32 template<
class CompType,
class ThermoType>
36 const word& compTypeName,
37 const word& thermoTypeName
40 CompType(mesh, thermoTypeName),
53 (this->
thermo()).speciesData()
57 nReaction_(reactions_.size()),
81 Info<<
"ODEChemistryModel: Number of species = " << nSpecie_
82 <<
" and reactions = " << nReaction_ <<
endl;
88 template<
class CompType,
class ThermoType>
95 template<
class CompType,
class ThermoType>
103 scalar pf, cf, pr, cr;
112 scalar omegai = omega
114 R, c, T, p, pf, cf, lRef, pr, cr, rRef
119 label si = R.
lhs()[s].index;
120 scalar sl = R.
lhs()[s].stoichCoeff;
126 label si = R.
rhs()[s].index;
127 scalar sr = R.
rhs()[s].stoichCoeff;
136 template<
class CompType,
class ThermoType>
152 for (label i=0; i<nSpecie_; i++)
154 c2[i] =
max(0.0, c[i]);
157 scalar kf = R.
kf(T, p, c2);
158 scalar kr = R.
kr(kf, T, p, c2);
167 lRef = R.
lhs()[slRef].index;
170 for (label s=1; s<Nl; s++)
172 label si = R.
lhs()[s].index;
176 scalar
exp = R.
lhs()[slRef].exponent;
177 pf *=
pow(
max(0.0, c[lRef]), exp);
183 scalar
exp = R.
lhs()[s].exponent;
184 pf *=
pow(
max(0.0, c[si]), exp);
187 cf =
max(0.0, c[lRef]);
190 scalar
exp = R.
lhs()[slRef].exponent;
195 pf *=
pow(cf, exp - 1.0);
204 pf *=
pow(cf, exp - 1.0);
209 rRef = R.
rhs()[srRef].index;
213 for (label s=1; s<Nr; s++)
215 label si = R.
rhs()[s].index;
218 scalar
exp = R.
rhs()[srRef].exponent;
219 pr *=
pow(
max(0.0, c[rRef]), exp);
225 scalar
exp = R.
rhs()[s].exponent;
226 pr *=
pow(
max(0.0, c[si]), exp);
229 cr =
max(0.0, c[rRef]);
232 scalar
exp = R.
rhs()[srRef].exponent;
237 pr *=
pow(cr, exp - 1.0);
246 pr *=
pow(cr, exp - 1.0);
250 return pf*cf - pr*cr;
254 template<
class CompType,
class ThermoType>
262 scalar T = c[nSpecie_];
263 scalar p = c[nSpecie_ + 1];
265 dcdt = omega(c, T, p);
271 for (label i=0; i<nSpecie_; i++)
273 scalar W = specieThermo_[i].W();
277 scalar mw = rho/cSum;
279 for (label i=0; i<nSpecie_; i++)
281 scalar cpi = specieThermo_[i].cp(T);
282 scalar
Xi = c[i]/
rho;
288 for (label i=0; i<nSpecie_; i++)
290 scalar hi = specieThermo_[i].h(T);
297 scalar dtMag =
min(500.0,
mag(dT));
298 dcdt[nSpecie_] = -dT*dtMag/(
mag(dT) + 1.0e-10);
301 dcdt[nSpecie_+1] = 0.0;
305 template<
class CompType,
class ThermoType>
314 scalar T = c[nSpecie_];
315 scalar p = c[nSpecie_ + 1];
318 for (label i=0; i<nSpecie_; i++)
320 c2[i] =
max(c[i], 0.0);
323 for (label i=0; i<nEqns(); i++)
325 for (label j=0; j<nEqns(); j++)
332 dcdt = omega(c2, T, p);
334 for (label ri=0; ri<reactions_.size(); ri++)
338 scalar kf0 = R.
kf(T, p, c2);
339 scalar kr0 = R.
kr(T, p, c2);
343 label sj = R.
lhs()[j].index;
347 label si = R.
lhs()[i].index;
348 scalar el = R.
lhs()[i].exponent;
355 kf *= el*
pow(c2[si] + VSMALL, el - 1.0);
364 kf *= el*
pow(c2[si], el - 1.0);
369 kf *=
pow(c2[si], el);
375 label si = R.
lhs()[i].index;
376 scalar sl = R.
lhs()[i].stoichCoeff;
377 dfdc[si][sj] -= sl*kf;
381 label si = R.
rhs()[i].index;
382 scalar sr = R.
rhs()[i].stoichCoeff;
383 dfdc[si][sj] += sr*kf;
389 label sj = R.
rhs()[j].index;
393 label si = R.
rhs()[i].index;
394 scalar er = R.
rhs()[i].exponent;
401 kr *= er*
pow(c2[si] + VSMALL, er - 1.0);
410 kr *= er*
pow(c2[si], er - 1.0);
415 kr *=
pow(c2[si], er);
421 label si = R.
lhs()[i].index;
422 scalar sl = R.
lhs()[i].stoichCoeff;
423 dfdc[si][sj] += sl*kr;
427 label si = R.
rhs()[i].index;
428 scalar sr = R.
rhs()[i].stoichCoeff;
429 dfdc[si][sj] -= sr*kr;
435 scalar delta = 1.0e-8;
439 for (label i=0; i<nEqns(); i++)
441 dfdc[i][nSpecie()] = 0.5*(dcdT1[i] - dcdT0[i])/delta;
447 template<
class CompType,
class ThermoType>
451 scalar pf, cf, pr, cr;
482 zeroGradientFvPatchScalarField::typeName
488 label nReaction = reactions_.size();
490 if (this->chemistry_)
494 scalar rhoi = rho[celli];
495 scalar Ti = this->
thermo().T()[celli];
496 scalar
pi = this->
thermo().p()[celli];
500 for (label i=0; i<nSpecie_; i++)
502 scalar Yi = Y_[i][celli];
503 c[i] = rhoi*Yi/specieThermo_[i].W();
513 R, c, Ti, pi, pf, cf, lRef, pr, cr, rRef
518 scalar sr = R.
rhs()[s].stoichCoeff;
519 tc[celli] += sr*pf*cf;
522 tc[celli] = nReaction*cSum/tc[celli];
527 ttc().correctBoundaryConditions();
533 template<
class CompType,
class ThermoType>
544 this->mesh_.time().timeName(),
552 zeroGradientFvPatchScalarField::typeName
556 if (this->chemistry_)
564 scalar hi = specieThermo_[i].Hc();
565 Sh[cellI] -= hi*RR_[i][cellI];
574 template<
class CompType,
class ThermoType>
585 this->mesh_.time().timeName(),
593 zeroGradientFvPatchScalarField::typeName
597 if (this->chemistry_)
607 template<
class CompType,
class ThermoType>
615 template<
class CompType,
class ThermoType>
632 for (label i=0; i<nSpecie_; i++)
634 RR_[i].setSize(rho.
size());
637 if (this->chemistry_)
641 for (label i=0; i<nSpecie_; i++)
646 scalar rhoi = rho[celli];
647 scalar Ti = this->
thermo().T()[celli];
648 scalar
pi = this->
thermo().p()[celli];
653 for (label i=0; i<nSpecie_; i++)
655 scalar Yi = Y_[i][celli];
656 c[i] = rhoi*Yi/specieThermo_[i].W();
659 dcdt = omega(c, Ti, pi);
661 for (label i=0; i<nSpecie_; i++)
663 RR_[i][celli] = dcdt[i]*specieThermo_[i].W();
670 template<
class CompType,
class ThermoType>
691 for (label i=0; i<nSpecie_; i++)
693 RR_[i].setSize(rho.
size());
696 if (!this->chemistry_)
701 scalar deltaTMin = GREAT;
708 for (label i=0; i<nSpecie_; i++)
713 scalar rhoi = rho[celli];
714 scalar Ti = this->
thermo().T()[celli];
715 scalar hi = this->
thermo().hs()[celli] + hc[celli];
716 scalar
pi = this->
thermo().p()[celli];
722 for (label i=0; i<nSpecie_; i++)
724 c[i] = rhoi*Y_[i][celli]/specieThermo_[i].W();
729 scalar tauC = this->deltaTChem_[celli];
730 scalar dt =
min(deltaT, tauC);
731 scalar timeLeft = deltaT;
736 while (timeLeft > SMALL)
738 tauC = solver().solve(c, Ti, pi, t, dt);
744 for (label i=0; i<nSpecie_; i++)
746 mixture += (c[i]/cTot)*specieThermo_[i];
748 Ti = mixture.TH(hi, Ti);
751 this->deltaTChem_[celli] = tauC;
752 dt =
min(timeLeft, tauC);
755 deltaTMin =
min(tauC, deltaTMin);
759 for (label i=0; i<nSpecie_; i++)
761 WTot += c[i]*specieThermo_[i].W();
765 for (label i=0; i<nSpecie_; i++)
767 RR_[i][celli] = dc[i]*specieThermo_[i].W()/deltaT;
772 deltaTMin =
min(deltaTMin, 2*deltaT);