multiphase and fireFoam: set phi based on the interpolated U before the pEqn construction for the p BCs

applications/solvers/combustion/fireFoam/pEqn.H

@@ -5,8 +5,8 @@ surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU));
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
 phi.boundaryField() =
-    fvc::interpolate(rho.boundaryField())
+    fvc::interpolate(rho.boundaryField()*U.boundaryField())
-   *(fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField());
+  & mesh.Sf().boundaryField();
 
 surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 

applications/solvers/multiphase/bubbleFoam/pEqn.H

@@ -14,12 +14,6 @@
     volVectorField HbyA2("HbyA2", U2);
     HbyA2 = rAU2*U2Eqn.H();
 
-    U.boundaryField() =
-        alpha1.boundaryField()*U1.boundaryField()
-      + alpha2.boundaryField()*U2.boundaryField();
-    phi.boundaryField() =
-        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
-
     surfaceScalarField phiDrag1
     (
         fvc::interpolate(alpha2/rho1*dragCoef*rAU1)*phi2
@@ -44,16 +38,18 @@
     (
         (fvc::interpolate(HbyA1) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU1, U1, phi1)
-      + phiDrag1
     );
 
     surfaceScalarField phiHbyA2
     (
         (fvc::interpolate(HbyA2) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU2, U2, phi2)
-      + phiDrag2
     );
 
+    phi = alpha1f*phiHbyA1 + alpha2f*phiHbyA2;
+
+    phiHbyA1 += phiDrag1;
+    phiHbyA2 += phiDrag2;
     surfaceScalarField phiHbyA("phiHbyA", alpha1f*phiHbyA1 + alpha2f*phiHbyA2);
 
     surfaceScalarField Dp

applications/solvers/multiphase/compressibleInterFoam/pEqn.H

@@ -29,8 +29,6 @@
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
-    phi.boundaryField() =
-        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
 
     surfaceScalarField phiHbyA
     (
@@ -38,6 +36,7 @@
         (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
+    phi = phiHbyA;
 
     surfaceScalarField phig
     (

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H

@@ -2,11 +2,6 @@
     rho1 = rho10 + psi1*p;
     rho2 = rho20 + psi2*p;
 
-    mrfZones.absoluteFlux(phi1.oldTime());
-    mrfZones.absoluteFlux(phi1);
-    mrfZones.absoluteFlux(phi2.oldTime());
-    mrfZones.absoluteFlux(phi2);
-
     surfaceScalarField alpha1f(fvc::interpolate(alpha1));
     surfaceScalarField alpha2f(scalar(1) - alpha1f);
 
@@ -22,11 +17,10 @@
     volVectorField HbyA2("HbyA2", U2);
     HbyA2 = rAU2*U2Eqn.H();
 
-    U.boundaryField() =
+    mrfZones.absoluteFlux(phi1.oldTime());
-        alpha1.boundaryField()*U1.boundaryField()
+    mrfZones.absoluteFlux(phi1);
-      + alpha2.boundaryField()*U2.boundaryField();
+    mrfZones.absoluteFlux(phi2.oldTime());
-    phi.boundaryField() =
+    mrfZones.absoluteFlux(phi2);
-        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
 
     surfaceScalarField ppDrag("ppDrag", 0.0*phi1);
 
@@ -47,18 +41,29 @@
         "phiHbyA1",
         (fvc::interpolate(HbyA1) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU1, alpha1, U1, phi1)
-      + fvc::interpolate((1.0/rho1)*rAU1*dragCoeff)*phi2
-      + ppDrag
-      + rAlphaAU1f*(g & mesh.Sf())
     );
-    mrfZones.relativeFlux(phiHbyA1);
 
     surfaceScalarField phiHbyA2
     (
         "phiHbyA2",
         (fvc::interpolate(HbyA2) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU2, alpha2, U2, phi2)
-      + fvc::interpolate((1.0/rho2)*rAU2*dragCoeff)*phi1
+    );
+
+    phi = alpha1f*phiHbyA1 + alpha2f*phiHbyA2;
+    mrfZones.relativeFlux(phi);
+
+    phiHbyA1 +=
+    (
+        fvc::interpolate((1.0/rho1)*rAU1*dragCoeff)*phi2
+      + ppDrag
+      + rAlphaAU1f*(g & mesh.Sf())
+    );
+    mrfZones.relativeFlux(phiHbyA1);
+
+    phiHbyA2 +=
+    (
+        fvc::interpolate((1.0/rho2)*rAU2*dragCoeff)*phi1
       + rAlphaAU2f*(g & mesh.Sf())
     );
     mrfZones.relativeFlux(phiHbyA2);

applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H

@@ -5,20 +5,15 @@
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
 
-    mrfZones.absoluteFlux(phi);
-    phi.boundaryField() =
-        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
-    mrfZones.relativeFlux(phi);
-
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
         (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
-    mrfZones.relativeFlux(phiHbyA);
-
     adjustPhi(phiHbyA, U, p_rgh);
+    mrfZones.relativeFlux(phiHbyA);
+    phi = phiHbyA;
 
     surfaceScalarField phig
     (

applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H

@@ -4,8 +4,6 @@
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
-    phi.boundaryField() =
-        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
 
     surfaceScalarField phiHbyA
     (
@@ -21,6 +19,8 @@
         fvc::makeAbsolute(phiHbyA, U);
     }
 
+    phiAbs = phiHbyA;
+
     surfaceScalarField phig
     (
         (

applications/solvers/multiphase/interFoam/pEqn.H

@@ -4,8 +4,6 @@
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
-    phi.boundaryField() =
-        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
 
     surfaceScalarField phiHbyA
     (
@@ -15,6 +13,7 @@
     );
 
     adjustPhi(phiHbyA, U, p_rgh);
+    phi = phiHbyA;
 
     surfaceScalarField phig
     (

applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H

@@ -4,8 +4,6 @@
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
-    phi.boundaryField() =
-        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
 
     surfaceScalarField phiHbyA
     (
@@ -13,8 +11,8 @@
         (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
-
     adjustPhi(phiHbyA, U, p_rgh);
+    phi = phiHbyA;
 
     surfaceScalarField phig
     (

applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H

@@ -58,6 +58,8 @@
         dimensionedScalar("phiHbyA", dimArea*dimVelocity, 0)
     );
 
+    phi = dimensionedScalar("phi", phi.dimensions(), 0);
+
     phasei = 0;
     forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
     {
@@ -74,14 +76,17 @@
         (
             (fvc::interpolate(HbyAs[phasei]) & mesh.Sf())
           + fvc::ddtPhiCorr(rAUs[phasei], alpha, phase.U(), phase.phi())
-          + rAlphaAUfs[phasei]
+        );
+        mrfZones.relativeFlux(phiHbyAs[phasei]);
+
+        phi += alphafs[phasei]*phiHbyAs[phasei];
+
+        phiHbyAs[phasei] +=
+            rAlphaAUfs[phasei]
            *(
                fluid.surfaceTension(phase)*mesh.magSf()/phase.rho()
              + (g & mesh.Sf())
-            )
+            );
-        );
-
-        mrfZones.relativeFlux(phiHbyAs[phasei]);
 
         multiphaseSystem::dragModelTable::const_iterator dmIter =
             fluid.dragModels().begin();

applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H

@@ -5,20 +5,15 @@
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
 
-    mrfZones.absoluteFlux(phi);
-    phi.boundaryField() =
-        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
-    mrfZones.relativeFlux(phi);
-
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
         (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
-    mrfZones.relativeFlux(phiHbyA);
-
     adjustPhi(phiHbyA, U, p_rgh);
+    mrfZones.relativeFlux(phiHbyA);
+    phi = phiHbyA;
 
     surfaceScalarField phig
     (

applications/solvers/multiphase/multiphaseInterFoam/pEqn.H

@@ -4,8 +4,6 @@
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
-    phi.boundaryField() =
-        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
 
     surfaceScalarField phiHbyA
     (
@@ -13,8 +11,8 @@
         (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
-
     adjustPhi(phiHbyA, U, p_rgh);
+    phi = phiHbyA;
 
     surfaceScalarField phig
     (

applications/solvers/multiphase/settlingFoam/pEqn.H

@@ -5,9 +5,6 @@
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
-    phi.boundaryField() =
-        fvc::interpolate(rho.boundaryField()*U.boundaryField())
-      & mesh.Sf().boundaryField();
 
     surfaceScalarField phiHbyA
     (
@@ -18,6 +15,7 @@
          + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
     );
+    phi = phiHbyA;
 
     surfaceScalarField phig
     (

applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H

@@ -4,8 +4,6 @@
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
-    phi.boundaryField() =
-        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
 
     surfaceScalarField phiHbyA
     (
@@ -13,8 +11,8 @@
         (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
-
     adjustPhi(phiHbyA, U, p_rgh);
+    phi = phiHbyA;
 
     surfaceScalarField phig
     (

applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H

@@ -18,13 +18,6 @@
     mrfZones.absoluteFlux(phi1);
     mrfZones.absoluteFlux(phi2.oldTime());
     mrfZones.absoluteFlux(phi2);
-    mrfZones.absoluteFlux(phi);
-
-    U.boundaryField() =
-        alpha1.boundaryField()*U1.boundaryField()
-      + alpha2.boundaryField()*U2.boundaryField();
-    phi.boundaryField() =
-        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
 
     surfaceScalarField ppDrag("ppDrag", 0.0*phi1);
 
@@ -43,18 +36,29 @@
         "phiHbyA1",
         (fvc::interpolate(HbyA1) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU1, U1, phi1)
-      + fvc::interpolate(alpha2/rho1*K*rAU1)*phi2
-      + ppDrag
-      + rAU1f*(g & mesh.Sf())
     );
-    mrfZones.relativeFlux(phiHbyA1);
 
     surfaceScalarField phiHbyA2
     (
         "phiHbyA2",
         (fvc::interpolate(HbyA2) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU2, U2, phi2)
-      + fvc::interpolate(alpha1/rho2*K*rAU2)*phi1
+    );
+
+    phi = alpha1f*phiHbyA1 + alpha2f*phiHbyA2;
+    mrfZones.relativeFlux(phi);
+
+    phiHbyA1 +=
+    (
+        fvc::interpolate(alpha2/rho1*K*rAU1)*phi2
+      + ppDrag
+      + rAU1f*(g & mesh.Sf())
+    );
+    mrfZones.relativeFlux(phiHbyA1);
+
+    phiHbyA2 +=
+    (
+        fvc::interpolate(alpha1/rho2*K*rAU2)*phi1
       + rAU2f*(g & mesh.Sf())
     );
     mrfZones.relativeFlux(phiHbyA2);
@@ -63,7 +67,6 @@
     mrfZones.relativeFlux(phi1);
     mrfZones.relativeFlux(phi2.oldTime());
     mrfZones.relativeFlux(phi2);
-    mrfZones.relativeFlux(phi);
 
     surfaceScalarField phiHbyA("phiHbyA", alpha1f*phiHbyA1 + alpha2f*phiHbyA2);