twoPhaseEulerFoam: Improvements to implicitPhasePressure

applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H

@@ -98,60 +98,3 @@
     Info<< "Creating field kinetic energy K\n" << endl;
     volScalarField K1(IOobject::groupName("K", phase1.name()), 0.5*magSqr(U1));
     volScalarField K2(IOobject::groupName("K", phase2.name()), 0.5*magSqr(U2));
-
-
-    volScalarField rAU1
-    (
-        IOobject
-        (
-            IOobject::groupName("rAU", phase1.name()),
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("zero", dimensionSet(-1, 3, 1, 0, 0), 0)
-    );
-
-    volScalarField rAU2
-    (
-        IOobject
-        (
-            IOobject::groupName("rAU", phase2.name()),
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("zero", dimensionSet(-1, 3, 1, 0, 0), 0)
-    );
-
-    surfaceScalarField rAU1f
-    (
-        IOobject
-        (
-            IOobject::groupName("rAUf", phase1.name()),
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("zero", dimensionSet(-1, 3, 1, 0, 0), 0)
-    );
-
-    surfaceScalarField rAU2f
-    (
-        IOobject
-        (
-            IOobject::groupName("rAUf", phase2.name()),
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("zero", dimensionSet(-1, 3, 1, 0, 0), 0)
-    );

applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H

@@ -1,11 +1,65 @@
-surfaceScalarField alpha1f("alpha1f", fvc::interpolate(alpha1));
+surfaceScalarField alphaf1("alphaf1", fvc::interpolate(alpha1));
-surfaceScalarField alpha2f("alpha2f", scalar(1) - alpha1f);
+surfaceScalarField alphaf2("alphaf2", scalar(1) - alphaf1);
 
-rAU1 = 1.0/U1Eqn.A();
+volScalarField rAU1(IOobject::groupName("rAU", phase1.name()), 1.0/U1Eqn.A());
-rAU2 = 1.0/U2Eqn.A();
+volScalarField rAU2(IOobject::groupName("rAU", phase2.name()), 1.0/U2Eqn.A());
+
+surfaceScalarField alpharAUf1(fvc::interpolate(alpha1*rAU1));
+surfaceScalarField alpharAUf2(fvc::interpolate(alpha2*rAU2));
+
+// Turbulent diffusion, particle-pressure, lift and wall-lubrication fluxes
+tmp<surfaceScalarField> phiF1;
+tmp<surfaceScalarField> phiF2;
+{
+    // Turbulent-dispersion diffusivity
+    volScalarField D(fluid.D());
+
+    // Phase-1 turbulent dispersion and particle-pressure flux
+    surfaceScalarField Df1
+    (
+        fvc::interpolate
+        (
+            rAU1*(D + phase1.turbulence().pPrime())
+        )
+    );
+
+    // Phase-2 turbulent dispersion and particle-pressure flux
+    surfaceScalarField Df2
+    (
+        fvc::interpolate
+        (
+            rAU2*(D + phase2.turbulence().pPrime())
+        )
+    );
+
+    // Cache the net diffusivity for implicit diffusion treatment in the
+    // phase-fraction equation
+    if (implicitPhasePressure)
+    {
+        fluid.pPrimeByA() = Df1 + Df2;
+    }
+
+    // Lift and wall-lubrication forces
+    volVectorField F(fluid.F());
+
+    // Phase-fraction face-gradient
+    surfaceScalarField snGradAlpha1(fvc::snGrad(alpha1)*mesh.magSf());
+
+    // Phase-1 dispersion, lift and wall-lubrication flux
+    phiF1 =
+    (
+        Df1*snGradAlpha1
+      + (fvc::interpolate(rAU1*F) & mesh.Sf())
+    );
+
+    // Phase-1 dispersion, lift and wall-lubrication flux
+    phiF2 =
+    (
+      - Df2*snGradAlpha1
+      - (fvc::interpolate(rAU2*F) & mesh.Sf())
+    );
+}
 
-surfaceScalarField alpharAU1f(fvc::interpolate(alpha1*rAU1));
-surfaceScalarField alpharAU2f(fvc::interpolate(alpha2*rAU2));
 
 // --- Pressure corrector loop
 while (pimple.correct())
@@ -33,61 +87,38 @@ while (pimple.correct())
     );
     HbyA2 = rAU2*U2Eqn.H();
 
-    // Face force fluxes
-    tmp<surfaceScalarField> phiF1;
-    tmp<surfaceScalarField> phiF2;
-
-    // Turbulent diffusion, particle-pressure lift and wall-lubrication fluxes
-    {
-        volScalarField D(fluid.D());
-        volVectorField F(fluid.F());
-        surfaceScalarField snGradAlpha1(fvc::snGrad(alpha1)*mesh.magSf());
-
-        phiF1 =
-        (
-            fvc::interpolate
-            (
-                rAU1*(D + phase1.turbulence().pPrime())
-            )*snGradAlpha1
-          + (fvc::interpolate(rAU1*F) & mesh.Sf())
-        );
-
-        phiF2 =
-        (
-          - fvc::interpolate
-            (
-                rAU2*(D + phase2.turbulence().pPrime())
-            )*snGradAlpha1
-          - (fvc::interpolate(rAU2*F) & mesh.Sf())
-        );
-    }
-
     // Mean density for buoyancy force and p_rgh -> p
     volScalarField rho("rho", fluid.rho());
 
-    // Add Buoyancy force
+    surfaceScalarField ghSnGradRho
-    {
+    (
-        surfaceScalarField ghSnGradRho(ghf*fvc::snGrad(rho)*mesh.magSf());
+        "ghSnGradRho",
+        ghf*fvc::snGrad(rho)*mesh.magSf()
+    );
 
-        phiF1() +=
+    surfaceScalarField phig1
-            alpharAU1f
+    (
-           *(
+        alpharAUf1
-                ghSnGradRho
+       *(
-              - alpha2f*fvc::interpolate(rho1 - rho2)*(g & mesh.Sf())
+           ghSnGradRho
-            );
+         - alphaf2*fvc::interpolate(rho1 - rho2)*(g & mesh.Sf())
+        )
+    );
+
+    surfaceScalarField phig2
+    (
+        alpharAUf2
+       *(
+           ghSnGradRho
+         - alphaf1*fvc::interpolate(rho2 - rho1)*(g & mesh.Sf())
+       )
+    );
 
-        phiF2() +=
-            alpharAU2f
-           *(
-                ghSnGradRho
-              - alpha1f*fvc::interpolate(rho2 - rho1)*(g & mesh.Sf())
-            );
-    }
 
     // ddtPhiCorr filter -- only apply in pure(ish) phases
-    surfaceScalarField alpha1fBar(fvc::interpolate(fvc::average(alpha1f)));
+    surfaceScalarField alphaf1Bar(fvc::interpolate(fvc::average(alphaf1)));
-    surfaceScalarField phiCorrCoeff1(pos(alpha1fBar - 0.99));
+    surfaceScalarField phiCorrCoeff1(pos(alphaf1Bar - 0.99));
-    surfaceScalarField phiCorrCoeff2(pos(0.01 - alpha1fBar));
+    surfaceScalarField phiCorrCoeff2(pos(0.01 - alphaf1Bar));
 
     forAll(mesh.boundary(), patchi)
     {
@@ -114,6 +145,7 @@ while (pimple.correct())
           - (fvc::interpolate(U1.oldTime()) & mesh.Sf())
         )/runTime.deltaT()
       - phiF1()
+      - phig1
     );
 
     // Phase-2 predicted flux
@@ -127,19 +159,20 @@ while (pimple.correct())
           - (fvc::interpolate(U2.oldTime()) & mesh.Sf())
         )/runTime.deltaT()
       - phiF2()
+      - phig2
     );
 
     // Face-drag coefficients
