From efad85b0e67b0e1b3f19ddcf93c17d9b7f5a720e Mon Sep 17 00:00:00 2001
From: Will Bainbridge <http://cfd.direct>
Date: Wed, 11 Apr 2018 14:50:13 +0100
Subject: [PATCH] reactingEulerFoam: Mixture properties for moving phases only

The calculations for mixture rho and U have been changed so that they
represent phase-averaged quantities over the moving phases only.

The mixture density is used as part of the pressure solution to
calculate buoyancy forces. The pressure within a stationary phase is
considered to be decoupled from the moving phases; i.e., it is
considered self-supporting. Therefore the stationary phase density
should not form a part of buoyancy calculations. This change to the
definition of mixture density ensures this.
---
 .../phaseSystems/phaseSystem/phaseSystem.C    | 50 +++++++++++++------
 1 file changed, 36 insertions(+), 14 deletions(-)

diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
index 8585f0f15..46eb6988e 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
@@ -237,33 +237,55 @@ Foam::phaseSystem::~phaseSystem()
 
 Foam::tmp<Foam::volScalarField> Foam::phaseSystem::rho() const
 {
-    tmp<volScalarField> tmpRho
-    (
-        phaseModels_[0]*phaseModels_[0].rho()
-    );
+    const label nMovingPhases = movingPhaseModels_.size();
 
-    for (label phasei=1; phasei<phaseModels_.size(); phasei++)
+    volScalarField rho(movingPhaseModels_[0]*movingPhaseModels_[0].rho());
+    for (label movingPhasei = 1; movingPhasei < nMovingPhases; ++ movingPhasei)
     {
-        tmpRho.ref() += phaseModels_[phasei]*phaseModels_[phasei].rho();
+        rho +=
+            movingPhaseModels_[movingPhasei]
+           *movingPhaseModels_[movingPhasei].rho();
     }
 
-    return tmpRho;
+    if (stationaryPhaseModels_.empty())
+    {
+        return rho;
+    }
+
+    volScalarField alpha(movingPhaseModels_[0]);
+    for (label movingPhasei = 1; movingPhasei < nMovingPhases; ++ movingPhasei)
+    {
+        alpha += movingPhaseModels_[movingPhasei];
+    }
+
+    return rho/alpha;
 }
 
 
 Foam::tmp<Foam::volVectorField> Foam::phaseSystem::U() const
 {
-    tmp<volVectorField> tmpU
-    (
-        phaseModels_[0]*phaseModels_[0].U()
-    );
+    const label nMovingPhases = movingPhaseModels_.size();
 
-    for (label phasei=1; phasei<phaseModels_.size(); phasei++)
+    volVectorField U(movingPhaseModels_[0]*movingPhaseModels_[0].U());
+    for (label movingPhasei = 1; movingPhasei < nMovingPhases; ++ movingPhasei)
     {
-        tmpU.ref() += phaseModels_[phasei]*phaseModels_[phasei].U();
+        U +=
+            movingPhaseModels_[movingPhasei]
+           *movingPhaseModels_[movingPhasei].U();
     }
 
-    return tmpU;
+    if (stationaryPhaseModels_.empty())
+    {
+        return U;
+    }
+
+    volScalarField alpha(movingPhaseModels_[0]);
+    for (label movingPhasei = 1; movingPhasei < nMovingPhases; ++ movingPhasei)
+    {
+        alpha += movingPhaseModels_[movingPhasei];
+    }
+
+    return U/alpha;
 }