Merge branch 'develop' of develop.openfoam.com:Development/OpenFOAM-plus into develop

applications/solvers/combustion/reactingFoam/createFields.H

@@ -52,27 +52,7 @@ volScalarField& p = thermo.p();
 
 #include "compressibleCreatePhi.H"
 
-dimensionedScalar rhoMax
+pressureControl pressureControl(p, rho, pimple.dict(), false);
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMax",
-        pimple.dict(),
-        dimDensity,
-        GREAT
-    )
-);
-
-dimensionedScalar rhoMin
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMin",
-        pimple.dict(),
-        dimDensity,
-        0
-    )
-);
 
 mesh.setFluxRequired(p.name());
 

applications/solvers/combustion/reactingFoam/pEqn.H

@@ -1,7 +1,4 @@
 rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
 
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
@@ -87,19 +84,17 @@ else
 // 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);
 K = 0.5*magSqr(U);
 
+if (pressureControl.limit(p))
+{
+    p.correctBoundaryConditions();
+    rho = thermo.rho();
+}
+
 if (thermo.dpdt())
 {
     dpdt = fvc::ddt(p);

applications/solvers/combustion/reactingFoam/pcEqn.H

@@ -1,7 +1,4 @@
 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()));
@@ -109,19 +106,13 @@ U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
+if (pressureControl.limit(p))
+{
+    p.correctBoundaryConditions();
+    rho = thermo.rho();
+}
+
 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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,6 +37,7 @@ Description
 #include "psiCombustionModel.H"
 #include "multivariateScheme.H"
 #include "pimpleControl.H"
+#include "pressureControl.H"
 #include "fvOptions.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
@@ -114,6 +115,8 @@ int main(int argc, char *argv[])
             }
         }
 
+        rho = thermo.rho();
+
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H

@@ -1,9 +1,8 @@
-{
 rho = thermo.rho();
 
 // Thermodynamic density needs to be updated by psi*d(p) after the
-    // pressure solution - done in 2 parts. Part 1:
+// pressure solution
-    thermo.rho() -= psi*p;
+const volScalarField psip0(psi*p);
 
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
@@ -63,8 +62,8 @@
 
 p = p_rgh + rho*gh;
 
-    // Second part of thermodynamic density update
+// Thermodynamic density update
-    thermo.rho() += psi*p;
+thermo.correctRho(psi*p - psip0);
 
 if (thermo.dpdt())
 {
@@ -73,4 +72,3 @@
 
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
-}

applications/solvers/combustion/reactingFoam/rhoReactingFoam/createFields.H

@@ -52,6 +52,8 @@ volScalarField& p = thermo.p();
 
 #include "compressibleCreatePhi.H"
 
+pressureControl pressureControl(p, rho, pimple.dict(), false);
+
 mesh.setFluxRequired(p.name());
 
 

applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H

@@ -1,109 +0,0 @@
-{
-    rho = thermo.rho();
-
-    // Thermodynamic density needs to be updated by psi*d(p) after the
-    // pressure solution - done in 2 parts. Part 1:
-    thermo.rho() -= psi*p;
-
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-
-    if (pimple.transonic())
-    {
-        surfaceScalarField phiHbyA
-        (
-            "phiHbyA",
-            (
-                fvc::flux(HbyA)
-              + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
-            )
-        );
-
-        MRF.makeRelative(phiHbyA);
-
-        surfaceScalarField phid("phid", fvc::interpolate(thermo.psi())*phiHbyA);
-
-        phiHbyA *= fvc::interpolate(rho);
-
-        fvScalarMatrix pDDtEqn
-        (
-            fvc::ddt(rho) + fvc::div(phiHbyA)
-          + correction(psi*fvm::ddt(p) + fvm::div(phid, p))
-        );
-
-        while (pimple.correctNonOrthogonal())
-        {
-            fvScalarMatrix pEqn
-            (
-                pDDtEqn
-              - fvm::laplacian(rhorAUf, p)
-             ==
-                fvOptions(psi, p, rho.name())
-            );
-
-            pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
-
-            if (pimple.finalNonOrthogonalIter())
-            {
-                phi = phiHbyA + pEqn.flux();
-            }
-        }
-    }
-    else
-    {
-        surfaceScalarField phiHbyA
-        (
-            "phiHbyA",
-            (
-                fvc::flux(rho*HbyA)
-              + rhorAUf*fvc::ddtCorr(rho, U, phi)
-            )
-        );
-
-        MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
-        // Update the pressure BCs to ensure flux consistency
-        constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
-
-        fvScalarMatrix pDDtEqn
-        (
-            fvc::ddt(rho) + psi*correction(fvm::ddt(p))
-          + fvc::div(phiHbyA)
-         ==
-            fvOptions(psi, p, rho.name())
-        );
-
-        while (pimple.correctNonOrthogonal())
-        {
-            fvScalarMatrix pEqn
-            (
-                pDDtEqn
-              - fvm::laplacian(rhorAUf, p)
-            );
-
-            pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
-
-            if (pimple.finalNonOrthogonalIter())
-            {
-                phi = phiHbyA + pEqn.flux();
-            }
-        }
-    }
-
-    // Second part of thermodynamic density update
-    thermo.rho() += psi*p;
-
-    #include "rhoEqn.H"
-    #include "compressibleContinuityErrs.H"
-
-    U = HbyA - rAU*fvc::grad(p);
-    U.correctBoundaryConditions();
-    fvOptions.correct(U);
-    K = 0.5*magSqr(U);
-
-    if (thermo.dpdt())
-    {
-        dpdt = fvc::ddt(p);
-    }
-}

applications/solvers/combustion/reactingFoam/rhoReactingFoam/rhoReactingFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -38,6 +38,7 @@ Description
 #include "turbulentFluidThermoModel.H"
 #include "multivariateScheme.H"
 #include "pimpleControl.H"
+#include "pressureControl.H"
 #include "fvOptions.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
@@ -100,7 +101,14 @@ int main(int argc, char *argv[])
             // --- Pressure corrector loop
             while (pimple.correct())
             {
-                #include "pEqn.H"
+                if (pimple.consistent())
+                {
+                    #include "../../../compressible/rhoPimpleFoam/pcEqn.H"
+                }
+                else
+                {
+                    #include "../../../compressible/rhoPimpleFoam/pEqn.H"
+                }
             }
 
             if (pimple.turbCorr())

applications/solvers/compressible/rhoPimpleFoam/createFields.H

@@ -42,6 +42,8 @@ volVectorField U
 
 pressureControl pressureControl(p, rho, pimple.dict(), false);
 
+mesh.setFluxRequired(p.name());
+
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence
 (
@@ -54,8 +56,6 @@ autoPtr<compressible::turbulenceModel> turbulence
     )
 );
 
-mesh.setFluxRequired(p.name());
-
 Info<< "Creating field dpdt\n" << endl;
 volScalarField dpdt
 (
@@ -73,3 +73,26 @@ Info<< "Creating field kinetic energy K\n" << endl;
 volScalarField K("K", 0.5*magSqr(U));
 
 #include "createMRF.H"
+
+
+dimensionedScalar rhoMax
+(
+    dimensionedScalar::lookupOrDefault
+    (
+        "rhoMax",
+        pimple.dict(),
+        dimDensity,
+        GREAT
+    )
+);
+
+dimensionedScalar rhoMin
+(
+    dimensionedScalar::lookupOrDefault
+    (
+        "rhoMin",
+        pimple.dict(),
+        dimDensity,
+        0
+    )
+);

applications/solvers/compressible/rhoPimpleFoam/pEqn.H

@@ -1,3 +1,12 @@
+if (!pimple.SIMPLErho())
+{
+    rho = thermo.rho();
+}
+
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
+
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@@ -10,7 +19,7 @@ if (pimple.nCorrPISO() <= 1)
 surfaceScalarField phiHbyA
 (
     "phiHbyA",
-    fvc::flux(rho*HbyA)
+    fvc::interpolate(rho)*fvc::flux(HbyA)
   + rhorAUf*fvc::ddtCorr(rho, U, phi)
 );
 
@@ -26,20 +35,21 @@ if (pimple.transonic())
         "phid",
         (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
     );
-    phiHbyA -= fvc::interpolate(p)*phid;
 
-    while (pimple.correctNonOrthogonal())
+    phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
-    {
+
-        fvScalarMatrix pEqn
+    fvScalarMatrix pDDtEqn
     (
-            fvm::ddt(psi, p)
+        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
-          + fvc::div(phiHbyA)
+      + fvc::div(phiHbyA) + fvm::div(phid, p)
-          + fvm::div(phid, p)
-          - fvm::laplacian(rhorAUf, p)
      ==
         fvOptions(psi, p, rho.name())
     );
 
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAUf, p));
+
         // Relax the pressure equation to ensure diagonal-dominance
         pEqn.relax();
 
@@ -53,17 +63,18 @@ if (pimple.transonic())
 }
 else
 {
-    while (pimple.correctNonOrthogonal())
+    fvScalarMatrix pDDtEqn
-    {
-        fvScalarMatrix pEqn
     (
-            fvm::ddt(psi, p)
+        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
       + fvc::div(phiHbyA)
-          - fvm::laplacian(rhorAUf, p)
      ==
         fvOptions(psi, p, rho.name())
     );
 
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAUf, p));
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
@@ -84,15 +95,14 @@ U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
-pressureControl.limit(p);
+if (pressureControl.limit(p))
-p.correctBoundaryConditions();
-rho = thermo.rho();
-
-if (!pimple.transonic())
 {
-    rho.relax();
+    p.correctBoundaryConditions();
 }
 
+thermo.correctRho(psi*p - psip0, rhoMin, rhoMax) ;
+rho = thermo.rho();
+
 if (thermo.dpdt())
 {
     dpdt = fvc::ddt(p);

applications/solvers/compressible/rhoPimpleFoam/pcEqn.H

@@ -1,3 +1,12 @@
+if (!pimple.SIMPLErho())
+{
+    rho = thermo.rho();
+}
+
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
+
 volScalarField rAU(1.0/UEqn.A());
 volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@@ -11,7 +20,7 @@ surfaceScalarField phiHbyA
 (
     "phiHbyA",
     (
-        fvc::flux(rho*HbyA)
+        fvc::interpolate(rho)*fvc::flux(HbyA)
       + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
     )
 );
@@ -33,22 +42,22 @@ if (pimple.transonic())
 
     phiHbyA +=
         fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf()
-      - fvc::interpolate(p)*phid;
+      - fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
 
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
-    while (pimple.correctNonOrthogonal())
+    fvScalarMatrix pDDtEqn
-    {
-        fvScalarMatrix pEqn
     (
-            fvm::ddt(psi, p)
+        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
-          + fvc::div(phiHbyA)
+      + fvc::div(phiHbyA) + fvm::div(phid, p)
-          + fvm::div(phid, p)
-          - fvm::laplacian(rhorAtU, p)
      ==
         fvOptions(psi, p, rho.name())
     );
 
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAtU, p));
+
         // Relax the pressure equation to ensure diagonal-dominance
         pEqn.relax();
 
@@ -65,17 +74,18 @@ else
     phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf();
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
-    while (pimple.correctNonOrthogonal())
+    fvScalarMatrix pDDtEqn
-    {
-        fvScalarMatrix pEqn
     (
-            fvm::ddt(psi, p)
+        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
       + fvc::div(phiHbyA)
-          - fvm::laplacian(rhorAtU, p)
      ==
         fvOptions(psi, p, rho.name())
     );
 
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAtU, p));
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
@@ -96,13 +106,20 @@ U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
-pressureControl.limit(p);
+if (pressureControl.limit(p))
-p.correctBoundaryConditions();
-rho = thermo.rho();
-
-if (!pimple.transonic())
 {
-    rho.relax();
+    p.correctBoundaryConditions();
+    thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
+    rho = thermo.rho();
+}
+else if (pimple.SIMPLErho())
+{
+    thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
+    rho = thermo.rho();
+}
+else
+{
+    thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
 }
 
 if (thermo.dpdt())

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H

@@ -1,3 +1,5 @@
+rho = thermo.rho();
+
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@@ -27,13 +29,14 @@ if (pimple.transonic())
         "phid",
         (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
     );
-    phiHbyA -= fvc::interpolate(p)*phid;
+
+    phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
-            fvm::ddt(psi, p)
+            fvc::ddt(rho) + psi*correction(fvm::ddt(p))
           + fvc::div(phiHbyA)
           + fvm::div(phid, p)
           - fvm::laplacian(rhorAUf, p)
