zhyang-dev/OpenFOAM-12
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

chemistryModel: Change of state variable to mass fraction

Browse Source 
The chemistry model now solves a system of mass fractions, temperature
and pressure, rather than a system of concentrations, temperature and
pressure.

The new form now accounts for the change in reaction rate associated
with thermal expansion. Thermal expansion (or contraction) can dilute
(or concentrate) the species concentrations, thereby reducing (or
increasing) the reaction rates. Previously it was not possible to
include this term because it was not computationally feasible to
evaluate it in a system in which the state variable was concentration.

The reaction rate interface has been simplified with respect to the
generation of derivatives. Reactions are defined as k*C, where k is the
reaction rate and C is the product of concentrations (raised to their
stoichiomentric coefficients and/or specified powers). Reaction rate
classes now provide two logical functions governing the derivatives of
k; i.e., ddT (derivative w.r.t. temperature) and ddc (derivative w.r.t.
concentration). Previously the reaction rate interface was closely
related to the form of third-body reactions, which made it inconvenient
to implement rates that were not very third-body-like.

It is now possible to verify the implementations of the jacobian methods
by comparison with finite-difference based evaluations of the rate
methods. This has been done and a number of bugs have been found and
fixed in the reaction rate classes.
This commit is contained in:
Will Bainbridge
2021-09-15 12:06:07 +01:00
parent 2d7a576d6b
commit 4f0dfc3bdf
57 changed files with 1595 additions and 1312 deletions
Download Patch File Download Diff File Expand all files Collapse all files

34
src/ODE/ODESystem/ODESystem.C
View File

@ -57,8 +57,15 @@ void Foam::ODESystem::check
jacobian(x, y, li, dfdx1, d2fdxdyAnalytic);
// Compare derivatives
Info<< "[derivatives] dfdx = " << dfdx0 << nl;
Info<< "[ jacobian] dfdx = " << dfdx1 << nl;
Info<< "[derivatives] dfdx = ( ";
forAll(dfdx0, i) { Info<< dfdx0[i] << ' '; }
Info<< ")" << nl;
Info<< "[ jacobian] dfdx = ( ";
forAll(dfdx1, i) { Info<< dfdx1[i] << ' '; }
Info<< ")" << nl;
Info<< "[ ratio] dfdx = ( ";
forAll(dfdx1, i) { Info<< dfdx1[i]/stabilise(dfdx0[i], rootVSmall) << ' '; }
Info<< ")" << nl;
// Construct a Jacobian using the finite differences and the derivatives
// method
@ -83,10 +90,25 @@ void Foam::ODESystem::check
for (label i = 0; i < nEqns(); ++ i)
{
Info<< "[derivatives] d2fdxdy[" << i << "] = "
<< UList<scalar>(d2fdxdyFiniteDifference[i], nEqns()) << nl;
Info<< "[ jacobian] d2fdxdy[" << i << "] = "
<< UList<scalar>(d2fdxdyAnalytic[i], nEqns()) << nl;
UList<scalar> FD(d2fdxdyFiniteDifference[i], nEqns());
Info<< "[derivatives] d2fdxdy[" << i << "] = ( ";
forAll(FD, i) { Info<< FD[i] << ' '; }
Info<< ")" << nl;
}
for (label i = 0; i < nEqns(); ++ i)
{
UList<scalar> A(d2fdxdyAnalytic[i], nEqns());
Info<< "[ jacobian] d2fdxdy[" << i << "] = ( ";
forAll(A, i) { Info<< A[i] << ' '; }
Info<< ")" << nl;
}
for (label i = 0; i < nEqns(); ++ i)
{
UList<scalar> FD(d2fdxdyFiniteDifference[i], nEqns());
UList<scalar> A(d2fdxdyAnalytic[i], nEqns());
Info<< "[ ratio] d2fdxdy[" << i << "] = ( ";
forAll(A, i) { Info<< A[i]/stabilise(FD[i], rootVSmall) << ' '; }
Info<< ")" << nl;
}
}

21
src/thermophysicalModels/chemistryModel/basicChemistryModel/basicChemistryModel.C
View File

@ -25,7 +25,25 @@ License
#include "basicChemistryModel.H"
/* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
/* * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * * */
namespace Foam
{
template<>
const char* NamedEnum<basicChemistryModel::jacobianType, 2>::names[] =
{
"fast",
"exact"
};
}
const Foam::NamedEnum
<
Foam::basicChemistryModel::jacobianType,
2
> Foam::basicChemistryModel::jacobianTypeNames_;
namespace Foam
{
@ -33,6 +51,7 @@ namespace Foam
defineRunTimeSelectionTable(basicChemistryModel, thermo);
}
// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
void Foam::basicChemistryModel::correct()

16
src/thermophysicalModels/chemistryModel/basicChemistryModel/basicChemistryModel.H
View File

