From d1b651533f9295aa4ac1656d87a5a5e45b79f37c Mon Sep 17 00:00:00 2001
From: sergio <sergio>
Date: Thu, 1 Jun 2017 12:39:28 -0700
Subject: [PATCH] Modification on rhoPimpleFoam pEq's for handling rho thermo
and incompressible EoS. Adding rho limiters if p is limited. This is
important when LTS stepping or large Co number is used.
Updating rhoBuoyantPimpleFoam to handle closed domain for rho thermo and incompressible Eos.
Consolidating chtMultiRegionSimpleFoam and chtMultiRegionFoam pEqs to use the same formulation as rhoBuoyantPimpleFoam and
rhoBuoyantSimpleFoam
---
.../compressible/rhoPimpleFoam/createFields.H | 23 ++++++++
.../solvers/compressible/rhoPimpleFoam/pEqn.H | 11 +++-
.../compressible/rhoPimpleFoam/pcEqn.H | 12 +++-
.../heatTransfer/buoyantPimpleFoam/pEqn.H | 46 +++++++++++++--
.../buoyantSimpleFoam/createFields.H | 23 ++++++++
.../heatTransfer/buoyantSimpleFoam/pEqn.H | 17 +++---
.../chtMultiRegionSimpleFoam/fluid/pEqn.H | 20 ++++++-
.../fluid/createFluidFields.H | 19 +++++++
.../chtMultiRegionFoam/fluid/pEqn.H | 57 ++++++++++++++-----
.../fluid/setRegionFluidFields.H | 3 +
.../basic/fluidThermo/fluidThermo.H | 7 ++-
.../basic/psiThermo/psiThermo.C | 8 ++-
.../basic/psiThermo/psiThermo.H | 7 ++-
.../basic/rhoThermo/rhoThermo.C | 10 +++-
.../basic/rhoThermo/rhoThermo.H | 8 ++-
15 files changed, 229 insertions(+), 42 deletions(-)
diff --git a/applications/solvers/compressible/rhoPimpleFoam/createFields.H b/applications/solvers/compressible/rhoPimpleFoam/createFields.H
index 3ec620dc22a..8079428225b 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/createFields.H
@@ -73,3 +73,26 @@ Info<< "Creating field kinetic energy K\n" << endl;
volScalarField K("K", 0.5*magSqr(U));
#include "createMRF.H"
+
+
+dimensionedScalar rhoMax
+(
+ dimensionedScalar::lookupOrDefault
+ (
+ "rhoMax",
+ pimple.dict(),
+ dimDensity,
+ GREAT
+ )
+);
+
+dimensionedScalar rhoMin
+(
+ dimensionedScalar::lookupOrDefault
+ (
+ "rhoMin",
+ pimple.dict(),
+ dimDensity,
+ 0
+ )
+);
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
index 83eeb8fc68b..3a0ca9d042a 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
@@ -84,9 +84,6 @@ else
}
}
-// Thermodynamic density update
-thermo.correctRho(psi*p - psip0);
-
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
@@ -101,12 +98,20 @@ K = 0.5*magSqr(U);
if (pressureControl.limit(p))
{
p.correctBoundaryConditions();
+ thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
rho = thermo.rho();
}
else if (pimple.SIMPLErho())
{
+ thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
rho = thermo.rho();
}
+else
+{
+ thermo.correctRho(psi*p - psip0, rhoMin, rhoMax) ;
+}
+
+Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
if (thermo.dpdt())
{
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
index e9b849bc05b..546dfaf251b 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
@@ -95,9 +95,6 @@ else
}
}
-// Thermodynamic density update
-thermo.correctRho(psi*p - psip0);
-
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
@@ -112,12 +109,21 @@ K = 0.5*magSqr(U);
if (pressureControl.limit(p))
{
p.correctBoundaryConditions();
+ thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
rho = thermo.rho();
}
+
else if (pimple.SIMPLErho())
{
+ thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
rho = thermo.rho();
}
+else
+{
+ thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
+}
+
+Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
if (thermo.dpdt())
{
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
index c0971187539..0f53379eef7 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
@@ -1,3 +1,7 @@
+
+dimensionedScalar compressibility = fvc::domainIntegrate(psi);
+bool compressible = (compressibility.value() > SMALL);
+
rho = thermo.rho();
// Thermodynamic density needs to be updated by psi*d(p) after the
@@ -41,6 +45,12 @@ while (pimple.correctNonOrthogonal())
- fvm::laplacian(rhorAUf, p_rgh)
);
+ p_rghEqn.setReference
+ (
+ pRefCell,
+ compressible ? getRefCellValue(p_rgh, pRefCell) : pRefValue
+ );
+
p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
if (pimple.finalNonOrthogonalIter())
@@ -63,12 +73,40 @@ while (pimple.correctNonOrthogonal())
p = p_rgh + rho*gh;
// Thermodynamic density update
-thermo.correctRho(psi*p - psip0);
+//thermo.correctRho(psi*p - psip0);
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+if (p_rgh.needReference())
+{
+ if (!compressible)
+ {
+ p += dimensionedScalar
+ (
+ "p",
+ p.dimensions(),
+ pRefValue - getRefCellValue(p, pRefCell)
+ );
+ }
+ else
+ {
+ p += (initialMass - fvc::domainIntegrate(psi*p))
+ /compressibility;
+ thermo.correctRho(psi*p - psip0);
+ rho = thermo.rho();
+
+ }
+ p_rgh = p - rho*gh;
+}
+else
+{
+ thermo.correctRho(psi*p - psip0);
+}
+
+rho = thermo.rho();
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);
}
-
-#include "rhoEqn.H"
-#include "compressibleContinuityErrs.H"
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
index b9f275407f4..64bc8558c8f 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
@@ -88,3 +88,26 @@ dimensionedScalar totalVolume = sum(mesh.V());
#include "createMRF.H"
#include "createRadiationModel.H"
+
+dimensionedScalar rhoMax
+(
+ dimensionedScalar::lookupOrDefault
+ (
+ "rhoMax",
+ simple.dict(),
+ dimDensity,
+ GREAT
+ )
+);
+
+dimensionedScalar rhoMin
+(
+ dimensionedScalar::lookupOrDefault
+ (
+ "rhoMin",
+ simple.dict(),
+ dimDensity,
+ 0
+ )
+);
+
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
index 647cf954bb2..aa3ad0f163e 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
@@ -1,7 +1,4 @@
{
- rho = thermo.rho();
- rho.relax();
-
volScalarField rAU("rAU", 1.0/UEqn.A());
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
@@ -24,6 +21,8 @@
// Update the pressure BCs to ensure flux consistency
constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
+ dimensionedScalar compressibility = fvc::domainIntegrate(psi);
+ bool compressible = (compressibility.value() > SMALL);
while (simple.correctNonOrthogonal())
{
@@ -32,7 +31,12 @@
fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phiHbyA)
);
- p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+ p_rghEqn.setReference
+ (
+ pRefCell,
+ compressible ? getRefCellValue(p_rgh, pRefCell) : pRefValue
+ );
+
p_rghEqn.solve();
if (simple.finalNonOrthogonalIter())
@@ -51,12 +55,9 @@
}
}
- #include "continuityErrs.H"
-
p = p_rgh + rho*gh;
- dimensionedScalar compressibility = fvc::domainIntegrate(psi);
- bool compressible = (compressibility.value() > SMALL);
+ #include "continuityErrs.H"
// For closed-volume cases adjust the pressure level
// to obey overall mass continuity
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
index 7fe544d6a57..5b8286ca4b8 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
@@ -1,8 +1,10 @@
{
+ /*
rho = thermo.rho();
rho = max(rho, rhoMin[i]);
rho = min(rho, rhoMax[i]);
rho.relax();
+ */
volScalarField rAU("rAU", 1.0/UEqn.A());
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
@@ -67,10 +69,22 @@
// For closed-volume cases adjust the pressure level
// to obey overall mass continuity
- if (closedVolume && compressible)
+ if (closedVolume)
{
- p += (initialMass - fvc::domainIntegrate(thermo.rho()))
- /compressibility;
+ if (!compressible)
+ {
+ p += dimensionedScalar
+ (
+ "p",
+ p.dimensions(),
+ pRefValue - getRefCellValue(p, pRefCell)
+ );
+ }
+ else
+ {
+ p += (initialMass - fvc::domainIntegrate(psi*p))
+ /compressibility;
+ }
p_rgh = p - rho*gh;
}
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
index be3ec79a2cc..87e372fcae3 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
@@ -19,6 +19,9 @@ List<bool> frozenFlowFluid(fluidRegions.size(), false);
PtrList<IOMRFZoneList> MRFfluid(fluidRegions.size());
PtrList<fv::options> fluidFvOptions(fluidRegions.size());
+List<label> refCellFluid(fluidRegions.size());
+List<scalar> refValueFluid(fluidRegions.size());
+
// Populate fluid field pointer lists
forAll(fluidRegions, i)
{
@@ -248,4 +251,20 @@ forAll(fluidRegions, i)
);
turbulence[i].validate();
+
+ refCellFluid[i] = 0;
+ refValueFluid[i] = 0.0;
+
+ if (p_rghFluid[i].needReference())
+ {
+ setRefCell
+ (
+ thermoFluid[i].p(),
+ p_rghFluid[i],
+ pimpleDict,
+ refCellFluid[i],
+ refValueFluid[i]
+ );
+ }
+
}
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
index d515c6a1838..0c9195ee9de 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
@@ -4,6 +4,10 @@ bool compressible = (compressibility.value() > SMALL);
rho = thermo.rho();
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
+
volScalarField rAU("rAU", 1.0/UEqn.A());
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
@@ -32,10 +36,6 @@ constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
+ fvc::div(phiHbyA)
);
- // 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
@@ -44,6 +44,12 @@ constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
- fvm::laplacian(rhorAUf, p_rgh)
);
+ p_rghEqn.setReference
+ (
+ pRefCell,
+ compressible ? getRefCellValue(p_rgh, pRefCell) : pRefValue
+ );
+
p_rghEqn.solve
(
mesh.solver
@@ -62,6 +68,9 @@ constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
if (nonOrth == nNonOrthCorr)
{
phi = phiHbyA + p_rghEqn.flux();
+
+ p_rgh.relax();
+
U = HbyA
+ rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
U.correctBoundaryConditions();
@@ -73,15 +82,10 @@ constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
p = p_rgh + rho*gh;
// Thermodynamic density update
- thermo.correctRho(psi*p - psip0);
+ //thermo.correctRho(psi*p - psip0);
}
-// Update pressure time derivative if needed
-if (thermo.dpdt())
-{
- dpdt = fvc::ddt(p);
-}
// Solve continuity
#include "rhoEqn.H"
@@ -91,10 +95,35 @@ if (thermo.dpdt())
// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
-if (closedVolume && compressible)
+if (closedVolume)
{
- p += (initialMass - fvc::domainIntegrate(thermo.rho()))
- /compressibility;
- rho = thermo.rho();
+ if (!compressible)
+ {
+ p += dimensionedScalar
+ (
+ "p",
+ p.dimensions(),
+ pRefValue - getRefCellValue(p, pRefCell)
+ );
+ }
+ else
+ {
+ p += (initialMass - fvc::domainIntegrate(psi*p))
+ /compressibility;
+ thermo.correctRho(psi*p - psip0);
+ rho = thermo.rho();
+ }
p_rgh = p - rho*gh;
}
+else
+{
+ thermo.correctRho(psi*p - psip0);
+}
+
+rho = thermo.rho();
+
+// Update pressure time derivative if needed
+if (thermo.dpdt())
+{
+ dpdt = fvc::ddt(p);
+}
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
index 8d3a291bef9..5edd27ac5e5 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
@@ -33,3 +33,6 @@
);
const bool frozenFlow = frozenFlowFluid[i];
+
+ const label pRefCell = refCellFluid[i];
+ const scalar pRefValue = refValueFluid[i];
diff --git a/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H b/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H
index 995b1e9ad57..bb0364dc799 100644
--- a/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H
+++ b/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H
@@ -105,7 +105,12 @@ public:
//- 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;
+ virtual void correctRho
+ (
+ const volScalarField& deltaRho,
+ const dimensionedScalar& rhoMin,
+ const dimensionedScalar& rhoMax
+ ) = 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 fff1f8dfd66..b3e65677d9f 100644
--- a/src/thermophysicalModels/basic/psiThermo/psiThermo.C
+++ b/src/thermophysicalModels/basic/psiThermo/psiThermo.C
@@ -102,10 +102,14 @@ Foam::tmp<Foam::scalarField> Foam::psiThermo::rho(const label patchi) const
}
-void Foam::psiThermo::correctRho(const Foam::volScalarField& deltaRho)
+void Foam::psiThermo::correctRho
+(
+ const Foam::volScalarField& deltaRho,
+ const dimensionedScalar& rhoMin,
+ const dimensionedScalar& rhoMax
+)
{}
-
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 454eb5674f2..b2e05a9e9e0 100644
--- a/src/thermophysicalModels/basic/psiThermo/psiThermo.H
+++ b/src/thermophysicalModels/basic/psiThermo/psiThermo.H
@@ -118,7 +118,12 @@ public:
//- 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);
+ virtual void correctRho
+ (
+ const volScalarField& deltaRho,
+ const dimensionedScalar& rhoMin,
+ const dimensionedScalar& rhoMax
+ );
//- 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 2f2c5d893f0..bbc1a64ba57 100644
--- a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
+++ b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
@@ -173,12 +173,18 @@ Foam::volScalarField& Foam::rhoThermo::rho()
}
-void Foam::rhoThermo::correctRho(const Foam::volScalarField& deltaRho)
+void Foam::rhoThermo::correctRho
+(
+ const Foam::volScalarField& deltaRho,
+ const dimensionedScalar& rhoMin,
+ const dimensionedScalar& rhoMax
+)
{
rho_ += deltaRho;
+ rho_ = max(rho_, rhoMin);
+ rho_ = min(rho_, rhoMax);
}
-
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 a4abe62a536..fd8bb6b81c5 100644
--- a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.H
+++ b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.H
@@ -138,7 +138,13 @@ public:
//- Add the given density correction to the density field.
// Used to update the density field following pressure solution
- virtual void correctRho(const volScalarField& deltaRho);
+ // Limit thermo rho between rhoMin and rhoMax
+ virtual void correctRho
+ (
+ const volScalarField& deltaRho,
+ const dimensionedScalar& rhoMin,
+ const dimensionedScalar& rhoMax
+ );
//- Compressibility [s^2/m^2]
virtual const volScalarField& psi() const;
--
GitLab