diff --git a/applications/solvers/combustion/reactingFoam/UEqn.H b/applications/solvers/combustion/reactingFoam/UEqn.H
index f493177d7616e5cfda0e0aa559d57fe5484aa099..1704ce4ab984743d03dc229c6d9ec10c4996c3d3 100644
--- a/applications/solvers/combustion/reactingFoam/UEqn.H
+++ b/applications/solvers/combustion/reactingFoam/UEqn.H
@@ -1,23 +1,25 @@
- MRF.correctBoundaryVelocity(U);
+// Solve the Momentum equation
- fvVectorMatrix UEqn
- (
- fvm::ddt(rho, U) + fvm::div(phi, U)
- + MRF.DDt(rho, U)
- + turbulence->divDevRhoReff(U)
- ==
- rho*g
- + fvOptions(rho, U)
- );
+MRF.correctBoundaryVelocity(U);
- UEqn.relax();
+tmp<fvVectorMatrix> UEqn
+(
+ fvm::ddt(rho, U) + fvm::div(phi, U)
+ + MRF.DDt(rho, U)
+ + turbulence->divDevRhoReff(U)
+ ==
+ rho*g
+ + fvOptions(rho, U)
+);
- fvOptions.constrain(UEqn);
+UEqn().relax();
- if (pimple.momentumPredictor())
- {
- solve(UEqn == -fvc::grad(p));
+fvOptions.constrain(UEqn());
- fvOptions.correct(U);
- K = 0.5*magSqr(U);
- }
+if (pimple.momentumPredictor())
+{
+ solve(UEqn() == -fvc::grad(p));
+
+ fvOptions.correct(U);
+ K = 0.5*magSqr(U);
+}
diff --git a/applications/solvers/combustion/reactingFoam/createFields.H b/applications/solvers/combustion/reactingFoam/createFields.H
index 19c115a8f2dfdceb4b749cefb062c7f7bacea427..71c805c43b83b21c660cd6e3cc11526c6e1c6a9b 100644
--- a/applications/solvers/combustion/reactingFoam/createFields.H
+++ b/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;
diff --git a/applications/solvers/combustion/reactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/pEqn.H
index cedab05eab6f72544e8b628a28715986da030a19..2f352db5da2523e04eedc4e48ff555c8a665dba6 100644
--- a/applications/solvers/combustion/reactingFoam/pEqn.H
+++ b/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);
diff --git a/applications/solvers/combustion/reactingFoam/pcEqn.H b/applications/solvers/combustion/reactingFoam/pcEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..3d7fffc459bdafe18405cf490478ec33578461fc
--- /dev/null
+++ b/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;
diff --git a/applications/solvers/combustion/reactingFoam/reactingFoam.C b/applications/solvers/combustion/reactingFoam/reactingFoam.C
index 79deedb763fbefe580b7c0c54b6a65bbf9528861..5ad452b52ba0b5fc3273f5c34e56e8790feaa511 100644
--- a/applications/solvers/combustion/reactingFoam/reactingFoam.C
+++ b/applications/solvers/combustion/reactingFoam/reactingFoam.C
@@ -95,7 +95,14 @@ int main(int argc, char *argv[])
// --- Pressure corrector loop
while (pimple.correct())
{
- #include "pEqn.H"
+ if (pimple.consistent())
+ {
+ #include "pcEqn.H"
+ }
+ else
+ {
+ #include "pEqn.H"
+ }
}
if (pimple.turbCorr())
diff --git a/applications/solvers/combustion/reactingFoam/setRDeltaT.H b/applications/solvers/combustion/reactingFoam/setRDeltaT.H
index a01f9bac3bec2726972c8aeacabf5898b7cd2aa3..b84be5a66d9356296b446d99dd902e1b73b04cce 100644
--- a/applications/solvers/combustion/reactingFoam/setRDeltaT.H
+++ b/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;