rhoPimpleFoam, rhoSimpleFoam, buoyantPimpleFoam, buoyantSimpleFoam: Consistency improvements

Various small changes to make comparison between pimple and simple
variants of the single-phase compressible solvers easier

applications/solvers/compressible/rhoPimpleFoam/createFields.H

@@ -40,7 +40,13 @@ volVectorField U
 
 #include "compressibleCreatePhi.H"
 
-pressureControl pressureControl(p, rho, pimple.dict(), false);
+pressureControl pressureControl
+(
+    p,
+    rho,
+    pimple.dict(),
+    thermo.incompressible()
+);
 
 mesh.setFluxRequired(p.name());
 

applications/solvers/compressible/rhoSimpleFoam/createFields.H

@@ -38,7 +38,13 @@ volVectorField U
 
 #include "compressibleCreatePhi.H"
 
-pressureControl pressureControl(p, rho, simple.dict());
+pressureControl pressureControl
+(
+    p,
+    rho,
+    simple.dict(),
+    thermo.incompressible()
+);
 
 mesh.setFluxRequired(p.name());
 

applications/solvers/compressible/rhoSimpleFoam/pEqn.H

@@ -27,11 +27,16 @@ volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
 
 tUEqn.clear();
 
-bool closedVolume = false;
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
+    fvc::interpolate(rho)*fvc::flux(HbyA)
+);
 
-surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho)*fvc::flux(HbyA));
 MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
+bool closedVolume = false;
+
 // Update the pressure BCs to ensure flux consistency
 constrainPressure(p, rho, U, phiHbyA, rhorAAtUf, MRF);
 
@@ -127,7 +132,7 @@ pressureControl.limit(p);
 
 // For closed-volume cases adjust the pressure and density levels
 // to obey overall mass continuity
-if (closedVolume)
+if (closedVolume && !thermo.incompressible())
 {
     p += (initialMass - fvc::domainIntegrate(psi*p))
         /fvc::domainIntegrate(psi);

applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H

@@ -7,8 +7,8 @@ if (!mesh.steady() && !pimple.simpleRho())
 // pressure solution
 const volScalarField psip0(psi*p);
 
-volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+const volScalarField rAU("rAU", 1.0/UEqn.A());
+const surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
 
 if (pimple.nCorrPiso() <= 1)
@@ -25,7 +25,7 @@ surfaceScalarField phiHbyA
 
 MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
-bool adjustMass = false;
+bool adjustMass = pimple.transonic() ? false : adjustPhi(phiHbyA, U, p_rgh);
 
 surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
@@ -172,6 +172,3 @@ if (thermo.dpdt())
         dpdt -= fvc::div(fvc::meshPhi(rho, U), p);
     }
 }
-
-Info<< "Min/max rho:" << min(rho).value() << ' '
-    << max(rho).value() << endl;

applications/solvers/heatTransfer/buoyantSimpleFoam/Make/options

@@ -1,5 +1,6 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/cfdTools \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \

applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,6 +35,7 @@ Description
 #include "turbulentFluidThermoModel.H"
 #include "radiationModel.H"
 #include "simpleControl.H"
+#include "pressureControl.H"
 #include "fvOptions.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H

@@ -70,21 +70,18 @@ volScalarField p_rgh
 // Force p_rgh to be consistent with p
 p_rgh = p - rho*gh;
 
-label pRefCell = 0;
-scalar pRefValue = 0.0;
-setRefCell
+pressureControl pressureControl
 (
     p,
     p_rgh,
+    rho,
     simple.dict(),
-    pRefCell,
-    pRefValue
+    thermo.incompressible()
 );
 
 mesh.setFluxRequired(p_rgh.name());
 
 dimensionedScalar initialMass = fvc::domainIntegrate(rho);
-dimensionedScalar totalVolume = sum(mesh.V());
 
 #include "createMRF.H"
 #include "createRadiationModel.H"

applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H

