From a1c8cde310f48925143781a940811a5fca8a87c2 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Fri, 24 Feb 2017 11:18:01 +0000
Subject: [PATCH] rhoSimpleFoam: added support for compressible liquid flows

rhoSimpleFoam now instantiates the lower-level fluidThermo which instantiates
either a psiThermo or rhoThermo according to the 'type' specification in
thermophysicalProperties, e.g.

thermoType
{
    type            hePsiThermo;
    mixture         pureMixture;
    transport       sutherland;
    thermo          janaf;
    equationOfState perfectGas;
    specie          specie;
    energy          sensibleInternalEnergy;
}

instantiates a psiThermo for a perfect gas with JANAF thermodynamics, whereas

thermoType
{
    type            heRhoThermo;
    mixture         pureMixture;
    properties      liquid;
    energy          sensibleInternalEnergy;
}

mixture
{
    H2O;
}

instantiates a rhoThermo for water, see new tutorial
compressible/rhoSimpleFoam/squareBendLiq.

In order to support complex equations of state the pressure can no longer be
unlimited and rhoSimpleFoam now limits the pressure rather than the density to
handle start-up more robustly.

For backward compatibility 'rhoMin' and 'rhoMax' can still be used in the SIMPLE
sub-dictionary of fvSolution which are converted into 'pMax' and 'pMin' but it
is better to set either 'pMax' and 'pMin' directly or use the more convenient
'pMinFactor' and 'pMinFactor' from which 'pMax' and 'pMin' are calculated using
the fixed boundary pressure or reference pressure e.g.

SIMPLE
{
    nNonOrthogonalCorrectors 0;

    pMinFactor      0.1;
    pMaxFactor      1.5;

    transonic       yes;
    consistent      yes;

    residualControl
    {
        p               1e-3;
        U               1e-4;
        e               1e-3;
        "(k|epsilon|omega)" 1e-3;
    }
}
---
 .../compressible/rhoSimpleFoam/createFields.H |  33 +---
 .../solvers/compressible/rhoSimpleFoam/pEqn.H |  22 ++-
 .../compressible/rhoSimpleFoam/pcEqn.H        |  20 ++-
 .../rhoPorousSimpleFoam/createFields.H        |  44 +----
 .../rhoSimpleFoam/rhoPorousSimpleFoam/pEqn.H  |  15 +-
 .../rhoPorousSimpleFoam/rhoPorousSimpleFoam.C |   7 +-
 .../rhoSimpleFoam/rhoSimpleFoam.C             |   5 +-
 .../system/fvSchemes                          |   4 +-
 .../system/fvSolution                         |   6 +-
 .../CompressibleTurbulenceModel.H             |  15 +-
 .../TurbulenceModel/TurbulenceModel.H         |   6 +-
 src/finiteVolume/Make/files                   |   1 +
 .../general/findRefCell/findRefCell.C         |  14 +-
 .../general/findRefCell/findRefCell.H         |  10 +-
 .../general/pressureControl/pressureControl.C | 155 ++++++++++++++++++
 .../general/pressureControl/pressureControl.H | 108 ++++++++++++
 .../pressureControl/pressureControlI.H        |  40 +++++
 .../angledDuctExplicit/system/fvSolution      |   4 +-
 .../angledDuctImplicit/system/fvSolution      |   4 +-
 .../system/fvSolution                         |   4 +-
 .../squareBend/system/fvSolution              |   5 +-
 .../rhoSimpleFoam/squareBendLiq/0/T           |  42 +++++
 .../rhoSimpleFoam/squareBendLiq/0/U           |  42 +++++
 .../rhoSimpleFoam/squareBendLiq/0/alphat      |  43 +++++
 .../rhoSimpleFoam/squareBendLiq/0/epsilon     |  47 ++++++
 .../rhoSimpleFoam/squareBendLiq/0/k           |  44 +++++
 .../rhoSimpleFoam/squareBendLiq/0/nut         |  45 +++++
 .../rhoSimpleFoam/squareBendLiq/0/p           |  39 +++++
 .../constant/thermophysicalProperties         |  31 ++++
 .../constant/turbulenceProperties             |  30 ++++
 .../squareBendLiq/system/blockMeshDict        | 127 ++++++++++++++
 .../squareBendLiq/system/controlDict          |  51 ++++++
 .../squareBendLiq/system/decomposeParDict     |  45 +++++
 .../squareBendLiq/system/fvSchemes            |  59 +++++++
 .../squareBendLiq/system/fvSolution           |  71 ++++++++
 35 files changed, 1115 insertions(+), 123 deletions(-)
 create mode 100644 src/finiteVolume/cfdTools/general/pressureControl/pressureControl.C
 create mode 100644 src/finiteVolume/cfdTools/general/pressureControl/pressureControl.H
 create mode 100644 src/finiteVolume/cfdTools/general/pressureControl/pressureControlI.H
 create mode 100644 tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/T
 create mode 100644 tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/U
 create mode 100644 tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/alphat
 create mode 100644 tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/epsilon
 create mode 100644 tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/k
 create mode 100644 tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/nut
 create mode 100644 tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/p
 create mode 100644 tutorials/compressible/rhoSimpleFoam/squareBendLiq/constant/thermophysicalProperties
 create mode 100644 tutorials/compressible/rhoSimpleFoam/squareBendLiq/constant/turbulenceProperties
 create mode 100644 tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/blockMeshDict
 create mode 100644 tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/controlDict
 create mode 100644 tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/decomposeParDict
 create mode 100644 tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/fvSchemes
 create mode 100644 tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/fvSolution