@@ -59,7 +62,7 @@ else
         // Pressure corrector
         fvScalarMatrix pEqn
         (
-            fvm::ddt(psi, p)
+            fvc::ddt(rho) + psi*correction(fvm::ddt(p))
           + fvc::div(phiHbyA)
           - fvm::laplacian(rhorAUf, p)
          ==
@@ -86,13 +89,10 @@ U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
-pressureControl.limit(p);
+if (pressureControl.limit(p))
+{
     p.correctBoundaryConditions();
     rho = thermo.rho();
-
-if (!pimple.transonic())
-{
-    rho.relax();
 }
 
 {

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C

@@ -122,6 +122,8 @@ int main(int argc, char *argv[])
             }
         }
 
+        rho = thermo.rho();
+
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/solvers/compressible/rhoSimpleFoam/pEqn.H

@@ -1,4 +1,3 @@
-{
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@@ -6,7 +5,7 @@
 
 bool closedVolume = false;
 
-    surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
+surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho)*fvc::flux(HbyA));
 MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
 // Update the pressure BCs to ensure flux consistency
@@ -19,7 +18,8 @@
         "phid",
         (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
     );
-        phiHbyA -= fvc::interpolate(p)*phid;
+
+    phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
 
     while (simple.correctNonOrthogonal())
     {
@@ -78,7 +78,6 @@
     }
 }
 
-
 #include "incompressible/continuityErrs.H"
 
 // Explicitly relax pressure for momentum corrector
@@ -88,7 +87,7 @@
 U.correctBoundaryConditions();
 fvOptions.correct(U);
 
-    pressureControl.limit(p);
+bool pLimited = pressureControl.limit(p);
 
 // For closed-volume cases adjust the pressure and density levels
 // to obey overall mass continuity
@@ -98,7 +97,10 @@
         /fvc::domainIntegrate(psi);
 }
 
+if (pLimited || closedVolume)
+{
     p.correctBoundaryConditions();
+}
 
 rho = thermo.rho();
 
@@ -106,4 +108,3 @@
 {
     rho.relax();
 }
-}

applications/solvers/compressible/rhoSimpleFoam/pcEqn.H

@@ -1,3 +1,5 @@
+rho = thermo.rho();
+
 volScalarField rAU(1.0/UEqn.A());
 volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@@ -5,7 +7,7 @@ tUEqn.clear();
 
 bool closedVolume = false;
 
-surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
+surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho)*fvc::flux(HbyA));
 MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
 volScalarField rhorAtU("rhorAtU", rho*rAtU);
@@ -23,7 +25,7 @@ if (simple.transonic())
 
     phiHbyA +=
         fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf()
-      - fvc::interpolate(p)*phid;
+      - fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
 
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
@@ -98,7 +100,7 @@ U = HbyA - rAtU*fvc::grad(p);
 U.correctBoundaryConditions();
 fvOptions.correct(U);
 
-pressureControl.limit(p);
+bool pLimited = pressureControl.limit(p);
 
 // For closed-volume cases adjust the pressure and density levels
 // to obey overall mass continuity
@@ -108,9 +110,11 @@ if (closedVolume)
         /fvc::domainIntegrate(psi);
 }
 
+if (pLimited || closedVolume)
+{
     p.correctBoundaryConditions();
+}
 
-// Recalculate density from the relaxed pressure
 rho = thermo.rho();
 
 if (!simple.transonic())

applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H

@@ -1,5 +1,3 @@
-{
-    bool closedVolume = p_rgh.needReference();
 
 dimensionedScalar compressibility = fvc::domainIntegrate(psi);
 bool compressible = (compressibility.value() > SMALL);
@@ -7,8 +5,8 @@
 rho = thermo.rho();
 
 // Thermodynamic density needs to be updated by psi*d(p) after the
-    // pressure solution - done in 2 parts. Part 1:
+// pressure solution
-    thermo.rho() -= psi*p_rgh;
+const volScalarField psip0(psi*p);
 
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
@@ -31,28 +29,26 @@
 // Update the pressure BCs to ensure flux consistency
 constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
 
-    tmp<fvScalarMatrix> p_rghDDtEqn
+fvScalarMatrix p_rghDDtEqn
-    (
-        new fvScalarMatrix(p_rgh, dimMass/dimTime)
-    );
-
-    if (compressible)
-    {
-        p_rghDDtEqn =
 (
     fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+  + fvc::div(phiHbyA)
+  ==
+    fvOptions(psi, p_rgh, rho.name())
 );
-    }
 
 while (pimple.correctNonOrthogonal())
 {
     fvScalarMatrix p_rghEqn
     (
-            p_rghDDtEqn()
+        p_rghDDtEqn
-          + fvc::div(phiHbyA)
       - fvm::laplacian(rhorAUf, p_rgh)
-        ==
+    );
-            fvOptions(psi, p_rgh, rho.name())
+
+    p_rghEqn.setReference
+    (
+        pRefCell,
+        compressible ? getRefCellValue(p_rgh, pRefCell) : pRefValue
     );
 
     p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
@@ -76,21 +72,10 @@
 
 p = p_rgh + rho*gh;
 
-    // Second part of thermodynamic density update
-    thermo.rho() += psi*p_rgh;
-
-    if (thermo.dpdt())
-    {
-        dpdt = fvc::ddt(p);
-    }
-
-    if (compressible)
-    {
 #include "rhoEqn.H"
-    }
 #include "compressibleContinuityErrs.H"
 
-    if (closedVolume)
+if (p_rgh.needReference())
 {
     if (!compressible)
     {
@@ -103,13 +88,21 @@
     }
     else
     {
-            p += (initialMass - fvc::domainIntegrate(thermo.rho()))
+        p += (initialMass - fvc::domainIntegrate(psi*p))
             /compressibility;
+        thermo.correctRho(psi*p - psip0);
         rho = thermo.rho();
-        }
         p_rgh = p - rho*gh;
     }
-
+}
-    Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value()
+else
-        << endl;
+{
+    thermo.correctRho(psi*p - psip0);
+}
+
+rho = thermo.rho();
+
+if (thermo.dpdt())
+{
+    dpdt = fvc::ddt(p);
 }

applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H

@@ -88,3 +88,26 @@ dimensionedScalar totalVolume = sum(mesh.V());
 
 #include "createMRF.H"
 #include "createRadiationModel.H"
+
+dimensionedScalar rhoMax
+(
+    dimensionedScalar::lookupOrDefault
+    (
+        "rhoMax",
+        simple.dict(),
+        dimDensity,
+        GREAT
+    )
+);
+
+dimensionedScalar rhoMin
+(
+    dimensionedScalar::lookupOrDefault
+    (
+        "rhoMin",
+        simple.dict(),
+        dimDensity,
+        0
+    )
+);
+

applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H

@@ -1,7 +1,4 @@
 {
-    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));
@@ -24,6 +21,8 @@
     // Update the pressure BCs to ensure flux consistency
     constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
 
+    dimensionedScalar compressibility = fvc::domainIntegrate(psi);
+    bool compressible = (compressibility.value() > SMALL);
 
     while (simple.correctNonOrthogonal())
     {
@@ -32,7 +31,12 @@
             fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
-        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+        p_rghEqn.setReference
+        (
+            pRefCell,
+            compressible ? getRefCellValue(p_rgh, pRefCell) : pRefValue
+        );
+
         p_rghEqn.solve();
 
         if (simple.finalNonOrthogonalIter())
@@ -51,12 +55,9 @@
         }
     }
 
-    #include "continuityErrs.H"
-
     p = p_rgh + rho*gh;
 
-    dimensionedScalar compressibility = fvc::domainIntegrate(psi);
+    #include "continuityErrs.H"
-    bool compressible = (compressibility.value() > SMALL);
 
     // For closed-volume cases adjust the pressure level
     // to obey overall mass continuity
@@ -75,9 +76,9 @@
         {
             p += (initialMass - fvc::domainIntegrate(psi*p))
                 /fvc::domainIntegrate(psi);
-        }
             p_rgh = p - rho*gh;
         }
+    }
 
     rho = thermo.rho();
     rho.relax();

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H

@@ -1,8 +1,10 @@
 {
+    /*
     rho = thermo.rho();
     rho = max(rho, rhoMin[i]);
     rho = min(rho, rhoMax[i]);
     rho.relax();
+    */
 
     volScalarField rAU("rAU", 1.0/UEqn.A());
     surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
@@ -67,12 +69,24 @@
 
     // For closed-volume cases adjust the pressure level
     // to obey overall mass continuity
-    if (closedVolume && compressible)
+    if (closedVolume)
     {
-        p += (initialMass - fvc::domainIntegrate(thermo.rho()))
+        if (!compressible)
+        {
+            p += dimensionedScalar
+            (
+                "p",
+                p.dimensions(),
+                pRefValue - getRefCellValue(p, pRefCell)
+            );
+        }
+        else
+        {
+            p += (initialMass - fvc::domainIntegrate(psi*p))
                 /compressibility;
             p_rgh = p - rho*gh;
         }
+    }
 
     rho = thermo.rho();
     rho = max(rho, rhoMin[i]);

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H

@@ -19,6 +19,9 @@ List<bool> frozenFlowFluid(fluidRegions.size(), false);
 PtrList<IOMRFZoneList> MRFfluid(fluidRegions.size());
 PtrList<fv::options> fluidFvOptions(fluidRegions.size());
 
+List<label> refCellFluid(fluidRegions.size());
+List<scalar> refValueFluid(fluidRegions.size());
+
 // Populate fluid field pointer lists
 forAll(fluidRegions, i)
 {
@@ -248,4 +251,20 @@ forAll(fluidRegions, i)
     );
 
     turbulence[i].validate();
+
+    refCellFluid[i] = 0;
+    refValueFluid[i] = 0.0;
+
+    if (p_rghFluid[i].needReference())
+    {
+        setRefCell
+        (
+            thermoFluid[i].p(),
+            p_rghFluid[i],
+            pimpleDict,
+            refCellFluid[i],
+            refValueFluid[i]
+        );
+    }
+
 }

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H

@@ -1,10 +1,13 @@
-{
 bool closedVolume = p_rgh.needReference();
 dimensionedScalar compressibility = fvc::domainIntegrate(psi);
 bool compressible = (compressibility.value() > SMALL);
 
 rho = thermo.rho();
 
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
+
 volScalarField rAU("rAU", 1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
@@ -26,35 +29,27 @@
 // Update the pressure BCs to ensure flux consistency
 constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
 
-    tmp<fvScalarMatrix> p_rghDDtEqn
-    (
-        new fvScalarMatrix(p_rgh, dimMass/dimTime)
-    );
-
 {
-        if (compressible)
+    fvScalarMatrix p_rghDDtEqn
-        {
-            p_rghDDtEqn =
     (
         fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
-             ==
+      + fvc::div(phiHbyA)
-                fvOptions(psi, p_rgh, rho.name())
     );
-        }
-
-        // Thermodynamic density needs to be updated by psi*d(p) after the
-        // pressure solution - done in 2 parts. Part 1:
-        thermo.rho() -= psi*p_rgh;
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
         fvScalarMatrix p_rghEqn
         (
-                p_rghDDtEqn()
+            p_rghDDtEqn
-              + fvc::div(phiHbyA)
           - fvm::laplacian(rhorAUf, p_rgh)
         );
 
+        p_rghEqn.setReference
+        (
+            pRefCell,
+            compressible ? getRefCellValue(p_rgh, pRefCell) : pRefValue
+        );
+
         p_rghEqn.solve
         (
             mesh.solver
@@ -73,6 +68,9 @@
         if (nonOrth == nNonOrthCorr)
         {
             phi = phiHbyA + p_rghEqn.flux();
+
+            p_rgh.relax();
+
             U = HbyA
               + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
             U.correctBoundaryConditions();
@@ -81,34 +79,51 @@
         }
     }
 
-        // Second part of thermodynamic density update
-        thermo.rho() += psi*p_rgh;
-    }
-
     p = p_rgh + rho*gh;
 
-    // Update pressure time derivative if needed
+    // Thermodynamic density update
-    if (thermo.dpdt())
+    //thermo.correctRho(psi*p - psip0);
-    {
-        dpdt = fvc::ddt(p);
 }
 
-    if (compressible)
+
-    {
+
 // Solve continuity
 #include "rhoEqn.H"
-    }
 
 // Update continuity errors
 #include "compressibleContinuityErrors.H"
 
 // For closed-volume cases adjust the pressure and density levels
 // to obey overall mass continuity
-    if (closedVolume && compressible)
+if (closedVolume)
 {
-        p += (initialMass - fvc::domainIntegrate(thermo.rho()))
+    if (!compressible)
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+    }
+    else
+    {
+        p += (initialMass - fvc::domainIntegrate(psi*p))
             /compressibility;
+        thermo.correctRho(psi*p - psip0);
         rho = thermo.rho();
         p_rgh = p - rho*gh;
     }
 }
+else
+{
+    thermo.correctRho(psi*p - psip0);
+}
+
+rho = thermo.rho();
+
+// Update pressure time derivative if needed
+if (thermo.dpdt())
+{
+    dpdt = fvc::ddt(p);
+}

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H

@@ -33,3 +33,6 @@
     );
 
     const bool frozenFlow = frozenFlowFluid[i];
