Merge branch 'lagrangianDev' into WBM

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

@@ -96,6 +96,10 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
     constProps_(particleProperties_),
     parcelTypeId_(readLabel(particleProperties_.lookup("parcelTypeId"))),
     coupled_(particleProperties_.lookup("coupled")),
+    cellValueSourceCorrection_
+    (
+        particleProperties_.lookup("cellValueSourceCorrection")
+    ),
     rndGen_(label(0)),
     rho_(rho),
     U_(U),

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

@@ -105,7 +105,11 @@ class KinematicCloud
 
         //- Flag to indicate whether parcels are coupled to the carrier phase
         //  i.e. whether or not to generate source terms for carrier phase
-        Switch coupled_;
+        const Switch coupled_;
+
+        //- Flag to correct cell values with latest transfer information
+        //  during the lagrangian timestep
+        const Switch cellValueSourceCorrection_;
 
         //- Random number generator - used by some injection routines
         Random rndGen_;
@@ -224,6 +228,9 @@ public:
             //- Return coupled flag
             inline const Switch coupled() const;
 
+            //- Return cell value correction flag
+            inline const Switch cellValueSourceCorrection() const;
+
             //- Return refernce to the random object
             inline Random& rndGen();
 

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H

@@ -72,6 +72,14 @@ inline const Foam::Switch Foam::KinematicCloud<ParcelType>::coupled() const
 }
 
 
+template <class ParcelType>
+inline const Foam::Switch
+Foam::KinematicCloud<ParcelType>::cellValueSourceCorrection() const
+{
+    return cellValueSourceCorrection_;
+}
+
+
 template<class ParcelType>
 inline const Foam::volScalarField&
 Foam::KinematicCloud<ParcelType>::rho() const

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

@@ -72,14 +72,14 @@ Foam::ReactingCloud<ParcelType>::ReactingCloud
     const volVectorField& U,
     const dimensionedVector& g,
     hCombustionThermo& thermo,
-    PtrList<specieReactingProperties>& gases
+    PtrList<specieReactingProperties>& carrierSpecies
 )
 :
     ThermoCloud<ParcelType>(cloudType, rho, U, g, thermo),
     reactingCloud(),
     constProps_(this->particleProperties()),
     carrierThermo_(thermo),
-    gases_(gases),
+    carrierSpecies_(carrierSpecies),
     compositionModel_
     (
         CompositionModel<ReactingCloud<ParcelType> >::New

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

@@ -80,7 +80,7 @@ class ReactingCloud
         hCombustionThermo& carrierThermo_;
 
         //- Gas phase properties
-        PtrList<specieReactingProperties>& gases_;
+        PtrList<specieReactingProperties>& carrierSpecies_;
 
 
         // References to the cloud sub-models
@@ -127,7 +127,7 @@ public:
             const volVectorField& U,
             const dimensionedVector& g,
             hCombustionThermo& thermo,
-            PtrList<specieReactingProperties>& gases
+            PtrList<specieReactingProperties>& carrierSpecies
         );
 
 
@@ -154,7 +154,7 @@ public:
             inline hCombustionThermo& carrierThermo();
 
             //- Gas phase properties
-            inline const PtrList<specieReactingProperties>& gases() const;
+            inline const PtrList<specieReactingProperties>& carrierSpecies() const;
 
 
             // Sub-models

src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudI.H

@@ -52,9 +52,9 @@ Foam::ReactingCloud<ParcelType>::carrierThermo()
 
 template<class ParcelType>
 inline const Foam::PtrList<Foam::specieReactingProperties>&
-Foam::ReactingCloud<ParcelType>::gases() const
+Foam::ReactingCloud<ParcelType>::carrierSpecies() const
 {
-    return gases_;
+    return carrierSpecies_;
 }
 
 

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

@@ -79,10 +79,10 @@ Foam::ReactingMultiphaseCloud<ParcelType>::ReactingMultiphaseCloud
     const volVectorField& U,
     const dimensionedVector& g,
     hCombustionThermo& thermo,
-    PtrList<specieReactingProperties>& gases
+    PtrList<specieReactingProperties>& carrierSpecies
 )
 :
