compressibleInterFoam: Improved mass conservation

using the continuity error correction formulation developed for
twoPhaseEulerFoam and reactingEulerFoam.

applications/solvers/multiphase/VoF/alphaEqn.H

@@ -113,6 +113,8 @@
                 phiCN,
                 upwind<scalar>(mesh, phiCN)
             ).fvmDiv(phiCN, alpha1)
+       // - fvm::Sp(fvc::ddt(dimensionedScalar("1", dimless, 1), mesh)
+       //           + fvc::div(phiCN), alpha1)
          ==
             Su + fvm::Sp(Sp + divU, alpha1)
         );
@@ -125,19 +127,19 @@
             << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
             << endl;
 
-        tmp<surfaceScalarField> talphaPhiUD(alpha1Eqn.flux());
-        alphaPhi = talphaPhiUD();
+        tmp<surfaceScalarField> talphaPhi1UD(alpha1Eqn.flux());
+        alphaPhi10 = talphaPhi1UD();
 
-        if (alphaApplyPrevCorr && talphaPhiCorr0.valid())
+        if (alphaApplyPrevCorr && talphaPhi1Corr0.valid())
         {
             Info<< "Applying the previous iteration compression flux" << endl;
-            MULES::correct(alpha1, alphaPhi, talphaPhiCorr0.ref(), 1, 0);
+            MULES::correct(alpha1, alphaPhi10, talphaPhi1Corr0.ref(), 1, 0);
 
-            alphaPhi += talphaPhiCorr0();
+            alphaPhi10 += talphaPhi1Corr0();
         }
 
         // Cache the upwind-flux
-        talphaPhiCorr0 = talphaPhiUD;
+        talphaPhi1Corr0 = talphaPhi1UD;
 
         alpha2 = 1.0 - alpha1;
 
@@ -151,7 +153,7 @@
 
         surfaceScalarField phir(phic*mixture.nHatf());
 
