CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

MULES::explicitSolve dimensions inconsistence

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree2Likes
  • 2 Post By jameswilson620

Reply
 
LinkBack Thread Tools Display Modes
Old   June 24, 2011, 19:29
Default MULES::explicitSolve dimensions inconsistence
  #1
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 430
Rep Power: 16
santiagomarquezd will become famous soon enough
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.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Research Scientist
Research Center for Computational Methods (CIMEC) - CONICET/UNL
Tel: 54-342-4511594 Int. 7032
Colectora Ruta Nac. 168 / Paraje El Pozo
(3000) Santa Fe - Argentina.
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   September 22, 2011, 04:09
Default
  #2
Senior Member
 
Join Date: Nov 2010
Posts: 113
Rep Power: 7
lindstroem is on a distinguished road
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!

Last edited by lindstroem; September 22, 2011 at 08:49.
lindstroem is offline   Reply With Quote

Old   November 7, 2011, 10:09
Default
  #3
Member
 
bojiezhang
Join Date: Jan 2010
Posts: 64
Rep Power: 8
bojiezhang is on a distinguished road
Quote:
Originally Posted by santiagomarquezd View Post
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
bojiezhang is offline   Reply With Quote

Old   August 1, 2016, 15:47
Default
  #4
Member
 
james wilson
Join Date: Aug 2014
Location: Orlando, Fl
Posts: 38
Rep Power: 4
jameswilson620 is on a distinguished road
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:

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
Attached Images
File Type: jpg New Doc 4_1.jpg (154.9 KB, 1 views)
File Type: jpg New Doc 4_2.jpg (165.7 KB, 0 views)
akidess and GerhardHolzinger like this.

Last edited by jameswilson620; August 1, 2016 at 19:09. Reason: Edited due to failure to attach images.
jameswilson620 is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Incompatible dimensions.... Amiga500 OpenFOAM Running, Solving & CFD 13 June 1, 2012 07:20
fields + dimensions santiagomarquezd OpenFOAM Programming & Development 3 March 7, 2011 13:47
Dimensions ! T.D. OpenFOAM Running, Solving & CFD 4 September 24, 2010 14:26
Different dimensions for FATAL ERROR retech OpenFOAM Running, Solving & CFD 2 August 14, 2007 10:17
Dimensions of laplacian in PISO loop kumar2 OpenFOAM Running, Solving & CFD 2 July 3, 2006 14:34


All times are GMT -4. The time now is 00:05.