ENH: Updated cloud coupling terms

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H

@@ -518,7 +518,7 @@ inline Foam::scalar Foam::KinematicParcel<ParcelType>::Re
     const scalar muc
 ) const
 {
-    return rhoc*mag(U - Uc_)*d/muc;
+    return rhoc*mag(U - Uc_)*d/(muc + ROOTVSMALL);
 }
 
 

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

@@ -63,7 +63,7 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::CpEff
 
 template<class ParcelType>
 template<class TrackData>
-Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::HEff
+Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::HsEff
 (
     TrackData& td,
     const scalar p,
@@ -74,9 +74,9 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::HEff
 ) const
 {
     return
-        this->Y_[GAS]*td.cloud().composition().H(idG, YGas_, p, T)
+        this->Y_[GAS]*td.cloud().composition().Hs(idG, YGas_, p, T)
-      + this->Y_[LIQ]*td.cloud().composition().H(idL, YLiquid_, p, T)
+      + this->Y_[LIQ]*td.cloud().composition().Hs(idL, YLiquid_, p, T)
-      + this->Y_[SLD]*td.cloud().composition().H(idS, YSolid_, p, T);
+      + this->Y_[SLD]*td.cloud().composition().Hs(idS, YSolid_, p, T);
 }
 
 
@@ -326,7 +326,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         updateMassFractions(mass0, dMassGas, dMassLiquid, dMassSolid);
 
 
-
     // Heat transfer
     // ~~~~~~~~~~~~~
 
@@ -383,25 +382,37 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         // Transfer mass lost from particle to carrier mass source
         forAll(YGas_, i)
         {
+            scalar dm = np0*dMassGas[i];
             label gid = composition.localToGlobalCarrierId(GAS, i);
-            td.cloud().rhoTrans(gid)[cellI] += np0*dMassGas[i];
+            scalar hs = composition.carrier().Hs(gid, 0.5*(T0 + T1));
+            td.cloud().rhoTrans(gid)[cellI] += dm;
+            td.cloud().hsTrans()[cellI] += dm*hs;
         }
         forAll(YLiquid_, i)
         {
+            scalar dm = np0*dMassLiquid[i];
             label gid = composition.localToGlobalCarrierId(LIQ, i);
-            td.cloud().rhoTrans(gid)[cellI] += np0*dMassLiquid[i];
+            scalar hs = composition.carrier().Hs(gid, 0.5*(T0 + T1));
+            td.cloud().rhoTrans(gid)[cellI] += dm;
+            td.cloud().hsTrans()[cellI] += dm*hs;
         }
 /*
         // No mapping between solid components and carrier phase
         forAll(YSolid_, i)
         {
+            scalar dm = np0*dMassSolid[i];
             label gid = composition.localToGlobalCarrierId(SLD, i);
-            td.cloud().rhoTrans(gid)[cellI] += np0*dMassSolid[i];
+            scalar hs = composition.carrier().Hs(gid, 0.5*(T0 + T1));
+            td.cloud().rhoTrans(gid)[cellI] += dm;
+            td.cloud().hsTrans()[cellI] += dm*hs;
         }
 */
         forAll(dMassSRCarrier, i)
         {
-            td.cloud().rhoTrans(i)[cellI] += np0*dMassSRCarrier[i];
+            scalar dm = np0*dMassSRCarrier[i];
+            scalar hs = composition.carrier().Hs(i, 0.5*(T0 + T1));
+            td.cloud().rhoTrans(i)[cellI] += dm;
+            td.cloud().hsTrans()[cellI] += dm*hs;
         }
 
         // Update momentum transfer
@@ -421,36 +432,38 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     // Remove the particle when mass falls below minimum threshold
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    if (mass1 < td.cloud().constProps().minParticleMass())
+    if (np0*mass1 < td.cloud().constProps().minParticleMass())
     {
         td.keepParticle = false;
 
         if (td.cloud().solution().coupled())
         {
+            scalar dm = np0*mass1;
+
             // Absorb parcel into carrier phase
             forAll(YGas_, i)
             {
                 label gid = composition.localToGlobalCarrierId(GAS, i);
-                td.cloud().rhoTrans(gid)[cellI] += np0*mass1*YMix[GAS]*YGas_[i];
+                td.cloud().rhoTrans(gid)[cellI] += dm*YMix[GAS]*YGas_[i];
             }
             forAll(YLiquid_, i)
             {
                 label gid = composition.localToGlobalCarrierId(LIQ, i);
-                td.cloud().rhoTrans(gid)[cellI] +=
+                td.cloud().rhoTrans(gid)[cellI] += dm*YMix[LIQ]*YLiquid_[i];
-                    np0*mass1*YMix[LIQ]*YLiquid_[i];
             }
 /*
             // No mapping between solid components and carrier phase
             forAll(YSolid_, i)
             {
                 label gid = composition.localToGlobalCarrierId(SLD, i);
-                td.cloud().rhoTrans(gid)[cellI] +=
+                td.cloud().rhoTrans(gid)[cellI] += dm*YMix[SLD]*YSolid_[i];
-                    np0*mass1*YMix[SLD]*YSolid_[i];
             }
 */
-            td.cloud().UTrans()[cellI] += np0*mass1*U1;
+            td.cloud().UTrans()[cellI] += dm*U1;
 
-            // enthalpy transfer accounted for via change in mass fractions
+            td.cloud().hsTrans()[cellI] += dm*HsEff(td, pc, T1, idG, idL, idS);
+
+            td.cloud().addToMassPhaseChange(dm);
         }
     }
 