@ -54,6 +54,22 @@ class basicChemistryModel
:
public IOdictionary
{
public:
// Public Enumerations
//- Enumeration for the type of Jacobian to be calculated by the
// chemistry model
enum class jacobianType
{
fast,
exact
};
//- Jacobian type names
static const NamedEnum<jacobianType, 2> jacobianTypeNames_;
protected:
// Protected data

481
src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.C
View File

@ -39,16 +39,22 @@ Foam::chemistryModel<ThermoType>::chemistryModel
basicChemistryModel(thermo),
ODESystem(),
log_(this->lookupOrDefault("log", false)),
Y_(this->thermo().composition().Y()),
jacobianType_
(
this->found("jacobian")
? jacobianTypeNames_.read(this->lookup("jacobian"))
: jacobianType::fast
),
Yvf_(this->thermo().composition().Y()),
mixture_(refCast<const multiComponentMixture<ThermoType>>(this->thermo())),
specieThermos_(mixture_.specieThermos()),
reactions_(mixture_.species(), specieThermos_, this->mesh(), *this),
nSpecie_(Y_.size()),
nReaction_(reactions_.size()),
Treact_(basicChemistryModel::template lookupOrDefault<scalar>("Treact", 0)),
nSpecie_(Yvf_.size()),
RR_(nSpecie_),
Y_(nSpecie_),
c_(nSpecie_),
dcdt_(nSpecie_),
YTpWork_(scalarField(nSpecie_ + 2)),
YTpYTpWork_(scalarSquareMatrix(nSpecie_ + 2)),
cTos_(nSpecie_, -1),
sToc_(nSpecie_),
mechRedPtr_
@ -81,7 +87,7 @@ Foam::chemistryModel<ThermoType>::chemistryModel
(
IOobject
(
"RR." + Y_[fieldi].name(),
"RR." + Yvf_[fieldi].name(),
this->mesh().time().timeName(),
this->mesh(),
IOobject::NO_READ,
@ -94,7 +100,7 @@ Foam::chemistryModel<ThermoType>::chemistryModel
}
Info<< "chemistryModel: Number of species = " << nSpecie_
<< " and reactions = " << nReaction_ << endl;
<< " and reactions = " << nReaction() << endl;
// When the mechanism reduction method is used, the 'active' flag for every
// species should be initialised (by default 'active' is true)
@ -102,11 +108,11 @@ Foam::chemistryModel<ThermoType>::chemistryModel
{
const basicSpecieMixture& composition = this->thermo().composition();
forAll(Y_, i)
forAll(Yvf_, i)
{
typeIOobject<volScalarField> header
(
Y_[i].name(),
Yvf_[i].name(),
this->mesh().time().timeName(),
this->mesh(),
IOobject::NO_READ
@ -137,101 +143,95 @@ Foam::chemistryModel<ThermoType>::~chemistryModel()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class ThermoType>
void Foam::chemistryModel<ThermoType>::omega
(
const scalar p,
const scalar T,
const scalarField& c, // Contains all species even when mechRed is active
const label li,
scalarField& dcdt
) const
{
if (mechRedActive_)
{
forAll(sToc_, si)
{
dcdt_[sToc_[si]] = 0;
}
forAll(reactions_, i)
{
if (!mechRed_.reactionDisabled(i))
{
const Reaction<ThermoType>& R = reactions_[i];
R.omega(p, T, c, li, dcdt_);
}
}
forAll(sToc_, si)
{
dcdt[si] = dcdt_[sToc_[si]];
}
}
else
{
dcdt = Zero;
forAll(reactions_, i)
{
const Reaction<ThermoType>& R = reactions_[i];
R.omega(p, T, c, li, dcdt);
}
}
}
template<class ThermoType>
void Foam::chemistryModel<ThermoType>::derivatives
(
const scalar time,
const scalarField& c,
const scalarField& YTp,
const label li,
scalarField& dcdt
scalarField& dYTpdt
) const
{
const scalar T = c[nSpecie_];
const scalar p = c[nSpecie_ + 1];
if (mechRedActive_)
{
forAll(sToc_, i)
{
c_[sToc_[i]] = max(c[i], 0);
Y_[sToc_[i]] = max(YTp[i], 0);
}
}
else
{
forAll(c_, i)
forAll(Y_, i)
{
c_[i] = max(c[i], 0);
Y_[i] = max(YTp[i], 0);
}
}
const scalar T = YTp[nSpecie_];
const scalar p = YTp[nSpecie_ + 1];
// Evaluate the mixture density
scalar rhoM = 0;
for (label i=0; i<Y_.size(); i++)
{
rhoM += Y_[i]/specieThermos_[i].rho(p, T);
}
rhoM = 1/rhoM;
// Evaluate the concentrations
for (label i=0; i<Y_.size(); i ++)
{
c_[i] = rhoM/specieThermos_[i].W()*Y_[i];
}
// Evaluate contributions from reactions
omega(p, T, c_, li, dcdt);
dYTpdt = Zero;
forAll(reactions_, ri)
{
if (!mechRed_.reactionDisabled(ri))
{
reactions_[ri].dNdtByV
(
p,
T,
c_,
li,
dYTpdt,
mechRedActive_,
cTos_,
0
);
}
}
// Reactions return dNdtByV, so we need to convert the result to dYdt
for (label i=0; i<nSpecie_; i++)
{
const scalar WiByrhoM = specieThermos_[sToc(i)].W()/rhoM;
scalar& dYidt = dYTpdt[i];
dYidt *= WiByrhoM;
}
// Evaluate the effect on the thermodynamic system ...
// c*Cp
scalar ccp = 0;
for (label i=0; i<c_.size(); i++)
// Evaluate the mixture Cp
scalar CpM = 0;
for (label i=0; i<Y_.size(); i++)
{
ccp += c_[i]*specieThermos_[i].cp(p, T);
CpM += Y_[i]*specieThermos_[i].Cp(p, T);
}
// dT/dt
scalar& dTdt = dcdt[nSpecie_];
dTdt = 0;
scalar& dTdt = dYTpdt[nSpecie_];
for (label i=0; i<nSpecie_; i++)
{
dTdt -= dcdt[i]*specieThermos_[sToc(i)].ha(p, T);
dTdt -= dYTpdt[i]*specieThermos_[sToc(i)].Ha(p, T);
}
dTdt /= ccp;
dTdt /= CpM;
// dp/dt = 0 (pressure is assumed constant)
dcdt[nSpecie_ + 1] = 0;
scalar& dpdt = dYTpdt[nSpecie_ + 1];
dpdt = 0;
}
@ -239,125 +239,217 @@ template<class ThermoType>
void Foam::chemistryModel<ThermoType>::jacobian
(
const scalar t,
const scalarField& c,
const scalarField& YTp,
const label li,
scalarField& dcdt,
scalarField& dYTpdt,
scalarSquareMatrix& J
) const
{
// If the mechanism reduction is active, the computed Jacobian
// is compact (size of the reduced set of species)
// but according to the information of the complete set
// (i.e. for the third-body efficiencies)
if (mechRedActive_)
{
forAll(sToc_, i)
{
c_[sToc_[i]] = max(c[i], 0);
Y_[sToc_[i]] = max(YTp[i], 0);
}
}
else
{
forAll(c_, i)
{
c_[i] = max(c[i], 0);
Y_[i] = max(YTp[i], 0);
}
}
const scalar T = c[nSpecie_];
const scalar p = c[nSpecie_ + 1];
const scalar T = YTp[nSpecie_];
const scalar p = YTp[nSpecie_ + 1];
dcdt = Zero;
J = Zero;
// Evaluate the specific volumes and mixture density
scalarField& v = YTpWork_[0];
for (label i=0; i<Y_.size(); i++)
{
v[i] = 1/specieThermos_[i].rho(p, T);
}
scalar rhoM = 0;
for (label i=0; i<Y_.size(); i++)
{
rhoM += Y_[i]*v[i];
}
rhoM = 1/rhoM;
// Evaluate the concentrations
for (label i=0; i<Y_.size(); i ++)
{
c_[i] = rhoM/specieThermos_[i].W()*Y_[i];
}
// Evaluate the derivatives of concentration w.r.t. mass fraction
scalarSquareMatrix& dcdY = YTpYTpWork_[0];
for (label i=0; i<nSpecie_; i++)
{
const scalar rhoMByWi = rhoM/specieThermos_[sToc(i)].W();
switch (jacobianType_)
{
case jacobianType::fast:
{
dcdY(i, i) = rhoMByWi;
}
break;
case jacobianType::exact:
for (label j=0; j<nSpecie_; j++)
{
dcdY(i, j) =
rhoMByWi*((i == j) - rhoM*v[sToc(j)]*Y_[sToc(i)]);
}
break;
}
}
// Evaluate the mixture thermal expansion coefficient
scalar alphavM = 0;
for (label i=0; i<Y_.size(); i++)
{
alphavM += Y_[i]*rhoM*v[i]*specieThermos_[i].alphav(p, T);
}
// Evaluate contributions from reactions
dYTpdt = Zero;
scalarSquareMatrix& ddNdtByVdcTp = YTpYTpWork_[1];
for (label i=0; i<nSpecie_ + 2; i++)
{
for (label j=0; j<nSpecie_ + 2; j++)
{
ddNdtByVdcTp[i][j] = 0;
}
}
forAll(reactions_, ri)
{
if (!mechRed_.reactionDisabled(ri))
{
const Reaction<ThermoType>& R = reactions_[ri];
scalar omegaI, kfwd, kbwd;
R.dwdc
reactions_[ri].ddNdtByVdcTp
(
p,
T,
c_,
li,
J,
dcdt,
omegaI,
kfwd,
kbwd,
mechRedActive_,
cTos_
);
R.dwdT
(
p,
T,
c_,
li,
omegaI,
kfwd,
kbwd,
J,
dYTpdt,
ddNdtByVdcTp,
mechRedActive_,
cTos_,
nSpecie_
0,
nSpecie_,
YTpWork_[1],
YTpWork_[2]
);
}
}
// Reactions return dNdtByV, so we need to convert the result to dYdt
for (label i=0; i<nSpecie_; i++)
{
const scalar WiByrhoM = specieThermos_[sToc(i)].W()/rhoM;
scalar& dYidt = dYTpdt[i];
dYidt *= WiByrhoM;
for (label j=0; j<nSpecie_; j++)
{
scalar ddNidtByVdYj = 0;
switch (jacobianType_)
{
case jacobianType::fast:
{
const scalar ddNidtByVdcj = ddNdtByVdcTp(i, j);
ddNidtByVdYj = ddNidtByVdcj*dcdY(j, j);
}
break;
case jacobianType::exact:
for (label k=0; k<nSpecie_; k++)
{
const scalar ddNidtByVdck = ddNdtByVdcTp(i, k);
ddNidtByVdYj += ddNidtByVdck*dcdY(k, j);
}
break;
}
scalar& ddYidtdYj = J(i, j);
ddYidtdYj = WiByrhoM*ddNidtByVdYj + rhoM*v[sToc(j)]*dYidt;
}
scalar ddNidtByVdT = ddNdtByVdcTp(i, nSpecie_);
for (label j=0; j<nSpecie_; j++)
{
const scalar ddNidtByVdcj = ddNdtByVdcTp(i, j);
ddNidtByVdT -= ddNidtByVdcj*c_[sToc(j)]*alphavM;
}
scalar& ddYidtdT = J(i, nSpecie_);
ddYidtdT = WiByrhoM*ddNidtByVdT + alphavM*dYidt;
scalar& ddYidtdp = J(i, nSpecie_ + 1);
ddYidtdp = 0;
}
// Evaluate the effect on the thermodynamic system ...
// c*Cp
scalar ccp = 0, dccpdT = 0;
forAll(c_, i)
// Evaluate the mixture Cp and its derivative
scalarField& Cp = YTpWork_[3];
scalar CpM = 0, dCpMdT = 0;
for (label i=0; i<Y_.size(); i++)
{
ccp += c_[i]*specieThermos_[i].cp(p, T);
dccpdT += c_[i]*specieThermos_[i].dcpdT(p, T);
Cp[i] = specieThermos_[i].Cp(p, T);
CpM += Y_[i]*Cp[i];
dCpMdT += Y_[i]*specieThermos_[i].dCpdT(p, T);
}
// dT/dt
scalar& dTdt = dcdt[nSpecie_];
scalarField& Ha = YTpWork_[4];
scalar& dTdt = dYTpdt[nSpecie_];
for (label i=0; i<nSpecie_; i++)
{
dTdt -= dcdt[i]*specieThermos_[sToc(i)].ha(p, T);
Ha[sToc(i)] = specieThermos_[sToc(i)].Ha(p, T);
dTdt -= dYTpdt[i]*Ha[sToc(i)];
}
dTdt /= ccp;
dTdt /= CpM;
// dp/dt = 0 (pressure is assumed constant)
scalar& dpdt = dYTpdt[nSpecie_ + 1];
dpdt = 0;
// d(dTdt)/dc
for (label i = 0; i < nSpecie_; i++)
// d(dTdt)/dY
for (label i=0; i<nSpecie_; i++)
{
scalar& d2Tdtdci = J(nSpecie_, i);
for (label j = 0; j < nSpecie_; j++)
scalar& ddTdtdYi = J(nSpecie_, i);
ddTdtdYi = 0;
for (label j=0; j<nSpecie_; j++)
{
const scalar d2cjdtdci = J(j, i);
d2Tdtdci -= d2cjdtdci*specieThermos_[sToc(j)].ha(p, T);
const scalar ddYjdtdYi = J(j, i);
ddTdtdYi -= ddYjdtdYi*Ha[sToc(j)];
}
d2Tdtdci -= specieThermos_[sToc(i)].cp(p, T)*dTdt;
d2Tdtdci /= ccp;
ddTdtdYi -= Cp[sToc(i)]*dTdt;
ddTdtdYi /= CpM;
}
// d(dTdt)/dT
scalar& d2TdtdT = J(nSpecie_, nSpecie_);
for (label i = 0; i < nSpecie_; i++)
scalar& ddTdtdT = J(nSpecie_, nSpecie_);
ddTdtdT = 0;
for (label i=0; i<nSpecie_; i++)
{
const scalar d2cidtdT = J(i, nSpecie_);
const label si = sToc(i);
d2TdtdT -=
dcdt[i]*specieThermos_[si].cp(p, T)
+ d2cidtdT*specieThermos_[si].ha(p, T);
const scalar dYidt = dYTpdt[i];
const scalar ddYidtdT = J(i, nSpecie_);
ddTdtdT -= dYidt*Cp[sToc(i)] + ddYidtdT*Ha[sToc(i)];
}
d2TdtdT -= dTdt*dccpdT;
d2TdtdT /= ccp;
ddTdtdT -= dTdt*dCpMdT;
ddTdtdT /= CpM;
// d(dpdt)/dc = 0 (pressure is assumed constant)
// d(dTdt)/dp = 0 (pressure is assumed constant)
scalar& ddTdtdp = J(nSpecie_, nSpecie_ + 1);
ddTdtdp = 0;
// d(dpdt)/dT = 0 (pressure is assumed constant)
// d(dpdt)/dYiTp = 0 (pressure is assumed constant)
for (label i=0; i<nSpecie_ + 2; i++)
{
scalar& ddpdtdYiTp = J(nSpecie_ + 1, i);
ddpdtdYiTp = 0;
}
}
@ -395,7 +487,7 @@ Foam::chemistryModel<ThermoType>::tc() const
for (label i=0; i<nSpecie_; i++)
{
c_[i] = rhoi*Y_[i][celli]/specieThermos_[i].W();
c_[i] = rhoi*Yvf_[i][celli]/specieThermos_[i].W();
}
// A reaction's rate scale is calculated as it's molar
@ -470,7 +562,7 @@ Foam::chemistryModel<ThermoType>::Qdot() const
scalarField& Qdot = tQdot.ref();
forAll(Y_, i)
forAll(Yvf_, i)
{
forAll(Qdot, celli)
{
@ -522,7 +614,7 @@ Foam::chemistryModel<ThermoType>::calculateRR
for (label i=0; i<nSpecie_; i++)
{
const scalar Yi = Y_[i][celli];
const scalar Yi = Yvf_[i][celli];
c_[i] = rhoi*Yi/specieThermos_[i].W();
}
@ -566,6 +658,8 @@ void Foam::chemistryModel<ThermoType>::calculate()
const scalarField& T = this->thermo().T();
const scalarField& p = this->thermo().p();
scalarField& dNdtByV = YTpWork_[0];
reactionEvaluationScope scope(*this);
forAll(rho, celli)
@ -576,15 +670,33 @@ void Foam::chemistryModel<ThermoType>::calculate()
for (label i=0; i<nSpecie_; i++)
{
const scalar Yi = Y_[i][celli];
const scalar Yi = Yvf_[i][celli];
c_[i] = rhoi*Yi/specieThermos_[i].W();
}
omega(pi, Ti, c_, celli, dcdt_);
dNdtByV = Zero;
for (label i=0; i<nSpecie_; i++)
forAll(reactions_, ri)
{
RR_[i][celli] = dcdt_[i]*specieThermos_[i].W();
if (!mechRed_.reactionDisabled(ri))
{
reactions_[ri].dNdtByV
(
pi,
Ti,
c_,
celli,
dNdtByV,
mechRedActive_,
cTos_,
0
);
}
}
for (label i=0; i<mechRed_.nActiveSpecies(); i++)
{
RR_[sToc(i)][celli] = dNdtByV[i]*specieThermos_[sToc(i)].W();
}
}
}
@ -615,35 +727,41 @@ Foam::scalar Foam::chemistryModel<ThermoType>::solve
return deltaTMin;
}
tmp<volScalarField> trho(this->thermo().rho0());
const scalarField& rho = trho();
tmp<volScalarField> trhovf(this->thermo().rho());
const volScalarField& rhovf = trhovf();
tmp<volScalarField> trho0vf(this->thermo().rho0());
const volScalarField& rho0vf = trho0vf();
const scalarField& T = this->thermo().T().oldTime();
const scalarField& p = this->thermo().p().oldTime();
const volScalarField& T0vf = this->thermo().T().oldTime();
const volScalarField& p0vf = this->thermo().p().oldTime();
reactionEvaluationScope scope(*this);
scalarField c0(nSpecie_);
scalarField Y0(nSpecie_);
// Composition vector (Yi, T, p, deltaT)
scalarField phiq(nEqns() + 1);
scalarField Rphiq(nEqns() + 1);
forAll(rho, celli)
forAll(rho0vf, celli)
{
const scalar rhoi = rho[celli];
scalar pi = p[celli];
scalar Ti = T[celli];
const scalar rho = rhovf[celli];
const scalar rho0 = rho0vf[celli];
scalar p = p0vf[celli];
scalar T = T0vf[celli];
for (label i=0; i<nSpecie_; i++)
{
c_[i] =
rhoi*Y_[i].oldTime()[celli]/specieThermos_[i].W();
c0[i] = c_[i];
phiq[i] = Y_[i].oldTime()[celli];
Y_[i] = Y0[i] = Yvf_[i].oldTime()[celli];
}
phiq[nSpecie()] = Ti;
phiq[nSpecie() + 1] = pi;
for (label i=0; i<nSpecie_; i++)
{
phiq[i] = Yvf_[i].oldTime()[celli];
}
phiq[nSpecie()] = T;
phiq[nSpecie() + 1] = p;
phiq[nSpecie() + 2] = deltaT[celli];
// Initialise time progress
@ -660,8 +778,10 @@ Foam::scalar Foam::chemistryModel<ThermoType>::solve
// Retrieved solution stored in Rphiq
for (label i=0; i<nSpecie(); i++)
{
c_[i] = rhoi*Rphiq[i]/specieThermos_[i].W();
Y_[i] = Rphiq[i];
}
T = Rphiq[nSpecie()];
p = Rphiq[nSpecie() + 1];
}
// This position is reached when tabulation is not used OR
// if the solution is not retrieved.
@ -671,17 +791,29 @@ Foam::scalar Foam::chemistryModel<ThermoType>::solve
{
if (mechRedActive_)
{
for (label i=0; i<nSpecie_; i++)
{
const scalar Yi = Yvf_[i][celli];
c_[i] = rho0*Yi/specieThermos_[i].W();
}
// Reduce mechanism change the number of species (only active)
mechRed_.reduceMechanism
(
pi,
Ti,
p,
T,
c_,
sc_,
cTos_,
sToc_,
celli
);
sY_.setSize(nSpecie_);
for (label i=0; i<nSpecie_; i++)
{
sY_[i] = specieThermos_[sToc(i)].W()*sc_[i]/rho0;
}
}
if (log_)
@ -699,9 +831,9 @@ Foam::scalar Foam::chemistryModel<ThermoType>::solve
// Solve the reduced set of ODE
solve
(
pi,
Ti,
sc_,
p,
T,
sY_,
celli,
dt,
deltaTChem_[celli]
@ -709,13 +841,12 @@ Foam::scalar Foam::chemistryModel<ThermoType>::solve
for (label i=0; i<mechRed_.nActiveSpecies(); i++)
{
c_[sToc_[i]] = sc_[i];
Y_[sToc_[i]] = sY_[i];
}
}
else
{
solve
(pi, Ti, c_, celli, dt, deltaTChem_[celli]);
solve(p, T, Y_, celli, dt, deltaTChem_[celli]);
}
timeLeft -= dt;
}
@ -729,15 +860,15 @@ Foam::scalar Foam::chemistryModel<ThermoType>::solve
// the stored points (either expand or add)
if (tabulation_.tabulates())
{
forAll(c_, i)
forAll(Y_, i)
{
Rphiq[i] = c_[i]/rhoi*specieThermos_[i].W();
Rphiq[i] = Y_[i];
}
Rphiq[Rphiq.size()-3] = Ti;
Rphiq[Rphiq.size()-2] = pi;
Rphiq[Rphiq.size()-3] = T;
Rphiq[Rphiq.size()-2] = p;
Rphiq[Rphiq.size()-1] = deltaT[celli];
tabulation_.add(phiq, Rphiq, celli, rhoi, deltaT[celli]);
tabulation_.add(phiq, Rphiq, celli, rho0, deltaT[celli]);
}
// When operations are done and if mechanism reduction is active,
@ -745,19 +876,17 @@ Foam::scalar Foam::chemistryModel<ThermoType>::solve
// to the total number of species (stored in the mechRed object)
if (mechRedActive_)
{
nSpecie_ = mechRed_.nSpecie();
setNSpecie(mechRed_.nSpecie());
}
deltaTMin = min(deltaTChem_[celli], deltaTMin);
deltaTChem_[celli] =
min(deltaTChem_[celli], deltaTChemMax_);
deltaTMin = min(deltaTChem_[celli], deltaTMin);
deltaTChem_[celli] = min(deltaTChem_[celli], deltaTChemMax_);
}
// Set the RR vector (used in the solver)
for (label i=0; i<nSpecie_; i++)
{
RR_[i][celli] =
(c_[i] - c0[i])*specieThermos_[i].W()/deltaT[celli];
RR_[i][celli] = (Y_[i]*rho - Y0[i]*rho0)/deltaT[celli];
}
}

54
src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.H
View File

@ -130,8 +130,11 @@ class chemistryModel
//- Switch to select performance logging
Switch log_;
//- Type of the Jacobian to be calculated
const jacobianType jacobianType_;
//- Reference to the field of specie mass fractions
const PtrList<volScalarField>& Y_;
const PtrList<volScalarField>& Yvf_;
//- Reference to the multi component mixture
const multiComponentMixture<ThermoType>& mixture_;
@ -145,24 +148,27 @@ class chemistryModel
//- Number of species
label nSpecie_;
//- Number of reactions
label nReaction_;
//- Temperature below which the reaction rates are assumed 0
scalar Treact_;
//- List of reaction rate per specie [kg/m^3/s]
PtrList<volScalarField::Internal> RR_;
//- Temporary mass fraction field
mutable scalarField Y_;
//- Temporary simplified mechanism mass fraction field
DynamicField<scalar> sY_;
//- Temporary concentration field
mutable scalarField c_;
//- Temporary rate-of-change of concentration field
mutable scalarField dcdt_;
//- Temporary simplified mechanism concentration field
DynamicField<scalar> sc_;
//- Specie-temperature-pressure workspace fields
mutable FixedList<scalarField, 5> YTpWork_;
//- Specie-temperature-pressure workspace matrices
mutable FixedList<scalarSquareMatrix, 2> YTpYTpWork_;
//- Temporary map from complete to simplified concentration fields
// c -> sc
List<label> cTos_;
@ -241,22 +247,6 @@ public:
//- The number of reactions
virtual inline label nReaction() const;
//- Temperature below which the reaction rates are assumed 0
inline scalar Treact() const;
//- Temperature below which the reaction rates are assumed 0
inline scalar& Treact();
//- dc/dt = omega, rate of change in concentration, for each species
virtual void omega
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcdt
) const;
//- Calculates the reaction rates
virtual void calculate();
@ -306,17 +296,17 @@ public:
virtual void derivatives
(
const scalar t,
const scalarField& c,
const scalarField& YTp,
const label li,
scalarField& dcdt
scalarField& dYTpdt
) const;
virtual void jacobian
(
const scalar t,
const scalarField& c,
const scalarField& YTp,
const label li,
scalarField& dcdt,
scalarField& dYTpdt,
scalarSquareMatrix& J
) const;
@ -324,7 +314,7 @@ public:
(
scalar& p,
scalar& T,
scalarField& c,
scalarField& Y,
const label li,
scalar& deltaT,
scalar& subDeltaT
@ -332,7 +322,7 @@ public:
//- Return a reference to the list of mass fraction fields
inline const PtrList<volScalarField>& Y();
inline const PtrList<volScalarField>& Y() const;
// Mechanism reduction access functions

20
src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModelI.H
View File

@ -91,21 +91,7 @@ inline Foam::label Foam::chemistryModel<ThermoType>::nSpecie() const
template<class ThermoType>
inline Foam::label Foam::chemistryModel<ThermoType>::nReaction() const
{
return nReaction_;
}
template<class ThermoType>
inline Foam::scalar Foam::chemistryModel<ThermoType>::Treact() const
{
return Treact_;
}
template<class ThermoType>
inline Foam::scalar& Foam::chemistryModel<ThermoType>::Treact()
{
return Treact_;
return reactions_.size();
}
@ -126,9 +112,9 @@ Foam::chemistryModel<ThermoType>::RR(const label i)
template<class ThermoType>
inline const Foam::PtrList<Foam::volScalarField>&
Foam::chemistryModel<ThermoType>::Y()
Foam::chemistryModel<ThermoType>::Y() const
{
return this->Y_;
return this->Yvf_;
}

70
src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/ISAT/ISAT.C
View File

@ -32,7 +32,7 @@ template<class ThermoType>
Foam::chemistryTabulationMethods::ISAT<ThermoType>::ISAT
(
const dictionary& chemistryProperties,
chemistryModel<ThermoType>& chemistry
const chemistryModel<ThermoType>& chemistry
)
:
chemistryTabulationMethod<ThermoType>
@ -352,16 +352,17 @@ void Foam::chemistryTabulationMethods::ISAT<ThermoType>::computeA
const scalar dt
)
{
const label speciesNumber = chemistry_.nSpecie();
scalarField Rcq(chemistry_.nEqns() + 1);
for (label i=0; i<speciesNumber; i++)
const label nSpecie = chemistry_.nSpecie();
scalarField Rphiqs(chemistry_.nEqns() + 1);
for (label i=0; i<nSpecie; i++)
{
const label si = chemistry_.sToc(i);
Rcq[i] = rhoi*Rphiq[si]/chemistry_.specieThermos()[si].W();
Rphiqs[i] = Rphiq[si];
}
Rcq[speciesNumber] = Rphiq[Rphiq.size() - 3];
Rcq[speciesNumber + 1] = Rphiq[Rphiq.size() - 2];
Rcq[speciesNumber + 2] = Rphiq[Rphiq.size() - 1];
Rphiqs[nSpecie] = Rphiq[Rphiq.size() - 3];
Rphiqs[nSpecie + 1] = Rphiq[Rphiq.size() - 2];
Rphiqs[nSpecie + 2] = Rphiq[Rphiq.size() - 1];
// Aaa is computed implicitly,
// A is given by A = C(psi0, t0+dt), where C is obtained through solving
@ -372,62 +373,29 @@ void Foam::chemistryTabulationMethods::ISAT<ThermoType>::computeA
// C(psi0,t0+dt)*(I-dt*J(psi(t0+dt))) = C(psi0, t0)
// A = C(psi0,t0)/(I-dt*J(psi(t0+dt)))
// where C(psi0,t0) = I
scalarField dcdt(speciesNumber + 2, Zero);
chemistry_.jacobian(runTime_.value(), Rcq, li, dcdt, A);
scalarField dYTpdt(nSpecie + 2, Zero);
chemistry_.jacobian(runTime_.value(), Rphiqs, li, dYTpdt, A);
// The jacobian is computed according to the molar concentration
// the following conversion allows the code to use A with mass fraction
for (label i=0; i<speciesNumber; i++)
// Inverse of I - dt*J(psi(t0 + dt))
for (label i=0; i<nSpecie + 2; i++)
{
const label si = chemistry_.sToc(i);
for (label j=0; j<speciesNumber; j++)
for (label j=0; j<nSpecie + 2; j++)
{
const label sj = chemistry_.sToc(j);
A(i, j) *=
-dt*chemistry_.specieThermos()[si].W()
/chemistry_.specieThermos()[sj].W();
A(i, j) *= -dt;
}
A(i, i) += 1;
// Columns for pressure and temperature
A(i, speciesNumber) *=
-dt*chemistry_.specieThermos()[si].W()/rhoi;
A(i, speciesNumber + 1) *=
-dt*chemistry_.specieThermos()[si].W()/rhoi;
}
// For the temperature and pressure lines, ddc(dTdt)
// should be converted in ddY(dTdt)
for (label i=0; i<speciesNumber; i++)
{
const label si = chemistry_.sToc(i);
A(speciesNumber, i) *=
-dt*rhoi/chemistry_.specieThermos()[si].W();
A(speciesNumber + 1, i) *=
-dt*rhoi/chemistry_.specieThermos()[si].W();
}
A(speciesNumber, speciesNumber) = -dt*A(speciesNumber, speciesNumber) + 1;
A(speciesNumber + 1, speciesNumber + 1) =
-dt*A(speciesNumber + 1, speciesNumber + 1) + 1;
A(speciesNumber + 2, speciesNumber + 2) = 1;
// Inverse of (I-dt*J(psi(t0+dt)))
A(nSpecie + 2, nSpecie + 2) = 1;
LUscalarMatrix LUA(A);
LUA.inv(A);
// After inversion, lines of p and T are set to 0 except diagonal. This
// avoid skewness of the ellipsoid of accuracy and potential issues in the
// binary tree.
for (label i=0; i<speciesNumber; i++)
for (label i=0; i<nSpecie; i++)
{
A(speciesNumber, i) = 0;
A(speciesNumber + 1, i) = 0;
A(nSpecie, i) = 0;
A(nSpecie + 1, i) = 0;
}
}

6
src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/ISAT/ISAT.H
View File

@ -63,7 +63,7 @@ class ISAT
const dictionary coeffsDict_;
chemistryModel<ThermoType>& chemistry_;
const chemistryModel<ThermoType>& chemistry_;
//- Switch to select performance logging
Switch log_;
@ -214,7 +214,7 @@ public:
ISAT
(
const dictionary& chemistryProperties,
chemistryModel<ThermoType>& chemistry
const chemistryModel<ThermoType>& chemistry
);
//- Disallow default bitwise copy construction
@ -233,7 +233,7 @@ public:
return true;
}
chemistryModel<ThermoType>& chemistry()
const chemistryModel<ThermoType>& chemistry()
{
return chemistry_;
}

