initial check-in with new composition structure

src/lagrangian/intermediate/Make/files

@@ -68,8 +68,12 @@ submodels/addOns/radiation/scatter/cloudScatter/cloudScatter.C
 /* integration schemes */
 IntegrationScheme/makeIntegrationSchemes.C
 
+/* phase properties */
+phaseProperties/phaseProperties/phaseProperties.C
+phaseProperties/phaseProperties/phasePropertiesIO.C
+phaseProperties/phasePropertiesList/phasePropertiesList.C
 
-/* Data entries */
+/* data entries */
 submodels/IO/DataEntry/makeDataEntries.C
 
 

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

@@ -52,9 +52,6 @@ void Foam::ReactingCloud<ParcelType>::addNewParcel
         d,
         U,
         nParticles,
-        composition().YGas0(),
-        composition().YLiquid0(),
-        composition().YSolid0(),
         composition().YMixture0(),
         constProps_
     );
@@ -120,12 +117,12 @@ Foam::ReactingCloud<ParcelType>::ReactingCloud
             (
                 IOobject
                 (
-                     this->name() + "rhoTrans" + Foam::name(i),
+                    this->name() + "rhoTrans" + Foam::name(i),
-                     this->db().time().timeName(),
+                    this->db().time().timeName(),
-                     this->db(),
+                    this->db(),
-                     IOobject::NO_READ,
+                    IOobject::NO_READ,
-                     IOobject::NO_WRITE,
+                    IOobject::NO_WRITE,
-                     false
+                    false
                 ),
                 this->mesh(),
                 dimensionedScalar("zero", dimMass, 0.0)

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

@@ -43,6 +43,10 @@ void Foam::ReactingMultiphaseCloud<ParcelType>::addNewParcel
     const scalar lagrangianDt
 )
 {
+    label idGas = this->composition().idGas();
+    label idLiquid = this->composition().idLiquid();
+    label idSolid = this->composition().idSolid();
+
     ParcelType* pPtr = new ParcelType
     (
         *this,
@@ -52,9 +56,9 @@ void Foam::ReactingMultiphaseCloud<ParcelType>::addNewParcel
         d,
         U,
         nParticles,
-        this->composition().YGas0(),
+        this->composition().Y0(idGas),
-        this->composition().YLiquid0(),
+        this->composition().Y0(idLiquid),
-        this->composition().YSolid0(),
+        this->composition().Y0(idSolid),
         this->composition().YMixture0(),
         constProps_
     );
@@ -112,9 +116,9 @@ void Foam::ReactingMultiphaseCloud<ParcelType>::resetSourceTerms()
 template<class ParcelType>
 void Foam::ReactingMultiphaseCloud<ParcelType>::evolve()
 {
-    const volScalarField& T = carrierThermo_.T();
+    const volScalarField& T = this->carrierThermo().T();
-    const volScalarField cp = carrierThermo_.Cp();
+    const volScalarField cp = this->carrierThermo().Cp();
-    const volScalarField& p = carrierThermo_.p();
+    const volScalarField& p = this->carrierThermo().p();
 
     autoPtr<interpolation<scalar> > rhoInterp = interpolation<scalar>::New
     (

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

@@ -68,12 +68,6 @@ class ReactingMultiphaseCloud
         //- Parcel constant properties
         typename ParcelType::constantProperties constProps_;
 
-        //- Thermodynamics package (combustion)
-        hCombustionThermo& carrierThermo_;
-
-        //- Gas phase properties
-        PtrList<specieReactingProperties>& gases_;
-
 
         // References to the cloud sub-models
 
@@ -85,13 +79,6 @@ class ReactingMultiphaseCloud
             devolatilisationModel_;
 
 
-
-        // Sources
-
-            //- Mass transfer fields - one per carrier phase specie
-            PtrList<DimensionedField<scalar, volMesh> > rhoTrans_;
-
-
     // Private Member Functions
 
         //- Disallow default bitwise copy construct
@@ -125,16 +112,6 @@ public:
 
         // Access
 
-            //- Return const access to carrier phase thermo package
-            inline const hCombustionThermo& carrierThermo() const;
-
-            //- Return access to carrier phase thermo package
-            inline hCombustionThermo& carrierThermo();
-
-            //- Gas phase properties
-            inline const PtrList<specieReactingProperties>& gases() const;
-
-
             // Sub-models
 
                 //- Return reference to devolatilisation model
@@ -145,29 +122,6 @@ public:
                 devolatilisation() const;
 
 
-            // Sources
-
-                //- Mass
-
-                    //- Return reference to mass source for field i
-                    inline DimensionedField<scalar, volMesh>&
-                        rhoTrans(const label i);
-
-                    //- Return reference to mass source fields
-                    inline PtrList<DimensionedField<scalar, volMesh> >&
-                        rhoTrans();
-
-                    //- Return tmp mass source for field i
-                    //  Fully explicit
-                    inline tmp<DimensionedField<scalar, volMesh> >
-                        Srho1(const label i) const;
-
-                    //- Return tmp total mass source for carrier phase
-                    //  Fully explicit
-                    inline tmp<DimensionedField<scalar, volMesh> >
-                        Srho1() const;
-
-
         // Cloud evolution functions
 
             //- Add new parcel

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

@@ -46,7 +46,7 @@ Foam::basicReactingMultiphaseCloud::basicReactingMultiphaseCloud
     PtrList<specieReactingProperties>& gases
 )
 :
-    ReactingCloud<basicReactingMultiphaseParcel>
+    ReactingMultiphaseCloud<basicReactingMultiphaseParcel>
     (
         cloudType,
         rho,

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

@@ -36,7 +36,7 @@ SourceFiles
 #ifndef basicReactingMultiphaseCloud_H
 #define basicReactingMultiphaseCloud_H
 
-#include "ReactingCloud.H"
+#include "ReactingMultiphaseCloud.H"
 #include "basicReactingMultiphaseParcel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -50,7 +50,7 @@ namespace Foam
 
 class basicReactingMultiphaseCloud
 :
-    public ReactingCloud<basicReactingMultiphaseParcel>
+    public ReactingMultiphaseCloud<basicReactingMultiphaseParcel>
 {
     // Private Member Functions
 

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

@@ -58,7 +58,11 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcCoupled
     const scalar cp0 = this->cp_;
     const scalar np0 = this->nParticle_;
     const scalar T0 = this->T_;
-    scalarList& YMix = this->YMixture_;
+    scalarField& YMix = this->YMixture_;
+
+    label idGas = td.cloud().composition().idGas();
+    label idLiquid = td.cloud().composition().idLiquid();
+    label idSolid = td.cloud().composition().idSolid();
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~
     // Initialise transfer terms
@@ -96,7 +100,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcCoupled
     // ~~~~~~~~~~~~~~~~~~~~~~
     // Calculate phase change
     // ~~~~~~~~~~~~~~~~~~~~~~
-    calcPhaseChange(td, dt, T0, T1, dMassMT);
+    calcPhaseChange(td, dt, T0, dMassMT);
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -104,7 +108,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcCoupled
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
     calcDevolatilisation(td, dt, T0, T1, dMassMT);
 
-
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Calculate surface reactions
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -132,12 +135,12 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcCoupled
     // Correct dhTrans to account for enthalpy of evolved volatiles
     dhTrans +=
         sum(dMassMT)
-       *td.cloud().composition().HGas(YGas_, 0.5*(T0 + T1));
+       *td.cloud().composition().H(idGas, YGas_, this->pc_, 0.5*(T0 + T1));
 
     // Correct dhTrans to account for enthalpy of consumed solids
     dhTrans +=
         sum(dMassSR)
-       *td.cloud().composition().HSolid(YSolid_, 0.5*(T0 + T1));
+       *td.cloud().composition().H(idSolid, YSolid_, this->pc_, 0.5*(T0 + T1));
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~
@@ -179,8 +182,10 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcCoupled
         td.cloud().hTrans()[celli] +=
             np0*mass1
            *(
-                YMix[0]*td.cloud().composition().HGas(YGas_, T1)
+                YMix[0]
-              + YMix[2]*td.cloud().composition().HSolid(YSolid_, T1)
+               *td.cloud().composition().H(idGas, YGas_, this->pc_, T1)
+              + YMix[2]
+               *td.cloud().composition().H(idSolid, YSolid_, this->pc_, T1)
             );
         td.cloud().UTrans()[celli] += np0*mass1*U1;
     }
@@ -192,9 +197,9 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcCoupled
         this->U_ = U1;
         this->T_ = T1;
         this->cp_ =
-            YMix[0]*td.cloud().composition().cpGas(YGas_, T1)
+            YMix[0]*td.cloud().composition().cp(idGas, YGas_, this->pc_, T1)
-          + YMix[1]*td.cloud().composition().cpLiquid(YLiquid_, this->pc_, T1)
+          + YMix[1]*td.cloud().composition().cp(idLiquid, YLiquid_, this->pc_, T1)
-          + YMix[2]*td.cloud().composition().cpSolid(YSolid_);
+          + YMix[2]*td.cloud().composition().cp(idSolid, YSolid_, this->pc_, T1);
 
         // Update particle density or diameter
         if (td.constProps().constantVolume())
@@ -224,7 +229,11 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcUncoupled
     const scalar T0 = this->T_;
     const scalar mass0 = this->mass();
     const scalar cp0 = this->cp_;
-    scalarList& YMix = this->YMixture_;
+    scalarField& YMix = this->YMixture_;
+
+    label idGas = td.cloud().composition().idGas();
+    label idLiquid = td.cloud().composition().idLiquid();
+    label idSolid = td.cloud().composition().idSolid();
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~
     // Initialise transfer terms
@@ -262,7 +271,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcUncoupled
     // ~~~~~~~~~~~~~~~~~~~~~~
     // Calculate phase change
     // ~~~~~~~~~~~~~~~~~~~~~~
-    calcPhaseChange(td, dt, T0, T1, dMassMT);
+    calcPhaseChange(td, dt, T0, dMassMT);
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -285,9 +294,9 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcUncoupled
 
     // New specific heat capacity
     const scalar cp1 =
-        YMix[0]*td.cloud().composition().cpGas(YGas_, T1)
+        YMix[0]*td.cloud().composition().cp(idGas, YGas_, this->pc_, T1)
-      + YMix[1]*td.cloud().composition().cpLiquid(YLiquid_, this->pc_, T1)
+      + YMix[1]*td.cloud().composition().cp(idLiquid, YLiquid_, this->pc_, T1)
-      + YMix[2]*td.cloud().composition().cpSolid(YSolid_);
+      + YMix[2]*td.cloud().composition().cp(idSolid, YSolid_, this->pc_, T1);
 
     // Add retained enthalpy to particle
     T1 += dhRet/(mass0*0.5*(cp0 + cp1));
@@ -364,7 +373,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
 
     // Determine mass to add to carrier phase
     const scalar mass = this->mass();
-    scalarList& YMix = this->YMixture_;
+    scalarField& YMix = this->YMixture_;
     const scalar dMassTot = td.cloud().devolatilisation().calculate
     (
         dt,
@@ -382,9 +391,10 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
     YMix[2] = 1.0 - YMix[0] - YMix[1];
 
     // Add to cummulative mass transfer
+    label idGas = td.cloud().composition().idGas();
     forAll (YGas_, i)
     {
-        label id = td.cloud().composition().gasGlobalIds()[i];
+        label id = td.cloud().composition().globalIds(idGas)[i];
 
         // Volatiles mass transfer
         scalar volatileMass = YGas_[i]*dMassTot;

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

@@ -40,6 +40,7 @@ SourceFiles
 #define ReactingMultiphaseParcel_H
 
 #include "ReactingParcel.H"
+#include "ReactingMultiphaseCloud.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -110,6 +111,9 @@ public:
     {
         // Private data
 
+            //- Reference to the cloud containing this particle
+            ReactingMultiphaseCloud<ParcelType>& cloud_;
+
             //- Particle constant properties
             const constantProperties& constProps_;
 
@@ -121,7 +125,7 @@ public:
             //- Construct from components
             inline trackData
             (
-                ReactingCloud<ParcelType>& cloud,
+                ReactingMultiphaseCloud<ParcelType>& cloud,
                 const constantProperties& constProps,
                 const interpolation<scalar>& rhoInterp,
                 const interpolation<vector>& UInterp,
@@ -135,6 +139,9 @@ public:
 
         // Member functions
 
+            //- Return access to the owner cloud
+            inline ReactingMultiphaseCloud<ParcelType>& cloud();
+
             //- Return const access to the constant properties
             inline const constantProperties& constProps() const;
     };
@@ -201,7 +208,7 @@ public:
         //- Construct from components
         inline ReactingMultiphaseParcel
         (
-            ReactingCloud<ParcelType>& owner,
+            ReactingMultiphaseCloud<ParcelType>& owner,
             const label typeId,
             const vector& position,
             const label celli,
@@ -292,10 +299,13 @@ public:
         // I-O
 
             //- Read
-            static void readFields(ReactingCloud<ParcelType>& c);
+            static void readFields(ReactingMultiphaseCloud<ParcelType>& c);
 
             //- Write
-            static void writeFields(const ReactingCloud<ParcelType>& c);
+            static void writeFields
+            (
+                const ReactingMultiphaseCloud<ParcelType>& c
+            );
 
 
     // Ostream Operator

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

@@ -43,7 +43,7 @@ constantProperties
 template<class ParcelType>
 inline Foam::ReactingMultiphaseParcel<ParcelType>::trackData::trackData
 (
-    ReactingCloud<ParcelType>& cloud,
+    ReactingMultiphaseCloud<ParcelType>& cloud,
     const constantProperties& constProps,
     const interpolation<scalar>& rhoInterp,
     const interpolation<vector>& UInterp,
@@ -65,14 +65,16 @@ inline Foam::ReactingMultiphaseParcel<ParcelType>::trackData::trackData
         CpInterp,
         pInterp,
         g
-    )
+    ),
+    cloud_(cloud),
+    constProps_(constProps)
 {}
 
 
 template<class ParcelType>
 inline Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
 (
-    ReactingCloud<ParcelType>& owner,
+    ReactingMultiphaseCloud<ParcelType>& owner,
     const label typeId,
     const vector& position,
     const label celli,
@@ -132,6 +134,14 @@ Foam::ReactingMultiphaseParcel<ParcelType>::constantProperties::Ldevol() const
 
 // * * * * * * * * * * * trackData Member Functions  * * * * * * * * * * * * //
 
+template<class ParcelType>
+inline Foam::ReactingMultiphaseCloud<ParcelType>&
+Foam::ReactingMultiphaseParcel<ParcelType>::trackData::cloud()
+{
+    return cloud_;
+}
+
+
 template<class ParcelType>
 inline const typename Foam::ReactingMultiphaseParcel<ParcelType>::
 constantProperties&

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

@@ -44,32 +44,34 @@ Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
 {
     if (readFields)
     {
-        const ReactingCloud<ParcelType>& cR =
+        const ReactingMultiphaseCloud<ParcelType>& cR =
-            dynamic_cast<const ReactingCloud<ParcelType>& >(cloud);
+            dynamic_cast<const ReactingMultiphaseCloud<ParcelType>& >(cloud);
 
-        const label nGas = cR.composition().gasNames().size();
+        const label idGas = cR.composition().idGas();
-        const label nLiquid = cR.composition().liquidNames().size();
+        const label nGas = cR.composition().componentNames(idGas).size();
-        const label nSolid = cR.composition().solidNames().size();
+        const label idLiquid = cR.composition().idLiquid();
+        const label nLiquid = cR.composition().componentNames(idLiquid).size();
+        const label idSolid = cR.composition().idGas();
+        const label nSolid = cR.composition().componentNames(idSolid).size();
 
         YGas_.setSize(nGas);
         YLiquid_.setSize(nLiquid);
         YSolid_.setSize(nSolid);
 
-        const scalarField& YMix = this->YMixture_;
         if (is.format() == IOstream::ASCII)
         {
             is >> YGas_ >> YLiquid_ >> YSolid_;
-            YGas_ /= YMix[0] + VSMALL;
-            YLiquid_ /= YMix[1] + VSMALL;
-            YSolid_ /= YMix[2] + VSMALL;
         }
         else
         {
             is >> YGas_ >> YLiquid_ >> YSolid_;
-            YGas_ /= YMix[0] + VSMALL;
-            YLiquid_ /= YMix[1] + VSMALL;
-            YSolid_ /= YMix[2] + VSMALL;
         }
+
+        // scale the mass fractions
+        const scalarField& YMix = this->YMixture_;
+        YGas_ /= YMix[0] + ROOTVSMALL;
+        YLiquid_ /= YMix[1] + ROOTVSMALL;
+        YSolid_ /= YMix[2] + ROOTVSMALL;
     }
 
     // Check state of Istream
@@ -88,7 +90,7 @@ Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
 template<class ParcelType>
 void Foam::ReactingMultiphaseParcel<ParcelType>::readFields
 (
-    ReactingCloud<ParcelType>& c
+    ReactingMultiphaseCloud<ParcelType>& c
 )
 {
     if (!c.size())
@@ -99,20 +101,20 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::readFields
     ReactingParcel<ParcelType>::readFields(c);
 
     // Get names and sizes for each Y...
-    const wordList gasNames = c.composition().gasNames();
+    const label idGas = c.composition().idGas();
-    const wordList liquidNames = c.composition().liquidNames();
+    const wordList gasNames = c.composition().componentNames(idGas);
-    const wordList solidNames = c.composition().solidNames();
+    const label idLiquid = c.composition().idLiquid();
-    const label nGas = gasNames.size();
+    const wordList liquidNames = c.composition().componentNames(idLiquid);
-    const label nLiquid = liquidNames.size();
+    const label idSolid = c.composition().idSolid();
-    const label nSolid = solidNames.size();
+    const wordList solidNames = c.composition().componentNames(idSolid);
 
     // Set storage for each Y... for each parcel
     forAllIter(typename Cloud<ParcelType>, c, iter)
     {
         ReactingMultiphaseParcel<ParcelType>& p = iter();
-        p.YGas_.setSize(nGas, 0.0);
+        p.YGas_.setSize(gasNames.size(), 0.0);
-        p.YLiquid_.setSize(nLiquid, 0.0);
+        p.YLiquid_.setSize(liquidNames.size(), 0.0);
-        p.YSolid_.setSize(nSolid, 0.0);
+        p.YSolid_.setSize(solidNames.size(), 0.0);
     }
 
     // Populate YGas for each parcel
@@ -127,7 +129,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::readFields
         forAllIter(typename Cloud<ParcelType>, c, iter)
         {
             ReactingMultiphaseParcel<ParcelType>& p = iter();
-            p.YGas_[j] = YGas[i++]/p.YMixture()[0];
+            p.YGas_[j] = YGas[i++]/(p.YMixture()[0] + ROOTVSMALL);
         }
     }
     // Populate YLiquid for each parcel
@@ -142,7 +144,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::readFields
         forAllIter(typename Cloud<ParcelType>, c, iter)
         {
             ReactingMultiphaseParcel<ParcelType>& p = iter();
-            p.YLiquid_[j] = YLiquid[i++]/p.YMixture()[1];
+            p.YLiquid_[j] = YLiquid[i++]/(p.YMixture()[1] + ROOTVSMALL);
         }
     }
     // Populate YSolid for each parcel
@@ -157,7 +159,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::readFields
         forAllIter(typename Cloud<ParcelType>, c, iter)
         {
             ReactingMultiphaseParcel<ParcelType>& p = iter();
-            p.YSolid_[j] = YSolid[i++]/p.YMixture()[2];
+            p.YSolid_[j] = YSolid[i++]/(p.YMixture()[2] + ROOTVSMALL);
         }
     }
 }
@@ -166,7 +168,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::readFields
 template<class ParcelType>
 void Foam::ReactingMultiphaseParcel<ParcelType>::writeFields
 (
-    const ReactingCloud<ParcelType>& c
+    const ReactingMultiphaseCloud<ParcelType>& c
 )
 {
     ReactingParcel<ParcelType>::writeFields(c);
@@ -176,7 +178,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::writeFields
     // Write the composition fractions
     if (np > 0)
     {
-        const wordList& gasNames = c.composition().gasNames();
+        const label idGas = c.composition().idGas();
+        const wordList gasNames = c.composition().componentNames(idGas);
         forAll(gasNames, j)
         {
             IOField<scalar> YGas
@@ -194,7 +197,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::writeFields
 
             YGas.write();
         }
-        const wordList& liquidNames = c.composition().liquidNames();
+        const label idLiquid = c.composition().idLiquid();
+        const wordList liquidNames = c.composition().componentNames(idLiquid);
         forAll(liquidNames, j)
         {
             IOField<scalar> YLiquid
@@ -212,7 +216,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::writeFields
 
             YLiquid.write();
         }
-        const wordList& solidNames = c.composition().solidNames();
+        const label idSolid = c.composition().idSolid();
+        const wordList solidNames = c.composition().componentNames(idSolid);
         forAll(solidNames, j)
         {
             IOField<scalar> YSolid
@@ -248,14 +253,14 @@ Foam::Ostream& Foam::operator<<
     scalarField YSolidLoc = p.YSolid()*p.YMixture()[2];
     if (os.format() == IOstream::ASCII)
     {
-        os  << static_cast<const ReactingParcel<ParcelType>& >(p)
+        os  << static_cast<const ReactingMultiphaseParcel<ParcelType>& >(p)
             << token::SPACE << YGasLoc
             << token::SPACE << YLiquidLoc
             << token::SPACE << YSolidLoc;
     }
     else
     {
-        os  << static_cast<const ReactingParcel<ParcelType>& >(p);
+        os  << static_cast<const ReactingMultiphaseParcel<ParcelType>& >(p);
         os << YGasLoc << YLiquidLoc << YSolidLoc;
     }
 

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

@@ -97,7 +97,7 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
     // ~~~~~~~~~~~~~~~~~~~~~~
     // Calculate phase change
     // ~~~~~~~~~~~~~~~~~~~~~~
-    calcPhaseChange(td, dt, T0, T1, dMassMT);
+    calcPhaseChange(td, dt, T, dMassMT);
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -117,6 +117,7 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
     T1 += dhRet/(0.5*(mass0 + mass1)*cp0);
     dhTrans -= dhRet;
 
+
     // ~~~~~~~~~~~~~~~~~~~~~~~
     // Accumulate source terms
     // ~~~~~~~~~~~~~~~~~~~~~~~
@@ -134,7 +135,6 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
     td.cloud().UCoeff()[celli] += np0*mass0*Cud;
 
     // Update enthalpy transfer
-    // - enthalpy of lost solids already accounted for
     td.cloud().hTrans()[celli] += np0*dhTrans;
 
     // Accumulate coefficient to be applied in carrier phase enthalpy coupling
@@ -229,7 +229,7 @@ void Foam::ReactingParcel<ParcelType>::calcUncoupled
     // ~~~~~~~~~~~~~~~~~~~~~~
     // Calculate phase change
     // ~~~~~~~~~~~~~~~~~~~~~~
-    calcPhaseChange(td, dt, T0, T1, dMassMT);
+    calcPhaseChange(td, dt, T, dMassMT);
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -285,16 +285,17 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
 (
     TrackData& td,
     const scalar dt,
-    const scalar T0,
+    const scalar T,
-    const scalar T1,
     scalarList& dMassMT
 )
 {
-    //    if (!td.cloud().phaseChange().active())
+    if (!td.cloud().phaseChange().active())
     {
         return;
     }
 
+    // TODO: separate treatment for boiling
+
     /*
     // Determine mass to add to carrier phase
     const scalar mass = this->mass();

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

@@ -174,8 +174,7 @@ protected:
         (
             TrackData& td,
             const scalar dt,
-            const scalar T0,
+            const scalar T,
-            const scalar T1,
             scalarList& dMassMT
         );
 

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

@@ -47,7 +47,7 @@ Foam::ReactingParcel<ParcelType>::ReactingParcel
         const ReactingCloud<ParcelType>& cR =
             dynamic_cast<const ReactingCloud<ParcelType>& >(cloud);
 
-        const label nMixture = cR.composition().compositionNames().size();
+        const label nMixture = cR.composition().phaseTypes().size();
         YMixture_.setSize(nMixture);
 
         if (is.format() == IOstream::ASCII)
@@ -102,22 +102,22 @@ void Foam::ReactingParcel<ParcelType>::readFields
     }
 
     // Get names and sizes for each Y...
-    const wordList compositionNames = c.composition().compositionNames();
+    const wordList phaseTypes = c.composition().phaseTypes();
-    const label nComposition = compositionNames.size();
+    const label nPhases = phaseTypes.size();
 
     // Set storage for each Y... for each parcel
     forAllIter(typename Cloud<ParcelType>, c, iter)
     {
         ReactingParcel<ParcelType>& p = iter();
-        p.YMixture_.setSize(nComposition, 0.0);
+        p.YMixture_.setSize(nPhases, 0.0);
     }
 
     // Populate YMixture for each parcel
-    forAll(compositionNames, j)
+    forAll(phaseTypes, j)
     {
         IOField<scalar> YMixture
         (
-            c.fieldIOobject("Y" + compositionNames[j], IOobject::MUST_READ)
+            c.fieldIOobject("Y" + phaseTypes[j], IOobject::MUST_READ)
         );
 
         label i = 0;
@@ -154,12 +154,12 @@ void Foam::ReactingParcel<ParcelType>::writeFields
     // Write the composition fractions
     if (np > 0)
     {
-        const wordList compositionNames = c.composition().compositionNames();
+        const wordList phaseTypes = c.composition().phaseTypes();
-        forAll(compositionNames, j)
+        forAll(phaseTypes, j)
         {
             IOField<scalar> YMixture
             (
-                c.fieldIOobject("Y" + compositionNames[j], IOobject::NO_READ),
+                c.fieldIOobject("Y" + phaseTypes[j], IOobject::NO_READ),
                 np
             );
 
@@ -196,6 +196,7 @@ Foam::Ostream& Foam::operator<<
         os  << static_cast<const ThermoParcel<ParcelType>& >(p);
         os.write
         (
+
             reinterpret_cast<const char*>(&p.mass0_),
             sizeof(p.mass0())
         );

src/lagrangian/intermediate/parcels/derived/basicReactingMultiphaseParcel/basicReactingMultiphaseParcel.C

@@ -40,7 +40,7 @@ namespace Foam
 
 Foam::basicReactingMultiphaseParcel::basicReactingMultiphaseParcel
 (
-    ReactingCloud<basicReactingMultiphaseParcel>& owner,
+    ReactingMultiphaseCloud<basicReactingMultiphaseParcel>& owner,
     const label typeId,
     const vector& position,
     const label celli,

src/lagrangian/intermediate/parcels/derived/basicReactingMultiphaseParcel/basicReactingMultiphaseParcel.H

@@ -63,7 +63,7 @@ public:
         //- Construct from components
         basicReactingMultiphaseParcel
         (
-             ReactingCloud<basicReactingMultiphaseParcel>& owner,
+             ReactingMultiphaseCloud<basicReactingMultiphaseParcel>& owner,
              const label typeId,
              const vector& position,
              const label celli,

src/lagrangian/intermediate/parcels/derived/basicReactingMultiphaseParcel/defineBasicReactingMultiphaseParcel.C

@@ -61,50 +61,37 @@ namespace Foam
         ReactingParcel<basicReactingMultiphaseParcel>,
         0
     );
+    defineParcelTypeNameAndDebug
+    (
+        ReactingMultiphaseParcel<basicReactingMultiphaseParcel>,
+        0
+    );
+    defineTemplateTypeNameAndDebug
+    (
+        ReactingMultiphaseParcel<basicReactingMultiphaseParcel>,
+        0
+    );
 
     defineParcelTypeNameAndDebug
     (
         KinematicCloud<basicReactingMultiphaseParcel>,
         0
     );
-//    defineTemplateTypeNameAndDebug
-//    (
-//        KinematicCloud<basicReactingMultiphaseParcel>,
-//        0
-//    );
-
     defineParcelTypeNameAndDebug
     (
         ThermoCloud<basicReactingMultiphaseParcel>,
         0
     );
-//    defineTemplateTypeNameAndDebug
-//    (
-//        ThermoCloud<basicReactingMultiphaseParcel>,
-//        0
-//    );
-
     defineParcelTypeNameAndDebug
     (
         ReactingCloud<basicReactingMultiphaseParcel>,
         0
     );
-//    defineTemplateTypeNameAndDebug
-//    (
-//        ReactingCloud<basicReactingMultiphaseParcel>,
-//        0
-//    );
-
     defineParcelTypeNameAndDebug
     (
         ReactingMultiphaseCloud<basicReactingMultiphaseParcel>,
         0
     );
-//    defineTemplateTypeNameAndDebug
-//    (
-//        ReactingMultiphaseCloud<basicReactingMultiphaseParcel>,
-//        0
-//    );
 
 };
 

src/lagrangian/intermediate/parcels/derived/basicReactingParcel/makeBasicReactingParcelCompositionModels.C

@@ -27,7 +27,7 @@ License
 #include "basicReactingParcel.H"
 #include "ReactingCloud.H"
 
-#include "SingleMixtureFraction.H"
+#include "SinglePhaseMixture.H"
 
 namespace Foam
 {
@@ -36,7 +36,7 @@ namespace Foam
     // Add instances of composition model to the table
     makeCompositionModelType
     (
-        SingleMixtureFraction,
+        SinglePhaseMixture,
         ReactingCloud,
         basicReactingParcel
     );

src/lagrangian/intermediate/phaseProperties/phaseProperties/phaseProperties.C

@@ -0,0 +1,305 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "phaseProperties.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+template<>
+const char* Foam::NamedEnum<Foam::phaseProperties::phaseType, 4>::names[] =
+{
+    "gas",
+    "liquid",
+    "solid",
+    "unknown"
+};
+
+
+const Foam::NamedEnum<Foam::phaseProperties::phaseType, 4>
+    Foam::phaseProperties::phaseTypeNames_;
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::phaseProperties::setGlobalGasIds
+(
+    const PtrList<volScalarField>& YGas
+)
+{
+    forAll(names_, i)
+    {
+        forAll (YGas, j)
+        {
+            word specieName = YGas[j].name();
+
+            if (specieName == names_[i])
+            {
+                globalIds_[i] = j;
+                break;
+            }
+        }
+        if (globalIds_[i] == -1)
+        {
+            wordList globalGasNames(YGas.size());
+
+            forAll (YGas, k)
+            {
+                globalGasNames[k] = YGas[k].name();
+            }
+
+            FatalErrorIn
+            (
+                "void phaseProperties::setGlobalGasIds\n"
+                "(\n"
+                "    const PtrList<volScalarField>&\n"
+                ")"
+            )   << "Could not find gas species " << names_[i]
+                << " in species list" <<  nl
+                << "Available species are: " << nl << globalGasNames << nl
+                << exit(FatalError);
+        }
+    }
+}
+
+
+void Foam::phaseProperties::setGlobalIds(const wordList& globalNames)
+{
+    forAll(names_, i)
+    {
+        forAll(globalNames, j)
+        {
+            if (globalNames[j] == names_[i])
+            {
+                globalIds_[i] = j;
+                break;
+            }
+        }
+        if (globalIds_[i] == -1)
+        {
+            FatalErrorIn
+            (
+                "void Foam::phaseProperties::setGlobalIds(const wordList&)"
+            )   << "Could not find specie " << names_[i]
+                << " in species list" <<  nl
+                << "Available species are: " << nl << globalNames << nl
+                << exit(FatalError);
+        }
+    }
+}
+
+
+void Foam::phaseProperties::checkTotalMassFraction() const
+{
+    scalar total = 0.0;
+    forAll(Y_, cmptI)
+    {
+        total += Y_[cmptI];
+    }
+
+    if (Y_.size() != 0 && mag(total - 1.0) > SMALL)
+    {
+        FatalErrorIn
+        (
+            "void Foam::phaseProperties::checkTotalMassFraction() const"
+        )   << "Component fractions must total to unity for phase "
+            << phaseTypeNames_[phase_] << nl
+            << "Components: " << nl << names_ << nl << exit(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::phaseProperties::phaseProperties()
+:
+    phase_(UNKNOWN),
+    names_(0),
+    Y_(0),
+    globalIds_(0)
+{}
+
+
+Foam::phaseProperties::phaseProperties(const phaseProperties& pp)
+:
+    phase_(pp.phase_),
+    names_(pp.names_),
+    Y_(pp.Y_),
+    globalIds_(pp.globalIds_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::phaseProperties::~phaseProperties()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::phaseProperties::initialiseGlobalIds
+(
+    const PtrList<volScalarField>& YGas,
+    const wordList& liquidNames,
+    const wordList& solidNames
+)
+{
+    // determine the addressing to map between components listed in the phase
+    // with those given in the (main) thermo properties
+    switch (phase_)
+    {
+        case GAS:
+        {
+            setGlobalGasIds(YGas);
+            break;
+        }
+        case LIQUID:
+        {
+            setGlobalIds(liquidNames);
+            break;
+        }
+        case SOLID:
+        {
+            setGlobalIds(solidNames);
+            break;
+        }
+        default:
+        {
+            FatalErrorIn
+            (
+                "Foam::phaseProperties::setGlobalIds\n"
+                "(\n"
+                "    const PtrList<volScalarField>&,\n"
+                "    const wordList&,\n"
+                "    const wordList&\n"
+                ")"
+            )   << "Invalid phase: " << phaseTypeNames_[phase_] << nl
+                << "    phase must be gas, liquid or solid" << nl
+                << exit(FatalError);
+        }
+    }
+}
+
+
+Foam::phaseProperties::phaseType Foam::phaseProperties::phase() const
+{
+    return phase_;
+}
+
+
+Foam::word Foam::phaseProperties::phaseTypeName() const
+{
+    return phaseTypeNames_[phase_];
+}
+
+
+const Foam::List<Foam::word>& Foam::phaseProperties::names() const
+{
+    return names_;
+}
+
+
+const Foam::word& Foam::phaseProperties::name(const label cmptI) const
+{
+    if (cmptI >= names_.size())
+    {
+        FatalErrorIn
+        (
+            "const Foam::word& Foam::phaseProperties::name"
+            "("
+                "const label cmptI"
+            ") const"
+        )   << "Requested component " << cmptI << "out of range" << nl
+            << "Available phase components:" << nl << names_ << nl
+            << exit(FatalError);
+    }
+
+    return names_[cmptI];
+}
+
+
+const Foam::scalarField& Foam::phaseProperties::Y() const
+{
+    return Y_;
+}
+
+
+Foam::scalar& Foam::phaseProperties::Y(const label cmptI)
+{
+    if (cmptI >= Y_.size())
+    {
+        FatalErrorIn
+        (
+            "const Foam::scalar& Foam::phaseProperties::Y"
+            "("
+                "const label cmptI"
+            ") const"
+        )   << "Requested component " << cmptI << "out of range" << nl
+            << "Available phase components:" << nl << names_ << nl
+            << exit(FatalError);
+    }
+
+    return Y_[cmptI];
+}
+
+
+Foam::label Foam::phaseProperties::globalId(const word& cmptName) const
+{
+    label id = this->id(cmptName);
+
+    if (id < 0)
+    {
+        return id;
+    }
+    else
+    {
+        return globalIds_[id];
+    }
+
+}
+
+
+const Foam::labelList& Foam::phaseProperties::globalIds() const
+{
+    return globalIds_;
+}
+
+
+Foam::label Foam::phaseProperties::id(const word& cmptName) const
+{
+    forAll(names_, cmptI)
+    {
+        if (names_[cmptI] == cmptName)
+        {
+            return cmptI;
+        }
+    }
+
+    return -1;
+}
+
+
+// ************************************************************************* //
+

src/lagrangian/intermediate/phaseProperties/phaseProperties/phaseProperties.C.bakNew

@@ -0,0 +1,357 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "phaseProperties.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+template<>
+const char* Foam::NamedEnum<Foam::phaseProperties::phaseType, 4>::names[] =
+{
+    "gas",
+    "liquid",
+    "solid",
+    "unknown"
+};
+
+
+const Foam::NamedEnum<Foam::phaseProperties::phaseType, 4>
+    Foam::phaseProperties::phaseTypeNames_;
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::phaseProperties::setGlobalGasIds
+(
+    const PtrList<volScalarField>& YGas
+)
+{
+    forAll(components_, i)
+    {
+        forAll (YGas, j)
+        {
+            word specieName = YGas[j].name();
+
+            if (specieName == components_[i].first())
+            {
+                globalIds_[i] = j;
+                break;
+            }
+        }
+        if (globalIds_[i] == -1)
+        {
+            wordList globalGasNames(YGas.size());
+
+            forAll (YGas, k)
+            {
+                globalGasNames[k] = YGas[k].name();
+            }
+
+            FatalErrorIn
+            (
+                "void phaseProperties::setGlobalGasIds"
+                "("
+                "    const hCombustionThermo&"
+                ")"
+            )   << "Could not find gas species " << components_[i].first()
+                << " in species list" <<  nl
+                << "Available species are: " << nl << globalGasNames << nl
+                << exit(FatalError);
+        }
+    }
+}
+
+
+void Foam::phaseProperties::setGlobalIds(const wordList& globalNames)
+{
+    forAll(components_, i)
+    {
+        forAll(globalNames, j)
+        {
+            if (globalNames[j] == components_[i].first())
+            {
+                globalIds_[i] = j;
+                break;
+            }
+        }
+        if (globalIds_[i] == -1)
+        {
+            FatalErrorIn
+            (
+                "void Foam::phaseProperties::setGlobalGasIds\n"
+                "(\n"
+                "    const PtrList<volScalarField>& YGas\n"
+                ")"
+            )   << "Could not find specie " << components_[i].first()
+                << " in species list" <<  nl
+                << "Available species are: " << nl << globalNames << nl
+                << exit(FatalError);
+        }
+    }
+}
+
+
+void Foam::phaseProperties::checkTotalMassFraction() const
+{
+    scalar total = 0.0;
+    forAll(fractions_, cmptI)
+    {
+        total += fractions_[cmptI];
+    }
+
+    if (mag(total - 1.0) > SMALL)
+    {
+        FatalErrorIn
+        (
+            "void Foam::phaseProperties::checkTotalMassFraction() const"
+        )   << "Component fractions must total to unity" << nl
+            << "Components: " << nl << components_ << nl << exit(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::phaseProperties::phaseProperties()
+:
+    phase_(UNKNOWN),
+    names_(0),
+    fractions_(0),
+    globalIds_(0)
+{}
+
+
+Foam::phaseProperties::phaseProperties(Istream& is)
+:
+    phase_(UNKNOWN),
+    names_(0),
+    fractions_(0),
+    globalIds_(0)
+{
+    is.check("Foam::phaseProperties::phaseProperties(Istream& is)");
+
+    Tuple2<word, List<Tuple2<word, scalar> > > components(is);
+
+    phase_ = phaseTypeNames_.read(components.first());
+
+    label nComponents = components.second().size();
+    names_.setSize(nComponents);
+    fractions_.setSize(nComponents);
+
+    forAll(components.second(), cmptI)
+    {
+        names_[cmptI] = components.second()[cmptI].first();
+        fractions_[cmptI] = components.second()[cmptI].second();
+    }
+
+    // initialise global ids to -1
+    globalIds_.setSize(nComponents, -1);
+
+    checkTotalMassFraction();
+}
+
+
+Foam::phaseProperties::phaseProperties(const phaseProperties& pp)
+:
+    phase_(pp.phase_),
+    names_(pp.names_),
+    fractions_(pp.fractions_),
+    globalIds_(pp.globalIds_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::phaseProperties::~phaseProperties()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::phaseProperties::initialiseGlobalIds
+(
+    const PtrList<volScalarField>& YGas,
+    const wordList& liquidNames,
+    const wordList& solidNames
+)
+{
+    // determine the addressing to map between components listed in the phase
+    // with those given in the (main) thermo properties
+    switch (phase_)
+    {
+        case GAS:
+        {
+            setGlobalGasIds(YGas);
+            break;
+        }
+        case LIQUID:
+        {
+            setGlobalIds(liquidNames);
+            break;
+        }
+        case SOLID:
+        {
+            setGlobalIds(solidNames);
+            break;
+        }
+        default:
+        {
+            FatalErrorIn
+            (
+                "Foam::phaseProperties::setGlobalIds\n"
+                "(\n"
+                "    const PtrList<specieReactingProperties>& gases,\n"
+                "    const wordList& liquidNames,\n"
+                "    const wordList& solidNames\n"
+                ")"
+            )   << "Invalid phase: " << phaseTypeNames_[phase_] << nl
+                << "    phase must be gas, liquid or solid" << nl
+                << exit(FatalError);
+        }
+}
+
+
+Foam::phaseProperties::phaseType Foam::phaseProperties::phase() const
+{
+    return phase_;
+}
+
+
+Foam::word Foam::phaseProperties::phaseTypeName() const
+{
+    return phaseTypeNames_[phase_];
+}
+
+
+const Foam::List<Foam::Tuple2<Foam::word, Foam::scalar> >&
+Foam::phaseProperties::components() const
+{
+    return components_;
+}
+
+
+const Foam::word& Foam::phaseProperties::name(const label cmptI) const
+{
+    if (cmptI >= components_.size())
+    {
+        FatalErrorIn
+        (
+            "const Foam::word& Foam::phaseProperties::name"
+            "("
+            "    const label cmptI"
+            ") const"
+        )   << "Requested component " << cmptI << "out of range" << nl
+            << "Available phase components:" << nl << components_ << nl
+            << exit(FatalError);
+    }
+
+    return components_[cmptI].first();
+}
+
+
+const Foam::wordList Foam::phaseProperties::names() const
+{
+    wordList cmptNames(components_.size());
+
+    forAll(components_, cmptI)
+    {
+        cmptNames[cmptI] = components_[cmptI].first();
+    }
+
+    return cmptNames;
+}
+
+
+Foam::label Foam::phaseProperties::id(const word& cmptName) const
+{
+    forAll(components_, cmptI)
+    {
+        if (components_[cmptI].first() == cmptName)
+        {
+            return cmptI;
+        }
+    }
+
+    return -1;
+}
+
+
+Foam::label Foam::phaseProperties::globalId(const word& cmptName) const
+{
+    label id = this->id(cmptName);
+
+    if (id < 0)
+    {
+        return id;
+    }
+    else
+    {
+        return globalIds_[id];
+    }
+
+}
+
+
+const Foam::labelList& Foam::phaseProperties::globalIds() const
+{
+    return globalIds_;
+}
+
+
+Foam::scalar& Foam::phaseProperties::Y(const label cmptI)
+{
+    if (cmptI >= components_.size())
+    {
+        FatalErrorIn
+        (
+            "const Foam::scalar& Foam::phaseProperties::Y"
+            "("
+            "    const label cmptI"
+            ") const"
+        )   << "Requested component " << cmptI << "out of range" << nl
+            << "Available phase components:" << nl << components_ << nl
+            << exit(FatalError);
+    }
+
+    return components_[cmptI].second();
+}
+
+
+const Foam::scalarList Foam::phaseProperties::Y() const
+{
+    scalarList cmptYs(components_.size(), 0.0);
+
+    forAll(cmptYs, i)
+    {
+        cmptYs[i] = components_[i].second();
+    }
+
+    return cmptYs;
+}
+
+
+// ************************************************************************* //
+

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

@@ -0,0 +1,177 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::phaseProperties
+
+Description
+    Helper class to manage multi-component phase properties
+
+SourceFiles
+    phaseProperties.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef phaseProperties_H
+#define phaseProperties_H
+
+#include "NamedEnum.H"
+#include "Tuple2.H"
+#include "PtrList.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class phaseProperties Declaration
+\*---------------------------------------------------------------------------*/
+
+class phaseProperties
+{
+public:
+
+    // Public data
+
+        //- Phase type enumeration
+        enum phaseType
+        {
+            GAS,
+            LIQUID,
+            SOLID,
+            UNKNOWN
+        };
+
+        //- Corresponding word representations for phase type enumerations
+        static const NamedEnum<phaseType, 4> phaseTypeNames_;
+
+
+private:
+
+   // Private data
+
+        //- Phase type
+        phaseType phase_;
+
+        //- List of component names
+        List<word> names_;
+
+        //- List of component mass fractions
+        scalarField Y_;
+
+        //- Global ids
+        labelList globalIds_;
+
+
+    // Private member functions
+
+        //- Set global ids - specialisation for carrier gas phases
+        void setGlobalGasIds(const PtrList<volScalarField>& YGas);
+
+        //- Set global ids - liquid and solid phases
+        void setGlobalIds(const wordList& globalNames);
+
+        //- Check the total mass fraction
+        void checkTotalMassFraction() const;
+
+
+public:
+
+    // Constructors
+
+        //- Null constructor
+        phaseProperties();
+
+        //- Construct from Istream
+        phaseProperties(Istream&);
+
+        //- Construct as copy
+        phaseProperties(const phaseProperties&);
+
+
+    //- Destructor
+    ~phaseProperties();
+
+
+    // Public member functions
+
+        //- Initialise the global ids
+        void initialiseGlobalIds
+        (
+            const PtrList<volScalarField>& YGas,
+            const wordList& liquidNames,
+            const wordList& solidNames
+        );
+
+
+        // Access
+
+            //- Return const access to the phase type
+            phaseType phase() const;
+
+            //- Return word representation of the phase type
+            word phaseTypeName() const;
+
+            //- Return the list of component names
+            const List<word>& names() const;
+
+            //- Return const access to a component name
+            const word& name(const label cmptI) const;
+
+            //- Return const access to all component mass fractions
+            const scalarField& Y() const;
+
+            //- Return non-const access to a component mass fraction
+            scalar& Y(const label cmptI);
+
+            //- Return const acccess to the global ids
+            const labelList& globalIds() const;
+
+            //- Return the global id of a component in the local list by name
+            //  Returns -1 if not found
+            label globalId(const word& cmptName) const;
+
+            //- Return the id of a component in the local list by name
+            //  Returns -1 if not found
+            label id(const word& cmptName) const;
+
+
+    // IOstream Operators
+
+        friend Istream& operator>>(Istream&, phaseProperties&);
+        friend Ostream& operator<<(Ostream&, const phaseProperties&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,177 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::phaseProperties
+
+Description
+    Helper class to manage multi-component phase properties
+
+SourceFiles
+    phaseProperties.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef phaseProperties_H
+#define phaseProperties_H
+
+#include "NamedEnum.H"
+#include "Tuple2.H"
+#include "PtrList.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class phaseProperties Declaration
+\*---------------------------------------------------------------------------*/
+
+class phaseProperties
+{
+public:
+
+    // Public data
+
+        //- Phase type enumeration
+        enum phaseType
+        {
+            GAS,
+            LIQUID,
+            SOLID,
+            UNKNOWN
+        };
+
+        //- Corresponding word representations for phase type enumerations
+        static const NamedEnum<phaseType, 4> phaseTypeNames_;
+
+
+private:
+
+   // Private data
+
+        //- Phase type
+        phaseType phase_;
+
+        //- List of component names
+        List<word> names_;
+
+        //- List of component mass fractions
+        List<scalar> fractions_;
+
+        //- Global ids
+        labelList globalIds_;
+
+
+    // Private member functions
+
+        //- Set global ids - specialisation for carrier gas phases
+        void setGlobalGasIds(const PtrList<volScalarField>& YGas);
+
+        //- Set global ids - liquid and solid phases
+        void setGlobalIds(const wordList& globalNames);
+
+        //- Check the total mass fraction
+        void checkTotalMassFraction() const;
+
+
+public:
+
+    // Constructors
+
+        //- Null constructor
+        phaseProperties();
+
+        //- Construct from Istream
+        phaseProperties(Istream&);
+
+        //- Construct as copy
+        phaseProperties(const phaseProperties&);
+
+
+    //- Destructor
+    ~phaseProperties();
+
+
+    // Public member functions
+
+        //- Initialise the global ids
+        void initialiseGlobalIds
+        (
+            const PtrList<volScalarField>& YGas,
+            const wordList& liquidNames,
+            const wordList& solidNames
+        );
+
+        //- Return const access to the phase type
+        phaseType phase() const;
+
+        //- Return word representation of the phase type
+        word phaseTypeName() const;
+
+        //- Return const access to the component name vs mass fraction list
+        const List<Tuple2<word, scalar> >& components() const;
+
+        //- Return const access to a component name
+        const word& name(const label cmptI) const;
+
+        //- Return a list of component names
+        const wordList names() const;
+
+        //- Return the id of a component in the local list by name
+        //  Returns -1 if not found
+        label id(const word& cmptName) const;
+
+        //- Return the global id of a component in the local list by name
+        //  Returns -1 if not found
+        label globalId(const word& cmptName) const;
+
+        //- Return const acccess to the global ids
+        const labelList& globalIds() const;
+
+        //- Return non-const access to a component mass fraction
+        scalar& Y(const label cmptI);
+
+        //- Return const access to all component mass fractions
+        const scalarList Y() const;
+
+
+    // IOstream Operators
+
+        friend Istream& operator>>(Istream&, phaseProperties&);
+        friend Ostream& operator<<(Ostream&, const phaseProperties&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/phaseProperties/phaseProperties/phasePropertiesIO.C

@@ -0,0 +1,145 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "phaseProperties.H"
+#include "dictionaryEntry.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+
+Foam::phaseProperties::phaseProperties(Istream& is)
+:
+    phase_(UNKNOWN),
+    names_(0),
+    Y_(0),
+    globalIds_(0)
+{
+    is.check("Foam::phaseProperties::phaseProperties(Istream& is)");
+
+    dictionaryEntry phaseInfo(dictionary::null, is);
+
+    if (!phaseInfo.isDict())
+    {
+        FatalErrorIn("Foam::phaseProperties::phaseProperties(Istream& is)")
+            << "Phase properties should be given in dictionary entries"
+            << nl << exit(FatalError);
+    }
+
+    label nComponents = phaseInfo.size();
+    names_.setSize(nComponents);
+    Y_.setSize(nComponents);
+
+    phase_ = phaseTypeNames_[phaseInfo.keyword()];
+
+    label cmptI = 0;
+    forAllConstIter(IDLList<entry>, phaseInfo, iter)
+    {
+        names_[cmptI] = iter().keyword();
+        Y_[cmptI] = readScalar(phaseInfo.lookup(names_[cmptI]));
+        cmptI++;
+    }
+
+    // initialise global ids to -1
+    globalIds_.setSize(nComponents, -1);
+
+    checkTotalMassFraction();
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Istream& Foam::operator>>(Istream& is, phaseProperties& pp)
+{
+    is.check
+    (
+        "Foam::Istream& Foam::operator>>"
+        "(Foam::Istream&, Foam::phaseProperties&)"
+    );
+
+    dictionaryEntry phaseInfo(dictionary::null, is);
+
+    if (!phaseInfo.isDict())
+    {
+        FatalErrorIn
+        (
+            "Foam::Istream& Foam::operator>>(Istream& is, phaseProperties& pp)"
+        )   << "Phase properties should be given in dictionary entries"
+            << nl << exit(FatalError);
+    }
+
+    label nComponents = phaseInfo.size();
+    pp.names_.setSize(nComponents);
+    pp.Y_.setSize(nComponents);
+
+    pp.phase_ = pp.phaseTypeNames_[phaseInfo.keyword()];
+
+    label cmptI = 0;
+    forAllConstIter(IDLList<entry>, phaseInfo, iter)
+    {
+        pp.names_[cmptI] = iter().keyword();
+        pp.Y_[cmptI] = readScalar(phaseInfo.lookup(pp.names_[cmptI]));
+        cmptI++;
+    }
+
+    // initialise global ids to -1
+    pp.globalIds_.setSize(nComponents, -1);
+
+    pp.checkTotalMassFraction();
+
+    return is;
+}
+
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const phaseProperties& pp)
+{
+    os.check
+    (
+        "Foam::Ostream& Foam::operator<<"
+        "(Foam::Ostream&, const Foam::phaseProperties&)"
+    );
+
+    os  << pp.phaseTypeNames_[pp.phase_] << nl << token::BEGIN_BLOCK << nl
+        << incrIndent;
+
+    forAll(pp.names_, cmptI)
+    {
+        os.writeKeyword(pp.names_[cmptI]) << pp.Y_[cmptI]
+            << token::END_STATEMENT << nl;
+    }
+
+    os  << decrIndent << token::END_BLOCK << nl;
+
+    os.check
+    (
+        "Foam::Ostream& Foam::operator<<"
+        "(Foam::Ostream&, const Foam::phaseProperties&)"
+    );
+
+    return os;
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/phaseProperties/phasePropertiesList/phasePropertiesList.C

@@ -22,44 +22,68 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Typedefs
-    Foam::cloudThermoTypes
-
-Description
-
 \*---------------------------------------------------------------------------*/
 
-#ifndef ReactingCloudThermoTypes_H
+#include "phasePropertiesList.H"
-#define ReactingCloudThermoTypes_H
 
-#include "sutherlandTransport.H"
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-#include "multiComponentMixture.H"
-#include "specie.H"
-#include "constTransport.H"
-#include "specieThermo.H"
-#include "hConstThermo.H"
-#include "janafThermo.H"
-#include "perfectGas.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+Foam::phasePropertiesList::phasePropertiesList
-
+(
-namespace Foam
+    Istream& is,
+    const PtrList<volScalarField>& YGas,
+    const wordList& liquidNames,
+    const wordList& solidNames
+)
+:
+    props_(is),
+    phaseTypeNames_()
 {
-//    typedef multiComponentMixture<constTransport<specieThermo<hConstThermo<perfectGas> > > > specieProperties;
+    forAll(props_, i)
-//    typedef hConstThermo<perfectGas> specieProperties;
+    {
+        props_[i].initialiseGlobalIds(YGas, liquidNames, solidNames);
+    }
 
-//    typedef sutherlandTransport<specieThermo<janafThermo<perfectGas> > >
+    phaseTypeNames_.setSize(props_.size());
-//        specieProperties;
+    forAll(phaseTypeNames_, i)
-
+    {
-    typedef sutherlandTransport<specieThermo<janafThermo<perfectGas> > >
+        phaseTypeNames_[i] = props_[i].phaseTypeName();
-        specieReactingProperties;
+    }
-
-    typedef constTransport<specieThermo<hConstThermo<perfectGas> > >
-        specieConstProperties;
 }
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#endif
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::phasePropertiesList::~phasePropertiesList()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::List<Foam::phaseProperties>&
+Foam::phasePropertiesList::props() const
+{
+    return props_;
+}
+
+
+const Foam::wordList& Foam::phasePropertiesList::phaseTypes() const
+{
+    return phaseTypeNames_;
+}
+
+
+Foam::label Foam::phasePropertiesList::size() const
+{
+    return props_.size();
+}
+
+
+const Foam::phaseProperties&
+Foam::phasePropertiesList::operator[](const label phaseI) const
+{
+    return props_[phaseI];
+}
+
 
 // ************************************************************************* //

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

@@ -0,0 +1,104 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::phasePropertiesList
+
+Description
+    Simple container for a list of phase properties
+
+SourceFiles
+    phasePropertiesList.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef phasePropertiesList_H
+#define phasePropertiesList_H
+
+#include "Istream.H"
+#include "volFields.H"
+#include "wordList.H"
+#include "phaseProperties.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class phasePropertiesList Declaration
+\*---------------------------------------------------------------------------*/
+
+class phasePropertiesList
+{
+    // Private data
+
+        //- List of phase properties
+        List<phaseProperties> props_;
+
+        //- List of word representation of phase types
+        wordList phaseTypeNames_;
+
+
+public:
+
+    //- Constructor
+    phasePropertiesList
+    (
+        Istream& is,
+        const PtrList<volScalarField>& YGas,
+        const wordList& liquidNames,
+        const wordList& solidNames
+    );
+
+    //- Destructor
+    ~phasePropertiesList();
+
+
+    // Public member functions
+
+        //- Return the list of phase properties
+        const List<phaseProperties>& props() const;
+
+        //- Return the list of word representation of phase types
+        const wordList& phaseTypes() const;
+
+        //- Return the size (number of phases)
+        label size() const;
+
+
+    // Member operators
+
+        const phaseProperties& operator[](const label) const;
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -46,7 +46,7 @@ Foam::CompositionModel<CloudType>::CompositionModel
         (
             owner.mesh().objectRegistry::lookupObject<dictionary>
             (
-                "thermophysicalProperties"
+                carrierThermo_.name()
             )
         )
     ),
@@ -56,9 +56,16 @@ Foam::CompositionModel<CloudType>::CompositionModel
         (
             owner.mesh().objectRegistry::lookupObject<dictionary>
             (
-                "thermophysicalProperties"
+                carrierThermo_.name()
             )
         )
+    ),
+    phaseProps_
+    (
+        coeffDict_.lookup("phases"),
+        carrierThermo_.composition().Y(),
+        liquids_().components(),
+        solids_().components()
     )
 {}
 
@@ -123,6 +130,240 @@ const Foam::solidMixture& Foam::CompositionModel<CloudType>::solids() const
 }
 
 
+template<class CloudType>
+const Foam::phasePropertiesList&
+Foam::CompositionModel<CloudType>::phaseProps() const
+{
+    return phaseProps_;
+}
+
+
+template<class CloudType>
+Foam::label Foam::CompositionModel<CloudType>::nPhase() const
+{
+    return phaseProps_.size();
+}
+
+
+template<class CloudType>
+const Foam::wordList&
+Foam::CompositionModel<CloudType>::phaseTypes() const
+{
+    // if only 1 phase, return the constituent component names
+    if (phaseProps_.size() == 1)
+    {
+        return phaseProps_[0].names();
+    }
+    else
+    {
+        return phaseProps_.phaseTypes();
+    }
+}
+
+
+template<class CloudType>
+const Foam::wordList&
+Foam::CompositionModel<CloudType>::componentNames(const label phaseI) const
+{
+    return phaseProps_[phaseI].names();
+}
+
+
+template<class CloudType>
+Foam::label Foam::CompositionModel<CloudType>::globalId
+(
+    const label phaseI,
+    const word& cmptName
+) const
+{
+    label id = phaseProps_[phaseI].globalId(cmptName);
+
+    if (id < 0)
+    {
+        FatalErrorIn
+        (
+            "Foam::label Foam::CompositionModel<CloudType>::globalId\n"
+            "(\n"
+            "    const label phaseI,\n"
+            "    const word& cmptName\n"
+            ") const"
+        )   << "Unable to determine global id for requested component "
+            << cmptName << nl << abort(FatalError);
+    }
+
+    return id;
+}
+
+
+template<class CloudType>
+const Foam::labelList& Foam::CompositionModel<CloudType>::globalIds
+(
+    const label phaseI
+) const
+{
+    return phaseProps_[phaseI].globalIds();
+}
+
+
+template<class CloudType>
+Foam::label Foam::CompositionModel<CloudType>::localId
+(
+    const label phaseI,
+    const word& cmptName
+) const
+{
+    label id = phaseProps_[phaseI].id(cmptName);
+
+    if (id < 0)
+    {
+        FatalErrorIn
+        (
+            "Foam::label Foam::CompositionModel<CloudType>::localId\n"
+            "(\n"
+            "    const label phaseI,\n"
+            "    const word& cmptName\n"
+            ") const"
+        )   << "Unable to determine local id for component " << cmptName
+            << nl << abort(FatalError);
+    }
+
+    return id;
+}
+
+
+template<class CloudType>
+const Foam::scalarField& Foam::CompositionModel<CloudType>::Y0
+(
+    const label phaseI
+) const
+{
+    return phaseProps_[phaseI].Y();
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::CompositionModel<CloudType>::H
+(
+    const label phaseI,
+    const scalarField& Y,
+    const scalar p,
+    const scalar T
+) const
+{
+    const phaseProperties& props = phaseProps_[phaseI];
+    scalar HMixture = 0.0;
+    switch (props.phase())
+    {
+        case phaseProperties::GAS:
+        {
+            forAll(Y, i)
+            {
+                label gid = props.globalIds()[i];
+                HMixture += Y[i]*this->gases()[gid].H(T);
+            }
+            break;
+        }
+        case phaseProperties::LIQUID:
+        {
+            forAll(Y, i)
+            {
+                label gid = props.globalIds()[i];
+                HMixture += Y[i]*this->liquids().properties()[gid].h(p, T);
+            }
+            break;
+        }
+        case phaseProperties::SOLID:
+        {
+            forAll(Y, i)
+            {
+                label gid = props.globalIds()[i];
+                HMixture +=
+                     Y[i]
+                    *(
+                        this->solids().properties()[gid].Hf()
+                      + this->solids().properties()[gid].cp()*T
+                     );
+            }
+            break;
+        }
+        default:
+        {
+            FatalErrorIn
+            (
+                "Foam::scalar Foam::CompositionModel<CloudType>::H\n"
+                "(\n"
+                "    const label phaseI,\n"
+                "    const scalarField& Y,"
+                "    const scalar p,\n"
+                "    const scalar T\n"
+                ") const"
+            )   << "Unknown phase enumeration" << nl << abort(FatalError);
+        }
+    }
+
+    return HMixture;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::CompositionModel<CloudType>::cp
+(
+    const label phaseI,
+    const scalarField& Y,
+    const scalar p,
+    const scalar T
+) const
+{
+    const phaseProperties& props = phaseProps_[phaseI];
+    scalar cpMixture = 0.0;
+    switch (props.phase())
+    {
+        case phaseProperties::GAS:
+        {
+            forAll(Y, i)
+            {
+                label gid = props.globalIds()[i];
+                cpMixture += Y[i]*this->gases()[gid].Cp(T);
+            }
+            break;
+        }
+        case phaseProperties::LIQUID:
+        {
+            forAll(Y, i)
+            {
+                label gid = props.globalIds()[i];
+                cpMixture += Y[i]*this->liquids().properties()[gid].cp(p, T);
+            }
+            break;
+        }
+        case phaseProperties::SOLID:
+        {
+            forAll(Y, i)
+            {
+                label gid = props.globalIds()[i];
+                cpMixture += Y[i]*this->solids().properties()[gid].cp();
+            }
+            break;
+        }
+        default:
+        {
+            FatalErrorIn
+            (
+                "Foam::scalar Foam::CompositionModel<CloudType>::cp\n"
+                "(\n"
+                "    const label phaseI,\n"
+                "    const scalarField& Y,"
+                "    const scalar p,\n"
+                "    const scalar T\n"
+                ") const"
+            )   << "Unknown phase enumeration" << nl << abort(FatalError);
+        }
+    }
+
+    return cpMixture;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #include "NewCompositionModel.C"

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

@@ -25,7 +25,6 @@ License
 Class
     Foam::CompositionModel
 
-
 Description
     Templated reacting parcel composition model class
     Consists of gases (via thermo package), liquids and solids
@@ -49,6 +48,8 @@ SourceFiles
 #include "liquidMixture.H"
 #include "solidMixture.H"
 
+#include "phasePropertiesList.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -84,6 +85,9 @@ class CompositionModel
         //- Global solid properties data
         autoPtr<solidMixture> solids_;
 
+        //- List of phase properties
+        phasePropertiesList phaseProps_;
+
 
 public:
 
@@ -143,103 +147,85 @@ public:
             //- Return the carrier phase thermo package
             const hCombustionThermo& carrierThermo() const;
 
-            //- Return the gases
-            const PtrList<specieReactingProperties>& gases() const;
 
-            //- Return the global liquids
+            // Composition lists
-            const liquidMixture& liquids() const;
 
-            //- Return the global solids
+                //- Return the gases
-            const solidMixture& solids() const;
+                const PtrList<specieReactingProperties>& gases() const;
 
-            //- Return the list of composition names
+                //- Return the global liquids
-            virtual const wordList compositionNames() const = 0;
+                const liquidMixture& liquids() const;
 
-            //- Return the list of gas names
+                //- Return the global solids
-            virtual const wordList& gasNames() const = 0;
+                const solidMixture& solids() const;
 
-            //- Return local id of gas based on gasName
+                //- Return the list of phase properties
-            virtual label gasLocalId(const word& gasName) const = 0;
+                const phasePropertiesList& phaseProps() const;
 
-            //- Return global id of gas based on gasName
+                //- Return the number of phases
-            virtual label gasGlobalId(const word& gasName) const = 0;
+                label nPhase() const;
 
-            //- Return the list indices of gases in global thermo list
-            virtual const labelList& gasGlobalIds() const = 0;
 
-            //- Return the list of gas mass fractions
+            // Phase properties
-            virtual const scalarField& YGas0() const = 0;
 
-            //- Return the list of liquid names
+                //- Return the list of phase type names
-            virtual const wordList& liquidNames() const = 0;
+                //  If only 1 phase, return the component names of that phase
+                const wordList& phaseTypes() const;
 
-            //- Return local id of liquid based on liquidName
+                //- Return the list of component names for phaseI
-            virtual label liquidLocalId(const word& liquidName) const = 0;
+                const wordList& componentNames(const label phaseI) const;
 
-            //- Return global id of liquid based on liquidName
+                //- Return global id of component cmptName in phase phaseI
-            virtual label liquidGlobalId(const word& liquidName) const = 0;
+                label globalId(const label phaseI, const word& cmptName) const;
 
-            //- Return the list indices of liquid in global thermo list
+                //- Return global ids of for phase phaseI
-            virtual const labelList& liquidGlobalIds() const = 0;
+                const labelList& globalIds(const label phaseI) const;
 
-            //- Return the list of liquid mass fractions
+                //- Return local id of component cmptName in phase phaseI
-            virtual const scalarField& YLiquid0() const = 0;
+                label localId(const label phaseI, const word& cmptName) const;
 
-            //- Return the list of solid names
+                //- Return the list of phase phaseI mass fractions
-            virtual const wordList& solidNames() const = 0;
+                const scalarField& Y0(const label phaseI) const;
 
-            //- Return local id of solid based on solidName
-            virtual label solidLocalId(const word& solidName) const = 0;
 
-            //- Return global id of solid based on solidName
+            // Mixture properties
-            virtual label solidGlobalId(const word& solidName) const = 0;
 
-            //- Return the list indices of solids in global thermo list
+                //- Return the list of mixture mass fractions
-            virtual const labelList& solidGlobalIds() const = 0;
+                //  If only 1 phase, return component fractions of that phase
+                virtual const scalarField& YMixture0() const = 0;
 
-            //- Return the list of solid mass fractions
+                // Indices of gas, liquid and solid phases in phase properties
-            virtual const scalarField& YSolid0() const = 0;
+                // list - returns -1 if not applicable
 
-            //- Return the list of mixture mass fractions
+                    //- Gas id
-            virtual const scalarField& YMixture0() const = 0;
+                    virtual label idGas() const = 0;
 
-            //- Return the gas constant for the gas mixture
+                    //- Liquid id
-            virtual scalar RGas(const scalarField& YGas) const = 0;
+                    virtual label idLiquid() const = 0;
 
-            //- Return enthalpy for the gas mixture [energy per unit mass]
+                    //- Solid id
-            virtual scalar HGas
+                    virtual label idSolid() const = 0;
+
+
+        // Evaluation
+
+            //- Return enthalpy for the phase phaseI
+            virtual scalar H
             (
-                 const scalarField& YGas,
+                const label phaseI,
-                 const scalar T
+                const scalarField& Y,
-            ) const = 0;
-
-            //- Return enthalpy for the solid mixture [energy per unit mass]
-            virtual scalar HSolid
-            (
-                const scalarField& YSolid,
-                const scalar T
-            ) const = 0;
-
-            //- Return specific heat caparcity for the gas mixture
-            virtual scalar cpGas
-            (
-                 const scalarField& YGas,
-                 const scalar T
-            ) const = 0;
-
-            //- Return specific heat caparcity for the liquid mixture
-            virtual scalar cpLiquid
-            (
-                const scalarField& YLiquid,
                 const scalar p,
                 const scalar T
-            ) const = 0;
+            ) const;
 
-            //- Return specific heat caparcity for the solid mixture
+            //- Return cp for the phase phaseI
-            virtual scalar cpSolid
+            virtual scalar cp
             (
-                const scalarField& YSolid
+                const label phaseI,
-            ) const = 0;
+                const scalarField& Y,
+                const scalar p,
+                const scalar T
+            ) const;
 };
 
 

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

@@ -36,34 +36,32 @@ Foam::CompositionModel<CloudType>::New
     CloudType& owner
 )
 {
-    word ReactingCompositionModelType
+    word CompositionModelType(dict.lookup("CompositionModel"));
-    (
-        dict.lookup("CompositionModel")
-    );
 
-    Info<< "Selecting CompositionModel "
+    Info<< "Selecting CompositionModel " << CompositionModelType << endl;
-        << ReactingCompositionModelType << endl;
 
     typename dictionaryConstructorTable::iterator cstrIter =
-        dictionaryConstructorTablePtr_->find(ReactingCompositionModelType);
+        dictionaryConstructorTablePtr_->find(CompositionModelType);
 
     if (cstrIter == dictionaryConstructorTablePtr_->end())
     {
         FatalErrorIn
         (
-            "CompositionModel<CloudType>::New"
+            "CompositionModel<CloudType>::New\n"
-            "(const dictionary&, CloudType&)"
+            "(\n"
+            "    const dictionary&,\n"
+            "    CloudType&\n"
+            ")"
         )
-            << "Unknown ReactingCompositionModelType type "
+            << "Unknown CompositionModelType type "
-            << ReactingCompositionModelType
+            << CompositionModelType
             << ", constructor not in hash table" << nl << nl
             << "    Valid CompositionModel types are :" << nl
-            << dictionaryConstructorTablePtr_->toc()
+            << dictionaryConstructorTablePtr_->toc() << nl
             << exit(FatalError);
     }
 
-    return autoPtr<CompositionModel<CloudType> >
+    return autoPtr<CompositionModel<CloudType> >(cstrIter()(dict, owner));
-        (cstrIter()(dict, owner));
 }
 
 

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

@@ -26,6 +26,62 @@ License
 
 #include "SingleMixtureFraction.H"
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::SingleMixtureFraction<CloudType>::constructIds()
+{
+    forAll(this->phaseProps(), phaseI)
+    {
+        switch (this->phaseProps()[phaseI].phase())
+        {
+            case phaseProperties::GAS:
+            {
+                idGas_ = phaseI;
+                break;
+            }
+            case phaseProperties::LIQUID:
+            {
+                idLiquid_ = phaseI;
+                break;
+            }
+            case phaseProperties::SOLID:
+            {
+                idSolid_ = phaseI;
+                break;
+            }
+            default:
+            {
+                FatalErrorIn
+                (
+                    "void Foam::SingleMixtureFraction<CloudType>::"
+                    "constructIds()"
+                )   << "Unknown phase enumeration" << nl << abort(FatalError);
+            }
+        }
+    }
+
+    if (idGas_ < 0)
+    {
+        FatalErrorIn("Foam::SingleMixtureFraction<CloudType>::constructIds()")
+            << "No gas phase found in phase list:" << nl
+            << this->phaseTypes() << nl << endl;
+    }
+    if (idLiquid_ < 0)
+    {
+        FatalErrorIn("Foam::SingleMixtureFraction<CloudType>::constructIds()")
+            << "No liquid phase found in phase list:" << nl
+            << this->phaseTypes() << nl << endl;
+    }
+    if (idSolid_ < 0)
+    {
+        FatalErrorIn("Foam::SingleMixtureFraction<CloudType>::constructIds()")
+            << "No solid phase found in phase list:" << nl
+            << this->phaseTypes() << nl << endl;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -37,181 +93,45 @@ Foam::SingleMixtureFraction<CloudType>::SingleMixtureFraction
 :
     CompositionModel<CloudType>(dict, owner, typeName),
 
-    gasNames_(this->coeffDict().lookup("gasNames")),
+    idGas_(-1),
-    gasGlobalIds_(gasNames_.size(), -1),
+    idLiquid_(-1),
-    YGas0_(this->coeffDict().lookup("YGas0")),
+    idSolid_(-1),
-    YGasTot0_(readScalar(this->coeffDict().lookup("YGasTot0"))),
-
-    liquidNames_(this->coeffDict().lookup("liquidNames")),
-    liquidGlobalIds_(liquidNames_.size(), -1),
-    YLiquid0_(this->coeffDict().lookup("YLiquid0")),
-    YLiquidTot0_(readScalar(this->coeffDict().lookup("YLiquidTot0"))),
-
-    solidNames_(this->coeffDict().lookup("solidNames")),
-    solidGlobalIds_(solidNames_.size(), -1),
-    YSolid0_(this->coeffDict().lookup("YSolid0")),
-    YSolidTot0_(readScalar(this->coeffDict().lookup("YSolidTot0"))),
 
     YMixture0_(3)
 {
-    // Construct gasGlobalIds_ list
+    constructIds();
-    forAll(gasNames_, i)
-    {
-        forAll (this->carrierThermo().composition().Y(), j)
-        {
-            word specieName(this->carrierThermo().composition().Y()[j].name());
 
-            if (specieName == gasNames_[i])
+    if (this->phaseProps().size() != 3)
-            {
-                gasGlobalIds_[i] = j;
-                break;
-            }
-        }
-        if (gasGlobalIds_[i] == -1)
-        {
-            Info<< "\nThermo package species composition comprises:" << endl;
-            forAll (this->carrierThermo().composition().Y(), k)
-            {
-                Info<< this->carrierThermo().composition().Y()[k].name()
-                    << endl;
-            }
-
-            FatalErrorIn
-            (
-                "Foam::SingleMixtureFraction<CloudType>::SingleMixtureFraction"
-                "(const dictionary&, CloudType&)"
-            )   << "Could not find gas species " << gasNames_[i]
-                << " in species list" <<  exit(FatalError);
-        }
-    }
-
-    // Construct liquidGlobalIds_ list
-    forAll(liquidNames_, i)
-    {
-        forAll (this->liquids().components(), j)
-        {
-            word specieName(this->liquids().components()[j]);
-
-            if (specieName == liquidNames_[i])
-            {
-                liquidGlobalIds_[i] = j;
-                break;
-            }
-        }
-        if (liquidGlobalIds_[i] == -1)
-        {
-            Info<< "\nLiquid mixture species composition comprises:" << endl;
-            forAll (this->liquids().components(), k)
-            {
-                Info<< this->liquids().components()[k] << endl;
-            }
-
-            FatalErrorIn
-            (
-                "Foam::SingleMixtureFraction<CloudType>::SingleMixtureFraction"
-                "(const dictionary&, CloudType&)"
-            )   << "Could not find liquid species " << liquidNames_[i]
-                << " in species list" <<  exit(FatalError);
-        }
-    }
-
-    // Construct solidGlobalIds_ list
-    forAll(solidNames_, i)
-    {
-        forAll (this->solids().components(), j)
-        {
-            word specieName(this->solids().components()[j]);
-
-            if (specieName == solidNames_[i])
-            {
-                solidGlobalIds_[i] = j;
-                break;
-            }
-        }
-        if (solidGlobalIds_[i] == -1)
-        {
-            Info<< "\nSolid mixture species composition comprises:" << endl;
-            forAll (this->solids().components(), k)
-            {
-                Info<< this->solids().components()[k] << endl;
-            }
-
-            FatalErrorIn
-            (
-                "Foam::SingleMixtureFraction<CloudType>::SingleMixtureFraction"
-                "(const dictionary&, CloudType&)"
-            )   << "Could not find solid species " << solidNames_[i]
-                << " in species list" <<  exit(FatalError);
-        }
-    }
-
-    // Set mixture fractions
-    YMixture0_[0] = YGasTot0_;
-    YMixture0_[1] = YLiquidTot0_;
-    YMixture0_[2] = YSolidTot0_;
-
-    // Check that total mass fractions = 1
-
-    if (YGas0_.size())
-    {
-        if (mag(sum(YGas0_) - 1) > SMALL)
-        {
-            FatalErrorIn
-            (
-                "Foam::SingleMixtureFraction<CloudType>::SingleMixtureFraction"
-                "(const dictionary&, CloudType&)"
-            )<< "Mass fractions of YGas0 must sum to unity"
-             <<  exit(FatalError);
-        }
-    }
-    else
-    {
-        YMixture0_[0] = 0.0;
-    }
-
-    if (YLiquid0_.size())
-    {
-        if (mag(sum(YLiquid0_) - 1) > SMALL)
-        {
-            FatalErrorIn
-            (
-                "Foam::SingleMixtureFraction<CloudType>::SingleMixtureFraction"
-                "(const dictionary&, CloudType&)"
-            )<< "Mass fractions of YLiquid0 must sum to unity"
-             <<  exit(FatalError);
-        }
-    }
-    else
-    {
-        YMixture0_[1] = 0.0;
-    }
-
-    if (YSolid0_.size())
-    {
-        if (mag(sum(YSolid0_) - 1) > SMALL)
-        {
-            FatalErrorIn
-            (
-                "Foam::SingleMixtureFraction<CloudType>::SingleMixtureFraction"
-                "(const dictionary&, CloudType&)"
-            )<< "Mass fractions of YSolid0 must sum to unity"
-             <<  exit(FatalError);
-        }
-    }
-    else
-    {
-        YMixture0_[2] = 0.0;
-    }
-
-    // Check total mixture fraction sums to 1
-    if (mag(sum(YMixture0_) - 1) > SMALL)
     {
         FatalErrorIn
         (
-            "Foam::SingleMixtureFraction<CloudType>::SingleMixtureFraction"
+            "Foam::SingleMixtureFraction<CloudType>::"
-            "(const dictionary&, CloudType&)"
+            "SingleMixtureFraction\n"
-        )   << "Mass fractions YGasTot0 + YSolidTot0 + YSolidTot0 must sum "
+            "(\n"
-            << "to unity" <<  exit(FatalError);
+            "    const dictionary&,\n"
+            "    CloudType&\n"
+            ")\n"
+        )   << "Incorrect numebr of phases: " << nl
+            << "    Please specify 1 gas, 1 liquid and 1 solid" << nl
+            << exit(FatalError);
+    }
+
+    this->coeffDict().lookup("YGasTot0") >> YMixture0_[idGas_];
+    this->coeffDict().lookup("YLiquidTot0") >> YMixture0_[idLiquid_];
+    this->coeffDict().lookup("YSolidTot0") >> YMixture0_[idSolid_];
+
+    if (mag(sum(YMixture0_) - 1.0) > SMALL)
+    {
+        FatalErrorIn
+        (
+            "Foam::SingleMixtureFraction<CloudType>::"
+            "SingleMixtureFraction\n"
+            "(\n"
+            "    const dictionary&,\n"
+            "    CloudType&\n"
+            ")\n"
+        )   << "Sum of phases should be 1. Phase fractions:" << nl
+            << YMixture0_ << nl << exit(FatalError);
     }
 }
 
@@ -225,157 +145,11 @@ Foam::SingleMixtureFraction<CloudType>::~SingleMixtureFraction()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class CloudType>
-const Foam::wordList
-Foam::SingleMixtureFraction<CloudType>::compositionNames() const
-{
-    wordList names(3);
-    names[0] = "Gas";
-    names[1] = "Liquid";
-    names[2] = "Solid";
-    return names;
-}
-
-template<class CloudType>
-const Foam::wordList&
-Foam::SingleMixtureFraction<CloudType>::gasNames() const
-{
-     return gasNames_;
-}
-
-
-template<class CloudType>
-Foam::label
-Foam::SingleMixtureFraction<CloudType>::gasLocalId(const word& gasName) const
-{
-    forAll(gasNames_, i)
-    {
-        if (gasName == gasNames_[i])
-        {
-            return i;
-        }
-    }
-
-    WarningIn
-    (
-        "Foam::label SingleMixtureFraction<CloudType>::"
-        "gasLocalId(const word& gasName) const"
-    )<< "Gas name " << gasName << " not found in gasNames_"
-     << endl;
-
-    return -1;
-}
-
-
-template<class CloudType>
-Foam::label
-Foam::SingleMixtureFraction<CloudType>::gasGlobalId(const word& gasName) const
-{
-    forAll(gasNames_, i)
-    {
-        if (gasName == gasNames_[i])
-        {
-            return gasGlobalIds_[i];
-        }
-    }
-
-    WarningIn
-    (
-        "Foam::label SingleMixtureFraction<CloudType>::"
-        "gasGlobalId(const word& gasName) const"
-    )<< "Gas name " << gasName << " not found in gasNames_"
-     << endl;
-
-    return -1;
-}
-
-
-template<class CloudType>
-const Foam::labelList&
-Foam::SingleMixtureFraction<CloudType>::gasGlobalIds() const
-{
-    return gasGlobalIds_;
-}
-
-
 template<class CloudType>
 const Foam::scalarField&
 Foam::SingleMixtureFraction<CloudType>::YGas0() const
 {
-    return YGas0_;
+    return this->phaseProps()[idGas_].Y();
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::SingleMixtureFraction<CloudType>::YGasTot0() const
-{
-    return YGasTot0_;
-}
-
-
-template<class CloudType>
-const Foam::wordList&
-Foam::SingleMixtureFraction<CloudType>::liquidNames() const
-{
-     return liquidNames_;
-}
-
-
-template<class CloudType>
-Foam::label Foam::SingleMixtureFraction<CloudType>::liquidLocalId
-(
-    const word& liquidName
-) const
-{
-    forAll(liquidNames_, i)
-    {
-        if (liquidName == liquidNames_[i])
-        {
-            return i;
-        }
-    }
-
-    WarningIn
-    (
-        "Foam::label SingleMixtureFraction<CloudType>::"
-        "liquidLocalId(const word& liquidName) const"
-    )<< "Liquid name " << liquidName << " not found in liquidNames_"
-     << endl;
-
-    return -1;
-}
-
-
-template<class CloudType>
-Foam::label Foam::SingleMixtureFraction<CloudType>::liquidGlobalId
-(
-    const word& liquidName
-) const
-{
-    forAll(liquidNames_, i)
-    {
-        if (liquidName == liquidNames_[i])
-        {
-            return liquidGlobalIds_[i];
-        }
-    }
-
-    WarningIn
-    (
-        "Foam::label SingleMixtureFraction<CloudType>::"
-        "liquidGlobalId(const word& liquidName) const"
-    )<< "Liquid name " << liquidName << " not found in liquidNames_"
-     << endl;
-
-    return -1;
-}
-
-
-template<class CloudType>
-const Foam::labelList&
-Foam::SingleMixtureFraction<CloudType>::liquidGlobalIds() const
-{
-    return liquidGlobalIds_;
 }
 
 
@@ -383,81 +157,7 @@ template<class CloudType>
 const Foam::scalarField&
 Foam::SingleMixtureFraction<CloudType>::YLiquid0() const
 {
-    return YLiquid0_;
+    return this->phaseProps()[idLiquid_].Y();
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::SingleMixtureFraction<CloudType>::YLiquidTot0() const
-{
-    return YLiquidTot0_;
-}
-
-
-template<class CloudType>
-const Foam::wordList&
-Foam::SingleMixtureFraction<CloudType>::solidNames() const
-{
-    return solidNames_;
-}
-
-
-template<class CloudType>
-Foam::label Foam::SingleMixtureFraction<CloudType>::solidLocalId
-(
-    const word& solidName
-) const
-{
-    forAll(solidNames_, i)
-    {
-        if (solidName == solidNames_[i])
-        {
-            return i;
-        }
-    }
-
-    WarningIn
-    (
-        "Foam::label SingleMixtureFraction<CloudType>::"
-        "SolididLocalId(const word& solidName) const"
-    )<< "Solid name " << solidName << " not found in solidNames_"
-     << endl;
-
-    return -1;
-}
-
-
-template<class CloudType>
-Foam::label
-Foam::SingleMixtureFraction<CloudType>::solidGlobalId
-(
-    const word& solidName
-) const
-{
-    forAll(solidNames_, i)
-    {
-        if (solidName == solidNames_[i])
-        {
-            return solidGlobalIds_[i];
-        }
-    }
-
-    WarningIn
-    (
-        "Foam::label SingleMixtureFraction<CloudType>::"
-        "solidGlobalId(const word& solidName) const"
-    )<< "Solid name " << solidName << " not found in solidNames_"
-     << endl;
-
-    return -1;
-}
-
-
-template<class CloudType>
-const Foam::labelList&
-Foam::SingleMixtureFraction<CloudType>::solidGlobalIds() const
-{
-    return solidGlobalIds_;
 }
 
 
@@ -465,15 +165,7 @@ template<class CloudType>
 const Foam::scalarField&
 Foam::SingleMixtureFraction<CloudType>::YSolid0() const
 {
-    return YSolid0_;
+    return this->phaseProps()[idSolid_].Y();
-}
-
-
-template<class CloudType>
-Foam::scalar
-Foam::SingleMixtureFraction<CloudType>::YSolidTot0() const
-{
-    return YSolidTot0_;
 }
 
 
@@ -486,115 +178,23 @@ Foam::SingleMixtureFraction<CloudType>::YMixture0() const
 
 
 template<class CloudType>
-Foam::scalar Foam::SingleMixtureFraction<CloudType>::RGas
+Foam::label Foam::SingleMixtureFraction<CloudType>::idGas() const
-(
-    const scalarField& YGas
-) const
 {
-    scalar RGasMixture = 0.0;
+    return idGas_;
-    forAll(YGas, i)
-    {
-        label id = gasGlobalIds_[i];
-        RGasMixture += YGas[i]*this->gases()[id].R();
-    }
-    return RGasMixture;
 }
 
 
 template<class CloudType>
-Foam::scalar Foam::SingleMixtureFraction<CloudType>::HGas
+Foam::label Foam::SingleMixtureFraction<CloudType>::idLiquid() const
-(
-    const scalarField& YGas,
-    const scalar T
-)
-const
 {
-    scalar HMixture = 0.0;
+    return idLiquid_;
-    forAll(YGas, i)
-    {
-        label id = gasGlobalIds_[i];
-        HMixture += YGas[i]*this->gases()[id].H(T);
-    }
-    return HMixture;
 }
 
 
 template<class CloudType>
-Foam::scalar Foam::SingleMixtureFraction<CloudType>::HSolid
+Foam::label Foam::SingleMixtureFraction<CloudType>::idSolid() const
-(
-    const scalarField& YSolid,
-    const scalar T
-)
-const
 {
-    scalar HMixture = 0.0;
+    return idSolid_;
-    forAll(YSolid, i)
-    {
-        label id = solidGlobalIds_[i];
-        HMixture +=
-            YSolid[i]
-           *(
-                this->solids().properties()[id].Hf()
-              + this->solids().properties()[id].cp()*T
-            );
-    }
-    return HMixture;
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::SingleMixtureFraction<CloudType>::cpGas
-(
-    const scalarField& YGas,
-    const scalar T
-)
-const
-{
-    scalar cpMixture = 0.0;
-    forAll(YGas, i)
-    {
-        label id = gasGlobalIds_[i];
-        cpMixture += YGas[i]*this->gases()[id].Cp(T);
-    }
-    return cpMixture;
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::SingleMixtureFraction<CloudType>::cpLiquid
-(
-    const scalarField& YLiquid,
-    const scalar p,
-    const scalar T
-)
-const
-{
-    scalar cpMixture = 0.0;
-    forAll(YLiquid, i)
-    {
-        label id = liquidGlobalIds_[i];
-        cpMixture += YLiquid[i]*this->liquids().properties()[id].cp(p, T);
-    }
-
-    return cpMixture;
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::SingleMixtureFraction<CloudType>::cpSolid
-(
-    const scalarField& YSolid
-)
-const
-{
-    scalar cpMixture = 0.0;
-    forAll(YSolid, i)
-    {
-        label id = solidGlobalIds_[i];
-        cpMixture += YSolid[i]*this->solids().properties()[id].cp();
-    }
-
-    return cpMixture;
 }
 
 

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

@@ -25,10 +25,8 @@ License
 Class
     Foam::SingleMixtureFraction
 
-
 Description
-    Templated parcel single mixture fraction class
+    Templated parcel multi-phase, multi-component class
-    - Each phase sums to a mass fraction of 1
 
 SourceFiles
     SingleMixtureFraction.C
@@ -56,67 +54,31 @@ class SingleMixtureFraction
 {
     // Private data
 
-        //- Parcel properties
+        // Indices of the phases
 
-            //- List of gas names
+            //- Gas
-            const wordList gasNames_;
+            label idGas_;
 
-            //- Return local id of gas based on gasName
+            //- Liquid
-            label gasLocalId(const word& gasName) const;
+            label idLiquid_;
 
-            //- Return global id of gas based on gasName
+            //- Solid
-            label gasGlobalId(const word& gasName) const;
+            label idSolid_;
 
-            //- List of IDs of gases in global list
-            labelList gasGlobalIds_;
 
-            //- Initial mass fractions of gases
+       // Mixture properties
-            const scalarField YGas0_;
 
-            //- Initial total mass fractions of gases
+            //- Phase component total fractions
-            const scalar YGasTot0_;
-
-            //- List of liquid names
-            const wordList liquidNames_;
-
-            //- Return local id of liquid based on liquidName
-            label liquidLocalId(const word& liquidName) const;
-
-            //- Return global id of liquid based on liquidName
-            label liquidGlobalId(const word& liquidName) const;
-
-            //- List of IDs of liquids in global list
-            labelList liquidGlobalIds_;
-
-            //- Initial mass fractions of liquids
-            const scalarField YLiquid0_;
-
-            //- Initial total mass fractions of liquids
-            const scalar YLiquidTot0_;
-
-            //- List of solids species
-            const wordList solidNames_;
-
-            //- Return local id of solid based on solidName
-            label solidLocalId(const word& solidName) const;
-
-            //- Return global id of solid based on solidName
-            label solidGlobalId(const word& solidName) const;
-
-            //- List of IDs of solids in global list
-            labelList solidGlobalIds_;
-
-            //- Initial mass fractions of solids
-            const scalarField YSolid0_;
-
-            //- Initial total mass fractions of solids
-            const scalar YSolidTot0_;
-
-            //- Initial total mass fractions of the mixture
-            //  Index 0 = gas, 1 = liquid, 2 = solid
             scalarField YMixture0_;
 
 
+    // Private member functions
+
+        //- Construct the indices and check correct specification of
+        //  1 gas, 1 liquid and 1 solid
+        void constructIds();
+
+
 public:
 
     //- Runtime type information
@@ -126,11 +88,7 @@ public:
     // Constructors
 
         //- Construct from dictionary
-        SingleMixtureFraction
+        SingleMixtureFraction(const dictionary& dict, CloudType& owner);
-        (
-            const dictionary& dict,
-            CloudType& owner
-        );
 
 
     //- Destructor
@@ -141,70 +99,49 @@ public:
 
         // Access
 
-            //- Return the list of composition names
+            // Gas properties
-            const wordList compositionNames() const;
 
-            //- Return the list of gas names
+                //- Return the list of gas mass fractions
-            const wordList& gasNames() const;
+                const scalarField& YGas0() const;
 
-            //- Return the list indices of gases in global thermo list
+                //- Return the total gas mass fraction
-            const labelList& gasGlobalIds() const;
+                scalar YGasTot0() const;
 
-            //- Return the list of gas mass fractions
-            const scalarField& YGas0() const;
 
-            //- Return the total gas mass fraction
+            // Liquid properties
-            scalar YGasTot0() const;
 
-            //- Return the list of liquid names
+                //- Return the list of liquid mass fractions
-            const wordList& liquidNames() const;
+                const scalarField& YLiquid0() const;
 
-            //- Return the list indices of liquid in global thermo list
+                //- Return the total liquid mass fraction
-            const labelList& liquidGlobalIds() const;
+                scalar YLiquidTot0() const;
 
-            //- Return the list of liquid mass fractions
-            const scalarField& YLiquid0() const;
 
-            //- Return the total liquid mass fraction
+            // Solid properties
-            scalar YLiquidTot0() const;
 
-             //- Return the list of solid names
+                //- Return the list of solid mass fractions
-            const wordList& solidNames() const;
+                const scalarField& YSolid0() const;
 
-            //- Return the list indices of solids in global thermo list
+                //- Return the total solid mass fraction
-            const labelList& solidGlobalIds() const;
+                scalar YSolidTot0() const;
 
-            //- Return the list of solid mass fractions
-            const scalarField& YSolid0() const;
 
-            //- Return the total solid mass fraction
+            // Mixture properties
-            scalar YSolidTot0() const;
 
-            //- Return the list of mixture mass fractions
+                //- Return the list of mixture mass fractions
-            const scalarField& YMixture0() const;
+                virtual const scalarField& YMixture0() const;
 
-            //- Return the gas constant for the gas mixture
+                // Indices of gas, liquid and solid phases in phase properties
-            scalar RGas(const scalarField& YGas) const;
+                // list
 
-            //- Return enthalpy for the gas mixture [energy per unit mass]
+                    //- Gas id
-            scalar HGas(const scalarField& YGas, const scalar T) const;
+                    virtual label idGas() const;
 
-            //- Return enthalpy for the solid mixture [energy per unit mass]
+                    //- Liquid id
-            scalar HSolid(const scalarField& YSolid, const scalar T) const;
+                    virtual label idLiquid() const;
 
-            //- Return specific heat caparcity for the gas mixture
+                    //- Solid id
-            scalar cpGas(const scalarField& YGas, const scalar T) const;
+                    virtual label idSolid() const;
-
-            //- Return specific heat caparcity for the liquid mixture
-            scalar cpLiquid
-            (
-                const scalarField& YLiquid,
-                const scalar p,
-                const scalar T
-            ) const;
-
-            //- Return specific heat caparcity for the solid mixture
-            scalar cpSolid(const scalarField& YSolid) const;
 };
 
 

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

@@ -0,0 +1,166 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "SinglePhaseMixture.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::SinglePhaseMixture<CloudType>::constructIds()
+{
+    if (this->phaseProps().size() != 1)
+    {
+        FatalErrorIn
+        (
+            "void Foam::SinglePhaseMixture<CloudType>::constructIds()"
+        )   << "Only one phase permitted" << nl << exit(FatalError);
+    }
+
+    switch (this->phaseProps()[0].phase())
+    {
+        case phaseProperties::GAS:
+        {
+            idGas_ = 0;
+            break;
+        }
+        case phaseProperties::LIQUID:
+        {
+            idLiquid_ = 0;
+            break;
+        }
+        case phaseProperties::SOLID:
+        {
+            idSolid_ = 0;
+            break;
+        }
+        default:
+        {
+            FatalErrorIn
+            (
+                "void Foam::SinglePhaseMixture<CloudType>::constructIds()"
+            )   << "Unknown phase enumeration" << nl << abort(FatalError);
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::SinglePhaseMixture<CloudType>::SinglePhaseMixture
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    CompositionModel<CloudType>(dict, owner, typeName),
+    idGas_(-1),
+    idLiquid_(-1),
+    idSolid_(-1)
+{
+    constructIds();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::SinglePhaseMixture<CloudType>::~SinglePhaseMixture()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+const Foam::scalarField&
+Foam::SinglePhaseMixture<CloudType>::YGas0() const
+{
+    notImplemented
+    (
+        "const Foam::scalarField& "
+        "Foam::SinglePhaseMixture<CloudType>::YGas0() const"
+    );
+    return this->phaseProps()[0].Y();
+}
+
+
+template<class CloudType>
+const Foam::scalarField&
+Foam::SinglePhaseMixture<CloudType>::YLiquid0() const
+{
+    notImplemented
+    (
+        "const Foam::scalarField& "
+        "Foam::SinglePhaseMixture<CloudType>::YLiquid0() const"
+    );
+    return this->phaseProps()[0].Y();
+}
+
+
+template<class CloudType>
+const Foam::scalarField&
+Foam::SinglePhaseMixture<CloudType>::YSolid0() const
+{
+    notImplemented
+    (
+        "const Foam::scalarField& "
+        "Foam::SinglePhaseMixture<CloudType>::YSolid0() const"
+    );
+    return this->phaseProps()[0].Y();
+}
+
+
+template<class CloudType>
+const Foam::scalarField&
+Foam::SinglePhaseMixture<CloudType>::YMixture0() const
+{
+    return this->phaseProps()[0].Y();
+}
+
+
+template<class CloudType>
+Foam::label Foam::SinglePhaseMixture<CloudType>::idGas() const
+{
+    return idGas_;
+}
+
+
+template<class CloudType>
+Foam::label Foam::SinglePhaseMixture<CloudType>::idLiquid() const
+{
+    return idLiquid_;
+}
+
+
+template<class CloudType>
+Foam::label Foam::SinglePhaseMixture<CloudType>::idSolid() const
+{
+    return idSolid_;
+}
+
+
+// ************************************************************************* //
+

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

@@ -0,0 +1,156 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::SinglePhaseMixture
+
+Description
+    Templated parcel single phase, multi-component class
+
+SourceFiles
+    SinglePhaseMixture.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef SinglePhaseMixture_H
+#define SinglePhaseMixture_H
+
+#include "CompositionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class SinglePhaseMixture Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class SinglePhaseMixture
+:
+    public CompositionModel<CloudType>
+{
+    // Private data
+
+        // Indices of the phases - only 1 will be set
+
+            //- Gas
+            label idGas_;
+
+            //- Liquid
+            label idLiquid_;
+
+            //- Solid
+            label idSolid_;
+
+
+    // Private member functions
+
+        //- Construct the indices and check correct specification of
+        //  1 gas or 1 liquid or 1 solid
+        void constructIds();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("SinglePhaseMixture");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        SinglePhaseMixture(const dictionary& dict, CloudType& owner);
+
+
+    //- Destructor
+    virtual ~SinglePhaseMixture();
+
+
+    // Member Functions
+
+        // Access
+
+            // Gas properties
+
+                //- Return the list of gas mass fractions
+                const scalarField& YGas0() const;
+
+                //- Return the total gas mass fraction
+                scalar YGasTot0() const;
+
+
+            // Liquid properties
+
+                //- Return the list of liquid mass fractions
+                const scalarField& YLiquid0() const;
+
+                //- Return the total liquid mass fraction
+                scalar YLiquidTot0() const;
+
+
+            // Solid properties
+
+                //- Return the list of solid mass fractions
+                const scalarField& YSolid0() const;
+
+                //- Return the total solid mass fraction
+                scalar YSolidTot0() const;
+
+
+            // Mixture properties
+
+                //- Return the list of mixture mass fractions
+                virtual const scalarField& YMixture0() const;
+
+                // Indices of gas, liquid and solid phases in phase properties
+                // list
+
+                    //- Gas id
+                    virtual label idGas() const;
+
+                    //- Liquid id
+                    virtual label idLiquid() const;
+
+                    //- Solid id
+                    virtual label idSolid() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "SinglePhaseMixture.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -58,8 +58,14 @@ bool Foam::NoPhaseChange<CloudType>::active() const
 template<class CloudType>
 Foam::scalar Foam::NoPhaseChange<CloudType>::calculate
 (
+    const scalar,
     const scalar,
     const scalarField&,
+    const scalarField&,
+    const vector,
+    const scalar,
+    const scalar,
+    const scalar,
     const scalar
 ) const
 {

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

@@ -26,7 +26,7 @@ Class
     Foam::NoPhaseChange
 
 Description
-    Dummy devolatilisation model for 'no phase change'
+    Dummy phase change model for 'no phase change'
 
 \*---------------------------------------------------------------------------*/
 
@@ -76,9 +76,15 @@ public:
         //- Update model
         scalar calculate
         (
-            const scalar,
+            const scalar T,
-            const scalarField&,
+            const scalar d,
-            const scalar
+            const scalarField& Xc,
+            const scalarField& dMassMT,
+            const vector Ur,
+            const scalar Tc,
+            const scalar pc,
+            const scalar nuc,
+            const scalar dt
         ) const;
 };
 

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

@@ -67,6 +67,12 @@ protected:
         const dictionary coeffDict_;
 
 
+    // Protected member functions
+
+        //- Sherwood number
+        scalar Sh() const;
+
+
 public:
 
     //- Runtime type information
@@ -129,9 +135,15 @@ public:
         //- Update model
         virtual scalar calculate
         (
-            const scalar dt,
+            const scalar T,
-            const scalarField& YMixture,
+            const scalar d,
-            const scalar T
+            const scalarField& Xc,
+            const scalarField& dMassMT,
+            const vector Ur,
+            const scalar Tc,
+            const scalar pc,
+            const scalar nuc,
+            const scalar dt
         ) const = 0;
 };
 

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

@@ -0,0 +1,201 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "liquidEvaporation.H"
+#include "specie.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::scalar Foam::liquidEvaporation<CloudType>::Sh
+(
+    const scalar Re,
+    const scalar Sc
+) const
+{
+    return 2.0 + 0.6*Foam::sqrt(Re)*pow(Sc, 0.333);
+}
+
+
+template <class CloudType>
+Foam::scalar Foam::liquidEvaporation<CloudType>::pSat
+(
+    const label i,
+    const scalar T
+) const
+{
+    const List<Tuple2<scalar, scalar> >& pSat = pSat_[i];
+
+    label id = 0;
+    label nT = pSat.size();
+    while ((id < nT) && (pSat[id].first() < T))
+    {
+        id++;
+    }
+
+    if (id == 0 || id == nT - 1)
+    {
+        return pSat[id].second();
+    }
+    else
+    {
+        return
+            (pSat[id].first() - T)
+           /(pSat[id].first() - pSat[id-1].first())
+           *(pSat[id].second() - pSat[id-1].second())
+          + pSat[id-1].second();
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::liquidEvaporation<CloudType>::liquidEvaporation
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    PhaseChangeModel<CloudType>(dict, owner, typeName),
+    gases_(owner.gases()),
+    liquids_
+    (
+        liquidMixture::New
+        (
+            owner.mesh().objectRegistry::lookupObject<dictionary>
+            (
+                owner.carrierThermo().name().name()
+            )
+        )
+    ),
+    Tvap_(readScalar(coeffDict().lookup("Tvap"))),
+    Dm_(coeffDict().lookup("DiffusionCoeffs")),
+    pSat_(coeffDict().lookup("pSat"))
+{
+    if (liquids_.size() != Dm_.size())
+    {
+        FatalErrorIn
+        (
+            "Foam::liquidEvaporation<CloudType>::liquidEvaporation\n"
+            "(\n"
+            "    const dictionary& dict,\n"
+            "    CloudType& cloud\n"
+            ")\n"
+        )   << "Size of diffusionCoeffs list not equal to the number of liquid "
+            << "species" << nl << exit(FatalError);
+    }
+
+    if (liquids_.size() != pSat_.size())
+    {
+        FatalErrorIn
+        (
+            "Foam::liquidEvaporation<CloudType>::liquidEvaporation\n"
+            "(\n"
+            "    const dictionary& dict,\n"
+            "    CloudType& cloud\n"
+            ")\n"
+        )   << "Size of saturation pressure lists not equal to the number of "
+            << "liquid species" << nl << exit(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::liquidEvaporation<CloudType>::~liquidEvaporation()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+bool Foam::liquidEvaporation<CloudType>::active() const
+{
+    return true;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::liquidEvaporation<CloudType>::calculate
+(
+    const scalar T,
+    const scalar d,
+    const scalarField& Xc,
+    const scalarField& dMassMT,
+    const vector Ur
+    const scalar Tc,
+    const scalar pc,
+    const scalar nuc
+    const scalar dt,
+) const
+{
+    if (T < Tvap_)
+    {
+        // not reached temperature threshold
+        return 0.0;
+    }
+    else
+    {
+        // droplet area
+        scalar A = mathematicalConstant::pi*sqr(d);
+
+        // universal gas constant
+        const scalar R = specie::RR.value();
+
+        // Reynolds number
+        scalar Re = mag(Ur)*d/(nuc + ROOTVSMALL);
+
+        // calculate mass transfer of each specie
+        forAll(dMassMT, i)
+        {
+            // Schmidt number
+            scalar Sc = nuc/(Dm_[i] + ROOTVSMALL);
+
+            // Sherwood number
+            scalar Sh = this->Sh(Re, Sc);
+
+            // mass transfer coefficient [m/s]
+            scalar kc = Sh*Dm_[i]/(d + ROOTVSMALL);
+
+            // vapour concentration at droplet surface [kgmol/m3]
+            scalar Cs = pSat(i, T)/(R*T);
+
+            // vapour concentration in bulk gas [kgmol/m3]
+            scalar Cinf = Xc[i]*pc/(R*Tc);
+
+            // molar flux of vapour [kgmol/m2/s]
+            scalar Ni = max(kc*(Cs - Cinf), 0.0);
+
+            // mass transfer
+            dMassMT[i] -= Ni*A*liquids_.properies()[i].W()*dt;
+        }
+    }
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,138 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::liquidEvaporation
+
+Description
+    Liquid evaporation model
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef liquidEvaporation_H
+#define liquidEvaporation_H
+
+#include "PhaseChangeModel.H"
+#include "liquidMixture.H"
+#include "Tuple2.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                       Class liquidEvaporation Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class liquidEvaporation
+:
+    public PhaseChangeModel<CloudType>
+{
+protected:
+
+    // Protected data
+
+        //- Global liquid properties data
+        autoPtr<liquidMixture> liquids_;
+
+        //- Vaporisation temperature [K]
+        //  Local droplet temperature must exceed Tvap before evaporation is
+        //  allowed
+        scalar Tvap_;
+
+        //- Component diffusion coefficients - one per specie
+        //  Note: need to be in same order as defined by liquids in thermo
+        //        props
+        scalarList Dm_;
+
+        //- Component saturation pressure vs temperature
+        List<List<Tuple2<scalar, scalar> > > pSat_;
+
+
+    // Protected member functions
+
+        //- Sherwood number as a function of Reynolds and Schmidt numbers
+        scalar Sh(const scalar Re, const scalar Sc) const;
+
+        //- Return the saturation pressure for species i at temperature, T
+        scalar pSat(const label i, const scalar T) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("liquidEvaporation");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        liquidEvaporation
+        (
+            const dictionary& dict,
+            CloudType& cloud
+        );
+
+
+    //- Destructor
+    virtual ~liquidEvaporation();
+
+
+    // Member Functions
+
+        //- Flag to indicate whether model activates phase change model
+        bool active() const;
+
+        //- Update model
+        scalar calculate
+        (
+            const scalar T,
+            const scalar d,
+            const scalarField& Xc,
+            const scalarField& dMassMT,
+            const vector Ur,
+            const scalar Tc,
+            const scalar pc,
+            const scalar nuc,
+            const scalar dt
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "liquidEvaporation.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -83,7 +83,7 @@ Foam::scalar Foam::HeatTransferModel<CloudType>::h
     const scalar muc
 ) const
 {
-    const scalar Re = rhoc*mag(Ur)*dp/(muc + SMALL);
+    const scalar Re = rhoc*mag(Ur)*dp/(muc + ROOTVSMALL);
 
 //    const scalar Pr = muc/alphac;
     const scalar Pr = this->Pr();