From 60cffb2fc3121609f13a3bbeea4b8b0da9449701 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Sat, 27 Jun 2015 18:56:24 +0100
Subject: [PATCH] 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.
---
 .../pimpleFoam/SRFPimpleFoam/pEqn.H           | 24 +++++++++----
 .../solvers/incompressible/pimpleFoam/pEqn.H  | 36 ++++++++++++-------
 .../pimpleFoam/pimpleDyMFoam/pEqn.H           | 28 ++++++++++-----
 3 files changed, 59 insertions(+), 29 deletions(-)

diff --git a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H
index d5c38e72ed..d305d1a261 100644
--- a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H
+++ b/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);
diff --git a/applications/solvers/incompressible/pimpleFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/pEqn.H
index 0abd31d143..13118296c5 100644
--- a/applications/solvers/incompressible/pimpleFoam/pEqn.H
+++ b/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
-(
-    "phiHbyA",
-    (fvc::interpolate(HbyA) & mesh.Sf())
-  + rAUf*fvc::ddtCorr(U, phi)
-);
-
-MRF.makeRelative(phiHbyA);
-
-adjustPhi(phiHbyA, U, p);
+surfaceScalarField rAUf("rAUf", fvc::interpolate(rAtU()));
 
 // 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);
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
index c133a95727..5ad9b80020 100644
--- a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
+++ b/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);