From d00c51e3b73636fae3a85c9f81942098f6e07629 Mon Sep 17 00:00:00 2001
From: "Lehnigk, Ronald (FWDC) - 10659" <r.lehnigk@hzdr.de>
Date: Fri, 4 Jan 2019 12:54:34 +0100
Subject: [PATCH] reactingMultiphaseEulerFoam: Fixes a bug where input order
 affects results

Prior to this commit, the order in which the phases involved are specified
through the entry "phases" in the phaseProperties dictionary affected the
results in case of an compressible simulation.

Patch contributed by Institute of Fluid Dynamics,
Helmholtz-Zentrum Dresden - Rossendorf (HZDR)
---
 .../multiphaseSystem/multiphaseSystem.C       | 67 ++++++++++++-------
 .../reactingMultiphaseEulerFoam/pU/pEqn.H     | 17 +++--
 2 files changed, 55 insertions(+), 29 deletions(-)

diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
index 62729c9b0..63f679db6 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -202,32 +202,38 @@ void Foam::multiphaseSystem::solveAlphas()
     }
     MULES::limitSum(alphafs, alphaPhiCorrs, fixedAlphaPhiCorrs);
 
-    // Solve for the moving phase alphas
+    // Compute alpha equation source terms for the moving phases
+    PtrList<volScalarField::Internal> Sps(movingPhases().size());
+    PtrList<volScalarField::Internal> Sus(movingPhases().size());
     forAll(movingPhases(), movingPhasei)
     {
         phaseModel& phase = movingPhases()[movingPhasei];
         volScalarField& alpha = phase;
 
-        surfaceScalarField& alphaPhi = alphaPhiCorrs[phase.index()];
-        alphaPhi += upwind<scalar>(mesh_, phi_).flux(phase);
-        phase.correctInflowOutflow(alphaPhi);
-
-        volScalarField::Internal Sp
+        Sps.set
         (
-            IOobject
+            phase.index(),
+            new volScalarField
             (
-                "Sp",
-                mesh_.time().timeName(),
-                mesh_
-            ),
-            mesh_,
-            dimensionedScalar(dimless/dimTime, 0)
+                IOobject
+                (
+                    "Sp",
+                    mesh_.time().timeName(),
+                    mesh_
+                ),
+                mesh_,
+                dimensionedScalar(dimless/dimTime, 0)
+            )
         );
 
-        volScalarField::Internal Su
+        Sus.set
         (
-            "Su",
-            min(alpha, scalar(1))*fvc::div(fvc::absolute(phi_, phase.U()))
+            phase.index(),
+            new volScalarField
+            (
+                "Su",
+                min(alpha, scalar(1))*fvc::div(fvc::absolute(phi_, phase.U()))
+            )
         );
 
         if (phase.divU().valid())
@@ -238,12 +244,12 @@ void Foam::multiphaseSystem::solveAlphas()
             {
                 if (dgdt[celli] > 0.0)
                 {
-                    Sp[celli] -= dgdt[celli];
-                    Su[celli] += dgdt[celli];
+                    Sps[phase.index()][celli] -= dgdt[celli];
+                    Sus[phase.index()][celli] += dgdt[celli];
                 }
                 else if (dgdt[celli] < 0.0)
                 {
-                    Sp[celli] +=
+                    Sps[phase.index()][celli] +=
                         dgdt[celli]
                        *(1 - alpha[celli])/max(alpha[celli], 1e-4);
                 }
@@ -265,29 +271,40 @@ void Foam::multiphaseSystem::solveAlphas()
                 {
                     if (dgdt2[celli] < 0.0)
                     {
-                        Sp[celli] +=
+                        Sps[phase.index()][celli] +=
                             dgdt2[celli]
                            *(1 - alpha2[celli])/max(alpha2[celli], 1e-4);
 
-                        Su[celli] -=
+                        Sus[phase.index()][celli] -=
                             dgdt2[celli]
                            *alpha[celli]/max(alpha2[celli], 1e-4);
                     }
                     else if (dgdt2[celli] > 0.0)
                     {
-                        Sp[celli] -= dgdt2[celli];
+                        Sps[phase.index()][celli] -= dgdt2[celli];
                     }
                 }
             }
         }
+    }
+
+    // Solve for the moving phase alphas
+    forAll(movingPhases(), movingPhasei)
+    {
+        phaseModel& phase = movingPhases()[movingPhasei];
+        volScalarField& alpha = phase;
+
+        surfaceScalarField& alphaPhi = alphaPhiCorrs[phase.index()];
+        alphaPhi += upwind<scalar>(mesh_, phi_).flux(phase);
+        phase.correctInflowOutflow(alphaPhi);
 
         MULES::explicitSolve
         (
             geometricOneField(),
             alpha,
             alphaPhi,
-            Sp,
-            Su
+            Sps[phase.index()],
+            Sus[phase.index()]
         );
 
         phase.alphaPhiRef() = alphaPhi;
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
index f50de4b7b..a92bb2811 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
@@ -361,13 +361,22 @@ while (pimple.correct())
 
         {
             fvScalarMatrix pEqn(pEqnIncomp);
-
-            forAll(phases, phasei)
             {
-                if (pEqnComps.set(phasei))
+                fvScalarMatrix pEqnComp
+                (
+                    p_rgh,
+                    pEqn.dimensions()
+                );
+
+                forAll(phases, phasei)
                 {
-                    pEqn += pEqnComps[phasei];
+                    if (pEqnComps.set(phasei))
+                    {
+                        pEqnComp += pEqnComps[phasei];
+                    }
                 }
+
+                pEqn += pEqnComp;
             }
 
             pEqn.solve();
-- 
2.20.1