2
src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.H
View File

@ -276,7 +276,7 @@ public:
// Member Functions
//- Access to the chemistryModel
inline chemistryModel<ThermoType>& chemistry()
inline const chemistryModel<ThermoType>& chemistry()
{
return table_.chemistry();
}

2
src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethod.C
View File

@ -32,7 +32,7 @@ template<class ThermoType>
Foam::chemistryTabulationMethod<ThermoType>::chemistryTabulationMethod
(
const dictionary& dict,
chemistryModel<ThermoType>& chemistry
const chemistryModel<ThermoType>& chemistry
)
{}

6
src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethod.H
View File

@ -71,7 +71,7 @@ public:
dictionary,
(
const dictionary& dict,
chemistryModel<ThermoType>& chemistry
const chemistryModel<ThermoType>& chemistry
),
(dict, chemistry)
);
@ -83,7 +83,7 @@ public:
chemistryTabulationMethod
(
const dictionary& dict,
chemistryModel<ThermoType>& chemistry
const chemistryModel<ThermoType>& chemistry
);
@ -92,7 +92,7 @@ public:
static autoPtr<chemistryTabulationMethod> New
(
const IOdictionary& dict,
chemistryModel<ThermoType>& chemistry
const chemistryModel<ThermoType>& chemistry
);

2
src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethodNew.C
View File

@ -34,7 +34,7 @@ Foam::autoPtr<Foam::chemistryTabulationMethod<ThermoType>>
Foam::chemistryTabulationMethod<ThermoType>::New
(
const IOdictionary& dict,
chemistryModel<ThermoType>& chemistry
const chemistryModel<ThermoType>& chemistry
)
{
if (dict.found("tabulation"))

2
src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/noChemistryTabulation/noChemistryTabulation.C
View File

@ -31,7 +31,7 @@ template<class ThermoType>
Foam::chemistryTabulationMethods::none<ThermoType>::none
(
const dictionary& chemistryProperties,
chemistryModel<ThermoType>& chemistry
const chemistryModel<ThermoType>& chemistry
)
:
chemistryTabulationMethod<ThermoType>

2
src/thermophysicalModels/chemistryModel/chemistryModel/tabulation/noChemistryTabulation/noChemistryTabulation.H
View File

@ -65,7 +65,7 @@ public:
none
(
const dictionary& chemistryProperties,
chemistryModel<ThermoType>& chemistry
const chemistryModel<ThermoType>& chemistry
);
//- Disallow default bitwise copy construction

4
src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C
View File

@ -77,8 +77,8 @@ void Foam::ode<ChemistryModel>::solve
if (debug)
{
scalarField dcTp(this->nEqns(), rootSmall);
dcTp[nSpecie] = 300*rootSmall;
dcTp[nSpecie+1] = 1e5*rootSmall;
dcTp[nSpecie] = T*rootSmall;
dcTp[nSpecie+1] = p*rootSmall;
this->check(0, cTp_, dcTp, li);
}

19
src/thermophysicalModels/chemistryModel/reactions/fluxLimitedLangmuirHinshelwood/fluxLimitedLangmuirHinshelwoodReactionRate.H
View File

@ -149,28 +149,15 @@ public:
const label li
) const;
//- Third-body efficiencies (beta = 1-alpha)
// non-empty only for third-body reactions
// with enhanced molecularity (alpha != 1)
inline const List<Tuple2<label, scalar>>& beta() const;
inline bool hasDdc() const;
//- Species concentration derivative of the pressure dependent term
inline void dcidc
inline void ddc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcidc
) const;
//- Temperature derivative of the pressure dependent term
inline scalar dcidT
(
const scalar p,
const scalar T,
const scalarField& c,
const label li
scalarField& ddc
) const;
//- Write to stream

21
src/thermophysicalModels/chemistryModel/reactions/fluxLimitedLangmuirHinshelwood/fluxLimitedLangmuirHinshelwoodReactionRateI.H
View File

@ -339,33 +339,22 @@ inline Foam::scalar Foam::fluxLimitedLangmuirHinshelwoodReactionRate::ddT
}
inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
Foam::fluxLimitedLangmuirHinshelwoodReactionRate::beta() const
inline bool Foam::fluxLimitedLangmuirHinshelwoodReactionRate::hasDdc() const
{
return NullObjectRef<List<Tuple2<label, scalar>>>();
return false;
}
inline void Foam::fluxLimitedLangmuirHinshelwoodReactionRate::dcidc
inline void Foam::fluxLimitedLangmuirHinshelwoodReactionRate::ddc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcidc
) const
{}
inline Foam::scalar Foam::fluxLimitedLangmuirHinshelwoodReactionRate::dcidT
(
const scalar p,
const scalar T,
const scalarField& c,
const label li
scalarField& ddc
) const
{
return 0;
ddc = 0;
}

20
src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.C
View File

@ -169,37 +169,39 @@ Foam::scalar Foam::IrreversibleReaction<ReactionThermo, ReactionRate>::dkrdT
template<class ReactionThermo, class ReactionRate>
const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
Foam::IrreversibleReaction<ReactionThermo, ReactionRate>::beta() const
bool Foam::IrreversibleReaction<ReactionThermo, ReactionRate>::hasDkdc() const
{
return k_.beta();
return k_.hasDdc();
}
template<class ReactionThermo, class ReactionRate>
void Foam::IrreversibleReaction<ReactionThermo, ReactionRate>::dcidc
void Foam::IrreversibleReaction<ReactionThermo, ReactionRate>::dkfdc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcidc
scalarField& dkfdc
) const
{
k_.dcidc(p, T, c, li, dcidc);
k_.ddc(p, T, c, li, dkfdc);
}
template<class ReactionThermo, class ReactionRate>
Foam::scalar Foam::IrreversibleReaction<ReactionThermo, ReactionRate>::dcidT
void Foam::IrreversibleReaction<ReactionThermo, ReactionRate>::dkrdc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li
const label li,
const scalarField& dkfdc,
const scalar kr,
scalarField& dkrdc
) const
{
return k_.dcidT(p, T, c, li);
dkrdc = 0;
}

28
src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.H
View File

@ -181,8 +181,7 @@ public:
const label li
) const;
//- Temperature derivative of reverse rate
// Returns 0
//- Temperature derivative of reverse rate. Returns zero.
virtual scalar dkrdT
(
const scalar p,
@ -193,35 +192,34 @@ public:
const scalar kr
) const;
//- Third-body efficiencies (beta = 1-alpha)
// non-empty only for third-body reactions
// with enhanced molecularity (alpha != 1)
virtual const List<Tuple2<label, scalar>>& beta() const;
//- Does this reaction have concentration-dependent rate constants?
virtual bool hasDkdc() const;
//- Species concentration derivative of the pressure dependent term
// By default this value is 1 as it multiplies the third-body term
virtual void dcidc
//- Concentration derivative of forward rate
void dkfdc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcidc
scalarField& dkfdc
) const;
//- Temperature derivative of the pressure dependent term
// By default this value is 0 since ddT of molecularity is approx.0
virtual scalar dcidT
//- Concentration derivative of reverse rate. Returns zero.
void dkrdc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li
const label li,
const scalarField& dkfdc,
const scalar kr,
scalarField& dkrdc
) const;
//- Write
virtual void write(Ostream&) const; // Private Member Functions
virtual void write(Ostream&) const;
// Member Operators

24
src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.C
View File

@ -181,40 +181,40 @@ Foam::NonEquilibriumReversibleReaction<ReactionThermo, ReactionRate>::dkrdT
template<class ReactionThermo, class ReactionRate>
const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
Foam::NonEquilibriumReversibleReaction<ReactionThermo, ReactionRate>::
beta() const
bool Foam::NonEquilibriumReversibleReaction<ReactionThermo, ReactionRate>::
hasDkdc() const
{
return fk_.beta();
return fk_.hasDdc() || rk_.hasDdc();
}
template<class ReactionThermo, class ReactionRate>
void
Foam::NonEquilibriumReversibleReaction<ReactionThermo, ReactionRate>::dcidc
void Foam::NonEquilibriumReversibleReaction<ReactionThermo, ReactionRate>::dkfdc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcidc
scalarField& dkfdc
) const
{
fk_.dcidc(p, T, c, li, dcidc);
fk_.ddc(p, T, c, li, dkfdc);
}
template<class ReactionThermo, class ReactionRate>
Foam::scalar
Foam::NonEquilibriumReversibleReaction<ReactionThermo, ReactionRate>::dcidT
void Foam::NonEquilibriumReversibleReaction<ReactionThermo, ReactionRate>::dkrdc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li
const label li,
const scalarField& dkfdc,
const scalar kr,
scalarField& dkrdc
) const
{
return fk_.dcidT(p, T, c, li);
rk_.ddc(p, T, c, li, dkrdc);
}

25
src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.H
View File

