Commit a1c8cde3 authored by Henry Weller's avatar Henry Weller
Browse files

rhoSimpleFoam: added support for compressible liquid flows

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

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

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

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

mixture
{
    H2O;
}

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

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

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

SIMPLE
{
    nNonOrthogonalCorrectors 0;

    pMinFactor      0.1;
    pMaxFactor      1.5;

    transonic       yes;
    consistent      yes;

    residualControl
    {
        p               1e-3;
        U               1e-4;
        e               1e-3;
        "(k|epsilon|omega)" 1e-3;
    }
}
parent c393b268
Info<< "Reading thermophysical properties\n" << endl;
autoPtr<psiThermo> pThermo
autoPtr<fluidThermo> pThermo
(
psiThermo::New(mesh)
fluidThermo::New(mesh)
);
psiThermo& thermo = pThermo();
fluidThermo& thermo = pThermo();
thermo.validate(args.executable(), "h", "e");
volScalarField rho
......@@ -38,35 +38,10 @@ volVectorField U
#include "compressibleCreatePhi.H"
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, simple.dict(), pRefCell, pRefValue);
pressureControl pressureControl(p, rho, simple.dict());
mesh.setFluxRequired(p.name());
dimensionedScalar rhoMax
(
dimensionedScalar::lookupOrDefault
(
"rhoMax",
simple.dict(),
dimDensity,
GREAT
)
);
dimensionedScalar rhoMin
(
dimensionedScalar::lookupOrDefault
(
"rhoMin",
simple.dict(),
dimDensity,
0
)
);
Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::turbulenceModel> turbulence
(
......
......@@ -30,7 +30,11 @@
// Relax the pressure equation to ensure diagonal-dominance
pEqn.relax();
pEqn.setReference(pRefCell, pRefValue);
pEqn.setReference
(
pressureControl.refCell(),
pressureControl.refValue()
);
pEqn.solve();
......@@ -60,7 +64,11 @@
fvOptions(psi, p, rho.name())
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.setReference
(
pressureControl.refCell(),
pressureControl.refValue()
);
pEqn.solve();
......@@ -81,6 +89,8 @@
U.correctBoundaryConditions();
fvOptions.correct(U);
pressureControl.limit(p);
// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
if (closedVolume)
......@@ -89,16 +99,12 @@
/fvc::domainIntegrate(psi);
}
p.correctBoundaryConditions();
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
if (!simple.transonic())
{
rho.relax();
}
Info<< "rho max/min : "
<< max(rho).value() << " "
<< min(rho).value() << endl;
}
......@@ -40,7 +40,11 @@ if (simple.transonic())
// Relax the pressure equation to maintain diagonal dominance
pEqn.relax();
pEqn.setReference(pRefCell, pRefValue);
pEqn.setReference
(
pressureControl.refCell(),
pressureControl.refValue()
);
pEqn.solve();
......@@ -75,7 +79,11 @@ else
fvOptions(psi, p, rho.name())
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.setReference
(
pressureControl.refCell(),
pressureControl.refValue()
);
pEqn.solve();
......@@ -97,6 +105,8 @@ U = HbyA - rAtU*fvc::grad(p);
U.correctBoundaryConditions();
fvOptions.correct(U);
pressureControl.limit(p);
// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
if (closedVolume)
......@@ -105,14 +115,12 @@ if (closedVolume)
/fvc::domainIntegrate(psi);
}
p.correctBoundaryConditions();
// Recalculate density from the relaxed pressure
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
if (!simple.transonic())
{
rho.relax();
}
Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
Info<< "Reading thermophysical properties\n" << endl;
autoPtr<rhoThermo> pThermo
autoPtr<fluidThermo> pThermo
(
rhoThermo::New(mesh)
fluidThermo::New(mesh)
);
rhoThermo& thermo = pThermo();
fluidThermo& thermo = pThermo();
thermo.validate(args.executable(), "h", "e");
volScalarField rho
......@@ -38,46 +38,10 @@ volVectorField U
#include "compressibleCreatePhi.H"
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, simple.dict(), pRefCell, pRefValue);
pressureControl pressureControl(p, rho, simple.dict());
mesh.setFluxRequired(p.name());
dimensionedScalar rhoMax
(
dimensionedScalar::lookupOrDefault
(
"rhoMax",
simple.dict(),
dimDensity,
GREAT
)
);
dimensionedScalar rhoMin
(
dimensionedScalar::lookupOrDefault
(
"rhoMin",
simple.dict(),
dimDensity,
0
)
);
dimensionedScalar pMin
(
dimensionedScalar::lookupOrDefault
(
"pMin",
simple.dict(),
dimPressure,
10
)
);
Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::turbulenceModel> turbulence
(
......
......@@ -48,7 +48,11 @@
fvScalarMatrix& pEqn = tpEqn.ref();
pEqn.setReference(pRefCell, pRefValue);
pEqn.setReference
(
pressureControl.refCell(),
pressureControl.refValue()
);
pEqn.solve();
......@@ -75,6 +79,8 @@
U.correctBoundaryConditions();
fvOptions.correct(U);
pressureControl.limit(p);
// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
if (closedVolume)
......@@ -83,13 +89,6 @@
/fvc::domainIntegrate(psi);
}
p = max(p, pMin);
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
rho.relax();
Info<< "rho max/min : "
<< max(rho).value() << " "
<< min(rho).value() << endl;
}
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -31,11 +31,12 @@ Description
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "rhoThermo.H"
#include "fluidThermo.H"
#include "turbulentFluidThermoModel.H"
#include "simpleControl.H"
#include "pressureControl.H"
#include "fvOptions.H"
#include "IOporosityModelList.H"
#include "simpleControl.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -30,9 +30,10 @@ Description
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "psiThermo.H"
#include "fluidThermo.H"
#include "turbulentFluidThermoModel.H"
#include "simpleControl.H"
#include "pressureControl.H"
#include "fvOptions.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......
......@@ -44,9 +44,9 @@ divSchemes
div(phi,e) $turbulence;
div(phi,h) $turbulence;
div(phi,K) $turbulence;
div(phi,Ekp) $turbulence;
div(phi,Ekp) $turbulence;
div(phid,p) bounded Gauss upwind;
div(phid,p) Gauss upwind;
div((phi|interpolate(rho)),p) bounded Gauss upwind;
div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
......
......@@ -43,8 +43,8 @@ SIMPLE
}
nNonOrthogonalCorrectors 0;
rhoMin rhoMin [ 1 -3 0 0 0 ] 0.1;
rhoMax rhoMax [ 1 -3 0 0 0 ] 1.5;
pMinFactor 0.1;
pMaxFactor 1.5;
}
relaxationFactors
......@@ -56,7 +56,7 @@ relaxationFactors
}
equations
{
U 0.7;
U 0.7;
"(e|h)" 0.7;
"(k|epsilon|omega)" 0.7;
}
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -115,6 +115,19 @@ public:
return this->transport_.mu(patchi);
}
//- Return the laminar viscosity
virtual tmp<volScalarField> nu() const
{
return this->transport_.mu()/this->rho_;
}
//- Return the laminar viscosity on patchi
virtual tmp<scalarField> nu(const label patchi) const
{
return
this->transport_.mu(patchi)/this->rho_.boundaryField()[patchi];
}
//- Return the turbulence dynamic viscosity
virtual tmp<volScalarField> mut() const
{
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -158,13 +158,13 @@ public:
}
//- Return the laminar viscosity
tmp<volScalarField> nu() const
virtual tmp<volScalarField> nu() const
{
return transport_.nu();
}
//- Return the laminar viscosity on patchi
tmp<scalarField> nu(const label patchi) const
virtual tmp<scalarField> nu(const label patchi) const
{
return transport_.nu(patchi);
}
......
......@@ -402,6 +402,7 @@ $(general)/constrainHbyA/constrainHbyA.C
$(general)/adjustPhi/adjustPhi.C
$(general)/bound/bound.C
$(general)/CorrectPhi/correctUphiBCs.C
$(general)/pressureControl/pressureControl.C
solutionControl = $(general)/solutionControl
$(solutionControl)/solutionControl/solutionControl.C
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -27,7 +27,7 @@ License
// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
void Foam::setRefCell
bool Foam::setRefCell
(
const volScalarField& field,
const volScalarField& fieldRef,
......@@ -109,11 +109,17 @@ void Foam::setRefCell
}
refValue = readScalar(dict.lookup(refValueName));
return true;
}
else
{
return false;
}
}
void Foam::setRefCell
bool Foam::setRefCell
(
const volScalarField& field,
const dictionary& dict,
......@@ -122,7 +128,7 @@ void Foam::setRefCell
const bool forceReference
)
{
setRefCell(field, field, dict, refCelli, refValue, forceReference);
return setRefCell(field, field, dict, refCelli, refValue, forceReference);
}
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -46,8 +46,8 @@ namespace Foam
//- If the field fieldRef needs referencing find the reference cell nearest
// (in index) to the given cell looked-up for field, but which is not on a
// cyclic, symmetry or processor patch.
void setRefCell
// cyclic, symmetry or processor patch and return true, otherwise return false.
bool setRefCell
(
const volScalarField& field,
const volScalarField& fieldRef,
......@@ -59,8 +59,8 @@ void setRefCell
//- If the field needs referencing find the reference cell nearest
// (in index) to the given cell looked-up for field, but which is not on a
// cyclic, symmetry or processor patch.
void setRefCell
// cyclic, symmetry or processor patch and return true, otherwise return false.
bool setRefCell
(
const volScalarField& field,
const dictionary& dict,
......
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "pressureControl.H"
#include "findRefCell.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::pressureControl::pressureControl
(
const volScalarField& p,
const volScalarField& rho,
const dictionary& dict
)
:
refCell_(0),
refValue_(0),
pMax_("pMax", dimPressure, 0),
pMin_("pMin", dimPressure, GREAT)
{
if (setRefCell(p, dict, refCell_, refValue_))
{
pMax_.value() = refValue_;
pMin_.value() = refValue_;
}
const volScalarField::Boundary& pbf = p.boundaryField();
const volScalarField::Boundary& rhobf = rho.boundaryField();
scalar rhoRefMax = -GREAT;
scalar rhoRefMin = GREAT;
bool rhoLimits = false;
forAll(pbf, patchi)
{
if (pbf[patchi].fixesValue())
{
rhoLimits = true;
pMax_.value() = max(pMax_.value(), max(pbf[patchi]));
pMin_.value() = min(pMin_.value(), min(pbf[patchi]));
rhoRefMax = max(rhoRefMax, max(rhobf[patchi]));
rhoRefMin = min(rhoRefMin, min(rhobf[patchi]));
}
}
if (dict.found("pMax"))
{
pMax_.value() = readScalar(dict.lookup("pMax"));
}
else if (dict.found("pMaxFactor"))
{
const scalar pMaxFactor(readScalar(dict.lookup("pMaxFactor")));
pMax_ *= pMaxFactor;
}
else if (dict.found("rhoMax"))
{
// For backward-compatibility infer the pMax from rhoMax
IOWarningInFunction(dict)
<< "'rhoMax' specified rather than 'pMax' or 'pMaxFactor'" << nl
<< " This is supported for backward-compatibility but "
"'pMax' or 'pMaxFactor' are more reliable." << endl;
if (!rhoLimits)
{
FatalIOErrorInFunction(dict)
<< "'rhoMax' specified rather than 'pMaxFactor'" << nl
<< " but the corresponding reference density cannot"
" be evaluated from the boundary conditions." << nl
<< "Please specify 'pMaxFactor' rather than 'rhoMax'"
<< exit(FatalError);
}
dimensionedScalar rhoMax("rhoMax", dimDensity, dict);
pMax_ *= max(rhoMax.value()/rhoRefMax, 1);
}
if (dict.found("pMin"))
{
pMin_.value() = readScalar(dict.lookup("pMin"));
}
else if (dict.found("pMinFactor"))
{
const scalar pMinFactor(readScalar(dict.lookup("pMinFactor")));
pMin_ *= pMinFactor;
}
else if (dict.found("rhoMin"))
{
// For backward-compatibility infer the pMin from rhoMin
IOWarningInFunction(dict)
<< "'rhoMin' specified rather than 'pMin' or 'pMinFactor'" << nl
<< " This is supported for backward-compatibility but"
"'pMin' or 'pMinFactor' are more reliable." << endl;
if (!rhoLimits)
{
FatalIOErrorInFunction(dict)