-    surfaceScalarField D1f(fvc::interpolate(rAU1*Kd));
+    surfaceScalarField rAUKd1(fvc::interpolate(rAU1*Kd));
-    surfaceScalarField D2f(fvc::interpolate(rAU2*Kd));
+    surfaceScalarField rAUKd2(fvc::interpolate(rAU2*Kd));
 
     // Construct the mean predicted flux
     // including explicit drag contributions based on absolute fluxes
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
-        alpha1f*(phiHbyA1 + D1f*mrfZones.absolute(phi2))
+        alphaf1*(phiHbyA1 + rAUKd1*mrfZones.absolute(phi2))
-      + alpha2f*(phiHbyA2 + D2f*mrfZones.absolute(phi1))
+      + alphaf2*(phiHbyA2 + rAUKd2*mrfZones.absolute(phi1))
     );
     mrfZones.makeRelative(phiHbyA);
 
@@ -147,7 +180,7 @@ while (pimple.correct())
     surfaceScalarField rAUf
     (
         "rAUf",
-        mag(alpha1f*alpharAU1f + alpha2f*alpharAU2f)
+        mag(alphaf1*alpharAUf1 + alphaf2*alpharAUf2)
     );
 
     // Update the fixedFluxPressure BCs to ensure flux consistency
@@ -158,9 +191,9 @@ while (pimple.correct())
             phiHbyA.boundaryField()
           - mrfZones.relative
             (
-                alpha1f.boundaryField()
+                alphaf1.boundaryField()
                *(mesh.Sf().boundaryField() & U1.boundaryField())
-              + alpha2f.boundaryField()
+              + alphaf2.boundaryField()
                *(mesh.Sf().boundaryField() & U2.boundaryField())
             )
         )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