@ -163,7 +163,7 @@ public:
const label li
) const;
//- Reverse rate constant from the given formard rate constant
//- Reverse rate constant from the given forward rate constant
virtual scalar kr
(
const scalar kfwd,
@ -207,30 +207,29 @@ public:
const scalar kr
) const;
//- Third-body efficiencies (beta = 1-alpha)
// non-empty only for third-body reactions
// with enhanced molecularity (alpha != 1)
virtual const List<Tuple2<label, scalar>>& beta() const;
//- Does this reaction have concentration-dependent rate constants?
virtual bool hasDkdc() const;
//- Species concentration derivative of the pressure dependent term
// By default this value is 1 as it multiplies the third-body term
virtual void dcidc
//- Concentration derivative of forward rate
void dkfdc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcidc
scalarField& dkfdc
) const;
//- Temperature derivative of the pressure dependent term
// By default this value is 0 since ddT of molecularity is approx.0
virtual scalar dcidT
//- Concentration derivative of reverse rate
void dkrdc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li
const label li,
const scalarField& dkfdc,
const scalar kr,
scalarField& dkrdc
) const;

354
src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C
View File

@ -286,30 +286,30 @@ void Foam::Reaction<ReactionThermo>::write(Ostream& os) const
template<class ReactionThermo>
void Foam::Reaction<ReactionThermo>::omega
void Foam::Reaction<ReactionThermo>::C
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcdt
scalar& Cf,
scalar& Cr
) const
{
scalar omegaf, omegar;
const scalar omegaI = omega(p, T, c, li, omegaf, omegar);
Cf = Cr = 1;
forAll(lhs(), i)
{
const label si = lhs()[i].index;
const scalar sl = lhs()[i].stoichCoeff;
dcdt[si] -= sl*omegaI;
const specieExponent& el = lhs()[i].exponent;
Cf *= c[si] >= small || el >= 1 ? pow(max(c[si], 0), el) : 0;
}
forAll(rhs(), i)
{
const label si = rhs()[i].index;
const scalar sr = rhs()[i].stoichCoeff;
dcdt[si] += sr*omegaI;
const specieExponent& er = rhs()[i].exponent;
Cr *= c[si] >= small || er >= 1 ? pow(max(c[si], 0), er) : 0;
}
}
@ -325,225 +325,229 @@ Foam::scalar Foam::Reaction<ReactionThermo>::omega
scalar& omegar
) const
{
scalar clippedT = min(max(T, this->Tlow()), this->Thigh());
const scalar clippedT = min(max(T, this->Tlow()), this->Thigh());
omegaf = this->kf(p, clippedT, c, li);
omegar = this->kr(omegaf, p, clippedT, c, li);
// Rate constants
const scalar kf = this->kf(p, clippedT, c, li);
const scalar kr = this->kr(kf, p, clippedT, c, li);
forAll(lhs(), i)
{
const label si = lhs()[i].index;
const scalar el = lhs()[i].exponent;
omegaf *= c[si] >= small || el >= 1 ? pow(max(c[si], 0), el) : 0;
}
forAll(rhs(), i)
{
const label si = rhs()[i].index;
const scalar er = rhs()[i].exponent;
omegar *= c[si] >= small || er >= 1 ? pow(max(c[si], 0), er) : 0;
}
// Concentration products
scalar Cf, Cr;
this->C(p, T, c, li, Cf, Cr);
omegaf = kf*Cf;
omegar = kr*Cr;
return omegaf - omegar;
}
template<class ReactionThermo>
void Foam::Reaction<ReactionThermo>::dwdc
void Foam::Reaction<ReactionThermo>::dNdtByV
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarSquareMatrix& J,
scalarField& dcdt,
scalar& omegaI,
scalar& kfwd,
scalar& kbwd,
scalarField& dNdtByV,
const bool reduced,
const List<label>& c2s
const List<label>& c2s,
const label Nsi0
) const
{
scalar omegaf, omegar;
omegaI = omega(p, T, c, li, omegaf, omegar);
const scalar omega = this->omega(p, T, c, li, omegaf, omegar);
forAll(lhs(), i)
{
const label si = reduced ? c2s[lhs()[i].index] : lhs()[i].index;
const scalar sl = lhs()[i].stoichCoeff;
dcdt[si] -= sl*omegaI;
dNdtByV[Nsi0 + si] -= sl*omega;
}
forAll(rhs(), i)
{
const label si = reduced ? c2s[rhs()[i].index] : rhs()[i].index;
const scalar sr = rhs()[i].stoichCoeff;
dcdt[si] += sr*omegaI;
}
kfwd = this->kf(p, T, c, li);
kbwd = this->kr(kfwd, p, T, c, li);
forAll(lhs(), j)
{
const label sj = reduced ? c2s[lhs()[j].index] : lhs()[j].index;
scalar kf = kfwd;
forAll(lhs(), i)
{
const label si = lhs()[i].index;
const scalar el = lhs()[i].exponent;
if (i == j)
{
kf *=
c[si] >= small || el >= 1
? el*pow(max(c[si], 0), el - 1)
: 0;
}
else
{
kf *=
c[si] >= small || el >= 1
? pow(max(c[si], 0), el)
: 0;
}
}
forAll(lhs(), i)
{
const label si = reduced ? c2s[lhs()[i].index] : lhs()[i].index;
const scalar sl = lhs()[i].stoichCoeff;
J(si, sj) -= sl*kf;
}
forAll(rhs(), i)
{
const label si = reduced ? c2s[rhs()[i].index] : rhs()[i].index;
const scalar sr = rhs()[i].stoichCoeff;
J(si, sj) += sr*kf;
}
}
forAll(rhs(), j)
{
const label sj = reduced ? c2s[rhs()[j].index] : rhs()[j].index;
scalar kr = kbwd;
forAll(rhs(), i)
{
const label si = rhs()[i].index;
const scalar er = rhs()[i].exponent;
if (i == j)
{
kr *=
c[si] >= small || er >= 1
? er*pow(max(c[si], 0), er - 1)
: 0;
}
else
{
kr *=
c[si] >= small || er >= 1
? pow(max(c[si], 0), er)
: 0;
}
}
forAll(lhs(), i)
{
const label si = reduced ? c2s[lhs()[i].index] : lhs()[i].index;
const scalar sl = lhs()[i].stoichCoeff;
J(si, sj) += sl*kr;
}
forAll(rhs(), i)
{
const label si = reduced ? c2s[rhs()[i].index] : rhs()[i].index;
const scalar sr = rhs()[i].stoichCoeff;
J(si, sj) -= sr*kr;
}
}
// When third-body species are involved, additional terms are added
// beta function returns an empty list when third-body are not involved
const List<Tuple2<label, scalar>>& beta = this->beta();
if (notNull(beta))
{
// This temporary array needs to be cached for efficiency
scalarField dcidc(beta.size());
this->dcidc(p, T, c, li, dcidc);
forAll(beta, j)
{
label sj = beta[j].first();
sj = reduced ? c2s[sj] : sj;
if (sj != -1)
{
forAll(lhs(), i)
{
const label si =
reduced ? c2s[lhs()[i].index] : lhs()[i].index;
const scalar sl = lhs()[i].stoichCoeff;
J(si, sj) -= sl*dcidc[j]*omegaI;
}
forAll(rhs(), i)
{
const label si =
reduced ? c2s[rhs()[i].index] : rhs()[i].index;
const scalar sr = rhs()[i].stoichCoeff;
J(si, sj) += sr*dcidc[j]*omegaI;
}
}
}
dNdtByV[Nsi0 + si] += sr*omega;
}
}
template<class ReactionThermo>
void Foam::Reaction<ReactionThermo>::dwdT
void Foam::Reaction<ReactionThermo>::ddNdtByVdcTp
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
const scalar omegaI,
const scalar kfwd,
const scalar kbwd,
scalarSquareMatrix& J,
scalarField& dNdtByV,
scalarSquareMatrix& ddNdtByVdcTp,
const bool reduced,
const List<label>& c2s,
const label indexT
const label Nsi0,
const label Tsi,
scalarField& cTpWork0,
scalarField& cTpWork1
) const
{
scalar dkfdT = this->dkfdT(p, T, c, li);
scalar dkrdT = this->dkrdT(p, T, c, li, dkfdT, kbwd);
// Rate constants
const scalar kf = this->kf(p, T, c, li);
const scalar kr = this->kr(kf, p, T, c, li);
forAll(lhs(), i)
{
const label si = lhs()[i].index;
const scalar el = lhs()[i].exponent;
dkfdT *= c[si] >= small || el >= 1 ? pow(max(c[si], 0), el) : 0;
}
forAll(rhs(), i)
{
const label si = rhs()[i].index;
const scalar er = rhs()[i].exponent;
dkrdT *= c[si] >= small || er >= 1 ? pow(max(c[si], 0), er) : 0;
}
// Concentration products
scalar Cf, Cr;
this->C(p, T, c, li, Cf, Cr);
const scalar dqidT = dkfdT - dkrdT;
// Overall reaction rate
const scalar omega = kf*Cf - kr*Cr;
// For reactions including third-body efficiencies or pressure dependent
// reaction, an additional term is needed
const scalar dcidT = omegaI*this->dcidT(p, T, c, li);
// J(i, indexT) = sum_reactions nu_i dqdT
// Specie reaction rates
forAll(lhs(), i)
{
const label si = reduced ? c2s[lhs()[i].index] : lhs()[i].index;
const scalar sl = lhs()[i].stoichCoeff;
J(si, indexT) -= sl*(dqidT + dcidT);
dNdtByV[Nsi0 + si] -= sl*omega;
}
forAll(rhs(), i)
{
const label si = reduced ? c2s[rhs()[i].index] : rhs()[i].index;
const scalar sr = rhs()[i].stoichCoeff;
J(si, indexT) += sr*(dqidT + dcidT);
dNdtByV[Nsi0 + si] += sr*omega;
}
// Jacobian contributions from the derivative of the concentration products
// w.r.t. concentration
{
forAll(lhs(), j)
{
const label sj = reduced ? c2s[lhs()[j].index] : lhs()[j].index;
scalar dCfdcj = 1;
forAll(lhs(), i)
{
const label si = lhs()[i].index;
const specieExponent& el = lhs()[i].exponent;
if (i == j)
{
dCfdcj *=
c[si] >= small || el >= 1
? el*pow(max(c[si], 0), el - specieExponent(1))
: 0;
}
else
{
dCfdcj *=
c[si] >= small || el >= 1
? pow(max(c[si], 0), el)
: 0;
}
}
forAll(lhs(), i)
{
const label si = reduced ? c2s[lhs()[i].index] : lhs()[i].index;
const scalar sl = lhs()[i].stoichCoeff;
ddNdtByVdcTp(Nsi0 + si, Nsi0 + sj) -= sl*kf*dCfdcj;
}
forAll(rhs(), i)
{
const label si = reduced ? c2s[rhs()[i].index] : rhs()[i].index;
const scalar sr = rhs()[i].stoichCoeff;
ddNdtByVdcTp(Nsi0 + si, Nsi0 + sj) += sr*kf*dCfdcj;
}
}
forAll(rhs(), j)
{
const label sj = reduced ? c2s[rhs()[j].index] : rhs()[j].index;
scalar dCrcj = 1;
forAll(rhs(), i)
{
const label si = rhs()[i].index;
const specieExponent& er = rhs()[i].exponent;
if (i == j)
{
dCrcj *=
c[si] >= small || er >= 1
? er*pow(max(c[si], 0), er - specieExponent(1))
: 0;
}
else
{
dCrcj *=
c[si] >= small || er >= 1
? pow(max(c[si], 0), er)
: 0;
}
}
forAll(lhs(), i)
{
const label si = reduced ? c2s[lhs()[i].index] : lhs()[i].index;
const scalar sl = lhs()[i].stoichCoeff;
ddNdtByVdcTp(Nsi0 + si, Nsi0 + sj) += sl*kr*dCrcj;
}
forAll(rhs(), i)
{
const label si = reduced ? c2s[rhs()[i].index] : rhs()[i].index;
const scalar sr = rhs()[i].stoichCoeff;
ddNdtByVdcTp(Nsi0 + si, Nsi0 + sj) -= sr*kr*dCrcj;
}
}
}
// Jacobian contributions from the derivative of the rate constants
// w.r.t. temperature
{
const scalar dkfdT = this->dkfdT(p, T, c, li);
const scalar dkrdT = this->dkrdT(p, T, c, li, dkfdT, kr);
const scalar dwdT = dkfdT*Cf - dkrdT*Cr;
forAll(lhs(), i)
{
const label si = reduced ? c2s[lhs()[i].index] : lhs()[i].index;
const scalar sl = lhs()[i].stoichCoeff;
ddNdtByVdcTp(Nsi0 + si, Tsi) -= sl*dwdT;
}
forAll(rhs(), i)
{
const label si = reduced ? c2s[rhs()[i].index] : rhs()[i].index;
const scalar sr = rhs()[i].stoichCoeff;
ddNdtByVdcTp(Nsi0 + si, Tsi) += sr*dwdT;
}
}
// Jacobian contributions from the derivative of the rate constants
// w.r.t. concentration
if (hasDkdc())
{
scalarField& dkfdc = cTpWork0;
scalarField& dkrdc = cTpWork1;
this->dkfdc(p, T, c, li, dkfdc);
this->dkrdc(p, T, c, li, dkfdc, kr, dkrdc);
forAll(c, j)
{
const label sj = reduced ? c2s[j] : j;
if (sj == -1) continue;
const scalar dwdc = dkfdc[j]*Cf - dkrdc[j]*Cr;
forAll(lhs(), i)
{
const label si = reduced ? c2s[lhs()[i].index] : lhs()[i].index;
const scalar sl = lhs()[i].stoichCoeff;
ddNdtByVdcTp(Nsi0 + si, Nsi0 + sj) -= sl*dwdc;
}
forAll(rhs(), i)
{
const label si = reduced ? c2s[rhs()[i].index] : rhs()[i].index;
const scalar sr = rhs()[i].stoichCoeff;
ddNdtByVdcTp(Nsi0 + si, Nsi0 + sj) += sr*dwdc;
}
}
}
}

96
src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H
View File

@ -224,14 +224,15 @@ public:
// Reaction rate coefficients
//- Net reaction rate for individual species
void omega
//- Concentration powers
void C
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcdt
scalar& Cf,
scalar& Cr
) const;
//- Net reaction rate
@ -245,6 +246,20 @@ public:
scalar& omegar
) const;
//- The net reaction rate for each species involved
void dNdtByV
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dNdtByV,
const bool reduced,
const List<label>& c2s,
const label Nsi0
) const;
// Reaction rate coefficients
//- Forward rate constant
@ -259,7 +274,7 @@ public:
//- Reverse rate constant from the given forward rate constant
virtual scalar kr
(
const scalar kfwd,
const scalar kf,
const scalar p,
const scalar T,
const scalarField& c,
@ -278,40 +293,6 @@ public:
// Jacobian coefficients
//- Derivative of the net reaction rate for each species involved
// w.r.t. the species concentration
void dwdc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarSquareMatrix& J,
scalarField& dcdt,
scalar& omegaI,
scalar& kfwd,
scalar& kbwd,
const bool reduced,
const List<label>& c2s
) const;
//- Derivative of the net reaction rate for each species involved
// w.r.t. the temperature
void dwdT
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
const scalar omegaI,
const scalar kfwd,
const scalar kbwd,
scalarSquareMatrix& J,
const bool reduced,
const List<label>& c2s,
const label indexT
) const;
//- Temperature derivative of forward rate
virtual scalar dkfdT
(
@ -332,30 +313,49 @@ public:
const scalar kr
) const = 0;
//- Third-body efficiencies (beta = 1-alpha)
// non-empty only for third-body reactions
// with enhanced molecularity (alpha != 1)
virtual const List<Tuple2<label, scalar>>& beta() const = 0;
//- Does this reaction have concentration-dependent rate constants?
virtual bool hasDkdc() const = 0;
//- Species concentration derivative of the pressure dependent term
virtual void dcidc
//- Concentration derivative of forward rate
virtual void dkfdc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcidc
scalarField& dkfdc
) const = 0;
//- Temperature derivative of the pressure dependent term
virtual scalar dcidT
//- Concentration derivative of reverse rate
virtual void dkrdc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li
const label li,
const scalarField& dkfdc,
const scalar kr,
scalarField& dkrdc
) const = 0;
//- Derivative of the net reaction rate for each species involved
// w.r.t. the concentration and temperature
void ddNdtByVdcTp
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dNdtByV,
scalarSquareMatrix& ddNdtByVdcTp,
const bool reduced,
const List<label>& c2s,
const label csi0,
const label Tsi,
scalarField& cTpWork0,
scalarField& cTpWork1
) const;
//- Write
virtual void write(Ostream&) const;

