Solvers: Updated to solve EEqn rather than h,hs,eEqn

This commit is contained in:
Henry
2012-09-22 16:34:20 +01:00
parent 790409ae2b
commit 78121bac4a
215 changed files with 13149 additions and 1094 deletions

View File

@ -5,6 +5,8 @@
psiuReactionThermo::New(mesh) psiuReactionThermo::New(mesh)
); );
psiuReactionThermo& thermo = pThermo(); psiuReactionThermo& thermo = pThermo();
thermo.validate(args.executable(), "ha");
basicMultiComponentMixture& composition = thermo.composition(); basicMultiComponentMixture& composition = thermo.composition();
volScalarField rho volScalarField rho

View File

@ -5,6 +5,8 @@
psiuReactionThermo::New(mesh) psiuReactionThermo::New(mesh)
); );
psiuReactionThermo& thermo = pThermo(); psiuReactionThermo& thermo = pThermo();
thermo.validate(args.executable(), "ha");
basicMultiComponentMixture& composition = thermo.composition(); basicMultiComponentMixture& composition = thermo.composition();
volScalarField rho volScalarField rho

View File

@ -30,6 +30,8 @@
scalar dtChem = refCast<const psiChemistryModel>(chemistry).deltaTChem()[0]; scalar dtChem = refCast<const psiChemistryModel>(chemistry).deltaTChem()[0];
psiReactionThermo& thermo = chemistry.thermo(); psiReactionThermo& thermo = chemistry.thermo();
thermo.validate(args.executable(), "h");
basicMultiComponentMixture& composition = thermo.composition(); basicMultiComponentMixture& composition = thermo.composition();
PtrList<volScalarField>& Y = composition.Y(); PtrList<volScalarField>& Y = composition.Y();
@ -47,7 +49,6 @@
); );
volScalarField& p = thermo.p(); volScalarField& p = thermo.p();
volScalarField& hs = thermo.he();
volScalarField Rspecific volScalarField Rspecific
( (

View File

@ -1,10 +1,12 @@
{ {
volScalarField& h = thermo.he();
if (constProp == "volume") if (constProp == "volume")
{ {
hs[0] = u0 + p[0]/rho[0] + integratedHeat; h[0] = u0 + p[0]/rho[0] + integratedHeat;
} }
else else
{ {
hs[0] = hs0 + integratedHeat; h[0] = h0 + integratedHeat;
} }
} }

View File

@ -83,20 +83,19 @@
} }
} }
scalar hs0 = 0.0; scalar h0 = 0.0;
forAll(Y, i) forAll(Y, i)
{ {
Y[i] = Y0[i]; Y[i] = Y0[i];
hs0 += Y0[i]*specieData[i].Hs(p[i], T0); h0 += Y0[i]*specieData[i].Hs(p[i], T0);
} }
hs = dimensionedScalar("h", dimEnergy/dimMass, hs0); thermo.he() = dimensionedScalar("h", dimEnergy/dimMass, h0);
thermo.correct(); thermo.correct();
rho = thermo.rho(); rho = thermo.rho();
scalar rho0 = rho[0]; scalar rho0 = rho[0];
scalar u0 = hs0 - p0/rho0; scalar u0 = h0 - p0/rho0;
scalar R0 = p0/(rho0*T0); scalar R0 = p0/(rho0*T0);
Rspecific[0] = R0; Rspecific[0] = R0;

View File

@ -1,6 +1,7 @@
EXE_INC = \ EXE_INC = \
-I../engineFoam \ -I../engineFoam \
-I../XiFoam \ -I../XiFoam \
-I../../compressible/rhoPimpleFoam \
-I$(LIB_SRC)/engine/lnInclude \ -I$(LIB_SRC)/engine/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \

View File

@ -81,7 +81,7 @@ int main(int argc, char *argv[])
// --- Pressure corrector loop // --- Pressure corrector loop
while (pimple.correct()) while (pimple.correct())
{ {
#include "hEqn.H" #include "EEqn.H"
#include "pEqn.H" #include "pEqn.H"
} }

View File

@ -5,6 +5,7 @@
psiThermo::New(mesh) psiThermo::New(mesh)
); );
psiThermo& thermo = pThermo(); psiThermo& thermo = pThermo();
thermo.validate(args.executable(), "h", "e");
volScalarField rho volScalarField rho
( (
@ -21,7 +22,6 @@
volScalarField& p = thermo.p(); volScalarField& p = thermo.p();
const volScalarField& psi = thermo.psi(); const volScalarField& psi = thermo.psi();
volScalarField& h = thermo.he();
const volScalarField& T = thermo.T(); const volScalarField& T = thermo.T();

View File

@ -1,12 +0,0 @@
{
solve
(
fvm::ddt(rho, h)
+ fvm::div(phi, h)
- fvm::laplacian(turbulence->alphaEff(), h)
==
- fvc::div(phi, 0.5*magSqr(U))
);
thermo.correct();
}

View File

@ -0,0 +1,26 @@
{
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(phi, volScalarField("Ep", p/rho))
: -dpdt
)
- fvm::laplacian(turbulence->alphaEff(), he)
// - fvm::laplacian(turbulence->muEff(), he) // unit lewis no.
==
reaction->Sh()
);
EEqn.relax();
EEqn.solve();
thermo.correct();
Info<< "min/max(T) = "
<< min(T).value() << ", " << max(T).value() << endl;
}

View File

