ENH: Updated lagrangian/intermediate library to operate using sensible enthalpy

- removed support for total enthalpy
- much cleaner/improved code for transfer terms

src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C

@@ -160,20 +160,6 @@ Foam::ThermoCloud<ParcelType>::ThermoCloud
         ),
         this->mesh(),
         dimensionedScalar("zero", dimEnergy, 0.0)
-    ),
-    hcTrans_
-    (
-        IOobject
-        (
-            this->name() + "hcTrans",
-            this->db().time().timeName(),
-            this->db(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        this->mesh(),
-        dimensionedScalar("zero", dimEnergy, 0.0)
     )
 {
     if (readFields)
@@ -220,7 +206,6 @@ void Foam::ThermoCloud<ParcelType>::resetSourceTerms()
 {
     KinematicCloud<ParcelType>::resetSourceTerms();
     hsTrans_.field() = 0.0;
-    hcTrans_.field() = 0.0;
 }
 
 

src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H

@@ -108,15 +108,9 @@ protected:
 
         // Sources
 
-            //- Sensible enthalpy transfer
+            //- Sensible enthalpy transfer [J/kg]
             DimensionedField<scalar, volMesh> hsTrans_;
 
-            //- Chemical enthalpy transfer
-            // - If solving for total enthalpy, the carrier phase enthalpy will
-            //   receive the full enthalpy of reaction via creation of reaction
-            //   products
-            DimensionedField<scalar, volMesh> hcTrans_;
-
 
     // Protected member functions
 
@@ -194,19 +188,10 @@ public:
 
                 // Enthalpy
 
-                    //- Return reference to sensible enthalpy source
+                    //- Sensible enthalpy transfer [J/kg]
                     inline DimensionedField<scalar, volMesh>& hsTrans();
 
-                    //- Return tmp total sensible enthalpy source term
+                    //- Return enthalpy source [J/kg/m3/s]
-                    inline tmp<DimensionedField<scalar, volMesh> > Shs() const;
-
-                    //- Return reference to chemical enthalpy source
-                    inline DimensionedField<scalar, volMesh>& hcTrans();
-
-                    //- Return tmp chemical enthalpy source term
-                    inline tmp<DimensionedField<scalar, volMesh> > Shc() const;
-
-                    //- Return tmp total enthalpy source term
                     inline tmp<DimensionedField<scalar, volMesh> > Sh() const;
 
 

src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H

@@ -85,82 +85,6 @@ Foam::ThermoCloud<ParcelType>::hsTrans()
 }
 
 
-template<class ParcelType>
-inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
-Foam::ThermoCloud<ParcelType>::Shs() const
-{
-    tmp<DimensionedField<scalar, volMesh> > tShs
-    (
-        new DimensionedField<scalar, volMesh>
-        (
-            IOobject
-            (
-                this->name() + "Shs",
-                this->db().time().timeName(),
-                this->mesh(),
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE,
-                false
-            ),
-            this->mesh(),
-            dimensionedScalar
-            (
-                "zero",
-                dimMass/dimLength/pow3(dimTime),
-                0.0
-            )
-        )
-    );
-
-    scalarField& Shs = tShs().field();
-    Shs = hsTrans_/(this->mesh().V()*this->db().time().deltaT());
-
-    return tShs;
-}
-
-
-template<class ParcelType>
-inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
-Foam::ThermoCloud<ParcelType>::hcTrans()
-{
-    return hcTrans_;
-}
-
-
-template<class ParcelType>
-inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
-Foam::ThermoCloud<ParcelType>::Shc() const
-{
-    tmp<DimensionedField<scalar, volMesh> > tShc
-    (
-        new DimensionedField<scalar, volMesh>
-        (
-            IOobject
-            (
-                this->name() + "Shc",
-                this->db().time().timeName(),
-                this->mesh(),
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE,
-                false
-            ),
-            this->mesh(),
-            dimensionedScalar
-            (
-                "zero",
-                dimMass/dimLength/pow3(dimTime),
-                0.0
-            )
-        )
-    );
-
-    scalarField& Shc = tShc().field();
-    Shc = hcTrans_/(this->mesh().V()*this->db().time().deltaT());
-
-    return tShc;
-}
-
-
 template<class ParcelType>
 inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
 Foam::ThermoCloud<ParcelType>::Sh() const
@@ -178,19 +102,10 @@ Foam::ThermoCloud<ParcelType>::Sh() const
                 IOobject::AUTO_WRITE,
                 false
             ),
-            this->mesh(),
+            hsTrans_/(this->mesh().V()*this->db().time().deltaT())
-            dimensionedScalar
-            (
-                "zero",
-                dimMass/dimLength/pow3(dimTime),
-                0.0
-            )
         )
     );
 
-    scalarField& Sh = tSh().field();
-    Sh = (hsTrans_ + hcTrans_)/(this->mesh().V()*this->db().time().deltaT());
-
     return tSh;
 }
 
