GIT: conflict resolution

This commit is contained in:
andy 2012-09-26 12:13:44 +01:00
commit f1deb4445e
636 changed files with 20768 additions and 5580 deletions

View File

@ -0,0 +1,17 @@
{
volScalarField& hea = thermo.he();
solve
(
betav*fvm::ddt(rho, hea) + mvConvection->fvmDiv(phi, hea)
+ betav*fvc::ddt(rho, K) + fvc::div(phi, K)
+ (
hea.name() == "ea"
? fvc::div(phi, volScalarField("Ep", p/rho))
: -betav*dpdt
)
- fvm::laplacian(Db, hea)
);
thermo.correct();
}

View File

@ -0,0 +1,22 @@
if (ign.ignited())
{
volScalarField& heau = thermo.heu();
solve
(
betav*fvm::ddt(rho, heau) + mvConvection->fvmDiv(phi, heau)
+ (betav*fvc::ddt(rho, K) + fvc::div(phi, K))*rho/thermo.rhou()
+ (
heau.name() == "eau"
? fvc::div(phi, volScalarField("Ep", p/rho))*rho/thermo.rhou()
: -betav*dpdt*rho/thermo.rhou()
)
- fvm::laplacian(Db, heau)
// These terms cannot be used in partially-premixed combustion due to
// the resultant inconsistency between ft and heau transport.
// A possible solution would be to solve for ftu as well as ft.
//- fvm::div(muEff*fvc::grad(b)/(b + 0.001), heau)
//+ fvm::Sp(fvc::div(muEff*fvc::grad(b)/(b + 0.001)), heau)
);
}

View File

@ -123,12 +123,12 @@ int main(int argc, char *argv[])
{
#include "bEqn.H"
#include "ftEqn.H"
#include "hauEqn.H"
#include "haEqn.H"
#include "EauEqn.H"
#include "EaEqn.H"
if (!ign.ignited())
{
hau == ha;
thermo.heu() == thermo.he();
}
#include "pEqn.H"

View File

@ -5,6 +5,8 @@
psiuReactionThermo::New(mesh)
);
psiuReactionThermo& thermo = pThermo();
thermo.validate(args.executable(), "ha", "ea");
basicMultiComponentMixture& composition = thermo.composition();
volScalarField rho
@ -22,14 +24,10 @@
volScalarField& p = thermo.p();
const volScalarField& psi = thermo.psi();
volScalarField& ha = thermo.he();
volScalarField& hau = thermo.heu();
volScalarField& b = composition.Y("b");
Info<< "min(b) = " << min(b).value() << endl;
//const volScalarField& T = thermo->T();
Info<< "\nReading field U\n" << endl;
volVectorField U
(
@ -196,6 +194,6 @@
}
fields.add(b);
fields.add(ha);
fields.add(hau);
fields.add(thermo.he());
fields.add(thermo.heu());
flameWrinkling->addXi(fields);

View File

@ -1,13 +0,0 @@
{
solve
(
betav*fvm::ddt(rho, ha)
+ mvConvection->fvmDiv(phi, ha)
- fvm::laplacian(Db, ha)
==
betav*dpdt
- betav*(fvc::ddt(rho, K) + fvc::div(phi, K))
);
thermo.correct();
}

View File

@ -1,18 +0,0 @@
if (ign.ignited())
{
solve
(
betav*fvm::ddt(rho, hau)
+ mvConvection->fvmDiv(phi, hau)
- fvm::laplacian(Db, hau)
// These terms cannot be used in partially-premixed combustion due to
// the resultant inconsistency between ft and hau transport.
// A possible solution would be to solve for ftu as well as ft.
//- fvm::div(muEff*fvc::grad(b)/(b + 0.001), hau)
//+ fvm::Sp(fvc::div(muEff*fvc::grad(b)/(b + 0.001)), hau)
==
betav*(dpdt - (fvc::ddt(rho, K) + fvc::div(phi, K)))*rho/thermo.rhou()
);
}

View File

