From f2ce1fa9ac5d8dab586b44abbf65413a68befb3d Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Fri, 28 Oct 2016 10:50:10 +0100
Subject: [PATCH] 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.
---
 .../twoPhaseSystem/phaseModel/phaseModel.C    | 26 +++++++++++++++++++
 .../twoPhaseSystem/phaseModel/phaseModel.H    |  3 +++
 .../twoPhaseSystem/twoPhaseSystem.C           | 24 ++---------------
 3 files changed, 31 insertions(+), 22 deletions(-)

diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.C
index a553a01b5e..cd60d2cbb0 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.C
+++ b/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];
+                }
+            }
+        }
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.H
index 5fb0245513..a5ae107c90 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.H
+++ b/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
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
index 7c2323a030..808fe255f9 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
+++ b/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();