-        tmp<surfaceScalarField> talphaPhiUn
+        tmp<surfaceScalarField> talphaPhi1Un
         (
             fvc::flux
             (
@@ -169,15 +171,15 @@
 
         if (MULESCorr)
         {
-            tmp<surfaceScalarField> talphaPhiCorr(talphaPhiUn() - alphaPhi);
+            tmp<surfaceScalarField> talphaPhi1Corr(talphaPhi1Un() - alphaPhi10);
             volScalarField alpha10("alpha10", alpha1);
 
             MULES::correct
             (
                 geometricOneField(),
                 alpha1,
-                talphaPhiUn(),
-                talphaPhiCorr.ref(),
+                talphaPhi1Un(),
+                talphaPhi1Corr.ref(),
                 Sp,
                 (-Sp*alpha1)(),
                 1,
@@ -187,24 +189,24 @@
             // Under-relax the correction for all but the 1st corrector
             if (aCorr == 0)
             {
-                alphaPhi += talphaPhiCorr();
+                alphaPhi10 += talphaPhi1Corr();
             }
             else
             {
                 alpha1 = 0.5*alpha1 + 0.5*alpha10;
-                alphaPhi += 0.5*talphaPhiCorr();
+                alphaPhi10 += 0.5*talphaPhi1Corr();
             }
         }
         else
         {
-            alphaPhi = talphaPhiUn;
+            alphaPhi10 = talphaPhi1Un;
 
             MULES::explicitSolve
             (
                 geometricOneField(),
                 alpha1,
                 phiCN,
-                alphaPhi,
+                alphaPhi10,
                 Sp,
                 (Su + divU*min(alpha1(), scalar(1)))(),
                 1,
@@ -219,34 +221,37 @@
 
     if (alphaApplyPrevCorr && MULESCorr)
     {
-        talphaPhiCorr0 = alphaPhi - talphaPhiCorr0;
-        talphaPhiCorr0.ref().rename("alphaPhiCorr0");
+        talphaPhi1Corr0 = alphaPhi10 - talphaPhi1Corr0;
+        talphaPhi1Corr0.ref().rename("alphaPhi1Corr0");
     }
     else
     {
-        talphaPhiCorr0.clear();
+        talphaPhi1Corr0.clear();
     }
 
+    #include "rhofs.H"
+
     if
     (
         word(mesh.ddtScheme("ddt(rho,U)"))
      == fv::EulerDdtScheme<vector>::typeName
+     || word(mesh.ddtScheme("ddt(rho,U)"))
+     == fv::localEulerDdtScheme<vector>::typeName
     )
     {
-        #include "rhofs.H"
-        rhoPhi = alphaPhi*(rho1f - rho2f) + phiCN*rho2f;
+        rhoPhi = alphaPhi10*(rho1f - rho2f) + phiCN*rho2f;
     }
     else
     {
         if (ocCoeff > 0)
         {
             // Calculate the end-of-time-step alpha flux
-            alphaPhi = (alphaPhi - (1.0 - cnCoeff)*alphaPhi.oldTime())/cnCoeff;
+            alphaPhi10 =
+                (alphaPhi10 - (1.0 - cnCoeff)*alphaPhi10.oldTime())/cnCoeff;
         }
 
         // Calculate the end-of-time-step mass flux
-        #include "rhofs.H"
-        rhoPhi = alphaPhi*(rho1f - rho2f) + phi*rho2f;
+        rhoPhi = alphaPhi10*(rho1f - rho2f) + phi*rho2f;
     }
 
     Info<< "Phase-1 volume fraction = "

applications/solvers/multiphase/VoF/createAlphaFluxes.H

-IOobject alphaPhiHeader
+IOobject alphaPhi10Header
 (
-    "alphaPhi",
+    "alphaPhi10",
     runTime.timeName(),
     mesh,
     IOobject::READ_IF_PRESENT,
     IOobject::AUTO_WRITE
 );
 
-const bool alphaRestart = alphaPhiHeader.typeHeaderOk<surfaceScalarField>(true);
+const bool alphaRestart =
+    alphaPhi10Header.typeHeaderOk<surfaceScalarField>(true);
 
 // MULES flux from previous time-step
-surfaceScalarField alphaPhi
+surfaceScalarField alphaPhi10
 (
-    alphaPhiHeader,
+    alphaPhi10Header,
     phi*fvc::interpolate(alpha1)
 );
 
 // MULES Correction
-tmp<surfaceScalarField> talphaPhiCorr0;
+tmp<surfaceScalarField> talphaPhi1Corr0;

applications/solvers/multiphase/compressibleInterFoam/TEqn.H

 {
     fvScalarMatrix TEqn
     (
-        fvm::ddt(rho, T)
-      + fvm::div(rhoPhi, T)
+        fvm::ddt(rho, T) + fvm::div(rhoPhi, T)
+      - fvm::Sp(contErr, T)
       - fvm::laplacian(mixture.alphaEff(turbulence->mut()), T)
       + (
             divU*p

applications/solvers/multiphase/compressibleInterFoam/UEqn.H

     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
+      - fvm::Sp(contErr, U)
       + MRF.DDt(rho, U)
       + turbulence->divDevRhoReff(U)
      ==

applications/solvers/multiphase/compressibleInterFoam/alphaSuSp.H

@@ -24,14 +24,14 @@ volScalarField::Internal Su
 
 forAll(dgdt, celli)
 {
-    if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
+    if (dgdt[celli] > 0.0)
     {
-        Sp[celli] -= dgdt[celli]*alpha1[celli];
-        Su[celli] += dgdt[celli]*alpha1[celli];
+        Sp[celli] -= dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
+        Su[celli] += dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
     }
-    else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
+    else if (dgdt[celli] < 0.0)
     {
-        Sp[celli] += dgdt[celli]*(1.0 - alpha1[celli]);
+        Sp[celli] += dgdt[celli]/max(alpha1[celli], 1e-4);
     }
 }
 

applications/solvers/multiphase/compressibleInterFoam/compressibleAlphaEqnSubCycle.H

+tmp<surfaceScalarField> talphaPhi1(alphaPhi10);
+
+if (nAlphaSubCycles > 1)
 {
-    #include "alphaEqnSubCycle.H"
+    dimensionedScalar totalDeltaT = runTime.deltaT();
+
+    talphaPhi1 = new surfaceScalarField
+    (
+        IOobject
+        (
+            "alphaPhi1",
+            runTime.timeName(),
+            mesh
+        ),
+        mesh,
+        dimensionedScalar("0", alphaPhi10.dimensions(), 0)
+    );
+
+    surfaceScalarField rhoPhiSum
+    (
+        IOobject
+        (
+            "rhoPhiSum",
+            runTime.timeName(),
+            mesh
+        ),
+        mesh,
+        dimensionedScalar("0", rhoPhi.dimensions(), 0)
+    );
+
+    tmp<volScalarField> trSubDeltaT;
+
+    if (LTS)
+    {
+        trSubDeltaT =
+            fv::localEulerDdt::localRSubDeltaT(mesh, nAlphaSubCycles);
+    }
+
+    for
+    (
+        subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
+        !(++alphaSubCycle).end();
+    )
+    {
+        #include "alphaEqn.H"
+        talphaPhi1.ref() += (runTime.deltaT()/totalDeltaT)*alphaPhi10;
+        rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
+    }
+
+    rhoPhi = rhoPhiSum;
 }
+else
+{
+    #include "alphaEqn.H"
+}
+
+rho == alpha1*rho1 + alpha2*rho2;
+
+const surfaceScalarField& alphaPhi1 = talphaPhi1();
+surfaceScalarField alphaPhi2("alphaPhi2", phi - alphaPhi1);
+
+volScalarField::Internal contErr
+(
+    (
+        fvc::ddt(rho) + fvc::div(rhoPhi)
+      - (fvOptions(alpha1, mixture.thermo1().rho())&rho1)
+      - (fvOptions(alpha2, mixture.thermo2().rho())&rho2)
+    )()
+);

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C

@@ -150,8 +150,6 @@ int main(int argc, char *argv[])
             #include "alphaControls.H"
             #include "compressibleAlphaEqnSubCycle.H"
 
-            solve(fvm::ddt(rho) + fvc::div(rhoPhi));
-
             #include "UEqn.H"
             #include "TEqn.H"
 

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H

@@ -56,13 +56,29 @@
     }
     else
     {
+        #include "rhofs.H"
+
         p_rghEqnComp1 =
-            fvc::ddt(rho1) + psi1*correction(fvm::ddt(p_rgh))
-          + fvc::div(phi, rho1) - fvc::Sp(fvc::div(phi), rho1);
+            pos(alpha1)
+           *(
+                (
+                    fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f)
+                  - (fvOptions(alpha1, mixture.thermo1().rho())&rho1)
+                )/rho1
+              - fvc::ddt(alpha1) - fvc::div(alphaPhi1)
+              + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh))
+            );
 
         p_rghEqnComp2 =
-            fvc::ddt(rho2) + psi2*correction(fvm::ddt(p_rgh))
-          + fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2);
+            pos(alpha2)
+           *(
+               (
+                   fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2*rho2f)
+                 - (fvOptions(alpha2, mixture.thermo2().rho())&rho2)
+               )/rho2
+             - fvc::ddt(alpha2) - fvc::div(alphaPhi2)
+             + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh))
+            );
     }
 
     // Cache p_rgh prior to solve for density update