@@ -256,21 +289,22 @@ while (pimple.correct())
             {
                 surfaceScalarField phi1s
                 (
-                    phiHbyA1 + alpharAU1f*mSfGradp
+                    phiHbyA1 + alpharAUf1*mSfGradp
                 );
 
                 surfaceScalarField phi2s
                 (
-                    phiHbyA2 + alpharAU2f*mSfGradp
+                    phiHbyA2 + alpharAUf2*mSfGradp
                 );
 
                 surfaceScalarField phir
                 (
-                    ((phi1s + D1f*phi2s) - (phi2s + D2f*phi1s))/(1 - D1f*D2f)
+                    ((phi1s + rAUKd1*phi2s) - (phi2s + rAUKd2*phi1s))
+                   /(1 - rAUKd1*rAUKd2)
                 );
 
-                phi1 = phi + alpha2f*phir;
+                phi1 = phi + alphaf2*phir;
-                phi2 = phi - alpha1f*phir;
+                phi2 = phi - alphaf1*phir;
             }
 
             // Compressibility correction for phase-fraction equations
@@ -287,21 +321,23 @@ while (pimple.correct())
 
             // Partial-elimination phase-velocity corrector
             {
-                volVectorField U1s
+                volVectorField Us1
                 (
-                    HbyA1 + fvc::reconstruct(alpharAU1f*mSfGradp - phiF1)
+                    HbyA1
+                  + fvc::reconstruct(alpharAUf1*mSfGradp - phiF1() - phig1)
                 );
 
-                volVectorField U2s
+                volVectorField Us2
                 (
-                    HbyA2 + fvc::reconstruct(alpharAU2f*mSfGradp - phiF2)
+                    HbyA2
+                  + fvc::reconstruct(alpharAUf2*mSfGradp - phiF2() - phig2)
                 );
 
                 volScalarField D1(rAU1*Kd);
                 volScalarField D2(rAU2*Kd);
 
-                U = alpha1*(U1s + D1*U2) + alpha2*(U2s + D2*U1);
+                U = alpha1*(Us1 + D1*U2) + alpha2*(Us2 + D2*U1);
-                volVectorField Ur(((1 - D2)*U1s - (1 - D1)*U2s)/(1 - D1*D2));
+                volVectorField Ur(((1 - D2)*Us1 - (1 - D1)*Us2)/(1 - D1*D2));
 
                 U1 = U + alpha2*Ur;
                 U1.correctBoundaryConditions();

applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H

@@ -1,9 +1,9 @@
-surfaceScalarField alpha1f("alpha1f", fvc::interpolate(alpha1));
+surfaceScalarField alphaf1("alphaf1", fvc::interpolate(alpha1));
-surfaceScalarField alpha2f("alpha2f", scalar(1) - alpha1f);
+surfaceScalarField alphaf2("alphaf2", scalar(1) - alphaf1);
 
-surfaceScalarField alphaRho1f0
+surfaceScalarField alphaRhof10
 (
-    "alphaRho1f0",
+    "alphaRhof10",
     fvc::interpolate
     (
         max(alpha1.oldTime(), fluid.residualAlpha(phase1))
@@ -11,9 +11,9 @@ surfaceScalarField alphaRho1f0
     )
 );
 
-surfaceScalarField alphaRho2f0
+surfaceScalarField alphaRhof20
 (
-    "alphaRho2f0",
+    "alphaRhof20",
     fvc::interpolate
     (
         max(alpha2.oldTime(), fluid.residualAlpha(phase2))
@@ -27,34 +27,76 @@ surfaceScalarField Kdf("Kdf", fluid.Kdf());
 // Virtual-mass coefficient
 surfaceScalarField Vmf("Vmf", fluid.Vmf());
 
-// Lift and wall-lubrication forces
+surfaceScalarField rAUf1
-surfaceScalarField Ff("Ff", fluid.Ff());
-
-tmp<surfaceScalarField> snGradAlpha1(fvc::snGrad(alpha1)*mesh.magSf());
-tmp<volScalarField> D(fluid.D());
-
-// Turbulent-dispersion force for phase 1
-surfaceScalarField Ftdf1
 (
-    IOobject::groupName("Ftdf", phase1.name()),
+    IOobject::groupName("rAUf", phase1.name()),
-    fvc::interpolate(D() + phase1.turbulence().pPrime())*snGradAlpha1()
+    1.0
+   /(
+        (alphaRhof10 + Vmf)/runTime.deltaT()
+      + fvc::interpolate(U1Eqn.A())
+      + Kdf
+   )
 );
 
-// Turbulent-dispersion force for phase 2
+surfaceScalarField rAUf2
-surfaceScalarField Ftdf2
 (
-    IOobject::groupName("Ftdf", phase2.name()),
+    IOobject::groupName("rAUf", phase2.name()),
-    fvc::interpolate(D + phase2.turbulence().pPrime())*snGradAlpha1
+    1.0
+   /(
+        (alphaRhof20 + Vmf)/runTime.deltaT()
+      + fvc::interpolate(U2Eqn.A())
+      + Kdf
+   )
 );
 
 
+// Turbulent dispersion, particle-pressure, lift and wall-lubrication forces
+tmp<surfaceScalarField> Ff1;
+tmp<surfaceScalarField> Ff2;
+{
+    // Turbulent-dispersion diffusivity
+    volScalarField D(fluid.D());
+
+    // Phase-1 turbulent dispersion and particle-pressure diffusivity
+    surfaceScalarField Df1
+    (
+        fvc::interpolate(D + phase1.turbulence().pPrime())
+    );
+
+    // Phase-2 turbulent dispersion and particle-pressure diffusivity
+    surfaceScalarField Df2
+    (
+        fvc::interpolate(D + phase2.turbulence().pPrime())
+    );
+
+    // Cache the net diffusivity for implicit diffusion treatment in the
+    // phase-fraction equation
+    if (implicitPhasePressure)
+    {
+        fluid.pPrimeByA() = rAUf1*Df1 + rAUf2*Df2;
+    }
+
+    // Lift and wall-lubrication forces
+    surfaceScalarField Ff(fluid.Ff());
+
+    // Phase-fraction face-gradient
+    surfaceScalarField snGradAlpha1(fvc::snGrad(alpha1)*mesh.magSf());
+
+    // Phase-1 dispersion, lift and wall-lubrication force
+    Ff1 = Df1*snGradAlpha1 + Ff;
+
+    // Phase-2 dispersion, lift and wall-lubrication force
+    Ff2 = -Df2*snGradAlpha1 - Ff;
+}
+
+
 while (pimple.correct())
 {
     // Update continuity errors due to temperature changes
     #include "correctContErrs.H"
 
-    surfaceScalarField rho1f(fvc::interpolate(rho1));
+    surfaceScalarField rhof1(fvc::interpolate(rho1));
-    surfaceScalarField rho2f(fvc::interpolate(rho2));
+    surfaceScalarField rhof2(fvc::interpolate(rho2));
 
     // Correct flux BCs to be consistent with the velocity BCs
     phi1.boundaryField() ==
@@ -62,41 +104,18 @@ while (pimple.correct())
     phi2.boundaryField() ==
         mrfZones.relative(mesh.Sf().boundaryField() & U2.boundaryField());
 
-    rAU1f =
+    surfaceScalarField alpharAUf1
     (
-        IOobject::groupName("rAUf", phase1.name()),
+        IOobject::groupName("alpharAUf", phase1.name()),
-        1.0
+        max(alphaf1, fluid.residualAlpha(phase1))*rAUf1
-       /(
-           (alphaRho1f0 + Vmf)/runTime.deltaT()
-         + fvc::interpolate(U1Eqn.A())
-         + Kdf
-        )
     );
 
-    rAU2f =
+    surfaceScalarField alpharAUf2
     (
-        IOobject::groupName("rAUf", phase2.name()),
+        IOobject::groupName("alpharAUf", phase2.name()),
-        1.0
+        max(alphaf2, fluid.residualAlpha(phase2))*rAUf2
-       /(
-           (alphaRho2f0 + Vmf)/runTime.deltaT()
-         + fvc::interpolate(U2Eqn.A())
-         + Kdf
-        )
     );
 
-    surfaceScalarField rAlphaAU1f
-    (
-        IOobject::groupName("rAlphaAUf", phase1.name()),
-        max(alpha1f, fluid.residualAlpha(phase1))*rAU1f
-    );
-
-    surfaceScalarField rAlphaAU2f
-    (
-        IOobject::groupName("rAlphaAUf", phase2.name()),
-        max(alpha2f, fluid.residualAlpha(phase2))*rAU2f
-    );
-
-
     volScalarField rho("rho", fluid.rho());
     surfaceScalarField ghSnGradRho
     (
@@ -104,25 +123,25 @@ while (pimple.correct())
         ghf*fvc::snGrad(rho)*mesh.magSf()
     );
 
-    // Add the phase-1 buoyancy force
+    // Phase-1 buoyancy flux
-    surfaceScalarField phiF1
+    surfaceScalarField phig1
     (
-        IOobject::groupName("phiF", phase1.name()),
+        IOobject::groupName("phig", phase1.name()),
-        rAlphaAU1f
+        alpharAUf1
        *(
             ghSnGradRho
-          - alpha2f*(rho1f - rho2f)*(g & mesh.Sf())
+          - alphaf2*(rhof1 - rhof2)*(g & mesh.Sf())
         )
     );
 
-    // Add the phase-2 buoyancy force
+    // Phase-2 buoyancy flux
-    surfaceScalarField phiF2
+    surfaceScalarField phig2
     (
-        IOobject::groupName("phiF", phase2.name()),
+        IOobject::groupName("phig", phase2.name()),
-        rAlphaAU2f
+        alpharAUf2
        *(
             ghSnGradRho
-          - alpha1f*(rho2f - rho1f)*(g & mesh.Sf())
+          - alphaf1*(rhof2 - rhof1)*(g & mesh.Sf())
         )
     );
 
@@ -135,15 +154,14 @@ while (pimple.correct())
     );
 
     phiHbyA1 =
-        rAU1f
+        rAUf1
        *(
-            (alphaRho1f0 + Vmf)
+            (alphaRhof10 + Vmf)
            *mrfZones.absolute(phi1.oldTime())/runTime.deltaT()
           + (fvc::interpolate(U1Eqn.H()) & mesh.Sf())
           + Vmf*ddtPhi2
           + Kdf*mrfZones.absolute(phi2)
-          - Ff
+          - Ff1()
-          - Ftdf1
         );
 
     // Phase-2 predicted flux
@@ -154,32 +172,31 @@ while (pimple.correct())
     );
 
     phiHbyA2 =
-        rAU2f
+        rAUf2
        *(
-            (alphaRho2f0 + Vmf)
+            (alphaRhof20 + Vmf)
            *mrfZones.absolute(phi2.oldTime())/runTime.deltaT()
           + (fvc::interpolate(U2Eqn.H()) & mesh.Sf())
           + Vmf*ddtPhi1
           + Kdf*mrfZones.absolute(phi1)
-          + Ff
+          - Ff2()
-          + Ftdf2
        );
 
 
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
-        alpha1f*(phiHbyA1 - phiF1) + alpha2f*(phiHbyA2 - phiF2)
+        alphaf1*(phiHbyA1 - phig1) + alphaf2*(phiHbyA2 - phig2)
     );
     mrfZones.makeRelative(phiHbyA);
 
-    phiHbyA1 -= phiF1;
+    phiHbyA1 -= phig1;
-    phiHbyA2 -= phiF2;
+    phiHbyA2 -= phig2;
 
     surfaceScalarField rAUf
     (
         "rAUf",
-        mag(alpha1f*rAlphaAU1f + alpha2f*rAlphaAU2f)
+        mag(alphaf1*alpharAUf1 + alphaf2*alpharAUf2)
     );
 
     // Update the fixedFluxPressure BCs to ensure flux consistency
@@ -189,8 +206,8 @@ while (pimple.correct())
         (
             phiHbyA.boundaryField()
           - (
-                alpha1f.boundaryField()*phi1.boundaryField()
+                alphaf1.boundaryField()*phi1.boundaryField()
-              + alpha2f.boundaryField()*phi2.boundaryField()
+              + alphaf2.boundaryField()*phi2.boundaryField()
             )
         )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
     );
