library updates + introduction of thermal shielding

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

@@ -103,7 +103,7 @@ Foam::scalar Foam::COxidationDiffusionLimitedRate<CloudType>::calculate
     const scalarField& YLiquid,
     const scalarField& YSolid,
     const scalarField& YMixture,
-    const scalarField& dMassVolatile,
+    const scalar N,
     scalarField& dMassGas,
     scalarField& dMassLiquid,
     scalarField& dMassSolid,

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

@@ -131,7 +131,7 @@ public:
             const scalarField& YLiquid,
             const scalarField& YSolid,
             const scalarField& YMixture,
-            const scalarField& dMassVolatile,
+            const scalar N,
             scalarField& dMassGas,
             scalarField& dMassLiquid,
             scalarField& dMassSolid,

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

@@ -111,7 +111,7 @@ Foam::scalar Foam::COxidationKineticDiffusionLimitedRate<CloudType>::calculate
     const scalarField& YLiquid,
     const scalarField& YSolid,
     const scalarField& YMixture,
-    const scalarField& dMassVolatile,
+    const scalar N,
     scalarField& dMassGas,
     scalarField& dMassLiquid,
     scalarField& dMassSolid,

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

@@ -139,7 +139,7 @@ public:
             const scalarField& YLiquid,
             const scalarField& YSolid,
             const scalarField& YMixture,
-            const scalarField& dMassVolatile,
+            const scalar N,
             scalarField& dMassGas,
             scalarField& dMassLiquid,
             scalarField& dMassSolid,

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

@@ -107,7 +107,7 @@ Foam::scalar Foam::COxidationMurphyShaddix<CloudType>::calculate
     const scalarField& YLiquid,
     const scalarField& YSolid,
     const scalarField& YMixture,
-    const scalarField& dMassVolatile,
+    const scalar N,
     scalarField& dMassGas,
     scalarField& dMassLiquid,
     scalarField& dMassSolid,
@@ -143,9 +143,6 @@ Foam::scalar Foam::COxidationMurphyShaddix<CloudType>::calculate
     // Far field partial pressure O2 [Pa]
     const scalar ppO2 = rhoO2/WO2_*specie::RR*Tc;
 
-    // Molar emission rate of volatiles per unit surface area
-    const scalar qVol = sum(dMassVolatile)/(WVol_*Ap);
-
     // Total molar concentration of the carrier phase [kmol/m^3]
     const scalar C = pc/(specie::RR*Tc);
 
@@ -174,7 +171,7 @@ Foam::scalar Foam::COxidationMurphyShaddix<CloudType>::calculate
     while ((mag(qCs - qCsOld)/qCs > tolerance_) && (iter <= maxIters_))
     {
         qCsOld = qCs;
-        const scalar PO2Surface = ppO2*exp(-(qCs + qVol)*d/(2*C*D));
+        const scalar PO2Surface = ppO2*exp(-(qCs + N)*d/(2*C*D));
         qCs = A_*exp(-E_/(specie::RR*T))*pow(PO2Surface, n_);
         qCs = (100.0*qCs + iter*qCsOld)/(100.0 + iter);
         qCs = min(qCs, qCsLim);

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

@@ -152,7 +152,7 @@ public:
             const scalarField& YLiquid,
             const scalarField& YSolid,
             const scalarField& YMixture,
-            const scalarField& dMassVolatile,
+            const scalar N,
             scalarField& dMassGas,
             scalarField& dMassLiquid,
             scalarField& dMassSolid,

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

@@ -101,6 +101,12 @@ void Foam::KinematicParcel<ParcelType>::calc
     const scalar rho0 = rho_;
     const scalar mass0 = mass();
 
+    // Reynolds number
+    const scalar Re = this->Re(U0, d0, muc_);
+
+
+    // Sources
+    //~~~~~~~~
 
     // Explicit momentum source for particle
     vector Su = vector::zero;
@@ -113,7 +119,8 @@ void Foam::KinematicParcel<ParcelType>::calc
     // ~~~~~~
 
     // Calculate new particle velocity
-    vector U1 = calcVelocity(td, dt, cellI, d0, U0, rho0, mass0, Su, dUTrans);
+    vector U1 =
+        calcVelocity(td, dt, cellI, Re, muc_, d0, U0, rho0, mass0, Su, dUTrans);
 
 
     // Accumulate carrier phase source terms
@@ -138,6 +145,8 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
     TrackData& td,
     const scalar dt,
     const label cellI,
+    const scalar Re,
+    const scalar mu,
     const scalar d,
     const vector& U,
     const scalar rho,
@@ -149,8 +158,7 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
     const polyMesh& mesh = this->cloud().pMesh();
 
     // Momentum transfer coefficient
-    const scalar utc =
+    const scalar utc = td.cloud().drag().utc(Re, d, mu) + ROOTVSMALL;
-        td.cloud().drag().utc(U - Uc_, d, rhoc_, muc_) + ROOTVSMALL;
 
     // Momentum source due to particle forces
     const vector FCoupled =

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

@@ -238,6 +238,8 @@ protected:
             TrackData& td,
             const scalar dt,           // timestep
             const label cellI,         // owner cell
+            const scalar Re,           // Reynolds number
+            const scalar mu,           // local carrier viscosity
             const scalar d,            // diameter
             const vector& U,           // velocity
             const scalar rho,          // density
@@ -387,6 +389,14 @@ public:
             //- Surface area for given diameter
             inline scalar areaS(const scalar d) const;
 
+            //- Reynolds number - particle properties input
+            inline scalar Re
+            (
+                const vector& U,
+                const scalar d,
+                const scalar mu
+            ) const;
+
 
         // Main calculation loop
 

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

@@ -385,4 +385,17 @@ Foam::KinematicParcel<ParcelType>::areaS(const scalar d) const
 }
 
 