@ -0,0 +1,20 @@
{
volScalarField& hea = thermo.he();
fvScalarMatrix EaEqn
(
fvm::ddt(rho, hea) + mvConvection->fvmDiv(phi, hea)
+ fvc::ddt(rho, K) + fvc::div(phi, K)
+ (
hea.name() == "ea"
? fvc::div(phi, volScalarField("Ep", p/rho))
: -dpdt
)
- fvm::laplacian(turbulence->alphaEff(), hea)
);
EaEqn.relax();
EaEqn.solve();
thermo.correct();
}

View File

@ -0,0 +1,22 @@
if (ign.ignited())
{
volScalarField& heau = thermo.heu();
solve
(
fvm::ddt(rho, heau) + mvConvection->fvmDiv(phi, heau)
+ (fvc::ddt(rho, K) + fvc::div(phi, K))*rho/thermo.rhou()
+ (
heau.name() == "eau"
? fvc::div(phi, volScalarField("Ep", p/rho))*rho/thermo.rhou()
: -dpdt*rho/thermo.rhou()
)
- fvm::laplacian(turbulence->alphaEff(), heau)
// These terms cannot be used in partially-premixed combustion due to
// the resultant inconsistency between ft and heau transport.
// A possible solution would be to solve for ftu as well as ft.
//- fvm::div(muEff*fvc::grad(b)/(b + 0.001), heau)
//+ fvm::Sp(fvc::div(muEff*fvc::grad(b)/(b + 0.001)), heau)
);
}

View File

@ -97,12 +97,12 @@ int main(int argc, char *argv[])
#include "ftEqn.H"
#include "bEqn.H"
#include "hauEqn.H"
#include "haEqn.H"
#include "EauEqn.H"
#include "EaEqn.H"
if (!ign.ignited())
{
hau == ha;
thermo.heu() == thermo.he();
}
// --- Pressure corrector loop

View File

@ -5,6 +5,8 @@
psiuReactionThermo::New(mesh)
);
psiuReactionThermo& thermo = pThermo();
thermo.validate(args.executable(), "ha", "ea");
basicMultiComponentMixture& composition = thermo.composition();
volScalarField rho
@ -22,14 +24,10 @@
volScalarField& p = thermo.p();
const volScalarField& psi = thermo.psi();
volScalarField& ha = thermo.he();
volScalarField& hau = thermo.heu();
volScalarField& b = composition.Y("b");
Info<< "min(b) = " << min(b).value() << endl;
const volScalarField& T = thermo.T();
Info<< "\nReading field U\n" << endl;
volVectorField U
@ -138,5 +136,5 @@
}
fields.add(b);
fields.add(ha);
fields.add(hau);
fields.add(thermo.he());
fields.add(thermo.heu());

View File

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

View File

@ -1,18 +0,0 @@
if (ign.ignited())
{
solve
(
fvm::ddt(rho, hau)
+ mvConvection->fvmDiv(phi, hau)
- fvm::laplacian(turbulence->alphaEff(), hau)
// These terms cannot be used in partially-premixed combustion due to
// the resultant inconsistency between ft and hau transport.
// A possible solution would be to solve for ftu as well as ft.
//- fvm::div(muEff*fvc::grad(b)/(b + 0.001), hau)
//+ fvm::Sp(fvc::div(muEff*fvc::grad(b)/(b + 0.001)), hau)
==
(dpdt - (fvc::ddt(rho, K) + fvc::div(phi, K)))*rho/thermo.rhou()
);
}

View File

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

View File

