buoyantReactingFoam: Added optional hydrostatic initialisation and replaced fireFoam

The fireFoam solver has solver has been replaced by the more general
buoyantReactingFoam solver, which supports buoyant compressible reacting flow
coupled to multiple run-time-selectable lagrangian clouds and surface film
modelling and optional hydrostatic initialisation of the pressure and p_rgh.

Hydrostatic initialisation of the pressure fields is useful for large fires in
open domains where the stability of the initial flow is dominated by the initial
pressure distribution in the domain and at the boundaries.  The optional
hydrostaticInitialization switch in fvSolution/PIMPLE with
nHydrostaticCorrectors enables hydrostatic initialisation, e.g.

PIMPLE
{
    momentumPredictor yes;
    nOuterCorrectors  1;
    nCorrectors       2;
    nNonOrthogonalCorrectors 0;

    hydrostaticInitialization yes;
    nHydrostaticCorrectors 5;
}

and the resulting ph_rgh field can be used with the prghTotalHydrostaticPressure
p_rgh boundary condition to apply this hydrostatic pressure distribution at the
boundaries throughout the simulation.

See the following cases for examples transferred from fireFoam:

    $FOAM_TUTORIALS/combustion/buoyantReactingFoam/RAS
This commit is contained in:
Henry Weller
2021-05-31 15:05:19 +01:00
parent 49ce8f6507
commit a997ddae5f
81 changed files with 56 additions and 593 deletions

View File

@ -1,3 +0,0 @@
fireFoam.C
EXE = $(FOAM_APPBIN)/fireFoam

View File

@ -1,33 +0,0 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
-I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
-I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude \
-I$(LIB_SRC)/ThermophysicalTransportModels/lnInclude \
-I$(LIB_SRC)/ThermophysicalTransportModels/fluidReactionThermo/lnInclude \
-I$(LIB_SRC)/transportModels/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
-I$(LIB_SRC)/combustionModels/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
-I$(LIB_SRC)/ODE/lnInclude
EXE_LIBS = \
-lfiniteVolume \
-lfvModels \
-lfvConstraints \
-lmeshTools \
-lsampling \
-lmomentumTransportModels \
-lcompressibleMomentumTransportModels \
-lthermophysicalTransportModels \
-lfluidReactionThermophysicalTransportModels \
-lspecie \
-lfluidThermophysicalModels \
-lreactionThermophysicalModels \
-lchemistryModel \
-lcombustionModels \
-lODE

View File

@ -1,33 +0,0 @@
MRF.correctBoundaryVelocity(U);
fvVectorMatrix UEqn
(
fvm::ddt(rho, U) + fvm::div(phi, U)
+ MRF.DDt(rho, U)
+ turbulence->divDevTau(U)
==
fvModels.source(rho, U)
);
UEqn.relax();
fvConstraints.constrain(UEqn);
if (pimple.momentumPredictor())
{
solve
(
UEqn
==
fvc::reconstruct
(
(
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
)*mesh.magSf()
)
);
fvConstraints.constrain(U);
K = 0.5*magSqr(U);
}

View File

@ -1,76 +0,0 @@
tmp<fv::convectionScheme<scalar>> mvConvection
(
fv::convectionScheme<scalar>::New
(
mesh,
fields,
phi,
mesh.divScheme("div(phi,Yi_h)")
)
);
{
combustion->correct();
forAll(Y, i)
{
if (composition.solve(i))
{
volScalarField& Yi = Y[i];
fvScalarMatrix YiEqn
(
fvm::ddt(rho, Yi)
+ mvConvection->fvmDiv(phi, Yi)
+ thermophysicalTransport->divj(Yi)
==
combustion->R(Yi)
+ fvModels.source(rho, Yi)
);
YiEqn.relax();
fvConstraints.constrain(YiEqn);
YiEqn.solve("Yi");
fvConstraints.constrain(Yi);
}
}
composition.normalise();
volScalarField& he = thermo.he();
fvScalarMatrix EEqn
(
fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he)
+ fvc::ddt(rho, K) + fvc::div(phi, K)
+ (
he.name() == "e"
? fvc::div
(
fvc::absolute(phi/fvc::interpolate(rho), U),
p,
"div(phiv,p)"
)
: -dpdt
)
+ thermophysicalTransport->divq(he)
==
combustion->Qdot()
+ fvModels.source(rho, he)
);
EEqn.relax();
fvConstraints.constrain(EEqn);
EEqn.solve();
fvConstraints.constrain(he);
thermo.correct();
Info<< "min/max(T) = "
<< min(T).value() << ", " << max(T).value() << endl;
}

View File

@ -1,2 +0,0 @@
const volScalarField& psi = thermo.psi();
const volScalarField& T = thermo.T();

View File

