31 template<
class CloudType>
34 template<
class CloudType>
40 template<
class CloudType>
62 O2GlobalId_(owner.composition().globalCarrierId(
"O2")),
63 CO2GlobalId_(owner.composition().globalCarrierId(
"CO2")),
68 label idSolid = owner.composition().idSolid();
69 CsLocalId_ = owner.composition().localId(idSolid,
"C");
72 WO2_ = owner.mcCarrierThermo().speciesData()[O2GlobalId_].W();
73 scalar WCO2 = owner.mcCarrierThermo().speciesData()[CO2GlobalId_].W();
80 template<
class CloudType>
87 template<
class CloudType>
94 template<
class CloudType>
117 const label idSolid = CloudType::parcelType::SLD;
118 const scalar fComb = YMixture[idSolid]*YSolid[CsLocalId_];
128 rhoc*this->owner().mcCarrierThermo().Y(O2GlobalId_)[cellI];
140 const scalar D = D0_*(rho0_/rhoc)*
pow(Tc/T0_, Dn_);
143 const scalar ppO2 = rhoO2/WO2_*specie::RR*Tc;
146 const scalar
C = pc/(specie::RR*Tc);
150 Pout<<
"mass = " << mass <<
nl
151 <<
"fComb = " << fComb <<
nl
152 <<
"Ap = " << Ap <<
nl
153 <<
"dt = " << dt <<
nl
162 const scalar qCsLim = mass*fComb/(WC_*Ap*dt);
166 Pout <<
"qCsLim = " << qCsLim <<
endl;
170 while ((
mag(qCs - qCsOld)/qCs > tolerance_) && (iter <= maxIters_))
173 const scalar PO2Surface = ppO2*
exp(-(qCs + N)*d/(2*C*D));
174 qCs = A_*
exp(-E_/(specie::RR*T))*
pow(PO2Surface, n_);
175 qCs = (100.0*qCs + iter*qCsOld)/(100.0 + iter);
176 qCs =
min(qCs, qCsLim);
180 Pout<<
"iter = " << iter
181 <<
", qCsOld = " << qCsOld
189 if (iter > maxIters_)
193 "scalar Foam::COxidationMurphyShaddix<CloudType>::calculate(...)"
194 ) <<
"iter limit reached (" << maxIters_ <<
")" <<
nl <<
endl;
198 scalar dOmega = qCs*Ap*dt;
201 dMassSRCarrier[O2GlobalId_] += -dOmega*WO2_;
202 dMassSRCarrier[CO2GlobalId_] += dOmega*(WC_ + WO2_);
205 dMassSolid[CsLocalId_] += dOmega*WC_;
208 this->owner().composition().solids().properties()[CsLocalId_].Hf()
209 + this->owner().composition().solids().properties()[CsLocalId_].cp()*
T;
211 this->owner().mcCarrierThermo().speciesData()[CO2GlobalId_].H(T);
213 this->owner().mcCarrierThermo().speciesData()[O2GlobalId_].H(T);
216 return dOmega*(WC_*HC + WO2_*HO2 - (WC_ + WO2_)*HCO2);