+template <class ParcelType>
+inline Foam::scalar
+Foam::KinematicParcel<ParcelType>::Re
+(
+    const vector& U,
+    const scalar d,
+    const scalar mu
+) const
+{
+    return rhoc_*mag(U - Uc_)*d/mu;
+}
+
+
 // ************************************************************************* //

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

@@ -206,6 +206,20 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     const label idL = td.cloud().composition().idLiquid();
     const label idS = td.cloud().composition().idSolid();
 
+
+    // Calc surface values
+    // ~~~~~~~~~~~~~~~~~~~
+    scalar Ts, rhos, mus, Pr, kappa;
+    ThermoParcel<ParcelType>::
+        calcSurfaceValues(td, cellI, T0, Ts, rhos, mus, Pr, kappa);
+
+    // Reynolds number
+    scalar Re = this->Re(U0, d0, mus);
+
+
+    // Sources
+    //~~~~~~~~
+
     // Explicit momentum source for particle
     vector Su = vector::zero;
 
@@ -219,6 +233,45 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     scalar dhsTrans = 0.0;
 
 
+    // Phase change in liquid phase
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Mass transfer due to phase change
+    scalarField dMassPC(YLiquid_.size(), 0.0);
+
+    // Molar flux of species emitted from the particle (kmol/m^2/s)
+    scalar Ne = 0.0;
+
+    // Sum Ni*Cpi*Wi of emission species
+    scalar NCpW = 0.0;
+
+    // Surface concentrations of emitted species
+    scalarField Cs(td.cloud().mcCarrierThermo().species().size(), 0.0);
+
+    // Calc mass and enthalpy transfer due to phase change
+    calcPhaseChange
+    (
+        td,
+        dt,
+        cellI,
+        Re,
+        Ts,
+        mus/rhos,
+        d0,
+        T0,
+        mass0,
+        idL,
+        YMix[LIQ],
+        YLiquid_,
+        dMassPC,
+        Sh,
+        dhsTrans,
+        Ne,
+        NCpW,
+        Cs
+    );
+
+
     // Devolatilisation
     // ~~~~~~~~~~~~~~~~
 
@@ -230,6 +283,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     (
         td,
         dt,
+        Ts,
+        d0,
         T0,
         mass0,
         this->mass0_,
@@ -239,9 +294,15 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         canCombust_,
         dMassDV,
         Sh,
-        dhsTrans
+        dhsTrans,
+        Ne,
+        NCpW,
+        Cs
     );
 
+    // Correct surface values due to emitted species
+    correctSurfaceValues(td, cellI, Ts, Cs, rhos, mus, Pr, kappa);
+
 
     // Surface reactions
     // ~~~~~~~~~~~~~~~~~
@@ -267,7 +328,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         T0,
         mass0,
         canCombust_,
-        dMassDV,    // assuming d(mass) due to phase change is non-volatile
+        Ne,
         YMix,
         YGas_,
         YLiquid_,
@@ -281,31 +342,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     );
 
 
-    // Phase change in liquid phase
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Mass transfer due to phase change
-    scalarField dMassPC(YLiquid_.size(), 0.0);
-
-    // Calc mass and enthalpy transfer due to phase change
-    calcPhaseChange
-    (
-        td,
-        dt,
-        cellI,
-        d0,
-        T0,
-        U0,
-        mass0,
-        idL,
-        YMix[LIQ],
-        YLiquid_,
-        dMassPC,
-        Sh,
-        dhsTrans
-    );
-
-
     // Update component mass fractions
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -322,14 +358,30 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
 
     // Calculate new particle temperature
     scalar T1 =
-        calcHeatTransfer(td, dt, cellI, d0, U0, rho0, T0, cp0, Sh, dhsTrans);
+        calcHeatTransfer
+        (
+            td,
+            dt,
+            cellI,
+            Re,
+            Pr,
+            kappa,
+            d0,
+            rho0,
+            T0,
+            cp0,
+            NCpW,
+            Sh,
+            dhsTrans
+        );
 
 
     // Motion
     // ~~~~~~
 
     // Calculate new particle velocity
-    vector U1 = calcVelocity(td, dt, cellI, d0, U0, rho0, mass0, Su, dUTrans);
+    vector U1 =
+        calcVelocity(td, dt, cellI, Re, mus, d0, U0, rho0, mass0, Su, dUTrans);
 
     dUTrans += 0.5*(mass0 - mass1)*(U0 + U1);
 
@@ -455,6 +507,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
 (
     TrackData& td,
     const scalar dt,
+    const scalar Ts,
+    const scalar d,
     const scalar T,
     const scalar mass,
     const scalar mass0,
@@ -464,7 +518,10 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
     bool& canCombust,
     scalarField& dMassDV,
     scalar& Sh,
-    scalar& dhsTrans
+    scalar& dhsTrans,
+    scalar& N,
+    scalar& NCpW,
+    scalarField& Cs
 ) const
 {
     // Check that model is active, and that the parcel temperature is
@@ -496,6 +553,30 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
     td.cloud().addToMassDevolatilisation(this->nParticle_*dMassTot);
 
     Sh -= dMassTot*td.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)
+    {
+        // Note: hardcoded gaseous diffusivities for now
+        // TODO: add to carrier thermo
+        const scalar beta = sqr(cbrt(15.0) + cbrt(15.0));
+        const label id =
+            td.cloud().composition().localToGlobalCarrierId(GAS, i);
+        const scalar Cp = td.cloud().mcCarrierThermo().speciesData()[id].Cp(Ts);
+        const scalar W = td.cloud().mcCarrierThermo().speciesData()[id].W();
+        const scalar Ni = dMassDV[i]/(this->areaS(d)*dt*W);
+
+        // Dab calc'd using API vapour mass diffusivity function
+        const scalar 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);
+     }
 }
 
 