-    ReactingCloud<ParcelType>(cloudType, rho, U, g, thermo, gases),
+    ReactingCloud<ParcelType>(cloudType, rho, U, g, thermo, carrierSpecies),
     reactingMultiphaseCloud(),
     constProps_(this->particleProperties()),
     devolatilisationModel_

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

@@ -121,7 +121,7 @@ public:
             const volVectorField& U,
             const dimensionedVector& g,
             hCombustionThermo& thermo,
-            PtrList<specieReactingProperties>& gases
+            PtrList<specieReactingProperties>& carrierSpecies
         );
 
 

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

@@ -102,11 +102,25 @@ Foam::ThermoCloud<ParcelType>::ThermoCloud
         )
     ),
     radiation_(this->particleProperties().lookup("radiation")),
-    hTrans_
+    hsTrans_
     (
         IOobject
         (
-            this->name() + "hTrans",
+            this->name() + "hsTrans",
+            this->db().time().timeName(),
+            this->db(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        this->mesh(),
+        dimensionedScalar("zero", dimensionSet(1, 2, -2, 0, 0), 0.0)
+    ),
+    hcTrans_
+    (
+        IOobject
+        (
+            this->name() + "hcTrans",
             this->db().time().timeName(),
             this->db(),
             IOobject::NO_READ,
@@ -132,7 +146,8 @@ template<class ParcelType>
 void Foam::ThermoCloud<ParcelType>::resetSourceTerms()
 {
     KinematicCloud<ParcelType>::resetSourceTerms();
-    hTrans_.field() = 0.0;
+    hsTrans_.field() = 0.0;
+    hcTrans_.field() = 0.0;
 }
 
 

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

@@ -99,8 +99,14 @@ class ThermoCloud
 
         // Sources
 
-            //- Enthalpy transfer
+            //- Sensible enthalpy transfer
-            DimensionedField<scalar, volMesh> hTrans_;
+            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_;
 
 
     // Private Member Functions
@@ -173,10 +179,19 @@ public:
 
                 // Enthalpy
 
-                    //- Return reference to enthalpy source
+                    //- Return reference to sensible enthalpy source
-                    inline DimensionedField<scalar, volMesh>& hTrans();
+                    inline DimensionedField<scalar, volMesh>& hsTrans();
 
-                    //- return tmp enthalpy source term - fully explicit
+                    //- Return tmp total sensible enthalpy source term
+                    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

@@ -77,9 +77,85 @@ inline bool Foam::ThermoCloud<ParcelType>::radiation() const
 
 template<class ParcelType>
 inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
-Foam::ThermoCloud<ParcelType>::hTrans()
+Foam::ThermoCloud<ParcelType>::hsTrans()
 {
-    return hTrans_;
+    return 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;
 }
 
 
@@ -111,7 +187,7 @@ Foam::ThermoCloud<ParcelType>::Sh() const
     );
 
     scalarField& Sh = tSh().field();
-    Sh = hTrans_/(this->mesh().V()*this->db().time().deltaT());
+    Sh = (hsTrans_ + hcTrans_)/(this->mesh().V()*this->db().time().deltaT());
 
     return tSh;
 }

src/lagrangian/intermediate/clouds/derived/basicReactingCloud/basicReactingCloud.C

@@ -43,7 +43,7 @@ Foam::basicReactingCloud::basicReactingCloud
     const volVectorField& U,
     const dimensionedVector& g,
     hCombustionThermo& thermo,
-    PtrList<specieReactingProperties>& gases
+    PtrList<specieReactingProperties>& carrierSpecies
 )
 :
     ReactingCloud<basicReactingParcel>
@@ -53,7 +53,7 @@ Foam::basicReactingCloud::basicReactingCloud
         U,
         g,
         thermo,
-        gases
+        carrierSpecies
     )
 {
     basicReactingParcel::readFields(*this);

src/lagrangian/intermediate/clouds/derived/basicReactingCloud/basicReactingCloud.H

@@ -77,7 +77,7 @@ public:
             const volVectorField& U,
             const dimensionedVector& g,
             hCombustionThermo& thermo,
-            PtrList<specieReactingProperties>& gases
+            PtrList<specieReactingProperties>& carrierSpecies
         );
 
 

src/lagrangian/intermediate/clouds/derived/basicReactingMultiphaseCloud/basicReactingMultiphaseCloud.C