@ -1,10 +1,12 @@
{
volScalarField& h = thermo.he();
if (constProp == "volume")
{
hs[0] = u0 + p[0]/rho[0] + integratedHeat;
h[0] = u0 + p[0]/rho[0] + integratedHeat;
}
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)
{
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();
rho = thermo.rho();
scalar rho0 = rho[0];
scalar u0 = hs0 - p0/rho0;
scalar u0 = h0 - p0/rho0;
scalar R0 = p0/(rho0*T0);
Rspecific[0] = R0;

View File

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

View File

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

View File

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

@ -103,12 +103,12 @@ int main(int argc, char *argv[])
#include "ftEqn.H"
#include "bEqn.H"
#include "hauEqn.H"
#include "haEqn.H"
#include "EauEqn.H"
#include "EaEqn.H"
if (!ign.ignited())
{
hau == ha;
thermo.heu() == thermo.he();
}
// --- Pressure corrector loop

View File

@ -1,5 +1,6 @@
Info<< "Mean pressure:" << p.weightedAverage(mesh.V()).value() << endl;
Info<< "Mean temperature:" << T.weightedAverage(mesh.V()).value() << endl;
Info<< "Mean temperature:" << thermo.T().weightedAverage(mesh.V()).value()
<< endl;
Info<< "Mean u':"
<< (sqrt((2.0/3.0)*turbulence->k()))().weightedAverage(mesh.V()).value()
<< endl;
@ -7,7 +8,7 @@ Info<< "Mean u':"
logSummaryFile
<< runTime.theta() << tab
<< p.weightedAverage(mesh.V()).value() << tab
<< T.weightedAverage(mesh.V()).value() << tab
<< thermo.T().weightedAverage(mesh.V()).value() << tab
<< (sqrt((2.0/3.0)*turbulence->k()))().weightedAverage(mesh.V()).value()
<< tab
<< 1 - b.weightedAverage(mesh.V()).value()

View File

@ -47,22 +47,27 @@ tmp<fv::convectionScheme<scalar> > mvConvection
Y[inertIndex] = scalar(1) - Yt;
Y[inertIndex].max(0.0);
fvScalarMatrix hsEqn
volScalarField& he = thermo.he();
fvScalarMatrix EEqn
(
fvm::ddt(rho, hs)
+ mvConvection->fvmDiv(phi, hs)
- fvm::laplacian(turbulence->alphaEff(), hs)
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)
==
dpdt
- (fvc::ddt(rho, K) + fvc::div(phi, K))
+ combustion->Sh()
combustion->Sh()
+ radiation->Sh(thermo)
+ parcels.Sh(hs)
+ parcels.Sh(he)
+ surfaceFilm.Sh()
);
hsEqn.relax();
hsEqn.solve();
EEqn.relax();
EEqn.solve();
thermo.correct();

View File