@@ -510,7 +591,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
     const scalar T,
     const scalar mass,
     const bool canCombust,
-    const scalarField& dMassVolatile,
+    const scalar N,
     const scalarField& YMix,
     const scalarField& YGas,
     const scalarField& YLiquid,
@@ -544,7 +625,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
         YLiquid,
         YSolid,
         YMix,
-        dMassVolatile,
+        N,
         dMassSRGas,
         dMassSRLiquid,
         dMassSRSolid,

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

@@ -227,6 +227,8 @@ protected:
         (
             TrackData& td,
             const scalar dt,           // timestep
+            const scalar Ts,           // Surface temperature
+            const scalar d,            // diameter
             const scalar T,            // temperature
             const scalar mass,         // mass
             const scalar mass0,        // mass (initial on injection)
@@ -236,7 +238,10 @@ protected:
             bool& canCombust,          // 'can combust' flag
             scalarField& dMassDV,      // mass transfer - local to particle
             scalar& Sh,                // explicit particle enthalpy source
-            scalar& dhsTrans           // sensible enthalpy transfer to carrier
+            scalar& dhsTrans,          // sensible enthalpy transfer to carrier
+            scalar& N,                 // flux of species emitted from particle
+            scalar& NCpW,              // sum of N*Cp*W of emission species
+            scalarField& Cs            // carrier conc. of emission species
         ) const;
 
         //- Calculate surface reactions
@@ -250,7 +255,7 @@ protected:
             const scalar T,            // temperature
             const scalar mass,         // mass
             const bool canCombust,     // 'can combust' flag
-            const scalarField& dMassVolatile, // mass transfer of volatiles
+            const scalar N,            // flux of species emitted from particle
             const scalarField& YMix,   // mixture mass fractions
             const scalarField& YGas,   // gas-phase mass fractions
             const scalarField& YLiquid,// liquid-phase mass fractions

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

@@ -26,6 +26,9 @@ License
 
 #include "ReactingParcel.H"
 #include "mathConstants.H"
+#include "specie.H"
+
+using namespace Foam::constant;
 
 // * * * * * * * * * * *  Protected Member Functions * * * * * * * * * * * * //
 
@@ -98,6 +101,92 @@ void Foam::ReactingParcel<ParcelType>::cellValueSourceCorrection
 }
 
 
+template<class ParcelType>
+template<class TrackData>
+void Foam::ReactingParcel<ParcelType>::correctSurfaceValues
+(
+    TrackData& td,
+    const label cellI,
+    const scalar T,
+    const scalarField& Cs,
+    scalar& rhos,
+    scalar& mus,
+    scalar& Pr,
+    scalar& kappa
+)
+{
+    // No correction if total concentration of emitted species is small
+    if (sum(Cs) < SMALL)
+    {
+        return;
+    }
+
+    // Far field gas molar fractions
+    scalarField Xinf(Y_.size());
+
+    forAll(Xinf, i)
+    {
+        Xinf[i] =
+            td.cloud().mcCarrierThermo().Y(i)[cellI]
+           /td.cloud().mcCarrierThermo().speciesData()[i].W();
+    }
+    Xinf /= sum(Xinf);
+
+    // Molar fraction of far field species at particle surface
+    const scalar Xsff = 1.0 - min(sum(Cs)*specie::RR*this->T_/pc_, 1.0);
+
+    // Surface gas total molar concentration
+    const scalar CsTot = pc_/(specie::RR*this->T_);
+
+    // Surface carrier composition (molar fraction)
+    scalarField Xs(Xinf.size());
+
+    // Surface carrier composition (mass fraction)
+    scalarField Ys(Xinf.size());
+
+    forAll(Xs, i)
+    {
+        // Molar concentration of species at particle surface
+        const scalar Csi = Cs[i] + Xsff*Xinf[i]*CsTot;
+
+        Xs[i] = (2.0*Csi + Xinf[i]*CsTot)/3.0;
+        Ys[i] = Xs[i]*td.cloud().mcCarrierThermo().speciesData()[i].W();
+    }
+    Xs /= sum(Xs);
+    Ys /= sum(Ys);
+
+
+    rhos = 0;
+    mus = 0;
+    kappa = 0;
+    scalar cps = 0.0;
+    scalar sumYiSqrtW = 0;
+    scalar sumYiCbrtW = 0;
+
+    forAll(Ys, i)
+    {
+        const scalar sqrtW =
+            sqrt(td.cloud().mcCarrierThermo().speciesData()[i].W());
+        const scalar cbrtW =
+            cbrt(td.cloud().mcCarrierThermo().speciesData()[i].W());
+
+        rhos += Xs[i]*td.cloud().mcCarrierThermo().speciesData()[i].W();
+        cps += Xs[i]*td.cloud().mcCarrierThermo().speciesData()[i].Cp(T);
+        mus += Ys[i]*sqrtW*td.cloud().mcCarrierThermo().speciesData()[i].mu(T);
+        kappa +=
+            Ys[i]*cbrtW*td.cloud().mcCarrierThermo().speciesData()[i].kappa(T);
+
+        sumYiSqrtW += Ys[i]*sqrtW;
+        sumYiCbrtW += Ys[i]*cbrtW;
+    }
+
+    rhos *= pc_/(specie::RR*T);
+    mus /= sumYiSqrtW;
+    kappa /= sumYiCbrtW;
+    Pr = cps*mus/kappa;
+}
+
+
 template<class ParcelType>
 Foam::scalar Foam::ReactingParcel<ParcelType>::updateMassFraction
 (
@@ -140,6 +229,20 @@ void Foam::ReactingParcel<ParcelType>::calc
     const scalar cp0 = this->cp_;
     const scalar mass0 = this->mass();
 
+
+    // Calc surface values
+    // ~~~~~~~~~~~~~~~~~~~
+    scalar Ts, rhos, mus, Pr, kappa;
+    ThermoParcel<ParcelType>::
+        calcSurfaceValues(td, cellI, T0, Ts, rhos, mus, Pr, kappa);
+
+    // Reynolds number
+    scalar Re = this->Re(U0, d0, mus);
+
+
+    // Sources
+    //~~~~~~~~
+
     // Explicit momentum source for particle
     vector Su = vector::zero;
 
@@ -159,24 +262,41 @@ void Foam::ReactingParcel<ParcelType>::calc
     // Mass transfer due to phase change
     scalarField dMassPC(Y_.size(), 0.0);
 
+    // Molar flux of species emitted from the particle (kmol/m^2/s)
+    scalar Ne = 0.0;
+
+    // Sum Ni*Cpi*Wi of emission species
+    scalar NCpW = 0.0;
+
+    // Surface concentrations of emitted species
+    scalarField Cs(td.cloud().mcCarrierThermo().species().size(), 0.0);
+
     // Calc mass and enthalpy transfer due to phase change
     calcPhaseChange
     (
         td,
         dt,
         cellI,
+        Ts,
+        mus/rhos,
+        Re,
         d0,
         T0,
-        U0,
         mass0,
         0,
         1.0,
         Y_,
         dMassPC,
         Sh,
-        dhsTrans
+        dhsTrans,
+        Ne,
+        NCpW,
+        Cs
     );
 
+    // Correct surface values due to emitted species
+    correctSurfaceValues(td, cellI, Ts, Cs, rhos, mus, Pr, kappa);
+
     // Update particle component mass and mass fractions
     scalar mass1 = updateMassFraction(mass0, dMassPC, Y_);
 
@@ -186,14 +306,30 @@ void Foam::ReactingParcel<ParcelType>::calc
 
     // Calculate new particle temperature
     scalar T1 =
-        calcHeatTransfer(td, dt, cellI, d0, U0, rho0, T0, cp0, Sh, dhsTrans);
+        calcHeatTransfer
+        (
+            td,
+            dt,
+            cellI,
+            Re,
+            Pr,
+            kappa,
+            d0,
+            rho0,
+            T0,
+            cp0,
+            NCpW,
+            Sh,
+            dhsTrans
+        );
 
 
     // Motion
     // ~~~~~~
 
     // Calculate new particle velocity
-    vector U1 = calcVelocity(td, dt, cellI, d0, U0, rho0, mass0, Su, dUTrans);
+    vector U1 =
+        calcVelocity(td, dt, cellI, Re, mus, d0, U0, rho0, mass0, Su, dUTrans);
 
     dUTrans += 0.5*(mass0 - mass1)*(U0 + U1);
 
@@ -258,7 +394,7 @@ void Foam::ReactingParcel<ParcelType>::calc
         }
         else
         {
-            this->d_ = cbrt(mass1/this->rho_*6.0/constant::math::pi);
+            this->d_ = cbrt(mass1/this->rho_*6.0/math::pi);
         }
     }
 }