17
src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.C
View File

@ -185,22 +185,21 @@ Foam::scalar Foam::ReactionProxy<ReactionThermo>::dkrdT
template<class ReactionThermo>
const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
Foam::ReactionProxy<ReactionThermo>::beta() const
bool Foam::ReactionProxy<ReactionThermo>::hasDkdc() const
{
NotImplemented;
return NullObjectRef<List<Tuple2<label, scalar>>>();
return false;
}
template<class ReactionThermo>
void Foam::ReactionProxy<ReactionThermo>::dcidc
void Foam::ReactionProxy<ReactionThermo>::dkfdc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcidc
scalarField& dkfdc
) const
{
NotImplemented;
@ -208,16 +207,18 @@ void Foam::ReactionProxy<ReactionThermo>::dcidc
template<class ReactionThermo>
Foam::scalar Foam::ReactionProxy<ReactionThermo>::dcidT
void Foam::ReactionProxy<ReactionThermo>::dkrdc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li
const label li,
const scalarField& dkfdc,
const scalar kr,
scalarField& dkrdc
) const
{
NotImplemented;
return 0;
}

21
src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.H
View File

@ -162,28 +162,29 @@ public:
const scalar kr
) const;
//- Third-body efficiencies (beta = 1-alpha)
// non-empty only for third-body reactions
// with enhanced molecularity (alpha != 1)
virtual const List<Tuple2<label, scalar>>& beta() const;
//- Does this reaction have concentration-dependent rate constants?
virtual bool hasDkdc() const;
//- Species concentration derivative of the pressure dependent term
virtual void dcidc
//- Concentration derivative of forward rate
void dkfdc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcidc
scalarField& dkfdc
) const;
//- Temperature derivative of the pressure dependent term
virtual scalar dcidT
//- Concentration derivative of reverse rate
void dkrdc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li
const label li,
const scalarField& dkfdc,
const scalar kr,
scalarField& dkrdc
) const;
};

26
src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.C
View File

@ -164,44 +164,48 @@ Foam::scalar Foam::ReversibleReaction<ReactionThermo, ReactionRate>::dkrdT
const scalar kr
) const
{
scalar Kc = max(this->Kc(p, T), rootSmall);
const scalar Kc = max(this->Kc(p, T), rootSmall);
return dkfdT/Kc - kr*this->dKcdTbyKc(p, T);
return dkfdT/Kc - (Kc > rootSmall ? kr*this->dKcdTbyKc(p, T) : 0);
}
template<class ReactionThermo, class ReactionRate>
const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
Foam::ReversibleReaction<ReactionThermo, ReactionRate>::beta() const
bool Foam::ReversibleReaction<ReactionThermo, ReactionRate>::hasDkdc() const
{
return k_.beta();
return k_.hasDdc();
}
template<class ReactionThermo, class ReactionRate>
void Foam::ReversibleReaction<ReactionThermo, ReactionRate>::dcidc
void Foam::ReversibleReaction<ReactionThermo, ReactionRate>::dkfdc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcidc
scalarField& dkfdc
) const
{
k_.dcidc(p, T, c, li, dcidc);
k_.ddc(p, T, c, li, dkfdc);
}
template<class ReactionThermo, class ReactionRate>
Foam::scalar Foam::ReversibleReaction<ReactionThermo, ReactionRate>::dcidT
void Foam::ReversibleReaction<ReactionThermo, ReactionRate>::dkrdc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li
const label li,
const scalarField& dkfdc,
const scalar kr,
scalarField& dkrdc
) const
{
return k_.dcidT(p, T, c, li);
const scalar Kc = max(this->Kc(p, T), rootSmall);
dkrdc = dkfdc/Kc;
}

25
src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.H
View File

@ -149,7 +149,7 @@ public:
const label li
) const;
//- Reverse rate constant from the given formard rate constant
//- Reverse rate constant from the given forward rate constant
virtual scalar kr
(
const scalar kfwd,
@ -193,30 +193,29 @@ public:
const scalar kr
) const;
//- Third-body efficiencies (beta = 1-alpha)
// non-empty only for third-body reactions
// with enhanced molecularity (alpha != 1)
virtual const List<Tuple2<label, scalar>>& beta() const;
//- Does this reaction have concentration-dependent rate constants?
virtual bool hasDkdc() const;
//- Species concentration derivative of the pressure dependent term
// By default this value is 1 as it multiplies the third-body term
virtual void dcidc
//- Concentration derivative of forward rate
void dkfdc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcidc
scalarField& dkfdc
) const;
//- Temperature derivative of the pressure dependent term
// By default this value is 0 since ddT of molecularity is approx.0
virtual scalar dcidT
//- Concentration derivative of reverse rate
void dkrdc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li
const label li,
const scalarField& dkfdc,
const scalar kr,
scalarField& dkrdc
) const;

25
src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRate.H
View File

@ -99,6 +99,7 @@ public:
//- Post-evaluation hook
inline void postEvaluate() const;
//- Return the rate
inline scalar operator()
(
const scalar p,
@ -107,6 +108,7 @@ public:
const label li
) const;
//- The derivative of the rate w.r.t. temperature
inline scalar ddT
(
const scalar p,
@ -115,30 +117,17 @@ public:
const label li
) const;
//- Third-body efficiencies (beta = 1-alpha)
// non-empty only for third-body reactions
// with enhanced molecularity (alpha != 1)
inline const List<Tuple2<label, scalar>>& beta() const;
//- Is the rate a function of concentration?
inline bool hasDdc() const;
//- Species concentration derivative of the pressure dependent term
// By default this value is 1 as it multiplies the third-body term
inline void dcidc
//- The derivative of the rate w.r.t. concentration
inline void ddc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcidc
) const;
//- Temperature derivative of the pressure dependent term
// By default this value is 0 since ddT of molecularity is approx.0
inline scalar dcidT
(
const scalar p,
const scalar T,
const scalarField& c,
const label li
scalarField& ddc
) const;
//- Write to stream

21
src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRateI.H
View File

@ -108,33 +108,22 @@ inline Foam::scalar Foam::ArrheniusReactionRate::ddT
}
inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
Foam::ArrheniusReactionRate::beta() const
inline bool Foam::ArrheniusReactionRate::hasDdc() const
{
return NullObjectRef<List<Tuple2<label, scalar>>>();
return false;
}
inline void Foam::ArrheniusReactionRate::dcidc
inline void Foam::ArrheniusReactionRate::ddc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcidc
) const
{}
inline Foam::scalar Foam::ArrheniusReactionRate::dcidT
(
const scalar p,
const scalar T,
const scalarField& c,
const label li
scalarField& ddc
) const
{
return 0;
ddc = 0;
}

29
src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRate.H
View File

@ -69,7 +69,6 @@ class ChemicallyActivatedReactionRate
ReactionRate kInf_;
ChemicallyActivationFunction F_;
thirdBodyEfficiencies thirdBodyEfficiencies_;
List<Tuple2<label, scalar>> beta_;
public:
@ -109,6 +108,7 @@ public:
//- Post-evaluation hook
inline void postEvaluate() const;
//- Return the rate
inline scalar operator()
(
const scalar p,
@ -117,6 +117,7 @@ public:
const label li
) const;
//- The derivative of the rate w.r.t. temperature
inline scalar ddT
(
const scalar p,
@ -125,33 +126,17 @@ public:
const label li
) const;
//- Third-body efficiencies (beta = 1-alpha)
// non-empty only for third-body reactions
// with enhanced molecularity (alpha != 1)
inline const List<Tuple2<label, scalar>>& beta() const
{
return beta_;
}
//- Is the rate a function of concentration?
inline bool hasDdc() const;
//- Species concentration derivative of the pressure dependent term
// By default this value is 1 as it multiplies the third-body term
inline void dcidc
//- The derivative of the rate w.r.t. concentration
inline void ddc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcidc
) const;
//- Temperature derivative of the pressure dependent term
// By default this value is 0 since ddT of molecularity is approx.0
inline scalar dcidT
(
const scalar p,
const scalar T,
const scalarField& c,
const label li
scalarField& ddc
) const;
//- Write to stream

170
src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRateI.H
View File

@ -23,6 +23,8 @@ License
\*---------------------------------------------------------------------------*/
#define CHEMICALLY_ACTIVATION_FUNCTION_JACOBIAN false
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class ReactionRate, class ChemicallyActivationFunction>
@ -42,12 +44,7 @@ inline Foam::ChemicallyActivatedReactionRate
kInf_(kInf),
F_(F),
thirdBodyEfficiencies_(tbes)
{
forAll(tbes, i)
{
beta_.append(Tuple2<label, scalar>(i, tbes[i]));
}
}
{}
template<class ReactionRate, class ChemicallyActivationFunction>
@ -61,23 +58,11 @@ inline Foam::ChemicallyActivatedReactionRate
const dictionary& dict
)
:
k0_(species, dict),
kInf_(species, dict),
F_(dict),
thirdBodyEfficiencies_(species, dict)
{
forAll(thirdBodyEfficiencies_, i)
{
beta_.append
(
Tuple2<label, scalar>
(
i,
thirdBodyEfficiencies_[i]
)
);
}
}
k0_(species, dict.subDict("k0")),
kInf_(species, dict.subDict("kInf")),
F_(dict.subDict("F")),
thirdBodyEfficiencies_(species, dict.subDict("thirdBodyEfficiencies"))
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -121,9 +106,11 @@ inline Foam::scalar Foam::ChemicallyActivatedReactionRate
{
const scalar k0 = k0_(p, T, c, li);
const scalar kInf = kInf_(p, T, c, li);
const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
const scalar M = thirdBodyEfficiencies_.M(c);
const scalar Pr = k0/kInf*M;
const scalar F = F_(T, Pr);
return k0*(1/(1 + Pr))*F_(T, Pr);
return k0/(1 + Pr)*F;
}
@ -142,9 +129,39 @@ inline Foam::scalar Foam::ChemicallyActivatedReactionRate
{
const scalar k0 = k0_(p, T, c, li);
const scalar kInf = kInf_(p, T, c, li);
const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
const scalar M = thirdBodyEfficiencies_.M(c);
const scalar Pr = k0/kInf*M;
const scalar F = F_(T, Pr);
return (1/(1 + Pr))*F_(T, Pr)*k0_.ddT(p, T, c, li);
const scalar dk0dT = k0_.ddT(p, T, c, li);
#if CHEMICALLY_ACTIVATION_FUNCTION_JACOBIAN
const scalar dkInfdT = kInf_.ddT(p, T, c, li);
const scalar dPrdT = (M*dk0dT - Pr*dkInfdT)/kInf;
const scalar dFdT = F_.ddT(T, Pr, F) + F_.ddPr(T, Pr, F)*dPrdT;
return
dk0dT/(1 + Pr)*F
- k0*dPrdT/sqr(1 + Pr)*F
+ k0/(1 + Pr)*dFdT;
#else
return dk0dT/(1 + Pr)*F;
#endif
}
template<class ReactionRate, class ChemicallyActivationFunction>
inline bool Foam::ChemicallyActivatedReactionRate
<
ReactionRate,
ChemicallyActivationFunction
>::hasDdc() const
{
return true;
}
@ -153,74 +170,44 @@ inline void Foam::ChemicallyActivatedReactionRate
<
ReactionRate,
ChemicallyActivationFunction
>::dcidc
>::ddc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcidc
scalarField& ddc
) const
{
const scalar k0 = k0_(p, T, c, li);
const scalar kInf = kInf_(p, T, c, li);
const scalar M = thirdBodyEfficiencies_.M(c);
const scalar Pr = k0/kInf*M;
const scalar F = F_(T, Pr);
if (M > small)
{
const scalar k0 = k0_(p, T, c, li);
const scalar kInf = kInf_(p, T, c, li);
#if CHEMICALLY_ACTIVATION_FUNCTION_JACOBIAN
const scalar Pr = k0*M/kInf;
const scalar F = F_(T, Pr);
scalarField dk0dc(c.size(), 0);
k0_.ddc(p, T, c, li, dk0dc);
scalarField dkInfdc(c.size(), 0);
kInf_.ddc(p, T, c, li, dkInfdc);
tmp<scalarField> tdMdc = thirdBodyEfficiencies_.dMdc(c);
const scalarField& dMdc = tdMdc();
scalarField dPrdc((M*dk0dc - Pr*dkInfdc + k0*dMdc)/kInf);
const scalar dFdPr = F_.ddPr(T, Pr, F);
forAll(beta_, i)
{
const scalar dPrdci = -beta_[i].second()*k0/kInf;
const scalar dFdci = F_.ddc(Pr, F, dPrdci, T);
dcidc[i] = (-dPrdci/(1 + Pr) + dFdci/F);
}
}
else
{
forAll(beta_, i)
{
dcidc[i] = 0;
}
}
}
ddc =
dk0dc/(1 + Pr)*F
- k0*dPrdc/sqr(1 + Pr)*F
+ k0/(1 + Pr)*dFdPr*dPrdc;
#else
template<class ReactionRate, class ChemicallyActivationFunction>
inline Foam::scalar Foam::ChemicallyActivatedReactionRate
<
ReactionRate,
ChemicallyActivationFunction
>::dcidT
(
const scalar p,
const scalar T,
const scalarField& c,
const label li
) const
{
const scalar M = thirdBodyEfficiencies_.M(c);
k0_.ddc(p, T, c, li, ddc);
if (M > small)
{
const scalar k0 = k0_(p, T, c, li);
const scalar kInf = kInf_(p, T, c, li);
ddc *= 1/(1 + Pr)*F;
const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
const scalar F = F_(T, Pr);
const scalar dPrdT =
Pr*(k0_.ddT(p, T, c, li)/k0 - kInf_.ddT(p, T, c, li)/kInf - 1/T);
const scalar dFdT = F_.ddT(Pr, F, dPrdT, T);
return (-dPrdT/(1 + Pr) + dFdT/F);
}
else
{
return 0;
}
#endif
}
@ -231,10 +218,33 @@ inline void Foam::ChemicallyActivatedReactionRate
ChemicallyActivationFunction
>::write(Ostream& os) const
{
os << indent << "k0" << nl;
os << indent << token::BEGIN_BLOCK << nl;
os << incrIndent;
k0_.write(os);
os << decrIndent;
os << indent << token::END_BLOCK << nl;
os << indent << "kInf" << nl;
os << indent << token::BEGIN_BLOCK << nl;
os << incrIndent;
kInf_.write(os);
os << decrIndent;
os << indent << token::END_BLOCK << nl;
os << indent << "F" << nl;
os << indent << token::BEGIN_BLOCK << nl;
os << incrIndent;
F_.write(os);
os << decrIndent;
os << indent << token::END_BLOCK << nl;
os << indent << "thirdBodyEfficiencies" << nl;
os << indent << token::BEGIN_BLOCK << nl;
os << incrIndent;
thirdBodyEfficiencies_.write(os);
os << decrIndent;
os << indent << token::END_BLOCK << nl;
}