@@ -43,7 +43,7 @@ Foam::basicReactingMultiphaseCloud::basicReactingMultiphaseCloud
     const volVectorField& U,
     const dimensionedVector& g,
     hCombustionThermo& thermo,
-    PtrList<specieReactingProperties>& gases
+    PtrList<specieReactingProperties>& carrierSpecies
 )
 :
     ReactingMultiphaseCloud<basicReactingMultiphaseParcel>
@@ -53,7 +53,7 @@ Foam::basicReactingMultiphaseCloud::basicReactingMultiphaseCloud
         U,
         g,
         thermo,
-        gases
+        carrierSpecies
     )
 {
     basicReactingMultiphaseParcel::readFields(*this);

src/lagrangian/intermediate/clouds/derived/basicReactingMultiphaseCloud/basicReactingMultiphaseCloud.H

@@ -77,7 +77,7 @@ public:
             const volVectorField& U,
             const dimensionedVector& g,
             hCombustionThermo& thermo,
-            PtrList<specieReactingProperties>& gases
+            PtrList<specieReactingProperties>& carrierSpecies
         );
 
 

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

@@ -31,7 +31,7 @@ License
 
 template<class ParcelType>
 template<class TrackData>
-void Foam::KinematicParcel<ParcelType>::updateCellQuantities
+void Foam::KinematicParcel<ParcelType>::setCellValues
 (
     TrackData& td,
     const scalar dt,
@@ -43,7 +43,7 @@ void Foam::KinematicParcel<ParcelType>::updateCellQuantities
     {
         WarningIn
         (
-            "void Foam::KinematicParcel<ParcelType>::updateCellQuantities"
+            "void Foam::KinematicParcel<ParcelType>::setCellValues"
             "("
                 "TrackData&, "
                 "const scalar, "
@@ -57,9 +57,6 @@ void Foam::KinematicParcel<ParcelType>::updateCellQuantities
 
     Uc_ = td.UInterp().interpolate(this->position(), cellI);
 
-    // Apply correction to cell velocity to account for momentum transfer
-    Uc_ += td.cloud().UTrans()[cellI]/(massCell(cellI));
-
     muc_ = td.muInterp().interpolate(this->position(), cellI);
 
     // Apply dispersion components to carrier phase velocity
@@ -75,6 +72,19 @@ void Foam::KinematicParcel<ParcelType>::updateCellQuantities
 }
 
 
+template<class ParcelType>
+template<class TrackData>
+void Foam::KinematicParcel<ParcelType>::cellValueSourceCorrection
+(
+    TrackData& td,
+    const scalar dt,
+    const label cellI
+)
+{
+    Uc_ += td.cloud().UTrans()[cellI]/massCell(cellI);
+}
+
+
 template<class ParcelType>
 template<class TrackData>
 void Foam::KinematicParcel<ParcelType>::calc
@@ -183,8 +193,8 @@ bool Foam::KinematicParcel<ParcelType>::move(TrackData& td)
         // Set the Lagrangian time-step
         scalar dt = min(dtMax, tEnd);
 
-        // Remember which cell the Parcel is in
+        // Remember which cell the Parcel is in since this will change if a
-        // since this will change if a face is hit
+        // face is hit
         label cellI = p.cell();
 
         dt *= p.trackToFace(p.position() + dt*U_, td);
@@ -192,12 +202,17 @@ bool Foam::KinematicParcel<ParcelType>::move(TrackData& td)
         tEnd -= dt;
         p.stepFraction() = 1.0 - tEnd/deltaT;
 
-        // Update cell based properties
-        p.updateCellQuantities(td, dt, cellI);
-
         // Avoid problems with extremely small timesteps
         if (dt > ROOTVSMALL)
         {
+            // Update cell based properties
+            p.setCellValues(td, dt, cellI);
+
+            if (td.cloud().cellValueSourceCorrection())
+            {
+                p.cellValueSourceCorrection(td, dt, cellI);
+            }
+
             p.calc(td, dt, cellI);
         }
 

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

@@ -372,9 +372,18 @@ public:
 
         // Main calculation loop
 
-            //- Update cell based quantities
+            //- Set cell values
             template<class TrackData>
-            void updateCellQuantities
+            void setCellValues
+            (
+                TrackData& td,
+                const scalar dt,
+                const label cellI
+            );
+
+            //- Correct cell values using latest transfer information
+            template<class TrackData>
+            void cellValueSourceCorrection
             (
                 TrackData& td,
                 const scalar dt,

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

@@ -129,14 +129,56 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::updateMassFractions
 
 template<class ParcelType>
 template<class TrackData>
-void Foam::ReactingMultiphaseParcel<ParcelType>::updateCellQuantities
+void Foam::ReactingMultiphaseParcel<ParcelType>::setCellValues
 (
     TrackData& td,
     const scalar dt,
     const label cellI
 )
 {
-    ReactingParcel<ParcelType>::updateCellQuantities(td, dt, cellI);
+    ReactingParcel<ParcelType>::setCellValues(td, dt, cellI);
+}
+
+
+template<class ParcelType>
+template<class TrackData>
+void Foam::ReactingMultiphaseParcel<ParcelType>::cellValueSourceCorrection
+(
+    TrackData& td,
+    const scalar dt,
+    const label cellI
+)
+{
+    scalar massCell = this->massCell(cellI);
+
+    scalar addedMass = 0.0;
+    forAll(td.cloud().rhoTrans(), i)
+    {
+        addedMass += td.cloud().rhoTrans(i)[cellI];
+    }
+
+    this->rhoc_ += addedMass/td.cloud().pMesh().cellVolumes()[cellI];
+
+    scalar massCellNew = massCell + addedMass;
+    this->Uc_ += td.cloud().UTrans()[cellI]/massCellNew;
+
+    scalar cpEff = 0;
+    if (addedMass > ROOTVSMALL)
+    {
+        forAll(td.cloud().rhoTrans(), i)
+        {
+            scalar Y = td.cloud().rhoTrans(i)[cellI]/addedMass;
+            cpEff += Y*td.cloud().carrierSpecies()[i].Cp(this->Tc_);
+        }
+    }
+    const scalar cpc = td.cpInterp().psi()[cellI];
+    this->cpc_ = (massCell*cpc + addedMass*cpEff)/massCellNew;
+
+    const scalar fCarrier = -1.0/td.constProps().hRetentionCoeff();
+    const scalar dh =
+        td.cloud().hsTrans()[cellI] + fCarrier*td.cloud().hcTrans[cellI];
+
+    this->Tc_ += dh/(this->cpc_*massCellNew);
 }
 
 
@@ -226,10 +268,10 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     scalarField dMassSRGas(YGas_.size(), 0.0);
     scalarField dMassSRLiquid(YLiquid_.size(), 0.0);
     scalarField dMassSRSolid(YSolid_.size(), 0.0);
-    scalarField dMassSRCarrier(td.cloud().gases().size(), 0.0);
+    scalarField dMassSRCarrier(td.cloud().carrierSpecies().size(), 0.0);
 
     // Return enthalpy source and calc mass transfer(s) due to surface reaction
-    scalar HReaction =
+    scalar hReaction =
         calcSurfaceReactions
         (
             td,
@@ -250,11 +292,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
             dMassSRCarrier
         );
 
-    // Heat of reaction split between component retained by particle
+    // Heat of reaction retained by particle
-    const scalar ShSR = td.constProps().hRetentionCoeff()*HReaction;
+    const scalar ShSR = td.constProps().hRetentionCoeff()*hReaction;
-
-    // ...and component added to the carrier phase
-    const scalar ShSRc = (1.0 - td.constProps().hRetentionCoeff())*HReaction;
 
 
     // Update component mass fractions
@@ -310,12 +349,14 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
             label id = td.cloud().composition().localToGlobalGasId(LIQ, i);
             td.cloud().rhoTrans(id)[cellI] += np0*dMassLiquid[i];
         }
-//        // No mapping between solid components and carrier phase
+/*
-//        forAll(YSolid_, i)
+        // No mapping between solid components and carrier phase
-//        {
+        forAll(YSolid_, i)
-//            label id = td.cloud().composition().localToGlobalGasId(SLD, i);
+        {
-//            td.cloud().rhoTrans(id)[cellI] += np0*dMassSolid[i];
+            label id = td.cloud().composition().localToGlobalGasId(SLD, i);
-//        }
+            td.cloud().rhoTrans(id)[cellI] += np0*dMassSolid[i];
+        }
+*/
         forAll(dMassSRCarrier, i)
         {
             td.cloud().rhoTrans(i)[cellI] += np0*dMassSRCarrier[i];
@@ -324,8 +365,11 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         // Update momentum transfer
         td.cloud().UTrans()[cellI] += np0*(mass0*U0 - mass1*U1);
 
-        // Update enthalpy transfer
+        // Update sensible enthalpy transfer
-        td.cloud().hTrans()[cellI] += np0*(mass0*H0 - (mass1*H1 + ShSRc));
+        td.cloud().hsTrans()[cellI] += np0*(mass0*H0 - mass1*H1);
+
+        // Update chemical enthalpy transfer
+        td.cloud().hcTrans()[cellI] -= np0*ShSR;
     }
 
 
@@ -346,21 +390,20 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
             }
             forAll(YLiquid_, i)
             {
-                label id =
+                label id = td.cloud().composition().localToGlobalGasId(LIQ, i);
-                    td.cloud().composition().localToGlobalGasId(LIQ, i);
                 td.cloud().rhoTrans(id)[cellI] +=
                     np0*mass1*YMix[LIQ]*YLiquid_[i];
             }
-//            // No mapping between solid components and carrier phase
+/*
-//            forAll(YSolid_, i)
+            // No mapping between solid components and carrier phase
-//            {
+            forAll(YSolid_, i)
-//                label id =
+            {
-//                    td.cloud().composition().localToGlobalGasId(SLD, i);
+                label id = td.cloud().composition().localToGlobalGasId(SLD, i);
-//                td.cloud().rhoTrans(id)[cellI] +=
+                td.cloud().rhoTrans(id)[cellI] +=
-//                    np0*mass1*YMix[SLD]*YSolid_[i];
+                    np0*mass1*YMix[SLD]*YSolid_[i];
-//            }
+            }
-
+*/
-            td.cloud().hTrans()[cellI] += np0*mass1*H1;
+            td.cloud().hsTrans()[cellI] += np0*mass1*H1;
             td.cloud().UTrans()[cellI] += np0*mass1*U1;
         }
     }