+
+    const label pRefCell = refCellFluid[i];
+    const scalar pRefValue = refValueFluid[i];

applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -42,6 +42,7 @@ Description
 #include "radiationModel.H"
 #include "SLGThermo.H"
 #include "pimpleControl.H"
+#include "pressureControl.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
 

applications/solvers/lagrangian/reactingParcelFoam/pEqn.H

@@ -1,9 +1,8 @@
-{
 rho = thermo.rho();
 
 // Thermodynamic density needs to be updated by psi*d(p) after the
-    // pressure solution - done in 2 parts. Part 1:
+// pressure solution
-    thermo.rho() -= psi*p;
+const volScalarField psip0(psi*p);
 
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
@@ -49,8 +48,8 @@
 
 p.relax();
 
-    // Second part of thermodynamic density update
+// Thermodynamic density update
-    thermo.rho() += psi*p;
+thermo.correctRho(psi*p - psip0);
 
 #include "rhoEqn.H" // NOTE: flux and time scales now inconsistent
 #include "compressibleContinuityErrs.H"
@@ -70,4 +69,3 @@
 rho = min(rho, rhoMax);
 
 Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
-}

applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H

@@ -1,7 +1,6 @@
-{
 // Thermodynamic density needs to be updated by psi*d(p) after the
-    // pressure solution - done in 2 parts. Part 1:
+// pressure solution
-    thermo.rho() -= psi*p;
+const volScalarField psip0(psi*p);
 
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
@@ -39,8 +38,8 @@
 
 p.relax();
 
-    // Second part of thermodynamic density update
+// Thermodynamic density update
-    thermo.rho() += psi*p;
+thermo.correctRho(psi*p - psip0);
 
 #include "compressibleContinuityErrs.H"
 
@@ -54,4 +53,3 @@
 rho.relax();
 
 Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
-}

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H

@@ -104,8 +104,8 @@
     }
 
     // Update densities from change in p_rgh
-    rho1 += psi1*(p_rgh - p_rgh_0);
+    mixture.thermo1().correctRho(psi1*(p_rgh - p_rgh_0));
-    rho2 += psi2*(p_rgh - p_rgh_0);
+    mixture.thermo2().correctRho(psi2*(p_rgh - p_rgh_0));
 
     rho = alpha1*rho1 + alpha2*rho2;
 

applications/solvers/multiphase/compressibleInterFoam/createFields.H

@@ -38,8 +38,8 @@ volScalarField& alpha2(mixture.alpha2());
 
 Info<< "Reading thermophysical properties\n" << endl;
 
-volScalarField& rho1 = mixture.thermo1().rho();
+const volScalarField& rho1 = mixture.thermo1().rho();
-volScalarField& rho2 = mixture.thermo2().rho();
+const volScalarField& rho2 = mixture.thermo2().rho();
 
 volScalarField rho
 (

applications/solvers/multiphase/compressibleInterFoam/pEqn.H

@@ -104,8 +104,8 @@
     }
 
     // Update densities from change in p_rgh
-    rho1 += psi1*(p_rgh - p_rgh_0);
+    mixture.thermo1().correctRho(psi1*(p_rgh - p_rgh_0));
-    rho2 += psi2*(p_rgh - p_rgh_0);
+    mixture.thermo2().correctRho(psi2*(p_rgh - p_rgh_0));
 
     rho = alpha1*rho1 + alpha2*rho2;
 

src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/fixedIncidentRadiation/fixedIncidentRadiationFvPatchScalarField.C

@@ -45,7 +45,7 @@ fixedIncidentRadiationFvPatchScalarField
 :
     fixedGradientFvPatchScalarField(p, iF),
     temperatureCoupledBase(patch(), "undefined", "undefined", "undefined-K"),
-    QrIncident_(p.size(), 0.0)
+    qrIncident_(p.size(), 0.0)
 {}
 
 
@@ -60,7 +60,7 @@ fixedIncidentRadiationFvPatchScalarField
 :
     fixedGradientFvPatchScalarField(psf, p, iF, mapper),
     temperatureCoupledBase(patch(), psf),
-    QrIncident_(psf.QrIncident_)
+    qrIncident_(psf.qrIncident_)
 {}
 
 
@@ -74,7 +74,7 @@ fixedIncidentRadiationFvPatchScalarField
 :
     fixedGradientFvPatchScalarField(p, iF),
     temperatureCoupledBase(patch(), dict),
-    QrIncident_("QrIncident", dict, p.size())
+    qrIncident_("qrIncident", dict, p.size())
 {
     if (dict.found("value") && dict.found("gradient"))
     {
@@ -99,7 +99,7 @@ fixedIncidentRadiationFvPatchScalarField
 :
     fixedGradientFvPatchScalarField(psf, iF),
     temperatureCoupledBase(patch(), psf),
-    QrIncident_(psf.QrIncident_)
+    qrIncident_(psf.qrIncident_)
 {}
 
 
@@ -111,7 +111,7 @@ fixedIncidentRadiationFvPatchScalarField
 :
     fixedGradientFvPatchScalarField(ptf),
     temperatureCoupledBase(patch(), ptf),
-    QrIncident_(ptf.QrIncident_)
+    qrIncident_(ptf.qrIncident_)
 {}
 
 
@@ -123,7 +123,7 @@ void Foam::radiation::fixedIncidentRadiationFvPatchScalarField::autoMap
 )
 {
     fixedGradientFvPatchScalarField::autoMap(m);
-    QrIncident_.autoMap(m);
+    qrIncident_.autoMap(m);
 }
 
 
@@ -141,7 +141,7 @@ void Foam::radiation::fixedIncidentRadiationFvPatchScalarField::rmap
             psf
         );
 
-    QrIncident_.rmap(thftpsf.QrIncident_, addr);
+    qrIncident_.rmap(thftpsf.qrIncident_, addr);
 }
 
 
@@ -165,7 +165,7 @@ void Foam::radiation::fixedIncidentRadiationFvPatchScalarField::updateCoeffs()
     gradient() =
         emissivity
        *(
-            QrIncident_
+            qrIncident_
           - physicoChemical::sigma.value()*pow4(*this)
         )/kappa(*this);
 
@@ -173,11 +173,11 @@ void Foam::radiation::fixedIncidentRadiationFvPatchScalarField::updateCoeffs()
 
     if (debug)
     {
-        scalar Qr = gSum(kappa(*this)*gradient()*patch().magSf());
+        scalar qr = gSum(kappa(*this)*gradient()*patch().magSf());
         Info<< patch().boundaryMesh().mesh().name() << ':'
             << patch().name() << ':'
             << this->internalField().name() << " -> "
-            << " radiativeFlux:" << Qr
+            << " radiativeFlux:" << qr
             << " walltemperature "
             << " min:" << gMin(*this)
             << " max:" << gMax(*this)
@@ -194,7 +194,7 @@ void Foam::radiation::fixedIncidentRadiationFvPatchScalarField::write
 {
     fixedGradientFvPatchScalarField::write(os);
     temperatureCoupledBase::write(os);
-    QrIncident_.writeEntry("QrIncident", os);
+    qrIncident_.writeEntry("qrIncident", os);
     writeEntry("value", os);
 }
 

src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/fixedIncidentRadiation/fixedIncidentRadiationFvPatchScalarField.H

@@ -34,19 +34,19 @@ Description
 
     the gradient heat flux is calculated as :
 
-        Qr = emissivity*(QrIncident - sigma_*T^4)
+        qr = emissivity*(qrIncident - sigma_*T^4)
 
     where:
 
     emissivity is the emissivity of the solid.
-    QrIncident is the specified fixed incident radiation.
+    qrIncident is the specified fixed incident radiation.
 
     Example usage:
 
     wall
     {
         type            fixedIncidentRadiation;
-        QrIncident      uniform 500;
+        qrIncident      uniform 500;
         kappa           solidThermo;
         KappaName       none;
     }
@@ -89,7 +89,7 @@ class fixedIncidentRadiationFvPatchScalarField
     // Private data
 
         //- Incident radiative heat flux
-        scalarField QrIncident_;
+        scalarField qrIncident_;
 
 
 public:

src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/lumpedMassWallTemperature/lumpedMassWallTemperatureFvPatchScalarField.C

@@ -150,18 +150,18 @@ void Foam::lumpedMassWallTemperatureFvPatchScalarField::updateCoeffs()
 
     if (debug)
     {
-        scalar qin(0);
+        scalar Qin(0);
-        scalar qout(0);
+        scalar Qout(0);
 
         forAll(q, facei)
         {
             if (q[facei] > 0.0) // out the wall
             {
-                qout += q[facei]*magSf[facei];
+                Qout += q[facei]*magSf[facei];
             }
             else if (q[facei] < 0.0) // into the wall
             {
-                qin += q[facei]*magSf[facei];
+                Qin += q[facei]*magSf[facei];
             }
         }
 
@@ -173,8 +173,8 @@ void Foam::lumpedMassWallTemperatureFvPatchScalarField::updateCoeffs()
             << " min:" << gMin(*this)
             << " max:" << gMax(*this)
             << " avg:" << gAverage(*this)
-            << " Qin [W]:" << qin
+            << " Qin [W]:" << Qin
-            << " Qout [W]:" << qout
+            << " Qout [W]:" << Qout
             << endl;
     }
 

src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/thermalBaffle1D/thermalBaffle1DFvPatchScalarField.C

@@ -50,7 +50,7 @@ thermalBaffle1DFvPatchScalarField
     TName_("T"),
     baffleActivated_(true),
     thickness_(p.size()),
-    Qs_(p.size()),
+    qs_(p.size()),
     solidDict_(),
     solidPtr_(nullptr),
     qrPrevious_(p.size()),
@@ -74,7 +74,7 @@ thermalBaffle1DFvPatchScalarField
     TName_(ptf.TName_),
     baffleActivated_(ptf.baffleActivated_),
     thickness_(ptf.thickness_, mapper),
-    Qs_(ptf.Qs_, mapper),
+    qs_(ptf.qs_, mapper),
     solidDict_(ptf.solidDict_),
     solidPtr_(ptf.solidPtr_),
     qrPrevious_(ptf.qrPrevious_, mapper),
@@ -97,7 +97,7 @@ thermalBaffle1DFvPatchScalarField
     TName_("T"),
     baffleActivated_(dict.lookupOrDefault<bool>("baffleActivated", true)),
     thickness_(),
-    Qs_(p.size(), 0),
+    qs_(p.size(), 0),
     solidDict_(dict),
     solidPtr_(),
     qrPrevious_(p.size(), 0.0),
@@ -111,9 +111,9 @@ thermalBaffle1DFvPatchScalarField
         thickness_ = scalarField("thickness", dict, p.size());
     }
 
-    if (dict.found("Qs"))
+    if (dict.found("qs"))
     {
-        Qs_ = scalarField("Qs", dict, p.size());
+        qs_ = scalarField("qs", dict, p.size());
     }
 
     if (dict.found("qrPrevious"))
@@ -151,7 +151,7 @@ thermalBaffle1DFvPatchScalarField
     TName_(ptf.TName_),
     baffleActivated_(ptf.baffleActivated_),
     thickness_(ptf.thickness_),
-    Qs_(ptf.Qs_),
+    qs_(ptf.qs_),
     solidDict_(ptf.solidDict_),
     solidPtr_(ptf.solidPtr_),
     qrPrevious_(ptf.qrPrevious_),
@@ -173,7 +173,7 @@ thermalBaffle1DFvPatchScalarField
     TName_(ptf.TName_),
     baffleActivated_(ptf.baffleActivated_),
     thickness_(ptf.thickness_),
-    Qs_(ptf.Qs_),
+    qs_(ptf.qs_),
     solidDict_(ptf.solidDict_),
     solidPtr_(ptf.solidPtr_),
     qrPrevious_(ptf.qrPrevious_),
@@ -264,11 +264,11 @@ baffleThickness() const
 
 
 template<class solidType>
-tmp<scalarField> thermalBaffle1DFvPatchScalarField<solidType>::Qs() const
+tmp<scalarField> thermalBaffle1DFvPatchScalarField<solidType>::qs() const
 {
     if (this->owner())
     {
-         return Qs_;
+         return qs_;
     }
     else
     {
@@ -283,10 +283,10 @@ tmp<scalarField> thermalBaffle1DFvPatchScalarField<solidType>::Qs() const
             nbrPatch.template lookupPatchField<volScalarField, scalar>(TName_)
         );
 
-        tmp<scalarField> tQs(new scalarField(nbrField.Qs()));
+        tmp<scalarField> tqs(new scalarField(nbrField.qs()));
-        scalarField& Qs = tQs.ref();
+        scalarField& qs = tqs.ref();
-        mapDist.distribute(Qs);
+        mapDist.distribute(qs);
-        return tQs;
+        return tqs;
     }
 }
 
