reactingFoam: Added support for PIMPLE-consistent and pressure relaxation

Pressure relaxation is useful with LTS to damp acoustic waves

applications/solvers/combustion/reactingFoam/UEqn.H

@@ -1,6 +1,8 @@
+// Solve the Momentum equation
+
 MRF.correctBoundaryVelocity(U);
 
-    fvVectorMatrix UEqn
+tmp<fvVectorMatrix> UEqn
 (
     fvm::ddt(rho, U) + fvm::div(phi, U)
   + MRF.DDt(rho, U)
@@ -10,13 +12,13 @@
   + fvOptions(rho, U)
 );
 
-    UEqn.relax();
+UEqn().relax();
 
-    fvOptions.constrain(UEqn);
+fvOptions.constrain(UEqn());
 
 if (pimple.momentumPredictor())
 {
-        solve(UEqn == -fvc::grad(p));
+    solve(UEqn() == -fvc::grad(p));
 
     fvOptions.correct(U);
     K = 0.5*magSqr(U);

applications/solvers/combustion/reactingFoam/createFields.H

@@ -45,6 +45,28 @@ const volScalarField& T = thermo.T();
 
 #include "compressibleCreatePhi.H"
 
+dimensionedScalar rhoMax
+(
+    dimensionedScalar::lookupOrDefault
+    (
+        "rhoMax",
+        pimple.dict(),
+        dimDensity,
+        GREAT
+    )
+);
+
+dimensionedScalar rhoMin
+(
+    dimensionedScalar::lookupOrDefault
+    (
+        "rhoMin",
+        pimple.dict(),
+        dimDensity,
+        0
+    )
+);
+
 mesh.setFluxRequired(p.name());
 
 Info << "Creating turbulence model.\n" << nl;

applications/solvers/combustion/reactingFoam/pEqn.H

@@ -1,10 +1,18 @@
 rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
 
-volScalarField rAU(1.0/UEqn.A());
+volScalarField rAU(1.0/UEqn().A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
 volVectorField HbyA("HbyA", U);
-HbyA = rAU*UEqn.H();
+HbyA = rAU*UEqn().H();
+
+if (pimple.nCorrPISO() <= 1)
+{
+    UEqn.clear();
+}
 
 if (pimple.transonic())
 {
@@ -26,7 +34,7 @@ if (pimple.transonic())
         (
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
-          - fvm::laplacian(rho*rAU, p)
+          - fvm::laplacian(rhorAUf, p)
          ==
             fvOptions(psi, p, rho.name())
         );
@@ -58,7 +66,7 @@ else
         (
             fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
-          - fvm::laplacian(rho*rAU, p)
+          - fvm::laplacian(rhorAUf, p)
          ==
             fvOptions(psi, p, rho.name())
         );
@@ -75,6 +83,17 @@ else
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
+// Explicitly relax pressure for momentum corrector
+p.relax();
+
+// Recalculate density from the relaxed pressure
+rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
+Info<< "rho max/min : " << max(rho).value()
+    << " " << min(rho).value() << endl;
+
 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 fvOptions.correct(U);

applications/solvers/combustion/reactingFoam/pcEqn.H

@@ -0,0 +1,125 @@
+rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
+
+volScalarField rAU(1.0/UEqn().A());
+volScalarField rAtU(1.0/(1.0/rAU - UEqn().H1()));
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*UEqn().H();
+
+if (pimple.nCorrPISO() <= 1)
+{
+    UEqn.clear();
+}
+
+if (pimple.transonic())
+{
+    surfaceScalarField phid
+    (
+        "phid",
+        fvc::interpolate(psi)
+       *(
+            (fvc::interpolate(HbyA) & mesh.Sf())
+          + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
+           /fvc::interpolate(rho)
+        )
+    );
+
+    MRF.makeRelative(fvc::interpolate(psi), phid);
+
+    surfaceScalarField phic
+    (
+        "phic",
+        fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf()
+    );
+
+    HbyA -= (rAU - rAtU)*fvc::grad(p);
+
+    volScalarField rhorAtU("rhorAtU", rho*rAtU);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvm::div(phid, p)
+          + fvc::div(phic)
+          - fvm::laplacian(rhorAtU, p)
+         ==
+            fvOptions(psi, p, rho.name())
+        );
+
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi == phic + pEqn.flux();
+        }
+    }
+}
+else
+{
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        (
+            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+          + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
+        )
+    );
+
+    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+    phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf();
+    HbyA -= (rAU - rAtU)*fvc::grad(p);
+
+    volScalarField rhorAtU("rhorAtU", rho*rAtU);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvc::div(phiHbyA)
+          - fvm::laplacian(rhorAtU, p)
+         ==
+            fvOptions(psi, p, rho.name())
+        );
+
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA + pEqn.flux();
+        }
+    }
+}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+// Explicitly relax pressure for momentum corrector
+p.relax();
+
+U = HbyA - rAtU*fvc::grad(p);
+U.correctBoundaryConditions();
+fvOptions.correct(U);
+K = 0.5*magSqr(U);
+
+if (thermo.dpdt())
+{
+    dpdt = fvc::ddt(p);
+}
+
+// Recalculate density from the relaxed pressure
+rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+
+if (!pimple.transonic())
+{
+    rho.relax();
+}
+
+Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;

applications/solvers/combustion/reactingFoam/reactingFoam.C

@@ -94,9 +94,16 @@ int main(int argc, char *argv[])
 
             // --- Pressure corrector loop
             while (pimple.correct())
+            {
+                if (pimple.consistent())
+                {
+                    #include "pcEqn.H"
+                }
+                else
                 {
                     #include "pEqn.H"
                 }
+            }
 
             if (pimple.turbCorr())
             {

applications/solvers/combustion/reactingFoam/setRDeltaT.H

@@ -116,6 +116,9 @@ License
         );
     }
 
+    // Update tho boundary values of the reciprocal time-step
+    rDeltaT.correctBoundaryConditions();
+
     Info<< "    Overall     = "
         << gMin(1/rDeltaT.internalField())
         << ", " << gMax(1/rDeltaT.internalField()) << endl;