30
src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRate.H
View File

@ -65,10 +65,12 @@ class FallOffReactionRate
// Private Data
ReactionRate k0_;
ReactionRate kInf_;
FallOffFunction F_;
thirdBodyEfficiencies thirdBodyEfficiencies_;
List<Tuple2<label, scalar>> beta_;
public:
@ -106,6 +108,7 @@ public:
//- Post-evaluation hook
inline void postEvaluate() const;
//- Return the rate
inline scalar operator()
(
const scalar p,
@ -114,6 +117,7 @@ public:
const label li
) const;
//- The derivative of the rate w.r.t. temperature
inline scalar ddT
(
const scalar p,
@ -122,31 +126,17 @@ public:
const label li
) const;
//- Third-body efficiencies (beta = 1-alpha)
// non-empty only for third-body reactions
// with enhanced molecularity (alpha != 1)
inline const List<Tuple2<label, scalar>>& beta() const
{
return beta_;
}
//- Is the rate a function of concentration?
inline bool hasDdc() const;
//- Species concentration derivative of the pressure dependent term
inline void dcidc
//- The derivative of the rate w.r.t. concentration
inline void ddc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcidc
) const;
//- Temperature derivative of the pressure dependent term
inline scalar dcidT
(
const scalar p,
const scalar T,
const scalarField& c,
const label li
scalarField& ddc
) const;
//- Write to stream

132
src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRateI.H
View File

@ -23,6 +23,8 @@ License
\*---------------------------------------------------------------------------*/
#define FALL_OFF_FUNCTION_JACOBIAN false
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class ReactionRate, class FallOffFunction>
@ -39,12 +41,7 @@ FallOffReactionRate
kInf_(kInf),
F_(F),
thirdBodyEfficiencies_(tbes)
{
forAll(tbes, i)
{
beta_.append(Tuple2<label, scalar>(i, tbes[i]));
}
}
{}
template<class ReactionRate, class FallOffFunction>
@ -59,19 +56,7 @@ FallOffReactionRate
kInf_(species, dict.subDict("kInf")),
F_(dict.subDict("F")),
thirdBodyEfficiencies_(species, dict.subDict("thirdBodyEfficiencies"))
{
forAll(thirdBodyEfficiencies_, i)
{
beta_.append
(
Tuple2<label, scalar>
(
i,
thirdBodyEfficiencies_[i]
)
);
}
}
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -106,9 +91,11 @@ Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::operator()
{
const scalar k0 = k0_(p, T, c, li);
const scalar kInf = kInf_(p, T, c, li);
const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
const scalar M = thirdBodyEfficiencies_.M(c);
const scalar Pr = k0/kInf*M;
const scalar F = F_(T, Pr);
return kInf*(Pr/(1 + Pr))*F_(T, Pr);
return kInf*(Pr/(1 + Pr))*F;
}
@ -124,77 +111,78 @@ Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::ddT
{
const scalar k0 = k0_(p, T, c, li);
const scalar kInf = kInf_(p, T, c, li);
const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
const scalar M = thirdBodyEfficiencies_.M(c);
const scalar Pr = k0/kInf*M;
const scalar F = F_(T, Pr);
return (Pr/(1 + Pr))*F_(T, Pr)*kInf_.ddT(p, T, c, li);
const scalar dkInfdT = kInf_.ddT(p, T, c, li);
#if FALL_OFF_FUNCTION_JACOBIAN
const scalar dk0dT = k0_.ddT(p, T, c, li);
const scalar dPrdT = (M*dk0dT - Pr*dkInfdT)/kInf;
const scalar dFdT = F_.ddT(T, Pr, F) + F_.ddPr(T, Pr, F)*dPrdT;
return
dkInfdT*(Pr/(1 + Pr))*F
+ kInf*dPrdT/sqr(1 + Pr)*F
+ kInf*(Pr/(1 + Pr))*dFdT;
#else
return dkInfdT*(Pr/(1 + Pr))*F;
#endif
}
template<class ReactionRate, class FallOffFunction>
inline void Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::dcidc
inline bool
Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::hasDdc() const
{
return true;
}
template<class ReactionRate, class FallOffFunction>
inline void Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::ddc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcidc
scalarField& ddc
) const
{
const scalar k0 = k0_(p, T, c, li);
const scalar kInf = kInf_(p, T, c, li);
const scalar M = thirdBodyEfficiencies_.M(c);
const scalar Pr = k0/kInf*M;
const scalar F = F_(T, Pr);
if (M > small)
{
const scalar k0 = k0_(p, T, c, li);
const scalar kInf = kInf_(p, T, c, li);
const scalar Pr = k0*M/kInf;
const scalar F = F_(T, Pr);
#if FALL_OFF_FUNCTION_JACOBIAN
forAll(beta_, i)
{
const scalar dPrdci = -beta_[i].second()*k0/kInf;
const scalar dFdci = F_.ddc(Pr, F, dPrdci, T);
dcidc[i] = (dPrdci/(Pr*(1 + Pr)) + dFdci/F);
}
}
else
{
forAll(beta_, i)
{
dcidc[i] = 0;
}
}
}
scalarField dk0dc(c.size(), 0);
k0_.ddc(p, T, c, li, dk0dc);
scalarField dkInfdc(c.size(), 0);
kInf_.ddc(p, T, c, li, dkInfdc);
tmp<scalarField> tdMdc = thirdBodyEfficiencies_.dMdc(c);
const scalarField& dMdc = tdMdc();
scalarField dPrdc((M*dk0dc - Pr*dkInfdc + k0*dMdc)/kInf);
const scalar dFdPr = F_.ddPr(T, Pr, F);
ddc =
dkInfdc*(Pr/(1 + Pr))*F
+ kInf*dPrdc/sqr(1 + Pr)*F
+ kInf*(Pr/(1 + Pr))*dFdPr*dPrdc;
template<class ReactionRate, class FallOffFunction>
inline Foam::scalar
Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::dcidT
(
const scalar p,
const scalar T,
const scalarField& c,
const label li
) const
{
const scalar M = thirdBodyEfficiencies_.M(c);
#else
if (M > small)
{
const scalar k0 = k0_(p, T, c, li);
const scalar kInf = kInf_(p, T, c, li);
kInf_.ddc(p, T, c, li, ddc);
const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
const scalar F = F_(T, Pr);
const scalar dPrdT =
Pr*(k0_.ddT(p, T, c, li)/k0 - kInf_.ddT(p, T, c, li)/kInf - 1/T);
const scalar dFdT = F_.ddT(Pr, F, dPrdT, T);
ddc *= (Pr/(1 + Pr))*F;
return (dPrdT/(Pr*(1 + Pr)) + dFdT/F);
}
else
{
return 0;
}
#endif
}

23
src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRate.H
View File

@ -102,6 +102,7 @@ public:
//- Post-evaluation hook
inline void postEvaluate() const;
//- Return the rate
inline scalar operator()
(
const scalar p,
@ -110,6 +111,7 @@ public:
const label li
) const;
//- The derivative of the rate w.r.t. temperature
inline scalar ddT
(
const scalar p,
@ -118,28 +120,17 @@ public:
const label li
) const;
//- Third-body efficiencies (beta = 1-alpha)
// non-empty only for third-body reactions
// with enhanced molecularity (alpha != 1)
inline const List<Tuple2<label, scalar>>& beta() const;
//- Is the rate a function of concentration?
inline bool hasDdc() const;
//- Species concentration derivative of the pressure dependent term
inline void dcidc
//- The derivative of the rate w.r.t. concentration
inline void ddc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcidc
) const;
//- Temperature derivative of the pressure dependent term
inline scalar dcidT
(
const scalar p,
const scalar T,
const scalarField& c,
const label li
scalarField& ddc
) const;
//- Write to stream

21
src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRateI.H
View File

@ -140,33 +140,22 @@ inline Foam::scalar Foam::JanevReactionRate::ddT
}
inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
Foam::JanevReactionRate::beta() const
inline bool Foam::JanevReactionRate::hasDdc() const
{
return NullObjectRef<List<Tuple2<label, scalar>>>();
return false;
}
inline void Foam::JanevReactionRate::dcidc
inline void Foam::JanevReactionRate::ddc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcidc
) const
{}
inline Foam::scalar Foam::JanevReactionRate::dcidT
(
const scalar p,
const scalar T,
const scalarField& c,
const label li
scalarField& ddc
) const
{
return 0;
ddc = 0;
}

23
src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRate.H
View File

@ -101,6 +101,7 @@ public:
//- Post-evaluation hook
inline void postEvaluate() const;
//- Return the rate
inline scalar operator()
(
const scalar p,
@ -109,6 +110,7 @@ public:
const label li
) const;
//- The derivative of the rate w.r.t. temperature
inline scalar ddT
(
const scalar p,
@ -117,28 +119,17 @@ public:
const label li
) const;
//- Third-body efficiencies (beta = 1-alpha)
// non-empty only for third-body reactions
// with enhanced molecularity (alpha != 1)
inline const List<Tuple2<label, scalar>>& beta() const;
//- Is the rate a function of concentration?
inline bool hasDdc() const;
//- Species concentration derivative of the pressure dependent term
inline void dcidc
//- The derivative of the rate w.r.t. concentration
inline void ddc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcidc
) const;
//- Temperature derivative of the pressure dependent term
inline scalar dcidT
(
const scalar p,
const scalar T,
const scalarField& c,
const label li
scalarField& ddc
) const;
//- Write to stream

21
src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRateI.H
View File

@ -155,33 +155,22 @@ inline Foam::scalar Foam::LandauTellerReactionRate::ddT
}
inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
Foam::LandauTellerReactionRate::beta() const
inline bool Foam::LandauTellerReactionRate::hasDdc() const
{
return NullObjectRef<List<Tuple2<label, scalar>>>();
return false;
}
inline void Foam::LandauTellerReactionRate::dcidc
inline void Foam::LandauTellerReactionRate::ddc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcidc
) const
{}
inline Foam::scalar Foam::LandauTellerReactionRate::dcidT
(
const scalar p,
const scalar T,
const scalarField& c,
const label li
scalarField& ddc
) const
{
return 0;
ddc = 0;
}

23
src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRate.H
View File

@ -104,6 +104,7 @@ public:
//- Post-evaluation hook
inline void postEvaluate() const;
//- Return the rate
inline scalar operator()
(
const scalar p,
@ -112,6 +113,7 @@ public:
const label li
) const;
//- The derivative of the rate w.r.t. temperature
inline scalar ddT
(
const scalar p,
@ -120,28 +122,17 @@ public:
const label li
) const;
//- Third-body efficiencies (beta = 1-alpha)
// non-empty only for third-body reactions
// with enhanced molecularity (alpha != 1)
inline const List<Tuple2<label, scalar>>& beta() const;
//- Is the rate a function of concentration?
inline bool hasDdc() const;
//- Species concentration derivative of the pressure dependent term
inline void dcidc
//- The derivative of the rate w.r.t. concentration
inline void ddc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcidc
) const;
//- Temperature derivative of the pressure dependent term
inline scalar dcidT
(
const scalar p,
const scalar T,
const scalarField& c,
const label li
scalarField& ddc
) const;
//- Write to stream

73
src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRateI.H
View File

@ -83,15 +83,13 @@ inline Foam::scalar Foam::LangmuirHinshelwoodReactionRate::operator()
const scalar c1m = pow(c[r_[0]], m_[1]);
const scalar c2m = pow(c[r_[1]], m_[2]);
const scalar TaByT0 = Ta_[0]/T;
const scalar TaByT1 = Ta_[1]/T;
const scalar TaByT2 = Ta_[2]/T;
const scalar k0 = A_[0]*pow(T, beta_[0])*exp(-Ta_[0]/T);
const scalar k1 = A_[1]*pow(T, beta_[1])*exp(-Ta_[1]/T);
const scalar k2 = A_[2]*pow(T, beta_[2])*exp(-Ta_[2]/T);
const scalar k0 = A_[0]*pow(T, beta_[0])*exp(-TaByT0);
const scalar k1 = A_[1]*pow(T, beta_[1])*exp(-TaByT1);
const scalar k2 = A_[2]*pow(T, beta_[2])*exp(-TaByT2);
const scalar b = a_ + k1*c1m + k2*c2m;
return k0/pow(a_ + k1*c1m + k2*c2m, m_[0]);
return k0/pow(b, m_[0]);
}
@ -106,50 +104,55 @@ inline Foam::scalar Foam::LangmuirHinshelwoodReactionRate::ddT
const scalar c1m = pow(c[r_[0]], m_[1]);
const scalar c2m = pow(c[r_[1]], m_[2]);
const scalar TaByT0 = Ta_[0]/T;
const scalar TaByT1 = Ta_[1]/T;
const scalar TaByT2 = Ta_[2]/T;
const scalar k0 = A_[0]*pow(T, beta_[0])*exp(-Ta_[0]/T);
const scalar k1 = A_[1]*pow(T, beta_[1])*exp(-Ta_[1]/T);
const scalar k2 = A_[2]*pow(T, beta_[2])*exp(-Ta_[2]/T);
const scalar k0 = A_[0]*pow(T, beta_[0])*exp(-TaByT0);
const scalar k1 = A_[1]*pow(T, beta_[1])*exp(-TaByT1);
const scalar k2 = A_[2]*pow(T, beta_[2])*exp(-TaByT2);
const scalar dk0dT = k0/T*(beta_[0] + Ta_[0]/T);
const scalar dk1dT = k1/T*(beta_[1] + Ta_[1]/T);
const scalar dk2dT = k2/T*(beta_[2] + Ta_[2]/T);
return
(
(beta_[0] + TaByT0)*k0
- m_[0]*k0*((beta_[1] + TaByT1)*k1*c1m + (beta_[2] + TaByT2)*k2*c2m)
/(a_ + k1*c1m + k2*c2m)
)/(pow(a_ + k1*c1m + k2*c2m, m_[0])*T);
const scalar b = a_ + k1*c1m + k2*c2m;
const scalar dbdT = dk1dT*c1m + dk2dT*c2m;
return (dk0dT - k0*m_[0]*dbdT/b)/pow(b, m_[0]);
}
inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
Foam::LangmuirHinshelwoodReactionRate::beta() const
inline bool Foam::LangmuirHinshelwoodReactionRate::hasDdc() const
{
return NullObjectRef<List<Tuple2<label, scalar>>>();
return true;
}
inline void Foam::LangmuirHinshelwoodReactionRate::dcidc
inline void Foam::LangmuirHinshelwoodReactionRate::ddc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcidc
) const
{}
inline Foam::scalar Foam::LangmuirHinshelwoodReactionRate::dcidT
(
const scalar p,
const scalar T,
const scalarField& c,
const label li
scalarField& ddc
) const
{
return 0;
const scalar c1m = pow(c[r_[0]], m_[1]);
const scalar c2m = pow(c[r_[1]], m_[2]);
const scalar dc1mdc1 = m_[1]*c1m/c[r_[0]];
const scalar dc2mdc2 = m_[2]*c2m/c[r_[1]];
const scalar k0 = A_[0]*pow(T, beta_[0])*exp(-Ta_[0]/T);
const scalar k1 = A_[1]*pow(T, beta_[1])*exp(-Ta_[1]/T);
const scalar k2 = A_[2]*pow(T, beta_[2])*exp(-Ta_[2]/T);
const scalar b = a_ + k1*c1m + k2*c2m;
const scalar dbdc1 = k1*dc1mdc1;
const scalar dbdc2 = k2*dc2mdc2;
const scalar k = k0/pow(b, m_[0]);
ddc = 0;
ddc[r_[0]] = - k*m_[0]/b*dbdc1;
ddc[r_[1]] = - k*m_[1]/b*dbdc2;
}