@@ -271,16 +407,21 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
     TrackData& td,
     const scalar dt,
     const label cellI,
+    const scalar Re,
+    const scalar Ts,
+    const scalar nus,
     const scalar d,
     const scalar T,
-    const vector& U,
     const scalar mass,
     const label idPhase,
     const scalar YPhase,
     const scalarField& YComponents,
     scalarField& dMassPC,
     scalar& Sh,
-    scalar& dhsTrans
+    scalar& dhsTrans,
+    scalar& N,
+    scalar& NCpW,
+    scalarField& Cs
 )
 {
     if
@@ -298,12 +439,12 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
     (
         dt,
         cellI,
+        Re,
         d,
-        min(T, td.constProps().Tbp()), // Limit to boiling temperature
+        nus,
+        T,
+        Ts,
         pc_,
-        this->Tc_,
-        this->muc_/(this->rhoc_ + ROOTVSMALL),
-        U - this->Uc_,
         dMassPC
     );
 
@@ -315,18 +456,38 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
     // Add to cumulative phase change mass
     td.cloud().addToMassPhaseChange(this->nParticle_*dMassTot);
 
-    // Enthalphy transfer to carrier phase
+    // Average molecular weight of carrier mix - assumes perfect gas
-    label id;
+    scalar Wc = this->rhoc_*specie::RR*this->Tc_/this->pc_;
+
     forAll(YComponents, i)
     {
-        id = td.cloud().composition().localToGlobalCarrierId(idPhase, i);
+        const label idc =
-        const scalar hv = td.cloud().mcCarrierThermo().speciesData()[id].H(T);
+            td.cloud().composition().localToGlobalCarrierId(idPhase, i);
+        const scalar hv = td.cloud().mcCarrierThermo().speciesData()[idc].H(T);
 
-        id = td.cloud().composition().globalIds(idPhase)[i];
+        const label idl = td.cloud().composition().globalIds(idPhase)[i];
         const scalar hl =
-            td.cloud().composition().liquids().properties()[id].h(pc_, T);
+            td.cloud().composition().liquids().properties()[idl].h(pc_, T);
 
+        // Enthalphy transfer to carrier phase
         Sh += dMassPC[i]*(hl - hv)/dt;
+
+        const scalar Dab =
+            td.cloud().composition().liquids().properties()[idl].D(pc_, Ts, Wc);
+
+        const scalar Cp =
+            td.cloud().mcCarrierThermo().speciesData()[idc].Cp(Ts);
+        const scalar W = td.cloud().mcCarrierThermo().speciesData()[idc].W();
+        const scalar Ni = dMassPC[i]/(this->areaS(d)*dt*W);
+
+        // 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/ReactingParcel/ReactingParcel.H

@@ -195,16 +195,21 @@ protected:
             TrackData& td,
             const scalar dt,           // timestep
             const label cellI,         // owner cell
+            const scalar Re,           // Reynolds number
+            const scalar Ts,           // Surface temperature
+            const scalar nus,          // Surface kinematic viscosity
             const scalar d,            // diameter
             const scalar T,            // temperature
-            const vector& U,           // velocity
             const scalar mass,         // mass
             const label idPhase,       // id of phase involved in phase change
             const scalar YPhase,       // total mass fraction
             const scalarField& YComponents, // component mass fractions
             scalarField& dMassPC,      // mass transfer - local to particle
             scalar& Sh,                // explicit particle enthalpy source
-            scalar& dhsTrans           // sensible enthalpy transfer to carrier
+            scalar& dhsTrans,          // sensible enthalpy transfer to carrier
+            scalar& N,                 // flux of species emitted from particle
+            scalar& NCpW,              // sum of N*Cp*W of emission species
+            scalarField& Cs            // carrier conc. of emission species
         );
 
         //- Update mass fraction
@@ -316,6 +321,20 @@ public:
                 const label cellI
             );
 