@@ -205,7 +120,7 @@ Foam::ThermoCloud<ParcelType>::Ep() const
         (
             IOobject
             (
-                this->name() + "radiationEp",
+                this->name() + "radiation::Ep",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
@@ -248,7 +163,7 @@ Foam::ThermoCloud<ParcelType>::ap() const
         (
             IOobject
             (
-                this->name() + "radiationAp",
+                this->name() + "radiation::ap",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
@@ -291,7 +206,7 @@ Foam::ThermoCloud<ParcelType>::sigmap() const
         (
             IOobject
             (
-                this->name() + "radiationSigmap",
+                this->name() + "radiation::sigmap",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,

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

@@ -267,7 +267,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         YLiquid_,
         dMassPC,
         Sh,
-        dhsTrans,
         Ne,
         NCpW,
         Cs
@@ -296,7 +295,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         canCombust_,
         dMassDV,
         Sh,
-        dhsTrans,
         Ne,
         NCpW,
         Cs
@@ -398,19 +396,11 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         {
             label gid = td.cloud().composition().localToGlobalCarrierId(GAS, i);
             td.cloud().rhoTrans(gid)[cellI] += np0*dMassGas[i];
-            td.cloud().hcTrans()[cellI] +=
-                np0
-               *dMassGas[i]
-               *td.cloud().mcCarrierThermo().speciesData()[gid].Hc();
         }
         forAll(YLiquid_, i)
         {
             label gid = td.cloud().composition().localToGlobalCarrierId(LIQ, i);
             td.cloud().rhoTrans(gid)[cellI] += np0*dMassLiquid[i];
-            td.cloud().hcTrans()[cellI] +=
-                np0
-               *dMassLiquid[i]
-               *td.cloud().mcCarrierThermo().speciesData()[gid].Hc();
         }
 /*
         // No mapping between solid components and carrier phase
@@ -418,19 +408,11 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         {
             label gid = td.cloud().composition().localToGlobalCarrierId(SLD, i);
             td.cloud().rhoTrans(gid)[cellI] += np0*dMassSolid[i];
-            td.cloud().hcTrans()[cellI] +=
-                np0
-               *dMassSolid[i]
-               *td.cloud().mcCarrierThermo().speciesData()[gid].Hc();
         }
 */
         forAll(dMassSRCarrier, i)
         {
             td.cloud().rhoTrans(i)[cellI] += np0*dMassSRCarrier[i];
-            td.cloud().hcTrans()[cellI] +=
-                np0
-               *dMassSRCarrier[i]
-               *td.cloud().mcCarrierThermo().speciesData()[i].Hc();
         }
 
         // Update momentum transfer
@@ -476,7 +458,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
 */
             td.cloud().UTrans()[cellI] += np0*mass1*U1;
             td.cloud().hsTrans()[cellI] +=
-                np0*mass1*HEff(td, pc, T1, idG, idL, idS);
+                np0*mass1*HEff(td, pc, T1, idG, idL, idS); // using total h
         }
     }
 
@@ -520,7 +502,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
     bool& canCombust,
     scalarField& dMassDV,
     scalar& Sh,
-    scalar& dhsTrans,
     scalar& N,
     scalar& NCpW,
     scalarField& Cs

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

@@ -238,7 +238,6 @@ 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& N,                 // flux of species emitted from particle
             scalar& NCpW,              // sum of N*Cp*W of emission species
             scalarField& Cs            // carrier conc. of emission species

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

@@ -287,7 +287,6 @@ void Foam::ReactingParcel<ParcelType>::calc
         Y_,
         dMassPC,
         Sh,
-        dhsTrans,
         Ne,
         NCpW,
         Cs
@@ -341,10 +340,6 @@ void Foam::ReactingParcel<ParcelType>::calc
         {
             label gid = td.cloud().composition().localToGlobalCarrierId(0, i);
             td.cloud().rhoTrans(gid)[cellI] += np0*dMassPC[i];
-            td.cloud().hcTrans()[cellI] +=
-                np0
-               *dMassPC[i]
-               *td.cloud().mcCarrierThermo().speciesData()[gid].Hc();
         }
 
         // Update momentum transfer
@@ -371,7 +366,7 @@ void Foam::ReactingParcel<ParcelType>::calc
                 td.cloud().rhoTrans(gid)[cellI] += np0*mass1*Y_[i];
             }
             td.cloud().UTrans()[cellI] += np0*mass1*U1;
-            td.cloud().hcTrans()[cellI] +=
+            td.cloud().hsTrans()[cellI] +=
                 np0*mass1*td.cloud().composition().H(0, Y_, pc_, T1);
         }
     }
@@ -417,7 +412,6 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
     const scalarField& YComponents,
     scalarField& dMassPC,
     scalar& Sh,
-    scalar& dhsTrans,               // TODO: not used
     scalar& N,
     scalar& NCpW,
     scalarField& Cs
@@ -469,6 +463,7 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
             td.cloud().composition().localToGlobalCarrierId(idPhase, i);
         const label idl = td.cloud().composition().globalIds(idPhase)[i];
 
+        // Calculate enthalpy transfer
         if
         (
             td.cloud().phaseChange().enthalpyTransfer()

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

@@ -209,7 +209,6 @@ protected:
             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& N,                 // flux of species emitted from particle
             scalar& NCpW,              // sum of N*Cp*W of emission species
             scalarField& Cs            // carrier conc. of emission species

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

@@ -410,6 +410,130 @@ Foam::scalar Foam::CompositionModel<CloudType>::H
 }
 
 
+template<class CloudType>
+Foam::scalar Foam::CompositionModel<CloudType>::Hs
+(
+    const label phaseI,
+    const scalarField& Y,
+    const scalar p,
+    const scalar T
+) const
+{
+    const phaseProperties& props = phaseProps_[phaseI];
+    scalar HsMixture = 0.0;
+    switch (props.phase())
+    {
+        case phaseProperties::GAS:
+        {
+            forAll(Y, i)
+            {
+                label gid = props.globalIds()[i];
+                HsMixture += Y[i]*mcCarrierThermo_.speciesData()[gid].Hs(T);
+            }
+            break;
+        }
+        case phaseProperties::LIQUID:
+        {
+            forAll(Y, i)
+            {
+                label gid = props.globalIds()[i];
+                HsMixture +=
+                    Y[i]
+                   *(
+                       this->liquids().properties()[gid].h(p, T)
+                     - this->liquids().properties()[gid].h(p, 298.25)
+                    );
+            }
+            break;
+        }
+        case phaseProperties::SOLID:
+        {
+            forAll(Y, i)
+            {
+                label gid = props.globalIds()[i];
+                HsMixture += Y[i]*this->solids().properties()[gid].cp()*T;
+            }
+            break;
+        }
+        default:
+        {
+            FatalErrorIn
+            (
+                "Foam::scalar Foam::CompositionModel<CloudType>::Hs"
+                "("
+                "    const label, "
+                "    const scalarField&, "
+                "    const scalar, "
+                "    const scalar"
+                ") const"
+            )   << "Unknown phase enumeration" << nl << abort(FatalError);
+        }
+    }
+
+    return HsMixture;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::CompositionModel<CloudType>::Hc
+(
+    const label phaseI,
+    const scalarField& Y,
+    const scalar p,
+    const scalar T
+) const
+{
+    const phaseProperties& props = phaseProps_[phaseI];
+    scalar HcMixture = 0.0;
+    switch (props.phase())
+    {
+        case phaseProperties::GAS:
+        {
+            forAll(Y, i)
+            {
+                label gid = props.globalIds()[i];
+                HcMixture += Y[i]*mcCarrierThermo_.speciesData()[gid].Hc();
+            }
+            break;
+        }
+        case phaseProperties::LIQUID:
+        {
+            forAll(Y, i)
+            {
+                label gid = props.globalIds()[i];
+                HcMixture +=
+                    Y[i]*this->liquids().properties()[gid].h(p, 298.15);
+            }
+            break;
+        }
+        case phaseProperties::SOLID:
+        {
+            forAll(Y, i)
+            {
+                label gid = props.globalIds()[i];
+                HcMixture += Y[i]*this->solids().properties()[gid].Hf();
+            }
+            break;
+        }
+        default:
+        {
+            FatalErrorIn
+            (
+                "Foam::scalar Foam::CompositionModel<CloudType>::Hc"
+                "("
+                "    const label, "
+                "    const scalarField&, "
+                "    const scalar, "
+                "    const scalar"
+                ") const"
+            )   << "Unknown phase enumeration" << nl << abort(FatalError);
+        }
+    }
+
+    return HcMixture;
+}
+
+
 template<class CloudType>
 Foam::scalar Foam::CompositionModel<CloudType>::cp
 (

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

@@ -226,7 +226,7 @@ public:
 
         // Evaluation
 
-            //- Return enthalpy for the phase phaseI
+            //- Return total enthalpy for the phase phaseI
             virtual scalar H
             (
                 const label phaseI,
@@ -235,6 +235,24 @@ public:
                 const scalar T
             ) const;
 
+            //- Return sensible enthalpy for the phase phaseI
+            virtual scalar Hs
+            (
+                const label phaseI,
+                const scalarField& Y,
+                const scalar p,
+                const scalar T
+            ) const;
+
+            //- Return chemical enthalpy for the phase phaseI
+            virtual scalar Hc
+            (
+                const label phaseI,
+                const scalarField& Y,
+                const scalar p,
+                const scalar T
+            ) const;
+
             //- Return specific heat caoacity for the phase phaseI
             virtual scalar cp
             (