25
src/thermophysicalModels/specie/reaction/reactionRate/MichaelisMenten/MichaelisMentenReactionRate.H
View File

@ -77,8 +77,6 @@ class MichaelisMentenReactionRate
//- The substrate specie index
label s_;
List<Tuple2<label, scalar>> beta_;
public:
@ -106,6 +104,7 @@ public:
//- Post-evaluation hook
inline void postEvaluate() const;
//- Return the rate
inline scalar operator()
(
const scalar p,
@ -114,6 +113,7 @@ public:
const label li
) const;
//- The derivative of the rate w.r.t. temperature
inline scalar ddT
(
const scalar p,
@ -122,28 +122,17 @@ public:
const label li
) const;
//- Third-body efficiencies (beta = 1-alpha)
// non-empty only for third-body reactions
// with enhanced molecularity (alpha != 1)
inline const List<Tuple2<label, scalar>>& beta() const;
//- Is the rate a function of concentration?
inline bool hasDdc() const;
//- Species concentration derivative of the pressure dependent term
inline void dcidc
//- The derivative of the rate w.r.t. concentration
inline void ddc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcidc
) const;
//- Temperature derivative of the pressure dependent term
inline scalar dcidT
(
const scalar p,
const scalar T,
const scalarField& c,
const label li
scalarField& ddc
) const;
//- Write to stream

28
src/thermophysicalModels/specie/reaction/reactionRate/MichaelisMenten/MichaelisMentenReactionRateI.H
View File

@ -35,9 +35,7 @@ inline Foam::MichaelisMentenReactionRate::MichaelisMentenReactionRate
Vmax_(dict.lookup<scalar>("Vmax")),
Km_(dict.lookup<scalar>("Km")),
s_(st[dict.lookup("S")])
{
beta_.append(Tuple2<label, scalar>(s_, 1.0));
}
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -74,35 +72,23 @@ inline Foam::scalar Foam::MichaelisMentenReactionRate::ddT
}
inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
Foam::MichaelisMentenReactionRate::beta() const
inline bool Foam::MichaelisMentenReactionRate::hasDdc() const
{
return beta_;
return true;
}
inline void Foam::MichaelisMentenReactionRate::dcidc
inline void Foam::MichaelisMentenReactionRate::ddc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcidc
scalarField& ddc
) const
{
dcidc[0] = -1.0/(Km_ + c[s_]);
}
inline Foam::scalar Foam::MichaelisMentenReactionRate::dcidT
(
const scalar p,
const scalar T,
const scalarField& c,
const label li
) const
{
return 0;
ddc = 0;
ddc[s_] = - Vmax_/sqr(Km_ + c[s_]);
}

14
src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/LindemannFallOffFunction/LindemannFallOffFunction.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -82,18 +82,16 @@ public:
inline scalar ddT
(
const scalar T,
const scalar Pr,
const scalar F,
const scalar dPrdT,
const scalar T
const scalar F
) const;
inline scalar ddc
inline scalar ddPr
(
const scalar T,
const scalar Pr,
const scalar F,
const scalar dPrdc,
const scalar T
const scalar F
) const;
//- Write to stream

6
src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/LindemannFallOffFunction/LindemannFallOffFunctionI.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -50,7 +50,6 @@ inline Foam::scalar Foam::LindemannFallOffFunction::operator()
inline Foam::scalar Foam::LindemannFallOffFunction::ddT
(
const scalar,
const scalar,
const scalar,
const scalar
@ -60,9 +59,8 @@ inline Foam::scalar Foam::LindemannFallOffFunction::ddT
}
inline Foam::scalar Foam::LindemannFallOffFunction::ddc
inline Foam::scalar Foam::LindemannFallOffFunction::ddPr
(
const scalar,
const scalar,
const scalar,
const scalar

14
src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/SRIFallOffFunction/SRIFallOffFunction.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -93,18 +93,16 @@ public:
inline scalar ddT
(
const scalar T,
const scalar Pr,
const scalar F,
const scalar dPrdT,
const scalar T
const scalar F
) const;
inline scalar ddc
inline scalar ddPr
(
const scalar T,
const scalar Pr,
const scalar F,
const scalar dPrdc,
const scalar T
const scalar F
) const;
//- Write to stream

57
src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/SRIFallOffFunction/SRIFallOffFunctionI.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -60,45 +60,52 @@ inline Foam::scalar Foam::SRIFallOffFunction::operator()
const scalar Pr
) const
{
scalar X = 1.0/(1 + sqr(log10(max(Pr, small))));
return d_*pow(a_*exp(-b_/T) + exp(-T/c_), X)*pow(T, e_);
const scalar logPr = log10(max(Pr, small));
const scalar X = 1/(1 + sqr(logPr));
const scalar psi = a_*exp(-b_/T) + exp(-T/c_);
return d_*pow(psi, X)*pow(T, e_);
}
inline Foam::scalar Foam::SRIFallOffFunction::ddT
(
const scalar T,
const scalar Pr,
const scalar F,
const scalar dPrdT,
const scalar T
const scalar F
) const
{
scalar X = 1.0/(1 + sqr(log10(max(Pr, small))));
scalar dXdPr = -X*X*2*log10(Pr)/Pr/log(10.0);
return
F
*(
e_/T
+ X
*(a_*b_*exp(-b_/T)/sqr(T) - exp(-T/c_)/c_)
/(a_*exp(-b_/T) + exp(-T/c_))
+ dXdPr*dPrdT*log(a_*exp(-b_/T) + exp(-T/c_))
);
const scalar logPr = log10(max(Pr, small));
const scalar X = 1/(1 + sqr(logPr));
const scalar psi = a_*exp(-b_/T) + exp(-T/c_);
const scalar dpsidT = a_*b_/sqr(T)*exp(-b_/T) - 1/c_*exp(-T/c_);
return F*(X/psi*dpsidT + e_/T);
}
inline Foam::scalar Foam::SRIFallOffFunction::ddc
inline Foam::scalar Foam::SRIFallOffFunction::ddPr
(
const scalar T,
const scalar Pr,
const scalar F,
const scalar dPrdc,
const scalar T
const scalar F
) const
{
scalar X = 1.0/(1 + sqr(log10(max(Pr, small))));
scalar dXdPr = -X*X*2*log10(Pr)/Pr/log(10.0);
scalar dXdc = dXdPr*dPrdc;
return F*dXdc*log(a_*exp(-b_/T) + exp(-T/c_));
static const scalar logTen = log(scalar(10));
const scalar logPr = log10(max(Pr, small));
const scalar dlogPrdPr = Pr >= small ? 1/(logTen*Pr) : 0;
const scalar X = 1/(1 + sqr(logPr));
const scalar dXdPr = -sqr(X)*2*logPr*dlogPrdPr;
const scalar psi = a_*exp(-b_/T) + exp(-T/c_);
return F*log(psi)*dXdPr;
}

14
src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/TroeFallOffFunction/TroeFallOffFunction.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -93,18 +93,16 @@ public:
inline scalar ddT
(
const scalar T,
const scalar Pr,
const scalar F,
const scalar dPrdT,
const scalar T
const scalar F
) const;
inline scalar ddc
inline scalar ddPr
(
const scalar T,
const scalar Pr,
const scalar F,
const scalar dPrdc,
const scalar T
const scalar F
) const;
//- Write to stream

157
src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/TroeFallOffFunction/TroeFallOffFunctionI.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -57,116 +57,103 @@ inline Foam::scalar Foam::TroeFallOffFunction::operator()
const scalar Pr
) const
{
scalar logFcent = log10
(
max
(
(1 - alpha_)*exp(-T/Tsss_) + alpha_*exp(-T/Ts_) + exp(-Tss_/T),
small
)
);
const scalar logPr = log10(max(Pr, small));
scalar c = -0.4 - 0.67*logFcent;
const scalar Fcent =
(1 - alpha_)*exp(-T/Tsss_)
+ alpha_*exp(-T/Ts_)
+ exp(-Tss_/T);
const scalar logFcent = log10(max(Fcent, small));
const scalar c = -0.4 - 0.67*logFcent;
static const scalar d = 0.14;
scalar n = 0.75 - 1.27*logFcent;
const scalar n = 0.75 - 1.27*logFcent;
scalar logPr = log10(max(Pr, small));
return pow(10, logFcent/(1 + sqr((logPr + c)/(n - d*(logPr + c)))));
const scalar x1 = n - d*(logPr + c);
const scalar x2 = (logPr + c)/x1;
const scalar x3 = 1 + sqr(x2);
const scalar x4 = logFcent/x3;
return pow(10, x4);
}
inline Foam::scalar Foam::TroeFallOffFunction::ddT
(
const scalar T,
const scalar Pr,
const scalar F,
const scalar dPrdT,
const scalar T
const scalar F
) const
{
scalar logPr = log10(max(Pr, small));
scalar logTen = log(10.0);
scalar Fcent =
(
max
(
(1 - alpha_)*exp(-T/Tsss_) + alpha_*exp(-T/Ts_) + exp(-Tss_/T),
small
)
);
scalar logFcent = log10(Fcent);
static const scalar logTen = log(scalar(10));
scalar dFcentdT =
(
(alpha_ - 1)*exp(-T/Tsss_)/Tsss_
- alpha_*exp(-T/Ts_)/Ts_
+ Tss_*exp(-Tss_/T)/sqr(T)
);
const scalar logPr = log10(max(Pr, small));
scalar d = 0.14;
scalar dlogFcentdT = dFcentdT/Fcent/logTen;
scalar c = -0.4 - 0.67*logFcent;
scalar dcdT = -0.67*dlogFcentdT;
scalar n = 0.75 - 1.27*logFcent;
scalar dndT = -1.27*dlogFcentdT;
const scalar Fcent =
(1 - alpha_)*exp(-T/Tsss_)
+ alpha_*exp(-T/Ts_)
+ exp(-Tss_/T);
const scalar dFcentdT =
- (1 - alpha_)/Tsss_*exp(-T/Tsss_)
- alpha_/Ts_*exp(-T/Ts_)
+ Tss_/sqr(T)*exp(-Tss_/T);
scalar dlogPrdT = dPrdT/Pr/logTen;
const scalar logFcent = log10(max(Fcent, small));
const scalar dlogFcentdT = Fcent >= small ? dFcentdT/Fcent/logTen : 0;
scalar dParentdT =
2.0*(logPr + c)/sqr(n - d*(logPr + c))
*(
(dlogPrdT + dcdT)
- (logPr + c)*(dndT - d*(dlogPrdT + dcdT))/(n - d*(logPr + c))
);
const scalar c = -0.4 - 0.67*logFcent;
const scalar dcdT = -0.67*dlogFcentdT;
static const scalar d = 0.14;
const scalar n = 0.75 - 1.27*logFcent;
const scalar dndT = - 1.27*dlogFcentdT;
return
(
F*logTen
*(
dlogFcentdT/(1.0 + sqr((logPr + c)/(n - d*(logPr + c))))
- logFcent*dParentdT/sqr(1.0 + sqr((logPr + c)/(n - d*(logPr + c))))
)
);
const scalar x1 = n - d*(logPr + c);
const scalar dx1dT = dndT - d*dcdT;
const scalar x2 = (logPr + c)/x1;
const scalar dx2dT = (dcdT - x2*dx1dT)/x1;
const scalar x3 = 1 + sqr(x2);
const scalar dx3dT = 2*x2*dx2dT;
const scalar x4 = logFcent/x3;
const scalar dx4dT = (dlogFcentdT - x4*dx3dT)/x3;
return logTen*F*dx4dT;
}
inline Foam::scalar Foam::TroeFallOffFunction::ddc
inline Foam::scalar Foam::TroeFallOffFunction::ddPr
(
const scalar T,
const scalar Pr,
const scalar F,
const scalar dPrdc,
const scalar T
const scalar F
) const
{
scalar logPr = log10(max(Pr, small));
scalar logTen = log(10.0);
scalar logFcent = log10
(
max
(
(1 - alpha_)*exp(-T/Tsss_) + alpha_*exp(-T/Ts_) + exp(-Tss_/T),
small
)
);
static const scalar logTen = log(scalar(10));
scalar dlogPrdc = dPrdc/Pr/logTen;
scalar d = 0.14;
scalar c = -0.4 - 0.67*logFcent;
scalar n = 0.75 - 1.27*logFcent;
const scalar logPr = log10(max(Pr, small));
const scalar dlogPrdPr = Pr >= small ? 1/(logTen*Pr) : 0;
scalar dParentdc =
2.0*(logPr + c)/sqr(n - d*(logPr + c))
*(
(dlogPrdc)
- (logPr + c)*(-d*(dlogPrdc))/(n - d*(logPr + c))
);
const scalar Fcent =
(1 - alpha_)*exp(-T/Tsss_)
+ alpha_*exp(-T/Ts_)
+ exp(-Tss_/T);
return
(
F*logTen
*(
- logFcent*dParentdc/sqr(1.0 + sqr((logPr + c)/(n - d*(logPr + c))))
)
);
const scalar logFcent = log10(max(Fcent, small));
const scalar c = -0.4 - 0.67*logFcent;
static const scalar d = 0.14;
const scalar n = 0.75 - 1.27*logFcent;
const scalar x1 = n - d*(logPr + c);
const scalar dx1dPr = -d*dlogPrdPr;
const scalar x2 = (logPr + c)/x1;
const scalar dx2dPr = (dlogPrdPr - x2*dx1dPr)/x1;
const scalar x3 = 1 + sqr(x2);
const scalar dx3dPr = 2*x2*dx2dPr;
const scalar x4 = logFcent/x3;
const scalar dx4dPr = -x4*dx3dPr/x3;
return logTen*F*dx4dPr;
}

23
src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRate.H
View File

@ -102,6 +102,7 @@ public:
//- Post-evaluation hook
inline void postEvaluate() const;
//- Return the rate
inline scalar operator()
(
const scalar p,
@ -110,6 +111,7 @@ public:
const label li
) const;
//- The derivative of the rate w.r.t. temperature
inline scalar ddT
(
const scalar p,
@ -118,28 +120,17 @@ public:
const label li
) const;
//- Third-body efficiencies (beta = 1-alpha)
// non-empty only for third-body reactions
// with enhanced molecularity (alpha != 1)
inline const List<Tuple2<label, scalar>>& beta() const;
//- Is the rate a function of concentration?
inline bool hasDdc() const;
//- Species concentration derivative of the pressure dependent term
inline void dcidc
//- The derivative of the rate w.r.t. concentration
inline void ddc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcidc
) const;
//- Temperature derivative of the pressure dependent term
inline scalar dcidT
(
const scalar p,
const scalar T,
const scalarField& c,
const label li
scalarField& ddc
) const;
//- Write to stream