@@ -466,7 +509,7 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
     }
 
     // Update surface reactions
-    const scalar HReaction = td.cloud().surfaceReaction().calculate
+    const scalar hReaction = td.cloud().surfaceReaction().calculate
     (
         dt,
         cellI,
@@ -493,7 +536,7 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
        *(sum(dMassSRGas) + sum(dMassSRLiquid) + sum(dMassSRSolid))
     );
 
-    return HReaction;
+    return hReaction;
 }
 
 
@@ -502,4 +545,3 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
 #include "ReactingMultiphaseParcelIO.C"
 
 // ************************************************************************* //
-

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

@@ -333,9 +333,18 @@ public:
 
         // Main calculation loop
 
-            //- Update cell based quantities
+            //- Set cell values
             template<class TrackData>
-            void updateCellQuantities
+            void setCellValues
+            (
+                TrackData& td,
+                const scalar dt,
+                const label cellI
+            );
+
+            //- Correct cell values using latest transfer information
+            template<class TrackData>
+            void cellValueSourceCorrection
             (
                 TrackData& td,
                 const scalar dt,
@@ -393,4 +402,3 @@ public:
 #endif
 
 // ************************************************************************* //
-

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

@@ -30,21 +30,21 @@ License
 
 template<class ParcelType>
 template<class TrackData>
-void Foam::ReactingParcel<ParcelType>::updateCellQuantities
+void Foam::ReactingParcel<ParcelType>::setCellValues
 (
     TrackData& td,
     const scalar dt,
     const label cellI
 )
 {
-    ThermoParcel<ParcelType>::updateCellQuantities(td, dt, cellI);
+    ThermoParcel<ParcelType>::setCellValues(td, dt, cellI);
 
     pc_ = td.pInterp().interpolate(this->position(), cellI);
     if (pc_ < td.constProps().pMin())
     {
         WarningIn
         (
-            "void Foam::ReactingParcel<ParcelType>::updateCellQuantities"
+            "void Foam::ReactingParcel<ParcelType>::setCellValues"
             "("
                 "TrackData&, "
                 "const scalar, "
@@ -55,14 +55,44 @@ void Foam::ReactingParcel<ParcelType>::updateCellQuantities
 
         pc_ = td.constProps().pMin();
     }
+}
+
+
+template<class ParcelType>
+template<class TrackData>
+void Foam::ReactingParcel<ParcelType>::cellValueSourceCorrection
+(
+    TrackData& td,
+    const scalar dt,
+    const label cellI
+)
+{
+    scalar massCell = this->massCell(cellI);
 
-    // Apply correction to cell density to account for mass transfer
     scalar addedMass = 0.0;
     forAll(td.cloud().rhoTrans(), i)
     {
         addedMass += td.cloud().rhoTrans(i)[cellI];
     }
+
     this->rhoc_ += addedMass/td.cloud().pMesh().cellVolumes()[cellI];
+
+    scalar massCellNew = massCell + addedMass;
+    this->Uc_ += td.cloud().UTrans()[cellI]/massCellNew;
+
+    scalar cpEff = 0;
+    if (addedMass > ROOTVSMALL)
+    {
+        forAll(td.cloud().rhoTrans(), i)
+        {
+            scalar Y = td.cloud().rhoTrans(i)[cellI]/addedMass;
+            cpEff += Y*td.cloud().carrierSpecies()[i].Cp(this->Tc_);
+        }
+    }
+    const scalar cpc = td.cpInterp().psi()[cellI];
+    this->cpc_ = (massCell*cpc + addedMass*cpEff)/massCellNew;
+
+    this->Tc_ += td.cloud().hsTrans()[cellI]/(this->cpc_*massCellNew);
 }
 
 
@@ -164,8 +194,8 @@ void Foam::ReactingParcel<ParcelType>::calc
         // Update momentum transfer
         td.cloud().UTrans()[cellI] += np0*(mass0*U0 - mass1*U1);
 
-        // Update enthalpy transfer
+        // Update sensible enthalpy transfer
-        td.cloud().hTrans()[cellI] += np0*(mass0*H0 - mass1*H1);
+        td.cloud().hsTrans()[cellI] += np0*(mass0*H0 - mass1*H1);
     }
 
 
@@ -184,7 +214,7 @@ void Foam::ReactingParcel<ParcelType>::calc
                 td.cloud().rhoTrans(id)[cellI] += np0*mass1*Y_[i];
             }
             td.cloud().UTrans()[cellI] += np0*mass1*U1;