@@ -78,16 +94,21 @@
 
         solve
         (
-            (
-                (max(alpha1, scalar(0))/rho1)*p_rghEqnComp1()
-              + (max(alpha2, scalar(0))/rho2)*p_rghEqnComp2()
-            )
-          + p_rghEqnIncomp,
+            p_rghEqnComp1() + p_rghEqnComp2() + p_rghEqnIncomp,
             mesh.solver(p_rgh.select(pimple.finalInnerIter()))
         );
 
         if (pimple.finalNonOrthogonalIter())
         {
+            p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
+            p_rgh = p - (alpha1*rho1 + alpha2*rho2)*gh;
+
+            dgdt =
+            (
+                alpha1*(p_rghEqnComp2 & p_rgh)
+              - alpha2*(p_rghEqnComp1 & p_rgh)
+            );
+
             phi = phiHbyA + p_rghEqnIncomp.flux();
 
             U = HbyA
@@ -109,15 +130,9 @@
 
     rho = alpha1*rho1 + alpha2*rho2;
 
-    p = max(p_rgh + rho*gh, pMin);
+    // Correct p_rgh for consistency with p and the updated densities
     p_rgh = p - rho*gh;
     p_rgh.correctBoundaryConditions();
 
-    dgdt =
-    (
-        pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2
-      - pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1
-    );
-
     K = 0.5*magSqr(U);
 }

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C

@@ -108,9 +108,7 @@ int main(int argc, char *argv[])
         while (pimple.loop())
         {
             #include "alphaControls.H"
-            #include "alphaEqnSubCycle.H"
-
-            solve(fvm::ddt(rho) + fvc::div(rhoPhi));
+            #include "compressibleAlphaEqnSubCycle.H"
 
             #include "UEqn.H"
             volScalarField divU(fvc::div(fvc::absolute(phi, U)));

applications/solvers/multiphase/compressibleInterFoam/createFields.H

@@ -89,7 +89,7 @@ surfaceScalarField rhoPhi
 
 volScalarField dgdt
 (
-    pos(alpha2)*fvc::div(phi)/max(alpha2, scalar(0.0001))
+    pos0(alpha2)*fvc::div(phi)/max(alpha2, scalar(0.0001))
 );
 
 // Construct compressible turbulence model

applications/solvers/multiphase/compressibleInterFoam/pEqn.H

@@ -53,13 +53,29 @@
     }
     else
     {
+        #include "rhofs.H"
+
         p_rghEqnComp1 =
-            fvc::ddt(rho1) + psi1*correction(fvm::ddt(p_rgh))
-          + fvc::div(phi, rho1) - fvc::Sp(fvc::div(phi), rho1);
+            pos(alpha1)
+           *(
+                (
+                    fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f)
+                  - (fvOptions(alpha1, mixture.thermo1().rho())&rho1)
+                )/rho1
+              - fvc::ddt(alpha1) - fvc::div(alphaPhi1)
+              + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh))
+            );
 
         p_rghEqnComp2 =
