Commit 2c6b4050 authored by Henry Weller's avatar Henry Weller
Browse files

fireFoam: Added optional hydrostatic initialization of the pressure and density

Also added the new prghTotalHydrostaticPressure p_rgh BC which uses the
hydrostatic pressure field as the reference state for the far-field
which provides much more accurate entrainment is large open domains
typical of many fire simulations.

The hydrostatic field solution is controlled by the optional entries in
the fvSolution.PIMPLE dictionary, e.g.

    hydrostaticInitialization yes;
    nHydrostaticCorrectors 5;

and the solver must also be specified for the hydrostatic p_rgh field
ph_rgh e.g.

    ph_rgh
    {
        $p_rgh;
    }

Suitable boundary conditions for ph_rgh cannot always be derived from
those for p_rgh and so the ph_rgh is read to provide them.

To avoid accuracy issues with IO, restart and post-processing the p_rgh
and ph_rgh the option to specify a suitable reference pressure is
provided via the optional pRef file in the constant directory, e.g.

    dimensions      [1 -1 -2 0 0 0 0];
    value           101325;

which is used in the relationship between p_rgh and p:

    p = p_rgh + rho*gh + pRef;

Note that if pRef is specified all pressure BC specifications in the
p_rgh and ph_rgh files are relative to the reference to avoid round-off
errors.

For examples of suitable BCs for p_rgh and ph_rgh for a range of
fireFoam cases please study the tutorials in
tutorials/combustion/fireFoam/les which have all been updated.