@ -11,6 +11,7 @@
Info<< "Reading thermophysical properties\n" << endl;
psiReactionThermo& thermo = combustion->thermo();
thermo.validate(args.executable(), "h", "e");
SLGThermo slgThermo(mesh, thermo);
@ -34,7 +35,6 @@
);
volScalarField& p = thermo.p();
volScalarField& hs = thermo.he();
const volScalarField& T = thermo.T();
const volScalarField& psi = thermo.psi();
@ -128,7 +128,7 @@
{
fields.add(Y[i]);
}
fields.add(hs);
fields.add(thermo.he());
IOdictionary additionalControlsDict
(

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -94,7 +94,7 @@ int main(int argc, char *argv[])
while (pimple.loop())
{
#include "UEqn.H"
#include "YhsEqn.H"
#include "YEEqn.H"
// --- Pressure corrector loop
while (pimple.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();
thermo.validate(args.executable(), "h", "e");
basicMultiComponentMixture& composition = thermo.composition();
PtrList<volScalarField>& Y = composition.Y();
@ -40,7 +41,6 @@ volVectorField U
volScalarField& p = thermo.p();
const volScalarField& psi = thermo.psi();
volScalarField& hs = thermo.he();
const volScalarField& T = thermo.T();
#include "compressibleCreatePhi.H"
@ -84,7 +84,7 @@ forAll(Y, i)
{
fields.add(Y[i]);
}
fields.add(hs);
fields.add(thermo.he());
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

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -70,7 +70,7 @@ int main(int argc, char *argv[])
{
#include "UEqn.H"
#include "YEqn.H"
#include "hsEqn.H"
#include "EEqn.H"
// --- Pressure corrector loop
while (pimple.correct())

View File

@ -1,4 +1,5 @@
EXE_INC = \
-I../reactingFoam \
-I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
-I$(LIB_SRC)/thermophysicalModels/specie/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();
thermo.validate(args.executable(), "h", "e");
basicMultiComponentMixture& composition = thermo.composition();
PtrList<volScalarField>& Y = composition.Y();
@ -40,7 +41,6 @@ volVectorField U
volScalarField& p = thermo.p();
const volScalarField& psi = thermo.psi();
volScalarField& hs = thermo.he();
const volScalarField& T = thermo.T();
@ -86,7 +86,7 @@ forAll(Y, i)
{
fields.add(Y[i]);
}
fields.add(hs);
fields.add(thermo.he());
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 "YEqn.H"
#include "hsEqn.H"
#include "EEqn.H"
// --- Pressure corrector loop
while (pimple.correct())

View File

@ -0,0 +1,20 @@
{
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(turbulence->alphaEff(), he)
);
EEqn.relax();
EEqn.solve();
thermo.correct();
}

View File

@ -4,7 +4,6 @@ tmp<fvVectorMatrix> UEqn
(
fvm::ddt(rho, U)
+ fvm::div(phi, U)
- fvm::Sp(fvc::ddt(rho) + fvc::div(phi), U)
+ turbulence->divDevRhoReff(U)
);

View File

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

View File

@ -1,20 +0,0 @@
{
fvScalarMatrix hEqn
(
fvm::ddt(rho, h)
+ fvm::div(phi, h)
- fvm::Sp(fvc::ddt(rho) + fvc::div(phi), h)
- fvm::laplacian(turbulence->alphaEff(), h)
==
dpdt
- (
fvc::ddt(rho, K) + fvc::div(phi, K)
- (fvc::ddt(rho) + 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())
{
#include "UEqn.H"
#include "hEqn.H"
#include "EEqn.H"
// --- Pressure corrector loop
while (pimple.correct())

View File

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

View File

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

View File

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

View File

@ -0,0 +1,19 @@
{
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(turbulence->alphaEff(), he)
);
EEqn.relax();
EEqn.solve();
thermo.correct();
}

View File

@ -3,7 +3,6 @@
tmp<fvVectorMatrix> UEqn
(
fvm::div(phi, U)
- fvm::Sp(fvc::div(phi), U)
+ turbulence->divDevRhoReff(U)
);

View File

@ -5,6 +5,7 @@
psiThermo::New(mesh)
);
psiThermo& thermo = pThermo();
thermo.validate(args.executable(), "h", "e");
volScalarField rho
(
@ -20,7 +21,6 @@
);
volScalarField& p = thermo.p();
volScalarField& e = thermo.he();
const volScalarField& psi = thermo.psi();
Info<< "Reading field U\n" << endl;

View File

@ -1,18 +0,0 @@
{
// Kinetic + pressure energy
volScalarField Ekp("Ekp", 0.5*magSqr(U) + p/rho);
fvScalarMatrix eEqn
(
fvm::div(phi, e)
- fvm::Sp(fvc::div(phi), e)
- fvm::laplacian(turbulence->alphaEff(), e)
==
fvc::div(phi)*Ekp - fvc::div(phi, Ekp)
);
eEqn.relax();
eEqn.solve();
thermo.correct();
}

View File

@ -0,0 +1,21 @@
{
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(turbulence->alphaEff(), he)
);
pZones.addEnergySource(thermo, rho, EEqn);
EEqn.relax();
EEqn.solve();
thermo.correct();
}

View File

@ -5,6 +5,7 @@
rhoThermo::New(mesh)
);
rhoThermo& thermo = pThermo();
thermo.validate(args.executable(), "h", "e");
volScalarField rho
(
@ -20,7 +21,6 @@
);
volScalarField& p = thermo.p();
volScalarField& e = thermo.he();
Info<< "Reading field U\n" << endl;
volVectorField U

View File

@ -1,20 +0,0 @@
{
// Kinetic + pressure energy
volScalarField Ekp("Ekp", 0.5*magSqr(U) + p/rho);
fvScalarMatrix eEqn
(
fvm::div(phi, e)
- fvm::Sp(fvc::div(phi), e)
- fvm::laplacian(turbulence->alphaEff(), e)
==
fvc::div(phi)*Ekp - fvc::div(phi, Ekp)
);
pZones.addEnergySource(thermo, rho, eEqn);
eEqn.relax();
eEqn.solve();
thermo.correct();
}

View File

@ -63,7 +63,7 @@ int main(int argc, char *argv[])
// Pressure-velocity SIMPLE corrector
{
#include "UEqn.H"
#include "eEqn.H"
#include "EEqn.H"
#include "pEqn.H"
}

View File

@ -59,7 +59,7 @@ int main(int argc, char *argv[])
// Pressure-velocity SIMPLE corrector
{
#include "UEqn.H"
#include "eEqn.H"
#include "EEqn.H"
#include "pEqn.H"
}

View File

@ -61,7 +61,7 @@ int main(int argc, char *argv[])
// Velocity-pressure-enthalpy SIMPLEC corrector
{
#include "UEqn.H"
#include "eEqn.H"
#include "EEqn.H"
#include "pEqn.H"
}

View File

@ -0,0 +1,10 @@
{
solve
(
fvm::ddt(rho, e) + fvm::div(phi, e)
+ fvc::ddt(rho, K) + fvc::div(phi, volScalarField("Ekp", K + p/rho))
- fvm::laplacian(turbulence->alphaEff(), e)
);
thermo.correct();
}

View File

@ -5,6 +5,7 @@
psiThermo::New(mesh)
);
psiThermo& thermo = pThermo();
thermo.validate(args.executable(), "e");
volScalarField& p = thermo.p();
volScalarField& e = thermo.he();

View File

@ -1,12 +0,0 @@
{
solve
(
fvm::ddt(rho, e)
+ fvm::div(phi, e)
- fvm::laplacian(turbulence->alphaEff(), e)
==
- (fvc::ddt(rho, K) + fvc::div(phi, volScalarField("Ekp", K + p/rho)))
);
thermo.correct();
}

View File

@ -17,7 +17,8 @@ surfaceScalarField phid
volScalarField Dp("Dp", rho*rAU);
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
// Non-orthogonal pressure corrector loop
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
@ -28,7 +29,7 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
pEqn.solve();
if (nonOrth == nNonOrthCorr)
if (pimple.finalNonOrthogonalIter())
{
phi = pEqn.flux();
}

View File

@ -0,0 +1,11 @@
{
solve
(
fvm::ddt(rho, e) + fvm::div(phi, e)
+ fvc::ddt(rho, K) + fvc::div(phi, K)
+ fvc::div(phi/fvc::interpolate(rho) + mesh.phi(), p, "div(phiv,p)")
- fvm::laplacian(turbulence->alphaEff(), e)
);
thermo.correct();
}

View File

@ -1,12 +0,0 @@
{
solve
(
fvm::ddt(rho, e)
+ fvm::div(phi, e)
- fvm::laplacian(turbulence->alphaEff(), e)
==
- p*fvc::div(phi/fvc::interpolate(rho) + mesh.phi())
);
thermo.correct();
}

View File

@ -34,6 +34,7 @@ Description
#include "psiThermo.H"
#include "turbulenceModel.H"
#include "motionSolver.H"
#include "pimpleControl.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -45,6 +46,8 @@ int main(int argc, char *argv[])
#include "createFields.H"
#include "initContinuityErrs.H"
pimpleControl pimple(mesh);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
@ -62,19 +65,23 @@ int main(int argc, char *argv[])
#include "rhoEqn.H"
#include "UEqn.H"
#include "eEqn.H"
// --- PISO loop
for (int corr=0; corr<nCorr; corr++)
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
#include "pEqn.H"
}
#include "UEqn.H"
#include "EEqn.H"
turbulence->correct();
// --- Pressure corrector loop
while (pimple.correct())
{
#include "pEqn.H"
}
if (pimple.turbCorr())
{
turbulence->correct();
}
}
rho = thermo.rho();