@ -6,6 +6,7 @@ autoPtr<combustionModels::psiCombustionModel> reaction
); );
psiReactionThermo& thermo = reaction->thermo(); psiReactionThermo& thermo = reaction->thermo();
thermo.validate(args.executable(), "h", "e");
basicMultiComponentMixture& composition = thermo.composition(); basicMultiComponentMixture& composition = thermo.composition();
PtrList<volScalarField>& Y = composition.Y(); PtrList<volScalarField>& Y = composition.Y();
@ -40,7 +41,6 @@ volVectorField U
volScalarField& p = thermo.p(); volScalarField& p = thermo.p();
const volScalarField& psi = thermo.psi(); const volScalarField& psi = thermo.psi();
volScalarField& hs = thermo.he();
const volScalarField& T = thermo.T(); const volScalarField& T = thermo.T();
#include "compressibleCreatePhi.H" #include "compressibleCreatePhi.H"
@ -84,7 +84,7 @@ forAll(Y, i)
{ {
fields.add(Y[i]); fields.add(Y[i]);
} }
fields.add(hs); fields.add(thermo.he());
volScalarField dQ volScalarField dQ
( (

View File

@ -1,21 +0,0 @@
{
fvScalarMatrix hsEqn
(
fvm::ddt(rho, hs)
+ mvConvection->fvmDiv(phi, hs)
- fvm::laplacian(turbulence->alphaEff(), hs)
// - fvm::laplacian(turbulence->muEff(), hs) // unit lewis no.
==
dpdt
- (fvc::ddt(rho, K) + fvc::div(phi, K))
+ reaction->Sh()
);
hsEqn.relax();
hsEqn.solve();
thermo.correct();
Info<< "T gas min/max = " << min(T).value() << ", "
<< max(T).value() << endl;
}

View File

@ -70,7 +70,7 @@ int main(int argc, char *argv[])
{ {
#include "UEqn.H" #include "UEqn.H"
#include "YEqn.H" #include "YEqn.H"
#include "hsEqn.H" #include "EEqn.H"
// --- Pressure corrector loop // --- Pressure corrector loop
while (pimple.correct()) while (pimple.correct())

View File

@ -1,4 +1,5 @@
EXE_INC = \ EXE_INC = \
-I../reactingFoam \
-I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \ -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \

View File

@ -1,16 +0,0 @@
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(phi, U)
+ turbulence->divDevRhoReff(U)
==
rho*g
);
UEqn.relax();
if (pimple.momentumPredictor())
{
solve(UEqn == -fvc::grad(p));
K = 0.5*magSqr(U);
}

View File

@ -1,47 +0,0 @@
tmp<fv::convectionScheme<scalar> > mvConvection
(
fv::convectionScheme<scalar>::New
(
mesh,
fields,
phi,
mesh.divScheme("div(phi,Yi_h)")
)
);
{
reaction->correct();
dQ = reaction->dQ();
label inertIndex = -1;
volScalarField Yt(0.0*Y[0]);
forAll(Y, i)
{
if (Y[i].name() != inertSpecie)
{
volScalarField& Yi = Y[i];
fvScalarMatrix YiEqn
(
fvm::ddt(rho, Yi)
+ mvConvection->fvmDiv(phi, Yi)
- fvm::laplacian(turbulence->muEff(), Yi)
==
reaction->R(Yi)
);
YiEqn.relax();
YiEqn.solve(mesh.solver("Yi"));
Yi.max(0.0);
Yt += Yi;
}
else
{
inertIndex = i;
}
}
Y[inertIndex] = scalar(1) - Yt;
Y[inertIndex].max(0.0);
}

View File

@ -6,6 +6,7 @@ autoPtr<combustionModels::rhoCombustionModel> reaction
); );
rhoReactionThermo& thermo = reaction->thermo(); rhoReactionThermo& thermo = reaction->thermo();
thermo.validate(args.executable(), "h", "e");
basicMultiComponentMixture& composition = thermo.composition(); basicMultiComponentMixture& composition = thermo.composition();
PtrList<volScalarField>& Y = composition.Y(); PtrList<volScalarField>& Y = composition.Y();
@ -40,7 +41,6 @@ volVectorField U
volScalarField& p = thermo.p(); volScalarField& p = thermo.p();
const volScalarField& psi = thermo.psi(); const volScalarField& psi = thermo.psi();
volScalarField& hs = thermo.he();
const volScalarField& T = thermo.T(); const volScalarField& T = thermo.T();
@ -86,7 +86,7 @@ forAll(Y, i)
{ {
fields.add(Y[i]); fields.add(Y[i]);
} }
fields.add(hs); fields.add(thermo.he());
volScalarField dQ volScalarField dQ
( (

View File

@ -1,21 +0,0 @@
{
fvScalarMatrix hsEqn
(
fvm::ddt(rho, hs)
+ mvConvection->fvmDiv(phi, hs)
- fvm::laplacian(turbulence->alphaEff(), hs)
// - fvm::laplacian(turbulence->muEff(), hs) // unit lewis no.
==
dpdt
- (fvc::ddt(rho, K) + fvc::div(phi, K))
+ reaction->Sh()
);
hsEqn.relax();
hsEqn.solve();
thermo.correct();
Info<< "min/max(T) = "
<< min(T).value() << ", " << max(T).value() << endl;
}

View File

@ -72,7 +72,7 @@ int main(int argc, char *argv[])
{ {
#include "UEqn.H" #include "UEqn.H"
#include "YEqn.H" #include "YEqn.H"
#include "hsEqn.H" #include "EEqn.H"
// --- Pressure corrector loop // --- Pressure corrector loop
while (pimple.correct()) while (pimple.correct())

View File

@ -5,9 +5,9 @@
psiThermo::New(mesh) psiThermo::New(mesh)
); );
psiThermo& thermo = pThermo(); psiThermo& thermo = pThermo();
thermo.validate(args.executable(), "h", "e");
volScalarField& p = thermo.p(); volScalarField& p = thermo.p();
volScalarField& h = thermo.he();
const volScalarField& psi = thermo.psi(); const volScalarField& psi = thermo.psi();
volScalarField rho volScalarField rho

View File

@ -1,16 +0,0 @@
{
fvScalarMatrix hEqn
(
fvm::ddt(rho, h)
+ fvm::div(phi, h)
- fvm::laplacian(turbulence->alphaEff(), h)
==
dpdt
- fvc::ddt(rho, K) + fvc::div(phi, K)
);
hEqn.relax();
hEqn.solve();
thermo.correct();
}

View File

@ -75,7 +75,7 @@ int main(int argc, char *argv[])
while (pimple.loop()) while (pimple.loop())
{ {
#include "UEqn.H" #include "UEqn.H"
#include "hEqn.H" #include "EEqn.H"
// --- Pressure corrector loop // --- Pressure corrector loop
while (pimple.correct()) while (pimple.correct())

View File

@ -75,7 +75,7 @@ int main(int argc, char *argv[])
while (pimple.loop()) while (pimple.loop())
{ {
#include "UEqn.H" #include "UEqn.H"
#include "hEqn.H" #include "EEqn.H"
// --- Pressure corrector loop // --- Pressure corrector loop
while (pimple.correct()) while (pimple.correct())

View File

@ -85,7 +85,7 @@ int main(int argc, char *argv[])
turbulence->correct(); turbulence->correct();
#include "UEqn.H" #include "UEqn.H"
#include "hEqn.H" #include "EEqn.H"
// --- Pressure corrector loop // --- Pressure corrector loop
while (pimple.correct()) while (pimple.correct())

View File

@ -78,7 +78,7 @@ int main(int argc, char *argv[])
while (pimple.loop()) while (pimple.loop())
{ {
#include "UEqn.H" #include "UEqn.H"
#include "hEqn.H" #include "EEqn.H"
// --- Pressure corrector loop // --- Pressure corrector loop
while (pimple.correct()) while (pimple.correct())

View File

@ -1,4 +1,5 @@
EXE_INC = \ EXE_INC = \
-I../../compressible/rhoPimpleFoam \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \ -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
-I$(LIB_SRC)/finiteVolume/lnInclude -I$(LIB_SRC)/finiteVolume/lnInclude

View File

@ -0,0 +1,26 @@
{
volScalarField& he = thermo.he();
fvScalarMatrix EEqn
(
fvm::div(phi, he)
+ (
he.name() == "e"
? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
: fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
)
- fvm::laplacian(turb.alphaEff(), he)
==
rad.Sh(thermo)
+ sources(rho, he)
);
EEqn.relax();
EEqn.solve();
thermo.correct();
rad.correct();
Info<< "Min/max T:" << min(thermo.T()).value() << ' '
<< max(thermo.T()).value() << endl;
}

View File

@ -1,22 +0,0 @@
{
fvScalarMatrix hEqn
(
fvm::div(phi, h)
- fvm::laplacian(turb.alphaEff(), h)
==
- fvc::div(phi, 0.5*magSqr(U), "div(phi,K)")
+ rad.Sh(thermo)
+ sources(rho, h)
);
hEqn.relax();
hEqn.solve();
thermo.correct();
rad.correct();
Info<< "Min/max T:" << min(thermo.T()).value() << ' '
<< max(thermo.T()).value() << endl;
}

View File

@ -1,6 +1,8 @@
const fvMesh& mesh = fluidRegions[i]; const fvMesh& mesh = fluidRegions[i];
rhoThermo& thermo = thermoFluid[i]; rhoThermo& thermo = thermoFluid[i];
thermo.validate(args.executable(), "h", "e");
volScalarField& rho = rhoFluid[i]; volScalarField& rho = rhoFluid[i];
volVectorField& U = UFluid[i]; volVectorField& U = UFluid[i];
surfaceScalarField& phi = phiFluid[i]; surfaceScalarField& phi = phiFluid[i];
@ -9,7 +11,6 @@
volScalarField& p = thermo.p(); volScalarField& p = thermo.p();
const volScalarField& psi = thermo.psi(); const volScalarField& psi = thermo.psi();
volScalarField& h = thermo.he();
IObasicSourceList& sources = heatSources[i]; IObasicSourceList& sources = heatSources[i];

View File

@ -4,7 +4,7 @@
rho.storePrevIter(); rho.storePrevIter();
{ {
#include "UEqn.H" #include "UEqn.H"
#include "hEqn.H" #include "EEqn.H"
#include "pEqn.H" #include "pEqn.H"
} }

View File

@ -0,0 +1,27 @@
{
volScalarField& he = thermo.he();
fvScalarMatrix EEqn
(
fvm::ddt(rho, he) + fvm::div(phi, he)
+ fvc::ddt(rho, K) + fvc::div(phi, K)
+ (
he.name() == "e"
? fvc::div(phi, volScalarField("Ep", p/rho))
: -dpdt
)
- fvm::laplacian(turb.alphaEff(), he)
==
rad.Sh(thermo)
+ sources(rho, he)
);
EEqn.relax();
EEqn.solve(mesh.solver(he.select(finalIter)));
thermo.correct();
rad.correct();
Info<< "Min/max T:" << min(thermo.T()).value() << ' '
<< max(thermo.T()).value() << endl;
}

View File

@ -1,24 +0,0 @@
{
fvScalarMatrix hEqn
(
fvm::ddt(rho, h)
+ fvm::div(phi, h)
- fvm::laplacian(turb.alphaEff(), h)
==
dpdt
- (fvc::ddt(rho, K) + fvc::div(phi, K))
+ rad.Sh(thermo)
+ sources(rho, h)
);
hEqn.relax();
hEqn.solve(mesh.solver(h.select(finalIter)));
thermo.correct();
rad.correct();
Info<< "Min/max T:" << min(thermo.T()).value() << ' '
<< max(thermo.T()).value() << endl;
}

View File

@ -1,6 +1,8 @@
fvMesh& mesh = fluidRegions[i]; fvMesh& mesh = fluidRegions[i];
rhoThermo& thermo = thermoFluid[i]; rhoThermo& thermo = thermoFluid[i];
thermo.validate(args.executable(), "h", "e");
volScalarField& rho = rhoFluid[i]; volScalarField& rho = rhoFluid[i];
volVectorField& U = UFluid[i]; volVectorField& U = UFluid[i];
surfaceScalarField& phi = phiFluid[i]; surfaceScalarField& phi = phiFluid[i];
@ -11,7 +13,6 @@
volScalarField& p = thermo.p(); volScalarField& p = thermo.p();
const volScalarField& psi = thermo.psi(); const volScalarField& psi = thermo.psi();
volScalarField& h = thermo.he();
volScalarField& p_rgh = p_rghFluid[i]; volScalarField& p_rgh = p_rghFluid[i];
const volScalarField& gh = ghFluid[i]; const volScalarField& gh = ghFluid[i];

View File

@ -9,8 +9,7 @@ if (oCorr == 0)
} }
#include "UEqn.H" #include "UEqn.H"
#include "EEqn.H"
#include "hEqn.H"
// --- PISO loop // --- PISO loop
for (int corr=0; corr<nCorr; corr++) for (int corr=0; corr<nCorr; corr++)

View File

@ -1,19 +0,0 @@
fvVectorMatrix UEqn
(
//pZones.ddt(rho, U)
fvm::ddt(rho, U)
+ fvm::div(phi, U)
+ turbulence->divDevRhoReff(U)
==
rho.dimensionedInternalField()*g
+ parcels.SU(U)
);
sources.constrain(UEqn);
pZones.addResistance(UEqn);
if (pimple.momentumPredictor())
{
solve(UEqn == -fvc::grad(p) + sources(rho, U));
}

View File

@ -1,52 +0,0 @@
tmp<fv::convectionScheme<scalar> > mvConvection
(
fv::convectionScheme<scalar>::New
(
mesh,
fields,
phi,
mesh.divScheme("div(phi,Yi_h)")
)
);
combustion->correct();
dQ = combustion->dQ();
if (solveSpecies)
{
label inertIndex = -1;
volScalarField Yt(0.0*Y[0]);
forAll(Y, i)
{
if (Y[i].name() != inertSpecie)
{
volScalarField& Yi = Y[i];
fvScalarMatrix YEqn
(
fvm::ddt(rho, Yi)
+ mvConvection->fvmDiv(phi, Yi)
- fvm::laplacian(turbulence->muEff(), Yi)
==
parcels.SYi(i, Yi)
+ combustion->R(Yi)
+ sources(rho, Yi)
);
sources.constrain(YEqn);
YEqn.solve(mesh.solver("Yi"));
Yi.max(0.0);
Yt += Yi;
}
else
{
inertIndex = i;
}
}
Y[inertIndex] = scalar(1) - Yt;
Y[inertIndex].max(0.0);
}

View File

@ -1,28 +0,0 @@
if (chemistry.chemistry())
{
Info<< "Solving chemistry" << endl;
// update reaction rates
chemistry.calculate();
// turbulent time scale
if (turbulentReaction)
{
typedef DimensionedField<scalar, volMesh> dsfType;
const dimensionedScalar e0("e0", sqr(dimLength)/pow3(dimTime), SMALL);
const dsfType tk =
Cmix*sqrt(turbulence->muEff()/rho/(turbulence->epsilon() + e0));
const dsfType tc = chemistry.tc()().dimensionedInternalField();
kappa = tc/(tc + tk);
}
else
{
kappa = 1.0;
}
chemistrySh = kappa*chemistry.Sh()();
}

View File

@ -1,9 +0,0 @@
Info<< "\nConstructing reacting cloud" << endl;
basicReactingMultiphaseCloud parcels
(
"reactingCloud1",
rho,
U,
g,
slgThermo
);

View File

@ -1,121 +0,0 @@
Info<< "Creating combustion model\n" << endl;
autoPtr<combustionModels::rhoCombustionModel> combustion
(
combustionModels::rhoCombustionModel::New(mesh)
);
rhoReactionThermo& thermo = combustion->thermo();
SLGThermo slgThermo(mesh, thermo);
basicMultiComponentMixture& composition = thermo.composition();
PtrList<volScalarField>& Y = composition.Y();
const word inertSpecie(thermo.lookup("inertSpecie"));
if (!composition.contains(inertSpecie))
{
FatalErrorIn(args.executable())
<< "Specified inert specie '" << inertSpecie << "' not found in "
<< "species list. Available species:" << composition.species()
<< exit(FatalError);
}
volScalarField& p = thermo.p();
volScalarField& hs = thermo.he();
const volScalarField& T = thermo.T();
const volScalarField& psi = thermo.psi();
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
thermo.rho()
);
Info<< "\nReading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
#include "compressibleCreatePhi.H"
dimensionedScalar rhoMax
(
mesh.solutionDict().subDict("PIMPLE").lookup("rhoMax")
);
dimensionedScalar rhoMin
(
mesh.solutionDict().subDict("PIMPLE").lookup("rhoMin")
);
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());
Info<< "Creating multi-variate interpolation scheme\n" << endl;
multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
forAll(Y, i)
{
fields.add(Y[i]);
}
fields.add(hs);
volScalarField dQ
(
IOobject
(
"dQ",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
);
volScalarField rDeltaT
(
IOobject
(
"rDeltaT",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("one", dimless/dimTime, 1),
zeroGradientFvPatchScalarField::typeName
);

View File

@ -1,25 +0,0 @@
{
fvScalarMatrix hsEqn
(
fvm::ddt(rho, hs)
+ mvConvection->fvmDiv(phi, hs)
- fvm::laplacian(turbulence->alphaEff(), hs)
==
- fvc::div(phi, 0.5*magSqr(U), "div(phi,K)")
+ parcels.Sh(hs)
+ radiation->Sh(thermo)
+ combustion->Sh()
+ sources(rho, hs)
);
sources.constrain(hsEqn);
hsEqn.solve();
thermo.correct();
radiation->correct();
Info<< "T gas min/max = " << min(T).value() << ", "
<< max(T).value() << endl;
}

View File

@ -1,66 +0,0 @@
{
rho = thermo.rho();
// Thermodynamic density needs to be updated by psi*d(p) after the
// pressure solution - done in 2 parts. Part 1:
thermo.rho() -= psi*p;
volScalarField rAU(1.0/UEqn.A());
volVectorField HbyA("HbyA", U);
HbyA = rAU*(UEqn == sources(rho, U))().H();
surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf());
if (pZones.size() == 0)
{
// ddtPhiCorr only used without porosity
phiHbyA += fvc::ddtPhiCorr(rAU, rho, U, phi);
}
phiHbyA *= fvc::interpolate(rho);
fvScalarMatrix pDDtEqn
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+ fvc::div(phiHbyA)
==
parcels.Srho()
+ sources(psi, p, rho.name())
);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
pDDtEqn
- fvm::laplacian(rho*rAU, p)
);
sources.constrain(pDDtEqn, rho.name());
pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA + pEqn.flux();
}
}
p.relax();
// Second part of thermodynamic density update
thermo.rho() += psi*p;
#include "rhoEqn.H" // NOTE: flux and time scales now inconsistent
#include "compressibleContinuityErrs.H"
U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();
sources.correct(U);
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
}

View File

@ -1,9 +0,0 @@
dictionary additional = mesh.solutionDict().subDict("additional");
// pressure work term for enthalpy equation
bool pressureWork = additional.lookupOrDefault("pressureWork", true);
bool pressureWorkTimeDerivative =
additional.lookupOrDefault("pressureWorkTimeDerivative", true);
// flag to activate solve transport for each specie (Y vector)
bool solveSpecies = additional.lookupOrDefault("solveSpecies", true);

View File

@ -1,50 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 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)
==
parcels.Srho(rho)
+ sources(rho)
);
sources.constrain(rhoEqn);
rhoEqn.solve();
Info<< "rho min/max = " << min(rho).value() << ", " << max(rho).value()
<< endl;
}
// ************************************************************************* //

