twoPhaseEulerFoam and compressibleTwoPhaseEulerFoam: Improved flux formulation for BCs

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H

@@ -2,6 +2,51 @@
     rho1 = rho10 + psi1*p;
     rho2 = rho20 + psi2*p;
 
+    surfaceScalarField alpha1f(fvc::interpolate(alpha1));
+    surfaceScalarField alpha2f(scalar(1) - alpha1f);
+
+    volScalarField rAU1(1.0/U1Eqn.A());
+    volScalarField rAU2(1.0/U2Eqn.A());
+
+    surfaceScalarField rAlphaAU1f = fvc::interpolate(alpha1*rAU1);
+    surfaceScalarField rAlphaAU2f = fvc::interpolate(alpha2*rAU2);
+
+    volVectorField HbyA1("HbyA1", U1);
+    HbyA1 = rAU1*U1Eqn.H();
+
+    volVectorField HbyA2("HbyA2", U2);
+    HbyA2 = rAU2*U2Eqn.H();
+
+    surfaceScalarField phiHbyA1
+    (
+        "phiHbyA1",
+        (fvc::interpolate(HbyA1) & mesh.Sf())
+      + fvc::ddtPhiCorr(rAU1, alpha1, U1, phi1)
+      + fvc::interpolate((1.0/rho1)*rAU1*dragCoeff)*phi2
+      + rAlphaAU1f*(g & mesh.Sf())
+    );
+
+    surfaceScalarField phiHbyA2
+    (
+        "phiHbyA2",
+        (fvc::interpolate(HbyA2) & mesh.Sf())
+      + fvc::ddtPhiCorr(rAU2, alpha2, U2, phi2)
+      + fvc::interpolate((1.0/rho2)*rAU2*dragCoeff)*phi1
+      + rAlphaAU2f*(g & mesh.Sf())
+    );
+
+    surfaceScalarField phiHbyA("phiHbyA", alpha1f*phiHbyA1 + alpha2f*phiHbyA2);
+
+    surfaceScalarField Dp
+    (
+        "Dp",
+        mag
+        (
+            alpha1f*rAlphaAU1f/fvc::interpolate(rho1)
+          + alpha2f*rAlphaAU2f/fvc::interpolate(rho2)
+        )
+    );
+
     tmp<fvScalarMatrix> pEqnComp1;
     tmp<fvScalarMatrix> pEqnComp2;
 
@@ -24,57 +69,11 @@
           - fvc::Sp(fvc::div(phid2), p);
     }
 
-    surfaceScalarField alpha1f(fvc::interpolate(alpha1));
-    surfaceScalarField alpha2f(scalar(1) - alpha1f);
-
-    volVectorField U10 = U1;
-    volVectorField U20 = U2;
-
-    volScalarField rAU1(1.0/U1Eqn.A());
-    volScalarField rAU2(1.0/U2Eqn.A());
-
-    surfaceScalarField rAlphaAU1f = fvc::interpolate(alpha1*rAU1);
-    surfaceScalarField rAlphaAU2f = fvc::interpolate(alpha2*rAU2);
-
-    U1 = rAU1*U1Eqn.H();
-    U2 = rAU2*U2Eqn.H();
-
-    {
-        surfaceScalarField phi10 = phi1;
-
-        phi1 =
-        (
-            (fvc::interpolate(U1) & mesh.Sf())
-          + fvc::ddtPhiCorr(rAU1, alpha1, U1, phi1)
-          + fvc::interpolate((1.0/rho1)*rAU1*dragCoeff)*phi2
-          + rAlphaAU1f*(g & mesh.Sf())
-        );
-
-        phi2 =
-        (
-            (fvc::interpolate(U2) & mesh.Sf())
-          + fvc::ddtPhiCorr(rAU2, alpha2, U2, phi2)
-          + fvc::interpolate((1.0/rho2)*rAU2*dragCoeff)*phi10
-          + rAlphaAU2f*(g & mesh.Sf())
-        );
-    }
-
-    surfaceScalarField phi0("phi0", alpha1f*phi1 + alpha2f*phi2);
-    phi = phi0;
-    adjustPhi(phi, U, p);
-
-    surfaceScalarField Dp
-    (
-        "Dp",
-        mag(alpha1f*rAlphaAU1f/fvc::interpolate(rho1)
-      + alpha2f*rAlphaAU2f/fvc::interpolate(rho2))
-    );
-
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqnIncomp
         (
-            fvc::div(phi)
+            fvc::div(phiHbyA)
           - fvm::laplacian(Dp, p)
         );
 