+            //- Correct surface values due to emitted species
+            template<class TrackData>
+            void correctSurfaceValues
+            (
+                TrackData& td,
+                const label cellI,
+                const scalar T,
+                const scalarField& Cs,
+                scalar& rhos,
+                scalar& mus,
+                scalar& Pr,
+                scalar& kappa
+            );
+
             //- Update parcel properties over the time interval
             template<class TrackData>
             void calc

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

@@ -27,6 +27,8 @@ License
 #include "ThermoParcel.H"
 #include "physicoChemicalConstants.H"
 
+using namespace Foam::constant;
+
 // * * * * * * * * * * *  Protected Member Functions * * * * * * * * * * * * //
 
 template<class ParcelType>
@@ -78,6 +80,33 @@ void Foam::ThermoParcel<ParcelType>::cellValueSourceCorrection
 }
 
 
+template<class ParcelType>
+template<class TrackData>
+void Foam::ThermoParcel<ParcelType>::calcSurfaceValues
+(
+    TrackData& td,
+    const label cellI,
+    const scalar T,
+    scalar& Ts,
+    scalar& rhos,
+    scalar& mus,
+    scalar& Pr,
+    scalar& kappa
+) const
+{
+    // Surface temperature using two thirds rule
+    Ts = (2.0*T + Tc_)/3.0;
+
+    // Assuming thermo props vary linearly with T for small dT
+    scalar factor = td.TInterp().interpolate(this->position(), cellI)/Ts;
+    rhos = this->rhoc_*factor;
+    mus = td.muInterp().interpolate(this->position(), cellI)/factor;
+
+    Pr = td.constProps().Pr();
+    kappa = cpc_*mus/Pr;
+}
+
+
 template<class ParcelType>
 template<class TrackData>
 void Foam::ThermoParcel<ParcelType>::calc
@@ -97,6 +126,19 @@ void Foam::ThermoParcel<ParcelType>::calc
     const scalar cp0 = this->cp_;
     const scalar mass0 = this->mass();
 
+
+    // Calc surface values
+    // ~~~~~~~~~~~~~~~~~~~
+    scalar Ts, rhos, mus, Pr, kappa;
+    calcSurfaceValues(td, cellI, T0, Ts, rhos, mus, Pr, kappa);
+
+    // Reynolds number
+    scalar Re = this->Re(U0, d0, mus);
+
+
+    // Sources
+    // ~~~~~~~
+
     // Explicit momentum source for particle
     vector Su = vector::zero;
 
@@ -109,19 +151,39 @@ void Foam::ThermoParcel<ParcelType>::calc
     // Sensible enthalpy transfer from the particle to the carrier phase
     scalar dhsTrans = 0.0;
 
+
     // Heat transfer
     // ~~~~~~~~~~~~~
 
+    // Sum Ni*Cpi*Wi of emission species
+    scalar NCpW = 0.0;
+
     // Calculate new particle velocity
     scalar T1 =
-        calcHeatTransfer(td, dt, cellI, d0, U0, rho0, T0, cp0, Sh, dhsTrans);
+        calcHeatTransfer
+        (
+            td,
+            dt,
+            cellI,
+            Re,
+            Pr,
+            kappa,
+            d0,
+            rho0,
+            T0,
+            cp0,
+            NCpW,
+            Sh,
+            dhsTrans
+        );
 
 
     // Motion
     // ~~~~~~
 
     // Calculate new particle velocity
-    vector U1 = calcVelocity(td, dt, cellI, d0, U0, rho0, mass0, Su, dUTrans);
+    vector U1 =
+        calcVelocity(td, dt, cellI, Re, mus, d0, U0, rho0, mass0, Su, dUTrans);
 
 
     //  Accumulate carrier phase source terms
@@ -149,11 +211,14 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
     TrackData& td,
     const scalar dt,
     const label cellI,
+    const scalar Re,
+    const scalar Pr,
+    const scalar kappa,
     const scalar d,
-    const vector& U,
     const scalar rho,
     const scalar T,
     const scalar cp,
+    const scalar NCpW,
     const scalar Sh,
     scalar& dhsTrans
 )
@@ -164,16 +229,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
     }
 
     // Calc heat transfer coefficient
-    scalar htc = td.cloud().heatTransfer().h
+    scalar htc = td.cloud().heatTransfer().htc(d, Re, Pr, kappa, NCpW);
-    (
-        d,
-        U - this->Uc_,
-        this->rhoc_,
-        rho,
-        cpc_,
-        cp,
-        this->muc_
-    );
 
     const scalar As = this->areaS(d);
 