diff --git a/applications/solvers/compressible/rhoSimpleFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/createFields.H
index e8c2ce6d548..4671347b66a 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/createFields.H
@@ -1,10 +1,10 @@
 Info<< "Reading thermophysical properties\n" << endl;
 
-autoPtr<psiThermo> pThermo
+autoPtr<fluidThermo> pThermo
 (
-    psiThermo::New(mesh)
+    fluidThermo::New(mesh)
 );
-psiThermo& thermo = pThermo();
+fluidThermo& thermo = pThermo();
 thermo.validate(args.executable(), "h", "e");
 
 volScalarField rho
@@ -38,35 +38,10 @@ volVectorField U
 
 #include "compressibleCreatePhi.H"
 
-
-label pRefCell = 0;
-scalar pRefValue = 0.0;
-setRefCell(p, simple.dict(), pRefCell, pRefValue);
+pressureControl pressureControl(p, rho, simple.dict());
 
 mesh.setFluxRequired(p.name());
 
-dimensionedScalar rhoMax
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMax",
-        simple.dict(),
-        dimDensity,
-        GREAT
-    )
-);
-
-dimensionedScalar rhoMin
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMin",
-        simple.dict(),
-        dimDensity,
-        0
-    )
-);
-
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence
 (
diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
index e46f2a66913..06fa890c057 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
@@ -30,7 +30,11 @@
             // Relax the pressure equation to ensure diagonal-dominance
             pEqn.relax();
 
-            pEqn.setReference(pRefCell, pRefValue);
+            pEqn.setReference
+            (
+                pressureControl.refCell(),
+                pressureControl.refValue()
+            );
 
             pEqn.solve();
 
@@ -60,7 +64,11 @@
                 fvOptions(psi, p, rho.name())
             );
 
-            pEqn.setReference(pRefCell, pRefValue);
+            pEqn.setReference
+            (
+                pressureControl.refCell(),
+                pressureControl.refValue()
+            );
 
             pEqn.solve();
 
@@ -81,6 +89,8 @@
     U.correctBoundaryConditions();
     fvOptions.correct(U);
 
+    pressureControl.limit(p);
+
     // For closed-volume cases adjust the pressure and density levels
     // to obey overall mass continuity
     if (closedVolume)
@@ -89,16 +99,12 @@
             /fvc::domainIntegrate(psi);
     }
 
+    p.correctBoundaryConditions();
+
     rho = thermo.rho();