@@ -304,7 +304,7 @@ void thermalBaffle1DFvPatchScalarField<solidType>::autoMap
     if (this->owner())
     {
         thickness_.autoMap(m);
-        Qs_.autoMap(m);
+        qs_.autoMap(m);
     }
 }
 
@@ -324,7 +324,7 @@ void thermalBaffle1DFvPatchScalarField<solidType>::rmap
     if (this->owner())
     {
         thickness_.rmap(tiptf.thickness_, addr);
-        Qs_.rmap(tiptf.Qs_, addr);
+        qs_.rmap(tiptf.qs_, addr);
     }
 }
 
@@ -397,7 +397,7 @@ void thermalBaffle1DFvPatchScalarField<solidType>::updateCoeffs()
 
         valueFraction() = alpha/(alpha + myKDelta);
 
-        refValue() = (KDeltaSolid*nbrTp + Qs()/2.0)/alpha;
+        refValue() = (KDeltaSolid*nbrTp + qs()/2.0)/alpha;
 
         if (debug)
         {
@@ -431,7 +431,7 @@ void thermalBaffle1DFvPatchScalarField<solidType>::write(Ostream& os) const
     if (this->owner())
     {
         baffleThickness()().writeEntry("thickness", os);
-        Qs()().writeEntry("Qs", os);
+        qs()().writeEntry("qs", os);
         solid().write(os);
     }
 

src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/thermalBaffle1D/thermalBaffle1DFvPatchScalarField.H

@@ -46,7 +46,7 @@ Usage
         samplePatch     <slavePatchName>;
 
         thickness       uniform 0.005;  // thickness [m]
-        Qs              uniform 100;    // heat flux [W/m2]
+        qs              uniform 100;    // heat flux [W/m2]
 
         qr              none;
         relaxation      1;
@@ -124,7 +124,7 @@ class thermalBaffle1DFvPatchScalarField
         mutable scalarField thickness_;
 
         //- Superficial heat source [W/m2]
-        mutable scalarField Qs_;
+        mutable scalarField qs_;
 
         //- Solid dictionary
         dictionary solidDict_;
@@ -147,8 +147,8 @@ class thermalBaffle1DFvPatchScalarField
         //- Return const solid thermo
         const solidType& solid() const;
 
-        //- Return Qs from master
+        //- Return qs from master
-        tmp<scalarField> Qs() const;
+        tmp<scalarField> qs() const;
 
         //- Return thickness from master
         tmp<scalarField> baffleThickness() const;

src/dynamicMesh/slidingInterface/slidingInterfaceProjectPoints.C

@@ -219,7 +219,7 @@ bool Foam::slidingInterface::projectPoints() const
             projectionAlgo_
         );
 
-//     Pout<< "USING N-SQAURED!!!" << endl;
+//     Pout<< "USING N-SQUARED!!!" << endl;
 //     List<objectHit> slavePointFaceHits =
 //         projectPointsNSquared<face, List, const pointField&>
 //         (

src/finiteVolume/cfdTools/general/pressureControl/pressureControl.C

@@ -38,24 +38,30 @@ Foam::pressureControl::pressureControl
 :
     refCell_(-1),
     refValue_(0),
-    pMax_("pMax", dimPressure, 0),
+    pMax_("pMax", dimPressure, GREAT),
-    pMin_("pMin", dimPressure, GREAT)
+    pMin_("pMin", dimPressure, 0),
+    limitMaxP_(false),
+    limitMinP_(false)
 {
     bool pLimits = false;
+    scalar pMax = -GREAT;
+    scalar pMin = GREAT;
 
     // Set the reference cell and value for closed domain simulations
     if (pRefRequired && setRefCell(p, dict, refCell_, refValue_))
     {
         pLimits = true;
 
-        pMax_.value() = refValue_;
+        pMax = refValue_;
-        pMin_.value() = refValue_;
+        pMin = refValue_;
     }
 
     if (dict.found("pMax") && dict.found("pMin"))
     {
         pMax_.value() = readScalar(dict.lookup("pMax"));
+        limitMaxP_ = true;
         pMin_.value() = readScalar(dict.lookup("pMin"));
+        limitMinP_ = true;
     }
     else
     {
@@ -73,8 +79,8 @@ Foam::pressureControl::pressureControl
                 pLimits = true;
                 rhoLimits = true;
 
-                pMax_.value() = max(pMax_.value(), max(pbf[patchi]));
+                pMax = max(pMax, max(pbf[patchi]));
-                pMin_.value() = min(pMin_.value(), min(pbf[patchi]));
+                pMin = min(pMin, min(pbf[patchi]));
 
                 rhoRefMax = max(rhoRefMax, max(rhobf[patchi]));
                 rhoRefMin = min(rhoRefMin, min(rhobf[patchi]));
@@ -84,8 +90,8 @@ Foam::pressureControl::pressureControl
         reduce(rhoLimits, andOp<bool>());
         if (rhoLimits)
         {
-            reduce(pMax_.value(), maxOp<scalar>());
+            reduce(pMax, maxOp<scalar>());
-            reduce(pMin_.value(), minOp<scalar>());
+            reduce(pMin, minOp<scalar>());
 
             reduce(rhoRefMax, maxOp<scalar>());
             reduce(rhoRefMin, minOp<scalar>());
@@ -94,6 +100,7 @@ Foam::pressureControl::pressureControl
         if (dict.found("pMax"))
         {
             pMax_.value() = readScalar(dict.lookup("pMax"));
+            limitMaxP_ = true;
         }
         else if (dict.found("pMaxFactor"))
         {
@@ -108,7 +115,8 @@ Foam::pressureControl::pressureControl
             }
 
             const scalar pMaxFactor(readScalar(dict.lookup("pMaxFactor")));
-            pMax_ *= pMaxFactor;
+            pMax_.value() = pMaxFactor*pMax;
+            limitMaxP_ = true;
         }
         else if (dict.found("rhoMax"))
         {
@@ -142,12 +150,14 @@ Foam::pressureControl::pressureControl
 
             dimensionedScalar rhoMax("rhoMax", dimDensity, dict);
 
-            pMax_ *= max(rhoMax.value()/rhoRefMax, 1);
+            pMax_.value() = max(rhoMax.value()/rhoRefMax, 1)*pMax;
+            limitMaxP_ = true;
         }
 
         if (dict.found("pMin"))
         {
             pMin_.value() = readScalar(dict.lookup("pMin"));
+            limitMinP_ = true;
         }
         else if (dict.found("pMinFactor"))
         {
@@ -162,7 +172,8 @@ Foam::pressureControl::pressureControl
             }
 
             const scalar pMinFactor(readScalar(dict.lookup("pMinFactor")));
-            pMin_ *= pMinFactor;
+            pMin_.value() = pMinFactor*pMin;
+            limitMinP_ = true;
         }
         else if (dict.found("rhoMin"))
         {
@@ -195,27 +206,65 @@ Foam::pressureControl::pressureControl
 
             dimensionedScalar rhoMin("rhoMin", dimDensity, dict);
 
-            pMin_ *= min(rhoMin.value()/rhoRefMin, 1);
+            pMin_.value() = min(rhoMin.value()/rhoRefMin, 1)*pMin;
+            limitMinP_ = true;
         }
     }
 
-    Info<< "pressureControl" << nl
+    if (limitMaxP_ || limitMinP_)
-        << "    pMax/pMin " << pMax_.value() << " " << pMin_.value()
+    {
-        << nl << endl;
+        Info<< "pressureControl" << nl;
+
+        if (limitMaxP_)
+        {
+            Info<< "    pMax " << pMax_.value() << nl;
+        }
+
+        if (limitMinP_)
+        {
+            Info<< "    pMin " << pMin_.value() << nl;
+        }
+
+        Info << endl;
+    }
 }
 
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-void Foam::pressureControl::limit(volScalarField& p) const
+bool Foam::pressureControl::limit(volScalarField& p) const
 {
-    Info<< "pressureControl: p max/min "
+    if (limitMaxP_ || limitMinP_)
-        << max(p).value() << " "
+    {
-        << min(p).value() << endl;
+        if (limitMaxP_)
+        {
+            const scalar pMax = max(p).value();
 
-    p = max(p, pMin_);
+            if (pMax > pMax_.value())
+            {
+                Info<< "pressureControl: p max " << pMax << endl;
                 p = min(p, pMax_);
             }
+        }
+
+        if (limitMinP_)
+        {
+            const scalar pMin = min(p).value();
+
+            if (pMin < pMin_.value())
+            {
+                Info<< "pressureControl: p min " << pMin << endl;
+                p = max(p, pMin_);
+            }
+        }
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
 
 
 // ************************************************************************* //

src/finiteVolume/cfdTools/general/pressureControl/pressureControl.H

@@ -66,6 +66,12 @@ class pressureControl
         //- Pressure upper-limit
         dimensionedScalar pMin_;
 
+        //- Pressure lower-limit
+        bool limitMaxP_;
+
+        //- Pressure upper-limit
+        bool limitMinP_;
+
 
 public:
 
@@ -89,8 +95,8 @@ public:
         //- Return the pressure reference level
         inline scalar refValue() const;
 
-        //- Limit the pressure
+        //- Limit the pressure if necessary and return true if so
-        void limit(volScalarField& p) const;
+        bool limit(volScalarField& p) const;
 };
 
 

src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -45,6 +45,7 @@ void Foam::pimpleControl::read()
     solveFlow_ = pimpleDict.lookupOrDefault<Switch>("solveFlow", true);
     nCorrPIMPLE_ = pimpleDict.lookupOrDefault<label>("nOuterCorrectors", 1);
     nCorrPISO_ = pimpleDict.lookupOrDefault<label>("nCorrectors", 1);
+    SIMPLErho_ = pimpleDict.lookupOrDefault<Switch>("SIMPLErho", false);
     turbOnFinalIterOnly_ =
         pimpleDict.lookupOrDefault<Switch>("turbOnFinalIterOnly", true);
 }
@@ -128,6 +129,7 @@ Foam::pimpleControl::pimpleControl(fvMesh& mesh, const word& dictName)
     nCorrPIMPLE_(0),
     nCorrPISO_(0),
     corrPISO_(0),
+    SIMPLErho_(false),
     turbOnFinalIterOnly_(true),
     converged_(false)
 {

src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -81,6 +81,10 @@ protected:
             //- Current PISO corrector
             label corrPISO_;
 
+            //- Flag to indicate whether to update density in SIMPLE
+            //  rather than PISO mode
+            bool SIMPLErho_;
+
             //- Flag to indicate whether to only solve turbulence on final iter
             bool turbOnFinalIterOnly_;
 
@@ -128,6 +132,10 @@ public:
             //- Current PISO corrector index
             inline label corrPISO() const;
 
+            //- Flag to indicate whether to update density in SIMPLE
+            //  rather than PISO mode
+            inline bool SIMPLErho() const;
+
 
         // Solution control
 

src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,6 +43,12 @@ inline Foam::label Foam::pimpleControl::corrPISO() const
 }
 
 
+inline bool Foam::pimpleControl::SIMPLErho() const
+{
+    return SIMPLErho_;
+}
+
+
 inline bool Foam::pimpleControl::correct()
 {
     corrPISO_++;

src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrad.C

@@ -92,6 +92,7 @@ Foam::fv::correctedSnGrad<Type>::correction
         )
     );
     GeometricField<Type, fvsPatchField, surfaceMesh>& ssf = tssf.ref();
+    ssf.setOriented();
 
     for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
     {

src/fvOptions/sources/derived/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.H

@@ -48,7 +48,7 @@ Description
     \endvartable
 
 
-    The distribution inside the hear exchanger is given by:
+    The distribution inside the heat exchanger is given by:
     \f[
         Q_c = \frac{V_c |U_c| (T_c - T_{ref})}{\sum(V_c |U_c| (T_c - T_{ref}))}
     \f]

src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/primaryRadiation/primaryRadiation.C

@@ -56,11 +56,11 @@ primaryRadiation::primaryRadiation
 )
 :
     filmRadiationModel(typeName, film, dict),
-    QinPrimary_
+    qinPrimary_
     (
         IOobject
         (
-            "Qin", // same name as Qin on primary region to enable mapping
+            "qin", // same name as qin on primary region to enable mapping
             film.time().timeName(),
             film.regionMesh(),
             IOobject::NO_READ,
@@ -83,8 +83,8 @@ primaryRadiation::~primaryRadiation()
 
 void primaryRadiation::correct()
 {
-    // Transfer Qin from primary region
+    // Transfer qin from primary region
-    QinPrimary_.correctBoundaryConditions();
+    qinPrimary_.correctBoundaryConditions();
 }
 
 
@@ -108,10 +108,10 @@ tmp<volScalarField> primaryRadiation::Shs()
     );
 
     scalarField& Shs = tShs.ref();
-    const scalarField& QinP = QinPrimary_;
+    const scalarField& qinP = qinPrimary_;
     const scalarField& alpha = filmModel_.alpha();
 
-    Shs = QinP*alpha;
+    Shs = qinP*alpha;
 
     return tShs;
 }

