twoPhaseEulerFoam::twoPhaseSystem: Ensure inlet flow of BOTH phases matches the BCs

Previously the inlet flow of phase 1 (the phase solved for) is corrected
to match the inlet specification for that phase.  However, if the second
phase is also constrained at inlets the inlet flux must also be
corrected to match the inlet specification.

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

@@ -248,4 +248,30 @@ bool Foam::phaseModel::read(const dictionary& phaseProperties)
 }
 
 
+void Foam::phaseModel::correctInflowFlux(surfaceScalarField& alphaPhi) const
+{
+    surfaceScalarField::Boundary& alphaPhiBf = alphaPhi.boundaryFieldRef();
+
+    // Ensure that the flux at inflow BCs is preserved
+    forAll(alphaPhiBf, patchi)
+    {
+        fvsPatchScalarField& alphaPhip = alphaPhiBf[patchi];
+
+        if (!alphaPhip.coupled())
+        {
+            const scalarField& phip = phi().boundaryField()[patchi];
+            const scalarField& alphap = boundaryField()[patchi];
+
+            forAll(alphaPhip, facei)
+            {
+                if (phip[facei] < SMALL)
+                {
+                    alphaPhip[facei] = alphap[facei]*phip[facei];
+                }
+            }
+        }
+    }
+}
+
+
 // ************************************************************************* //

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

@@ -319,6 +319,9 @@ public:
             return alphaRhoPhi_;
         }
 
+        //- Ensure that the flux at inflow BCs is preserved
+        void correctInflowFlux(surfaceScalarField& alphaPhi) const;
+
         //- Correct the phase properties
         //  other than the thermodynamics and turbulence
         //  which have special treatment

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

@@ -444,28 +444,7 @@ void Foam::twoPhaseSystem::solve()
             )
         );
 
-        surfaceScalarField::Boundary& alphaPhic1Bf =
-            alphaPhic1.boundaryFieldRef();
-
-        // Ensure that the flux at inflow BCs is preserved
-        forAll(alphaPhic1Bf, patchi)
-        {
-            fvsPatchScalarField& alphaPhic1p = alphaPhic1Bf[patchi];
-
-            if (!alphaPhic1p.coupled())
-            {
-                const scalarField& phi1p = phi1.boundaryField()[patchi];
-                const scalarField& alpha1p = alpha1.boundaryField()[patchi];
-
-                forAll(alphaPhic1p, facei)
-                {
-                    if (phi1p[facei] < 0)
-                    {
-                        alphaPhic1p[facei] = alpha1p[facei]*phi1p[facei];
-                    }
-                }
-            }
-        }
+        phase1_.correctInflowFlux(alphaPhic1);
 
         if (nAlphaSubCycles > 1)
         {
@@ -537,6 +516,7 @@ void Foam::twoPhaseSystem::solve()
 
         phase2_.alphaPhi() = phi_ - phase1_.alphaPhi();
         alpha2 = scalar(1) - alpha1;
+        phase2_.correctInflowFlux(phase2_.alphaPhi());
         phase2_.alphaRhoPhi() =
             fvc::interpolate(phase2_.rho())*phase2_.alphaPhi();