@@ -189,7 +245,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
     {
         const scalarField& G =
             td.cloud().mesh().objectRegistry::lookupObject<volScalarField>("G");
-        const scalar sigma = constant::physicoChemical::sigma.value();
+        const scalar sigma = physicoChemical::sigma.value();
         const scalar epsilon = td.constProps().epsilon0();
 
         ap =

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

@@ -97,6 +97,9 @@ public:
             //- Particle scattering factor [] (radiation)
             const scalar f0_;
 
+            //- Default carrier Prandtl number []
+            const scalar Pr_;
+
 
     public:
 
@@ -124,6 +127,9 @@ public:
                 //- Return const access to the particle scattering factor []
                 //  Active for radiation only
                 inline scalar f0() const;
+
+                //- Return const access to the default carrier Prandtl number []
+                inline scalar Pr() const;
     };
 
 
@@ -217,11 +223,14 @@ protected:
             TrackData& td,
             const scalar dt,           // timestep
             const label cellI,         // owner cell
+            const scalar Re,           // Reynolds number
+            const scalar Pr,           // Prandtl number - surface
+            const scalar kappa,        // Thermal conductivity - surface
             const scalar d,            // diameter
-            const vector& U,           // velocity
             const scalar rho,          // density
             const scalar T,            // temperature
             const scalar cp,           // specific heat capacity
+            const scalar NCpW,         // Sum of N*Cp*W of emission species
             const scalar Sh,           // explicit particle enthalpy source
             scalar& dhsTrans           // sensible enthalpy transfer to carrier
         );
@@ -323,6 +332,20 @@ public:
                 const label cellI
             );
 
+            //- Calculate surface thermo properties
+            template<class TrackData>
+            void calcSurfaceValues
+            (
+                TrackData& td,
+                const label cellI,
+                const scalar T,
+                scalar& Ts,
+                scalar& rhos,
+                scalar& mus,
+                scalar& Pr,
+                scalar& kappa
+            ) const;
+
             //- Update parcel properties over the time interval
             template<class TrackData>
             void calc

src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelI.H

@@ -37,7 +37,8 @@ inline Foam::ThermoParcel<ParcelType>::constantProperties::constantProperties
     TMin_(dimensionedScalar(this->dict().lookup("TMin")).value()),
     cp0_(dimensionedScalar(this->dict().lookup("cp0")).value()),
     epsilon0_(dimensionedScalar(this->dict().lookup("epsilon0")).value()),
-    f0_(dimensionedScalar(this->dict().lookup("f0")).value())
+    f0_(dimensionedScalar(this->dict().lookup("f0")).value()),
+    Pr_(dimensionedScalar(this->dict().lookup("Pr")).value())
 {}
 
 
@@ -159,6 +160,14 @@ Foam::ThermoParcel<ParcelType>::constantProperties::f0() const
 }
 
 
+template <class ParcelType>
+inline Foam::scalar
+Foam::ThermoParcel<ParcelType>::constantProperties::Pr() const
+{
+    return Pr_;
+}
+
+
 // * * * * * * * * * * * trackData Member Functions  * * * * * * * * * * * * //
 
 template<class ParcelType>

src/lagrangian/intermediate/parcels/include/makeParcelHeatTransferModels.H

@@ -33,6 +33,7 @@ License
 
 #include "NoHeatTransfer.H"
 #include "RanzMarshall.H"
+#include "RanzMarshallBirdCorrection.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -52,6 +53,12 @@ License
         ThermoCloud,                                                          \
         ParcelType                                                            \
     );                                                                        \
+    makeHeatTransferModelType                                                 \
+    (                                                                         \
+        RanzMarshallBirdCorrection,                                           \
+        ThermoCloud,                                                          \
+        ParcelType                                                            \
+    );
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/lagrangian/intermediate/parcels/include/makeReactingParcelHeatTransferModels.H

@@ -34,6 +34,7 @@ License
 
 #include "NoHeatTransfer.H"
 #include "RanzMarshall.H"
+#include "RanzMarshallBirdCorrection.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -74,6 +75,13 @@ License
         ParcelType,                                                           \
         ThermoType                                                            \
     );                                                                        \
+    makeHeatTransferModelThermoType                                           \
+    (                                                                         \
+        RanzMarshallBirdCorrection,                                           \
+        ThermoCloud,                                                          \
+        ParcelType,                                                           \
+        ThermoType                                                            \
+    );
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/lagrangian/intermediate/submodels/Kinematic/DragModel/DragModel/DragModel.C

@@ -66,19 +66,14 @@ const Foam::dictionary& Foam::DragModel<CloudType>::dict() const
 template<class CloudType>
 Foam::scalar Foam::DragModel<CloudType>::utc
 (
-    const vector& Ur,
+    const scalar Re,
     const scalar d,
-    const scalar rhoc,
     const scalar mu
 ) const
 {
-    const scalar magUr = mag(Ur);
-
-    const scalar Re = rhoc*magUr*d/(mu + ROOTVSMALL);
-
     const scalar Cd = this->Cd(Re);
 
-    return Cd*rhoc*magUr/8.0;
+    return Cd*Re/d*mu/8.0;
 }
 
 

src/lagrangian/intermediate/submodels/Kinematic/DragModel/DragModel/DragModel.H

@@ -122,13 +122,7 @@ public:
 
         //- Return momentum transfer coefficient
         //  Drag force per unit particle surface area = utc(U - Up)
-        scalar utc
+        scalar utc(const scalar Re, const scalar d, const scalar mu) const;
-        (
-            const vector& Ur,
-            const scalar d,
-            const scalar rhoc,
-            const scalar mu
-        ) const;
 };
 
 

src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C