@@ -280,25 +297,25 @@ while (pimple.correct())
             surfaceScalarField phi1s
             (
                 phiHbyA1
-              + rAlphaAU1f*mSfGradp
+              + alpharAUf1*mSfGradp
-              - rAU1f*Kdf*mrfZones.absolute(phi2)
+              - rAUf1*Kdf*mrfZones.absolute(phi2)
             );
 
             surfaceScalarField phi2s
             (
                 phiHbyA2
-              + rAlphaAU2f*mSfGradp
+              + alpharAUf2*mSfGradp
-              - rAU2f*Kdf*mrfZones.absolute(phi1)
+              - rAUf2*Kdf*mrfZones.absolute(phi1)
             );
 
             surfaceScalarField phir
             (
-                ((phi2s + rAU2f*Kdf*phi1s) - (phi1s + rAU1f*Kdf*phi2s))
+                ((phi2s + rAUf2*Kdf*phi1s) - (phi1s + rAUf1*Kdf*phi2s))
-               /(1.0 - rAU1f*rAU2f*sqr(Kdf))
+               /(1.0 - rAUf1*rAUf2*sqr(Kdf))
             );
 
-            phi1 = phi - alpha2f*phir;
+            phi1 = phi - alphaf2*phir;
-            phi2 = phi + alpha1f*phir;
+            phi2 = phi + alphaf1*phir;
 
             U1 = fvc::reconstruct(mrfZones.absolute(phi1));
             U1.correctBoundaryConditions();

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C

