Commit 791e1ca2 authored by Henry Weller's avatar Henry Weller Committed by Andrew Heather
Browse files

Merged reactingParcelFilmFoam into reactingParcelFoam

The combined solver includes the most advanced and general functionality from
each solver including:

    Continuous phase
    Lagrangian multiphase parcels
    Optional film
    Continuous and Lagrangian phase reactions
    Radiation
    Strong buoyancy force support by solving for p_rgh

The reactingParcelFoam and reactingParcelFilmFoam tutorials have been combined
and updated.
parent 5cb3bb7b
const volScalarField& psi = thermo.psi();
const volScalarField& T = thermo.T();
filmModelType& surfaceFilm = tsurfaceFilm();
regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm();
const label inertIndex(composition.species()[inertSpecie]);
Info<< "\nConstructing surface film model" << endl;
typedef regionModels::surfaceFilmModels::surfaceFilmModel filmModelType;
autoPtr<filmModelType> tsurfaceFilm(filmModelType::New(mesh, g));
autoPtr<regionModels::surfaceFilmModel> tsurfaceFilm
(
regionModels::surfaceFilmModel::New(mesh, g)
);
reactingParcelFilmFoam.C
EXE = $(FOAM_APPBIN)/reactingParcelFilmFoam
EXE_INC = \
-I. \
-I../reactingParcelFoam \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I${LIB_SRC}/sampling/lnInclude \
-I${LIB_SRC}/meshTools/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-I$(LIB_SRC)/transportModels/compressible/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
-I$(LIB_SRC)/regionModels/regionModel/lnInclude \
-I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
-I$(LIB_SRC)/lagrangian/basic/lnInclude \
-I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
-I$(LIB_SRC)/ODE/lnInclude \
-I$(LIB_SRC)/combustionModels/lnInclude \
-I$(FOAM_SOLVERS)/combustion/reactingFoam
EXE_LIBS = \
-lfiniteVolume \
-lfvOptions \
-lsampling \
-lmeshTools \
-lturbulenceModels \
-lcompressibleTurbulenceModels \
-lspecie \
-lcompressibleTransportModels \
-lfluidThermophysicalModels \
-lreactionThermophysicalModels \
-lSLGThermo \
-lchemistryModel \
-lregionModels \
-lradiationModels \
-lsurfaceFilmModels \
-lsurfaceFilmDerivedFvPatchFields \
-llagrangianIntermediate \
-llagrangianTurbulence \
-lODE \
-lcombustionModels
MRF.correctBoundaryVelocity(U);
fvVectorMatrix UEqn
(
fvm::ddt(rho, U) + fvm::div(phi, U)
+ MRF.DDt(rho, U)
+ turbulence->divDevRhoReff(U)
==
parcels.SU(U)
+ fvOptions(rho, U)
);
UEqn.relax();
fvOptions.constrain(UEqn);
if (pimple.momentumPredictor())
{
solve
(
UEqn
==
fvc::reconstruct
(
(
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
)*mesh.magSf()
)
);
fvOptions.correct(U);
K = 0.5*magSqr(U);
}
Info<< "\nConstructing reacting cloud" << endl;
basicReactingCloud parcels
(
"reactingCloud1",
rho,
U,
g,
slgThermo
);
#include "createRDeltaT.H"
Info<< "Creating combustion model\n" << endl;
autoPtr<combustionModels::psiCombustionModel> combustion
(
combustionModels::psiCombustionModel::New(mesh)
);
psiReactionThermo& thermo = combustion->thermo();
thermo.validate(args.executable(), "h", "e");
SLGThermo slgThermo(mesh, thermo);
basicSpecieMixture& composition = thermo.composition();
PtrList<volScalarField>& Y = composition.Y();
const word inertSpecie(thermo.lookup("inertSpecie"));
if (!composition.species().found(inertSpecie))
{
FatalIOErrorIn(args.executable().c_str(), thermo)
<< "Inert specie " << inertSpecie << " not found in available species "
<< composition.species()
<< exit(FatalIOError);
}
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"
Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::turbulenceModel> turbulence
(
compressible::turbulenceModel::New
(
rho,
U,
phi,
thermo
)
);
// Set the turbulence into the combustion model
combustion->setTurbulence(turbulence());
#include "readGravitationalAcceleration.H"
#include "readhRef.H"
#include "gh.H"
volScalarField p_rgh
(
IOobject
(
"p_rgh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
// Force p_rgh to be consistent with p
p_rgh = p - rho*gh;
mesh.setFluxRequired(p_rgh.name());
multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
forAll(Y, i)
{
fields.add(Y[i]);
}
fields.add(thermo.he());
IOdictionary additionalControlsDict
(
IOobject
(
"additionalControls",
runTime.constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
Switch solvePrimaryRegion
(
additionalControlsDict.lookup("solvePrimaryRegion")
);
volScalarField Qdot
(
IOobject
(
"Qdot",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
);
#include "createDpdt.H"
#include "createK.H"
#include "createMRF.H"
#include "createClouds.H"
#include "createRadiationModel.H"
#include "createSurfaceFilmModel.H"
Info<< "\nConstructing surface film model" << endl;
typedef regionModels::surfaceFilmModels::surfaceFilmModel filmModelType;
autoPtr<filmModelType> tsurfaceFilm(filmModelType::New(mesh, g));
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(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
surfaceScalarField 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);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
fvc::ddt(psi, rho)*gh
+ fvc::div(phiHbyA)
+ fvm::ddt(psi, p_rgh)
- fvm::laplacian(rhorAUf, p_rgh)
==
parcels.Srho()
+ surfaceFilm.Srho()
+ fvOptions(psi, p_rgh, rho.name())
);
p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA + p_rghEqn.flux();
U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/rhorAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
}
}
p = p_rgh + rho*gh;
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
K = 0.5*magSqr(U);
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);
}
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-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/>.
Application
reactingParcelFilmFoam
Group
grpLagrangianSolvers
Description
Transient solver for compressible, turbulent flow with a reacting,
multiphase particle cloud, and surface film modelling.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "turbulentFluidThermoModel.H"
#include "basicReactingCloud.H"
#include "surfaceFilmModel.H"
#include "psiCombustionModel.H"
#include "radiationModel.H"
#include "SLGThermo.H"
#include "fvOptions.H"
#include "pimpleControl.H"
#include "localEulerDdtScheme.H"
#include "fvcSmooth.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "postProcess.H"
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "createControl.H"
#include "createTimeControls.H"
#include "createFields.H"
#include "createFieldRefs.H"
#include "createFvOptions.H"
#include "initContinuityErrs.H"
turbulence->validate();
if (!LTS)
{
#include "compressibleCourantNo.H"
#include "setInitialDeltaT.H"
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
#include "readTimeControls.H"
if (!LTS)
{
#include "compressibleCourantNo.H"
#include "setMultiRegionDeltaT.H"
#include "setDeltaT.H"
}
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
parcels.evolve();
surfaceFilm.evolve();
if (LTS)
{
#include "setRDeltaT.H"
}
if (solvePrimaryRegion)
{
#include "rhoEqn.H"
// --- PIMPLE loop
while (pimple.loop())
{
#include "UEqn.H"
#include "YEqn.H"
#include "EEqn.H"
// --- Pressure corrector loop
while (pimple.correct())
{
#include "pEqn.H"
}
if (pimple.turbCorr())
{
turbulence->correct();
}
}
rho = thermo.rho();
}
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End" << endl;
return 0;
}
// ************************************************************************* //
......@@ -19,6 +19,7 @@
==
rho*(U&g)
+ parcels.Sh(he)
+ surfaceFilm.Sh()
+ radiation->Sh(thermo, he)
+ Qdot
+ fvOptions(rho, he)
......@@ -35,6 +36,6 @@
thermo.correct();
radiation->correct();
Info<< "T gas min/max " << min(T).value() << ", "
Info<< "T gas min/max = " << min(T).value() << ", "
<< max(T).value() << endl;
}
EXE_INC = \
-I. \
-I../reactingParcelFoam \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I${LIB_SRC}/sampling/lnInclude \
-I${LIB_SRC}/meshTools/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-I$(LIB_SRC)/lagrangian/basic/lnInclude \
-I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
-I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \
-I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-I$(LIB_SRC)/transportModels/compressible/lnInclude \
......@@ -17,33 +16,33 @@ EXE_INC = \
-I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
-I$(LIB_SRC)/ODE/lnInclude \
-I$(LIB_SRC)/regionModels/regionModel/lnInclude \
-I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
-I$(LIB_SRC)/lagrangian/basic/lnInclude \
-I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
-I$(LIB_SRC)/ODE/lnInclude \
-I$(LIB_SRC)/combustionModels/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
-I$(FOAM_SOLVERS)/combustion/reactingFoam
EXE_LIBS = \
-lfiniteVolume \
-lfvOptions \
-lsampling \
-lmeshTools \
-lturbulenceModels \
-lcompressibleTurbulenceModels \
-llagrangian \
-llagrangianIntermediate \
-llagrangianTurbulence \
-lspecie \
-lcompressibleTransportModels \
-lfluidThermophysicalModels \
-lthermophysicalProperties \
-lreactionThermophysicalModels \
-lSLGThermo \
-lchemistryModel \
-lradiationModels \
-lODE \
-lregionModels \
-lradiationModels \
-lsurfaceFilmModels \
-lcombustionModels \
-lfvOptions \
-lsampling
-lsurfaceFilmDerivedFvPatchFields \
-llagrangian \
-llagrangianIntermediate \
-llagrangianTurbulence \
-lODE \
-lcombustionModels
......@@ -6,8 +6,7 @@
+ MRF.DDt(rho, U)
+ turbulence->divDevRhoReff(U)
==
rho()*g
+ parcels.SU(U)
parcels.SU(U)
+ fvOptions(rho, U)
);
......@@ -17,7 +16,18 @@
if (pimple.momentumPredictor())
{
solve(UEqn == -fvc::grad(p));
solve
(
UEqn
==
fvc::reconstruct
(
(
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
)*mesh.magSf()
)
);
fvOptions.correct(U);
K = 0.5*magSqr(U);
......
......@@ -9,6 +9,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
)
);
{
combustion->correct();
Qdot = combustion->Qdot();
......@@ -24,11 +25,12 @@ tmp<fv::convectionScheme<scalar>> mvConvection
(
fvm::ddt(rho, Yi)
+ mvConvection->fvmDiv(phi, Yi)
- fvm::laplacian(turbulence->muEff(), Yi)
==
- fvm::laplacian(turbulence->alphaEff(), Yi)
==
parcels.SYi(i, Yi)
+ combustion->R(Yi)
+ fvOptions(rho, Yi)
+ combustion->R(Yi)
+ surfaceFilm.Srho(i)
);
YEqn.relax();
......