@@ -38,15 +38,14 @@ Foam::scalarField Foam::LiquidEvaporation<CloudType>::calcXc
 {
     scalarField Xc(this->owner().mcCarrierThermo().Y().size());
 
-    scalar Winv = 0.0;
     forAll(Xc, i)
     {
-        scalar Y = this->owner().mcCarrierThermo().Y()[i][cellI];
+        Xc[i] =
-        Winv += Y/this->owner().mcCarrierThermo().speciesData()[i].W();
+            this->owner().mcCarrierThermo().Y()[i][cellI]
-        Xc[i] = Y/this->owner().mcCarrierThermo().speciesData()[i].W();
+           /this->owner().mcCarrierThermo().speciesData()[i].W();
     }
 
-    return Xc/Winv;
+    return Xc/sum(Xc);
 }
 
 
@@ -136,12 +135,12 @@ void Foam::LiquidEvaporation<CloudType>::calculate
 (
     const scalar dt,
     const label cellI,
+    const scalar Re,
     const scalar d,
+    const scalar nu,
     const scalar T,
+    const scalar Ts,
     const scalar pc,
-    const scalar Tc,
-    const scalar nuc,
-    const vector& Ur,
     scalarField& dMassPC
 ) const
 {
@@ -151,29 +150,25 @@ void Foam::LiquidEvaporation<CloudType>::calculate
     // droplet surface area
     scalar A = constant::math::pi*sqr(d);
 
-    // Reynolds number
-    scalar Re = mag(Ur)*d/(nuc + ROOTVSMALL);
-
-    // film temperature evaluated using the 2/3 rule
-    scalar Tf = (2.0*T + Tc)/3.0;
-
     // calculate mass transfer of each specie in liquid
     forAll(activeLiquids_, i)
     {
         label gid = liqToCarrierMap_[i];
         label lid = liqToLiqMap_[i];
 
-        // vapour diffusivity at film temperature and cell pressure [m2/s]
+        // vapour diffusivity [m2/s]
-        scalar Dab = liquids_->properties()[lid].D(pc, Tf);
+        scalar Dab = liquids_->properties()[lid].D(pc, Ts);
 
-        // saturation pressure for species i at film temperature and cell
+        // Saturation pressure for species i [pa]
-        // pressure [pa] - carrier phase pressure assumed equal to the liquid
+        // - carrier phase pressure assumed equal to the liquid vapour pressure
-        // vapour pressure close to the surface
+        //   close to the surface
-        // - limited to pc if pSat > pc
+        // NOTE: if pSat > pc then particle is superheated
-        scalar pSat = min(liquids_->properties()[lid].pv(pc, Tf), pc);
+        // calculated evaporation rate will be greater than that of a particle
+        // at boiling point, but this is not a boiling model
+        scalar pSat = liquids_->properties()[lid].pv(pc, T);
 
         // Schmidt number
-        scalar Sc = nuc/(Dab + ROOTVSMALL);
+        scalar Sc = nu/(Dab + ROOTVSMALL);
 
         // Sherwood number
         scalar Sh = this->Sh(Re, Sc);
@@ -181,11 +176,11 @@ void Foam::LiquidEvaporation<CloudType>::calculate
         // mass transfer coefficient [m/s]
         scalar kc = Sh*Dab/(d + ROOTVSMALL);
 
-        // vapour concentration at droplet surface at film temperature [kmol/m3]
+        // vapour concentration at droplet surface [kmol/m3] at film temperature
-        scalar Cs = pSat/(specie::RR*Tf);
+        scalar Cs = pSat/(specie::RR*Ts);
 
-        // vapour concentration in bulk gas at film temperature [kmol/m3]
+        // vapour concentration in bulk gas [kmol/m3] at film temperature
-        scalar Cinf = Xc[gid]*pc/(specie::RR*Tf);
+        scalar Cinf = Xc[gid]*pc/(specie::RR*Ts);
 
         // molar flux of vapour [kmol/m2/s]
         scalar Ni = max(kc*(Cs - Cinf), 0.0);

src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.H

@@ -106,12 +106,12 @@ public:
         (
             const scalar dt,
             const label cellI,
+            const scalar Re,
             const scalar d,
+            const scalar nu,
             const scalar T,
+            const scalar Ts,
             const scalar pc,
-            const scalar Tc,
-            const scalar nuc,
-            const vector& Ur,
             scalarField& dMassPC
         ) const;
 };

src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/NoPhaseChange/NoPhaseChange.C

@@ -65,7 +65,7 @@ void Foam::NoPhaseChange<CloudType>::calculate
     const scalar,
     const scalar,
     const scalar,
-    const vector&,
+    const scalar,
     scalarField&
 ) const
 {

src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/NoPhaseChange/NoPhaseChange.H

@@ -74,12 +74,12 @@ public:
         (
             const scalar dt,
             const label cellI,
+            const scalar Re,
             const scalar d,
+            const scalar nu,
             const scalar T,
+            const scalar Ts,
             const scalar pc,
-            const scalar Tc,
-            const scalar nuc,
-            const vector& Ur,
             scalarField& dMassPC
         ) const;
 };

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

@@ -140,12 +140,12 @@ public:
         (
             const scalar dt,
             const label cellI,
+            const scalar Re,
             const scalar d,
+            const scalar nu,
             const scalar T,
+            const scalar Ts,
             const scalar pc,
-            const scalar Tc,
-            const scalar nuc,
-            const vector& Ur,
             scalarField& dMassPC
         ) const = 0;
 };

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

@@ -70,7 +70,7 @@ Foam::scalar Foam::NoSurfaceReaction<CloudType>::calculate
     const scalarField&,
     const scalarField&,
     const scalarField&,
-    const scalarField&,
+    const scalar,
     scalarField&,
     scalarField&,
     scalarField&,

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

@@ -88,7 +88,7 @@ public:
             const scalarField& YLiquid,
             const scalarField& YSolid,
             const scalarField& YMixture,
