reactingEulerFoam: Added fvOption support for incompressible phases

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 =
-            (
-                phase1.continuityError()
-              - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
-            )/rho1
-          + correction
+            surfaceScalarField phid1
             (
-                (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 = fvm::Su(-(fvOptions(alpha1, rho1)&rho1)/rho1, p_rgh);
+    }
+
+    if (phase2.compressible())
+    {
+        if (pimple.transonic())
         {
-            pEqnComp1 =
+            surfaceScalarField phid2
             (
-                phase1.continuityError()
-              - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
-            )/rho1
-          + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh));
-        }
+                IOobject::groupName("phid", phase2.name()),
+                fvc::interpolate(psi2)*phi2
+            );
 
-        if (phase2.compressible())
+            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())
     {