@@ -1,70 +1,110 @@
+rho = thermo.rho();
+
+const volScalarField rAU("rAU", 1.0/UEqn.A());
+const surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
+
+tUEqn.clear();
+
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
+    fvc::interpolate(rho)*fvc::flux(HbyA)
+);
+
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+bool closedVolume = simple.transonic() ? false : adjustPhi(phiHbyA, U, p_rgh);
+
+surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+
+phiHbyA += phig;
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
+
+fvScalarMatrix p_rghEqn(p_rgh, dimMass/dimTime);
+
+if (simple.transonic())
 {
-    rho = thermo.rho();
-    rho.relax();
-
-    volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
-    tUEqn.clear();
-
-    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
-
-    surfaceScalarField phiHbyA
+    surfaceScalarField phid
     (
-        "phiHbyA",
-        fvc::interpolate(rho)*fvc::flux(HbyA)
+        "phid",
+        (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
     );
 
-    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
-    bool closedVolume = adjustPhi(phiHbyA, U, p_rgh);
-
-    phiHbyA += phig;
-
-    // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
+    phiHbyA -= fvc::interpolate(psi*p_rgh)*phiHbyA/fvc::interpolate(rho);
 
     while (simple.correctNonOrthogonal())
     {
-        fvScalarMatrix p_rghEqn
+        p_rghEqn =
+            fvc::div(phiHbyA) + fvm::div(phid, p_rgh)
+          - fvm::laplacian(rhorAUf, p_rgh)
+         ==
+            fvOptions(psi, p_rgh, rho.name());
+
+        // Relax the pressure equation to ensure diagonal-dominance
+        p_rghEqn.relax();
+
+        p_rghEqn.setReference
         (
-            fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phiHbyA)
+            pressureControl.refCell(),
+            pressureControl.refValue()
         );
 
-        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
         p_rghEqn.solve();
-
-        if (simple.finalNonOrthogonalIter())
-        {
-            // Calculate the conservative fluxes
-            phi = phiHbyA - p_rghEqn.flux();
-
-            // Explicitly relax pressure for momentum corrector
-            p_rgh.relax();
-
-            // Correct the momentum source with the pressure gradient flux
-            // calculated from the relaxed pressure
-            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rhorAUf);
-            U.correctBoundaryConditions();
-            fvOptions.correct(U);
-        }
     }
-
-    #include "continuityErrs.H"
-
-    p = p_rgh + rho*gh;
-
-    // For closed-volume cases adjust the pressure level
-    // to obey overall mass continuity
-    if (!thermo.incompressible() && closedVolume)
-    {
-        p += (initialMass - fvc::domainIntegrate(psi*p))
-            /fvc::domainIntegrate(psi);
-        p_rgh = p - rho*gh;
-    }
-
-    rho = thermo.rho();
-    rho.relax();
-    Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value()
-        << endl;
+}
+else
+{
+    while (simple.correctNonOrthogonal())
+    {
+        p_rghEqn =
+            fvc::div(phiHbyA)
+          - fvm::laplacian(rhorAUf, p_rgh)
+          ==
+            fvOptions(psi, p_rgh, rho.name());
+
+        p_rghEqn.setReference
+        (
+            pressureControl.refCell(),
+            pressureControl.refValue()
+        );
+
+        p_rghEqn.solve();
+    }
+}
+
+phi = phiHbyA + p_rghEqn.flux();
+
+p = p_rgh + rho*gh;
+
+#include "incompressible/continuityErrs.H"
+
+// Explicitly relax pressure for momentum corrector
+p_rgh.relax();
+
+// Correct the momentum source with the pressure gradient flux
+// calculated from the relaxed pressure
+U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
+U.correctBoundaryConditions();
+fvOptions.correct(U);
+
+pressureControl.limit(p);
+
+// For closed-volume compressible cases adjust the pressure level
+// to obey overall mass continuity
+if (closedVolume && !thermo.incompressible())
+{
+    p += (initialMass - fvc::domainIntegrate(thermo.rho()))
+        /fvc::domainIntegrate(psi);
+    p_rgh = p - rho*gh;
+    p.correctBoundaryConditions();
+}
+
+rho = thermo.rho();
+
+if (!simple.transonic())
+{
+    rho.relax();
 }