updates to the reacting parcel/cloud

src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C

@@ -59,7 +59,6 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
     const scalar mass0 = this->mass();
     const scalar np0 = this->nParticle_;
     const scalar T0 = this->T_;
-    const scalar cp0 = this->cp_;
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~
     // Initialise transfer terms
@@ -79,9 +78,12 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
     // - components do not necessarily exist in particle already
     scalarList dMassSR(td.cloud().gases().size(), 0.0);
 
-    // Total mass lost from particle due to surface reactions
+
-    // - sub-model will adjust component mass fractions
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    scalar dMassMTSR = 0.0;
+    // Calculate velocity - update U
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    scalar Cud = 0.0;
+    const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -91,13 +93,6 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
     scalar T1 = calcHeatTransfer(td, dt, celli, htc, dhTrans);
 
 
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Calculate velocity - update U
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    scalar Cud = 0.0;
-    const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
-
-
     // ~~~~~~~~~~~~~~~~~~~~~~~
     // Calculate mass transfer
     // ~~~~~~~~~~~~~~~~~~~~~~~
@@ -107,27 +102,36 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Calculate surface reactions
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    calcSurfaceReactions(td, dt, celli, T0, T1, dMassMTSR, dMassSR);
+
+    // Initialise enthalpy retention to zero
+    scalar dhRet = 0.0;
+
+    calcSurfaceReactions(td, dt, celli, T0, T1, dMassMT, dMassSR, dhRet);
 
     // New total mass
-    const scalar mass1 = mass0 - sum(dMassMT) - dMassMTSR;
+    const scalar mass1 = mass0 - sum(dMassMT) - sum(dMassSR);
 
-    // Ratio of mass devolatilised to the total volatile mass of the particle
+    // Correct particle temperature to account for latent heat of
-    const scalar fVol = 1 -
+    // devolatilisation
-        (YMixture_[0]*mass1)
+    T1 -=
-       /(td.cloud().composition().YMixture0()[0]*mass0_);
+        td.constProps().Ldevol()
+       *sum(dMassMT)
+       /(0.5*(mass0 + mass1)*this->cp_);
 
-    // Specific heat capacity of non-volatile components
+    // Add retained enthalpy from surface reaction to particle and remove
-    const scalar cpNonVolatile =
+    // from gas
-        (
+    T1 += dhRet/(0.5*(mass0 + mass1)*this->cp_);
-            YMixture_[1]*td.cloud().composition().cpLiquid(YLiquid_, pc_, this->Tc_)
+    dhTrans -= dhRet;
-          + YMixture_[2]*td.cloud().composition().cpSolid(YSolid_)
-        )/(YMixture_[1] + YMixture_[2]);
 
-    // New specific heat capacity - linear variation until volatiles
+    // Correct dhTrans to account for enthalpy of evolved volatiles
-    // have evolved
+    dhTrans +=
-    const scalar cp1 = (cpNonVolatile - td.constProps().cp0())*fVol
+        sum(dMassMT)
-       + td.constProps().cp0();
+       *td.cloud().composition().HGas(YGas_, 0.5*(T0 + T1));
+
+    // Correct dhTrans to account for enthalpy of consumed solids
+    dhTrans +=
+        sum(dMassSR)
+       *td.cloud().composition().HSolid(YSolid_, 0.5*(T0 + T1));
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~
@@ -147,8 +151,8 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
     td.cloud().UCoeff()[celli] += np0*mass0*Cud;
 
     // Update enthalpy transfer
-//    td.cloud().hTrans()[celli] += np0*(mass0*cp0*T0 - mass1*cp1*T1);
+    // - enthalpy of lost solids already accounted for
-    td.cloud().hTrans()[celli] += np0*((mass0*cp0 - mass1*cp1)*T0 + dhTrans);
+    td.cloud().hTrans()[celli] += np0*dhTrans;
 
     // Accumulate coefficient to be applied in carrier phase enthalpy coupling
     td.cloud().hCoeff()[celli] += np0*htc*this->areaS();
@@ -166,7 +170,12 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
         {
             td.cloud().rhoTrans(i)[celli] += np0*dMassMT[i];
         }
