MRG: Integrated Foundation code to commit 19e602b

applications/solvers/multiphase/VoF/alphaEqn.H

@@ -2,46 +2,57 @@
     word alphaScheme("div(phi,alpha)");
     word alpharScheme("div(phirb,alpha)");
 
-    tmp<fv::ddtScheme<scalar>> ddtAlpha
-    (
-        fv::ddtScheme<scalar>::New
-        (
-            mesh,
-            mesh.ddtScheme("ddt(alpha)")
-        )
-    );
-
     // Set the off-centering coefficient according to ddt scheme
     scalar ocCoeff = 0;
-    if
-    (
-        isType<fv::EulerDdtScheme<scalar>>(ddtAlpha())
-     || isType<fv::localEulerDdtScheme<scalar>>(ddtAlpha())
-    )
     {
-        ocCoeff = 0;
-    }
-    else if (isType<fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha()))
-    {
-        if (nAlphaSubCycles > 1)
+        tmp<fv::ddtScheme<scalar>> tddtAlpha
+        (
+            fv::ddtScheme<scalar>::New
+            (
+                mesh,
+                mesh.ddtScheme("ddt(alpha)")
+            )
+        );
+        const fv::ddtScheme<scalar>& ddtAlpha = tddtAlpha();
+
+        if
+        (
+            isType<fv::EulerDdtScheme<scalar>>(ddtAlpha)
+         || isType<fv::localEulerDdtScheme<scalar>>(ddtAlpha)
+        )
+        {
+            ocCoeff = 0;
+        }
+        else if (isType<fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha))
+        {
+            if (nAlphaSubCycles > 1)
+            {
+                FatalErrorInFunction
+                    << "Sub-cycling is not supported "
+                       "with the CrankNicolson ddt scheme"
+                    << exit(FatalError);
+            }
+
+            if
+            (
+                alphaRestart
+             || mesh.time().timeIndex() > mesh.time().startTimeIndex() + 1
+            )
+            {
+                ocCoeff =
+                    refCast<const fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha)
+                   .ocCoeff();
+            }
+        }
+        else
         {
             FatalErrorInFunction
-                << "Sub-cycling is not supported "
-                   "with the CrankNicolson ddt scheme"
+                << "Only Euler and CrankNicolson ddt schemes are supported"
                 << exit(FatalError);
         }
-
-        ocCoeff =
-            refCast<const fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha())
-           .ocCoeff();
-    }
-    else
-    {
-        FatalErrorInFunction
-            << "Only Euler and CrankNicolson ddt schemes are supported"
-            << exit(FatalError);
     }
 
+    // Set the time blending factor, 1 for Euler
     scalar cnCoeff = 1.0/(1.0 + ocCoeff);
 
     // Standard face-flux compression coefficient
@@ -136,8 +147,8 @@
         (
             fvc::flux
             (
-                phi,
-                alpha1,
+                phiCN(),
+                cnCoeff*alpha1 + (1.0 - cnCoeff)*alpha1.oldTime(),
                 alphaScheme
             )
           + fvc::flux
@@ -148,13 +159,6 @@
             )
         );
 