-            td.cloud().hTrans()[cellI] += np0*mass1*H1;
+            td.cloud().hsTrans()[cellI] += np0*mass1*H1;
         }
     }
 

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

@@ -276,9 +276,18 @@ public:
 
         // Main calculation loop
 
-            //- Update cell based quantities
+            //- Set cell values
             template<class TrackData>
-            void updateCellQuantities
+            void setCellValues
+            (
+                TrackData& td,
+                const scalar dt,
+                const label cellI
+            );
+
+            //- Correct cell values using latest transfer information
+            template<class TrackData>
+            void cellValueSourceCorrection
             (
                 TrackData& td,
                 const scalar dt,
@@ -333,4 +342,3 @@ public:
 #endif
 
 // ************************************************************************* //
-

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

@@ -30,28 +30,24 @@ License
 
 template<class ParcelType>
 template<class TrackData>
-void Foam::ThermoParcel<ParcelType>::updateCellQuantities
+void Foam::ThermoParcel<ParcelType>::setCellValues
 (
     TrackData& td,
     const scalar dt,
     const label cellI
 )
 {
-    KinematicParcel<ParcelType>::updateCellQuantities(td, dt, cellI);
+    KinematicParcel<ParcelType>::setCellValues(td, dt, cellI);
 
     cpc_ = td.cpInterp().interpolate(this->position(), cellI);
 
     Tc_ = td.TInterp().interpolate(this->position(), cellI);
 
-    // Apply correction to cell temperature to account for enthalpy transfer
-    scalar cpMean = td.cpInterp().psi()[cellI];
-    Tc_ += td.cloud().hTrans()[cellI]/(cpMean*this->massCell(cellI));
-
     if (Tc_ < td.constProps().TMin())
     {
         WarningIn
         (
-            "void Foam::ThermoParcel<ParcelType>::updateCellQuantities"
+            "void Foam::ThermoParcel<ParcelType>::setCellValues"
             "("
                 "TrackData&, "
                 "const scalar, "
@@ -65,6 +61,22 @@ void Foam::ThermoParcel<ParcelType>::updateCellQuantities
 }
 
 
+template<class ParcelType>
+template<class TrackData>
+void Foam::ThermoParcel<ParcelType>::cellValueSourceCorrection
+(
+    TrackData& td,
+    const scalar dt,
+    const label cellI
+)
+{
+    this->Uc_ += td.cloud().UTrans()[cellI]/this->massCell(cellI);
+
+    scalar cpMean = td.cpInterp().psi()[cellI];
+    Tc_ += td.cloud().hsTrans()[cellI]/(cpMean*this->massCell(cellI));
+}
+
+
 template<class ParcelType>
 template<class TrackData>
 void Foam::ThermoParcel<ParcelType>::calc
@@ -118,8 +130,8 @@ void Foam::ThermoParcel<ParcelType>::calc
         // Update momentum transfer
         td.cloud().UTrans()[cellI] += np0*mass0*(U0 - U1);
 
-        // Update enthalpy transfer
+        // Update sensible enthalpy transfer
-        td.cloud().hTrans()[cellI] += np0*mass0*(H0 - H1);
+        td.cloud().hsTrans()[cellI] += np0*mass0*(H0 - H1);
     }
 
     // Set new particle properties
@@ -198,4 +210,3 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
 #include "ThermoParcelIO.C"
 
 // ************************************************************************* //
-

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

@@ -286,9 +286,18 @@ public:
 
         // Main calculation loop
 
-            //- Update cell based quantities
+            //- Set cell values
             template<class TrackData>
-            void updateCellQuantities
+            void setCellValues
+            (
+                TrackData& td,
+                const scalar dt,
+                const label cellI
+            );
+
+            //- Correct cell values using latest transfer information
+            template<class TrackData>
+            void cellValueSourceCorrection
             (
                 TrackData& td,
                 const scalar dt,
@@ -343,4 +352,3 @@ public:
 #endif
 
 // ************************************************************************* //
-

src/lagrangian/intermediate/phaseProperties/phaseProperties/phaseProperties.H

@@ -100,7 +100,8 @@ private:
         //- Set global ids - liquid and solid phases
         void setGlobalIds(const wordList& globalNames);
 
-        //- Set global gas ids - attempts to map component names to global gases
+        //- Set global gas ids - attempts to map component names to global
+        //  carrier species
         void setGlobalGasIds(const PtrList<volScalarField>& YGas);
 
         //- Check the total mass fraction

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

@@ -40,7 +40,7 @@ Foam::CompositionModel<CloudType>::CompositionModel
     owner_(owner),
     coeffDict_(dict.subDict(type + "Coeffs")),
     carrierThermo_(owner.carrierThermo()),
-    gases_(owner.gases()),
+    carrierSpecies_(owner.carrierSpecies()),
     liquids_
     (
         liquidMixture::New
@@ -111,9 +111,9 @@ Foam::CompositionModel<CloudType>::carrierThermo() const
 
 template<class CloudType>
 const Foam::PtrList<Foam::specieReactingProperties>&
-Foam::CompositionModel<CloudType>::gases() const
+Foam::CompositionModel<CloudType>::carrierSpecies() const
 {
-    return gases_;
+    return carrierSpecies_;
 }
 
 
@@ -294,8 +294,8 @@ Foam::scalarField Foam::CompositionModel<CloudType>::X
             forAll(Y, i)
             {
                 label gid = props.globalIds()[i];
-                WInv += Y[i]/this->gases()[gid].W();
+                WInv += Y[i]/this->carrierSpecies()[gid].W();
-                X[i] = Y[i]/this->gases()[gid].W();
+                X[i] = Y[i]/this->carrierSpecies()[gid].W();
             }
             break;
         }
@@ -345,7 +345,7 @@ Foam::scalar Foam::CompositionModel<CloudType>::H
             forAll(Y, i)
             {
                 label gid = props.globalIds()[i];
-                HMixture += Y[i]*this->gases()[gid].H(T);
+                HMixture += Y[i]*this->carrierSpecies()[gid].H(T);
             }
             break;
         }
@@ -409,7 +409,7 @@ Foam::scalar Foam::CompositionModel<CloudType>::cp
             forAll(Y, i)
             {
                 label gid = props.globalIds()[i];
-                cpMixture += Y[i]*this->gases()[gid].Cp(T);
+                cpMixture += Y[i]*this->carrierSpecies()[gid].Cp(T);
             }
             break;
         }

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

@@ -27,7 +27,8 @@ Class
 
 Description
     Templated reacting parcel composition model class
-    Consists of gases (via thermo package), liquids and solids
+    Consists of carrier species (via thermo package), and additional liquids
+    and solids
 
 SourceFiles
     CompositionModel.C
@@ -77,12 +78,12 @@ class CompositionModel
         hCombustionThermo& carrierThermo_;
 
         //- Reference to the carrier phase species
-        const PtrList<specieReactingProperties>& gases_;
+        const PtrList<specieReactingProperties>& carrierSpecies_;
 
-        //- Global liquid properties data
+        //- Global (additional) liquid properties data
         autoPtr<liquidMixture> liquids_;
 
-        //- Global solid properties data
+        //- Global (additional) solid properties data
         autoPtr<solidMixture> solids_;
 
         //- List of phase properties
@@ -150,13 +151,13 @@ public:
 
             // Composition lists
 
-                //- Return the gases
+                //- Return the carrier species
-                const PtrList<specieReactingProperties>& gases() const;
+                const PtrList<specieReactingProperties>& carrierSpecies() const;
 
-                //- Return the global liquids
+                //- Return the global (additional) liquids
                 const liquidMixture& liquids() const;
 
-                //- Return the global solids
+                //- Return the global (additional) solids
                 const solidMixture& solids() const;
 
                 //- Return the list of phase properties

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

@@ -41,8 +41,8 @@ Foam::scalarField Foam::LiquidEvaporation<CloudType>::calcXc
     forAll(Xc, i)
     {
         scalar Y = this->owner().carrierThermo().composition().Y()[i][cellI];
-        Winv += Y/this->owner().gases()[i].W();
+        Winv += Y/this->owner().carrierSpecies()[i].W();
-        Xc[i] = Y/this->owner().gases()[i].W();
+        Xc[i] = Y/this->owner().carrierSpecies()[i].W();
     }
 
     return Xc/Winv;
@@ -153,6 +153,9 @@ void Foam::LiquidEvaporation<CloudType>::calculate
     // Reynolds number
     scalar Re = mag(Ur)*d/(nuc + ROOTVSMALL);
 
+    // film temperature evaluated using the 1/3 rule
+    scalar Tf = (2.0*T + Tc)/3.0;
+
     // calculate mass transfer of each specie in liquid
     forAll(activeLiquids_, i)
     {
@@ -160,10 +163,12 @@ void Foam::LiquidEvaporation<CloudType>::calculate
         label lid = liqToLiqMap_[i];
 
         // vapour diffusivity [m2/s]
-        scalar Dab = liquids_->properties()[lid].D(pc, T);
+        scalar Dab = liquids_->properties()[lid].D(pc, Tf);
 
         // saturation pressure for species i [pa]
-        scalar pSat = liquids_->properties()[lid].pv(pc, T);
+        // - carrier phase pressure assumed equal to the liquid vapour pressure
+        //   close to the surface
+        scalar pSat = liquids_->properties()[lid].pv(pc, Tf);
 
         // Schmidt number
         scalar Sc = nuc/(Dab + ROOTVSMALL);
@@ -175,7 +180,7 @@ void Foam::LiquidEvaporation<CloudType>::calculate
         scalar kc = Sh*Dab/(d + ROOTVSMALL);
 
         // vapour concentration at droplet surface [kgmol/m3]
-        scalar Cs = pSat/(specie::RR*T);
+        scalar Cs = pSat/(specie::RR*Tf);
 
         // vapour concentration in bulk gas [kgmol/m3]
         scalar Cinf = Xc[gid]*pc/(specie::RR*Tc);

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

@@ -27,6 +27,7 @@ Class
 
 Description
     Liquid evaporation model
+    - uses ideal gas assumption
 
 \*---------------------------------------------------------------------------*/