-        td.cloud().hTrans()[celli] += np0*mass1*cp1*T1;
+        td.cloud().hTrans()[celli] +=
+            np0*mass1
+           *(
+                YMixture_[0]*td.cloud().composition().HGas(YGas_, T1)
+              + YMixture_[2]*td.cloud().composition().HSolid(YSolid_, T1)
+            );
         td.cloud().UTrans()[celli] += np0*mass1*U1;
     }
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -176,7 +185,10 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
     {
         this->U_ = U1;
         this->T_ = T1;
-        this->cp_ = cp1;
+        this->cp_ =
+            YMixture_[0]*td.cloud().composition().cpGas(YGas_, T1)
+          + YMixture_[1]*td.cloud().composition().cpLiquid(YLiquid_, pc_, T1)
+          + YMixture_[2]*td.cloud().composition().cpSolid(YSolid_);
 
         // Update particle density or diameter
         if (td.cloud().massTransfer().changesVolume())
@@ -205,7 +217,7 @@ void Foam::ReactingParcel<ParcelType>::calcUncoupled
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     const scalar T0 = this->T_;
     const scalar mass0 = this->mass();
-//    const scalar cp0 = this->cp();
+    const scalar cp0 = this->cp_;
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~
     // Initialise transfer terms
@@ -225,9 +237,12 @@ void Foam::ReactingParcel<ParcelType>::calcUncoupled
     // - components do not necessarily exist in particle already
     scalarList dMassSR(td.cloud().gases().size(), 0.0);
 
-    // Total mass lost from particle due to surface reactions
+
-    // - sub-model will adjust component mass fractions
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    scalar dMassMTSR = 0.0;
+    // Calculate velocity - update U
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    scalar Cud = 0.0;
+    const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -240,13 +255,6 @@ void Foam::ReactingParcel<ParcelType>::calcUncoupled
     T1 = min(td.constProps().Tvap(), T1);
 
 
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Calculate velocity - update U
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    scalar Cud = 0.0;
-    const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
-
-
     // ~~~~~~~~~~~~~~~~~~~~~~~
     // Calculate mass transfer
     // ~~~~~~~~~~~~~~~~~~~~~~~
@@ -256,37 +264,23 @@ void Foam::ReactingParcel<ParcelType>::calcUncoupled
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Calculate surface reactions
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    calcSurfaceReactions
+
-    (
+    // Initialise enthalpy retention to zero
-        td,
+    scalar dhRet = 0.0;
-        dt,
+
-        celli,
+    calcSurfaceReactions(td, dt, celli, T0, T1, dMassMT, dMassSR, dhRet);
-        T0,
-        T1,
-        dMassMTSR,
-        dMassSR
-    );
 
     // New total mass
-    const scalar mass1 = mass0 - sum(dMassMT) - dMassMTSR;
+    const scalar mass1 = mass0 - sum(dMassMT) - sum(dMassSR);
 
-    // Ratio of mass devolatilised to the total volatile mass of the particle
+    // New specific heat capacity
-    const scalar fVol = 1 -
+    const scalar cp1 =
-        (YMixture_[0]*mass1)
+        YMixture_[0]*td.cloud().composition().cpGas(YGas_, T1)
-       /(td.cloud().composition().YMixture0()[0]*mass0_);
+      + YMixture_[1]*td.cloud().composition().cpLiquid(YLiquid_, pc_, T1)
-
+      + YMixture_[2]*td.cloud().composition().cpSolid(YSolid_);
-    // Specific heat capacity of non-volatile components
-    const scalar cpNonVolatile =
-        (
-            YMixture_[1]*td.cloud().composition().cpLiquid(YLiquid_, pc_, this->Tc_)
-          + YMixture_[2]*td.cloud().composition().cpSolid(YSolid_)
-        )/(YMixture_[1] + YMixture_[2]);
-
-    // New specific heat capacity - linear variation until volatiles
-    // have evolved
-    const scalar cp1 = (cpNonVolatile - td.constProps().cp0())*fVol
-       + td.constProps().cp0();
 
+    // Add retained enthalpy to particle
+    T1 += dhRet/(mass0*0.5*(cp0 + cp1));
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Remove the particle when mass falls below minimum threshold
@@ -389,8 +383,9 @@ void Foam::ReactingParcel<ParcelType>::calcSurfaceReactions
     const label celli,
     const scalar T0,
     const scalar T1,
-    scalar& dMassMTSR,
+    const scalarList& dMassMT,
-    scalarList& dMassMT
+    scalarList& dMassSR,
+    scalar& dhRet
 )
 {
     // Check that model is active
@@ -409,21 +404,20 @@ void Foam::ReactingParcel<ParcelType>::calcSurfaceReactions
         T0,
         T1,
         this->Tc_,
+        pc_,
         this->rhoc_,
         this->mass(),
+        dMassMT,
         YGas_,
         YLiquid_,
         YSolid_,
         YMixture_,
-        dMassMTSR,
+        dMassSR,
-        dMassMT
+        dhRet
     );
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * *  IOStream operators * * * * * * * * * * * //
 
 #include "ReactingParcelIO.C"

src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H

@@ -89,6 +89,9 @@ public:
             //- Boiling point [K]
             const scalar Tbp_;
 
+            //- Latent heat of devolatilisation [J/kg]
+            const scalar Ldevol_;
+
 
     public:
 
@@ -102,6 +105,9 @@ public:
 
             //- Return const access to the boiling point
             inline scalar Tbp() const;
+
+            //- Return const access to the latent heat of devolatilisation
+            inline scalar Ldevol() const;
     };
 
 
@@ -210,8 +216,9 @@ protected:
             const label celli,
             const scalar T0,
             const scalar T1,
-            scalar& dMassMTSR,
+            const scalarList& dMassMT,
-            scalarList& dMassMT
+            scalarList& dMassSR,
+            scalar& dhRet
         );
 
 

src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H

@@ -34,7 +34,8 @@ inline Foam::ReactingParcel<ParcelType>::constantProperties::constantProperties
 :
     ThermoParcel<ParcelType>::constantProperties(dict),
     Tvap_(dimensionedScalar(dict.lookup("Tvap")).value()),