View File

@ -0,0 +1,31 @@
{
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(phi, volScalarField("Ep", p/rho))
: -dpdt
)
- fvm::laplacian(turbulence->alphaEff(), he)
==
combustion->Sh()
+ coalParcels.Sh(he)
+ limestoneParcels.Sh(he)
+ radiation->Sh(thermo)
+ sources(rho, he)
);
EEqn.relax();
sources.constrain(EEqn);
EEqn.solve();
thermo.correct();
radiation->correct();
Info<< "T gas min/max = " << min(T).value() << ", "
<< max(T).value() << endl;
}

View File

@ -92,7 +92,7 @@ int main(int argc, char *argv[])
{ {
#include "UEqn.H" #include "UEqn.H"
#include "YEqn.H" #include "YEqn.H"
#include "hsEqn.H" #include "EEqn.H"
// --- Pressure corrector loop // --- Pressure corrector loop
while (pimple.correct()) while (pimple.correct())

View File

@ -6,6 +6,7 @@
); );
psiReactionThermo& thermo = combustion->thermo(); psiReactionThermo& thermo = combustion->thermo();
thermo.validate(args.executable(), "h", "e");
SLGThermo slgThermo(mesh, thermo); SLGThermo slgThermo(mesh, thermo);
@ -23,7 +24,6 @@
} }
volScalarField& p = thermo.p(); volScalarField& p = thermo.p();
volScalarField& hs = thermo.he();
const volScalarField& T = thermo.T(); const volScalarField& T = thermo.T();
const volScalarField& psi = thermo.psi(); const volScalarField& psi = thermo.psi();
@ -33,7 +33,7 @@
{ {
fields.add(Y[i]); fields.add(Y[i]);
} }
fields.add(hs); fields.add(thermo.he());
volScalarField rho volScalarField rho
( (

View File

@ -1,29 +0,0 @@
{
fvScalarMatrix hsEqn
(
fvm::ddt(rho, hs)
+ mvConvection->fvmDiv(phi, hs)
- fvm::laplacian(turbulence->alphaEff(), hs)
==
dpdt
- (fvc::ddt(rho, K) + fvc::div(phi, K))
+ combustion->Sh()
+ coalParcels.Sh(hs)
+ limestoneParcels.Sh(hs)
+ radiation->Sh(thermo)
+ sources(rho, hs)
);
hsEqn.relax();
sources.constrain(hsEqn);
hsEqn.solve();
thermo.correct();
radiation->correct();
Info<< "T gas min/max = " << min(T).value() << ", "
<< max(T).value() << endl;
}

View File

@ -0,0 +1,8 @@
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
set -x
wmake
wmake icoUncoupledKinematicParcelDyMFoam
# ----------------------------------------------------------------- end-of-file

View File

@ -1,5 +1,5 @@
EXE_INC = \ EXE_INC = \
-I../icoUncoupledKinematicParcelFoam \ -I.. \
-I$(LIB_SRC)/lagrangian/basic/lnInclude \ -I$(LIB_SRC)/lagrangian/basic/lnInclude \
-I$(LIB_SRC)/lagrangian/intermediate/lnInclude \ -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \

View File

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

View File

@ -1,54 +0,0 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I${LIB_SRC}/meshTools/lnInclude \
-I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
-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)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/properties/liquidMixtureProperties/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/properties/solidMixtureProperties/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/radiationModels/lnInclude \
-I$(LIB_SRC)/ODE/lnInclude \
-I$(LIB_SRC)/regionModels/regionModel/lnInclude \
-I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
-I$(LIB_SRC)/combustionModels/lnInclude \
-I$(LIB_SRC)/fieldSources/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
-I$(FOAM_SOLVERS)/combustion/reactingFoam
EXE_LIBS = \
-lfiniteVolume \
-lmeshTools \
-lcompressibleTurbulenceModel \
-lcompressibleRASModels \
-lcompressibleLESModels \
-llagrangian \
-llagrangianIntermediate \
-lspecie \
-lfluidThermophysicalModels \
-lliquidProperties \
-lliquidMixtureProperties \
-lsolidProperties \
-lsolidMixtureProperties \
-lthermophysicalFunctions \
-lreactionThermophysicalModels \
-lSLGThermo \
-lchemistryModel \
-lradiationModels \
-lODE \
-lregionModels \
-lsurfaceFilmModels \
-lcombustionModels \
-lfieldSources \
-lsampling

View File

@ -1,9 +0,0 @@
Info<< "\nConstructing reacting cloud" << endl;
basicReactingMultiphaseCloud parcels
(
"reactingCloud1",
rho,
U,
g,
slgThermo
);

View File

@ -1,2 +0,0 @@
Info<< "Creating sources\n" << endl;
IObasicSourceList sources(mesh);

View File

@ -1,3 +0,0 @@
Info<< "Creating porous zones" << nl << endl;
porousZones pZones(mesh);

View File

@ -1,32 +0,0 @@
{
fvScalarMatrix hsEqn
(
fvm::ddt(rho, hs)
+ mvConvection->fvmDiv(phi, hs)
- fvm::laplacian(turbulence->alphaEff(), hs)
==
- (fvc::ddt(rho, K) + fvc::div(phi, K))
+ parcels.Sh(hs)
+ radiation->Sh(thermo)
+ combustion->Sh()
+ sources(rho, hs)
);
if (pressureWorkTimeDerivative)
{
hsEqn -= dpdt;
}
hsEqn.relax();
sources.constrain(hsEqn);
hsEqn.solve();
thermo.correct();
radiation->correct();
Info<< "T gas min/max = " << min(T).value() << ", "
<< max(T).value() << endl;
}

View File

@ -1,66 +0,0 @@
{
rho = thermo.rho();
// Thermodynamic density needs to be updated by psi*d(p) after the
// pressure solution - done in 2 parts. Part 1:
thermo.rho() -= psi*p;
volScalarField rAU(1.0/UEqn.A());
volVectorField HbyA("HbyA", U);
HbyA = rAU*(UEqn == sources(rho, U))().H();
surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf());
if (pZones.size() == 0)
{
// ddtPhiCorr only used without porosity
phiHbyA += fvc::ddtPhiCorr(rAU, rho, U, phi);
}
phiHbyA *= fvc::interpolate(rho);
fvScalarMatrix pDDtEqn
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+ fvc::div(phiHbyA)
==
parcels.Srho()
+ sources(psi, p, rho.name())
);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
pDDtEqn
- fvm::laplacian(rho*rAU, p)
);
sources.constrain(pDDtEqn, rho.name());
pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA + pEqn.flux();
}
}
// Second part of thermodynamic density update
thermo.rho() += psi*p;
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();
sources.correct(U);
K = 0.5*magSqr(U);
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);
}
}