View File

@ -33,6 +33,7 @@ Description
#include "fvCFD.H"
#include "psiThermo.H"
#include "turbulenceModel.H"
#include "pimpleControl.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -44,6 +45,8 @@ int main(int argc, char *argv[])
#include "createFields.H"
#include "initContinuityErrs.H"
pimpleControl pimple(mesh);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
@ -52,20 +55,28 @@ int main(int argc, char *argv[])
{
Info<< "Time = " << runTime.timeName() << nl << endl;
#include "readPISOControls.H"
#include "readTimeControls.H"
#include "compressibleCourantNo.H"
#include "rhoEqn.H"
#include "UEqn.H"
// --- PISO loop
for (int corr=0; corr<nCorr; corr++)
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
#include "eEqn.H"
#include "pEqn.H"
}
#include "UEqn.H"
#include "EEqn.H"
turbulence->correct();
// --- Pressure corrector loop
while (pimple.correct())
{
#include "pEqn.H"
}
if (pimple.turbCorr())
{
turbulence->correct();
}
}
rho = thermo.rho();

View File

@ -31,6 +31,7 @@ Description
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "pimpleControl.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -44,6 +45,8 @@ int main(int argc, char *argv[])
#include "createFields.H"
#include "initContinuityErrs.H"
pimpleControl pimple(mesh);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
@ -52,57 +55,59 @@ int main(int argc, char *argv[])
{
Info<< "Time = " << runTime.timeName() << nl << endl;
#include "readPISOControls.H"
#include "readTimeControls.H"
#include "compressibleCourantNo.H"
#include "rhoEqn.H"
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(phi, U)
- fvm::laplacian(mu, U)
);
solve(UEqn == -fvc::grad(p));
// --- PISO loop
for (int corr=0; corr<nCorr; corr++)
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
volScalarField rAU(1.0/UEqn.A());
U = rAU*UEqn.H();
surfaceScalarField phid
fvVectorMatrix UEqn
(
"phid",
psi
*(
(fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rAU, rho, U, phi)
)
fvm::ddt(rho, U)
+ fvm::div(phi, U)
- fvm::laplacian(mu, U)
);
phi = (rhoO/psi)*phid;
volScalarField Dp("Dp", rho*rAU);
solve(UEqn == -fvc::grad(p));
fvScalarMatrix pEqn
(
fvm::ddt(psi, p)
+ fvc::div(phi)
+ fvm::div(phid, p)
- fvm::laplacian(Dp, p)
);
// --- Pressure corrector loop
while (pimple.correct())
{
volScalarField rAU(1.0/UEqn.A());
U = rAU*UEqn.H();
pEqn.solve();
surfaceScalarField phid
(
"phid",
psi
*(
(fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rAU, rho, U, phi)
)
);
phi += pEqn.flux();
phi = (rhoO/psi)*phid;
volScalarField Dp("Dp", rho*rAU);
#include "compressibleContinuityErrs.H"
fvScalarMatrix pEqn
(
fvm::ddt(psi, p)
+ fvc::div(phi)
+ fvm::div(phid, p)
- fvm::laplacian(Dp, p)
);
U -= rAU*fvc::grad(p);
U.correctBoundaryConditions();
pEqn.solve();
phi += pEqn.flux();
#include "compressibleContinuityErrs.H"
U -= rAU*fvc::grad(p);
U.correctBoundaryConditions();
}
}
rho = rhoO + psi*p;