-    rho = max(rho, rhoMin);
-    rho = min(rho, rhoMax);
 
     if (!simple.transonic())
     {
         rho.relax();
     }
-
-    Info<< "rho max/min : "
-        << max(rho).value() << " "
-        << min(rho).value() << endl;
 }
diff --git a/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H
index 67b962f8050..0dbde829608 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H
@@ -40,7 +40,11 @@ if (simple.transonic())
         // Relax the pressure equation to maintain diagonal dominance
         pEqn.relax();
 
-        pEqn.setReference(pRefCell, pRefValue);
+        pEqn.setReference
+        (
+            pressureControl.refCell(),
+            pressureControl.refValue()
+        );
 
         pEqn.solve();
 
@@ -75,7 +79,11 @@ else
             fvOptions(psi, p, rho.name())
         );
 
-        pEqn.setReference(pRefCell, pRefValue);
+        pEqn.setReference
+        (
+            pressureControl.refCell(),
+            pressureControl.refValue()
+        );
 
         pEqn.solve();
 
@@ -97,6 +105,8 @@ U = HbyA - rAtU*fvc::grad(p);
 U.correctBoundaryConditions();
 fvOptions.correct(U);
 
+pressureControl.limit(p);
+
 // For closed-volume cases adjust the pressure and density levels
 // to obey overall mass continuity
 if (closedVolume)
@@ -105,14 +115,12 @@ if (closedVolume)
         /fvc::domainIntegrate(psi);
 }
 
+p.correctBoundaryConditions();
+
 // Recalculate density from the relaxed pressure
 rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
 
 if (!simple.transonic())
 {
     rho.relax();
 }
-
-Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/createFields.H
index 5daef0c2de9..4671347b66a 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/createFields.H
@@ -1,10 +1,10 @@
 Info<< "Reading thermophysical properties\n" << endl;
 
-autoPtr<rhoThermo> pThermo
+autoPtr<fluidThermo> pThermo
 (
-    rhoThermo::New(mesh)
+    fluidThermo::New(mesh)
 );
-rhoThermo& thermo = pThermo();
+fluidThermo& thermo = pThermo();
 thermo.validate(args.executable(), "h", "e");
 
 volScalarField rho
@@ -38,46 +38,10 @@ volVectorField U
 
 #include "compressibleCreatePhi.H"
 
-
-label pRefCell = 0;
-scalar pRefValue = 0.0;
-setRefCell(p, simple.dict(), pRefCell, pRefValue);
+pressureControl pressureControl(p, rho, simple.dict());
 
 mesh.setFluxRequired(p.name());
 
-dimensionedScalar rhoMax
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMax",
-        simple.dict(),
-        dimDensity,
-        GREAT
-    )
-);
-
-dimensionedScalar rhoMin
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMin",
-        simple.dict(),
-        dimDensity,
-        0
-    )
-);
-
-dimensionedScalar pMin
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "pMin",
-        simple.dict(),
-        dimPressure,
-        10
-    )
-);
-
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence
 (
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/pEqn.H
index ca8725b4994..4167f80e599 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/pEqn.H
@@ -48,7 +48,11 @@
 
         fvScalarMatrix& pEqn = tpEqn.ref();
 
-        pEqn.setReference(pRefCell, pRefValue);
+        pEqn.setReference
+        (
+            pressureControl.refCell(),
+            pressureControl.refValue()
+        );
 
         pEqn.solve();
 
@@ -75,6 +79,8 @@
     U.correctBoundaryConditions();
     fvOptions.correct(U);
 
+    pressureControl.limit(p);
+
     // For closed-volume cases adjust the pressure and density levels
     // to obey overall mass continuity
     if (closedVolume)
@@ -83,13 +89,6 @@
             /fvc::domainIntegrate(psi);
     }
 
-    p = max(p, pMin);
-
     rho = thermo.rho();
-    rho = max(rho, rhoMin);
-    rho = min(rho, rhoMax);
     rho.relax();
