CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (http://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   MULES::explicitSolve dimensions inconsistence (http://www.cfd-online.com/Forums/openfoam-programming-development/89898-mules-explicitsolve-dimensions-inconsistence.html)

santiagomarquezd June 24, 2011 19:29

MULES::explicitSolve dimensions inconsistence
 
Hi FOAMers I have a question about following lines in MULESTemplates.C

Code:

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 September 22, 2011 04:09

Hello Santiago!

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

Info << Su.field() << endl;
The error is sth like
Code:

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

bojiezhang November 7, 2011 10:09

Quote:

Originally Posted by santiagomarquezd (Post 313487)
Hi FOAMers I have a question about following lines in MULESTemplates.C

Code:

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 August 1, 2016 15:47

2 Attachment(s)
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 :cool:

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

Code:

{
// 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.

Code:

...
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:

Code:

...
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:

Code:

...
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


All times are GMT -4. The time now is 07:42.