pimpleFoam family: Added PIMPLEC (PIMPLE-consistent) support

Select using the optional

    consistent true;

setting in the PIMPLE dictionary of fvSolution.

This option is generally only beneficial for cases run in PIMPLE-mode
with a large maximum Courant number.

applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H

@@ -2,11 +2,6 @@ volScalarField rAUrel(1.0/UrelEqn().A());
 volVectorField HbyA("HbyA", Urel);
 HbyA = rAUrel*UrelEqn().H();
 
-if (pimple.nCorrPISO() <= 1)
-{
-    UrelEqn.clear();
-}
-
 surfaceScalarField phiHbyA
 (
     "phiHbyA",
@@ -16,13 +11,28 @@ surfaceScalarField phiHbyA
 
 adjustPhi(phiHbyA, Urel, p);
 
+tmp<volScalarField> rAtUrel(rAUrel);
+
+if (pimple.consistent())
+{
+    rAtUrel = 1.0/max(1.0/rAUrel - UrelEqn().H1(), 0.1/rAUrel);
+    phiHbyA +=
+        fvc::interpolate(rAtUrel() - rAUrel)*fvc::snGrad(p)*mesh.magSf();
+    HbyA -= (rAUrel - rAtUrel())*fvc::grad(p);
+}
+
+if (pimple.nCorrPISO() <= 1)
+{
+    UrelEqn.clear();
+}
+
 // Non-orthogonal pressure corrector loop
 while (pimple.correctNonOrthogonal())
 {
     // Pressure corrector
     fvScalarMatrix pEqn
     (
-        fvm::laplacian(rAUrel, p) == fvc::div(phiHbyA)
+        fvm::laplacian(rAtUrel(), p) == fvc::div(phiHbyA)
     );
 
     pEqn.setReference(pRefCell, pRefValue);
@@ -40,6 +50,6 @@ while (pimple.correctNonOrthogonal())
 p.relax();
 
 // Momentum corrector
-Urel = HbyA - rAUrel*fvc::grad(p);
+Urel = HbyA - rAtUrel()*fvc::grad(p);
 Urel.correctBoundaryConditions();
 fvOptions.correct(Urel);

applications/solvers/incompressible/pimpleFoam/pEqn.H

@@ -1,23 +1,33 @@
-surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
-
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn().H();
 
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
+    (fvc::interpolate(HbyA) & mesh.Sf())
+  + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
+);
+
+MRF.makeRelative(phiHbyA);
+
+adjustPhi(phiHbyA, U, p);
+
+tmp<volScalarField> rAtU(rAU);
+
+if (pimple.consistent())
+{
+    rAtU = 1.0/max(1.0/rAU - UEqn().H1(), 0.1/rAU);
+    phiHbyA +=
+        fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
+    HbyA -= (rAU - rAtU())*fvc::grad(p);
+}
+
 if (pimple.nCorrPISO() <= 1)
 {
     UEqn.clear();
 }
 
-surfaceScalarField phiHbyA
+surfaceScalarField rAUf("rAUf", fvc::interpolate(rAtU()));
-(
-    "phiHbyA",
-    (fvc::interpolate(HbyA) & mesh.Sf())
-  + rAUf*fvc::ddtCorr(U, phi)
-);
-
-MRF.makeRelative(phiHbyA);
-
-adjustPhi(phiHbyA, U, p);
 
 // Update the fixedFluxPressure BCs to ensure flux consistency
 setSnGrad<fixedFluxPressureFvPatchScalarField>
@@ -53,6 +63,6 @@ while (pimple.correctNonOrthogonal())
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
-U = HbyA - rAU*fvc::grad(p);
+U = HbyA - rAtU()*fvc::grad(p);
 U.correctBoundaryConditions();
 fvOptions.correct(U);

applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H

@@ -1,18 +1,11 @@
-surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
-
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn().H();
 
-if (pimple.nCorrPISO() <= 1)
-{
-    UEqn.clear();
-}
-
 surfaceScalarField phiHbyA
 (
     "phiHbyA",
     (fvc::interpolate(HbyA) & mesh.Sf())
-  + rAUf*fvc::ddtCorr(U, Uf)
+  + fvc::interpolate(rAU)*fvc::ddtCorr(U, Uf)
 );
 
 MRF.makeRelative(phiHbyA);
@@ -24,6 +17,23 @@ if (p.needReference())
     fvc::makeAbsolute(phiHbyA, U);
 }
 
+tmp<volScalarField> rAtU(rAU);
+
+if (pimple.consistent())
+{
+    rAtU = 1.0/max(1.0/rAU - UEqn().H1(), 0.1/rAU);
+    phiHbyA +=
+        fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
+    HbyA -= (rAU - rAtU())*fvc::grad(p);
+}
+
+if (pimple.nCorrPISO() <= 1)
+{
+    UEqn.clear();
+}
+
+surfaceScalarField rAUf("rAUf", fvc::interpolate(rAtU()));
+
 // Update the fixedFluxPressure BCs to ensure flux consistency
 setSnGrad<fixedFluxPressureFvPatchScalarField>
 (
@@ -57,7 +67,7 @@ while (pimple.correctNonOrthogonal())
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
-U = HbyA - rAU*fvc::grad(p);
+U = HbyA - rAtU*fvc::grad(p);
 U.correctBoundaryConditions();
 fvOptions.correct(U);