[santiagomarquezd]

Hi FOAMers I have a question about following lines in MULESTemplates.C

00117     {
00118         psiIf =
00119         (
00120             rho.oldTime().field()*psi0/deltaT
00121           + Su.field()
00122           - psiIf
00123         )/(rho.field()/deltaT - Sp.field());
00124     }

they look like not dimensionally consistent, particularly:

rho.oldTime().field()*psi0/deltaT-psiIf

psiIf was previously used in and integration, nevertheless it wasn't multiplied by rho ever.

So, what's the trick?

Regards.

[lindstroem]

Hello Santiago!

I was considering the same question, have you found out sth? Do you know why i cant get an output via

Info << Su.field() << endl;

The error is sth like

MULESTemplates.C:132:5: error: no match for ‘operator<<’ in ‘Foam::operator<<(((F...

Greetings!

[bojiezhang]

Originally Posted by santiagomarquezd

Hi FOAMers I have a question about following lines in MULESTemplates.C

00117     {
00118         psiIf =
00119         (
00120             rho.oldTime().field()*psi0/deltaT
00121           + Su.field()
00122           - psiIf
00123         )/(rho.field()/deltaT - Sp.field());
00124     }

they look like not dimensionally consistent, particularly:

rho.oldTime().field()*psi0/deltaT-psiIf

psiIf was previously used in and integration, nevertheless it wasn't multiplied by rho ever.

So, what's the trick?

Regards.

hi:
I am also confused by this problem. Have you solved it?
bojiezhang

[jameswilson620]

Better late than never I suppose.

See attached.

The first image explains how for an incompressible flow using the FV method, we can use the divergence operator identity since div(u)=0 and surface integrals from Gauss divergence theorem to calculate the net accumulation/removal of alpha per unit time resulting from volume fluxes into the control volume (cell) enclosed by faces, dS.

actual, physical convective term:
u&grad(alpha) (1) (finite difference-esque expression, does not lend itself to local mass conservationfor the finite volume frame work)

The assumption that MULES::explicitSolve makes
div(alpha U) = u&grad(alpha) + alpha div(u) (2) (last term is zero <-> incompressible)

div(alpha U) is then used in place of Eq(1) for reasons that are explained below.

Gauss Divergence:
Vol_Int(div(alpha U))dV -> Surf_Int(alpha U&ns)dS = SUM_Over_faces_i(alpha_fi*U_fi*Area_fi) =
= SUM_Over_faces_i(alpha_fi*Phi_fi)

The last term above is calculated using surfaceIntegrate in MULES::explicitSolve where the units after integration are [m^3/s] (U&area ->[m^3/s]) and then become [1/s] after dividing the local cell result by the local cell volume.

Given the above explanation, and the equations in attached image 1, it can be deduced that if the flow is not incompressible, MULES::explicitSolve(..) will not return the correct result, i.e. there will be unboundedness/error in regions where div(phi)=/=0. If the user considers and understands the identities used to develop MULES::explicitSolve for an incompressible flow using surface integrals, then the user can add an implicit source term to correct for the inclusion of compressibility (div(phi)=/0).

In the first image, Eq(1) describes how a scalar transport equation in a compressible flow field SHOULD look while using div(phiAlpha) as the convective term, with no other external sources and including (-alphadiv(phi)) to correct for the conversion to div() for the convective term.

Eq(2) shows you what MULES::explicitSolve IS ACTUALLY SOLVING provided you supply:

MULES::explicitSolve(geometricOneField(),alpha,phi Alpha,Sp,Su,1,0)

where:
geometricOneField()->rho in image 1

rho can be anything that attenuates(>1)/advances(<1) the rate of accumulation/efflux of alpha in each cell. An example here would be isotropic porous media where the volume averaged darcian flux, phi_darcian, differs from the interstitial pore velocity, phi_interstitial, by the Dupuit-Forscheim relation phi_interstitial = phi_darcian/porosity. Since the interstitial pore velocity, phi_darcian/porosity, advects an interface and not phi_darcian, each convective term and source term can be modified by:
(convective and source terms)/porosity

-or- the transient term can be modified:
d(porosity*alpha)/dt

can be modified since porosity is assumed to be constant. This is the purpose of rho, it is constant and multiplies ONLY the transient term MULES::explicitSolve().

Finally
Sp->fvc::div(phi) (corrects the addition of alpha*div(phi) when converting convective term to div() in the presence of an incompressible flow in MULES::explicitSolve)

Su=0

It goes without saying that there are a number of 'incompressible' flows that might require the use of sources alpha*Sp and Su such as interPhaseChangeFoam where liquid cells are: 1 selected for removal as the pressure drops below Psat and 2: corrected as there is a local expansion/contraction of vapor resulting from an intended source/sink in the continuity equation, i.e. div(phi)=/=0.

With regard to the specific question asked here, look at image 1, Eq.2, then look at image 2. Connect the dots between Eq.1 in image 1, the manipulation of the explicit solution for the convection of alpha involving alpha*Sp and Su and the file mulesTemplates.C.

psiIf serves two purposes:

1. it is what returns to the solver the updated value of psi, the field representing alpha.

2. its a scalarField that is used to temporarily store the local values of surfaceIntegrate, i.e. the convective term for the VOF equation.

For example we want the new value of A, which is related to other values/variables by the relation:
A = 1+|b|;

where: b = -1;

To compute mag(b), we need to define a container to store it since the operation mag() isn't something that can be placed in the definition of A (theoretically); since A and b are of the same type, we'll temporarily use A to store our result for |b| to eliminate the need to initialize more variables just to store values:
A = |b|;

So finally, we have:
A = 1+|b| = 1+A

This is the same as what we see in MULES::explicitSolve() for psiIf.

EDITED BELOW: To add more detail and add the attachments I failed to attath the first time

Lets start with what we're all familiar with, alphaEqn.H:

{
// Add label for reference in fvSchemes
    word alphaScheme("div(phi,alpha)");                  
    word alpharScheme("div(phirb,alpha)");

// Define compressive velocity

// phic is defined as face velocities[m/s]
    surfaceScalarField phic(mag(phi/mesh.magSf()));

// Ensures that the compressive velocity is <= the max face velocity in the domain [m/s]
    phic = min(interface.cAlpha()*phic, max(phic));

// Takes phic and directs its components orthogonal to the interface using the interface normal nHatf [m/s]
    surfaceScalarField phir(phic*interface.nHatf());

// Couldnt tell you.. Can you tell me? Im guessing phiAlpha is updated at each loop since fvc::flux() is using old values of phi with new values of alpha1 when nAlphaCorr>1. Not sure why using phiAlpha calculated for the current time step solution is used to re-advect the old alpha1 is beneficial..
    for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
    {
    // Construct phiAlpha using the discretization schemes labeled above
    // Results in something like: alpha1*phi + alpha1*(1-alpha1)*phir
        surfaceScalarField phiAlpha
        (
            fvc::flux
            (
                phi,
                alpha1,
                alphaScheme
            )
          + fvc::flux
            (
                -fvc::flux(-phir, alpha2, alpharScheme),
                alpha1,
                alpharScheme
            )
        );
        
        // Store divergence. 0 if velocity field is solenoidal
        volScalarField divU = fvc::div(phi); 

        // Define source term fields
        volScalarField::DimensionedInternalField Su
        (
            IOobject
            (
                "Su",
                runTime.timeName(),
                mesh,
        	IOobject::NO_READ,
        	IOobject::NO_WRITE
            ),
	    mesh,
            dimensionedScalar("Su", dimensionSet(0,0,-1,0,0,0,0), 0.0)       
	);

        // Assign coefficient portion for the implicit source term field
        // Does nothing if flow field is solenoidal but corrects the MULES solver when there is a divergence in phi
        volScalarField::DimensionedInternalField Sp
        (
           IOobject
           (
               "Sp",
               runTime.timeName(),
               mesh,
               IOobject::NO_READ,
               IOobject::NO_WRITE
    	   ),
    	   divU
        );

     // CALL EXPLICIT SOLVER
         // See the snipett of code below for MULES::explicitSolve(...)
        MULES::explicitSolve(geometricOneField(), alpha1, phi, phiAlpha, Sp, Su, 1, 0); 
    // SEE BELOW FOR EXPLANATION
//      MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);

        alpha2 = 1.0 - alpha1;

        // Define the mixture mass flux for the momentum equation from the updated interface location
        rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
    }

    Info<< "Phase-1 volume fraction = "
        << alpha1.weightedAverage(mesh.Vsc()).value()
        << "  Min(alpha1) = " << min(alpha1).value()
        << "  Max(alpha1) = " << max(alpha1).value()
        << endl;
}

And for the function call, MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0), a call to ANOTHER function in MULES.C is made.

...
void Foam::MULES::explicitSolve  // we called this in alphaEqn.H 
(
    volScalarField& psi,
    const surfaceScalarField& phi,
    surfaceScalarField& phiPsi,
    const scalar psiMax,
    const scalar psiMin
)
{
    explicitSolve // This then gets called with the additional information required to get an explicit solution
    (
        geometricOneField(),  // ADDED since the transient term is present, but un-modified. hence multiplication by ones
        psi,
        phi,
        phiPsi,
        zeroField(), zeroField(),  // ADDED,ADDED since there are no source terms provided, i.e. Sp=Su=0
        psiMax, psiMin
    );
}
...

The first call to MULES::explicitSolve() filled in missing data required for an explicit solution and then called another MULES::explicitSolve() template while including the added information (above). The call to MULES::explicitSolve() now refers to a function in MULESTemplates.C:

...
template<class RhoType, class SpType, class SuType>
void Foam::MULES::explicitSolve
(
    const RhoType& rho,
    volScalarField& psi,
    const surfaceScalarField& phi,
    surfaceScalarField& phiPsi,  // This is where phiPsi updates phiAlpha in alphaEqn.H as limit() takes phiPsi and updates it
    const SpType& Sp,
    const SuType& Su,
    const scalar psiMax,
    const scalar psiMin
)
{
    const fvMesh& mesh = psi.mesh();
    const scalar rDeltaT = 1.0/mesh.time().deltaTValue(); // one over the time step
    psi.correctBoundaryConditions(); // Enforces boundary conditions (this is done after the solution too..:confused:)
    limit(rDeltaT, rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, 3, false);  // Flux corrected transport that makes MULES, MULES (function and constituents found later in MULESTemplates.C
    explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su);  // Final call to the explicit solver where we update the interface location (SEE BELOW)
}
...

Im omitting a discussion of limit() and limiter(). All i know is this is where the Flux Correction takes place and 'phiCorr' is defined. The surfaceScalarField phiPsi gets updated from these functions and results in a bounded solution when the sharp interface is advected using this field.

And finally, the last call to MULES::explicitSolve() in MULESTemplates.C:

...
template<class RhoType, class SpType, class SuType>
void Foam::MULES::explicitSolve
(
    const RhoType& rho,                  // Multiplies transient term in VOF equation (see image 2 attached)
    volScalarField& psi,                   // this is alpha. psiIf (defined below) updates psi, which based on the definition here, updates alpha in alphaEqn.H 
    const surfaceScalarField& phiPsi, // this is phiAlpha. Based on this definition, 
    const SpType& Sp,                    // This is the implicit source term added so the user can add whatever source or sink desired. 
                                                      // NOTE: when fvc::div(phi)=/=0 for your simulation, this term should be included in your call
                                                      // to: MULES::explicitSolve(....,Sp,..). Here Sp=fvc::div(phi) and NOT Sp= alpha*fvc::div(phi) 
                                                      // since the implicit part of the source term, i.e. the unknown dependent variable 'alpha', is 
                                                      // alreadyincluded in the MULES::explicitSolve() framework below
    const SuType& Su                     // This is similar to the above but is explicit, i.e. if you: 1. know how 
                                                      // much 'alpha'/second is appearing or disappearing and 2. is uncoupled 
                                                      // with alpha, i.e. its structure looks like 
)
{
    Info<< "MULES: Solving for " << psi.name() << endl; // This is the terminal output that exclaims: "MULES: Solving for alpha.water"

    const fvMesh& mesh = psi.mesh();

    scalarField& psiIf = psi;
    const scalarField& psi0 = psi.oldTime();
    const scalar deltaT = mesh.time().deltaTValue();

    psiIf = 0.0;

    // calculates the convective term [psiIf = surfaceIntegrate(phiPsi)] where
    // the first argument (psiIf) is the scalarField that stores the result and the 
    // second argument is the surfaceScalarField that is being integrated at each face for each cell in the domain
    fvc::surfaceIntegrate(psiIf,phiPsi) 

    if (mesh.moving())            // MOVING MESHES 
    {
        psiIf =
        (
            mesh.Vsc0()().field()*rho.oldTime().field()
           *psi0/(deltaT*mesh.Vsc()().field())
          + Su.field()
          - psiIf
        )/(rho.field()/deltaT - Sp.field());
    }
    else                               // STATIC MESH 
    {                                   // SEE ATTACHED IMAGE 2 TO SHOW HOW THIS WAS DERIVED 
        psiIf =                       // psiIf updates psi since: 'scalarField& psiIf = psi;'
        (
            rho.oldTime().field()*psi0/deltaT
          + Su.field()
          - psiIf                      // This is the convective term
        )/(rho.field()/deltaT - Sp.field());
    }

    psi.correctBoundaryConditions();
}

Hope this helps! -J

[mszeto715]

Hi,

I am wondering why the equation in image one starts as

ddt(rho alpha) + grad(alpha U) = Sp alpha + Su.

Shouldn't there be a rho in the second term? grad(rho alpha U)


Thank you.

-Mimi

[saturn_K]

Hello everyone

I have a question about the following lines in MULESTemplates.C

345         else
346         {
347             // Add the optional additional allowed boundary extrema
348             if (boundaryDeltaExtremaCoeff > 0)
349             {
350                 forAll(phiCorrPf, pFacei)
351                 {
352                     const label pfCelli = pFaceCells[pFacei];
353 
354                     const scalar extrema =
355                         boundaryDeltaExtremaCoeff
356                        *(psiMax[pfCelli] - psiMin[pfCelli]);
357 
358                     psiMaxn[pfCelli] += extrema;
359                     psiMinn[pfCelli] -= extrema;
360                 }
361             }
362         }

As you can see in the above lines, for psiMax and psiMin the subscript operator [] is used, and as far as I understand from other lines in this class both psiMax and psiMin are pointers, so basically we should be able to use the subscript operator for them but when I try to compile my solver which is included this class, I get the following error!

/opt/openfoam7/src/finiteVolume/lnInclude/MULESTemplates.C:356:32: error: subscripted value is neither array nor pointer
                        *(psiMax[pfCelli] - psiMin[pfCelli]);
                          ~~~~~~^
/opt/openfoam7/src/finiteVolume/lnInclude/MULESTemplates.C:356:50: error: subscripted value is neither array nor pointer
                        *(psiMax[pfCelli] - psiMin[pfCelli]);

Is there anybody who can explain to me exactly what is going on?!!! I am really confused and I will really appreciate it if you try to make it clear.

Regards,
Keivan