-    Info<< "rho max/min : "
-        << max(rho).value() << " "
-        << min(rho).value() << endl;
 }
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/rhoPorousSimpleFoam.C b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/rhoPorousSimpleFoam.C
index 6cca40e2c0d..29839274e8d 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/rhoPorousSimpleFoam.C
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/rhoPorousSimpleFoam.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
@@ -31,11 +31,12 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "rhoThermo.H"
+#include "fluidThermo.H"
 #include "turbulentFluidThermoModel.H"
+#include "simpleControl.H"
+#include "pressureControl.H"
 #include "fvOptions.H"
 #include "IOporosityModelList.H"
-#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C b/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C
index 4a3cdf8e01a..200e620eef4 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.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
@@ -30,9 +30,10 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "psiThermo.H"
+#include "fluidThermo.H"
 #include "turbulentFluidThermoModel.H"
 #include "simpleControl.H"
+#include "pressureControl.H"
 #include "fvOptions.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/etc/templates/compressibleInflowOutflow/system/fvSchemes b/etc/templates/compressibleInflowOutflow/system/fvSchemes
index 9986456cf79..ad7ca94ca6e 100644
--- a/etc/templates/compressibleInflowOutflow/system/fvSchemes
+++ b/etc/templates/compressibleInflowOutflow/system/fvSchemes
@@ -44,9 +44,9 @@ divSchemes
     div(phi,e)      $turbulence;
     div(phi,h)      $turbulence;
     div(phi,K)      $turbulence;
-    div(phi,Ekp)    $turbulence; 
+    div(phi,Ekp)    $turbulence;
 
-    div(phid,p)     bounded Gauss upwind;
+    div(phid,p)     Gauss upwind;
     div((phi|interpolate(rho)),p)  bounded Gauss upwind;
 
     div(((rho*nuEff)*dev2(T(grad(U)))))    Gauss linear;
diff --git a/etc/templates/compressibleInflowOutflow/system/fvSolution b/etc/templates/compressibleInflowOutflow/system/fvSolution
index 80579478d3b..d20f03588b6 100644
--- a/etc/templates/compressibleInflowOutflow/system/fvSolution
+++ b/etc/templates/compressibleInflowOutflow/system/fvSolution
@@ -43,8 +43,8 @@ SIMPLE
     }
 
     nNonOrthogonalCorrectors 0;
-    rhoMin          rhoMin [ 1 -3 0 0 0 ] 0.1;
-    rhoMax          rhoMax [ 1 -3 0 0 0 ] 1.5;
+    pMinFactor      0.1;
+    pMaxFactor      1.5;
 }
 
 relaxationFactors
@@ -56,7 +56,7 @@ relaxationFactors
     }
     equations
     {
-     	U               0.7;
+        U               0.7;
         "(e|h)"         0.7;
         "(k|epsilon|omega)" 0.7;
     }
diff --git a/src/TurbulenceModels/compressible/CompressibleTurbulenceModel/CompressibleTurbulenceModel.H b/src/TurbulenceModels/compressible/CompressibleTurbulenceModel/CompressibleTurbulenceModel.H
index 763303447ec..07b32950d78 100644
--- a/src/TurbulenceModels/compressible/CompressibleTurbulenceModel/CompressibleTurbulenceModel.H
+++ b/src/TurbulenceModels/compressible/CompressibleTurbulenceModel/CompressibleTurbulenceModel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -115,6 +115,19 @@ public:
             return this->transport_.mu(patchi);
         }
 
+        //- Return the laminar viscosity
+        virtual tmp<volScalarField> nu() const
+        {
+            return this->transport_.mu()/this->rho_;
+        }
+
+        //- Return the laminar viscosity on patchi
+        virtual tmp<scalarField> nu(const label patchi) const
+        {
+            return
+                this->transport_.mu(patchi)/this->rho_.boundaryField()[patchi];
+        }
+
         //- Return the turbulence dynamic viscosity
         virtual tmp<volScalarField> mut() const
         {
diff --git a/src/TurbulenceModels/turbulenceModels/TurbulenceModel/TurbulenceModel.H b/src/TurbulenceModels/turbulenceModels/TurbulenceModel/TurbulenceModel.H
index bff236d4864..61e0ac3d7f7 100644
--- a/src/TurbulenceModels/turbulenceModels/TurbulenceModel/TurbulenceModel.H
+++ b/src/TurbulenceModels/turbulenceModels/TurbulenceModel/TurbulenceModel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -158,13 +158,13 @@ public:
         }
 
         //- Return the laminar viscosity
