reacting*EulerFoam: Added pressure referencing for incompressible phase systems

Patch contributed by Institute of Fluid Dynamics,
Helmholtz-Zentrum Dresden - Rossendorf (HZDR)

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C

@@ -407,6 +407,20 @@ Foam::PtrList<Foam::volScalarField> Foam::phaseSystem::dmdts() const
 }
 
 
+bool Foam::phaseSystem::incompressible() const
+{
+    forAll(phaseModels_, phasei)
+    {
+        if (!phaseModels_[phasei].incompressible())
+        {
+            return false;
+        }
+    }
+
+    return true;
+}
+
+
 bool Foam::phaseSystem::implicitPhasePressure(const phaseModel& phase) const
 {
     return false;

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -421,6 +421,9 @@ public:
             //- Return the mass transfer rates for each phase
             virtual PtrList<volScalarField> dmdts() const;
 
+            //- Return incompressibility
+            bool incompressible() const;
+
 
         // Transfers
 

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/createFields.H

@@ -38,6 +38,12 @@ volScalarField p_rgh
 
 label pRefCell = 0;
 scalar pRefValue = 0.0;
+if (fluid.incompressible())
+{
+    p = max(p_rgh + fluid.rho()*gh, pMin);
+
+    if (p_rgh.needReference())
+    {
         setRefCell
         (
             p,
@@ -46,6 +52,16 @@ setRefCell
             pRefCell,
             pRefValue
         );
+
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+        p_rgh = p - fluid.rho()*gh;
+    }
+}
 mesh.setFluxRequired(p_rgh.name());
 
 PtrList<volScalarField> rAUs;

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H

@@ -335,6 +335,11 @@ while (pimple.correct())
                 pEqn += pEqnComps[phasei];
             }
 
+            if (fluid.incompressible())
+            {
+                pEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+            }
+
             pEqn.solve();
         }
 
@@ -393,6 +398,17 @@ while (pimple.correct())
     // Update and limit the static pressure
     p = max(p_rgh + rho*gh, pMin);
 
+    // Account for static pressure reference
+    if (p_rgh.needReference() && fluid.incompressible())
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+    }
+
     // Limit p_rgh
     p_rgh = p - rho*gh;
 

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/pEqn.H

@@ -324,6 +324,11 @@ while (pimple.correct())
                 pEqn += pEqnComps[phasei];
             }
 
+            if (fluid.incompressible())
+            {
+                pEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+            }
+
             pEqn.solve();
         }
 
@@ -370,6 +375,17 @@ while (pimple.correct())
     // Update and limit the static pressure
     p = max(p_rgh + rho*gh, pMin);
 
+    // Account for static pressure reference
+    if (p_rgh.needReference() && fluid.incompressible())
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+    }
+
     // Limit p_rgh
     p_rgh = p - rho*gh;
 

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/createFields.H

@@ -37,6 +37,12 @@ volScalarField p_rgh
 
 label pRefCell = 0;
 scalar pRefValue = 0.0;
+if (fluid.incompressible())
+{
+    p = max(p_rgh + fluid.rho()*gh, pMin);
+
+    if (p_rgh.needReference())
+    {
         setRefCell
         (
             p,
@@ -45,6 +51,16 @@ setRefCell
             pRefCell,
             pRefValue
         );
+
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+        p_rgh = p - fluid.rho()*gh;
+    }
+}
 mesh.setFluxRequired(p_rgh.name());
 
 PtrList<volScalarField> rAUs;

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H

@@ -286,6 +286,11 @@ while (pimple.correct())
           - fvm::laplacian(rAUf, p_rgh)
         );
 
+        if (fluid.incompressible())
+        {
+            pEqnIncomp.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+        }
+
         // Solve
         solve(pEqnIncomp + pEqnComp1 + pEqnComp2);
 
@@ -365,6 +370,17 @@ while (pimple.correct())
     // Update and limit the static pressure
     p = max(p_rgh + rho*gh, pMin);
 
+    // Account for static pressure reference
+    if (p_rgh.needReference() && fluid.incompressible())
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+    }
+
     // Limit p_rgh
     p_rgh = p - rho*gh;
 

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H

@@ -298,6 +298,11 @@ while (pimple.correct())
           - fvm::laplacian(rAUf, p_rgh)
         );
 
+        if (fluid.incompressible())
+        {
+            pEqnIncomp.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+        }
+
         // Solve
         solve(pEqnIncomp + pEqnComp1 + pEqnComp2);
 
@@ -347,6 +352,17 @@ while (pimple.correct())
     // Update and limit the static pressure
     p = max(p_rgh + rho*gh, pMin);
 
+    // Account for static pressure reference
+    if (p_rgh.needReference() && fluid.incompressible())
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+    }
+
     // Limit p_rgh
     p_rgh = p - rho*gh;
 

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/trickleBed/system/fvSolution

@@ -61,8 +61,6 @@ PIMPLE
     nOuterCorrectors 3;
     nCorrectors      1;
     nNonOrthogonalCorrectors 0;
-    pRefCell         0;
-    pRefValue        1e5;
     partialElimination true;
 }