View File

@ -1,126 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 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
porousExplicitSourceReactingParcelFoam
Description
Transient PIMPLE solver for compressible, laminar or turbulent flow with
reacting multiphase Lagrangian parcels for porous media, including explicit
sources for mass, momentum and energy
The solver includes:
- reacting multiphase parcel cloud
- porous media
- mass, momentum and energy sources
Note: ddtPhiCorr not used here when porous zones are active
- not well defined for porous calculations
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "turbulenceModel.H"
#include "basicReactingMultiphaseCloud.H"
#include "rhoCombustionModel.H"
#include "radiationModel.H"
#include "porousZones.H"
#include "IObasicSourceList.H"
#include "SLGThermo.H"
#include "pimpleControl.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "readGravitationalAcceleration.H"
#include "createFields.H"
#include "createRadiationModel.H"
#include "createClouds.H"
#include "createExplicitSources.H"
#include "createPorousZones.H"
#include "initContinuityErrs.H"
#include "readTimeControls.H"
#include "compressibleCourantNo.H"
#include "setInitialDeltaT.H"
pimpleControl pimple(mesh);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
#include "readTimeControls.H"
#include "readAdditionalSolutionControls.H"
#include "compressibleCourantNo.H"
#include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
parcels.evolve();
#include "rhoEqn.H"
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
#include "UEqn.H"
#include "YEqn.H"
#include "hsEqn.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\n" << endl;
return(0);
}
// ************************************************************************* //