src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/primaryRadiation/primaryRadiation.H

@@ -61,7 +61,7 @@ private:
     // Private data
 
         //- Incident radiative flux mapped from  the primary region / [kg/s3]
-        volScalarField QinPrimary_;
+        volScalarField qinPrimary_;
 
 
     // Private member functions

src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.C

@@ -57,11 +57,11 @@ standardRadiation::standardRadiation
 )
 :
     filmRadiationModel(typeName, film, dict),
-    QinPrimary_
+    qinPrimary_
     (
         IOobject
         (
-            "Qin", // same name as Qin on primary region to enable mapping
+            "qin", // same name as qin on primary region to enable mapping
             film.time().timeName(),
             film.regionMesh(),
             IOobject::NO_READ,
@@ -101,7 +101,7 @@ standardRadiation::~standardRadiation()
 void standardRadiation::correct()
 {
     // Transfer qr from primary region
-    QinPrimary_.correctBoundaryConditions();
+    qinPrimary_.correctBoundaryConditions();
 }
 
 
@@ -125,14 +125,14 @@ tmp<volScalarField> standardRadiation::Shs()
     );
 
     scalarField& Shs = tShs.ref();
-    const scalarField& QinP = QinPrimary_;
+    const scalarField& qinP = qinPrimary_;
     const scalarField& delta = filmModel_.delta();
     const scalarField& alpha = filmModel_.alpha();
 
-    Shs = beta_*QinP*alpha*(1.0 - exp(-kappaBar_*delta));
+    Shs = beta_*qinP*alpha*(1.0 - exp(-kappaBar_*delta));
 
     // Update net qr on local region
-    qrNet_.primitiveFieldRef() = QinP - Shs;
+    qrNet_.primitiveFieldRef() = qinP - Shs;
     qrNet_.correctBoundaryConditions();
 
     return tShs;

src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.H

@@ -60,7 +60,7 @@ private:
     // Private data
 
         //- Radiative incident flux mapped from  the primary region / [kg/s3]
-        volScalarField QinPrimary_;
+        volScalarField qinPrimary_;
 
         //- Remaining radiative flux after removing local contribution
         volScalarField qrNet_;

src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.H

@@ -367,10 +367,10 @@ public:
         // Derived fields (calculated on-the-fly)
 
             //- Return the convective heat energy from film to wall
-            inline tmp<scalarField> Qconvw(const label patchi) const;
+            inline tmp<scalarField> qconvw(const label patchi) const;
 
             //- Return the convective heat energy from primary region to film
-            inline tmp<scalarField> Qconvp(const label patchi) const;
+            inline tmp<scalarField> qconvp(const label patchi) const;
 
 
         // Evolution

src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayerI.H

@@ -155,7 +155,7 @@ inline const filmRadiationModel& thermoSingleLayer::radiation() const
 }
 
 
-inline tmp<scalarField> thermoSingleLayer::Qconvw(const label patchi) const
+inline tmp<scalarField> thermoSingleLayer::qconvw(const label patchi) const
 {
     const scalarField htc(htcw_->h()().boundaryField()[patchi]);
     const scalarField& Tp = T_.boundaryField()[patchi];
@@ -165,7 +165,7 @@ inline tmp<scalarField> thermoSingleLayer::Qconvw(const label patchi) const
 }
 
 
-inline tmp<scalarField> thermoSingleLayer::Qconvp(const label patchi) const
+inline tmp<scalarField> thermoSingleLayer::qconvp(const label patchi) const
 {
     const scalarField htc(htcs_->h()().boundaryField()[patchi]);
     const scalarField& Tp = T_.boundaryField()[patchi];

src/regionModels/thermalBaffleModels/thermalBaffle/thermalBaffle.C

@@ -112,7 +112,7 @@ void thermalBaffle::solveEnergy()
                 const label cellId = cells[i];
 
                 Q[cellId] =
-                    Qs_.boundaryField()[patchi][localFacei]
+                    qs_.boundaryField()[patchi][localFacei]
                    /thickness_[localFacei];
 
                 rho[cellId] *= delta_.value()/thickness_[localFacei];
@@ -167,11 +167,11 @@ thermalBaffle::thermalBaffle
     nNonOrthCorr_(readLabel(solution().lookup("nNonOrthCorr"))),
     thermo_(solidThermo::New(regionMesh(), dict)),
     h_(thermo_->he()),
-    Qs_
+    qs_
     (
         IOobject
         (
-            "Qs",
+            "qs",
             regionMesh().time().timeName(),
             regionMesh(),
             IOobject::READ_IF_PRESENT,
@@ -227,11 +227,11 @@ thermalBaffle::thermalBaffle
     nNonOrthCorr_(readLabel(solution().lookup("nNonOrthCorr"))),
     thermo_(solidThermo::New(regionMesh())),
     h_(thermo_->he()),
-    Qs_
+    qs_
     (
         IOobject
         (
-            "Qs",
+            "qs",
             regionMesh().time().timeName(),
             regionMesh(),
             IOobject::READ_IF_PRESENT,
@@ -289,13 +289,13 @@ void thermalBaffle::init()
     if (oneD_ && !constantThickness_)
     {
         label patchi = intCoupledPatchIDs_[0];
-        const label Qsb = Qs_.boundaryField()[patchi].size();
+        const label qsb = qs_.boundaryField()[patchi].size();
 
-        if (Qsb!= thickness_.size())
+        if (qsb!= thickness_.size())
         {
             FatalErrorInFunction
-                << "the boundary field of Qs is "
+                << "the boundary field of qs is "
-                << Qsb << " and " << nl
+                << qsb << " and " << nl
                 << "the field 'thickness' is " << thickness_.size() << nl
                 << exit(FatalError);
         }

src/regionModels/thermalBaffleModels/thermalBaffle/thermalBaffle.H

@@ -94,7 +94,7 @@ protected:
         // Source term fields
 
             //- Surface energy source  / [J/m2/s]
-            volScalarField Qs_;
+            volScalarField qs_;
 
             //- Volumetric energy source  / [J/m3/s]
             volScalarField Q_;

src/thermophysicalModels/basic/fluidThermo/fluidThermo.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -103,6 +103,24 @@ public:
 
         // Access to thermodynamic state variables
 
+            //- Add the given density correction to the density field.
+            //  Used to update the density field following pressure solution
+            //  Limit thermo rho between rhoMin and rhoMax
+            virtual void correctRho
+            (
+                const volScalarField& deltaRho,
+                const dimensionedScalar& rhoMin,
+                const dimensionedScalar& rhoMax
+            ) = 0;
+
+
+            //- Add the given density correction to the density field.
+            //  Used to update the density field following pressure solution
+            virtual void correctRho
+            (
+                const volScalarField& deltaRho
+            ) = 0;
+
             //- Compressibility [s^2/m^2]
             virtual const volScalarField& psi() const = 0;
 

src/thermophysicalModels/basic/psiThermo/psiThermo.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -102,6 +102,20 @@ Foam::tmp<Foam::scalarField> Foam::psiThermo::rho(const label patchi) const
 }
 
 
+void Foam::psiThermo::correctRho
+(
+    const Foam::volScalarField& deltaRho,
+    const dimensionedScalar& rhoMin,
+    const dimensionedScalar& rhoMax
+)
+{}
+
+void Foam::psiThermo::correctRho
+(
+    const Foam::volScalarField& deltaRho
+)
+{}
+
 const Foam::volScalarField& Foam::psiThermo::psi() const
 {
     return psi_;

src/thermophysicalModels/basic/psiThermo/psiThermo.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -115,6 +115,23 @@ public:
 
         // Fields derived from thermodynamic state variables
 
+            //- Add the given density correction to the density field.
+            //  Used to update the density field following pressure solution.
+            //  For psiThermo does nothing.
+            virtual void correctRho
+            (
+                const volScalarField& deltaRho,
+                const dimensionedScalar& rhoMin,
+                const dimensionedScalar& rhoMax
+            );
+
+            //- Add the given density correction to the density field.
+            //  Used to update the density field following pressure solution
+            virtual void correctRho
+            (
+                const volScalarField& deltaRho
+            );
+
             //- Density [kg/m^3] - uses current value of pressure
             virtual tmp<volScalarField> rho() const;
 

src/thermophysicalModels/basic/rhoThermo/rhoThermo.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -173,6 +173,24 @@ Foam::volScalarField& Foam::rhoThermo::rho()
 }
 
 
+void Foam::rhoThermo::correctRho
+(
+    const Foam::volScalarField& deltaRho,
+    const dimensionedScalar& rhoMin,
+    const dimensionedScalar& rhoMax
+)
+{
+    rho_ += deltaRho;
+    rho_ = max(rho_, rhoMin);
+    rho_ = min(rho_, rhoMax);
+}
+
+void Foam::rhoThermo::correctRho(const Foam::volScalarField& deltaRho)
+{
+    rho_ += deltaRho;
+}
+
+
 const Foam::volScalarField& Foam::rhoThermo::psi() const
 {
     return psi_;

src/thermophysicalModels/basic/rhoThermo/rhoThermo.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -136,6 +136,20 @@ public:
             //- Return non-const access to the local density field [kg/m^3]
             virtual volScalarField& rho();
 
+            //- Add the given density correction to the density field.
+            //  Used to update the density field following pressure solution
+            //  Limit thermo rho between rhoMin and rhoMax
+            virtual void correctRho
+            (
+                const volScalarField& deltaRho,
+                const dimensionedScalar& rhoMin,
+                const dimensionedScalar& rhoMax
+            );
+
+            //- Add the given density correction to the density field.
+            //  Used to update the density field following pressure solution
+            virtual void correctRho(const volScalarField& deltaRho);
+
             //- Compressibility [s^2/m^2]
             virtual const volScalarField& psi() const;
 

src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C

@@ -184,18 +184,18 @@ updateCoeffs()
 
     const scalarField& emissivity = temissivity();
 
-    scalarField& Qem = ray.Qem().boundaryFieldRef()[patchi];
+    scalarField& qem = ray.qem().boundaryFieldRef()[patchi];
-    scalarField& Qin = ray.Qin().boundaryFieldRef()[patchi];
+    scalarField& qin = ray.qin().boundaryFieldRef()[patchi];
 
     const vector& myRayId = dom.IRay(rayId).d();
 
     // Use updated Ir while iterating over rays
-    // avoids to used lagged Qin
+    // avoids to used lagged qin
-    scalarField Ir = dom.IRay(0).Qin().boundaryField()[patchi];
+    scalarField Ir = dom.IRay(0).qin().boundaryField()[patchi];
 
     for (label rayI=1; rayI < dom.nRay(); rayI++)
     {
-        Ir += dom.IRay(rayI).Qin().boundaryField()[patchi];
+        Ir += dom.IRay(rayI).qin().boundaryField()[patchi];
     }
 
     if (solarLoad_)
@@ -221,7 +221,7 @@ updateCoeffs()
                 )/pi;
 
             // Emmited heat flux from this ray direction
-            Qem[faceI] = refValue()[faceI]*nAve[faceI];
+            qem[faceI] = refValue()[faceI]*nAve[faceI];
         }
         else
         {
@@ -231,7 +231,7 @@ updateCoeffs()
             refValue()[faceI] = 0.0; //not used
 
             // Incident heat flux on this ray direction
-            Qin[faceI] = Iw[faceI]*nAve[faceI];
+            qin[faceI] = Iw[faceI]*nAve[faceI];
         }
     }
 

src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C

@@ -177,17 +177,17 @@ updateCoeffs()
 
     const scalarField& emissivity = temissivity();
 
-    scalarField& Qem = ray.Qem().boundaryFieldRef()[patchi];
+    scalarField& qem = ray.qem().boundaryFieldRef()[patchi];
-    scalarField& Qin = ray.Qin().boundaryFieldRef()[patchi];
+    scalarField& qin = ray.qin().boundaryFieldRef()[patchi];
 
     // Use updated Ir while iterating over rays
-    // avoids to used lagged Qin
+    // avoids to used lagged qin
     /*
-    scalarField Ir = dom.IRay(0).Qin().boundaryField()[patchi];
+    scalarField Ir = dom.IRay(0).qin().boundaryField()[patchi];
 
     for (label rayI=1; rayI < dom.nRay(); rayI++)
     {
-        Ir += dom.IRay(rayI).Qin().boundaryField()[patchi];
+        Ir += dom.IRay(rayI).qin().boundaryField()[patchi];
     }
     */
 
@@ -227,7 +227,7 @@ updateCoeffs()
                 )/pi;
 
             // Emmited heat flux from this ray direction (sum over lambdaId)
-            Qem[facei] += refValue()[facei]*nAve[facei];
+            qem[facei] += refValue()[facei]*nAve[facei];
         }
         else
         {
@@ -237,7 +237,7 @@ updateCoeffs()
             refValue()[facei] = 0.0; //not used
 
             // Incident heat flux on this ray direction (sum over lambdaId)
-            Qin[facei] += Iw[facei]*nAve[facei];
+            qin[facei] += Iw[facei]*nAve[facei];
         }
     }
 

