From 69d994e794d7506f4af0be3885e8ca2e6ef8c623 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Thu, 25 Jun 2015 16:08:21 +0100
Subject: [PATCH] twoPhaseEulerFoam: Change the implicit particle-pressure and
 turbulence dispersion to be phase-symmetric so that the results are
 independent of which phase-fraction is solved.

---
 .../twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C    | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
index 2c686bd32d9..76af4a310a3 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
@@ -374,16 +374,19 @@ void Foam::twoPhaseSystem::solve()
     surfaceScalarField phic("phic", phi_);
     surfaceScalarField phir("phir", phi1 - phi2);
 
-    surfaceScalarField alpha1f(fvc::interpolate(max(alpha1, scalar(0))));
+    tmp<surfaceScalarField> alpha1alpha2f;
 
     if (pPrimeByA_.valid())
     {
+        alpha1alpha2f =
+            fvc::interpolate(max(alpha1, scalar(0)))
+           *fvc::interpolate(max(alpha2, scalar(0)));
+
         surfaceScalarField phiP
         (
             pPrimeByA_()*fvc::snGrad(alpha1, "bounded")*mesh_.magSf()
         );
 
-        phic += alpha1f*phiP;
         phir += phiP;
     }
 
@@ -520,7 +523,7 @@ void Foam::twoPhaseSystem::solve()
             fvScalarMatrix alpha1Eqn
             (
                 fvm::ddt(alpha1) - fvc::ddt(alpha1)
-              - fvm::laplacian(alpha1f*pPrimeByA_(), alpha1, "bounded")
+              - fvm::laplacian(alpha1alpha2f*pPrimeByA_(), alpha1, "bounded")
             );
 
             alpha1Eqn.relax();
-- 
GitLab