@ -1,117 +0,0 @@
Info<< "Reading thermophysical properties\n" << endl;
autoPtr<fluidReactionThermo> pThermo(fluidReactionThermo::New(mesh));
fluidReactionThermo& thermo = pThermo();
thermo.validate(args.executable(), "h", "e");
basicSpecieMixture& composition = thermo.composition();
PtrList<volScalarField>& Y = composition.Y();
Info<< "Creating field rho\n" << endl;
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
thermo.rho()
);
volScalarField& p = thermo.p();
Info<< "\nReading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
#include "compressibleCreatePhi.H"
#include "createMRF.H"
Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::momentumTransportModel> turbulence
(
compressible::momentumTransportModel::New
(
rho,
U,
phi,
thermo
)
);
Info<< "Creating thermophysical transport model\n" << endl;
autoPtr<fluidReactionThermophysicalTransportModel> thermophysicalTransport
(
fluidReactionThermophysicalTransportModel::New(turbulence(), thermo)
);
Info<< "Creating combustion model\n" << endl;
autoPtr<combustionModel> combustion
(
combustionModel::New(thermo, turbulence())
);
#include "readGravitationalAcceleration.H"
#include "readhRef.H"
#include "gh.H"
#include "readpRef.H"
volScalarField p_rgh
(
IOobject
(
"p_rgh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
mesh.setFluxRequired(p_rgh.name());
#include "phrghEqn.H"
multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
forAll(Y, i)
{
fields.add(Y[i]);
}
fields.add(thermo.he());
Info<< "Creating field dpdt\n" << endl;
volScalarField dpdt
(
IOobject
(
"dpdt",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar(p.dimensions()/dimTime, 0)
);
Info<< "Creating field kinetic energy K\n" << endl;
volScalarField K("K", 0.5*magSqr(U));
#include "createFvModels.H"
#include "createFvConstraints.H"

View File

@ -1,139 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2021 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/>.
Application
fireFoam
Description
Transient solver for fires and turbulent diffusion flames with optional
reacting particle clouds and surface film modelling via fvModels.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "dynamicMomentumTransportModel.H"
#include "fluidReactionThermophysicalTransportModel.H"
#include "fluidReactionThermo.H"
#include "combustionModel.H"
#include "pimpleControl.H"
#include "fvModels.H"
#include "fvConstraints.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "postProcess.H"
#include "setRootCaseLists.H"
#include "createTime.H"
#include "createMesh.H"
#include "createControl.H"
#include "createFields.H"
#include "createFieldRefs.H"
#include "initContinuityErrs.H"
#include "createTimeControls.H"
#include "compressibleCourantNo.H"
#include "setInitialDeltaT.H"
turbulence->validate();
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (pimple.run(runTime))
{
#include "readTimeControls.H"
#include "compressibleCourantNo.H"
#include "setMultiRegionDeltaT.H"
#include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
// --- PIMPLE loop
while (pimple.loop())
{
if (!pimple.flow())
{
if (pimple.models())
{
fvModels.correct();
}
if (pimple.thermophysics())
{
#include "YEEqn.H"
}
}
else
{
if (pimple.firstPimpleIter() && !pimple.simpleRho())
{
#include "rhoEqn.H"
}
if (pimple.models())
{
fvModels.correct();
}
#include "UEqn.H"
if (pimple.thermophysics())
{
#include "YEEqn.H"
}
// --- Pressure corrector loop
while (pimple.correct())
{
#include "pEqn.H"
}
if (pimple.turbCorr())
{
turbulence->correct();
thermophysicalTransport->correct();
}
}
rho = thermo.rho();
}
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -1,58 +0,0 @@
rho = thermo.rho();
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
surfaceScalarField phig("phig", -rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
surfaceScalarField phiHbyA
(
"phiHbyA",
(
fvc::flux(rho*HbyA)
+ MRF.zeroFilter(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);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
fvm::ddt(psi, p_rgh)
+ fvc::ddt(psi, rho)*gh
+ fvc::ddt(psi)*pRef
+ fvc::div(phiHbyA)
- fvm::laplacian(rhorAUf, p_rgh)
==
fvModels.source(psi, p_rgh, rho.name())
);
p_rghEqn.solve();
if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA + p_rghEqn.flux();
U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/rhorAUf);
U.correctBoundaryConditions();
fvConstraints.constrain(U);
}
}
p = p_rgh + rho*gh + pRef;
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
K = 0.5*magSqr(U);
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);
}

View File

@ -1,62 +0,0 @@
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;
}
}

View File

@ -1,46 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2021 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/>.
Global
rhoEqn
Description
Solve the continuity for density.
\*---------------------------------------------------------------------------*/
{
fvScalarMatrix rhoEqn
(
fvm::ddt(rho)
+ fvc::div(phi)
==
fvModels.source(rho)
);
rhoEqn.solve();
fvConstraints.constrain(rho);
}
// ************************************************************************* //

View File

@ -1,56 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2021 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/>.
Global
setMultiRegionDeltaT
Description
Reset the timestep to maintain a constant maximum Courant numbers.
Reduction of time-step is immediate, but increase is damped to avoid
unstable oscillations.
\*---------------------------------------------------------------------------*/
if (adjustTimeStep)
{
if (CoNum == -great)
{
CoNum = small;
}
const scalar TFactorFluid = maxCo/(CoNum + small);
const scalar dt0 = runTime.deltaTValue();
runTime.setDeltaT
(
min
(
dt0*min(TFactorFluid, 1.2),
maxDeltaT
)
);
}
// ************************************************************************* //

View File

@ -86,9 +86,6 @@ volScalarField p_rgh
mesh
);
// Force p_rgh to be consistent with p
p_rgh = p - rho*gh - pRef;
pressureReference pressureReference
(
p,
@ -99,6 +96,8 @@ pressureReference pressureReference
mesh.setFluxRequired(p_rgh.name());
#include "hydrostaticInitialisation.H"
Info<< "Creating field dpdt\n" << endl;
volScalarField dpdt
(