21
src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRateI.H
View File

@ -122,33 +122,22 @@ inline Foam::scalar Foam::powerSeriesReactionRate::ddT
}
inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
Foam::powerSeriesReactionRate::beta() const
inline bool Foam::powerSeriesReactionRate::hasDdc() const
{
return NullObjectRef<List<Tuple2<label, scalar>>>();
return false;
}
inline void Foam::powerSeriesReactionRate::dcidc
inline void Foam::powerSeriesReactionRate::ddc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcidc
) const
{}
inline Foam::scalar Foam::powerSeriesReactionRate::dcidT
(
const scalar p,
const scalar T,
const scalarField& c,
const label li
scalarField& ddc
) const
{
return 0;
ddc = 0;
}

27
src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRate.H
View File

@ -61,7 +61,6 @@ class thirdBodyArrheniusReactionRate
// Private Data
thirdBodyEfficiencies thirdBodyEfficiencies_;
List<Tuple2<label, scalar>> beta_;
public:
@ -99,6 +98,7 @@ public:
//- Post-evaluation hook
inline void postEvaluate() const;
//- Return the rate
inline scalar operator()
(
const scalar p,
@ -107,6 +107,7 @@ public:
const label li
) const;
//- The derivative of the rate w.r.t. temperature
inline scalar ddT
(
const scalar p,
@ -115,31 +116,17 @@ public:
const label li
) const;
//- Third-body efficiencies (beta = 1-alpha)
// non-empty only for third-body reactions
// with enhanced molecularity (alpha != 1)
inline const List<Tuple2<label, scalar>>& beta() const
{
return beta_;
}
//- Is the rate a function of concentration?
inline bool hasDdc() const;
//- Species concentration derivative of the pressure dependent term
inline void dcidc
//- The derivative of the rate w.r.t. concentration
inline void ddc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcidc
) const;
//- Temperature derivative of the pressure dependent term
inline scalar dcidT
(
const scalar p,
const scalar T,
const scalarField& c,
const label li
scalarField& ddc
) const;
//- Write to stream

50
src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRateI.H
View File

@ -35,12 +35,7 @@ inline Foam::thirdBodyArrheniusReactionRate::thirdBodyArrheniusReactionRate
:
ArrheniusReactionRate(A, beta, Ta),
thirdBodyEfficiencies_(tbes)
{
forAll(tbes, i)
{
beta_.append(Tuple2<label, scalar>(i, tbes[i]));
}
}
{}
inline Foam::thirdBodyArrheniusReactionRate::thirdBodyArrheniusReactionRate
@ -55,19 +50,7 @@ inline Foam::thirdBodyArrheniusReactionRate::thirdBodyArrheniusReactionRate
dict
),
thirdBodyEfficiencies_(species, dict)
{
forAll(thirdBodyEfficiencies_, i)
{
beta_.append
(
Tuple2<label, scalar>
(
i,
thirdBodyEfficiencies_[i]
)
);
}
}
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -112,32 +95,25 @@ inline Foam::scalar Foam::thirdBodyArrheniusReactionRate::ddT
}
inline void Foam::thirdBodyArrheniusReactionRate::dcidc
inline bool Foam::thirdBodyArrheniusReactionRate::hasDdc() const
{
return true;
}
inline void Foam::thirdBodyArrheniusReactionRate::ddc
(
const scalar p,
const scalar T,
const scalarField& c,
const label li,
scalarField& dcidc
scalarField& ddc
) const
{
scalar M = thirdBodyEfficiencies_.M(c);
forAll(beta_, i)
{
dcidc[i] = beta_[i].second()/max(M, small);
}
}
const scalar k = ArrheniusReactionRate::operator()(p, T, c, li);
inline Foam::scalar Foam::thirdBodyArrheniusReactionRate::dcidT
(
const scalar p,
const scalar T,
const scalarField& c,
const label li
) const
{
return -1.0/T;
ddc = thirdBodyEfficiencies_.dMdc(c);
ddc *= k;
}

12
src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyEfficiencies/thirdBodyEfficiencies.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -35,7 +35,7 @@ SourceFiles
#ifndef thirdBodyEfficiencies_H
#define thirdBodyEfficiencies_H
#include "scalarList.H"
#include "scalarField.H"
#include "speciesTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -55,7 +55,7 @@ Ostream& operator<<(Ostream&, const thirdBodyEfficiencies&);
class thirdBodyEfficiencies
:
public scalarList
public scalarField
{
// Private Data
@ -84,7 +84,11 @@ public:
// Member Functions
//- Calculate and return M, the concentration of the third-bodies
inline scalar M(const scalarList& c) const;
inline scalar M(const scalarField& c) const;
//- Calculate and return the derivative of M, w.r.t. the species
// concentrations
inline tmp<scalarField> dMdc(const scalarField& c) const;
//- Write to stream
inline void write(Ostream& os) const;

17
src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyEfficiencies/thirdBodyEfficienciesI.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -33,7 +33,7 @@ inline Foam::thirdBodyEfficiencies::thirdBodyEfficiencies
const scalarList& efficiencies
)
:
scalarList(efficiencies),
scalarField(efficiencies),
species_(species)
{
if (size() != species_.size())
@ -52,7 +52,7 @@ inline Foam::thirdBodyEfficiencies::thirdBodyEfficiencies
const dictionary& dict
)
:
scalarList(species.size()),
scalarField(species.size()),
species_(species)
{
if (dict.found("coeffs"))
@ -74,14 +74,14 @@ inline Foam::thirdBodyEfficiencies::thirdBodyEfficiencies
else
{
scalar defaultEff = dict.lookup<scalar>("defaultEfficiency");
scalarList::operator=(defaultEff);
scalarField::operator=(defaultEff);
}
}
// * * * * * * * * * * * * * * * Member functions * * * * * * * * * * * * * //
inline Foam::scalar Foam::thirdBodyEfficiencies::M(const scalarList& c) const
inline Foam::scalar Foam::thirdBodyEfficiencies::M(const scalarField& c) const
{
scalar M = 0;
forAll(*this, i)
@ -93,6 +93,13 @@ inline Foam::scalar Foam::thirdBodyEfficiencies::M(const scalarList& c) const
}
inline Foam::tmp<Foam::scalarField>
Foam::thirdBodyEfficiencies::dMdc(const scalarField& c) const
{
return *this;
}
inline void Foam::thirdBodyEfficiencies::write(Ostream& os) const
{
List<Tuple2<word, scalar>> coeffs(species_.size());

4
src/thermophysicalModels/specie/reaction/specieCoeffs/specieCoeffs.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2019-2020 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2019-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -191,7 +191,7 @@ void Foam::specieCoeffs::reactionStr
reaction << scs[i].stoichCoeff;
}
reaction << species[scs[i].index];
if (mag(scs[i].exponent - scs[i].stoichCoeff) > small)
if (mag(scs[i].exponent.operator scalar() - scs[i].stoichCoeff) > small)
{
reaction << "^" << scs[i].exponent;
}

121
src/thermophysicalModels/specie/reaction/specieCoeffs/specieCoeffs.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2019-2020 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2019-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -36,6 +36,7 @@ SourceFiles
#ifndef specieCoeffs_H
#define specieCoeffs_H
#include "specieExponent.H"
#include "speciesTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -47,7 +48,6 @@ namespace Foam
class specieCoeffs;
Ostream& operator<<(Ostream&, const specieCoeffs&);
/*---------------------------------------------------------------------------*\
Class specieCoeffs Declaration
\*---------------------------------------------------------------------------*/
@ -56,62 +56,83 @@ class specieCoeffs
{
public:
label index;
scalar stoichCoeff;
scalar exponent;
// Public Data
specieCoeffs()
:
index(-1),
stoichCoeff(0),
exponent(1)
{}
//- Index of the specie
label index;
specieCoeffs(const speciesTable& species, Istream& is);
//- Stoichiometric coefficient
scalar stoichCoeff;
bool operator==(const specieCoeffs& sc) const
{
return index == sc.index;
}
//- Exponent of the specie concentration
specieExponent exponent;
bool operator!=(const specieCoeffs& sc) const
{
return index != sc.index;
}
friend Ostream& operator<<(Ostream& os, const specieCoeffs& sc)
{
os << sc.index << token::SPACE
<< sc.stoichCoeff << token::SPACE
<< sc.exponent;
return os;
}
// Constructors
//- Construct the left- and right-hand-side reaction coefficients
static void setLRhs
(
Istream&,
const speciesTable&,
List<specieCoeffs>& lhs,
List<specieCoeffs>& rhs
);
//- Construct null
specieCoeffs()
:
index(-labelMax),
stoichCoeff(NaN),
exponent()
{}
//- Write the string representation of the specieCoeffs list
static void reactionStr
(
OStringStream& reaction,
const speciesTable&,
const List<specieCoeffs>& scs
);
//- Construct from species table and input stream
specieCoeffs(const speciesTable& species, Istream& is);
//- Return string representation of reaction
static string reactionStr
(
OStringStream& reaction,
const speciesTable&,
const List<specieCoeffs>& lhs,
const List<specieCoeffs>& rhs
);
// Member Functions
//- Construct the left- and right-hand-side reaction coefficients
static void setLRhs
(
Istream&,
const speciesTable&,
List<specieCoeffs>& lhs,
List<specieCoeffs>& rhs
);
//- Write the string representation of the specieCoeffs list
static void reactionStr
(
OStringStream& reaction,
const speciesTable&,
const List<specieCoeffs>& scs
);
//- Return string representation of reaction
static string reactionStr
(
OStringStream& reaction,
const speciesTable&,
const List<specieCoeffs>& lhs,
const List<specieCoeffs>& rhs
);
// Member Operators
//- Equality comparison
bool operator==(const specieCoeffs& sc) const
{
return index == sc.index;
}
//- Inequality comparison
bool operator!=(const specieCoeffs& sc) const
{
return index != sc.index;
}
//- Write to output stream
friend Ostream& operator<<(Ostream& os, const specieCoeffs& sc)
{
os << sc.index << token::SPACE
<< sc.stoichCoeff << token::SPACE
<< sc.exponent;
return os;
}
};

154
src/thermophysicalModels/specie/reaction/specieExponent/specieExponent.H Normal file
View File

@ -0,0 +1,154 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2021 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/>.
Class
Foam::specieExponent
Description
SourceFiles
specieExponentI.H
specieExponent.C
specieExponentIO.C
\*---------------------------------------------------------------------------*/
#ifndef specieExponent_H
#define specieExponent_H
#include "label.H"
#include "scalar.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
class Ostream;
// Forward declaration of friend functions and operators
class specieExponent;
scalar pow(scalar x, const specieExponent& e);
specieExponent operator+(const specieExponent& a, const specieExponent& b);
specieExponent operator-(const specieExponent& a, const specieExponent& b);
Ostream& operator<<(Ostream&, const specieExponent&);
/*---------------------------------------------------------------------------*\
Class specieExponent Declaration
\*---------------------------------------------------------------------------*/
class specieExponent
{
// Private Data
//- Marker value to indicate that there is no possible integer
// representation of this exponent
static const label noIntegerExponent_ = labelMax;
//- Integer exponent
label integerExponent_;
//- Scalar exponent
scalar scalarExponent_;
// Private Member Functions
//- Return whether or not this exponent has an integer representation
inline bool hasIntegerExponent() const;
public:
// Constructors
//- Construct null
inline specieExponent();
//- Construct from integer
inline specieExponent(const label integerExponent);
//- Construct from scalar
inline specieExponent(const scalar scalarExponent);
// Member Operators
//- Cast to scalar
inline operator scalar() const;
//- Assign to integer
inline specieExponent& operator=(const label integerExponent);
//- Assign to scalar
inline specieExponent& operator=(const scalar scalarExponent);
//- Negate a specie exponent
inline specieExponent operator-() const;
// Friend Functions
//- Compute the power of a number to a specie exponent
inline friend scalar pow(const scalar x, const specieExponent& e);
// Friend Operators
//- Sum two specie exponents
inline friend specieExponent operator+
(
const specieExponent& a,
const specieExponent& b
);
//- Subtract two specie exponents
inline friend specieExponent operator-
(
const specieExponent& a,
const specieExponent& b
);
// IOstream Operators
//- Write to output stream
inline friend Ostream& operator<<(Ostream& os, const specieExponent& e);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "specieExponentI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

174
src/thermophysicalModels/specie/reaction/specieExponent/specieExponentI.H Normal file
View File

@ -0,0 +1,174 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2021 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/>.
\*---------------------------------------------------------------------------*/
#include "specieExponent.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
inline bool Foam::specieExponent::hasIntegerExponent() const
{
return integerExponent_ != noIntegerExponent_;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
inline Foam::specieExponent::specieExponent()
:
integerExponent_(noIntegerExponent_),
scalarExponent_(NaN)
{}
inline Foam::specieExponent::specieExponent(const label integerExponent)
:
integerExponent_(integerExponent),
scalarExponent_(integerExponent)
{}
inline Foam::specieExponent::specieExponent(const scalar scalarExponent)
:
integerExponent_(labelMax),
scalarExponent_(scalarExponent)
{
const label integerExponent = floor(scalarExponent);
if (integerExponent == scalarExponent)
{
integerExponent_ = integerExponent;
}
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
inline Foam::specieExponent::operator scalar() const
{
return scalarExponent_;
}
inline Foam::specieExponent&
Foam::specieExponent::operator=(const label integerExponent)
{
return *this = specieExponent(integerExponent);
}
inline Foam::specieExponent&
Foam::specieExponent::operator=(const scalar scalarExponent)
{
return *this = specieExponent(scalarExponent);
}
inline Foam::specieExponent Foam::specieExponent::operator-() const
{
if (hasIntegerExponent())
{
return -integerExponent_;
}
else
{
return -scalarExponent_;
}
}
// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
inline Foam::scalar Foam::pow(const scalar x, const specieExponent& e)
{
if (e.hasIntegerExponent())
{
if (e.integerExponent_ == 0)
{
return 1;
}
else
{
scalar xx = e.integerExponent_ > 0 ? x : 1/x;
scalar y = 1;
for (label i = e.integerExponent_; i != 0; i /= 2)
{
if (i % 2)
{
y *= xx;
}
xx *= xx;
}
return y;
}
}
else
{
return pow(x, e.scalarExponent_);
}
}
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
inline Foam::specieExponent Foam::operator+
(
const specieExponent& a,
const specieExponent& b
)
{
if (a.hasIntegerExponent() && b.hasIntegerExponent())
{
return a.integerExponent_ + b.integerExponent_;
}
else
{
return a.scalarExponent_ + b.scalarExponent_;
}
}
inline Foam::specieExponent Foam::operator-
(
const specieExponent& a,
const specieExponent& b
)
{
return a + (-b);
}
// * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * * //
inline Foam::Ostream& Foam::operator<<(Ostream& os, const specieExponent& e)
{
os << e.operator scalar();
return os;
}
// ************************************************************************* //