-        // Calculate the Crank-Nicolson off-centred alpha flux
-        if (ocCoeff > 0)
-        {
-            talphaPhiUn =
-                cnCoeff*talphaPhiUn + (1.0 - cnCoeff)*alphaPhi.oldTime();
-        }
-
         if (MULESCorr)
         {
             tmp<surfaceScalarField> talphaPhiCorr(talphaPhiUn() - alphaPhi);

applications/solvers/multiphase/VoF/createAlphaFluxes.H

@@ -0,0 +1,20 @@
+IOobject alphaPhiHeader
+(
+    "alphaPhi",
+    runTime.timeName(),
+    mesh,
+    IOobject::READ_IF_PRESENT,
+    IOobject::AUTO_WRITE
+);
+
+const bool alphaRestart = alphaPhiHeader.headerOk();
+
+// MULES flux from previous time-step
+surfaceScalarField alphaPhi
+(
+    alphaPhiHeader,
+    phi*fvc::interpolate(alpha1)
+);
+
+// MULES Correction
+tmp<surfaceScalarField> talphaPhiCorr0;

applications/solvers/multiphase/compressibleInterFoam/TEqn.H

@@ -5,7 +5,7 @@
       + fvm::div(rhoPhi, T)
       - fvm::laplacian(mixture.alphaEff(turbulence->mut()), T)
       + (
-            divU*p
+            fvc::div(fvc::absolute(phi, U), p)
           + fvc::ddt(rho, K) + fvc::div(rhoPhi, K)
         )
        *(

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

@@ -64,6 +64,7 @@ int main(int argc, char *argv[])
     #include "initContinuityErrs.H"
     #include "createControl.H"
     #include "createFields.H"
+    #include "createAlphaFluxes.H"
     #include "createFvOptions.H"
     #include "createUf.H"
     #include "createControls.H"

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C

@@ -64,6 +64,7 @@ int main(int argc, char *argv[])
     #include "createControl.H"
     #include "createTimeControls.H"
     #include "createFields.H"
+    #include "createAlphaFluxes.H"
     #include "createFvOptions.H"
 
     volScalarField& p = mixture.p();
@@ -112,7 +113,6 @@ int main(int argc, char *argv[])
             solve(fvm::ddt(rho) + fvc::div(rhoPhi));
 
             #include "UEqn.H"
-            volScalarField divU(fvc::div(fvc::absolute(phi, U)));
             #include "TEqn.H"
 
             // --- Pressure corrector loop

applications/solvers/multiphase/compressibleInterFoam/createFields.H

@@ -101,21 +101,4 @@ autoPtr<compressible::turbulenceModel> turbulence
 Info<< "Creating field kinetic energy K\n" << endl;
 volScalarField K("K", 0.5*magSqr(U));
 
-// MULES flux from previous time-step
-surfaceScalarField alphaPhi
-(
-    IOobject
-    (
-        "alphaPhi",
-        runTime.timeName(),
-        mesh,
-        IOobject::READ_IF_PRESENT,
-        IOobject::AUTO_WRITE
-    ),
-    phi*fvc::interpolate(alpha1)
-);
-
-// MULES Correction
-tmp<surfaceScalarField> talphaPhiCorr0;
-
 #include "createMRF.H"

applications/solvers/multiphase/interFoam/createFields.H

@@ -121,20 +121,6 @@ if (p_rgh.needReference())
 mesh.setFluxRequired(p_rgh.name());
 mesh.setFluxRequired(alpha1.name());
 
-// MULES flux from previous time-step
-surfaceScalarField alphaPhi
-(
-    IOobject
-    (
-        "alphaPhi",
-        runTime.timeName(),
-        mesh,
-        IOobject::READ_IF_PRESENT,
-        IOobject::AUTO_WRITE
-    ),
-    phi*fvc::interpolate(alpha1)
-);
-
 // MULES compressed flux is registered in case scalarTransport FO needs it.
 surfaceScalarField alphaPhiUn
 (
@@ -150,7 +136,4 @@ surfaceScalarField alphaPhiUn
     dimensionedScalar("zero", phi.dimensions(), 0.0)
 );
 
-// MULES Correction
-tmp<surfaceScalarField> talphaPhiCorr0;
-
 #include "createMRF.H"

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

@@ -63,6 +63,7 @@ int main(int argc, char *argv[])
     #include "createTimeControls.H"
     #include "createDyMControls.H"
     #include "createFields.H"
+    #include "createAlphaFluxes.H"
     #include "createFvOptions.H"
 
     volScalarField rAU

applications/solvers/multiphase/interFoam/interFoam.C

@@ -66,6 +66,7 @@ int main(int argc, char *argv[])
     #include "createTimeControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
+    #include "createAlphaFluxes.H"
     #include "createFvOptions.H"
     #include "correctPhi.H"
 

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H

@@ -294,12 +294,11 @@ while (pimple.correct())
     forAll(phases, phasei)
     {
         phaseModel& phase = phases[phasei];
+        const volScalarField& alpha = phase;
+        volScalarField& rho = phase.thermo().rho();
 
         if (phase.compressible())
         {
-            const volScalarField& alpha = phase;
-            const volScalarField& rho = phase.rho();
-
             if (pimple.transonic())
             {
                 surfaceScalarField phid
@@ -357,21 +356,29 @@ while (pimple.correct())
                     ).ptr()
                 );
             }
+        }
+        else
+        {
+            pEqnComps.set
+            (
+                phasei,
+                fvm::Su(-(fvOptions(alpha, rho)&rho)/rho, p_rgh).ptr()
+            );
+        }
 
-            if (fluid.transfersMass(phase))
+        if (fluid.transfersMass(phase))
+        {
+            if (pEqnComps.set(phasei))
             {
-                if (pEqnComps.set(phasei))
-                {
-                    pEqnComps[phasei] -= fluid.dmdt(phase)/rho;
-                }
-                else
-                {
-                    pEqnComps.set
-                    (
-                        phasei,
-                        fvm::Su(-fluid.dmdt(phase)/rho, p_rgh)
-                    );
-                }
+                pEqnComps[phasei] -= fluid.dmdt(phase)/rho;
+            }
+            else
+            {
+                pEqnComps.set
+                (
+                    phasei,
+                    fvm::Su(-fluid.dmdt(phase)/rho, p_rgh)
+                );
             }
         }
     }
