diff --git a/applications/solvers/combustion/reactingFoam/createFields.H b/applications/solvers/combustion/reactingFoam/createFields.H
index 6fcf83624c621c96b621d78bb06f8597ff86d165..bc093581930e6d3a6ba3109a8e95d8709eb2cdef 100644
--- a/applications/solvers/combustion/reactingFoam/createFields.H
+++ b/applications/solvers/combustion/reactingFoam/createFields.H
@@ -52,27 +52,7 @@ volScalarField& p = thermo.p();
#include "compressibleCreatePhi.H"
-dimensionedScalar rhoMax
-(
- dimensionedScalar::lookupOrDefault
- (
- "rhoMax",
- pimple.dict(),
- dimDensity,
- GREAT
- )
-);
-
-dimensionedScalar rhoMin
-(
- dimensionedScalar::lookupOrDefault
- (
- "rhoMin",
- pimple.dict(),
- dimDensity,
- 0
- )
-);
+pressureControl pressureControl(p, rho, pimple.dict(), false);
mesh.setFluxRequired(p.name());
diff --git a/applications/solvers/combustion/reactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/pEqn.H
index ac7107acf066722b286275569b0a8b2b39af3e05..81ab6424b4c1a7bfac43df5538ff4d089da5dca0 100644
--- a/applications/solvers/combustion/reactingFoam/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/pEqn.H
@@ -1,7 +1,4 @@
rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
@@ -87,19 +84,17 @@ else
// Explicitly relax pressure for momentum corrector
p.relax();
-// Recalculate density from the relaxed pressure
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
-Info<< "rho max/min : " << max(rho).value()
- << " " << min(rho).value() << endl;
-
U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();
fvOptions.correct(U);
K = 0.5*magSqr(U);
+if (pressureControl.limit(p))
+{
+ p.correctBoundaryConditions();
+ rho = thermo.rho();
+}
+
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);
diff --git a/applications/solvers/combustion/reactingFoam/pcEqn.H b/applications/solvers/combustion/reactingFoam/pcEqn.H
index 713f443fc5dbed1e31266ae52d845518954f36e9..a1564b0f204ce1d212ae69aa039b29a38ef6676f 100644
--- a/applications/solvers/combustion/reactingFoam/pcEqn.H
+++ b/applications/solvers/combustion/reactingFoam/pcEqn.H
@@ -1,7 +1,4 @@
rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
volScalarField rAU(1.0/UEqn.A());
volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
@@ -109,19 +106,13 @@ U.correctBoundaryConditions();
fvOptions.correct(U);
K = 0.5*magSqr(U);
-if (thermo.dpdt())
+if (pressureControl.limit(p))
{
- dpdt = fvc::ddt(p);
+ p.correctBoundaryConditions();
+ rho = thermo.rho();
}
-// Recalculate density from the relaxed pressure
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-
-if (!pimple.transonic())
+if (thermo.dpdt())
{
- rho.relax();
+ dpdt = fvc::ddt(p);
}
-
-Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
diff --git a/applications/solvers/combustion/reactingFoam/reactingFoam.C b/applications/solvers/combustion/reactingFoam/reactingFoam.C
index a06eef279c3e9e0627a221eefaeb370502aa21fc..d3e1b19ad38b37b22bcef3b81fd336e0c9570bc7 100644
--- a/applications/solvers/combustion/reactingFoam/reactingFoam.C
+++ b/applications/solvers/combustion/reactingFoam/reactingFoam.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -37,6 +37,7 @@ Description
#include "psiCombustionModel.H"
#include "multivariateScheme.H"
#include "pimpleControl.H"
+#include "pressureControl.H"
#include "fvOptions.H"
#include "localEulerDdtScheme.H"
#include "fvcSmooth.H"
@@ -114,6 +115,8 @@ int main(int argc, char *argv[])
}
}
+ rho = thermo.rho();
+
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H
index 926e69f7c1f7933c402e0289b2e7b0f22cc25bb7..04296ddfc2c53d13dd5fa4f5433e91f6bcb7557f 100644
--- a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H
@@ -1,76 +1,74 @@
-{
- rho = thermo.rho();
+rho = thermo.rho();
- // Thermodynamic density needs to be updated by psi*d(p) after the
- // pressure solution - done in 2 parts. Part 1:
- thermo.rho() -= psi*p;
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
- volScalarField rAU(1.0/UEqn.A());
- surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
- volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
- surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
- surfaceScalarField phiHbyA
+surfaceScalarField phiHbyA
+(
+ "phiHbyA",
(
- "phiHbyA",
- (
- fvc::flux(rho*HbyA)
- + rhorAUf*fvc::ddtCorr(rho, U, phi)
- )
- + phig
- );
-
- MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
- // Update the pressure BCs to ensure flux consistency
- constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
-
- fvScalarMatrix p_rghDDtEqn
+ fvc::flux(rho*HbyA)
+ + rhorAUf*fvc::ddtCorr(rho, U, phi)
+ )
+ + phig
+);
+
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
+
+fvScalarMatrix p_rghDDtEqn
+(
+ fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+ + fvc::div(phiHbyA)
+ ==
+ fvOptions(psi, p_rgh, rho.name())
+);
+
+while (pimple.correctNonOrthogonal())
+{
+ fvScalarMatrix p_rghEqn
(
- fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
- + fvc::div(phiHbyA)
- ==
- fvOptions(psi, p_rgh, rho.name())
+ p_rghDDtEqn
+ - fvm::laplacian(rhorAUf, p_rgh)
);
- while (pimple.correctNonOrthogonal())
+ p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+
+ if (pimple.finalNonOrthogonalIter())
{
- fvScalarMatrix p_rghEqn
- (
- p_rghDDtEqn
- - fvm::laplacian(rhorAUf, p_rgh)
- );
-
- p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
-
- if (pimple.finalNonOrthogonalIter())
- {
- // Calculate the conservative fluxes
- phi = phiHbyA + p_rghEqn.flux();
-
- // Explicitly relax pressure for momentum corrector
- p_rgh.relax();
-
- // Correct the momentum source with the pressure gradient flux
- // calculated from the relaxed pressure
- U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
- U.correctBoundaryConditions();
- fvOptions.correct(U);
- K = 0.5*magSqr(U);
- }
+ // Calculate the conservative fluxes
+ phi = phiHbyA + p_rghEqn.flux();
+
+ // Explicitly relax pressure for momentum corrector
+ p_rgh.relax();
+
+ // Correct the momentum source with the pressure gradient flux
+ // calculated from the relaxed pressure
+ U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
+ U.correctBoundaryConditions();
+ fvOptions.correct(U);
+ K = 0.5*magSqr(U);
}
+}
- p = p_rgh + rho*gh;
-
- // Second part of thermodynamic density update
- thermo.rho() += psi*p;
+p = p_rgh + rho*gh;
- if (thermo.dpdt())
- {
- dpdt = fvc::ddt(p);
- }
+// Thermodynamic density update
+thermo.correctRho(psi*p - psip0);
- #include "rhoEqn.H"
- #include "compressibleContinuityErrs.H"
+if (thermo.dpdt())
+{
+ dpdt = fvc::ddt(p);
}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/createFields.H b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/createFields.H
index 687dffcc13db38d9b916c07df8cee4ee1b11c928..e56035f24a703420b00d2e539ee89f3a7af56090 100644
--- a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/createFields.H
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/createFields.H
@@ -52,6 +52,8 @@ volScalarField& p = thermo.p();
#include "compressibleCreatePhi.H"
+pressureControl pressureControl(p, rho, pimple.dict(), false);
+
mesh.setFluxRequired(p.name());
diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H
deleted file mode 100644
index eff014cf43d847917101add0d0be78cf34b64628..0000000000000000000000000000000000000000
--- a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H
+++ /dev/null
@@ -1,109 +0,0 @@
-{
- rho = thermo.rho();
-
- // Thermodynamic density needs to be updated by psi*d(p) after the
- // pressure solution - done in 2 parts. Part 1:
- thermo.rho() -= psi*p;
-
- volScalarField rAU(1.0/UEqn.A());
- surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
- volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-
- if (pimple.transonic())
- {
- surfaceScalarField phiHbyA
- (
- "phiHbyA",
- (
- fvc::flux(HbyA)
- + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
- )
- );
-
- MRF.makeRelative(phiHbyA);
-
- surfaceScalarField phid("phid", fvc::interpolate(thermo.psi())*phiHbyA);
-
- phiHbyA *= fvc::interpolate(rho);
-
- fvScalarMatrix pDDtEqn
- (
- fvc::ddt(rho) + fvc::div(phiHbyA)
- + correction(psi*fvm::ddt(p) + fvm::div(phid, p))
- );
-
- while (pimple.correctNonOrthogonal())
- {
- fvScalarMatrix pEqn
- (
- pDDtEqn
- - fvm::laplacian(rhorAUf, p)
- ==
- fvOptions(psi, p, rho.name())
- );
-
- pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
-
- if (pimple.finalNonOrthogonalIter())
- {
- phi = phiHbyA + pEqn.flux();
- }
- }
- }
- else
- {
- surfaceScalarField phiHbyA
- (
- "phiHbyA",
- (
- fvc::flux(rho*HbyA)
- + rhorAUf*fvc::ddtCorr(rho, U, phi)
- )
- );
-
- MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
- // Update the pressure BCs to ensure flux consistency
- constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
-
- fvScalarMatrix pDDtEqn
- (
- fvc::ddt(rho) + psi*correction(fvm::ddt(p))
- + fvc::div(phiHbyA)
- ==
- fvOptions(psi, p, rho.name())
- );
-
- while (pimple.correctNonOrthogonal())
- {
- fvScalarMatrix pEqn
- (
- pDDtEqn
- - fvm::laplacian(rhorAUf, p)
- );
-
- pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
-
- if (pimple.finalNonOrthogonalIter())
- {
- phi = phiHbyA + pEqn.flux();
- }
- }
- }
-
- // Second part of thermodynamic density update
- thermo.rho() += psi*p;
-
- #include "rhoEqn.H"
- #include "compressibleContinuityErrs.H"
-
- U = HbyA - rAU*fvc::grad(p);
- U.correctBoundaryConditions();
- fvOptions.correct(U);
- K = 0.5*magSqr(U);
-
- if (thermo.dpdt())
- {
- dpdt = fvc::ddt(p);
- }
-}
diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/rhoReactingFoam.C b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/rhoReactingFoam.C
index 55de274d839e85cbd63f37423b72a967b12bad17..5be897c6c7f4903ea4bf33d561c7b0a1ca48db46 100644
--- a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/rhoReactingFoam.C
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/rhoReactingFoam.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -38,6 +38,7 @@ Description
#include "turbulentFluidThermoModel.H"
#include "multivariateScheme.H"
#include "pimpleControl.H"
+#include "pressureControl.H"
#include "fvOptions.H"
#include "localEulerDdtScheme.H"
#include "fvcSmooth.H"
@@ -100,7 +101,14 @@ int main(int argc, char *argv[])
// --- Pressure corrector loop
while (pimple.correct())
{
- #include "pEqn.H"
+ if (pimple.consistent())
+ {
+ #include "../../../compressible/rhoPimpleFoam/pcEqn.H"
+ }
+ else
+ {
+ #include "../../../compressible/rhoPimpleFoam/pEqn.H"
+ }
}
if (pimple.turbCorr())
diff --git a/applications/solvers/compressible/rhoPimpleFoam/createFields.H b/applications/solvers/compressible/rhoPimpleFoam/createFields.H
index 84f70f2cd02c8c7c0c8a406ea5324cf0dc480a0b..3ec620dc22a0e043433cc9f091b1dfa47a42bfab 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/createFields.H
@@ -42,6 +42,8 @@ volVectorField U
pressureControl pressureControl(p, rho, pimple.dict(), false);
+mesh.setFluxRequired(p.name());
+
Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::turbulenceModel> turbulence
(
@@ -54,8 +56,6 @@ autoPtr<compressible::turbulenceModel> turbulence
)
);
-mesh.setFluxRequired(p.name());
-
Info<< "Creating field dpdt\n" << endl;
volScalarField dpdt
(
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
index 9eecd01d99098f69f2b91358617a0d5a36523cc3..83eeb8fc68b115e1f38a58f73adf6afddf4d5343 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
@@ -1,4 +1,11 @@
-rho = thermo.rho();
+if (!pimple.SIMPLErho())
+{
+ rho = thermo.rho();
+}
+
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
@@ -31,17 +38,17 @@ if (pimple.transonic())
phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
+ fvScalarMatrix pDDtEqn
+ (
+ fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+ + fvc::div(phiHbyA) + fvm::div(phid, p)
+ ==
+ fvOptions(psi, p, rho.name())
+ );
+
while (pimple.correctNonOrthogonal())
{
- fvScalarMatrix pEqn
- (
- fvm::ddt(psi, p)
- + fvc::div(phiHbyA)
- + fvm::div(phid, p)
- - fvm::laplacian(rhorAUf, p)
- ==
- fvOptions(psi, p, rho.name())
- );
+ fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAUf, p));
// Relax the pressure equation to ensure diagonal-dominance
pEqn.relax();
@@ -56,16 +63,17 @@ if (pimple.transonic())
}
else
{
+ fvScalarMatrix pDDtEqn
+ (
+ fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+ + fvc::div(phiHbyA)
+ ==
+ fvOptions(psi, p, rho.name())
+ );
+
while (pimple.correctNonOrthogonal())
{
- fvScalarMatrix pEqn
- (
- fvm::ddt(psi, p)
- + fvc::div(phiHbyA)
- - fvm::laplacian(rhorAUf, p)
- ==
- fvOptions(psi, p, rho.name())
- );
+ fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAUf, p));
pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
@@ -76,6 +84,9 @@ else
}
}
+// Thermodynamic density update
+thermo.correctRho(psi*p - psip0);
+
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
@@ -92,10 +103,9 @@ if (pressureControl.limit(p))
p.correctBoundaryConditions();
rho = thermo.rho();
}
-
-if (!pimple.transonic())
+else if (pimple.SIMPLErho())
{
- rho.relax();
+ rho = thermo.rho();
}
if (thermo.dpdt())
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
index 34657a189287c3fa28cdf4c4cd100f72838399dc..e9b849bc05bbea6e075765e3b68ae4c445fe04b9 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
@@ -1,4 +1,11 @@
-rho = thermo.rho();
+if (!pimple.SIMPLErho())
+{
+ rho = thermo.rho();
+}
+
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
volScalarField rAU(1.0/UEqn.A());
volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
@@ -39,17 +46,17 @@ if (pimple.transonic())
HbyA -= (rAU - rAtU)*fvc::grad(p);
+ fvScalarMatrix pDDtEqn
+ (
+ fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+ + fvc::div(phiHbyA) + fvm::div(phid, p)
+ ==
+ fvOptions(psi, p, rho.name())
+ );
+
while (pimple.correctNonOrthogonal())
{
- fvScalarMatrix pEqn
- (
- fvm::ddt(psi, p)
- + fvc::div(phiHbyA)
- + fvm::div(phid, p)
- - fvm::laplacian(rhorAtU, p)
- ==
- fvOptions(psi, p, rho.name())
- );
+ fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAtU, p));
// Relax the pressure equation to ensure diagonal-dominance
pEqn.relax();
@@ -67,16 +74,17 @@ else
phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf();
HbyA -= (rAU - rAtU)*fvc::grad(p);
+ fvScalarMatrix pDDtEqn
+ (
+ fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+ + fvc::div(phiHbyA)
+ ==
+ fvOptions(psi, p, rho.name())
+ );
+
while (pimple.correctNonOrthogonal())
{
- fvScalarMatrix pEqn
- (
- fvm::ddt(psi, p)
- + fvc::div(phiHbyA)
- - fvm::laplacian(rhorAtU, p)
- ==
- fvOptions(psi, p, rho.name())
- );
+ fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAtU, p));
pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
@@ -87,6 +95,9 @@ else
}
}
+// Thermodynamic density update
+thermo.correctRho(psi*p - psip0);
+
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
@@ -103,10 +114,9 @@ if (pressureControl.limit(p))
p.correctBoundaryConditions();
rho = thermo.rho();
}
-
-if (!pimple.transonic())
+else if (pimple.SIMPLErho())
{
- rho.relax();
+ rho = thermo.rho();
}
if (thermo.dpdt())
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
index 1540f601c128a65f47225a45721d5fb8a46dfd65..ed802ba6d2cd197ff21bd5248d4ef1b2ab9c5b37 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
@@ -36,7 +36,7 @@ if (pimple.transonic())
{
fvScalarMatrix pEqn
(
- fvm::ddt(psi, p)
+ fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+ fvc::div(phiHbyA)
+ fvm::div(phid, p)
- fvm::laplacian(rhorAUf, p)
@@ -62,7 +62,7 @@ else
// Pressure corrector
fvScalarMatrix pEqn
(
- fvm::ddt(psi, p)
+ fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+ fvc::div(phiHbyA)
- fvm::laplacian(rhorAUf, p)
==
@@ -95,11 +95,6 @@ if (pressureControl.limit(p))
rho = thermo.rho();
}
-if (!pimple.transonic())
-{
- rho.relax();
-}
-
{
rhoUf = fvc::interpolate(rho*U);
surfaceVectorField n(mesh.Sf()/mesh.magSf());
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
index 0fbe7717e3569bec9026d98bd30d3ad51dc6a1ba..5d5ed3f5f4c3961b3e2e941e4a42278d73038bb5 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
@@ -122,6 +122,8 @@ int main(int argc, char *argv[])
}
}
+ rho = thermo.rho();
+
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
index 15c87f9c85eecf51d1a6d83e0633646ca1b14758..c0971187539f55678f521ef16e6b7de65ffcb00b 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
@@ -1,115 +1,74 @@
-{
- bool closedVolume = p_rgh.needReference();
-
- dimensionedScalar compressibility = fvc::domainIntegrate(psi);
- bool compressible = (compressibility.value() > SMALL);
-
- rho = thermo.rho();
+rho = thermo.rho();
- // Thermodynamic density needs to be updated by psi*d(p) after the
- // pressure solution - done in 2 parts. Part 1:
- thermo.rho() -= psi*p_rgh;
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
- volScalarField rAU(1.0/UEqn.A());
- surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
- volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
+volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
- surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
- surfaceScalarField phiHbyA
+surfaceScalarField phiHbyA
+(
+ "phiHbyA",
(
- "phiHbyA",
- (
- fvc::flux(rho*HbyA)
- + rhorAUf*fvc::ddtCorr(rho, U, phi)
- )
- + phig
- );
-
- MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
- // Update the pressure BCs to ensure flux consistency
- constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
-
- tmp<fvScalarMatrix> p_rghDDtEqn
+ fvc::flux(rho*HbyA)
+ + rhorAUf*fvc::ddtCorr(rho, U, phi)
+ )
+ + phig
+);
+
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
+
+fvScalarMatrix p_rghDDtEqn
+(
+ fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+ + fvc::div(phiHbyA)
+ ==
+ fvOptions(psi, p_rgh, rho.name())
+);
+
+while (pimple.correctNonOrthogonal())
+{
+ fvScalarMatrix p_rghEqn
(
- new fvScalarMatrix(p_rgh, dimMass/dimTime)
+ p_rghDDtEqn
+ - fvm::laplacian(rhorAUf, p_rgh)
);
- if (compressible)
- {
- p_rghDDtEqn =
- (
- fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
- );
- }
-
- while (pimple.correctNonOrthogonal())
- {
- fvScalarMatrix p_rghEqn
- (
- p_rghDDtEqn()
- + fvc::div(phiHbyA)
- - fvm::laplacian(rhorAUf, p_rgh)
- ==
- fvOptions(psi, p_rgh, rho.name())
- );
-
- p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
-
- if (pimple.finalNonOrthogonalIter())
- {
- // Calculate the conservative fluxes
- phi = phiHbyA + p_rghEqn.flux();
-
- // Explicitly relax pressure for momentum corrector
- p_rgh.relax();
-
- // Correct the momentum source with the pressure gradient flux
- // calculated from the relaxed pressure
- U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
- U.correctBoundaryConditions();
- fvOptions.correct(U);
- K = 0.5*magSqr(U);
- }
- }
-
- p = p_rgh + rho*gh;
+ p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
- // Second part of thermodynamic density update
- thermo.rho() += psi*p_rgh;
-
- if (thermo.dpdt())
+ if (pimple.finalNonOrthogonalIter())
{
- dpdt = fvc::ddt(p);
+ // Calculate the conservative fluxes
+ phi = phiHbyA + p_rghEqn.flux();
+
+ // Explicitly relax pressure for momentum corrector
+ p_rgh.relax();
+
+ // Correct the momentum source with the pressure gradient flux
+ // calculated from the relaxed pressure
+ U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
+ U.correctBoundaryConditions();
+ fvOptions.correct(U);
+ K = 0.5*magSqr(U);
}
+}
- if (compressible)
- {
- #include "rhoEqn.H"
- }
- #include "compressibleContinuityErrs.H"
+p = p_rgh + rho*gh;
- if (closedVolume)
- {
- if (!compressible)
- {
- p += dimensionedScalar
- (
- "p",
- p.dimensions(),
- pRefValue - getRefCellValue(p, pRefCell)
- );
- }
- else
- {
- p += (initialMass - fvc::domainIntegrate(thermo.rho()))
- /compressibility;
- rho = thermo.rho();
- }
- p_rgh = p - rho*gh;
- }
+// Thermodynamic density update
+thermo.correctRho(psi*p - psip0);
- Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value()
- << endl;
+if (thermo.dpdt())
+{
+ dpdt = fvc::ddt(p);
}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
index 09d90906c0237c79dc632492b6e3675db3e23f15..bd955882b11e86a4c53254bd8672be8c9552d97a 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
@@ -1,114 +1,100 @@
-{
- bool closedVolume = p_rgh.needReference();
- dimensionedScalar compressibility = fvc::domainIntegrate(psi);
- bool compressible = (compressibility.value() > SMALL);
+bool closedVolume = p_rgh.needReference();
+dimensionedScalar compressibility = fvc::domainIntegrate(psi);
+bool compressible = (compressibility.value() > SMALL);
- rho = thermo.rho();
+rho = thermo.rho();
- volScalarField rAU("rAU", 1.0/UEqn.A());
- surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
- volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
+volScalarField rAU("rAU", 1.0/UEqn.A());
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
- surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
- surfaceScalarField phiHbyA
+surfaceScalarField phiHbyA
+(
+ "phiHbyA",
(
- "phiHbyA",
- (
- fvc::flux(rho*HbyA)
- + rhorAUf*fvc::ddtCorr(rho, U, phi)
- )
- + phig
- );
+ fvc::flux(rho*HbyA)
+ + rhorAUf*fvc::ddtCorr(rho, U, phi)
+ )
+ + phig
+);
- MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
- // Update the pressure BCs to ensure flux consistency
- constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
- tmp<fvScalarMatrix> p_rghDDtEqn
+{
+ fvScalarMatrix p_rghDDtEqn
(
- new fvScalarMatrix(p_rgh, dimMass/dimTime)
+ fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+ + fvc::div(phiHbyA)
);
- {
- if (compressible)
- {
- p_rghDDtEqn =
- (
- fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
- ==
- fvOptions(psi, p_rgh, rho.name())
- );
- }
-
- // Thermodynamic density needs to be updated by psi*d(p) after the
- // pressure solution - done in 2 parts. Part 1:
- thermo.rho() -= psi*p_rgh;
+ // Thermodynamic density needs to be updated by psi*d(p) after the
+ // pressure solution
+ const volScalarField psip0(psi*p);
- for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
- {
- fvScalarMatrix p_rghEqn
- (
- p_rghDDtEqn()
- + fvc::div(phiHbyA)
- - fvm::laplacian(rhorAUf, p_rgh)
- );
+ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+ {
+ fvScalarMatrix p_rghEqn
+ (
+ p_rghDDtEqn
+ - fvm::laplacian(rhorAUf, p_rgh)
+ );
- p_rghEqn.solve
+ solvPerfp_rgh = p_rghEqn.solve
+ (
+ mesh.solver
(
- mesh.solver
+ p_rgh.select
(
- p_rgh.select
(
- (
- oCorr == nOuterCorr-1
- && corr == nCorr-1
- && nonOrth == nNonOrthCorr
- )
+ oCorr == nOuterCorr-1
+ && corr == nCorr-1
+ && nonOrth == nNonOrthCorr
)
)
- );
-
- if (nonOrth == nNonOrthCorr)
- {
- phi = phiHbyA + p_rghEqn.flux();
- U = HbyA
- + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
- U.correctBoundaryConditions();
- fvOptions.correct(U);
- K = 0.5*magSqr(U);
- }
- }
+ )
+ );
- // Second part of thermodynamic density update
- thermo.rho() += psi*p_rgh;
+ if (nonOrth == nNonOrthCorr)
+ {
+ phi = phiHbyA + p_rghEqn.flux();
+ U = HbyA
+ + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
+ U.correctBoundaryConditions();
+ fvOptions.correct(U);
+ K = 0.5*magSqr(U);
+ }
}
p = p_rgh + rho*gh;
- // Update pressure time derivative if needed
- if (thermo.dpdt())
- {
- dpdt = fvc::ddt(p);
- }
+ // Thermodynamic density update
+ thermo.correctRho(psi*p - psip0);
+}
- if (compressible)
- {
- // Solve continuity
- #include "rhoEqn.H"
- }
- // Update continuity errors
- #include "compressibleContinuityErrors.H"
+// Update pressure time derivative if needed
+if (thermo.dpdt())
+{
+ dpdt = fvc::ddt(p);
+}
- // For closed-volume cases adjust the pressure and density levels
- // to obey overall mass continuity
- if (closedVolume && compressible)
- {
- p += (initialMass - fvc::domainIntegrate(thermo.rho()))
- /compressibility;
- rho = thermo.rho();
- p_rgh = p - rho*gh;
- }
+// Solve continuity
+#include "rhoEqn.H"
+
+// Update continuity errors
+#include "compressibleContinuityErrors.H"
+
+// For closed-volume cases adjust the pressure and density levels
+// to obey overall mass continuity
+if (closedVolume && compressible)
+{
+ p += (initialMass - fvc::domainIntegrate(thermo.rho()))
+ /compressibility;
+ rho = thermo.rho();
+ p_rgh = p - rho*gh;
}
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C b/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
index 544cccdae2ee420c8991c9894cc9ad6b5d1e492c..c7ea95929c950877d6c3f1ce512aa3659988672a 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
+++ b/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -42,6 +42,7 @@ Description
#include "radiationModel.H"
#include "SLGThermo.H"
#include "pimpleControl.H"
+#include "pressureControl.H"
#include "localEulerDdtScheme.H"
#include "fvcSmooth.H"
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
index f6c8ceec21d28dc6064fee4c129d1d875d724b02..f2cb2947186adfcee7a90c18cbdcc63c2e8ff7fb 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
@@ -1,73 +1,71 @@
-{
- rho = thermo.rho();
+rho = thermo.rho();
- // Thermodynamic density needs to be updated by psi*d(p) after the
- // pressure solution - done in 2 parts. Part 1:
- thermo.rho() -= psi*p;
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
- volScalarField rAU(1.0/UEqn.A());
- surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
- volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
- surfaceScalarField phiHbyA
+volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+surfaceScalarField phiHbyA
+(
+ "phiHbyA",
(
- "phiHbyA",
- (
- fvc::flux(rho*HbyA)
- + rhorAUf*fvc::ddtCorr(rho, U, phi)
- )
- );
+ fvc::flux(rho*HbyA)
+ + rhorAUf*fvc::ddtCorr(rho, U, phi)
+ )
+);
- MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
- // Update the pressure BCs to ensure flux consistency
- constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
- fvScalarMatrix pDDtEqn
+fvScalarMatrix pDDtEqn
+(
+ fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+ + fvc::div(phiHbyA)
+ ==
+ parcels.Srho()
+ + fvOptions(psi, p, rho.name())
+);
+
+while (pimple.correctNonOrthogonal())
+{
+ fvScalarMatrix pEqn
(
- fvc::ddt(rho) + psi*correction(fvm::ddt(p))
- + fvc::div(phiHbyA)
- ==
- parcels.Srho()
- + fvOptions(psi, p, rho.name())
+ pDDtEqn
+ - fvm::laplacian(rhorAUf, p)
);
- while (pimple.correctNonOrthogonal())
- {
- fvScalarMatrix pEqn
- (
- pDDtEqn
- - fvm::laplacian(rhorAUf, p)
- );
-
- pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+ pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
- if (pimple.finalNonOrthogonalIter())
- {
- phi = phiHbyA + pEqn.flux();
- }
+ if (pimple.finalNonOrthogonalIter())
+ {
+ phi = phiHbyA + pEqn.flux();
}
+}
- p.relax();
+p.relax();
- // Second part of thermodynamic density update
- thermo.rho() += psi*p;
+// Thermodynamic density update
+thermo.correctRho(psi*p - psip0);
- #include "rhoEqn.H" // NOTE: flux and time scales now inconsistent
- #include "compressibleContinuityErrs.H"
+#include "rhoEqn.H" // NOTE: flux and time scales now inconsistent
+#include "compressibleContinuityErrs.H"
- U = HbyA - rAU*fvc::grad(p);
- U.correctBoundaryConditions();
- fvOptions.correct(U);
- K = 0.5*magSqr(U);
+U = HbyA - rAU*fvc::grad(p);
+U.correctBoundaryConditions();
+fvOptions.correct(U);
+K = 0.5*magSqr(U);
- if (thermo.dpdt())
- {
- dpdt = fvc::ddt(p);
- }
+if (thermo.dpdt())
+{
+ dpdt = fvc::ddt(p);
+}
- rho = thermo.rho();
- rho = max(rho, rhoMin);
- rho = min(rho, rhoMax);
+rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
- Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
-}
+Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H
index 055eff6f0b34ec7e30fb3da5d8ac9f6f96386d25..1ce5db0ec8f1fbd4758d126be362ebb0dd2d3354 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H
@@ -1,57 +1,55 @@
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
+
+volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+tUEqn.clear();
+surfaceScalarField phiHbyA
+(
+ "phiHbyA",
+ fvc::interpolate(rho)*fvc::flux(HbyA)
+);
+
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+
+while (simple.correctNonOrthogonal())
{
- // Thermodynamic density needs to be updated by psi*d(p) after the
- // pressure solution - done in 2 parts. Part 1:
- thermo.rho() -= psi*p;
-
- volScalarField rAU(1.0/UEqn.A());
- surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
- volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
- tUEqn.clear();
- surfaceScalarField phiHbyA
+ fvScalarMatrix pEqn
(
- "phiHbyA",
- fvc::interpolate(rho)*fvc::flux(HbyA)
+ fvc::div(phiHbyA)
+ - fvm::laplacian(rhorAUf, p)
+ ==
+ parcels.Srho()
+ + fvOptions(psi, p, rho.name())
);
- MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+ pEqn.solve();
- // Update the pressure BCs to ensure flux consistency
- constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
-
- while (simple.correctNonOrthogonal())
+ if (simple.finalNonOrthogonalIter())
{
- fvScalarMatrix pEqn
- (
- fvc::div(phiHbyA)
- - fvm::laplacian(rhorAUf, p)
- ==
- parcels.Srho()
- + fvOptions(psi, p, rho.name())
- );
-
- pEqn.solve();
-
- if (simple.finalNonOrthogonalIter())
- {
- phi = phiHbyA + pEqn.flux();
- }
+ phi = phiHbyA + pEqn.flux();
}
+}
- p.relax();
+p.relax();
- // Second part of thermodynamic density update
- thermo.rho() += psi*p;
+// Thermodynamic density update
+thermo.correctRho(psi*p - psip0);
- #include "compressibleContinuityErrs.H"
+#include "compressibleContinuityErrs.H"
- U = HbyA - rAU*fvc::grad(p);
- U.correctBoundaryConditions();
- fvOptions.correct(U);
+U = HbyA - rAU*fvc::grad(p);
+U.correctBoundaryConditions();
+fvOptions.correct(U);
- rho = thermo.rho();
- rho = max(rho, rhoMin);
- rho = min(rho, rhoMax);
- rho.relax();
+rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
- Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
-}
+Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
index 7cd3875e85c6ed1c72b5e155ce687d3634cec71c..e84c19453ffcdfbb87c71099f8b536096a55fba3 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
@@ -104,8 +104,8 @@
}
// Update densities from change in p_rgh
- rho1 += psi1*(p_rgh - p_rgh_0);
- rho2 += psi2*(p_rgh - p_rgh_0);
+ mixture.thermo1().correctRho(psi1*(p_rgh - p_rgh_0));
+ mixture.thermo2().correctRho(psi2*(p_rgh - p_rgh_0));
rho = alpha1*rho1 + alpha2*rho2;
diff --git a/applications/solvers/multiphase/compressibleInterFoam/createFields.H b/applications/solvers/multiphase/compressibleInterFoam/createFields.H
index 46168e6b87f7bfb44bba6c1313f37bf2e72b1787..0616b80ea817e02d219348a6db3de57f73ac41d8 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/createFields.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/createFields.H
@@ -38,8 +38,8 @@ volScalarField& alpha2(mixture.alpha2());
Info<< "Reading thermophysical properties\n" << endl;
-volScalarField& rho1 = mixture.thermo1().rho();
-volScalarField& rho2 = mixture.thermo2().rho();
+const volScalarField& rho1 = mixture.thermo1().rho();
+const volScalarField& rho2 = mixture.thermo2().rho();
volScalarField rho
(
diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
index 51c58ba53689d0b3582c9c0407c81d0136b90f01..93412c9dbaca55742f28d77bd15f74e354a15a93 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
@@ -104,8 +104,8 @@
}
// Update densities from change in p_rgh
- rho1 += psi1*(p_rgh - p_rgh_0);
- rho2 += psi2*(p_rgh - p_rgh_0);
+ mixture.thermo1().correctRho(psi1*(p_rgh - p_rgh_0));
+ mixture.thermo2().correctRho(psi2*(p_rgh - p_rgh_0));
rho = alpha1*rho1 + alpha2*rho2;
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C
index e500ab75872fa4676acf2ca686e929150e1b644f..6d031f05b6dacf5e9a87b44158710b6a79c00e90 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C
+++ b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -45,6 +45,7 @@ void Foam::pimpleControl::read()
solveFlow_ = pimpleDict.lookupOrDefault<Switch>("solveFlow", true);
nCorrPIMPLE_ = pimpleDict.lookupOrDefault<label>("nOuterCorrectors", 1);
nCorrPISO_ = pimpleDict.lookupOrDefault<label>("nCorrectors", 1);
+ SIMPLErho_ = pimpleDict.lookupOrDefault<Switch>("SIMPLErho", false);
turbOnFinalIterOnly_ =
pimpleDict.lookupOrDefault<Switch>("turbOnFinalIterOnly", true);
}
@@ -128,6 +129,7 @@ Foam::pimpleControl::pimpleControl(fvMesh& mesh, const word& dictName)
nCorrPIMPLE_(0),
nCorrPISO_(0),
corrPISO_(0),
+ SIMPLErho_(false),
turbOnFinalIterOnly_(true),
converged_(false)
{
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H
index 033ce84d2f54a6029712f74f7c49c584fc96402d..5e87b20ea33832a19914a830f794341584ae1a9a 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H
+++ b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -81,6 +81,10 @@ protected:
//- Current PISO corrector
label corrPISO_;
+ //- Flag to indicate whether to update density in SIMPLE
+ // rather than PISO mode
+ bool SIMPLErho_;
+
//- Flag to indicate whether to only solve turbulence on final iter
bool turbOnFinalIterOnly_;
@@ -128,6 +132,10 @@ public:
//- Current PISO corrector index
inline label corrPISO() const;
+ //- Flag to indicate whether to update density in SIMPLE
+ // rather than PISO mode
+ inline bool SIMPLErho() const;
+
// Solution control
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H
index f6c6e0ce35fbd7bd100455ed0035f46ba28cb34b..deed2958c0e4ba023fd4f626e1dd9717ee198a83 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H
+++ b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -43,6 +43,12 @@ inline Foam::label Foam::pimpleControl::corrPISO() const
}
+inline bool Foam::pimpleControl::SIMPLErho() const
+{
+ return SIMPLErho_;
+}
+
+
inline bool Foam::pimpleControl::correct()
{
corrPISO_++;
diff --git a/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H b/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H
index e18cda73431f6485264520edd12ffb4e872f8acc..995b1e9ad5713b0b16b67def9d72406560589bb4 100644
--- a/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H
+++ b/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2012-2015 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2012-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -103,6 +103,10 @@ public:
// Access to thermodynamic state variables
+ //- Add the given density correction to the density field.
+ // Used to update the density field following pressure solution
+ virtual void correctRho(const volScalarField& deltaRho) = 0;
+
//- Compressibility [s^2/m^2]
virtual const volScalarField& psi() const = 0;
diff --git a/src/thermophysicalModels/basic/psiThermo/psiThermo.C b/src/thermophysicalModels/basic/psiThermo/psiThermo.C
index 630402280a448259cd8ca3d233d95d2cdd28222d..fff1f8dfd66b243cc4405c406dfb21852a7bad84 100644
--- a/src/thermophysicalModels/basic/psiThermo/psiThermo.C
+++ b/src/thermophysicalModels/basic/psiThermo/psiThermo.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -102,6 +102,10 @@ Foam::tmp<Foam::scalarField> Foam::psiThermo::rho(const label patchi) const
}
+void Foam::psiThermo::correctRho(const Foam::volScalarField& deltaRho)
+{}
+
+
const Foam::volScalarField& Foam::psiThermo::psi() const
{
return psi_;
diff --git a/src/thermophysicalModels/basic/psiThermo/psiThermo.H b/src/thermophysicalModels/basic/psiThermo/psiThermo.H
index 76404a026b36f3668d94da12296b0d0ab93b8c7e..454eb5674f20a1ecdf357e1c1b4803929bcd1254 100644
--- a/src/thermophysicalModels/basic/psiThermo/psiThermo.H
+++ b/src/thermophysicalModels/basic/psiThermo/psiThermo.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -115,6 +115,11 @@ public:
// Fields derived from thermodynamic state variables
+ //- Add the given density correction to the density field.
+ // Used to update the density field following pressure solution.
+ // For psiThermo does nothing.
+ virtual void correctRho(const volScalarField& deltaRho);
+
//- Density [kg/m^3] - uses current value of pressure
virtual tmp<volScalarField> rho() const;
diff --git a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
index ada7716c4a8ee0b69c8bd722b2c36e3eec1b5e7d..2f2c5d893f052de20ec4f76144f742f84e96a590 100644
--- a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
+++ b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -173,6 +173,12 @@ Foam::volScalarField& Foam::rhoThermo::rho()
}
+void Foam::rhoThermo::correctRho(const Foam::volScalarField& deltaRho)
+{
+ rho_ += deltaRho;
+}
+
+
const Foam::volScalarField& Foam::rhoThermo::psi() const
{
return psi_;
diff --git a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.H b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.H
index 680094b90f0a6e28b6903961e8507c20489272a0..a4abe62a5360503c869ca4f423458a262e1578c0 100644
--- a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.H
+++ b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -136,6 +136,10 @@ public:
//- Return non-const access to the local density field [kg/m^3]
virtual volScalarField& rho();
+ //- Add the given density correction to the density field.
+ // Used to update the density field following pressure solution
+ virtual void correctRho(const volScalarField& deltaRho);
+
//- Compressibility [s^2/m^2]
virtual const volScalarField& psi() const;
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSolution b/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSolution
index 5dc858a4d40859793b8ee4a2c409e71c0379f93b..828b16db4fc689ea030fa89f9f96831967670f43 100644
--- a/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSolution
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSolution
@@ -54,6 +54,7 @@ PIMPLE
nCorrectors 1;
nNonOrthogonalCorrectors 0;
consistent yes;
+ SIMPLErho yes;
pMaxFactor 1.5;
pMinFactor 0.9;
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/T b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/T
new file mode 100644
index 0000000000000000000000000000000000000000..2df93ba37eda80c3a2c2c94e2a1a55659c5aa4f6
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/T
@@ -0,0 +1,42 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: dev |
+| \\ / A nd | Web: www.OpenFOAM.org |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class volScalarField;
+ object T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [0 0 0 1 0 0 0];
+
+internalField uniform 300;
+
+boundaryField
+{
+ Default_Boundary_Region
+ {
+ type fixedValue;
+ value uniform 350;
+ }
+
+ inlet
+ {
+ type fixedValue;
+ value $internalField;
+ }
+
+ outlet
+ {
+ type inletOutlet;
+ inletValue $internalField;
+ }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/U b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/U
new file mode 100644
index 0000000000000000000000000000000000000000..03d1d6a8d856e001773f425f707425c765f90b21
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/U
@@ -0,0 +1,43 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: dev |
+| \\ / A nd | Web: www.OpenFOAM.org |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class volVectorField;
+ object U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [0 1 -1 0 0 0 0];
+
+internalField uniform (0 0 0);
+
+boundaryField
+{
+ Default_Boundary_Region
+ {
+ type noSlip;
+ }
+
+ inlet
+ {
+ type flowRateInletVelocity;
+ massFlowRate constant 5;
+ rhoInlet 1000; // Guess for rho
+ }
+
+ outlet
+ {
+ type pressureInletOutletVelocity;
+ value $internalField;
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/alphat b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/alphat
new file mode 100644
index 0000000000000000000000000000000000000000..61043c15b176c51a24630b5f6f8f22b9921eda30
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/alphat
@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: dev |
+| \\ / A nd | Web: www.OpenFOAM.org |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class volScalarField;
+ location "0";
+ object alphat;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [1 -1 -1 0 0 0 0];
+
+internalField uniform 0;
+
+boundaryField
+{
+ Default_Boundary_Region
+ {
+ type compressible::alphatWallFunction;
+ Prt 0.85;
+ value uniform 0;
+ }
+
+ inlet
+ {
+ type calculated;
+ value uniform 0;
+ }
+
+ outlet
+ {
+ type calculated;
+ value uniform 0;
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/epsilon b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/epsilon
new file mode 100644
index 0000000000000000000000000000000000000000..20d9e01d3fd428b66e421cc471ac3f25a6dfe970
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/epsilon
@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: dev |
+| \\ / A nd | Web: www.OpenFOAM.org |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class volScalarField;
+ location "0";
+ object epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [0 2 -3 0 0 0 0];
+
+internalField uniform 200;
+
+boundaryField
+{
+ Default_Boundary_Region
+ {
+ type epsilonWallFunction;
+ Cmu 0.09;
+ kappa 0.41;
+ E 9.8;
+ value uniform 200;
+ }
+
+ inlet
+ {
+ type turbulentMixingLengthDissipationRateInlet;
+ mixingLength 0.005;
+ value uniform 200;
+ }
+
+ outlet
+ {
+ type inletOutlet;
+ inletValue uniform 200;
+ value uniform 200;
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/k b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/k
new file mode 100644
index 0000000000000000000000000000000000000000..4bddfd6596532883ae70970018ab4c36036bc17d
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/k
@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: dev |
+| \\ / A nd | Web: www.OpenFOAM.org |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class volScalarField;
+ location "0";
+ object k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [0 2 -2 0 0 0 0];
+
+internalField uniform 1;
+
+boundaryField
+{
+ Default_Boundary_Region
+ {
+ type kqRWallFunction;
+ value uniform 1;
+ }
+
+ inlet
+ {
+ type turbulentIntensityKineticEnergyInlet;
+ intensity 0.05;
+ value uniform 1;
+ }
+
+ outlet
+ {
+ type inletOutlet;
+ inletValue uniform 1;
+ value uniform 1;
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/nut b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/nut
new file mode 100644
index 0000000000000000000000000000000000000000..a7ffddb820aee299c106ee28be41ab84c21c18c9
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/nut
@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: dev |
+| \\ / A nd | Web: www.OpenFOAM.org |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class volScalarField;
+ location "0";
+ object nut;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [0 2 -1 0 0 0 0];
+
+internalField uniform 0;
+
+boundaryField
+{
+ Default_Boundary_Region
+ {
+ type nutkWallFunction;
+ Cmu 0.09;
+ kappa 0.41;
+ E 9.8;
+ value uniform 0;
+ }
+
+ inlet
+ {
+ type calculated;
+ value uniform 0;
+ }
+
+ outlet
+ {
+ type calculated;
+ value uniform 0;
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/p b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/p
new file mode 100644
index 0000000000000000000000000000000000000000..b4c015f22217d0f122bd94673345805f084a6201
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/p
@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: dev |
+| \\ / A nd | Web: www.OpenFOAM.org |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class volScalarField;
+ object p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [1 -1 -2 0 0 0 0];
+
+internalField uniform 1e5;
+
+boundaryField
+{
+ Default_Boundary_Region
+ {
+ type zeroGradient;
+ }
+
+ inlet
+ {
+ type zeroGradient;
+ }
+
+ outlet
+ {
+ type fixedValue;
+ value $internalField;
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/constant/thermophysicalProperties b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/constant/thermophysicalProperties
new file mode 100644
index 0000000000000000000000000000000000000000..426ab721160866a4dbd6eda166dfb8a5e083df74
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/constant/thermophysicalProperties
@@ -0,0 +1,31 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: dev |
+| \\ / A nd | Web: www.OpenFOAM.org |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class dictionary;
+ location "constant";
+ object thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+ type heRhoThermo;
+ mixture pureMixture;
+ properties liquid;
+ energy sensibleInternalEnergy;
+}
+
+mixture
+{
+ H2O;
+}
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/constant/turbulenceProperties b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/constant/turbulenceProperties
new file mode 100644
index 0000000000000000000000000000000000000000..cd2daf8229ba0b2be3dca97cab3a5c08f20b0e8a
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/constant/turbulenceProperties
@@ -0,0 +1,30 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: dev |
+| \\ / A nd | Web: www.OpenFOAM.org |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class dictionary;
+ location "constant";
+ object turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType RAS;
+
+RAS
+{
+ RASModel kEpsilon;
+
+ turbulence on;
+
+ printCoeffs on;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/blockMeshDict b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/blockMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..9bec6a9547f35be865bbf60a0a4d8c18b0258980
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/blockMeshDict
@@ -0,0 +1,127 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: dev |
+| \\ / A nd | Web: www.OpenFOAM.org |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class dictionary;
+ object blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 0.001;
+
+vertices
+(
+ // front-plane: z = +25mm
+ // inlet region
+ ( -50 25 25) // pt 0
+ ( 0 25 25) // pt 1
+ ( -50 75 25) // pt 2
+ ( 0 75 25) // pt 3
+ // outlet region
+ ( -500 -75 25) // pt 4
+ ( 0 -75 25) // pt 5
+ ( -500 -25 25) // pt 6
+ ( 0 -25 25) // pt 7
+ // bend mid-points
+ ( 25 0 25) // pt 8
+ ( 75 0 25) // pt 9
+ // back-plane: z = -25mm
+ // inlet region
+ ( -50 25 -25) // pt 0 + 10
+ ( 0 25 -25) // pt 1 + 10
+ ( -50 75 -25) // pt 2 + 10
+ ( 0 75 -25) // pt 3 + 10
+ // outlet region
+ ( -500 -75 -25) // pt 4 + 10
+ ( 0 -75 -25) // pt 5 + 10
+ ( -500 -25 -25) // pt 7 + 10
+ ( 0 -25 -25) // pt 8 + 10
+ // bend mid-points
+ ( 25 0 -25) // pt 8 + 10
+ ( 75 0 -25) // pt 9 + 10
+);
+
+blocks
+(
+ hex (0 1 11 10 2 3 13 12) inlet ( 20 20 20) simpleGrading (1 1 1)
+ hex (4 5 15 14 6 7 17 16) outlet (200 20 20) simpleGrading (1 1 1)
+
+ hex (1 8 18 11 3 9 19 13) bend1 ( 30 20 20) simpleGrading (1 1 1)
+ hex (5 9 19 15 7 8 18 17) bend2 ( 30 20 20) simpleGrading (1 1 1)
+);
+
+edges
+(
+ // block 2
+ arc 1 8 ( 17.678 17.678 25)
+ arc 11 18 ( 17.678 17.678 -25)
+ arc 3 9 ( 53.033 53.033 25)
+ arc 13 19 ( 53.033 53.033 -25)
+ // block 3
+ arc 7 8 ( 17.678 -17.678 25)
+ arc 17 18 ( 17.678 -17.678 -25)
+ arc 5 9 ( 53.033 -53.033 25)
+ arc 15 19 ( 53.033 -53.033 -25)
+);
+
+boundary
+(
+ // is there no way of defining all my 'defaultFaces' to be 'wall'?
+ Default_Boundary_Region
+ {
+ type wall;
+ faces
+ (
+ // block0
+ ( 0 1 3 2 )
+ ( 11 10 12 13 )
+ ( 0 10 11 1 )
+ ( 2 3 13 12 )
+ // block1
+ ( 4 5 7 6 )
+ ( 15 14 16 17 )
+ ( 4 14 15 5 )
+ ( 6 7 17 16 )
+ // block2
+ ( 1 8 9 3 )
+ ( 18 11 13 19 )
+ ( 3 9 19 13 )
+ ( 1 11 18 8 )
+ // block3
+ ( 5 9 8 7 )
+ ( 19 15 17 18 )
+ ( 5 15 19 9 )
+ ( 7 8 18 17 )
+ );
+ }
+ inlet
+ {
+ type patch;
+ faces
+ (
+ (0 2 12 10)
+ );
+ }
+
+ outlet
+ {
+ type patch;
+ faces
+ (
+ (4 6 16 14)
+ );
+ }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/controlDict b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/controlDict
new file mode 100644
index 0000000000000000000000000000000000000000..07384c85ac8ed84b73dc54c7c8848c466b5037f9
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/controlDict
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: dev |
+| \\ / A nd | Web: www.OpenFOAM.org |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class dictionary;
+ location "system";
+ object controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application rhoPimpleFoam;
+
+startFrom startTime;
+
+startTime 0;
+
+stopAt endTime;
+
+endTime 0.5;
+
+deltaT 2e-3;
+
+writeControl timeStep;
+
+writeInterval 10;
+
+purgeWrite 0;
+
+writeFormat ascii;
+
+writePrecision 6;
+
+writeCompression off;
+
+timeFormat general;
+
+timePrecision 6;
+
+graphFormat raw;
+
+runTimeModifiable true;
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/decomposeParDict b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/decomposeParDict
new file mode 100644
index 0000000000000000000000000000000000000000..9e844a62ba7901b16f2d0c8c470ae7b5acb2504a
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/decomposeParDict
@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: dev |
+| \\ / A nd | Web: www.OpenFOAM.org |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class dictionary;
+ location "system";
+ object decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 8;
+
+method hierarchical;
+
+simpleCoeffs
+{
+ n (8 1 1);
+ delta 0.001;
+}
+
+hierarchicalCoeffs
+{
+ n (4 2 1);
+ delta 0.001;
+ order xyz;
+}
+
+manualCoeffs
+{
+ dataFile "";
+}
+
+distributed no;
+
+roots ( );
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/fvSchemes b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/fvSchemes
new file mode 100644
index 0000000000000000000000000000000000000000..8c9c2b631b2a4c7272f8a23eae7edb61bc216d4f
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/fvSchemes
@@ -0,0 +1,60 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: dev |
+| \\ / A nd | Web: www.OpenFOAM.org |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class dictionary;
+ location "system";
+ object fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+ default Euler;
+}
+
+gradSchemes
+{
+ default Gauss linear;
+
+ limited cellLimited Gauss linear 1;
+}
+
+divSchemes
+{
+ default none;
+
+ div(phi,U) Gauss linearUpwind limited;
+ div(phi,e) Gauss linearUpwind limited;
+ div(phi,epsilon) Gauss linearUpwind limited;
+ div(phi,k) Gauss linearUpwind limited;
+ div(phi,K) Gauss linearUpwind limited;
+ div(phiv,p) Gauss linearUpwind limited;
+
+ div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
+}
+
+laplacianSchemes
+{
+ default Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+ default linear;
+}
+
+snGradSchemes
+{
+ default corrected;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/fvSolution b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/fvSolution
new file mode 100644
index 0000000000000000000000000000000000000000..173adf6c59d86acbb1d11971622a1487082f3aec
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/fvSolution
@@ -0,0 +1,63 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: dev |
+| \\ / A nd | Web: www.OpenFOAM.org |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class dictionary;
+ location "system";
+ object fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+ "rho.*"
+ {
+ solver nthn;
+ }
+
+ "p.*"
+ {
+ solver GAMG;
+ smoother GaussSeidel;
+
+ tolerance 1e-7;
+ relTol 0;
+ }
+
+ "(U|e|k|epsilon).*"
+ {
+ solver PBiCGStab;
+ preconditioner DILU;
+
+ tolerance 1e-7;
+ relTol 0;
+ }
+}
+
+PIMPLE
+{
+ nCorrectors 3;
+ nNonOrthogonalCorrectors 0;
+ pMinFactor 0.1;
+ pMaxFactor 1.5;
+
+ transonic no;
+ consistent no;
+}
+
+relaxationFactors
+{
+ equations
+ {
+ ".*" 1;
+ }
+}
+
+// ************************************************************************* //