View File

@ -7,7 +7,6 @@
fvScalarMatrix TEqn
(
fvm::div(phi, T)
- fvm::Sp(fvc::div(phi), T)
- fvm::laplacian(kappaEff, T)
);

View File

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

View File

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

View File

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

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

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

View File

@ -0,0 +1,19 @@
{
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(turbulence->alphaEff(), he)
);
EEqn.relax();
EEqn.solve();
thermo.correct();
}

View File

@ -59,7 +59,7 @@ int main(int argc, char *argv[])
// Pressure-velocity SIMPLE corrector
{
#include "UEqn.H"
#include "hEqn.H"
#include "EEqn.H"
#include "pEqn.H"
}

View File

@ -0,0 +1,22 @@
{
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(turbulence->alphaEff(), he)
==
radiation->Sh(thermo)
);
EEqn.relax();
EEqn.solve();
thermo.correct();
radiation->correct();
}

View File

@ -1,5 +1,5 @@
EXE_INC = \
-I../buoyantSimpleFoam \
-I.. \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \
-I$(LIB_SRC)/turbulenceModels \

View File

@ -62,7 +62,7 @@ int main(int argc, char *argv[])
// Pressure-velocity SIMPLE corrector
{
#include "UEqn.H"
#include "hEqn.H"
#include "EEqn.H"
#include "pEqn.H"
}