@@ -421,7 +428,7 @@ while (pimple.correct())
                 phase.phi() = phiHbyAs[phasei] + alpharAUfs[phasei]*mSfGradp;
 
                 // Set the phase dilatation rates
-                if (phase.compressible())
+                if (pEqnComps.set(phasei))
                 {
                     phase.divU(-pEqnComps[phasei] & p_rgh);
                 }

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H

@@ -230,9 +230,10 @@ while (pimple.correct())
     tmp<fvScalarMatrix> pEqnComp2;
 
     // Construct the compressibility parts of the pressure equation
-    if (pimple.transonic())
+
+    if (phase1.compressible())
     {
-        if (phase1.compressible())
+        if (pimple.transonic())
         {
             surfaceScalarField phid1
             (
@@ -257,8 +258,24 @@ while (pimple.correct())
             deleteDemandDrivenData(pEqnComp1.ref().faceFluxCorrectionPtr());
             pEqnComp1.ref().relax();
         }
+        else
+        {
+            pEqnComp1 =
+                (
+                    phase1.continuityError()
+                  - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
+                )/rho1
+              + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh));
+        }
+    }
+    else
+    {
+        pEqnComp1 = fvm::Su(-(fvOptions(alpha1, rho1)&rho1)/rho1, p_rgh);
+    }
 
-        if (phase2.compressible())
+    if (phase2.compressible())
+    {
+        if (pimple.transonic())
         {
             surfaceScalarField phid2
             (
@@ -279,23 +296,11 @@ while (pimple.correct())
                       + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
                     )
                 );
+
             deleteDemandDrivenData(pEqnComp2.ref().faceFluxCorrectionPtr());
             pEqnComp2.ref().relax();
         }
-    }
-    else
-    {
-        if (phase1.compressible())
-        {
-            pEqnComp1 =
-                (
-                    phase1.continuityError()
-                  - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
-                )/rho1
-              + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh));
-        }
-
-        if (phase2.compressible())
+        else
         {
             pEqnComp2 =
                 (
@@ -305,6 +310,10 @@ while (pimple.correct())
               + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh));
         }
     }
+    else
+    {
+        pEqnComp2 = fvm::Su(-(fvOptions(alpha2, rho2)&rho2)/rho2, p_rgh);
+    }
 
     if (fluid.transfersMass())
     {
@@ -390,11 +399,11 @@ while (pimple.correct())
             }
 
             // Set the phase dilatation rates
-            if (phase1.compressible())
+            if (pEqnComp1.valid())
             {
                 phase1.divU(-pEqnComp1 & p_rgh);
             }
-            if (phase2.compressible())
+            if (pEqnComp2.valid())
             {
                 phase2.divU(-pEqnComp2 & p_rgh);
             }

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H

@@ -215,79 +215,91 @@ while (pimple.correct())
     tmp<fvScalarMatrix> pEqnComp1;
     tmp<fvScalarMatrix> pEqnComp2;
 
