From 06eb5b74777ead763ddffa52064395ad775e2462 Mon Sep 17 00:00:00 2001
From: henry <Henry Weller h.weller@opencfd.co.uk>
Date: Wed, 25 Nov 2009 11:59:19 +0000
Subject: [PATCH] Added preliminary version of the rhoSimplecFoam code and
 tutorial case. This currently only runs in serial, parallelisation is worked
 on.

---
 .../solvers/compressible/rhoSimpleFoam/pEqn.H |   2 +-
 .../compressible/rhoSimplecFoam/Make/files    |   3 +
 .../compressible/rhoSimplecFoam/Make/options  |  14 ++
 .../compressible/rhoSimplecFoam/UEqn.H        |  17 ++
 .../rhoSimplecFoam/createFields.H             |  63 ++++++
 .../compressible/rhoSimplecFoam/hEqn.H        |  29 +++
 .../compressible/rhoSimplecFoam/pEqn.H        | 123 +++++++++++
 .../rhoSimplecFoam/rhoSimplecFoam.C           |  92 ++++++++
 .../rhoSimplecFoam/squareBend/0/T             |  43 ++++
 .../rhoSimplecFoam/squareBend/0/U             |  43 ++++
 .../rhoSimplecFoam/squareBend/0/alphat        |  43 ++++
 .../rhoSimplecFoam/squareBend/0/epsilon       |  47 +++++
 .../rhoSimplecFoam/squareBend/0/k             |  44 ++++
 .../rhoSimplecFoam/squareBend/0/mut           |  45 ++++
 .../rhoSimplecFoam/squareBend/0/p             |  54 +++++
 .../squareBend/constant/RASProperties         |  25 +++
 .../constant/polyMesh/blockMeshDict           | 115 ++++++++++
 .../squareBend/constant/polyMesh/boundary     |  40 ++++
 .../constant/thermophysicalProperties         |  23 ++
 .../squareBend/constant/turbulenceProperties  |  21 ++
 .../squareBend/system/controlDict             |  51 +++++
 .../squareBend/system/fvSchemes               |  71 +++++++
 .../squareBend/system/fvSolution              | 199 ++++++++++++++++++
 23 files changed, 1206 insertions(+), 1 deletion(-)
 create mode 100644 applications/solvers/compressible/rhoSimplecFoam/Make/files
 create mode 100644 applications/solvers/compressible/rhoSimplecFoam/Make/options
 create mode 100644 applications/solvers/compressible/rhoSimplecFoam/UEqn.H
 create mode 100644 applications/solvers/compressible/rhoSimplecFoam/createFields.H
 create mode 100644 applications/solvers/compressible/rhoSimplecFoam/hEqn.H
 create mode 100644 applications/solvers/compressible/rhoSimplecFoam/pEqn.H
 create mode 100644 applications/solvers/compressible/rhoSimplecFoam/rhoSimplecFoam.C
 create mode 100644 tutorials/compressible/rhoSimplecFoam/squareBend/0/T
 create mode 100644 tutorials/compressible/rhoSimplecFoam/squareBend/0/U
 create mode 100644 tutorials/compressible/rhoSimplecFoam/squareBend/0/alphat
 create mode 100644 tutorials/compressible/rhoSimplecFoam/squareBend/0/epsilon
 create mode 100644 tutorials/compressible/rhoSimplecFoam/squareBend/0/k
 create mode 100644 tutorials/compressible/rhoSimplecFoam/squareBend/0/mut
 create mode 100644 tutorials/compressible/rhoSimplecFoam/squareBend/0/p
 create mode 100644 tutorials/compressible/rhoSimplecFoam/squareBend/constant/RASProperties
 create mode 100644 tutorials/compressible/rhoSimplecFoam/squareBend/constant/polyMesh/blockMeshDict
 create mode 100644 tutorials/compressible/rhoSimplecFoam/squareBend/constant/polyMesh/boundary
 create mode 100644 tutorials/compressible/rhoSimplecFoam/squareBend/constant/thermophysicalProperties
 create mode 100644 tutorials/compressible/rhoSimplecFoam/squareBend/constant/turbulenceProperties
 create mode 100644 tutorials/compressible/rhoSimplecFoam/squareBend/system/controlDict
 create mode 100644 tutorials/compressible/rhoSimplecFoam/squareBend/system/fvSchemes
 create mode 100644 tutorials/compressible/rhoSimplecFoam/squareBend/system/fvSolution

diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
index f6a433fd616..38922c99c9c 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
@@ -27,7 +27,7 @@ if (transonic)
 
         pEqn.setReference(pRefCell, pRefValue);
 
-        // retain the residual from the first iteration
+        // Retain the residual from the first iteration
         if (nonOrth == 0)
         {
             eqnResidual = pEqn.solve().initialResidual();
diff --git a/applications/solvers/compressible/rhoSimplecFoam/Make/files b/applications/solvers/compressible/rhoSimplecFoam/Make/files
new file mode 100644
index 00000000000..6637e49aa3b
--- /dev/null
+++ b/applications/solvers/compressible/rhoSimplecFoam/Make/files
@@ -0,0 +1,3 @@
+rhoSimplecFoam.C
+
+EXE = $(FOAM_APPBIN)/rhoSimplecFoam
diff --git a/applications/solvers/compressible/rhoSimplecFoam/Make/options b/applications/solvers/compressible/rhoSimplecFoam/Make/options
new file mode 100644
index 00000000000..9d578f011a7
--- /dev/null
+++ b/applications/solvers/compressible/rhoSimplecFoam/Make/options
@@ -0,0 +1,14 @@
+EXE_INC = \
+    -I../rhoSimpleFoam \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels \
+    -I$(LIB_SRC)/turbulenceModels/compressible/RAS/RASModel \
+    -I$(LIB_SRC)/finiteVolume/cfdTools \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -lbasicThermophysicalModels \
+    -lspecie \
+    -lcompressibleRASModels \
+    -lfiniteVolume \
+    -lmeshTools
diff --git a/applications/solvers/compressible/rhoSimplecFoam/UEqn.H b/applications/solvers/compressible/rhoSimplecFoam/UEqn.H
new file mode 100644
index 00000000000..c41bc9b6c7b
--- /dev/null
+++ b/applications/solvers/compressible/rhoSimplecFoam/UEqn.H
@@ -0,0 +1,17 @@
+    // Solve the Momentum equation
+
+    tmp<fvVectorMatrix> UEqn
+    (
+        fvm::div(phi, U)
+      - fvm::Sp(fvc::div(phi), U)
+      + turbulence->divDevRhoReff(U)
+    );
+
+    UEqn().relax();
+
+    eqnResidual = solve
+    (
+        UEqn() == -fvc::grad(p)
+    ).initialResidual();
+
+    maxResidual = max(eqnResidual, maxResidual);
diff --git a/applications/solvers/compressible/rhoSimplecFoam/createFields.H b/applications/solvers/compressible/rhoSimplecFoam/createFields.H
new file mode 100644
index 00000000000..d97ee4705b0
--- /dev/null
+++ b/applications/solvers/compressible/rhoSimplecFoam/createFields.H
@@ -0,0 +1,63 @@
+    Info<< "Reading thermophysical properties\n" << endl;
+
+    autoPtr<basicPsiThermo> pThermo
+    (
+        basicPsiThermo::New(mesh)
+    );
+    basicPsiThermo& thermo = pThermo();
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        thermo.rho()
+    );
+
+    volScalarField& p = thermo.p();
+    volScalarField& h = thermo.h();
+    const volScalarField& psi = thermo.psi();
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    #include "compressibleCreatePhi.H"
+
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);
+
+    dimensionedScalar pMin
+    (
+        mesh.solutionDict().subDict("SIMPLE").lookup("pMin")
+    );
+
+    Info<< "Creating turbulence model\n" << endl;
+    autoPtr<compressible::RASModel> turbulence
+    (
+        compressible::RASModel::New
+        (
+            rho,
+            U,
+            phi,
+            thermo
+        )
+    );
+
+    dimensionedScalar initialMass = fvc::domainIntegrate(rho);
diff --git a/applications/solvers/compressible/rhoSimplecFoam/hEqn.H b/applications/solvers/compressible/rhoSimplecFoam/hEqn.H
new file mode 100644
index 00000000000..8ac1c9f5109
--- /dev/null
+++ b/applications/solvers/compressible/rhoSimplecFoam/hEqn.H
@@ -0,0 +1,29 @@
+{
+    fvScalarMatrix hEqn
+    (
+        fvm::div(phi, h)
+      - fvm::Sp(fvc::div(phi), h)
+      - fvm::laplacian(turbulence->alphaEff(), h)
+     ==
+        fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p, "div(U,p)"))
+      - p*fvc::div(phi/fvc::interpolate(rho))
+    );
+
+    hEqn.relax();
+
+    eqnResidual = hEqn.solve().initialResidual();
+    maxResidual = max(eqnResidual, maxResidual);
+
+    thermo.correct();
+
+    rho = thermo.rho();
+
+    if (!transonic)
+    {
+        rho.relax();
+    }
+
+    Info<< "rho max/min : "
+        << max(rho).value() << " "
+        << min(rho).value() << endl;
+}
diff --git a/applications/solvers/compressible/rhoSimplecFoam/pEqn.H b/applications/solvers/compressible/rhoSimplecFoam/pEqn.H
new file mode 100644
index 00000000000..cd3ae5ff1f4
--- /dev/null
+++ b/applications/solvers/compressible/rhoSimplecFoam/pEqn.H
@@ -0,0 +1,123 @@
+volScalarField p0 = p;
+
+volScalarField AU = UEqn().A();
+volScalarField AtU = AU - UEqn().H1();
+U = UEqn().H()/AU;
+UEqn.clear();
+
+bool closedVolume = false;
+
+if (transonic)
+{
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        surfaceScalarField phid
+        (
+            "phid",
+            fvc::interpolate(psi*U) & mesh.Sf()
+        );
+
+        surfaceScalarField phic
+        (
+            "phic",
+            fvc::interpolate(rho/AtU - rho/AU)*fvc::snGrad(p)*mesh.magSf()
+          + phid*(fvc::interpolate(p) - fvc::interpolate(p, "UD"))
+        );
+
+        refCast<mixedFvPatchScalarField>(p.boundaryField()[1]).refValue()
+            = p.boundaryField()[1];
+
+        fvScalarMatrix pEqn
+        (
+            fvm::div(phid, p)
+          + fvc::div(phic)
+          - fvm::Sp(fvc::div(phid), p)
+          + fvc::div(phid)*p
+          - fvm::laplacian(rho/AtU, p)
+        );
+
+        pEqn.setReference(pRefCell, pRefValue);
+
+        // Retain the residual from the first iteration
+        if (nonOrth == 0)
+        {
+            eqnResidual = pEqn.solve().initialResidual();
+            maxResidual = max(eqnResidual, maxResidual);
+        }
+        else
+        {
+            pEqn.solve();
+        }
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi == phic + pEqn.flux();
+        }
+    }
+}
+else
+{
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        phi = fvc::interpolate(rho*U) & mesh.Sf();
+        closedVolume = adjustPhi(phi, U, p);
+        phi += fvc::interpolate(rho/AtU - rho/AU)*fvc::snGrad(p)*mesh.magSf();
+
+        fvScalarMatrix pEqn
+        (
+            fvc::div(phi)
+        //- fvm::laplacian(rho/AU, p)
+          - fvm::laplacian(rho/AtU, p)
+        );
+
+        pEqn.setReference(pRefCell, pRefValue);
+
+        // Retain the residual from the first iteration
+        if (nonOrth == 0)
+        {
+            eqnResidual = pEqn.solve().initialResidual();
+            maxResidual = max(eqnResidual, maxResidual);
+        }
+        else
+        {
+            pEqn.solve();
+        }
+
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi += pEqn.flux();
+        }
+    }
+}
+
+// The incompressibe for of the continuity error check is appropriate for
+// steady-state compressible also.
+#include "incompressible/continuityErrs.H"
+
+// Explicitly relax pressure for momentum corrector
+p.relax();
+
+U -= (fvc::grad(p0)*(1.0/AU - 1.0/AtU) + fvc::grad(p)/AtU);
+//U -= fvc::grad(p)/AU;
+
+U.correctBoundaryConditions();
+
+bound(p, pMin);
+
+// For closed-volume cases adjust the pressure and density levels
+// to obey overall mass continuity
+if (closedVolume)
+{
+    p += (initialMass - fvc::domainIntegrate(psi*p))
+        /fvc::domainIntegrate(psi);
+}
+
+rho = thermo.rho();
+
+if (!transonic)
+{
+    rho.relax();
+}
+
+Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
diff --git a/applications/solvers/compressible/rhoSimplecFoam/rhoSimplecFoam.C b/applications/solvers/compressible/rhoSimplecFoam/rhoSimplecFoam.C
new file mode 100644
index 00000000000..8f706182b56
--- /dev/null
+++ b/applications/solvers/compressible/rhoSimplecFoam/rhoSimplecFoam.C
@@ -0,0 +1,92 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    rhoSimplecFoam
+
+Description
+    Steady-state SIMPLEC solver for laminar or turbulent RANS flow of
+    compressible fluids.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "basicPsiThermo.H"
+#include "RASModel.H"
+#include "mixedFvPatchFields.H"
+#include "bound.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    for (runTime++; !runTime.end(); runTime++)
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        #include "readSIMPLEControls.H"
+        #include "initConvergenceCheck.H"
+
+        p.storePrevIter();
+
+        if (!transonic)
+        {
+            rho.storePrevIter();
+        }
+
+        // Velocity-pressure-enthalpy SIMPLEC corrector
+        {
+            #include "UEqn.H"
+            #include "pEqn.H"
+            #include "hEqn.H"
+        }
+
+        turbulence->correct();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+
+        #include "convergenceCheck.H"
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimplecFoam/squareBend/0/T b/tutorials/compressible/rhoSimplecFoam/squareBend/0/T
new file mode 100644
index 00000000000..c8a138aa587
--- /dev/null
+++ b/tutorials/compressible/rhoSimplecFoam/squareBend/0/T
@@ -0,0 +1,43 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://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 1000;
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            zeroGradient;
+    }
+
+    inlet
+    {
+        type		fixedValue;
+        value		uniform 1000;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        //type            zeroGradient;
+        value           uniform 1000;
+        inletValue      uniform 1000;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimplecFoam/squareBend/0/U b/tutorials/compressible/rhoSimplecFoam/squareBend/0/U
new file mode 100644
index 00000000000..971e760d3de
--- /dev/null
+++ b/tutorials/compressible/rhoSimplecFoam/squareBend/0/U
@@ -0,0 +1,43 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://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            fixedValue;
+        value           uniform (0 0 0);
+    }
+    inlet
+    {
+        type            flowRateInletVelocity;
+        flowRate        0.5; //0.75;
+        value           uniform (0 0 0);
+    }
+    outlet
+    {
+        type            inletOutlet;
+        value           uniform (0 0 0);
+        inletValue      uniform (0 0 0);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimplecFoam/squareBend/0/alphat b/tutorials/compressible/rhoSimplecFoam/squareBend/0/alphat
new file mode 100644
index 00000000000..7b530fa40d6
--- /dev/null
+++ b/tutorials/compressible/rhoSimplecFoam/squareBend/0/alphat
@@ -0,0 +1,43 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    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            alphatWallFunction;
+        Prt             0.85;
+        value           uniform 0;
+    }
+    inlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    outlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimplecFoam/squareBend/0/epsilon b/tutorials/compressible/rhoSimplecFoam/squareBend/0/epsilon
new file mode 100644
index 00000000000..c63316a51f1
--- /dev/null
+++ b/tutorials/compressible/rhoSimplecFoam/squareBend/0/epsilon
@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    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            compressible::epsilonWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 200;
+    }
+    inlet
+    {
+        type            compressible::turbulentMixingLengthDissipationRateInlet;
+        mixingLength    0.005;
+        value           uniform 200;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 200;
+        value           uniform 200;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimplecFoam/squareBend/0/k b/tutorials/compressible/rhoSimplecFoam/squareBend/0/k
new file mode 100644
index 00000000000..b885efd9029
--- /dev/null
+++ b/tutorials/compressible/rhoSimplecFoam/squareBend/0/k
@@ -0,0 +1,44 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    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            compressible::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/rhoSimplecFoam/squareBend/0/mut b/tutorials/compressible/rhoSimplecFoam/squareBend/0/mut
new file mode 100644
index 00000000000..4d70cca517b
--- /dev/null
+++ b/tutorials/compressible/rhoSimplecFoam/squareBend/0/mut
@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      mut;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            mutkWallFunction;
+        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/rhoSimplecFoam/squareBend/0/p b/tutorials/compressible/rhoSimplecFoam/squareBend/0/p
new file mode 100644
index 00000000000..fa45b63e2ce
--- /dev/null
+++ b/tutorials/compressible/rhoSimplecFoam/squareBend/0/p
@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://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 110000;
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            zeroGradient;
+    }
+    inlet
+    {
+        //type            zeroGradient;
+        type            mixed;
+        refValue        uniform 110000;
+        refGradient     uniform 0;
+        valueFraction   uniform 0.3;
+    }
+    outlet
+    {
+        type            fixedValue;
+        value           uniform 110000;
+
+        //type            mixed;
+        //refValue        uniform 110000;
+        //refGradient     uniform 0;
+        //valueFraction   uniform 1;
+        //type            transonicOutletPressure;
+        //U               U;
+        //phi             phi;
+        //gamma           1.4;
+        //psi             psi;
+        //pInf            uniform 110000;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimplecFoam/squareBend/constant/RASProperties b/tutorials/compressible/rhoSimplecFoam/squareBend/constant/RASProperties
new file mode 100644
index 00000000000..81b1ec91152
--- /dev/null
+++ b/tutorials/compressible/rhoSimplecFoam/squareBend/constant/RASProperties
@@ -0,0 +1,25 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      RASProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+RASModel        kEpsilon;
+
+turbulence      on;
+
+printCoeffs     on;
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimplecFoam/squareBend/constant/polyMesh/blockMeshDict b/tutorials/compressible/rhoSimplecFoam/squareBend/constant/polyMesh/blockMeshDict
new file mode 100644
index 00000000000..f8b06e3446c
--- /dev/null
+++ b/tutorials/compressible/rhoSimplecFoam/squareBend/constant/polyMesh/blockMeshDict
@@ -0,0 +1,115 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://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)
+);
+
+patches
+(
+    // is there no way of defining all my 'defaultFaces' to be 'wall'?
+    wall Default_Boundary_Region
+    (
+        // 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 )
+    )
+    patch inlet
+    (
+        (0 2 12 10)
+    )
+
+    patch outlet
+    (
+        (4 6 16 14)
+    )
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimplecFoam/squareBend/constant/polyMesh/boundary b/tutorials/compressible/rhoSimplecFoam/squareBend/constant/polyMesh/boundary
new file mode 100644
index 00000000000..ae1e146725b
--- /dev/null
+++ b/tutorials/compressible/rhoSimplecFoam/squareBend/constant/polyMesh/boundary
@@ -0,0 +1,40 @@
+/*--------------------------------*- 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       polyBoundaryMesh;
+    location    "constant/polyMesh";
+    object      boundary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+3
+(
+    Default_Boundary_Region
+    {
+        type            wall;
+        nFaces          22400;
+        startFace       324400;
+    }
+    inlet
+    {
+        type            patch;
+        nFaces          400;
+        startFace       346800;
+    }
+    outlet
+    {
+        type            patch;
+        nFaces          400;
+        startFace       347200;
+    }
+)
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimplecFoam/squareBend/constant/thermophysicalProperties b/tutorials/compressible/rhoSimplecFoam/squareBend/constant/thermophysicalProperties
new file mode 100644
index 00000000000..a2c84a89d38
--- /dev/null
+++ b/tutorials/compressible/rhoSimplecFoam/squareBend/constant/thermophysicalProperties
@@ -0,0 +1,23 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType      hPsiThermo<pureMixture<sutherlandTransport<specieThermo<hConstThermo<perfectGas>>>>>;
+
+mixture         air 1 28.9 1007 0 1.4792e-06 116;
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimplecFoam/squareBend/constant/turbulenceProperties b/tutorials/compressible/rhoSimplecFoam/squareBend/constant/turbulenceProperties
new file mode 100644
index 00000000000..f6753662e27
--- /dev/null
+++ b/tutorials/compressible/rhoSimplecFoam/squareBend/constant/turbulenceProperties
@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  RASModel;
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimplecFoam/squareBend/system/controlDict b/tutorials/compressible/rhoSimplecFoam/squareBend/system/controlDict
new file mode 100644
index 00000000000..49ae3f8a17d
--- /dev/null
+++ b/tutorials/compressible/rhoSimplecFoam/squareBend/system/controlDict
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     rhoSimplecFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         500;
+
+deltaT          1;
+
+writeControl    timeStep;
+
+writeInterval   100;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+graphFormat     raw;
+
+runTimeModifiable yes;
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimplecFoam/squareBend/system/fvSchemes b/tutorials/compressible/rhoSimplecFoam/squareBend/system/fvSchemes
new file mode 100644
index 00000000000..e84e32f437f
--- /dev/null
+++ b/tutorials/compressible/rhoSimplecFoam/squareBend/system/fvSchemes
@@ -0,0 +1,71 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    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;
+}
+
+divSchemes
+{
+    default             none;
+
+    div(phi,U)          Gauss upwind; //limitedLinearV 1; //upwind;
+    div((muEff*dev2(grad(U).T())))      Gauss linear;
+    div(phi,h)          Gauss upwind; //limitedLinear 1; //upwind;
+    div(phi,epsilon)    Gauss upwind;
+    div(phi,k)          Gauss upwind;
+
+    div(phid,p)         Gauss upwind;
+    div(U,p)            Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default                 none;
+    interpolate(rho)        linear;
+    div(U,p)                upwind phi;
+    interpolate((psi*U))    linear;
+    interpolate(U)          linear;
+    UD                      upwind phid;
+    interpolate(p)          linear;
+    interpolate(((rho|(A(U)-H(1)))-(rho|A(U)))) linear;
+}
+
+snGradSchemes
+{
+    default                 corrected;
+}
+
+fluxRequired
+{
+    default                 no;
+    p;
+    pCorr;
+}
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimplecFoam/squareBend/system/fvSolution b/tutorials/compressible/rhoSimplecFoam/squareBend/system/fvSolution
new file mode 100644
index 00000000000..d114a855712
--- /dev/null
+++ b/tutorials/compressible/rhoSimplecFoam/squareBend/system/fvSolution
@@ -0,0 +1,199 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p0
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-08;
+        relTol          0.01;
+    }
+
+    p
+    {
+        solver          GAMG;
+        tolerance       1e-08;
+        relTol          0.1;
+        smoother        GaussSeidel;
+        nPreSweeps      0;
+        nPostSweeps     2;
+        nFinestSweeps   2;
+        cacheAgglomeration true;
+        nCellsInCoarsestLevel 20;
+        agglomerator    faceAreaPair;
+        mergeLevels     1;
+    }
+
+    U0
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-08;
+        relTol          0.1;
+    }
+
+    U1
+    {
+        solver          smoothSolver;
+        smoother        GaussSeidel;
+        nSweeps         1;
+        tolerance       1e-08;
+        relTol          0.1;
+    }
+
+    U
+    {
+        solver          GAMG;
+        tolerance       1e-08;
+        relTol          0.1;
+        smoother        GaussSeidel;
+        nPreSweeps      0;
+        nPostSweeps     2;
+        nFinestSweeps   2;
+        cacheAgglomeration true;
+        nCellsInCoarsestLevel 20;
+        agglomerator    faceAreaPair;
+        mergeLevels     1;
+    }
+
+    h0
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-08;
+        relTol          0.1;
+    }
+
+    h1
+    {
+        solver          smoothSolver;
+        smoother        GaussSeidel;
+        nSweeps         1;
+        tolerance       1e-08;
+        relTol          0.1;
+    }
+
+    h
+    {
+        solver          GAMG;
+        tolerance       1e-08;
+        relTol          0.1;
+        smoother        GaussSeidel;
+        nPreSweeps      0;
+        nPostSweeps     2;
+        nFinestSweeps   2;
+        cacheAgglomeration true;
+        nCellsInCoarsestLevel 20;
+        agglomerator    faceAreaPair;
+        mergeLevels     1;
+    }
+
+    k0
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-08;
+        relTol          0.1;
+    }
+
+    k1
+    {
+        solver          smoothSolver;
+        smoother        GaussSeidel;
+        nSweeps         1;
+        tolerance       1e-08;
+        relTol          0.1;
+    }
+
+    k
+    {
+        solver          GAMG;
+        tolerance       1e-08;
+        relTol          0.1;
+        smoother        GaussSeidel;
+        nPreSweeps      0;
+        nPostSweeps     2;
+        nFinestSweeps   2;
+        cacheAgglomeration true;
+        nCellsInCoarsestLevel 20;
+        agglomerator    faceAreaPair;
+        mergeLevels     1;
+    }
+
+    epsilon0
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-08;
+        relTol          0.1;
+    }
+
+    epsilon1
+    {
+        solver          smoothSolver;
+        smoother        GaussSeidel;
+        nSweeps         1;
+        tolerance       1e-08;
+        relTol          0.1;
+    }
+
+    epsilon
+    {
+        solver          GAMG;
+        tolerance       1e-08;
+        relTol          0.1;
+        smoother        GaussSeidel;
+        nPreSweeps      0;
+        nPostSweeps     2;
+        nFinestSweeps   2;
+        cacheAgglomeration true;
+        nCellsInCoarsestLevel 20;
+        agglomerator    faceAreaPair;
+        mergeLevels     1;
+    }
+}
+
+SIMPLE
+{
+    nNonOrthogonalCorrectors        0;
+    pMin                            pMin [1 -1 -2 0 0 0 0] 10000;
+    transonic                       yes;
+}
+
+relaxationFactors
+{
+    p           1;
+    rho         1; //0.1;
+    U           0.9;
+    h           0.95;
+    k           0.9;
+    epsilon     0.9;
+}
+
+relaxationFactors0
+{
+    p           0.3;
+    rho         0.1;
+    U           0.7;
+    h           0.7;
+    k           0.7;
+    epsilon     0.7;
+}
+
+// ************************************************************************* //
-- 
GitLab