-        tmp<volScalarField> nu() const
+        virtual tmp<volScalarField> nu() const
         {
             return transport_.nu();
         }
 
         //- Return the laminar viscosity on patchi
-        tmp<scalarField> nu(const label patchi) const
+        virtual tmp<scalarField> nu(const label patchi) const
         {
             return transport_.nu(patchi);
         }
diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files
index 2d7fbce1304..8f9122b58eb 100644
--- a/src/finiteVolume/Make/files
+++ b/src/finiteVolume/Make/files
@@ -402,6 +402,7 @@ $(general)/constrainHbyA/constrainHbyA.C
 $(general)/adjustPhi/adjustPhi.C
 $(general)/bound/bound.C
 $(general)/CorrectPhi/correctUphiBCs.C
+$(general)/pressureControl/pressureControl.C
 
 solutionControl = $(general)/solutionControl
 $(solutionControl)/solutionControl/solutionControl.C
diff --git a/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C b/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C
index 0d24e83a7bf..136fdda5959 100644
--- a/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C
+++ b/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.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
@@ -27,7 +27,7 @@ License
 
 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 
-void Foam::setRefCell
+bool Foam::setRefCell
 (
     const volScalarField& field,
     const volScalarField& fieldRef,
@@ -109,11 +109,17 @@ void Foam::setRefCell
         }
 
         refValue = readScalar(dict.lookup(refValueName));
+
+        return true;
+    }
+    else
+    {
+        return false;
     }
 }
 
 
-void Foam::setRefCell
+bool Foam::setRefCell
 (
     const volScalarField& field,
     const dictionary& dict,
@@ -122,7 +128,7 @@ void Foam::setRefCell
     const bool forceReference
 )
 {
-    setRefCell(field, field, dict, refCelli, refValue, forceReference);
+    return setRefCell(field, field, dict, refCelli, refValue, forceReference);
 }
 
 
diff --git a/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.H b/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.H
index 8a9d3d100a0..1532c718e78 100644
--- a/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.H
+++ b/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -46,8 +46,8 @@ namespace Foam
 
 //- If the field fieldRef needs referencing find the reference cell nearest
 //  (in index) to the given cell looked-up for field, but which is not on a