@@ -531,28 +544,35 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
 
     Sh -= dMassTot*td.cloud().constProps().LDevol()/dt;
 
-    // Molar average molecular weight of carrier mix
-    const scalar Wc = this->rhoc_*specie::RR*this->Tc_/this->pc_;
-
     // Update molar emissions
-    forAll(dMassDV, i)
+    if (td.cloud().heatTransfer().BirdCorrection())
     {
+        // Molar average molecular weight of carrier mix
+        const scalar Wc =
+            max(SMALL, this->rhoc_*specie::RR*this->Tc_/this->pc_);
+
         // Note: hardcoded gaseous diffusivities for now
         // TODO: add to carrier thermo
         const scalar beta = sqr(cbrt(15.0) + cbrt(15.0));
-        const label id = composition.localToGlobalCarrierId(GAS, i);
-        const scalar Cp = composition.carrier().Cp(id, Ts);
-        const scalar W = composition.carrier().W(id);
-        const scalar Ni = dMassDV[i]/(this->areaS(d)*dt*W);
 
-        // Dab calc'd using API vapour mass diffusivity function
+        forAll(dMassDV, i)
-        const scalar Dab =
+        {
-            3.6059e-3*(pow(1.8*Ts, 1.75))*sqrt(1.0/W + 1.0/Wc)/(this->pc_*beta);
+            const label id = composition.localToGlobalCarrierId(GAS, i);
+            const scalar Cp = composition.carrier().Cp(id, Ts);
+            const scalar W = composition.carrier().W(id);
+            const scalar Ni = dMassDV[i]/(this->areaS(d)*dt*W);
 
-        N += Ni;
+            // Dab calc'd using API vapour mass diffusivity function
-        NCpW += Ni*Cp*W;
+            const scalar Dab =
-        Cs[id] += Ni*d/(2.0*Dab);
+                3.6059e-3*(pow(1.8*Ts, 1.75))
-     }
+               *sqrt(1.0/W + 1.0/Wc)
+               /(this->pc_*beta);
+
+            N += Ni;
+            NCpW += Ni*Cp*W;
+            Cs[id] += Ni*d/(2.0*Dab);
+        }
+    }
 }
 
 

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

@@ -134,9 +134,9 @@ private:
             const label idS
         ) const;
 
-        //- Return the mixture effective enthalpy
+        //- Return the mixture effective sensible enthalpy
         template<class TrackData>
-        scalar HEff
+        scalar HsEff
         (
             TrackData& td,
             const scalar p,

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

@@ -213,9 +213,17 @@ void Foam::ReactingParcel<ParcelType>::correctSurfaceValues
         sumYiCbrtW += Ys[i]*cbrtW;
     }
 
+    Cps = max(Cps, ROOTVSMALL);
+
     rhos *= pc_/(specie::RR*T);
+    rhos = max(rhos, ROOTVSMALL);
+
     mus /= sumYiSqrtW;
+    mus = max(mus, ROOTVSMALL);
+
     kappas /= sumYiCbrtW;
+    kappas = max(kappas, ROOTVSMALL);
+
     Prs = Cps*mus/kappas;
 }
 
@@ -335,7 +343,9 @@ void Foam::ReactingParcel<ParcelType>::calc
     Res = this->Re(U0, d0, rhos, mus);
 
     // Update particle component mass and mass fractions
-    scalar mass1 = updateMassFraction(mass0, dMassPC, Y_);
+    scalarField dMass(dMassPC);
+
+    scalar mass1 = updateMassFraction(mass0, dMass, Y_);
 
 
     // Heat transfer