-            fvc::ddt(rho2) + psi2*correction(fvm::ddt(p_rgh))
-          + fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2);
+            pos(alpha2)
+           *(
+               (
+                   fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2*rho2f)
+                 - (fvOptions(alpha2, mixture.thermo2().rho())&rho2)
+               )/rho2
+             - fvc::ddt(alpha2) - fvc::div(alphaPhi2)
+             + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh))
+            );
     }
 
     // Cache p_rgh prior to solve for density update
@@ -75,11 +91,7 @@
 
         solve
         (
-            (
-                (max(alpha1, scalar(0))/rho1)*p_rghEqnComp1()
-              + (max(alpha2, scalar(0))/rho2)*p_rghEqnComp2()
-            )
-          + p_rghEqnIncomp,
+            p_rghEqnComp1() + p_rghEqnComp2() + p_rghEqnIncomp,
             mesh.solver(p_rgh.select(pimple.finalInnerIter()))
         );
 
@@ -90,8 +102,8 @@
 
             dgdt =
             (
-                pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2
-              - pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1
+                alpha1*(p_rghEqnComp2 & p_rgh)
+              - alpha2*(p_rghEqnComp1 & p_rgh)
             );
 
             phi = phiHbyA + p_rghEqnIncomp.flux();

applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -950,7 +950,7 @@ Foam::multiphaseMixtureThermo::nearInterface() const
     forAllConstIter(PtrDictionary<phaseModel>, phases_, phase)
     {
         tnearInt.ref() =
-            max(tnearInt(), pos(phase() - 0.01)*pos(0.99 - phase()));
+            max(tnearInt(), pos0(phase() - 0.01)*pos0(0.99 - phase()));
     }
 
     return tnearInt;

applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H

@@ -104,7 +104,7 @@
             )
             {
                 phase().dgdt() =
-                    pos(phase())
+                    pos0(phase())
                   *(p_rghEqnComps[phasei] & p_rgh)/phase().thermo().rho();
             }
 

applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C

@@ -132,7 +132,7 @@ int main(int argc, char *argv[])
                     // if the mesh topology changed
                     if (mesh.topoChanging())
                     {
-                        talphaPhiCorr0.clear();
+                        talphaPhi1Corr0.clear();
                     }
 
                     gh = (g & mesh.C()) - ghRef;

src/finiteVolume/fvMesh/fvMeshGeometry.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -434,7 +434,7 @@ const Foam::surfaceScalarField& Foam::fvMesh::phi() const
 
     // Set zero current time
     // mesh motion fluxes if the time has been incremented
-    if (phiPtr_->timeIndex() != time().timeIndex())
+    if (!time().subCycling() && phiPtr_->timeIndex() != time().timeIndex())
     {
         (*phiPtr_) = dimensionedScalar("0", dimVolume/dimTime, 0.0);
     }

tutorials/multiphase/compressibleInterDyMFoam/laminar/sloshingTank2D/system/fvSchemes

@@ -31,8 +31,6 @@ divSchemes
     div(phirb,alpha) Gauss linear;
 
     div(rhoPhi,U)  Gauss vanLeerV;
-    div(phi,thermo:rho.water) Gauss linear;
-    div(phi,thermo:rho.air) Gauss linear;
     div(rhoPhi,T)  Gauss linear;
     div(rhoPhi,K)  Gauss linear;
     div((phi+meshPhi),p)  Gauss linear;

tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/system/fvSchemes

@@ -31,8 +31,6 @@ divSchemes
     div(phirb,alpha) Gauss linear;
 
     div(rhoPhi,U)  Gauss upwind;
-    div(phi,thermo:rho.water) Gauss upwind;
-    div(phi,thermo:rho.air) Gauss upwind;
     div(rhoPhi,T)  Gauss upwind;
     div(rhoPhi,K)  Gauss upwind;
     div(phi,p)      Gauss upwind;

tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/system/fvSchemes

@@ -31,8 +31,6 @@ divSchemes
     div(phirb,alpha) Gauss linear;
 
     div(rhoPhi,U)  Gauss upwind;
-    div(phi,thermo:rho.water) Gauss upwind;
-    div(phi,thermo:rho.air) Gauss upwind;
     div(rhoPhi,T)  Gauss upwind;
     div(rhoPhi,K)  Gauss upwind;
     div(phi,p)      Gauss upwind;