-//  cyclic, symmetry or processor patch.
-void setRefCell
+//  cyclic, symmetry or processor patch and return true, otherwise return false.
+bool setRefCell
 (
     const volScalarField& field,
     const volScalarField& fieldRef,
@@ -59,8 +59,8 @@ void setRefCell
 
 //- If the field needs referencing find the reference cell nearest
 //  (in index) to the given cell looked-up for field, but which is not on a
-//  cyclic, symmetry or processor patch.
-void setRefCell
+//  cyclic, symmetry or processor patch and return true, otherwise return false.
+bool setRefCell
 (
     const volScalarField& field,
     const dictionary& dict,
diff --git a/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.C b/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.C
new file mode 100644
index 00000000000..101126a5e87
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.C
@@ -0,0 +1,155 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "pressureControl.H"
+#include "findRefCell.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::pressureControl::pressureControl
+(
+    const volScalarField& p,
+    const volScalarField& rho,
+    const dictionary& dict
+)
+:
+    refCell_(0),
+    refValue_(0),
+    pMax_("pMax", dimPressure, 0),
+    pMin_("pMin", dimPressure, GREAT)
+{
+    if (setRefCell(p, dict, refCell_, refValue_))
+    {
+        pMax_.value() = refValue_;
+        pMin_.value() = refValue_;
+    }
+
+    const volScalarField::Boundary& pbf = p.boundaryField();
+    const volScalarField::Boundary& rhobf = rho.boundaryField();
+
+    scalar rhoRefMax = -GREAT;
+    scalar rhoRefMin = GREAT;
+    bool rhoLimits = false;
+
+    forAll(pbf, patchi)
+    {
+        if (pbf[patchi].fixesValue())
+        {
+            rhoLimits = true;
+
+            pMax_.value() = max(pMax_.value(), max(pbf[patchi]));
+            pMin_.value() = min(pMin_.value(), min(pbf[patchi]));
+
+            rhoRefMax = max(rhoRefMax, max(rhobf[patchi]));
+            rhoRefMin = min(rhoRefMin, min(rhobf[patchi]));
+        }
+    }
+
+    if (dict.found("pMax"))
+    {
+        pMax_.value() = readScalar(dict.lookup("pMax"));
+    }
+    else if (dict.found("pMaxFactor"))
+    {
+        const scalar pMaxFactor(readScalar(dict.lookup("pMaxFactor")));
+        pMax_ *= pMaxFactor;
+    }
+    else if (dict.found("rhoMax"))
+    {
+        // For backward-compatibility infer the pMax from rhoMax
+
+        IOWarningInFunction(dict)
+            << "'rhoMax' specified rather than 'pMax' or 'pMaxFactor'" << nl
+            << "    This is supported for backward-compatibility but "
+               "'pMax' or 'pMaxFactor' are more reliable." << endl;
+
+        if (!rhoLimits)
+        {
+            FatalIOErrorInFunction(dict)
+                << "'rhoMax' specified rather than 'pMaxFactor'" << nl
+                << "    but the corresponding reference density cannot"
+                   " be evaluated from the boundary conditions." << nl
+                << "Please specify 'pMaxFactor' rather than 'rhoMax'"
+                << exit(FatalError);
+        }
+
+        dimensionedScalar rhoMax("rhoMax", dimDensity, dict);
+
+        pMax_ *= max(rhoMax.value()/rhoRefMax, 1);
+    }
+
+    if (dict.found("pMin"))
+    {
+        pMin_.value() = readScalar(dict.lookup("pMin"));
+    }
+    else if (dict.found("pMinFactor"))
+    {
+        const scalar pMinFactor(readScalar(dict.lookup("pMinFactor")));
+        pMin_ *= pMinFactor;
+    }
+    else if (dict.found("rhoMin"))
+    {
+        // For backward-compatibility infer the pMin from rhoMin
+
+        IOWarningInFunction(dict)
+            << "'rhoMin' specified rather than 'pMin' or 'pMinFactor'" << nl
+            << "    This is supported for backward-compatibility but"
+               "'pMin' or 'pMinFactor' are more reliable." << endl;
+
+        if (!rhoLimits)
+        {
+            FatalIOErrorInFunction(dict)
+                << "'rhoMin' specified rather than 'pMinFactor'" << nl
+                << "    but the corresponding reference density cannot"
+                   " be evaluated from the boundary conditions." << nl
+                << "Please specify 'pMinFactor' rather than 'rhoMin'"
+                << exit(FatalError);
+        }
+
+        dimensionedScalar rhoMin("rhoMin", dimDensity, dict);
+
+        pMin_ *= min(rhoMin.value()/rhoRefMin, 1);
+    }
+
+    Info<< "pressureControl" << nl
+        << "    pMax/pMin " << pMax_.value() << " " << pMin_.value()
+        << nl << endl;
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::pressureControl::limit(volScalarField& p) const
+{
+    Info<< "pressureControl: p max/min "
+        << max(p).value() << " "
+        << min(p).value() << endl;
+
+    p = max(p, pMin_);
+    p = min(p, pMax_);
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.H b/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.H
new file mode 100644
index 00000000000..19bd053d430
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.H
@@ -0,0 +1,108 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::pressureControl
+
+Description
+    Provides controls for the pressure reference is closed-volume simulations
+    and a general method for limiting the pressure during the startup of
+    steady-state simulations.
+
+SourceFiles
+    pressureControlI.H
+    pressureControl.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef pressureControl_H
+#define pressureControl_H
+
+#include "dimensionedScalar.H"
+#include "volFieldsFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class pressureControl Declaration
+\*---------------------------------------------------------------------------*/
+
+class pressureControl
+{
+    // Private data
+
+        //- Optional cell in which the reference pressure is set
+        label refCell_;
+
+        //- Optional pressure reference level
+        scalar refValue_;
+
+        //- Pressure lower-limit
+        dimensionedScalar pMax_;
+
+        //- Pressure upper-limit
+        dimensionedScalar pMin_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from the SIMPLE/PIMPLE sub-dictionary
+        pressureControl
+        (
+            const volScalarField& p,
+            const volScalarField& rho,
+            const dictionary& dict
+        );
+
+
+    // Member Functions
+
+        //- Return the cell in which the reference pressure is set
+        inline label refCell() const;
+
+        //- Return the pressure reference level
+        inline scalar refValue() const;
+
+        //- Limit the pressure
+        void limit(volScalarField& p) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "pressureControlI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/pressureControl/pressureControlI.H b/src/finiteVolume/cfdTools/general/pressureControl/pressureControlI.H
new file mode 100644
index 00000000000..6548d43455d
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/pressureControl/pressureControlI.H
@@ -0,0 +1,40 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline Foam::label Foam::pressureControl::refCell() const
+{
+    return refCell_;
+}
+
+
+inline Foam::scalar Foam::pressureControl::refValue() const
+{
+    return refValue_;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctExplicit/system/fvSolution b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctExplicit/system/fvSolution
index 9eee4c23fb0..4b7b15b9276 100644
--- a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctExplicit/system/fvSolution
+++ b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctExplicit/system/fvSolution
@@ -54,8 +54,8 @@ solvers
 SIMPLE
 {
     nNonOrthogonalCorrectors 0;
-    rhoMin          0.5;
-    rhoMax          1.5;
+    pMinFactor      0.5;
+    pMaxFactor      1.5;
 
     residualControl
     {
diff --git a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/system/fvSolution b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/system/fvSolution
index 8699fd70954..421197bf139 100644
--- a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/system/fvSolution
+++ b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/system/fvSolution
@@ -48,8 +48,8 @@ SIMPLE
 {
     nUCorrectors        2;
     nNonOrthogonalCorrectors 0;
-    rhoMin              0.5;
-    rhoMax              2.0;
+    pMinFactor          0.5;
+    pMaxFactor          2.0;
 
     residualControl
     {
diff --git a/tutorials/compressible/rhoSimpleFoam/angledDuctExplicitFixedCoeff/system/fvSolution b/tutorials/compressible/rhoSimpleFoam/angledDuctExplicitFixedCoeff/system/fvSolution
index a054853c3ee..25619f1f3ff 100644
--- a/tutorials/compressible/rhoSimpleFoam/angledDuctExplicitFixedCoeff/system/fvSolution
+++ b/tutorials/compressible/rhoSimpleFoam/angledDuctExplicitFixedCoeff/system/fvSolution
@@ -54,8 +54,8 @@ solvers
 SIMPLE
 {
     nNonOrthogonalCorrectors 0;
-    rhoMin          0.4;
-    rhoMax          1.5;
+    pMinFactor      0.4;
+    pMaxFactor      1.5;
 
     residualControl
     {
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBend/system/fvSolution b/tutorials/compressible/rhoSimpleFoam/squareBend/system/fvSolution
index 3855d4cb054..2ffe1fe4be9 100644
--- a/tutorials/compressible/rhoSimpleFoam/squareBend/system/fvSolution
+++ b/tutorials/compressible/rhoSimpleFoam/squareBend/system/fvSolution
@@ -39,8 +39,8 @@ solvers
 SIMPLE
 {
     nNonOrthogonalCorrectors 0;
-    rhoMin          0.1;
-    rhoMax          1.0;
+    pMinFactor      0.1;
+    pMaxFactor      1.5;
     transonic       yes;
     consistent      yes;
 
@@ -60,7 +60,6 @@ relaxationFactors
     fields
     {
         p               1;
-        rho             1;
     }
     equations
     {
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/T b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/T
new file mode 100644
index 00000000000..2df93ba37ed
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/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;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/U b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/U
new file mode 100644
index 00000000000..99e3c1d6354
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/U
@@ -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       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            inletOutlet;
+        value           uniform (0 0 0);
+        inletValue      uniform (0 0 0);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/alphat b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/alphat
new file mode 100644
index 00000000000..3a6625042cc
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/alphat
@@ -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       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;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/epsilon b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/epsilon
new file mode 100644
index 00000000000..71b9c782d16
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/epsilon
@@ -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      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;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/k b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/k
new file mode 100644
index 00000000000..7ee52f8dbfb
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/k
@@ -0,0 +1,44 @@
+/*--------------------------------*- 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;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/nut b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/nut
new file mode 100644
index 00000000000..ac8d562f40e
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/nut
@@ -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      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;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/p b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/p
new file mode 100644
index 00000000000..5e0d23f65fd
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/p
@@ -0,0 +1,39 @@
+/*--------------------------------*- 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;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/constant/thermophysicalProperties b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/constant/thermophysicalProperties
new file mode 100644
index 00000000000..426ab721160
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/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;
+}
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/constant/turbulenceProperties b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/constant/turbulenceProperties
new file mode 100644
index 00000000000..cd2daf8229b
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/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;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/blockMeshDict b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/blockMeshDict
new file mode 100644
index 00000000000..9bec6a9547f
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/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
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/controlDict b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/controlDict
new file mode 100644
index 00000000000..edddd719e6e
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/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     rhoSimpleFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         500;
+
+deltaT          1;
+
+writeControl    timeStep;
+
+writeInterval   100;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+graphFormat     raw;
+
+runTimeModifiable true;
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/decomposeParDict b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/decomposeParDict
new file mode 100644
index 00000000000..9e844a62ba7
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/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           ( );
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/fvSchemes b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/fvSchemes
new file mode 100644
index 00000000000..c68f548c95f
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/fvSchemes
@@ -0,0 +1,59 @@
+/*--------------------------------*- 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             steadyState;
+}
+
+gradSchemes
+{
+    default             Gauss linear;
+
+    limited             cellLimited Gauss linear 1;
+}
+
+divSchemes
+{
+    default             none;
+
+    div(phi,U)          bounded Gauss linearUpwind limited;
+    div(phi,e)          bounded Gauss linearUpwind limited;
+    div(phi,epsilon)    bounded Gauss linearUpwind limited;
+    div(phi,k)          bounded Gauss linearUpwind limited;
+    div(phi,Ekp)        bounded Gauss linearUpwind limited;
+
+    div(((rho*nuEff)*dev2(T(grad(U)))))      Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/fvSolution b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/fvSolution
new file mode 100644
index 00000000000..3b30d4ebe5b
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/fvSolution
@@ -0,0 +1,71 @@
+/*--------------------------------*- 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
+{
+    p
+    {
+        solver          GAMG;
+        smoother        GaussSeidel;
+
+        tolerance       1e-7;
+        relTol          0.1;
+    }
+
+    "(U|e|k|epsilon)"
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+
+        tolerance       1e-7;
+        relTol          0.1;
+    }
+}
+
+SIMPLE
+{
+    nNonOrthogonalCorrectors 0;
+    pMinFactor      0.1;
+    pMaxFactor      1.5;
+
+    transonic       no;
+    consistent      no;
+
+    residualControl
+    {
+        p               1e-4;
+        U               1e-4;
+        "(k|omega|epsilon|e|h)" 1e-3;
+    }
+}
+
+relaxationFactors
+{
+    fields
+    {
+        p               0.3;
+    }
+    equations
+    {
+        U               0.7;
+        e               0.7;
+        k               0.7;
+        epsilon         0.7;
+    }
+}
+
+// ************************************************************************* //
-- 
GitLab