zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

ENH: Updates to ODEChemistryModel.[CH]

Browse Source
This commit is contained in:
andy
2010-08-02 15:19:37 +01:00
parent 9c5ecc1a90
commit 1c8fd73a79
2 changed files with 91 additions and 101 deletions
Download Patch File Download Diff File Expand all files Collapse all files

186
src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.C
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) 2009-2009 OpenCFD Ltd. \\ / A nd | Copyright (C) 2009-2010 OpenCFD Ltd.
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -93,7 +93,8 @@ Foam::ODEChemistryModel<CompType, ThermoType>::~ODEChemistryModel()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class CompType, class ThermoType> template<class CompType, class ThermoType>
Foam::scalarField Foam::ODEChemistryModel<CompType, ThermoType>::omega Foam::tmp<Foam::scalarField>
Foam::ODEChemistryModel<CompType, ThermoType>::omega
( (
const scalarField& c, const scalarField& c,
const scalar T, const scalar T,
@ -103,7 +104,8 @@ Foam::scalarField Foam::ODEChemistryModel<CompType, ThermoType>::omega
scalar pf, cf, pr, cr; scalar pf, cf, pr, cr;
label lRef, rRef; label lRef, rRef;
scalarField om(nEqns(), 0.0); tmp<scalarField> tom(new scalarField(nEqns(), 0.0));
scalarField& om = tom();
forAll(reactions_, i) forAll(reactions_, i)
{ {
@ -116,20 +118,20 @@ Foam::scalarField Foam::ODEChemistryModel<CompType, ThermoType>::omega
forAll(R.lhs(), s) forAll(R.lhs(), s)
{ {
label si = R.lhs()[s].index; const label si = R.lhs()[s].index;
scalar sl = R.lhs()[s].stoichCoeff; const scalar sl = R.lhs()[s].stoichCoeff;
om[si] -= sl*omegai; om[si] -= sl*omegai;
} }
forAll(R.rhs(), s) forAll(R.rhs(), s)
{ {
label si = R.rhs()[s].index; const label si = R.rhs()[s].index;
scalar sr = R.rhs()[s].stoichCoeff; const scalar sr = R.rhs()[s].stoichCoeff;
om[si] += sr*omegai; om[si] += sr*omegai;
} }
} }
return om; return tom;
} }
@ -154,8 +156,8 @@ Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::omega
c2[i] = max(0.0, c[i]); c2[i] = max(0.0, c[i]);
} }
scalar kf = R.kf(T, p, c2); const scalar kf = R.kf(T, p, c2);
scalar kr = R.kr(kf, T, p, c2); const scalar kr = R.kr(kf, T, p, c2);
pf = 1.0; pf = 1.0;
pr = 1.0; pr = 1.0;
@ -169,26 +171,26 @@ Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::omega
pf = kf; pf = kf;
for (label s=1; s<Nl; s++) for (label s=1; s<Nl; s++)
{ {
label si = R.lhs()[s].index; const label si = R.lhs()[s].index;
if (c[si] < c[lRef]) if (c[si] < c[lRef])
{ {
scalar exp = R.lhs()[slRef].exponent; const scalar exp = R.lhs()[slRef].exponent;
pf *= pow(max(0.0, c[lRef]), exp); pf *= pow(max(0.0, c[lRef]), exp);
lRef = si; lRef = si;
slRef = s; slRef = s;
} }
else else
{ {
scalar exp = R.lhs()[s].exponent; const scalar exp = R.lhs()[s].exponent;
pf *= pow(max(0.0, c[si]), exp); pf *= pow(max(0.0, c[si]), exp);
} }
} }
cf = max(0.0, c[lRef]); cf = max(0.0, c[lRef]);
{ {
scalar exp = R.lhs()[slRef].exponent; const scalar exp = R.lhs()[slRef].exponent;
if (exp<1.0) if (exp < 1.0)
{ {
if (cf > SMALL) if (cf > SMALL)
{ {
@ -212,25 +214,25 @@ Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::omega
pr = kr; pr = kr;
for (label s=1; s<Nr; s++) for (label s=1; s<Nr; s++)
{ {
label si = R.rhs()[s].index; const label si = R.rhs()[s].index;
if (c[si] < c[rRef]) if (c[si] < c[rRef])
{ {
scalar exp = R.rhs()[srRef].exponent; const scalar exp = R.rhs()[srRef].exponent;
pr *= pow(max(0.0, c[rRef]), exp); pr *= pow(max(0.0, c[rRef]), exp);
rRef = si; rRef = si;
srRef = s; srRef = s;
} }
else else
{ {
scalar exp = R.rhs()[s].exponent; const scalar exp = R.rhs()[s].exponent;
pr *= pow(max(0.0, c[si]), exp); pr *= pow(max(0.0, c[si]), exp);
} }
} }
cr = max(0.0, c[rRef]); cr = max(0.0, c[rRef]);
{ {
scalar exp = R.rhs()[srRef].exponent; const scalar exp = R.rhs()[srRef].exponent;
if (exp<1.0) if (exp < 1.0)
{ {
if (cr>SMALL) if (cr>SMALL)
{ {
@ -259,8 +261,8 @@ void Foam::ODEChemistryModel<CompType, ThermoType>::derivatives
scalarField& dcdt scalarField& dcdt
) const ) const
{ {
scalar T = c[nSpecie_]; const scalar T = c[nSpecie_];
scalar p = c[nSpecie_ + 1]; const scalar p = c[nSpecie_ + 1];
dcdt = omega(c, T, p); dcdt = omega(c, T, p);
@ -270,16 +272,16 @@ void Foam::ODEChemistryModel<CompType, ThermoType>::derivatives
scalar cSum = 0.0; scalar cSum = 0.0;
for (label i=0; i<nSpecie_; i++) for (label i=0; i<nSpecie_; i++)
{ {
scalar W = specieThermo_[i].W(); const scalar W = specieThermo_[i].W();
cSum += c[i]; cSum += c[i];
rho += W*c[i]; rho += W*c[i];
} }
scalar mw = rho/cSum; const scalar mw = rho/cSum;
scalar cp = 0.0; scalar cp = 0.0;
for (label i=0; i<nSpecie_; i++) for (label i=0; i<nSpecie_; i++)
{ {
scalar cpi = specieThermo_[i].cp(T); const scalar cpi = specieThermo_[i].cp(T);
scalar Xi = c[i]/rho; const scalar Xi = c[i]/rho;
cp += Xi*cpi; cp += Xi*cpi;
} }
cp /= mw; cp /= mw;
@ -287,18 +289,18 @@ void Foam::ODEChemistryModel<CompType, ThermoType>::derivatives
scalar dT = 0.0; scalar dT = 0.0;
for (label i=0; i<nSpecie_; i++) for (label i=0; i<nSpecie_; i++)
{ {
scalar hi = specieThermo_[i].h(T); const scalar hi = specieThermo_[i].h(T);
dT += hi*dcdt[i]; dT += hi*dcdt[i];
} }
dT /= rho*cp; dT /= rho*cp;
// limit the time-derivative, this is more stable for the ODE // limit the time-derivative, this is more stable for the ODE
// solver when calculating the allowed time step // solver when calculating the allowed time step
scalar dtMag = min(500.0, mag(dT)); const scalar dTLimited = min(500.0, mag(dT));
dcdt[nSpecie_] = -dT*dtMag/(mag(dT) + 1.0e-10); dcdt[nSpecie_] = -dT*dTLimited/(mag(dT) + 1.0e-10);
// dp/dt = ... // dp/dt = ...
dcdt[nSpecie_+1] = 0.0; dcdt[nSpecie_ + 1] = 0.0;
} }
@ -311,11 +313,11 @@ void Foam::ODEChemistryModel<CompType, ThermoType>::jacobian
scalarSquareMatrix& dfdc scalarSquareMatrix& dfdc
) const ) const
{ {
scalar T = c[nSpecie_]; const scalar T = c[nSpecie_];
scalar p = c[nSpecie_ + 1]; const scalar p = c[nSpecie_ + 1];
scalarField c2(nSpecie_, 0.0); scalarField c2(nSpecie_, 0.0);
for (label i=0; i<nSpecie_; i++) forAll(c2, i)
{ {
c2[i] = max(c[i], 0.0); c2[i] = max(c[i], 0.0);
} }
@ -335,22 +337,22 @@ void Foam::ODEChemistryModel<CompType, ThermoType>::jacobian
{ {
const Reaction<ThermoType>& R = reactions_[ri]; const Reaction<ThermoType>& R = reactions_[ri];
scalar kf0 = R.kf(T, p, c2); const scalar kf0 = R.kf(T, p, c2);
scalar kr0 = R.kr(T, p, c2); const scalar kr0 = R.kr(T, p, c2);
forAll(R.lhs(), j) forAll(R.lhs(), j)
{ {
label sj = R.lhs()[j].index; const label sj = R.lhs()[j].index;
scalar kf = kf0; scalar kf = kf0;
forAll(R.lhs(), i) forAll(R.lhs(), i)
{ {
label si = R.lhs()[i].index; const label si = R.lhs()[i].index;
scalar el = R.lhs()[i].exponent; const scalar el = R.lhs()[i].exponent;
if (i == j) if (i == j)
{ {
if (el < 1.0) if (el < 1.0)
{ {
if (c2[si]>SMALL) if (c2[si] > SMALL)
{ {
kf *= el*pow(c2[si] + VSMALL, el - 1.0); kf *= el*pow(c2[si] + VSMALL, el - 1.0);
} }
@ -372,31 +374,31 @@ void Foam::ODEChemistryModel<CompType, ThermoType>::jacobian
forAll(R.lhs(), i) forAll(R.lhs(), i)
{ {
label si = R.lhs()[i].index; const label si = R.lhs()[i].index;
scalar sl = R.lhs()[i].stoichCoeff; const scalar sl = R.lhs()[i].stoichCoeff;
dfdc[si][sj] -= sl*kf; dfdc[si][sj] -= sl*kf;
} }
forAll(R.rhs(), i) forAll(R.rhs(), i)
{ {
label si = R.rhs()[i].index; const label si = R.rhs()[i].index;
scalar sr = R.rhs()[i].stoichCoeff; const scalar sr = R.rhs()[i].stoichCoeff;
dfdc[si][sj] += sr*kf; dfdc[si][sj] += sr*kf;
} }
} }
forAll(R.rhs(), j) forAll(R.rhs(), j)
{ {
label sj = R.rhs()[j].index; const label sj = R.rhs()[j].index;
scalar kr = kr0; scalar kr = kr0;
forAll(R.rhs(), i) forAll(R.rhs(), i)
{ {
label si = R.rhs()[i].index; const label si = R.rhs()[i].index;
scalar er = R.rhs()[i].exponent; const scalar er = R.rhs()[i].exponent;
if (i==j) if (i == j)
{ {
if (er<1.0) if (er < 1.0)
{ {
if (c2[si]>SMALL) if (c2[si] > SMALL)
{ {
kr *= er*pow(c2[si] + VSMALL, er - 1.0); kr *= er*pow(c2[si] + VSMALL, er - 1.0);
} }
@ -418,25 +420,25 @@ void Foam::ODEChemistryModel<CompType, ThermoType>::jacobian
forAll(R.lhs(), i) forAll(R.lhs(), i)
{ {
label si = R.lhs()[i].index; const label si = R.lhs()[i].index;
scalar sl = R.lhs()[i].stoichCoeff; const scalar sl = R.lhs()[i].stoichCoeff;
dfdc[si][sj] += sl*kr; dfdc[si][sj] += sl*kr;
} }
forAll(R.rhs(), i) forAll(R.rhs(), i)
{ {
label si = R.rhs()[i].index; const label si = R.rhs()[i].index;
scalar sr = R.rhs()[i].stoichCoeff; const scalar sr = R.rhs()[i].stoichCoeff;
dfdc[si][sj] -= sr*kr; dfdc[si][sj] -= sr*kr;
} }
} }
} }
// calculate the dcdT elements numerically // calculate the dcdT elements numerically
scalar delta = 1.0e-8; const scalar delta = 1.0e-8;
scalarField dcdT0 = omega(c2, T - delta, p); const scalarField dcdT0 = omega(c2, T - delta, p);
scalarField dcdT1 = omega(c2, T + delta, p); const scalarField dcdT1 = omega(c2, T + delta, p);
for (label i=0; i<nEqns(); i++) for (label i = 0; i < nEqns(); i++)
{ {
dfdc[i][nSpecie()] = 0.5*(dcdT1[i] - dcdT0[i])/delta; dfdc[i][nSpecie()] = 0.5*(dcdT1[i] - dcdT0[i])/delta;
} }
@ -485,7 +487,7 @@ Foam::ODEChemistryModel<CompType, ThermoType>::tc() const
scalarField& tc = ttc(); scalarField& tc = ttc();
label nReaction = reactions_.size(); const label nReaction = reactions_.size();
if (this->chemistry_) if (this->chemistry_)
{ {
@ -626,34 +628,31 @@ void Foam::ODEChemistryModel<CompType, ThermoType>::calculate()
this->thermo().rho() this->thermo().rho()
); );
if (this->mesh().changing())
{
for (label i=0; i<nSpecie_; i++) for (label i=0; i<nSpecie_; i++)
{ {
RR_[i].setSize(rho.size()); RR_[i].setSize(rho.size());
RR_[i] = 0.0;
}
} }
if (this->chemistry_) if (this->chemistry_)
{ {
forAll(rho, celli) forAll(rho, celli)
{ {
const scalar rhoi = rho[celli];
const scalar Ti = this->thermo().T()[celli];
const scalar pi = this->thermo().p()[celli];
scalarField c(nSpecie_, 0.0);
for (label i=0; i<nSpecie_; i++) for (label i=0; i<nSpecie_; i++)
{ {
RR_[i][celli] = 0.0; const scalar Yi = Y_[i][celli];
}
scalar rhoi = rho[celli];
scalar Ti = this->thermo().T()[celli];
scalar pi = this->thermo().p()[celli];
scalarField c(nSpecie_);
scalarField dcdt(nEqns(), 0.0);
for (label i=0; i<nSpecie_; i++)
{
scalar Yi = Y_[i][celli];
c[i] = rhoi*Yi/specieThermo_[i].W(); c[i] = rhoi*Yi/specieThermo_[i].W();
} }
dcdt = omega(c, Ti, pi); const scalarField dcdt = omega(c, Ti, pi);
for (label i=0; i<nSpecie_; i++) for (label i=0; i<nSpecie_; i++)
{ {
@ -671,6 +670,8 @@ Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::solve
const scalar deltaT const scalar deltaT
) )
{ {
scalar deltaTMin = GREAT;
const volScalarField rho const volScalarField rho
( (
IOobject IOobject
@ -685,35 +686,33 @@ Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::solve
this->thermo().rho() this->thermo().rho()
); );
for (label i=0; i<nSpecie_; i++) if (this->mesh().changing())
{ {
RR_[i].setSize(rho.size()); for (label i = 0; i < nSpecie_; i++)
{
RR_[i].setSize(this->mesh().nCells());
RR_[i] = 0.0;
}
} }
if (!this->chemistry_) if (!this->chemistry_)
{ {
return GREAT; return deltaTMin;
} }
scalar deltaTMin = GREAT;
tmp<volScalarField> thc = this->thermo().hc(); tmp<volScalarField> thc = this->thermo().hc();
const scalarField& hc = thc(); const scalarField& hc = thc();
forAll(rho, celli) forAll(rho, celli)
{ {
for (label i=0; i<nSpecie_; i++) const scalar rhoi = rho[celli];
{ const scalar hi = this->thermo().hs()[celli] + hc[celli];
RR_[i][celli] = 0.0; const scalar pi = this->thermo().p()[celli];
}
scalar rhoi = rho[celli];
scalar Ti = this->thermo().T()[celli]; scalar Ti = this->thermo().T()[celli];
scalar hi = this->thermo().hs()[celli] + hc[celli];
scalar pi = this->thermo().p()[celli];
scalarField c(nSpecie_); scalarField c(nSpecie_, 0.0);
scalarField c0(nSpecie_); scalarField c0(nSpecie_, 0.0);
scalarField dc(nSpecie_, 0.0); scalarField dc(nSpecie_, 0.0);
for (label i=0; i<nSpecie_; i++) for (label i=0; i<nSpecie_; i++)
@ -722,21 +721,20 @@ Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::solve
} }
c0 = c; c0 = c;
// initialise timing parameters
scalar t = t0; scalar t = t0;
scalar tauC = this->deltaTChem_[celli]; scalar tauC = this->deltaTChem_[celli];
scalar dt = min(deltaT, tauC); scalar dt = min(deltaT, tauC);
scalar timeLeft = deltaT; scalar timeLeft = deltaT;
// calculate the chemical source terms // calculate the chemical source terms
scalar cTot = 0.0;
while (timeLeft > SMALL) while (timeLeft > SMALL)
{ {
tauC = solver().solve(c, Ti, pi, t, dt); tauC = solver().solve(c, Ti, pi, t, dt);
t += dt; t += dt;
// update the temperature // update the temperature
cTot = sum(c); scalar cTot = sum(c);
ThermoType mixture(0.0*specieThermo_[0]); ThermoType mixture(0.0*specieThermo_[0]);
for (label i=0; i<nSpecie_; i++) for (label i=0; i<nSpecie_; i++)
{ {
@ -746,19 +744,11 @@ Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::solve
timeLeft -= dt; timeLeft -= dt;
this->deltaTChem_[celli] = tauC; this->deltaTChem_[celli] = tauC;
dt = min(timeLeft, tauC); dt = max(SMALL, min(timeLeft, tauC));
dt = max(dt, SMALL);
} }
deltaTMin = min(tauC, deltaTMin); deltaTMin = min(tauC, deltaTMin);
dc = c - c0; dc = c - c0;
scalar WTot = 0.0;
for (label i=0; i<nSpecie_; i++)
{
WTot += c[i]*specieThermo_[i].W();
}
WTot /= cTot;
for (label i=0; i<nSpecie_; i++) for (label i=0; i<nSpecie_; i++)
{ {
RR_[i][celli] = dc[i]*specieThermo_[i].W()/deltaT; RR_[i][celli] = dc[i]*specieThermo_[i].W()/deltaT;

2
src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.H
View File

@ -141,7 +141,7 @@ public:
inline const chemistrySolver<CompType, ThermoType>& solver() const; inline const chemistrySolver<CompType, ThermoType>& solver() const;
//- dc/dt = omega, rate of change in concentration, for each species //- dc/dt = omega, rate of change in concentration, for each species
virtual scalarField omega virtual tmp<scalarField> omega
( (
const scalarField& c, const scalarField& c,
const scalar T, const scalar T,