@@ -50,11 +50,6 @@ int main(int argc, char *argv[])
 
     pimpleControl pimple(mesh);
 
-    Switch faceMomentum
-    (
-        pimple.dict().lookupOrDefault<Switch>("faceMomentum", false)
-    );
-
     #include "createFields.H"
     #include "createMRFZones.H"
     #include "createFvOptions.H"
@@ -63,6 +58,19 @@ int main(int argc, char *argv[])
     #include "CourantNos.H"
     #include "setInitialDeltaT.H"
 
+    Switch faceMomentum
+    (
+        pimple.dict().lookupOrDefault<Switch>("faceMomentum", false)
+    );
+
+    Switch implicitPhasePressure
+    (
+        mesh.solverDict(alpha1.name()).lookupOrDefault<Switch>
+        (
+            "implicitPhasePressure", false
+        )
+    );
+
     #include "pUf/createDDtU.H"
     #include "pU/createDDtU.H"
 

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C

@@ -365,13 +365,6 @@ void Foam::twoPhaseSystem::solve()
     const surfaceScalarField& phi1 = phase1_.phi();
     const surfaceScalarField& phi2 = phase2_.phi();
 
-    Switch faceMomentum
-    (
-        //pimple.dict().lookupOrDefault<Switch>("faceMomentum", false)
-        mesh_.solutionDict().subDict("PIMPLE")
-           .lookupOrDefault<Switch>("faceMomentum", false)
-    );
-
     const dictionary& alphaControls = mesh_.solverDict
     (
         alpha1.name()
@@ -379,10 +372,6 @@ void Foam::twoPhaseSystem::solve()
 
     label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
     label nAlphaCorr(readLabel(alphaControls.lookup("nAlphaCorr")));
-    Switch implicitPhasePressure
-    (
-        alphaControls.lookupOrDefault<Switch>("implicitPhasePressure", false)
-    );
 
     word alphaScheme("div(phi," + alpha1.name() + ')');
     word alpharScheme("div(phir," + alpha1.name() + ')');
@@ -395,48 +384,11 @@ void Foam::twoPhaseSystem::solve()
 
     surfaceScalarField alpha1f(fvc::interpolate(max(alpha1, scalar(0))));
 
-    tmp<surfaceScalarField> pPrimeByA;
+    if (pPrimeByA_.valid())
-
-    if (implicitPhasePressure)
     {
-        if (faceMomentum)
-        {
-            const surfaceScalarField& rAU1f =
-                mesh_.lookupObject<surfaceScalarField>
-                (
-                    IOobject::groupName("rAUf", phase1_.name())
-                );
-            const surfaceScalarField& rAU2f =
-                mesh_.lookupObject<surfaceScalarField>
-                (
-                    IOobject::groupName("rAUf", phase2_.name())
-                );
-
-            volScalarField D(this->D());
-
-            pPrimeByA =
-                rAU1f*fvc::interpolate(D + phase1_.turbulence().pPrime())
-              + rAU2f*fvc::interpolate(D + phase2_.turbulence().pPrime());
-        }
-        else
-        {
-            const volScalarField& rAU1 = mesh_.lookupObject<volScalarField>
-            (
-                IOobject::groupName("rAU", phase1_.name())
-            );
-            const volScalarField& rAU2 = mesh_.lookupObject<volScalarField>
-            (
-                IOobject::groupName("rAU", phase2_.name())
-            );
-
-            pPrimeByA =
-                fvc::interpolate(rAU1*phase1_.turbulence().pPrime())
-              + fvc::interpolate(rAU2*phase2_.turbulence().pPrime());
-        }
-
         surfaceScalarField phiP
         (
-            pPrimeByA()*fvc::snGrad(alpha1, "bounded")*mesh_.magSf()
+            pPrimeByA_()*fvc::snGrad(alpha1, "bounded")*mesh_.magSf()
         );
 
         phic += alpha1f*phiP;
@@ -571,12 +523,12 @@ void Foam::twoPhaseSystem::solve()
             phase1_.alphaPhi() = alphaPhic1;
         }
 
-        if (implicitPhasePressure)
+        if (pPrimeByA_.valid())
         {
             fvScalarMatrix alpha1Eqn
             (
                 fvm::ddt(alpha1) - fvc::ddt(alpha1)
-              - fvm::laplacian(alpha1f*pPrimeByA(), alpha1, "bounded")
+              - fvm::laplacian(alpha1f*pPrimeByA_(), alpha1, "bounded")
             );
 
             alpha1Eqn.relax();

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H

@@ -78,9 +78,12 @@ class twoPhaseSystem
         //- Total volumetric flux
         surfaceScalarField phi_;
 
-        //-  Dilatation term
+        //- Dilatation term
         volScalarField dgdt_;
 
+        //- Optional dispersion diffusivity
+        tmp<surfaceScalarField> pPrimeByA_;
+
         //- Unordered phase pair
         autoPtr<phasePair> pair_;
 
@@ -113,7 +116,7 @@ class twoPhaseSystem
         autoPtr<BlendedInterfacialModel<turbulentDispersionModel> >
             turbulentDispersion_;
 
-        //-
+        //- Default residual alpha (0)
         static dimensionedScalar zeroResidualAlpha_;
 
 
@@ -227,6 +230,9 @@ public:
 
             //- Return non-const access to the dilatation parameter
             inline volScalarField& dgdt();
+
+            //- Return non-const access to the dispersion diffusivity
+            inline tmp<surfaceScalarField>& pPrimeByA();
 };
 
 

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystemI.H

@@ -95,4 +95,10 @@ inline Foam::volScalarField& Foam::twoPhaseSystem::dgdt()
 }
 
 
+inline Foam::tmp<Foam::surfaceScalarField>& Foam::twoPhaseSystem::pPrimeByA()
+{
+    return pPrimeByA_;
+}
+
+
 // ************************************************************************* //