-    if (pimple.transonic())
-    {
-        surfaceScalarField phid1
-        (
-            IOobject::groupName("phid", phase1.name()),
-            fvc::interpolate(psi1)*phi1
-        );
-        surfaceScalarField phid2
-        (
-            IOobject::groupName("phid", phase2.name()),
-            fvc::interpolate(psi2)*phi2
-        );
+    // Construct the compressibility parts of the pressure equation
 
-        if (phase1.compressible())
+    if (phase1.compressible())
+    {
+        if (pimple.transonic())
         {
-            pEqnComp1 =
+            surfaceScalarField phid1
             (
-                phase1.continuityError()
-              - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
-            )/rho1
-          + correction
-            (
-                (alpha1/rho1)*
-                (
-                    psi1*fvm::ddt(p_rgh)
-                  + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
-                )
+                IOobject::groupName("phid", phase1.name()),
+                fvc::interpolate(psi1)*phi1
             );
+
+            pEqnComp1 =
+                (
+                    phase1.continuityError()
+                  - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
+                )/rho1
+              + correction
+                (
+                    (alpha1/rho1)*
+                    (
+                        psi1*fvm::ddt(p_rgh)
+                      + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
+                    )
+                );
+
             deleteDemandDrivenData(pEqnComp1.ref().faceFluxCorrectionPtr());
             pEqnComp1.ref().relax();
         }
-
-        if (phase2.compressible())
+        else
         {
-            pEqnComp2 =
-            (
-                phase2.continuityError()
-              - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
-            )/rho2
-          + correction
-            (
-                (alpha2/rho2)*
+            pEqnComp1 =
                 (
-                    psi2*fvm::ddt(p_rgh)
-                  + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
-                )
-            );
-            deleteDemandDrivenData(pEqnComp2.ref().faceFluxCorrectionPtr());
-            pEqnComp2.ref().relax();
+                    phase1.continuityError()
+                  - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
+                )/rho1
+              + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh));
         }
     }
     else
     {
-        if (phase1.compressible())
-        {
-            pEqnComp1 =
-            (
-                phase1.continuityError()
-              - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
-            )/rho1
-          + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh));
-        }
+        pEqnComp1 = fvm::Su(-(fvOptions(alpha1, rho1)&rho1)/rho1, p_rgh);
+    }
 
-        if (phase2.compressible())
+    if (phase2.compressible())
+    {
+        if (pimple.transonic())
+        {
+            surfaceScalarField phid2
+            (
+                IOobject::groupName("phid", phase2.name()),
+                fvc::interpolate(psi2)*phi2
+            );
+
+            pEqnComp2 =
+                (
+                    phase2.continuityError()
+                  - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
+                )/rho2
+              + correction
+                (
+                    (alpha2/rho2)*
+                    (
+                        psi2*fvm::ddt(p_rgh)
+                      + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
+                    )
+                );
+
+            deleteDemandDrivenData(pEqnComp2.ref().faceFluxCorrectionPtr());
+            pEqnComp2.ref().relax();
+        }
+        else
         {
             pEqnComp2 =
-            (
-                phase2.continuityError()
-              - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
-            )/rho2
-          + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh));
+                (
+                    phase2.continuityError()
+                  - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
+                )/rho2
+              + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh));
         }
     }
+    else
+    {
+        pEqnComp2 = fvm::Su(-(fvOptions(alpha2, rho2)&rho2)/rho2, p_rgh);
+    }
 
     if (fluid.transfersMass())
     {

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/derivedFvPatchFields/wallBoilingSubModels/partitioningModels/Lavieville/Lavieville.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -73,13 +73,13 @@ Lavieville::fLiquid
 ) const
 {
     return
-        pos(alphaLiquid-alphaCrit_)
+        pos(alphaLiquid - alphaCrit_)
        *(
             1 - 0.5*exp(-20*(alphaLiquid - alphaCrit_))
         )
       + neg(alphaLiquid - alphaCrit_)
        *(
-            pow(0.5*(alphaLiquid/alphaCrit_), 20*alphaCrit_)
+            0.5*pow(alphaLiquid/alphaCrit_, 20*alphaCrit_)
         );
 }