View File

@ -5,6 +5,7 @@
psiThermo::New(mesh)
);
psiThermo& thermo = pThermo();
thermo.validate(args.executable(), "h", "e");
volScalarField rho
(
@ -20,7 +21,6 @@
);
volScalarField& p = thermo.p();
volScalarField& h = thermo.he();
const volScalarField& psi = thermo.psi();
Info<< "Reading field U\n" << endl;

View File

@ -1,15 +0,0 @@
{
fvScalarMatrix hEqn
(
fvm::div(phi, h)
- fvm::Sp(fvc::div(phi), h)
- fvm::laplacian(turbulence->alphaEff(), h)
==
- fvc::div(phi, 0.5*magSqr(U), "div(phi,K)")
);
hEqn.relax();
hEqn.solve();
thermo.correct();
}

View File

@ -1,19 +0,0 @@
{
fvScalarMatrix hEqn
(
fvm::div(phi, h)
- fvm::Sp(fvc::div(phi), h)
- fvm::laplacian(turbulence->alphaEff(), h)
==
- fvc::div(phi, 0.5*magSqr(U), "div(phi,K)")
+ radiation->Sh(thermo)
);
hEqn.relax();
hEqn.solve();
thermo.correct();
radiation->correct();
}

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,23 +0,0 @@
{
fvScalarMatrix hEqn
(
fvm::div(phi, h)
- fvm::Sp(fvc::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

@ -7,23 +7,31 @@
volScalarField rAU(1.0/UEqn().A());
surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU));
U = rAU*UEqn().H();
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn().H();
UEqn.clear();
phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
bool closedVolume = adjustPhi(phi, U, p_rgh);
surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf())
);
bool closedVolume = adjustPhi(phiHbyA, U, p_rgh);
phiHbyA += phig;
dimensionedScalar compressibility = fvc::domainIntegrate(psi);
bool compressible = (compressibility.value() > SMALL);
surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
phi += phig;
// Solve pressure
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix p_rghEqn
(
fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phi)
fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phiHbyA)
);
p_rghEqn.setReference
@ -37,14 +45,14 @@
if (nonOrth == nNonOrthCorr)
{
// Calculate the conservative fluxes
phi -= p_rghEqn.flux();
phi = phiHbyA - p_rghEqn.flux();
// Explicitly relax pressure for momentum corrector
p_rgh.relax();
// Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure
U += rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rhorAUf);
U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rhorAUf);
U.correctBoundaryConditions();
}
}

View File

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

View File

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

View File

@ -2,7 +2,6 @@
fvScalarMatrix hPorousEqn
(
fvm::div(porousPhi, porousH)
- fvm::Sp(fvc::div(porousPhi), porousH)
- fvm::laplacian(turbPorous.alphaEff(), porousH)
==
- fvc::div(porousPhi, 0.5*magSqr(porousU), "div(phi,K)")

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];
rhoThermo& thermo = thermoFluid[i];
thermo.validate(args.executable(), "h", "e");
volScalarField& rho = rhoFluid[i];
volVectorField& U = UFluid[i];
surfaceScalarField& phi = phiFluid[i];
@ -11,7 +13,6 @@
volScalarField& p = thermo.p();
const volScalarField& psi = thermo.psi();
volScalarField& h = thermo.he();
volScalarField& p_rgh = p_rghFluid[i];
const volScalarField& gh = ghFluid[i];

View File

@ -9,8 +9,7 @@ if (oCorr == 0)
}
#include "UEqn.H"
#include "hEqn.H"
#include "EEqn.H"
// --- PISO loop
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

@ -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 "YEqn.H"
#include "hsEqn.H"
#include "EEqn.H"
// --- Pressure corrector loop
while (pimple.correct())

View File

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

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