-            const scalarField& dMassVolatile,
+            const scalar N,
             scalarField& dMassGas,
             scalarField& dMassLiquid,
             scalarField& dMassSolid,

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

@@ -146,7 +146,7 @@ public:
             const scalarField& YLiquid,
             const scalarField& YSolid,
             const scalarField& YMixture,
-            const scalarField& dMassVolatile,
+            const scalar N,
             scalarField& dMassGas,
             scalarField& dMassLiquid,
             scalarField& dMassSolid,

src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/HeatTransferModel/HeatTransferModel.C

@@ -82,57 +82,21 @@ const Foam::dictionary& Foam::HeatTransferModel<CloudType>::coeffDict() const
 
 
 template<class CloudType>
-Foam::scalar Foam::HeatTransferModel<CloudType>::h
+Foam::scalar Foam::HeatTransferModel<CloudType>::htc
 (
     const scalar dp,
-    const vector& Ur,
+    const scalar Re,
-    const scalar rhoc,
+    const scalar Pr,
-    const scalar rhop,
+    const scalar kappa,
-    const scalar cpc,
+    const scalar NCpW
-    const scalar cpp,
-    const scalar muc
 ) const
 {
-    const scalar Re = rhoc*mag(Ur)*dp/(muc + ROOTVSMALL);
-
-//    const scalar Pr = muc/alphac;
-    const scalar Pr = this->Pr();
-
     const scalar Nu = this->Nu(Re, Pr);
 
-    const scalar kappa = cpc*muc/Pr;
-
     return Nu*kappa/dp;
 }
 
 
-template<class CloudType>
-Foam::scalar Foam::HeatTransferModel<CloudType>::Cu
-(
-    const scalar dp,
-    const vector& Ur,
-    const scalar rhoc,
-    const scalar rhop,
-    const scalar cpc,
-    const scalar cpp,
-    const scalar muc
-) const
-{
-    const scalar Re = rhoc*mag(Ur)*dp/(muc + ROOTVSMALL);
-
-//    const scalar Pr = muc/alphac;
-    const scalar Pr = this->Pr();
-
-    const scalar Nu = this->Nu(Re, Pr);
-
-    const scalar kappa = cpc*muc/Pr;
-
-    const scalar htc = Nu*kappa/dp;
-
-    return 6.0*htc/(dp*rhop*cpp);
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #include "NewHeatTransferModel.C"

src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/HeatTransferModel/HeatTransferModel.H

@@ -138,27 +138,13 @@ public:
         virtual scalar Pr() const = 0;
 
         //- Return heat transfer coefficient
-        virtual scalar h
+        virtual scalar htc
         (
             const scalar dp,
-            const vector& Ur,
+            const scalar Re,
-            const scalar rhoc,
+            const scalar Pr,
-            const scalar rhop,
+            const scalar kappa,
-            const scalar cpc,
+            const scalar NCpW
-            const scalar cpp,
-            const scalar muc
-        ) const;
-
-        //- Return linearised coefficient for temperature equation
-        virtual scalar Cu
-        (
-            const scalar dp,
-            const vector& Ur,
-            const scalar rhoc,
-            const scalar rhop,
-            const scalar cpc,
-            const scalar cpp,
-            const scalar muc
         ) const;
 };
 

src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/RanzMarshallBirdCorrection/RanzMarshallBirdCorrection.C

@@ -0,0 +1,79 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "RanzMarshallBirdCorrection.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::RanzMarshallBirdCorrection<CloudType>::RanzMarshallBirdCorrection
+(
+    const dictionary& dict,
+    CloudType& cloud
+)
+:
+    RanzMarshall<CloudType>(dict, cloud)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::RanzMarshallBirdCorrection<CloudType>::~RanzMarshallBirdCorrection()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::scalar Foam::RanzMarshallBirdCorrection<CloudType>::htc
+(
+    const scalar dp,
+    const scalar Re,
+    const scalar Pr,
+    const scalar kappa,
+    const scalar NCpW
+) const
+{
+    scalar htc = RanzMarshall<CloudType>::htc(dp, Re, Pr, kappa, NCpW);
+
+    // Bird correction
+    if (mag(htc) > ROOTVSMALL && mag(NCpW) > ROOTVSMALL)
+    {
+        const scalar phit = min(NCpW/htc, 50);
+        scalar fBird = 1.0;
+        if (phit > 0.001)
+        {
+            fBird = phit/(exp(phit) - 1.0);
+        }
+        htc *= fBird;
+    }
+
+    return htc;
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/RanzMarshallBirdCorrection/RanzMarshallBirdCorrection.H

@@ -0,0 +1,101 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::RanzMarshallBirdCorrection
+
+Description
+    The Ranz-Marshall correlation for heat transfer with the Bird correction
+    to account for the local shielding effect due to emitted species.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef RanzMarshallBirdCorrection_H
+#define RanzMarshallBirdCorrection_H
+
+#include "RanzMarshall.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                 Class RanzMarshallBirdCorrection Declaration
+\*---------------------------------------------------------------------------*/
+
+template <class CloudType>
+class RanzMarshallBirdCorrection
+:
+    public RanzMarshall<CloudType>
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("RanzMarshallBirdCorrection");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        RanzMarshallBirdCorrection
+        (
+            const dictionary& dict,
+            CloudType& cloud
+        );
+
+
+    //- Destructor
+    virtual ~RanzMarshallBirdCorrection();
+
+
+    // Member Functions
+
+        //- Return heat transfer coefficient
+        virtual scalar htc
+        (
+            const scalar dp,
+            const scalar Re,
+            const scalar Pr,
+            const scalar kappa,
+            const scalar NCpW
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "RanzMarshallBirdCorrection.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //