
[Sponsors] 
June 24, 2011, 19:29 
MULES::explicitSolve dimensions inconsistence

#1 
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 438
Rep Power: 17 
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 } rho.oldTime().field()*psi0/deltaTpsiIf 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: 543424511594 Int. 7032 Colectora Ruta Nac. 168 / Paraje El Pozo (3000) Santa Fe  Argentina. http://www.cimec.org.ar 

September 22, 2011, 04:09 

#2 
Senior Member
Join Date: Nov 2010
Posts: 113
Rep Power: 9 
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; Code:
MULESTemplates.C:132:5: error: no match for ‘operator<<’ in ‘Foam::operator<<(((F... Last edited by lindstroem; September 22, 2011 at 08:49. 

November 7, 2011, 10:09 

#3  
Member
bojiezhang
Join Date: Jan 2010
Posts: 64
Rep Power: 9 
Quote:
hi: I am also confused by this problem. Have you solved it? bojiezhang 

August 1, 2016, 15:47 

#4 
Member
james wilson
Join Date: Aug 2014
Location: Orlando, Fl
Posts: 38
Rep Power: 5 
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 differenceesque 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 DupuitForscheim 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 readvect 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*(1alpha1)*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<< "Phase1 volume fraction = " << alpha1.weightedAverage(mesh.Vsc()).value() << " Min(alpha1) = " << min(alpha1).value() << " Max(alpha1) = " << max(alpha1).value() << endl; } 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 unmodified. 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 ); } ... 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) } ... 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(); } Last edited by jameswilson620; August 1, 2016 at 19:09. Reason: Edited due to failure to attach images. 

Thread Tools  
Display Modes  


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 