View File

@ -0,0 +1,29 @@
{
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(phi, volScalarField("Ep", p/rho))
: -dpdt
)
- fvm::laplacian(turbulence->alphaEff(), he)
==
parcels.Sh(he)
+ surfaceFilm.Sh()
+ radiation->Sh(thermo)
+ combustion->Sh()
);
EEqn.relax();
EEqn.solve();
thermo.correct();
radiation->correct();
Info<< "T gas min/max = " << min(T).value() << ", "
<< max(T).value() << endl;
}

View File

@ -6,6 +6,7 @@
); );
psiReactionThermo& thermo = combustion->thermo(); psiReactionThermo& thermo = combustion->thermo();
thermo.validate(args.executable(), "h", "e");
SLGThermo slgThermo(mesh, thermo); SLGThermo slgThermo(mesh, thermo);
@ -29,7 +30,6 @@
); );
volScalarField& p = thermo.p(); volScalarField& p = thermo.p();
volScalarField& hs = thermo.he();
const volScalarField& T = thermo.T(); const volScalarField& T = thermo.T();
const volScalarField& psi = thermo.psi(); const volScalarField& psi = thermo.psi();
@ -107,7 +107,7 @@
{ {
fields.add(Y[i]); fields.add(Y[i]);
} }
fields.add(hs); fields.add(thermo.he());
IOdictionary additionalControlsDict IOdictionary additionalControlsDict
( (

View File

@ -1,25 +0,0 @@
{
fvScalarMatrix hsEqn
(
fvm::ddt(rho, hs)
+ mvConvection->fvmDiv(phi, hs)
- fvm::laplacian(turbulence->alphaEff(), hs)
==
dpdt
- (fvc::ddt(rho, K) + fvc::div(phi, K))
+ parcels.Sh(hs)
+ surfaceFilm.Sh()
+ radiation->Sh(thermo)
+ combustion->Sh()
);
hsEqn.relax();
hsEqn.solve();
thermo.correct();
radiation->correct();
Info<< "min/max(T) = " << min(T).value() << ", " << max(T).value() << endl;
}

View File

@ -87,7 +87,7 @@ int main(int argc, char *argv[])
{ {
#include "UEqn.H" #include "UEqn.H"
#include "YEqn.H" #include "YEqn.H"
#include "hsEqn.H" #include "EEqn.H"
// --- Pressure corrector loop // --- Pressure corrector loop
while (pimple.correct()) while (pimple.correct())

View File

@ -0,0 +1,8 @@
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
set -x
wmake
wmake LTSReactingParcelFoam
# ----------------------------------------------------------------- end-of-file

View File

@ -0,0 +1,30 @@
{
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(phi, volScalarField("Ep", p/rho))
: -dpdt
)
- fvm::laplacian(turbulence->alphaEff(), he)
==
parcels.Sh(he)
+ radiation->Sh(thermo)
+ combustion->Sh()
+ sources(rho, he)
);
EEqn.relax();
sources.constrain(EEqn);
EEqn.solve();
thermo.correct();
radiation->correct();
Info<< "T gas min/max = " << min(T).value() << ", "
<< max(T).value() << endl;
}

View File

@ -61,6 +61,7 @@ int main(int argc, char *argv[])
#include "readTimeControls.H" #include "readTimeControls.H"
#include "readAdditionalSolutionControls.H" #include "readAdditionalSolutionControls.H"
#include "createFields.H" #include "createFields.H"
#include "createRDeltaT.H"
#include "createRadiationModel.H" #include "createRadiationModel.H"
#include "createClouds.H" #include "createClouds.H"
#include "createExplicitSources.H" #include "createExplicitSources.H"
@ -93,7 +94,7 @@ int main(int argc, char *argv[])
#include "UEqn.H" #include "UEqn.H"
#include "YEqn.H" #include "YEqn.H"
#include "hsEqn.H" #include "EEqn.H"
// --- Pressure corrector loop // --- Pressure corrector loop
while (pimple.correct()) while (pimple.correct())

View File

@ -1,4 +1,5 @@
EXE_INC = \ EXE_INC = \
-I.. \
-I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/finiteVolume/lnInclude \
-I${LIB_SRC}/meshTools/lnInclude \ -I${LIB_SRC}/meshTools/lnInclude \
-I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \ -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \

View File

@ -0,0 +1,14 @@
volScalarField rDeltaT
(
IOobject
(
"rDeltaT",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("one", dimless/dimTime, 1),
zeroGradientFvPatchScalarField::typeName
);

View File

@ -1,10 +1,10 @@
EXE_INC = \ EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/finiteVolume/lnInclude \
-I${LIB_SRC}/meshTools/lnInclude \ -I${LIB_SRC}/meshTools/lnInclude \
-I${LIB_SRC}/sampling/lnInclude \
-I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \ -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
-I$(LIB_SRC)/lagrangian/basic/lnInclude \ -I$(LIB_SRC)/lagrangian/basic/lnInclude \
-I$(LIB_SRC)/lagrangian/intermediate/lnInclude \ -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
-I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \
-I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \ -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
@ -22,13 +22,13 @@ EXE_INC = \
-I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \ -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
-I$(LIB_SRC)/combustionModels/lnInclude \ -I$(LIB_SRC)/combustionModels/lnInclude \
-I$(LIB_SRC)/fieldSources/lnInclude \ -I$(LIB_SRC)/fieldSources/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
-I$(FOAM_SOLVERS)/combustion/reactingFoam -I$(FOAM_SOLVERS)/combustion/reactingFoam
EXE_LIBS = \ EXE_LIBS = \
-lfiniteVolume \ -lfiniteVolume \
-lmeshTools \ -lmeshTools \
-lsampling \
-lcompressibleTurbulenceModel \ -lcompressibleTurbulenceModel \
-lcompressibleRASModels \ -lcompressibleRASModels \
-lcompressibleLESModels \ -lcompressibleLESModels \
@ -48,5 +48,6 @@ EXE_LIBS = \
-lODE \ -lODE \
-lregionModels \ -lregionModels \
-lsurfaceFilmModels \ -lsurfaceFilmModels \
-lcombustionModels \
-lfieldSources \ -lfieldSources \
-lcombustionModels -lsampling

View File

@ -1,5 +1,6 @@
fvVectorMatrix UEqn fvVectorMatrix UEqn
( (
//pZones.ddt(rho, U)
fvm::ddt(rho, U) fvm::ddt(rho, U)
+ fvm::div(phi, U) + fvm::div(phi, U)
+ turbulence->divDevRhoReff(U) + turbulence->divDevRhoReff(U)
@ -12,6 +13,8 @@
sources.constrain(UEqn); sources.constrain(UEqn);
pZones.addResistance(UEqn);
if (pimple.momentumPredictor()) if (pimple.momentumPredictor())
{ {
solve(UEqn == -fvc::grad(p) + sources(rho, U)); solve(UEqn == -fvc::grad(p) + sources(rho, U));

View File

@ -5,14 +5,15 @@ tmp<fv::convectionScheme<scalar> > mvConvection
mesh, mesh,
fields, fields,
phi, phi,
mesh.divScheme("div(phi,Yi_hs)") mesh.divScheme("div(phi,Yi_h)")
) )
); );
combustion->correct();
dQ = combustion->dQ();
if (solveSpecies)
{ {
combustion->correct();
dQ = combustion->dQ();
label inertIndex = -1; label inertIndex = -1;
volScalarField Yt(0.0*Y[0]); volScalarField Yt(0.0*Y[0]);
@ -22,7 +23,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
{ {
volScalarField& Yi = Y[i]; volScalarField& Yi = Y[i];
fvScalarMatrix YiEqn fvScalarMatrix YEqn
( (
fvm::ddt(rho, Yi) fvm::ddt(rho, Yi)
+ mvConvection->fvmDiv(phi, Yi) + mvConvection->fvmDiv(phi, Yi)
@ -33,11 +34,11 @@ tmp<fv::convectionScheme<scalar> > mvConvection
+ sources(rho, Yi) + sources(rho, Yi)
); );
YiEqn.relax(); YEqn.relax();
sources.constrain(YiEqn); sources.constrain(YEqn);
YiEqn.solve(mesh.solver("Yi")); YEqn.solve(mesh.solver("Yi"));
Yi.max(0.0); Yi.max(0.0);
Yt += Yi; Yt += Yi;

View File

@ -1,5 +1,5 @@
Info<< "\nConstructing reacting cloud" << endl; Info<< "\nConstructing reacting cloud" << endl;
basicReactingCloud parcels basicReactingMultiphaseCloud parcels
( (
"reactingCloud1", "reactingCloud1",
rho, rho,

View File

@ -1,11 +1,12 @@
Info<< "Creating combustion model\n" << endl; Info<< "Creating combustion model\n" << endl;
autoPtr<combustionModels::psiCombustionModel> combustion autoPtr<combustionModels::rhoCombustionModel> combustion
( (
combustionModels::psiCombustionModel::New(mesh) combustionModels::rhoCombustionModel::New(mesh)
); );
psiReactionThermo& thermo = combustion->thermo(); rhoReactionThermo& thermo = combustion->thermo();
thermo.validate(args.executable(), "h", "e");
SLGThermo slgThermo(mesh, thermo); SLGThermo slgThermo(mesh, thermo);
@ -23,7 +24,6 @@
} }
volScalarField& p = thermo.p(); volScalarField& p = thermo.p();
volScalarField& hs = thermo.he();
const volScalarField& T = thermo.T(); const volScalarField& T = thermo.T();
const volScalarField& psi = thermo.psi(); const volScalarField& psi = thermo.psi();
@ -56,6 +56,22 @@
#include "compressibleCreatePhi.H" #include "compressibleCreatePhi.H"
dimensionedScalar rhoMax
(
"rhoMax",
dimDensity,
mesh.solutionDict().subDict("PIMPLE")
.lookupOrDefault<scalar>("rhoMax", GREAT)
);
dimensionedScalar rhoMin
(
"rhoMin",
dimDensity,
mesh.solutionDict().subDict("PIMPLE")
.lookupOrDefault<scalar>("rhoMin", -GREAT)
);
Info<< "Creating turbulence model\n" << endl; Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::turbulenceModel> turbulence autoPtr<compressible::turbulenceModel> turbulence
( (
@ -86,13 +102,15 @@
Info<< "Creating field kinetic energy K\n" << endl; Info<< "Creating field kinetic energy K\n" << endl;
volScalarField K("K", 0.5*magSqr(U)); volScalarField K("K", 0.5*magSqr(U));
Info<< "Creating multi-variate interpolation scheme\n" << endl;
multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields; multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
forAll(Y, i) forAll(Y, i)
{ {
fields.add(Y[i]); fields.add(Y[i]);
} }
fields.add(hs); fields.add(thermo.he());
volScalarField dQ volScalarField dQ
( (
@ -107,6 +125,3 @@
mesh, mesh,
dimensionedScalar("dQ", dimEnergy/dimTime, 0.0) dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
); );
Info<< "Creating sources\n" << endl;
IObasicSourceList sources(mesh);

View File

@ -1,65 +1,43 @@
rho = thermo.rho();
volScalarField rAU(1.0/UEqn.A());
volVectorField HbyA("HbyA", U);
HbyA = rAU*(UEqn == sources(rho, U))().H();
if (pimple.transonic())
{ {
surfaceScalarField phid rho = thermo.rho();
(
"phid",
fvc::interpolate(psi)
*(
(fvc::interpolate(HbyA) & mesh.Sf())
+ fvc::ddtPhiCorr(rAU, rho, U, phi)
)
);
while (pimple.correctNonOrthogonal()) // Thermodynamic density needs to be updated by psi*d(p) after the
// pressure solution - done in 2 parts. Part 1:
thermo.rho() -= psi*p;
volScalarField rAU(1.0/UEqn.A());
volVectorField HbyA("HbyA", U);
HbyA = rAU*(UEqn == sources(rho, U))().H();
surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf());
if (pZones.size() == 0)
{ {
fvScalarMatrix pEqn // ddtPhiCorr only used without porosity
( phiHbyA += fvc::ddtPhiCorr(rAU, rho, U, phi);
fvm::ddt(psi, p)
+ fvm::div(phid, p)
- fvm::laplacian(rho*rAU, p)
==
parcels.Srho()
+ sources(psi, p, rho.name())
);
pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
if (pimple.finalNonOrthogonalIter())
{
phi == pEqn.flux();
} }
}
}
else
{
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::interpolate(rho)
*(
(fvc::interpolate(HbyA) & mesh.Sf())
+ fvc::ddtPhiCorr(rAU, rho, U, phi)
)
);
while (pimple.correctNonOrthogonal()) phiHbyA *= fvc::interpolate(rho);
{
fvScalarMatrix pEqn
fvScalarMatrix pDDtEqn
( (
fvm::ddt(psi, p) fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+ fvc::div(phiHbyA) + fvc::div(phiHbyA)
- fvm::laplacian(rho*rAU, p)
== ==
parcels.Srho() parcels.Srho()
+ sources(psi, p, rho.name()) + sources(psi, p, rho.name())
); );
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
pDDtEqn
- fvm::laplacian(rho*rAU, p)
);
sources.constrain(pDDtEqn, rho.name());
pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
if (pimple.finalNonOrthogonalIter()) if (pimple.finalNonOrthogonalIter())
@ -67,17 +45,28 @@ else
phi = phiHbyA + pEqn.flux(); phi = phiHbyA + pEqn.flux();
} }
} }
}
#include "rhoEqn.H" p.relax();
#include "compressibleContinuityErrs.H"
U = HbyA - rAU*fvc::grad(p); // Second part of thermodynamic density update
U.correctBoundaryConditions(); thermo.rho() += psi*p;
sources.correct(U);
K = 0.5*magSqr(U);
if (thermo.dpdt()) #include "rhoEqn.H" // NOTE: flux and time scales now inconsistent
{ #include "compressibleContinuityErrs.H"
U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();
sources.correct(U);
K = 0.5*magSqr(U);
if (thermo.dpdt())
{
dpdt = fvc::ddt(p); dpdt = fvc::ddt(p);
}
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
} }

View File

@ -22,22 +22,32 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application Application
reactingParcelFoam porousExplicitSourceReactingParcelFoam
Description Description
Transient PIMPLE solver for compressible, laminar or turbulent flow with Transient PIMPLE solver for compressible, laminar or turbulent flow with
reacting Lagrangian parcels. reacting multiphase Lagrangian parcels for porous media, including explicit
sources for mass, momentum and energy
The solver includes:
- reacting multiphase parcel cloud
- porous media
- mass, momentum and energy sources
Note: ddtPhiCorr not used here when porous zones are active
- not well defined for porous calculations
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#include "fvCFD.H" #include "fvCFD.H"
#include "turbulenceModel.H" #include "turbulenceModel.H"
#include "basicReactingCloud.H" #include "basicReactingMultiphaseCloud.H"
#include "psiCombustionModel.H" #include "rhoCombustionModel.H"
#include "radiationModel.H" #include "radiationModel.H"
#include "porousZones.H"
#include "IObasicSourceList.H"
#include "SLGThermo.H" #include "SLGThermo.H"
#include "pimpleControl.H" #include "pimpleControl.H"
#include "IObasicSourceList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -48,16 +58,19 @@ int main(int argc, char *argv[])
#include "createTime.H" #include "createTime.H"
#include "createMesh.H" #include "createMesh.H"
#include "readGravitationalAcceleration.H" #include "readGravitationalAcceleration.H"
pimpleControl pimple(mesh);
#include "createFields.H" #include "createFields.H"
#include "createClouds.H"
#include "createRadiationModel.H" #include "createRadiationModel.H"
#include "createClouds.H"
#include "createExplicitSources.H"
#include "createPorousZones.H"
#include "initContinuityErrs.H" #include "initContinuityErrs.H"
#include "readTimeControls.H" #include "readTimeControls.H"
#include "compressibleCourantNo.H" #include "compressibleCourantNo.H"
#include "setInitialDeltaT.H" #include "setInitialDeltaT.H"
pimpleControl pimple(mesh);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl; Info<< "\nStarting time loop\n" << endl;
@ -65,6 +78,7 @@ int main(int argc, char *argv[])
while (runTime.run()) while (runTime.run())
{ {
#include "readTimeControls.H" #include "readTimeControls.H"
#include "readAdditionalSolutionControls.H"
#include "compressibleCourantNo.H" #include "compressibleCourantNo.H"
#include "setDeltaT.H" #include "setDeltaT.H"
@ -81,7 +95,7 @@ int main(int argc, char *argv[])
{ {
#include "UEqn.H" #include "UEqn.H"
#include "YEqn.H" #include "YEqn.H"
#include "hsEqn.H" #include "EEqn.H"
// --- Pressure corrector loop // --- Pressure corrector loop
while (pimple.correct()) while (pimple.correct())

View File

@ -1,8 +1,4 @@
dictionary additional = mesh.solutionDict().subDict("additional"); dictionary additional = mesh.solutionDict().subDict("additional");
// pressure work term for enthalpy equation
bool pressureWorkTimeDerivative =
additional.lookupOrDefault("pressureWorkTimeDerivative", false);
// flag to activate solve transport for each specie (Y vector) // flag to activate solve transport for each specie (Y vector)
bool solveSpecies = additional.lookupOrDefault("solveSpecies", true); bool solveSpecies = additional.lookupOrDefault("solveSpecies", true);

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -42,6 +42,9 @@ Description
sources.constrain(rhoEqn); sources.constrain(rhoEqn);
rhoEqn.solve(); rhoEqn.solve();
Info<< "rho min/max = " << min(rho).value() << ", " << max(rho).value()
<< endl;
} }
// ************************************************************************* // // ************************************************************************* //

View File

@ -1,5 +1,4 @@
EXE_INC = \ EXE_INC = \
-I$(FOAM_SOLVERS)/lagrangian/reactingParcelFoam \
-I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/finiteVolume/lnInclude \
-I${LIB_SRC}/meshTools/lnInclude \ -I${LIB_SRC}/meshTools/lnInclude \
-I${LIB_SRC}/sampling/lnInclude \ -I${LIB_SRC}/sampling/lnInclude \

View File

@ -1,6 +1,5 @@
fvVectorMatrix UEqn fvVectorMatrix UEqn
( (
//pZones.ddt(rho, U)
fvm::ddt(rho, U) fvm::ddt(rho, U)
+ fvm::div(phi, U) + fvm::div(phi, U)
+ turbulence->divDevRhoReff(U) + turbulence->divDevRhoReff(U)
@ -13,11 +12,8 @@
sources.constrain(UEqn); sources.constrain(UEqn);
pZones.addResistance(UEqn);
if (pimple.momentumPredictor()) if (pimple.momentumPredictor())
{ {
solve(UEqn == -fvc::grad(p) + sources(rho, U)); solve(UEqn == -fvc::grad(p) + sources(rho, U));
K = 0.5*magSqr(U); K = 0.5*magSqr(U);
} }

View File

@ -1,4 +1,3 @@
tmp<fv::convectionScheme<scalar> > mvConvection tmp<fv::convectionScheme<scalar> > mvConvection
( (
fv::convectionScheme<scalar>::New fv::convectionScheme<scalar>::New
@ -6,15 +5,14 @@ tmp<fv::convectionScheme<scalar> > mvConvection
mesh, mesh,
fields, fields,
phi, phi,
mesh.divScheme("div(phi,Yi_h)") mesh.divScheme("div(phi,Yi_hs)")
) )
); );
combustion->correct();
dQ = combustion->dQ();
if (solveSpecies)
{ {
combustion->correct();
dQ = combustion->dQ();
label inertIndex = -1; label inertIndex = -1;
volScalarField Yt(0.0*Y[0]); volScalarField Yt(0.0*Y[0]);

View File

@ -1,11 +1,11 @@
Info<< "Creating combustion model\n" << endl; Info<< "Creating combustion model\n" << endl;
autoPtr<combustionModels::rhoCombustionModel> combustion autoPtr<combustionModels::psiCombustionModel> combustion
( (
combustionModels::rhoCombustionModel::New(mesh) combustionModels::psiCombustionModel::New(mesh)
); );
rhoReactionThermo& thermo = combustion->thermo(); psiReactionThermo& thermo = combustion->thermo();
SLGThermo slgThermo(mesh, thermo); SLGThermo slgThermo(mesh, thermo);
@ -71,7 +71,21 @@
// Set the turbulence into the combustion model // Set the turbulence into the combustion model
combustion->setTurbulence(turbulence()); combustion->setTurbulence(turbulence());
Info<< "Creating multi-variate interpolation scheme\n" << endl; 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));
multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields; multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
forAll(Y, i) forAll(Y, i)
@ -94,18 +108,5 @@
dimensionedScalar("dQ", dimEnergy/dimTime, 0.0) dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
); );
Info<< "Creating field dpdt\n" << endl; Info<< "Creating sources\n" << endl;
volScalarField dpdt IObasicSourceList sources(mesh);
(
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));

View File

@ -0,0 +1,83 @@
rho = thermo.rho();
volScalarField rAU(1.0/UEqn.A());
volVectorField HbyA("HbyA", U);
HbyA = rAU*(UEqn == sources(rho, U))().H();
if (pimple.transonic())
{
surfaceScalarField phid
(
"phid",
fvc::interpolate(psi)
*(
(fvc::interpolate(HbyA) & mesh.Sf())
+ fvc::ddtPhiCorr(rAU, rho, U, phi)
)
);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvm::ddt(psi, p)
+ fvm::div(phid, p)
- fvm::laplacian(rho*rAU, p)
==
parcels.Srho()
+ sources(psi, p, rho.name())
);
pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
if (pimple.finalNonOrthogonalIter())
{
phi == pEqn.flux();
}
}
}
else
{
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::interpolate(rho)
*(
(fvc::interpolate(HbyA) & mesh.Sf())
+ fvc::ddtPhiCorr(rAU, rho, U, phi)
)
);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvm::ddt(psi, p)
+ fvc::div(phiHbyA)
- fvm::laplacian(rho*rAU, p)
==
parcels.Srho()
+ sources(psi, p, rho.name())
);
pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA + pEqn.flux();
}
}
}
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();
sources.correct(U);
K = 0.5*magSqr(U);
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);
}

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License

View File

@ -1,6 +1,5 @@
EXE_INC = \ EXE_INC = \
-I$(FOAM_SOLVERS)/lagrangian/sprayFoam \ -I.. \
-I$(FOAM_SOLVERS)/lagrangian/reactingParcelFoam \
-I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/finiteVolume/lnInclude \
-I${LIB_SRC}/meshTools/lnInclude \ -I${LIB_SRC}/meshTools/lnInclude \
-I${LIB_SRC}/sampling/lnInclude \ -I${LIB_SRC}/sampling/lnInclude \

View File

@ -38,6 +38,7 @@ FoamFile
frontAndBack frontAndBack
{ {
type empty; type empty;
inGroups 1(empty);
nFaces 8000; nFaces 8000;
startFace 8140; startFace 8140;
} }

View File

@ -0,0 +1,58 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class polyBoundaryMesh;
location "constant/polyMesh";
object boundary;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
6
(
front
{
type wall;
nFaces 700;
startFace 63400;
}
back
{
type wall;
nFaces 700;
startFace 64100;
}
wall
{
type wall;
nFaces 1400;
startFace 64800;
}
porosityWall
{
type wall;
nFaces 1600;
startFace 66200;
}
inlet
{
type patch;
nFaces 400;
startFace 67800;
}
outlet
{
type patch;
nFaces 400;
startFace 68200;
}
)
// ************************************************************************* //

View File

@ -32,6 +32,7 @@ FoamFile
frontAndBack frontAndBack
{ {
type empty; type empty;
inGroups 1(empty);
nFaces 800; nFaces 800;
startFace 840; startFace 840;
} }

View File

@ -38,6 +38,7 @@ FoamFile
frontAndBack frontAndBack
{ {
type empty; type empty;
inGroups 1(empty);
nFaces 8000; nFaces 8000;
startFace 8140; startFace 8140;
} }

View File

@ -20,12 +20,14 @@ FoamFile
back back
{ {
type symmetryPlane; type symmetryPlane;
inGroups 1(symmetryPlane);
nFaces 9340; nFaces 9340;
startFace 265900; startFace 265900;
} }
front front
{ {
type symmetryPlane; type symmetryPlane;
inGroups 1(symmetryPlane);
nFaces 9340; nFaces 9340;
startFace 275240; startFace 275240;
} }

View File

@ -15,7 +15,7 @@ FoamFile
} }
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
psiChemistryModel ODEChemistryModel<gasThermoPhysics>; rhoChemistryModel ODEChemistryModel<gasThermoPhysics>;
chemistry off; chemistry off;

View File

@ -15,7 +15,7 @@ FoamFile
} }
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
combustionModel PaSR<psiChemistryCombustion>; combustionModel PaSR<rhoChemistryCombustion>;
active false; active false;

View File

@ -26,6 +26,7 @@ FoamFile
frontAndBack frontAndBack
{ {
type empty; type empty;
inGroups 1(empty);
nFaces 800; nFaces 800;
startFace 840; startFace 840;
} }

Some files were not shown because too many files have changed in this diff Show More