src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C

@@ -243,31 +243,31 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
         mesh_,
         dimensionedScalar("qr", dimMass/pow3(dimTime), 0.0)
     ),
-    Qem_
+    qem_
     (
         IOobject
         (
-            "Qem",
+            "qem",
             mesh_.time().timeName(),
             mesh_,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
         mesh_,
-        dimensionedScalar("Qem", dimMass/pow3(dimTime), 0.0)
+        dimensionedScalar("qem", dimMass/pow3(dimTime), 0.0)
     ),
-    Qin_
+    qin_
     (
         IOobject
         (
-            "Qin",
+            "qin",
             mesh_.time().timeName(),
             mesh_,
             IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
         mesh_,
-        dimensionedScalar("Qin", dimMass/pow3(dimTime), 0.0)
+        dimensionedScalar("qin", dimMass/pow3(dimTime), 0.0)
     ),
     a_
     (
@@ -336,31 +336,31 @@ Foam::radiation::fvDOM::fvDOM
         mesh_,
         dimensionedScalar("qr", dimMass/pow3(dimTime), 0.0)
     ),
-    Qem_
+    qem_
     (
         IOobject
         (
-            "Qem",
+            "qem",
             mesh_.time().timeName(),
             mesh_,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
         mesh_,
-        dimensionedScalar("Qem", dimMass/pow3(dimTime), 0.0)
+        dimensionedScalar("qem", dimMass/pow3(dimTime), 0.0)
     ),
-    Qin_
+    qin_
     (
         IOobject
         (
-            "Qin",
+            "qin",
             mesh_.time().timeName(),
             mesh_,
             IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
         mesh_,
-        dimensionedScalar("Qin", dimMass/pow3(dimTime), 0.0)
+        dimensionedScalar("qin", dimMass/pow3(dimTime), 0.0)
     ),
     a_
     (
@@ -509,16 +509,16 @@ void Foam::radiation::fvDOM::updateG()
 {
     G_ = dimensionedScalar("zero",dimMass/pow3(dimTime), 0.0);
     qr_ = dimensionedScalar("zero",dimMass/pow3(dimTime), 0.0);
-    Qem_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
+    qem_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
-    Qin_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
+    qin_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
 
     forAll(IRay_, rayI)
     {
         IRay_[rayI].addIntensity();
         G_ += IRay_[rayI].I()*IRay_[rayI].omega();
         qr_.boundaryFieldRef() += IRay_[rayI].qr().boundaryField();
-        Qem_.boundaryFieldRef() += IRay_[rayI].Qem().boundaryField();
+        qem_.boundaryFieldRef() += IRay_[rayI].qem().boundaryField();
-        Qin_.boundaryFieldRef() += IRay_[rayI].Qin().boundaryField();
+        qin_.boundaryFieldRef() += IRay_[rayI].qin().boundaryField();
     }
 }
 

src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.H

@@ -107,10 +107,10 @@ class fvDOM
         volScalarField qr_;
 
          //- Emmited radiative heat flux [W/m2]
-        volScalarField Qem_;
+        volScalarField qem_;
 
         //- Incidet radiative heat flux [W/m2]
-        volScalarField Qin_;
+        volScalarField qin_;
 
         //- Total absorption coefficient [1/m]
         volScalarField a_;
@@ -255,10 +255,10 @@ public:
             inline const volScalarField& qr() const;
 
             //- Const access to incident radiative heat flux field
-            inline const volScalarField& Qin() const;
+            inline const volScalarField& qin() const;
 
             //- Const access to emitted radiative heat flux field
-            inline const volScalarField& Qem() const;
+            inline const volScalarField& qem() const;
 
             //- Const access to black body
             inline const blackBodyEmission& blackBody() const;

src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOMI.H

@@ -91,15 +91,15 @@ inline const Foam::volScalarField& Foam::radiation::fvDOM::qr() const
     return qr_;
 }
 
-inline const Foam::volScalarField& Foam::radiation::fvDOM::Qin() const
+inline const Foam::volScalarField& Foam::radiation::fvDOM::qin() const
 {
-    return Qin_;
+    return qin_;
 }
 
 
-inline const Foam::volScalarField& Foam::radiation::fvDOM::Qem() const
+inline const Foam::volScalarField& Foam::radiation::fvDOM::qem() const
 {
-    return Qem_;
+    return qem_;
 }
 
 

src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C

@@ -80,31 +80,31 @@ Foam::radiation::radiativeIntensityRay::radiativeIntensityRay
         mesh_,
         dimensionedScalar("qr", dimMass/pow3(dimTime), 0.0)
     ),
-    Qin_
+    qin_
     (
         IOobject
         (
-            "Qin" + name(rayId),
+            "qin" + name(rayId),
             mesh_.time().timeName(),
             mesh_,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
         mesh_,
-        dimensionedScalar("Qin", dimMass/pow3(dimTime), 0.0)
+        dimensionedScalar("qin", dimMass/pow3(dimTime), 0.0)
     ),
-    Qem_
+    qem_
     (
         IOobject
         (
-            "Qem" + name(rayId),
+            "qem" + name(rayId),
             mesh_.time().timeName(),
             mesh_,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
         mesh_,
-        dimensionedScalar("Qem", dimMass/pow3(dimTime), 0.0)
+        dimensionedScalar("qem", dimMass/pow3(dimTime), 0.0)
     ),
     d_(Zero),
     dAve_(Zero),

src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.H

@@ -83,10 +83,10 @@ private:
         volScalarField qr_;
 
         //- Incident radiative heat flux on boundary
-        volScalarField Qin_;
+        volScalarField qin_;
 
         //- Emitted radiative heat flux on boundary
-        volScalarField Qem_;
+        volScalarField qem_;
 
         //- Direction
         vector d_;
@@ -182,16 +182,16 @@ public:
             inline volScalarField& qr();
 
             //- Return non-const access to the boundary incident heat flux
-            inline volScalarField& Qin();
+            inline volScalarField& qin();
 
             //- Return non-const access to the boundary emmited heat flux
-            inline volScalarField& Qem();
+            inline volScalarField& qem();
 
             //- Return const access to the boundary incident heat flux
-            inline const volScalarField& Qin() const;
+            inline const volScalarField& qin() const;
 
             //- Return const access to the boundary emmited heat flux
-            inline const volScalarField& Qem() const;
+            inline const volScalarField& qem() const;
 
             //- Return direction
             inline const vector& d() const;

src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRayI.H

@@ -43,28 +43,28 @@ inline Foam::volScalarField& Foam::radiation::radiativeIntensityRay::qr()
 }
 
 inline const Foam::volScalarField& Foam::radiation::
-radiativeIntensityRay::Qin() const
+radiativeIntensityRay::qin() const
 {
-    return Qin_;
+    return qin_;
 }
 
 
-inline Foam::volScalarField& Foam::radiation::radiativeIntensityRay::Qin()
+inline Foam::volScalarField& Foam::radiation::radiativeIntensityRay::qin()
 {
-    return Qin_;
+    return qin_;
 }
 
 
 inline const Foam::volScalarField& Foam::radiation::
-radiativeIntensityRay::Qem() const
+radiativeIntensityRay::qem() const
 {
-    return Qem_;
+    return qem_;
 }
 
 
-inline Foam::volScalarField& Foam::radiation::radiativeIntensityRay::Qem()
+inline Foam::volScalarField& Foam::radiation::radiativeIntensityRay::qem()
 {
-    return Qem_;
+    return qem_;
 }
 
 

src/thermophysicalModels/radiation/radiationModels/radiationModel/radiationModel.C

@@ -43,7 +43,7 @@ namespace Foam
 }
 
 const Foam::word Foam::radiation::radiationModel::externalRadHeatFieldName_ =
-    "QrExt";
+    "qrExt";
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 

src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.C

@@ -119,7 +119,7 @@ void Foam::radiation::solarLoad::updateDirectHitRadiation
 {
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
     const scalarField& V = mesh_.V();
-    volScalarField::Boundary& QrBf = Qr_.boundaryFieldRef();
+    volScalarField::Boundary& qrBf = qr_.boundaryFieldRef();
 
     forAll(hitFacesId, i)
     {
@@ -135,7 +135,7 @@ void Foam::radiation::solarLoad::updateDirectHitRadiation
 
             for (label bandI = 0; bandI < nBands_; bandI++)
             {
-                QrBf[patchID][localFaceI] +=
+                qrBf[patchID][localFaceI] +=
                     (qPrim & n[localFaceI])
                   * spectralDistribution_[bandI]
                   * absorptivity_[patchID][bandI]()[localFaceI];
@@ -167,7 +167,7 @@ void Foam::radiation::solarLoad::updateSkyDiffusiveRadiation
 {
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
     const scalarField& V = mesh_.V();
-    volScalarField::Boundary& QrBf = Qr_.boundaryFieldRef();
+    volScalarField::Boundary& qrBf = qr_.boundaryFieldRef();
 
     switch(solarCalc_.sunLoadModel())
     {
@@ -231,7 +231,7 @@ void Foam::radiation::solarLoad::updateSkyDiffusiveRadiation
                     {
                         for (label bandI = 0; bandI < nBands_; bandI++)
                         {
-                            QrBf[patchID][faceI] +=
+                            qrBf[patchID][faceI] +=
                                 (Ed + Er)
                               * spectralDistribution_[bandI]
                               * absorptivity_[patchID][bandI]()[faceI];
@@ -269,7 +269,7 @@ void Foam::radiation::solarLoad::updateSkyDiffusiveRadiation
                     {
                         for (label bandI = 0; bandI < nBands_; bandI++)
                         {
-                            QrBf[patchID][faceI] +=
+                            qrBf[patchID][faceI] +=
                                 solarCalc_.diffuseSolarRad()
                               * spectralDistribution_[bandI]
                               * absorptivity_[patchID][bandI]()[faceI];
@@ -555,9 +555,9 @@ void Foam::radiation::solarLoad::calculateQdiff
         }
     }
 
-    volScalarField::Boundary& QsBf = QsecondRad_.boundaryFieldRef();
+    volScalarField::Boundary& qsBf = qsecondRad_.boundaryFieldRef();
 
-    // Fill QsecondRad_
+    // Fill qsecondRad_
     label compactId = 0;
     forAll(includePatches_, i)
     {
@@ -566,7 +566,7 @@ void Foam::radiation::solarLoad::calculateQdiff
 
         if (pp.size() > 0)
         {
-            scalarField& Qrp = QsBf[patchID];
+            scalarField& qrp = qsBf[patchID];
 
             const labelList& coarsePatchFace =
                 coarseMesh_->patchFaceMap()[patchID];
@@ -578,7 +578,7 @@ void Foam::radiation::solarLoad::calculateQdiff
                 forAll(fineFaces, k)
                 {
                     label faceI = fineFaces[k];
-                    Qrp[faceI] = localqDiffusive[compactId];
+                    qrp[faceI] = localqDiffusive[compactId];
                 }
                 compactId ++;
             }
@@ -587,15 +587,15 @@ void Foam::radiation::solarLoad::calculateQdiff
 
     const scalarField& V = mesh_.V();
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-    volScalarField::Boundary& QrBf = Qr_.boundaryFieldRef();
+    volScalarField::Boundary& qrBf = qr_.boundaryFieldRef();
 
     forAllConstIter(labelHashSet, includePatches, iter)
     {
         const label patchID = iter.key();
-        const scalarField& qSecond = QsecondRad_.boundaryField()[patchID];
+        const scalarField& qSecond = qsecondRad_.boundaryField()[patchID];
         if (includeMappedPatchBasePatches[patchID])
         {
-            QrBf[patchID] += qSecond;
+            qrBf[patchID] += qSecond;
         }
         else
         {
@@ -630,31 +630,31 @@ Foam::radiation::solarLoad::solarLoad(const volScalarField& T)
         )
     ),
     coarseMesh_(),
-    Qr_
+    qr_
     (
         IOobject
         (
-            "Qr",
+            "qr",
             mesh_.time().timeName(),
             mesh_,
             IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
         mesh_,
-        dimensionedScalar("Qr", dimMass/pow3(dimTime), 0.0)
+        dimensionedScalar("qr", dimMass/pow3(dimTime), 0.0)
     ),
-    QsecondRad_
+    qsecondRad_
     (
         IOobject
         (
-            "QsecondRad",
+            "qsecondRad",
             mesh_.time().timeName(),
             mesh_,
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
         mesh_,
-        dimensionedScalar("QsecondRad", dimMass/pow3(dimTime), 0.0)
+        dimensionedScalar("qsecondRad", dimMass/pow3(dimTime), 0.0)
     ),
     hitFaces_(),
     Ru_
@@ -720,31 +720,31 @@ Foam::radiation::solarLoad::solarLoad
         )
     ),
     coarseMesh_(),
-    Qr_
+    qr_
     (
         IOobject
         (
-            "Qr",
+            "qr",
             mesh_.time().timeName(),
             mesh_,
             IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
         mesh_,
-        dimensionedScalar("Qr", dimMass/pow3(dimTime), 0.0)
+        dimensionedScalar("qr", dimMass/pow3(dimTime), 0.0)
     ),
-    QsecondRad_
+    qsecondRad_
     (
         IOobject
         (
-            "QsecondRad",
+            "qsecondRad",
             mesh_.time().timeName(),
             mesh_,
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
         mesh_,
-        dimensionedScalar("QsecondRad", dimMass/pow3(dimTime), 0.0)
+        dimensionedScalar("qsecondRad", dimMass/pow3(dimTime), 0.0)
     ),
     hitFaces_(),
     Ru_
@@ -812,7 +812,7 @@ Foam::radiation::solarLoad::solarLoad
         )
     ),
     coarseMesh_(),
-    Qr_
+    qr_
     (
         IOobject
         (
@@ -823,20 +823,20 @@ Foam::radiation::solarLoad::solarLoad
             IOobject::AUTO_WRITE
         ),
         mesh_,
-        dimensionedScalar("Qr", dimMass/pow3(dimTime), 0.0)
+        dimensionedScalar("qr", dimMass/pow3(dimTime), 0.0)
     ),
-    QsecondRad_
+    qsecondRad_
     (
         IOobject
         (
-            "QsecondRad",
+            "qsecondRad",
             mesh_.time().timeName(),
             mesh_,
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
         mesh_,
-        dimensionedScalar("QsecondRad", dimMass/pow3(dimTime), 0.0)
+        dimensionedScalar("qsecondRad", dimMass/pow3(dimTime), 0.0)
     ),
     hitFaces_(),
     Ru_
@@ -942,13 +942,13 @@ void Foam::radiation::solarLoad::calculate()
     }
 
     bool facesChanged = updateHitFaces();
-    volScalarField::Boundary& QrBf = Qr_.boundaryFieldRef();
+    volScalarField::Boundary& qrBf = qr_.boundaryFieldRef();
 
     if (facesChanged)
     {
-        // Reset Ru and Qr
+        // Reset Ru and qr
         Ru_ = dimensionedScalar("Ru", dimMass/dimLength/pow3(dimTime), 0.0);
-        QrBf = 0.0;
+        qrBf = 0.0;
 
         // Add direct hit radation
         const labelList& hitFacesId = hitFaces_->rayStartFaces();

src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.H

@@ -105,10 +105,10 @@ private:
         autoPtr<singleCellFvMesh> coarseMesh_;
 
         //- Net radiative heat flux [W/m2]
-        volScalarField Qr_;
+        volScalarField qr_;
 
         //- Secondary solar radiative heat flux [W/m2]
-        volScalarField QsecondRad_;
+        volScalarField qsecondRad_;
 
         //- Direct hit faces Ids
         autoPtr<faceShading> hitFaces_;
@@ -143,10 +143,10 @@ private:
         //- Face-compact map
         labelListIOList visibleFaceFaces_;
 
-        //- Couple solids through mapped boundary patch using Qr (default:true)
+        //- Couple solids through mapped boundary patch using qr (default:true)
         bool solidCoupled_;
 
-        //- Couple wall patches using Qr (default:false)
+        //- Couple wall patches using qr (default:false)
         bool wallCoupled_;
 
         //- Absorptivity list

src/thermophysicalModels/thermophysicalPropertiesFvPatchFields/liquidProperties/humidityTemperatureCoupledMixed/humidityTemperatureCoupledMixedFvPatchScalarField.C

@@ -149,8 +149,8 @@ humidityTemperatureCoupledMixedFvPatchScalarField
     rhoName_("rho"),
     muName_("thermo:mu"),
     TnbrName_("T"),
-    QrNbrName_("none"),
+    qrNbrName_("none"),
-    QrName_("none"),
+    qrName_("none"),
     specieName_("none"),
     liquid_(nullptr),
     liquidDict_(nullptr),
@@ -189,8 +189,8 @@ humidityTemperatureCoupledMixedFvPatchScalarField
     rhoName_(psf.rhoName_),
     muName_(psf.muName_),
     TnbrName_(psf.TnbrName_),
-    QrNbrName_(psf.QrNbrName_),
+    qrNbrName_(psf.qrNbrName_),
-    QrName_(psf.QrName_),
+    qrName_(psf.qrName_),
     specieName_(psf.specieName_),
     liquid_(psf.liquid_, false),
     liquidDict_(psf.liquidDict_),
@@ -224,8 +224,8 @@ humidityTemperatureCoupledMixedFvPatchScalarField
     rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
     muName_(dict.lookupOrDefault<word>("mu", "thermo:mu")),
     TnbrName_(dict.lookupOrDefault<word>("Tnbr", "T")),
-    QrNbrName_(dict.lookupOrDefault<word>("QrNbr", "none")),
+    qrNbrName_(dict.lookupOrDefault<word>("qrNbr", "none")),
-    QrName_(dict.lookupOrDefault<word>("Qr", "none")),
+    qrName_(dict.lookupOrDefault<word>("qr", "none")),
     specieName_(dict.lookupOrDefault<word>("specie", "none")),
     liquid_(nullptr),
     liquidDict_(),
@@ -348,8 +348,8 @@ humidityTemperatureCoupledMixedFvPatchScalarField
     rhoName_(psf.rhoName_),
     muName_(psf.muName_),
     TnbrName_(psf.TnbrName_),
-    QrNbrName_(psf.QrNbrName_),
+    qrNbrName_(psf.qrNbrName_),
-    QrName_(psf.QrName_),
+    qrName_(psf.qrName_),
     specieName_(psf.specieName_),
     liquid_(psf.liquid_, false),
     liquidDict_(psf.liquidDict_),
@@ -675,26 +675,26 @@ void Foam::humidityTemperatureCoupledMixedFvPatchScalarField::updateCoeffs()
         mpp.distribute(dmHfgNbr);
     }
 
-    // Obtain Rad heat (Qr)
+    // Obtain Rad heat (qr)
-    scalarField Qr(Tp.size(), 0.0);
+    scalarField qr(Tp.size(), 0.0);
-    if (QrName_ != "none")
+    if (qrName_ != "none")
     {
-        Qr = patch().lookupPatchField<volScalarField, scalar>(QrName_);
+        qr = patch().lookupPatchField<volScalarField, scalar>(qrName_);
     }
 
-    scalarField QrNbr(Tp.size(), 0.0);
+    scalarField qrNbr(Tp.size(), 0.0);
-    if (QrNbrName_ != "none")
+    if (qrNbrName_ != "none")
     {
-        QrNbr = nbrPatch.lookupPatchField<volScalarField, scalar>(QrNbrName_);
+        qrNbr = nbrPatch.lookupPatchField<volScalarField, scalar>(qrNbrName_);
-        mpp.distribute(QrNbr);
+        mpp.distribute(qrNbr);
     }
 
     const scalarField dmHfg(dmHfgNbr + dmHfg_);
 
     const scalarField mpCpdt(mpCpTpNbr + mpCpTp_);
 
-    // Qr > 0 (heat up the wall)
+    // qr > 0 (heat up the wall)
-    scalarField alpha(KDeltaNbr + mpCpdt - (Qr + QrNbr)/Tp);
+    scalarField alpha(KDeltaNbr + mpCpdt - (qr + qrNbr)/Tp);
 
     valueFraction() = alpha/(alpha + myKDelta_);
 
@@ -742,8 +742,8 @@ void Foam::humidityTemperatureCoupledMixedFvPatchScalarField::write
     writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
     writeEntryIfDifferent<word>(os, "mu", "thermo:mu", muName_);
     writeEntryIfDifferent<word>(os, "Tnbr", "T", TnbrName_);
-    writeEntryIfDifferent<word>(os, "QrNbr", "none", QrNbrName_);
+    writeEntryIfDifferent<word>(os, "qrNbr", "none", qrNbrName_);
-    writeEntryIfDifferent<word>(os, "Qr", "none", QrName_);
+    writeEntryIfDifferent<word>(os, "qr", "none", qrName_);
 
     if (fluid_)
     {

src/thermophysicalModels/thermophysicalPropertiesFvPatchFields/liquidProperties/humidityTemperatureCoupledMixed/humidityTemperatureCoupledMixedFvPatchScalarField.H

@@ -212,10 +212,10 @@ private:
             const word TnbrName_;
 
             //- Name of the radiative heat flux in the neighbout region
-            const word QrNbrName_;
+            const word qrNbrName_;
 
             //- Name of the radiative heat flux field
-            const word QrName_;
+            const word qrName_;
 
             //- Name of the species on which the mass transfered (default H2O)
             const word specieName_;

tutorials/combustion/fireFoam/LES/compartmentFire/0/T

@@ -43,8 +43,8 @@ boundaryField
         type            compressible::turbulentTemperatureRadCoupledMixed;
         value           uniform 294.75;
         Tnbr            T;
-        QrNbr           none;
+        qrNbr           none;
-        Qr              Qr;
+        qr              qr;
         kappaMethod     fluidThermo;
     }
     region0_to_panelRegion_internalWallPanel_top
@@ -52,8 +52,8 @@ boundaryField
         type            compressible::turbulentTemperatureRadCoupledMixed;
         value           uniform 294.75;
         Tnbr            T;
-        QrNbr           none;
+        qrNbr           none;
-        Qr              Qr;
+        qr              qr;
         kappaMethod     fluidThermo;
     }
     region0_to_panelRegion_internalWallPanel_bottom
@@ -61,8 +61,8 @@ boundaryField
         type            compressible::turbulentTemperatureRadCoupledMixed;
         value           uniform 294.75;
         Tnbr            T;
-        QrNbr           none;
+        qrNbr           none;
-        Qr              Qr;
+        qr              qr;
         kappaMethod     fluidThermo;
     }
 }

tutorials/combustion/fireFoam/LES/compartmentFire/0/panelRegion/T

@@ -47,8 +47,8 @@ boundaryField
         type                compressible::turbulentTemperatureRadCoupledMixed;
         neighbourField      T;
         kappaMethod         solidThermo;
-        QrNbr               Qr;
+        qrNbr               qr;
-        Qr                  none;
+        qr                  none;
         value               $internalField;
     }
 }

tutorials/combustion/fireFoam/LES/compartmentFire/system/controlDict

@@ -70,7 +70,7 @@ functions
         fields              (phi);
     }
 
-    wallPanel_Qin
+    wallPanel_qin
     {
         type                patchProbes;
         libs                ("libsampling.so");
@@ -84,10 +84,10 @@ functions
             (0.2  0.2 0.01)  // HF3
             (0.0  0.4 0.01)  // HF4
         );
-        fields              (Qin);
+        fields              (qin);
     }
 
-    inletQr_Qin
+    inletqr_qin
     {
         type                patchProbes;
         libs                ("libsampling.so");
@@ -103,7 +103,7 @@ functions
             (-0.02 0.0  0.02) // HF4
             (-0.02 0.0 -0.02) // HF5
         );
-        fields    (Qr Qin);
+        fields    (qr qin);
     }
 
     thermoCouple

tutorials/combustion/fireFoam/LES/simplePMMApanel/0/panelRegion/T

@@ -38,7 +38,7 @@ boundaryField
         type                fixedIncidentRadiation;
         kappaMethod         solidThermo;
         kappa               none;
-        QrIncident          uniform 60000.0;   //W
+        qrIncident          uniform 60000.0;   //W
         value               uniform 298.15;
     }
 

tutorials/combustion/fireFoam/LES/simplePMMApanel/constant/pyrolysisZones

@@ -25,7 +25,7 @@ FoamFile
 
         reactingOneDimCoeffs
         {
-            QrHSource       no; //Energy source term due in depht radiation
+            qrHSource       no; //Energy source term due in depht radiation
 
             filmCoupled     false;
 

tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/system/controlDict

@@ -7,7 +7,7 @@
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
-    version         2;
+    version         2.0;
     format          ascii;
     class           dictionary;
     location        "system";
@@ -23,7 +23,7 @@ startTime       0;
 
 stopAt          endTime;
 
-endTime         7000;
+endTime         6000;
 
 deltaT          1;
 

tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSolution

@@ -54,6 +54,7 @@ PIMPLE
     nCorrectors         1;
     nNonOrthogonalCorrectors 0;
     consistent          yes;
+    SIMPLErho           yes;
 
     pMaxFactor          1.5;
     pMinFactor          0.9;

tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/T

@@ -0,0 +1,42 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 300;
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            fixedValue;
+        value           uniform 350;
+    }
+
+    inlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      $internalField;
+    }
+}
+
+// ************************************************************************* //

tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/U

@@ -0,0 +1,43 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            noSlip;
+    }
+
+    inlet
+    {
+        type            flowRateInletVelocity;
+        massFlowRate    constant 5;
+        rhoInlet        1000;    // Guess for rho
+    }
+
+    outlet
+    {
+        type            pressureInletOutletVelocity;
+        value           $internalField;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/alphat

@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alphat;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            compressible::alphatWallFunction;
+        Prt             0.85;
+        value           uniform 0;
+    }
+
+    inlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+
+    outlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/epsilon

@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -3 0 0 0 0];
+
+internalField   uniform 200;
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            epsilonWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 200;
+    }
+
+    inlet
+    {
+        type            turbulentMixingLengthDissipationRateInlet;
+        mixingLength    0.005;
+        value           uniform 200;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 200;
+        value           uniform 200;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/k

@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 1;
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            kqRWallFunction;
+        value           uniform 1;
+    }
+
+    inlet
+    {
+        type            turbulentIntensityKineticEnergyInlet;
+        intensity       0.05;
+        value           uniform 1;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 1;
+        value           uniform 1;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/nut

@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      nut;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            nutkWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0;
+    }
+
+    inlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+
+    outlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/p

@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 1e5;
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            zeroGradient;
+    }
+
+    inlet
+    {
+        type            zeroGradient;
+    }
+
+    outlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/constant/thermophysicalProperties