@@ -85,13 +84,21 @@
               + (alpha2/rho2)*pEqnComp2()
             )
           + pEqnIncomp,
-            mesh.solver(p.select(pimple.finalInnerIter())));
+            mesh.solver(p.select(pimple.finalInnerIter()))
+        );
 
         if (pimple.finalNonOrthogonalIter())
         {
             surfaceScalarField mSfGradp = pEqnIncomp.flux()/Dp;
-            phi1 += rAlphaAU1f*mSfGradp/fvc::interpolate(rho1);
-            phi2 += rAlphaAU2f*mSfGradp/fvc::interpolate(rho2);
+
+            phi1.boundaryField() ==
+                (fvc::interpolate(U1) & mesh.Sf())().boundaryField();
+            phi1 = phiHbyA1 + rAlphaAU1f*mSfGradp/fvc::interpolate(rho1);
+
+            phi2.boundaryField() ==
+                (fvc::interpolate(U2) & mesh.Sf())().boundaryField();
+            phi2 = phiHbyA2 + rAlphaAU2f*mSfGradp/fvc::interpolate(rho2);
+
             phi = alpha1f*phi1 + alpha2f*phi2;
 
             dgdt =
@@ -103,7 +110,10 @@
             p.relax();
             mSfGradp = pEqnIncomp.flux()/Dp;
 
-            U1 += (1.0/rho1)*rAU1*dragCoeff*U20
+            volVectorField U10 = U1;
+
+            U1 = HbyA1
+               + (1.0/rho1)*rAU1*dragCoeff*U2
                + fvc::reconstruct
                  (
                      rAlphaAU1f*(g & mesh.Sf())
@@ -111,7 +121,8 @@
                  );
             U1.correctBoundaryConditions();
 
-            U2 += (1.0/rho2)*rAU2*dragCoeff*U10
+            U2 = HbyA2
+               + (1.0/rho2)*rAU2*dragCoeff*U10
                + fvc::reconstruct
                  (
                      rAlphaAU2f*(g & mesh.Sf())

applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H

@@ -5,11 +5,6 @@
     volScalarField rUaA(1.0/UaEqn.A());
     volScalarField rUbA(1.0/UbEqn.A());
 
-    phia.boundaryField() ==
-        (fvc::interpolate(Ua) & mesh.Sf())().boundaryField();
-    phib.boundaryField() ==
-        (fvc::interpolate(Ub) & mesh.Sf())().boundaryField();
-
     rUaAf = fvc::interpolate(rUaA);
     surfaceScalarField rUbAf(fvc::interpolate(rUbA));
 
@@ -51,6 +46,7 @@
 
     surfaceScalarField phiHbyAa
     (
+        "phiHbyAa",
         (fvc::interpolate(HbyAa) & mesh.Sf())
       + fvc::ddtPhiCorr(rUaA, Ua, phia)
       + phiDraga
@@ -58,6 +54,7 @@
 
     surfaceScalarField phiHbyAb
     (
+        "phiHbyAb",
         (fvc::interpolate(HbyAb) & mesh.Sf())
       + fvc::ddtPhiCorr(rUbA, Ub, phib)
       + phiDragb
@@ -86,8 +83,14 @@
         {
             surfaceScalarField SfGradp(pEqn.flux()/Dp);
 
+            phia.boundaryField() ==
+                (fvc::interpolate(Ua) & mesh.Sf())().boundaryField();
             phia = phiHbyAa - rUaAf*SfGradp/rhoa;
+
+            phib.boundaryField() ==
+                (fvc::interpolate(Ub) & mesh.Sf())().boundaryField();
             phib = phiHbyAb - rUbAf*SfGradp/rhob;
+
             phi = alphaf*phia + betaf*phib;
 
             p.relax();