Henry G. Weller
CFD Direct Ltd.
parent c79f0607
......@@ -54,6 +54,9 @@ volVectorField U
#include "compressibleCreatePhi.H"
#include "createMRF.H"
Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::turbulenceModel> turbulence
(
......@@ -69,42 +72,11 @@ autoPtr<compressible::turbulenceModel> turbulence
// Set the turbulence into the combustion model
combustion->setTurbulence(turbulence());
volScalarField dQ
(
IOobject
(
"dQ",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
);
Info<< "Creating field dpdt\n" << endl;
volScalarField dpdt
(
IOobject
(
"dpdt",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
);
Info<< "Creating field kinetic energy K\n" << endl;
volScalarField K("K", 0.5*magSqr(U));
#include "readGravitationalAcceleration.H"
#include "readhRef.H"
#include "gh.H"
#include "readpRef.H"
volScalarField p_rgh
(
......@@ -119,11 +91,11 @@ volScalarField p_rgh
mesh
);
// Force p_rgh to be consistent with p
p_rgh = p - rho*gh;
mesh.setFluxRequired(p_rgh.name());
#include "phrghEqn.H"
multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
forAll(Y, i)
......@@ -148,3 +120,34 @@ Switch solvePrimaryRegion
(
additionalControlsDict.lookup("solvePrimaryRegion")
);
volScalarField dQ
(
IOobject
(
"dQ",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
);
Info<< "Creating field dpdt\n" << endl;
volScalarField dpdt
(
IOobject
(
"dpdt",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
);
Info<< "Creating field kinetic energy K\n" << endl;
volScalarField K("K", 0.5*magSqr(U));
......@@ -54,7 +54,6 @@ int main(int argc, char *argv[])
pimpleControl pimple(mesh);
#include "createFields.H"
#include "createMRF.H"
#include "createFvOptions.H"
#include "createClouds.H"
#include "createSurfaceFilmModel.H"
......
......@@ -4,7 +4,7 @@ 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("phig", -rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
surfaceScalarField phiHbyA
(
......@@ -25,9 +25,10 @@ while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
fvc::ddt(psi, rho)*gh
fvm::ddt(psi, p_rgh)
+ fvc::ddt(psi, rho)*gh
+ fvc::ddt(psi)*pRef
+ fvc::div(phiHbyA)
+ fvm::ddt(psi, p_rgh)
- fvm::laplacian(rhorAUf, p_rgh)
==
parcels.Srho()
......@@ -46,7 +47,7 @@ while (pimple.correctNonOrthogonal())
}
}
p = p_rgh + rho*gh;
p = p_rgh + rho*gh + pRef;
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
......
if (pimple.dict().lookupOrDefault<bool>("hydrostaticInitialization", false))
{
volScalarField& ph_rgh = regIOobject::store
(
new volScalarField
(
IOobject
(
"ph_rgh",
"0",
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
),
mesh
)
);
if (equal(runTime.value(), 0))
{
p = ph_rgh + rho*gh + pRef;
thermo.correct();
rho = thermo.rho();
label nCorr
(
pimple.dict().lookupOrDefault<label>("nHydrostaticCorrectors", 5)
);
for (label i=0; i<nCorr; i++)
{
surfaceScalarField rhof("rhof", fvc::interpolate(rho));
surfaceScalarField phig
(
"phig",
-rhof*ghf*fvc::snGrad(rho)*mesh.magSf()
);
// Update the pressure BCs to ensure flux consistency
constrainPressure(ph_rgh, rho, U, phig, rhof);
fvScalarMatrix ph_rghEqn
(
fvm::laplacian(rhof, ph_rgh) == fvc::div(phig)
);
ph_rghEqn.solve();
p = ph_rgh + rho*gh + pRef;
thermo.correct();
rho = thermo.rho();
Info<< "Hydrostatic pressure variation "
<< (max(ph_rgh) - min(ph_rgh)).value() << endl;
}
ph_rgh.write();
p_rgh = ph_rgh;
}
}
......@@ -202,6 +202,7 @@ $(derivedFvPatchFields)/waveSurfacePressure/waveSurfacePressureFvPatchScalarFiel
$(derivedFvPatchFields)/interstitialInletVelocity/interstitialInletVelocityFvPatchVectorField.C
$(derivedFvPatchFields)/prghPressure/prghPressureFvPatchScalarField.C
$(derivedFvPatchFields)/prghTotalPressure/prghTotalPressureFvPatchScalarField.C
$(derivedFvPatchFields)/prghTotalHydrostaticPressure/prghTotalHydrostaticPressureFvPatchScalarField.C
$(derivedFvPatchFields)/fixedProfile/fixedProfileFvPatchFields.C
fvsPatchFields = fields/fvsPatchFields
......
Info<< "\nReading pRef" << endl;
uniformDimensionedScalarField pRef
(
IOobject
(
"pRef",
runTime.constant(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
),
dimensionedScalar("pRef", dimPressure, 0)
);
......@@ -43,7 +43,7 @@ Description
myPatch
{
type fixedValue;
value uniform 0; // example for scalar field usage
value uniform 0; // Example for scalar field usage
}
\endverbatim
......
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016 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 "prghTotalHydrostaticPressureFvPatchScalarField.H"
#include "addToRunTimeSelectionTable.H"
#include "fvPatchFieldMapper.H"
#include "volFields.H"
#include "surfaceFields.H"
#include "uniformDimensionedFields.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::prghTotalHydrostaticPressureFvPatchScalarField::
prghTotalHydrostaticPressureFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF
)
:
fixedValueFvPatchScalarField(p, iF),
UName_("U"),
phiName_("phi"),
rhoName_("rho"),
ph_rghName_("ph_rgh")
{}
Foam::prghTotalHydrostaticPressureFvPatchScalarField::
prghTotalHydrostaticPressureFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
const dictionary& dict
)
:
fixedValueFvPatchScalarField(p, iF, dict),
UName_(dict.lookupOrDefault<word>("U", "U")),
phiName_(dict.lookupOrDefault<word>("phi", "phi")),
rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
ph_rghName_(dict.lookupOrDefault<word>("ph_rgh", "ph_rgh"))
{}
Foam::prghTotalHydrostaticPressureFvPatchScalarField::
prghTotalHydrostaticPressureFvPatchScalarField
(
const prghTotalHydrostaticPressureFvPatchScalarField& ptf,
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
const fvPatchFieldMapper& mapper
)
:
fixedValueFvPatchScalarField(ptf, p, iF, mapper),
UName_(ptf.UName_),
phiName_(ptf.phiName_),
rhoName_(ptf.rhoName_),
ph_rghName_(ptf.ph_rghName_)
{}
Foam::prghTotalHydrostaticPressureFvPatchScalarField::
prghTotalHydrostaticPressureFvPatchScalarField
(
const prghTotalHydrostaticPressureFvPatchScalarField& ptf
)
:
fixedValueFvPatchScalarField(ptf),
UName_(ptf.UName_),
phiName_(ptf.phiName_),
rhoName_(ptf.rhoName_),
ph_rghName_(ptf.ph_rghName_)
{}
Foam::prghTotalHydrostaticPressureFvPatchScalarField::
prghTotalHydrostaticPressureFvPatchScalarField
(
const prghTotalHydrostaticPressureFvPatchScalarField& ptf,
const DimensionedField<scalar, volMesh>& iF
)
:
fixedValueFvPatchScalarField(ptf, iF),
UName_(ptf.UName_),
phiName_(ptf.phiName_),
rhoName_(ptf.rhoName_),
ph_rghName_(ptf.ph_rghName_)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::prghTotalHydrostaticPressureFvPatchScalarField::updateCoeffs()
{
if (updated())
{
return;
}
const scalarField& rhop =
patch().lookupPatchField<volScalarField, scalar>(rhoName_);
const scalarField& ph_rghp =
patch().lookupPatchField<volScalarField, scalar>(ph_rghName_);
const scalarField& phip =
patch().lookupPatchField<surfaceScalarField, scalar>(phiName_);
const vectorField& Up =
patch().lookupPatchField<volVectorField, vector>(UName_);
operator==
(
ph_rghp
- 0.5*rhop*(1.0 - pos(phip))*magSqr(Up)
);
fixedValueFvPatchScalarField::updateCoeffs();
}
void Foam::prghTotalHydrostaticPressureFvPatchScalarField::write
(
Ostream& os
) const
{
fvPatchScalarField::write(os);
writeEntryIfDifferent<word>(os, "U", "U", UName_);
writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
writeEntryIfDifferent<word>(os, "ph_rghName", "ph_rghName", ph_rghName_);
writeEntry("value", os);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
makePatchTypeField
(
fvPatchScalarField,
prghTotalHydrostaticPressureFvPatchScalarField
);
}
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::prghTotalHydrostaticPressureFvPatchScalarField
Group
grpGenericBoundaryConditions
Description
This boundary condition provides static pressure condition for p_rgh,
calculated as:
\f[
p_rgh = ph_rgh - 0.5 \rho |U|^2
\f]
where
\vartable
p_rgh | Pressure - \rho g.(h - hRef) [Pa]
ph_rgh | Hydrostatic pressure - \rho g.(h - hRef) [Pa]
h | Height in the opposite direction to gravity
hRef | Reference height in the opposite direction to gravity
\rho | Density
g | Acceleration due to gravity [m/s^2]
\endtable
\heading Patch usage
\table
Property | Description | Required | Default value
U | Velocity field name | no | U
phi | Flux field name | no | phi
rho | Density field name | no | rho
ph_rgh | ph_rgh field name | no | ph_rgh
value | Patch face values | yes |
\endtable
Example of the boundary condition specification:
\verbatim
myPatch
{
type prghTotalHydrostaticPressure;
value uniform 0;
}
\endverbatim
SeeAlso
Foam::fixedValueFvPatchScalarField
Foam::prghTotalPressureFvPatchScalarField
SourceFiles
prghTotalHydrostaticPressureFvPatchScalarField.C
\*---------------------------------------------------------------------------*/
#ifndef prghTotalHydrostaticPressureFvPatchScalarField_H
#define prghTotalHydrostaticPressureFvPatchScalarField_H
#include "fixedValueFvPatchFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class prghTotalHydrostaticPressureFvPatchScalarField Declaration
\*---------------------------------------------------------------------------*/
class prghTotalHydrostaticPressureFvPatchScalarField
:
public fixedValueFvPatchScalarField
{
protected:
// Protected data
//- Name of the velocity field
word UName_;
//- Name of the flux transporting the field
word phiName_;
//- Name of density field
word rhoName_;
//- Name of hydrostatic pressure field
word ph_rghName_;
public:
//- Runtime type information
TypeName("prghTotalHydrostaticPressure");
// Constructors
//- Construct from patch and internal field
prghTotalHydrostaticPressureFvPatchScalarField
(
const fvPatch&,
const DimensionedField<scalar, volMesh>&
);
//- Construct from patch, internal field and dictionary
prghTotalHydrostaticPressureFvPatchScalarField
(
const fvPatch&,
const DimensionedField<scalar, volMesh>&,
const dictionary&
);
//- Construct by mapping given
// prghTotalHydrostaticPressureFvPatchScalarField onto a new patch
prghTotalHydrostaticPressureFvPatchScalarField
(
const prghTotalHydrostaticPressureFvPatchScalarField&,
const fvPatch&,
const DimensionedField<scalar, volMesh>&,
const fvPatchFieldMapper&
);
//- Construct as copy
prghTotalHydrostaticPressureFvPatchScalarField
(
const prghTotalHydrostaticPressureFvPatchScalarField&
);
//- Construct and return a clone
virtual tmp<fvPatchScalarField> clone() const
{
return tmp<fvPatchScalarField >
(
new prghTotalHydrostaticPressureFvPatchScalarField(*this)
);
}
//- Construct as copy setting internal field reference
prghTotalHydrostaticPressureFvPatchScalarField
(
const prghTotalHydrostaticPressureFvPatchScalarField&,
const DimensionedField<scalar, volMesh>&
);
//- Construct and return a clone setting internal field reference
virtual tmp