-    Tbp_(dimensionedScalar(dict.lookup("Tbp")).value())
+    Tbp_(dimensionedScalar(dict.lookup("Tbp")).value()),
+    Ldevol_(dimensionedScalar(dict.lookup("Ldevol")).value())
 {}
 
 
@@ -127,6 +128,14 @@ Foam::ReactingParcel<ParcelType>::constantProperties::Tbp() const
 }
 
 
+template<class ParcelType>
+inline Foam::scalar
+Foam::ReactingParcel<ParcelType>::constantProperties::Ldevol() const
+{
+    return Ldevol_;
+}
+
+
 // * * * * * * * * * * * trackData Member Functions  * * * * * * * * * * * * //
 
 template<class ParcelType>

src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H

@@ -209,13 +209,20 @@ public:
             //- Return the gas constant for the gas mixture
             virtual scalar RGas(const scalarField& YGas) const = 0;
 
-            //- Return enthalpy for the gas mixture
+            //- Return enthalpy for the gas mixture [energy per unit mass]
             virtual scalar HGas
             (
                  const scalarField& YGas,
                  const scalar T
             ) const = 0;
 
+            //- Return enthalpy for the solid mixture [energy per unit mass]
+            virtual scalar HSolid
+            (
+                const scalarField& YSolid,
+                const scalar T
+            ) const = 0;
+
             //- Return specific heat caparcity for the gas mixture
             virtual scalar cpGas
             (

src/lagrangian/intermediate/submodels/Reacting/CompositionModel/SingleMixtureFraction/SingleMixtureFraction.C

@@ -519,6 +519,29 @@ const
 }
 
 
+template<class CloudType>
+Foam::scalar Foam::SingleMixtureFraction<CloudType>::HSolid
+(
+    const scalarField& YSolid,
+    const scalar T
+)
+const
+{
+    scalar HMixture = 0.0;
+    forAll(YSolid, i)
+    {
+        label id = solidGlobalIds_[i];
+        HMixture +=
+            YSolid[i]
+           *(
+                this->solids().properties()[id].Hf()
+              + this->solids().properties()[id].cp()*T
+            );
+    }
+    return HMixture;
+}
+
+
 template<class CloudType>
 Foam::scalar Foam::SingleMixtureFraction<CloudType>::cpGas
 (

src/lagrangian/intermediate/submodels/Reacting/CompositionModel/SingleMixtureFraction/SingleMixtureFraction.H

@@ -188,9 +188,12 @@ public:
             //- Return the gas constant for the gas mixture
             scalar RGas(const scalarField& YGas) const;
 
-            //- Return enthalpy for the gas mixture
+            //- Return enthalpy for the gas mixture [energy per unit mass]
             scalar HGas(const scalarField& YGas, const scalar T) const;
 
+            //- Return enthalpy for the solid mixture [energy per unit mass]
+            scalar HSolid(const scalarField& YSolid, const scalar T) const;
+
             //- Return specific heat caparcity for the gas mixture
             scalar cpGas(const scalarField& YGas, const scalar T) const;
 

src/lagrangian/intermediate/submodels/Reacting/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.C

@@ -64,16 +64,20 @@ void Foam::NoSurfaceReaction<CloudType>::calculate
     const scalar T0,
     const scalar T1,
     const scalar Tc,
+    const scalar pc,
     const scalar rhoc,
     const scalar massp,
+    const scalarList& dMassMT,
     scalarField& YGas,
     scalarField& YLiquid,
     scalarField& YSolid,
     scalarField& YMixture,
-    scalar& dMassMTSR,
+    scalarList& dMassSR,
-    scalarList& dMassSR
+    scalar& dhRet
 ) const
-{}
+{
+    // do nothing
+}
 
 
 // ************************************************************************* //

src/lagrangian/intermediate/submodels/Reacting/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.H

@@ -84,14 +84,16 @@ public:
             const scalar T0,
             const scalar T1,
             const scalar Tc,
+            const scalar pc,
             const scalar rhoc,
             const scalar massp,
+            const scalarList& dMassMT,
             scalarField& YGas,
             scalarField& YLiquid,
             scalarField& YSolid,
             scalarField& YMixture,
-            scalar& dMassMTSR,
+            scalarList& dMassSR,
-            scalarList& dMassSR
+            scalar& dhRet
         ) const;
 };
 

src/lagrangian/intermediate/submodels/Reacting/SurfaceReactionModel/SurfaceReactionModel/SurfaceReactionModel.H

@@ -141,14 +141,16 @@ public:
             const scalar T0,
             const scalar T1,
             const scalar Tc,
+            const scalar pc,
             const scalar rhoc,
             const scalar massp,
+            const scalarList& dMassMT,
             scalarField& YGas,
             scalarField& YLiquid,
             scalarField& YSolid,
             scalarField& YMixture,
-            scalar& dMassMTSR,
+            scalarList& dMassSR,
-            scalarList& dMassSR
+            scalar& dhRet
         ) const = 0;
 };