ENH: Updated coal surface reaction models

src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationDiffusionLimitedRate/COxidationDiffusionLimitedRate.C

@@ -26,6 +26,8 @@ License
 #include "COxidationDiffusionLimitedRate.H"
 #include "mathematicalConstants.H"
 
+using namespace Foam::constant;
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -130,12 +132,13 @@ Foam::scalar Foam::COxidationDiffusionLimitedRate<CloudType>::calculate
         return 0.0;
     }
 
+    const SLGThermo& thermo = this->owner().thermo();
+
     // Local mass fraction of O2 in the carrier phase
-    const scalar YO2 = this->owner().thermo().carrier().Y(O2GlobalId_)[cellI];
+    const scalar YO2 = thermo.carrier().Y(O2GlobalId_)[cellI];
 
     // Change in C mass [kg]
-    scalar dmC =
-        4.0*constant::mathematical::pi*d*D_*YO2*Tc*rhoc/(Sb_*(T + Tc))*dt;
+    scalar dmC = 4.0*mathematical::pi*d*D_*YO2*Tc*rhoc/(Sb_*(T + Tc))*dt;
 
     // Limit mass transfer by availability of C
     dmC = min(mass*fComb, dmC);
@@ -153,8 +156,12 @@ Foam::scalar Foam::COxidationDiffusionLimitedRate<CloudType>::calculate
     dMassSRCarrier[O2GlobalId_] -= dmO2;
     dMassSRCarrier[CO2GlobalId_] += dmCO2;
 
+    const scalar HC = thermo.solids().properties()[CsLocalId_].H(T);
+    const scalar HCO2 = thermo.carrier().H(CO2GlobalId_, T);
+    const scalar HO2 = thermo.carrier().H(O2GlobalId_, T);
+
     // Heat of reaction [J]
-    return -HcCO2_*dmCO2;
+    return dmC*HC + dmO2*HO2 - dmCO2*HCO2;
 }
 
 

src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationKineticDiffusionLimitedRate/COxidationKineticDiffusionLimitedRate.C

@@ -42,7 +42,6 @@ COxidationKineticDiffusionLimitedRate
         owner,
         typeName
     ),
-    Sb_(dimensionedScalar(this->coeffDict().lookup("Sb")).value()),
     C1_(dimensionedScalar(this->coeffDict().lookup("C1")).value()),
     C2_(dimensionedScalar(this->coeffDict().lookup("C2")).value()),
     E_(dimensionedScalar(this->coeffDict().lookup("E")).value()),
@@ -62,19 +61,6 @@ COxidationKineticDiffusionLimitedRate
     const scalar WCO2 = owner.thermo().carrier().W(CO2GlobalId_);
     WC_ = WCO2 - WO2_;
     HcCO2_ = owner.thermo().carrier().Hc(CO2GlobalId_);
-
-    if (Sb_ < 0)
-    {
-        FatalErrorIn
-        (
-            "COxidationKineticDiffusionLimitedRate"
-            "("
-                "const dictionary&, "
-                "CloudType&"
-            ")"
-        )   << "Stoichiometry of reaction, Sb, must be greater than zero" << nl
-            << exit(FatalError);
-    }
 }
 
 
@@ -86,7 +72,6 @@ COxidationKineticDiffusionLimitedRate
 )
 :
     SurfaceReactionModel<CloudType>(srm),
-    Sb_(srm.Sb_),
     C1_(srm.C1_),
     C2_(srm.C2_),
     E_(srm.E_),
@@ -140,8 +125,10 @@ Foam::scalar Foam::COxidationKineticDiffusionLimitedRate<CloudType>::calculate
         return 0.0;
     }
 
+    const SLGThermo& thermo = this->owner().thermo();
+
     // Local mass fraction of O2 in the carrier phase
-    const scalar YO2 = this->owner().thermo().carrier().Y(O2GlobalId_)[cellI];
+    const scalar YO2 = thermo.carrier().Y(O2GlobalId_)[cellI];
 
     // Diffusion rate coefficient
     const scalar D0 = C1_/d*pow(0.5*(T + Tc), 0.75);
@@ -158,21 +145,28 @@ Foam::scalar Foam::COxidationKineticDiffusionLimitedRate<CloudType>::calculate
     // Limit mass transfer by availability of C
     dmC = min(mass*fComb, dmC);
 
+    // Molar consumption
+    const scalar dOmega = dmC/WC_;
+
     // Change in O2 mass [kg]
-    const scalar dmO2 = dmC/WC_*Sb_*WO2_;
+    const scalar dmO2 = dOmega*WO2_;
 
     // Mass of newly created CO2 [kg]
-    const scalar dmCO2 = dmC + dmO2;
+    const scalar dmCO2 = dOmega*(WC_ + WO2_);
 
     // Update local particle C mass
-    dMassSolid[CsLocalId_] += dmC;
+    dMassSolid[CsLocalId_] += dOmega*WC_;
 
     // Update carrier O2 and CO2 mass
     dMassSRCarrier[O2GlobalId_] -= dmO2;
     dMassSRCarrier[CO2GlobalId_] += dmCO2;
 
+    const scalar HC = thermo.solids().properties()[CsLocalId_].H(T);
+    const scalar HCO2 = thermo.carrier().H(CO2GlobalId_, T);
+    const scalar HO2 = thermo.carrier().H(O2GlobalId_, T);
+
     // Heat of reaction [J]
-    return -HcCO2_*dmCO2;
+    return dOmega*(WC_*HC + WO2_*HO2 - (WC_ + WO2_)*HCO2);
 }
 
 

src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationKineticDiffusionLimitedRate/COxidationKineticDiffusionLimitedRate.H

@@ -61,9 +61,6 @@ class COxidationKineticDiffusionLimitedRate
 
         // Model constants
 
-            //- Stoichiometry of reaction
-            const scalar Sb_;
-
             //- Mass diffusion limited rate constant, C1
             const scalar C1_;
 

src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationMurphyShaddix/COxidationMurphyShaddix.C

@@ -139,9 +139,10 @@ Foam::scalar Foam::COxidationMurphyShaddix<CloudType>::calculate
         return 0.0;
     }
 
+    const SLGThermo& thermo = this->owner().thermo();
+
     // Cell carrier phase O2 species density [kg/m^3]
-    const scalar rhoO2 =
-        rhoc*this->owner().thermo().carrier().Y(O2GlobalId_)[cellI];
+    const scalar rhoO2 = rhoc*thermo.carrier().Y(O2GlobalId_)[cellI];
 
     if (rhoO2 < SMALL)
     {
@@ -220,11 +221,9 @@ Foam::scalar Foam::COxidationMurphyShaddix<CloudType>::calculate
     // Add to particle mass transfer
     dMassSolid[CsLocalId_] += dOmega*WC_;
 
-    const scalar HC =
-        this->owner().composition().solids().properties()[CsLocalId_].Hf()
-      + this->owner().composition().solids().properties()[CsLocalId_].Cp()*T;
-    const scalar HCO2 = this->owner().thermo().carrier().H(CO2GlobalId_, T);
-    const scalar HO2 = this->owner().thermo().carrier().H(O2GlobalId_, T);
+    const scalar HC = thermo.solids().properties()[CsLocalId_].H(T);
+    const scalar HCO2 = thermo.carrier().H(CO2GlobalId_, T);
+    const scalar HO2 = thermo.carrier().H(O2GlobalId_, T);
 
     // Heat of reaction
     return dOmega*(WC_*HC + WO2_*HO2 - (WC_ + WO2_)*HCO2);