@@ -0,0 +1,31 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         pureMixture;
+    properties      liquid;
+    energy          sensibleInternalEnergy;
+}
+
+mixture
+{
+    H2O;
+}
+
+// ************************************************************************* //

tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/constant/turbulenceProperties

@@ -0,0 +1,30 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  RAS;
+
+RAS
+{
+    RASModel        kEpsilon;
+
+    turbulence      on;
+
+    printCoeffs     on;
+}
+
+
+// ************************************************************************* //

tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/blockMeshDict

@@ -0,0 +1,127 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 0.001;
+
+vertices
+(
+    // front-plane: z = +25mm
+    // inlet region
+    (   -50  25   25)           // pt 0
+    (     0  25   25)           // pt 1
+    (   -50  75   25)           // pt 2
+    (     0  75   25)           // pt 3
+    // outlet region
+    (  -500 -75   25)           // pt 4
+    (     0 -75   25)           // pt 5
+    (  -500 -25   25)           // pt 6
+    (     0 -25   25)           // pt 7
+    // bend mid-points
+    (    25   0   25)           // pt 8
+    (    75   0   25)           // pt 9
+    // back-plane: z = -25mm
+    // inlet region
+    (   -50  25   -25)          // pt 0 + 10
+    (     0  25   -25)          // pt 1 + 10
+    (   -50  75   -25)          // pt 2 + 10
+    (     0  75   -25)          // pt 3 + 10
+    // outlet region
+    (  -500 -75   -25)          // pt 4 + 10
+    (     0 -75   -25)          // pt 5 + 10
+    (  -500 -25   -25)          // pt 7 + 10
+    (     0 -25   -25)          // pt 8 + 10
+    // bend mid-points
+    (    25   0   -25)          // pt 8 + 10
+    (    75   0   -25)          // pt 9 + 10
+);
+
+blocks
+(
+    hex (0 1 11 10  2 3 13 12) inlet  ( 20 20 20)  simpleGrading (1 1 1)
+    hex (4 5 15 14  6 7 17 16) outlet (200 20 20)  simpleGrading (1 1 1)
+
+    hex (1 8 18 11  3 9 19 13) bend1  ( 30 20 20)  simpleGrading (1 1 1)
+    hex (5 9 19 15  7 8 18 17) bend2  ( 30 20 20)  simpleGrading (1 1 1)
+);
+
+edges
+(
+   // block 2
+   arc  1  8  ( 17.678  17.678  25)
+   arc 11 18  ( 17.678  17.678 -25)
+   arc  3  9  ( 53.033  53.033  25)
+   arc 13 19  ( 53.033  53.033 -25)
+   // block 3
+   arc  7  8  ( 17.678  -17.678  25)
+   arc 17 18  ( 17.678  -17.678 -25)
+   arc  5  9  ( 53.033  -53.033  25)
+   arc 15 19  ( 53.033  -53.033 -25)
+);
+
+boundary
+(
+        // is there no way of defining all my 'defaultFaces' to be 'wall'?
+    Default_Boundary_Region
+    {
+        type wall;
+        faces
+        (
+            // block0
+            ( 0 1 3 2 )
+            ( 11 10 12 13 )
+            ( 0 10 11 1 )
+            ( 2 3 13 12 )
+            // block1
+            ( 4 5 7 6 )
+            ( 15 14 16 17 )
+            ( 4 14 15 5 )
+            ( 6 7 17 16 )
+            // block2
+            ( 1 8 9 3 )
+            ( 18 11 13 19 )
+            ( 3 9 19 13 )
+            ( 1 11 18 8 )
+            // block3
+            ( 5 9 8 7 )
+            ( 19 15 17 18 )
+            ( 5 15 19 9 )
+            ( 7 8 18 17 )
+        );
+    }
+    inlet
+    {
+        type patch;
+        faces
+        (
+            (0 2 12 10)
+        );
+    }
+
+    outlet
+    {
+        type patch;
+        faces
+        (
+            (4 6 16 14)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //

tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/controlDict

@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     rhoPimpleFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         0.5;
+
+deltaT          2e-3;
+
+writeControl    timeStep;
+
+writeInterval   10;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+graphFormat     raw;
+
+runTimeModifiable true;
+
+
+// ************************************************************************* //

tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/decomposeParDict

@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 8;
+
+method          hierarchical;
+
+simpleCoeffs
+{
+    n               (8 1 1);
+    delta           0.001;
+}
+
+hierarchicalCoeffs
+{
+    n               (4 2 1);
+    delta           0.001;
+    order           xyz;
+}
+
+manualCoeffs
+{
+    dataFile        "";
+}
+
+distributed     no;
+
+roots           ( );
+
+
+// ************************************************************************* //

tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/fvSchemes

@@ -0,0 +1,60 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default             Euler;
+}
+
+gradSchemes
+{
+    default             Gauss linear;
+
+    limited             cellLimited Gauss linear 1;
+}
+
+divSchemes
+{
+    default             none;
+
+    div(phi,U)          Gauss linearUpwind limited;
+    div(phi,e)          Gauss linearUpwind limited;
+    div(phi,epsilon)    Gauss linearUpwind limited;
+    div(phi,k)          Gauss linearUpwind limited;
+    div(phi,K)          Gauss linearUpwind limited;
+    div(phiv,p)         Gauss linearUpwind limited;
+
+    div(((rho*nuEff)*dev2(T(grad(U)))))      Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+
+// ************************************************************************* //

tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/fvSolution

@@ -0,0 +1,63 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    "rho.*"
+    {
+        solver nthn;
+    }
+
+    "p.*"
+    {
+        solver          GAMG;
+        smoother        GaussSeidel;
+
+        tolerance       1e-7;
+        relTol          0;
+    }
+
+    "(U|e|k|epsilon).*"
+    {
+        solver          PBiCGStab;
+        preconditioner  DILU;
+
+        tolerance       1e-7;
+        relTol          0;
+    }
+}
+
+PIMPLE
+{
+    nCorrectors     3;
+    nNonOrthogonalCorrectors 0;
+    pMinFactor      0.1;
+    pMaxFactor      1.5;
+
+    transonic       no;
+    consistent      no;
+}
+
+relaxationFactors
+{
+    equations
+    {
+        ".*"            1;
+    }
+}
+
+// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.orig/include/1DBaffle/1DTemperatureMasterBafflePatches

@@ -11,7 +11,7 @@ T
     type   compressible::thermalBaffle1D<hConstSolidThermoPhysics>;
 
     thickness       uniform 0.005;  // thickness [m]
-    Qs              uniform 100;    // heat flux [W/m2]
+    qs              uniform 100;    // heat flux [W/m2]
 
     # include "1DbaffleSolidThermo"
 

tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/0/air/T

@@ -58,8 +58,8 @@ boundaryField
         inletValue      uniform 300;
         Tnbr            T;
         kappaMethod     fluidThermo;
-        QrNbr           none;
+        qrNbr           none;
-        Qr              Qr;
+        qr              qr;
         kappa           none;
     }
     air_to_solid
@@ -69,8 +69,8 @@ boundaryField
         inletValue      uniform 300;
         Tnbr            T;
         kappaMethod     fluidThermo;
-        QrNbr           none;
+        qrNbr           none;
-        Qr              Qr;
+        qr              qr;
         kappa           none;
     }
 }

tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/0/air/qr

@@ -11,7 +11,7 @@ FoamFile
     format      ascii;
     class       volScalarField;
     location    "0/air";
-    object      Qr;
+    object      qr;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/0/floor/T

@@ -64,8 +64,8 @@ boundaryField
         value           uniform 300;
         Tnbr            T;
         kappaMethod     solidThermo;
-        QrNbr           Qr;
+        qrNbr           qr;
-        Qr              none;
+        qr              none;
         kappa           none;
     }
     floor_to_solid
@@ -74,8 +74,8 @@ boundaryField
         value           uniform 300;
         Tnbr            T;
         kappaMethod     solidThermo;
-        QrNbr           none;
+        qrNbr           none;
-        Qr              none;
+        qr              none;
         kappa           none;
     }
 }

tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/0/solid/T

@@ -32,8 +32,8 @@ boundaryField
         value           uniform 300;
         Tnbr            T;
         kappaMethod     solidThermo;
-        QrNbr           Qr;
+        qrNbr           qr;
-        Qr              none;
+        qr              none;
         kappa           none;
     }
     solid_to_floor
@@ -42,8 +42,8 @@ boundaryField
         value           uniform 300;
         Tnbr            T;
         kappaMethod     solidThermo;
-        QrNbr           none;
+        qrNbr           none;
-        Qr              none;
+        qr              none;
         kappa           none;
     }
 }

tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/Allrun.pre

@@ -16,7 +16,7 @@ rm -r system/domain3
 # remove fluid fields from solid regions (important for post-processing)
 for i in solid floor
 do
-   rm -f 0*/$i/{rho,mut,alphat,epsilon,k,U,p_rgh,Qr,G,IDefault}
+   rm -f 0*/$i/{rho,mut,alphat,epsilon,k,U,p_rgh,qr,G,IDefault}
 done
 
 for i in air solid floor

tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/constant/air/radiationProperties

@@ -72,8 +72,8 @@ solarLoadCoeffs
         C   0.058;  // Model constant
 
     // Radiative flux coupling flags
-    solidCoupled    true;  //Couple through Qr the solid regions (default true)
+    solidCoupled    true;  //Couple through qr the solid regions (default true)
-    wallCoupled     false; //Couple through Qr wall patches (default false)
+    wallCoupled     false; //Couple through qr wall patches (default false)
 }
 
 

tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/system/air/changeDictionaryDict

@@ -69,8 +69,8 @@ dictionaryReplacement
                 type            compressible::turbulentTemperatureRadCoupledMixed;
                 Tnbr            T;
                 kappaMethod     fluidThermo;
-                QrNbr           none;
+                qrNbr           none;
-                Qr              Qr;
+                qr              qr;
                 kappa           none;
                 value           uniform 300;
             }

tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/system/floor/changeDictionaryDict

@@ -40,8 +40,8 @@ dictionaryReplacement
                 type            compressible::turbulentTemperatureRadCoupledMixed;
                 Tnbr            T;
                 kappaMethod     solidThermo;
-                QrNbr           none;
+                qrNbr           none;
-                Qr              none;
+                qr              none;
                 kappa           none;
                 value           uniform 300;
             }
@@ -51,8 +51,8 @@ dictionaryReplacement
                 type            compressible::turbulentTemperatureRadCoupledMixed;
                 Tnbr            T;
                 kappaMethod     solidThermo;
-                QrNbr           Qr;
+                qrNbr           qr;
-                Qr              none;
+                qr              none;
                 kappa           none;
                 value           uniform 300;
             }

tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/system/solid/changeDictionaryDict

@@ -34,8 +34,8 @@ dictionaryReplacement
                 type            compressible::turbulentTemperatureRadCoupledMixed;
                 Tnbr            T;
                 kappaMethod     solidThermo;
-                QrNbr           Qr;
+                qrNbr           qr;
-                Qr              none;
+                qr              none;
                 kappa           none;
                 value           uniform 300;
             }
@@ -45,8 +45,8 @@ dictionaryReplacement
                 type            compressible::turbulentTemperatureRadCoupledMixed;
                 Tnbr            T;
                 kappaMethod     solidThermo;
-                QrNbr           none;
+                qrNbr           none;
-                Qr              none;
+                qr              none;
                 kappa           none;
                 value           uniform 300;
             }