@@ -390,11 +400,15 @@ void Foam::ReactingParcel<ParcelType>::calc
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     if (td.cloud().solution().coupled())
     {
-        // Transfer mass lost from particle to carrier mass source
+        // Transfer mass lost to carrier mass and enthalpy sources
-        forAll(dMassPC, i)
+        forAll(dMass, i)
         {
+            scalar dm = np0*dMass[i];
             label gid = composition.localToGlobalCarrierId(0, i);
-            td.cloud().rhoTrans(gid)[cellI] += np0*dMassPC[i];
+            scalar hs = composition.carrier().Hs(gid, 0.5*(T0 + T1));
+
+            td.cloud().rhoTrans(gid)[cellI] += dm;
+            td.cloud().hsTrans()[cellI] += dm*hs;
         }
 
         // Update momentum transfer
@@ -413,21 +427,27 @@ void Foam::ReactingParcel<ParcelType>::calc
 
     // Remove the particle when mass falls below minimum threshold
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    if (mass1 < td.cloud().constProps().minParticleMass())
+    if (np0*mass1 < td.cloud().constProps().minParticleMass())
     {
         td.keepParticle = false;
 
         if (td.cloud().solution().coupled())
         {
+            scalar dm = np0*mass1;
+
             // Absorb parcel into carrier phase
             forAll(Y_, i)
             {
+                scalar dmi = dm*Y_[i];
                 label gid = composition.localToGlobalCarrierId(0, i);
-                td.cloud().rhoTrans(gid)[cellI] += np0*mass1*Y_[i];
+                scalar hs = composition.carrier().Hs(gid, T1);
-            }
-            td.cloud().UTrans()[cellI] += np0*mass1*U1;
 
-            // enthalpy transfer accounted for via change in mass fractions
+                td.cloud().rhoTrans(gid)[cellI] += dmi;
+                td.cloud().hsTrans()[cellI] += dmi*hs;
+            }
+            td.cloud().UTrans()[cellI] += dm*U1;
+
+            td.cloud().addToMassPhaseChange(dm);
         }
     }
 
@@ -514,33 +534,44 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
     // Add to cumulative phase change mass
     td.cloud().addToMassPhaseChange(this->nParticle_*dMassTot);
 
-    // Average molecular weight of carrier mix - assumes perfect gas
+    forAll(dMassPC, i)
-    const scalar Wc = this->rhoc_*specie::RR*this->Tc_/this->pc_;
-
-    forAll(YComponents, i)
     {
         const label idc = composition.localToGlobalCarrierId(idPhase, i);
         const label idl = composition.globalIds(idPhase)[i];
 
         const scalar dh = td.cloud().phaseChange().dh(idc, idl, pc_, T);
         Sh -= dMassPC[i]*dh/dt;
+    }
 
-        // Update particle surface thermo properties
-        const scalar Dab =
-            composition.liquids().properties()[idl].D(pc_, Ts, Wc);
 
-        const scalar Cp = composition.carrier().Cp(idc, Ts);
+    // Update molar emissions
-        const scalar W = composition.carrier().W(idc);
+    if (td.cloud().heatTransfer().BirdCorrection())
-        const scalar Ni = dMassPC[i]/(this->areaS(d)*dt*W);
+    {
+        // Average molecular weight of carrier mix - assumes perfect gas
+        const scalar Wc = this->rhoc_*specie::RR*this->Tc_/this->pc_;
 
-        // Molar flux of species coming from the particle (kmol/m^2/s)
-        N += Ni;
 
-        // Sum of Ni*Cpi*Wi of emission species
+        forAll(dMassPC, i)
-        NCpW += Ni*Cp*W;
+        {
+            const label idc = composition.localToGlobalCarrierId(idPhase, i);
+            const label idl = composition.globalIds(idPhase)[i];
 
-        // Concentrations of emission species
+            const scalar Cp = composition.carrier().Cp(idc, Ts);
-        Cs[idc] += Ni*d/(2.0*Dab);
+            const scalar W = composition.carrier().W(idc);
+            const scalar Ni = dMassPC[i]/(this->areaS(d)*dt*W);
+
+            const scalar Dab =
+                composition.liquids().properties()[idl].D(pc_, Ts, Wc);
+
+            // Molar flux of species coming from the particle (kmol/m^2/s)
+            N += Ni;
+
+            // Sum of Ni*Cpi*Wi of emission species
+            NCpW += Ni*Cp*W;
+
+            // Concentrations of emission species
+            Cs[idc] += Ni*d/(2.0*Dab);
+        }
     }
 }
 

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

@@ -153,7 +153,10 @@ void Foam::ThermoParcel<ParcelType>::calcSurfaceValues
     mus = td.muInterp().interpolate(this->position(), tetIs)/TRatio;
 
     Pr = td.cloud().constProps().Pr();
+    Pr = max(ROOTVSMALL, Pr);
+
     kappas = Cpc_*mus/Pr;
+    kappas = max(ROOTVSMALL, kappas);
 }