diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
index f35a0cd11f..1072ae273a 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
@@ -562,7 +562,11 @@ void Foam::KinematicCloud<CloudType>::setParcelThermoProperties
     const scalar lagrangianDt
 )
 {
-    parcel.rho() = constProps_.rho0();
+    // If rho0 is given in the const properties
+    if (constProps_.rho0() != -1)
+    {
+        parcel.rho() = constProps_.rho0();
+    }
 }
 
 
@@ -581,6 +585,14 @@ void Foam::KinematicCloud<CloudType>::checkParcelProperties
     {
         parcel.typeId() = constProps_.parcelTypeId();
     }
+
+    if (parcel.rho() == -1)
+    {
+        FatalErrorInFunction
+            << "The kinematic cloud needs rho0 in the constantProperties "
+            << " dictionary. " << nl
+            << abort(FatalError);
+    }
 }
 
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingHeterogeneousCloud/ReactingHeterogeneousCloud.C b/src/lagrangian/intermediate/clouds/Templates/ReactingHeterogeneousCloud/ReactingHeterogeneousCloud.C
index d8800ea264..e11627b45a 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingHeterogeneousCloud/ReactingHeterogeneousCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingHeterogeneousCloud/ReactingHeterogeneousCloud.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2019 OpenCFD Ltd.
+    Copyright (C) 2018-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -136,6 +136,25 @@ void Foam::ReactingHeterogeneousCloud<CloudType>::setParcelThermoProperties
 
     // Set the parcel to combust
     parcel.canCombust() = 1;
+
+    // If rho0 was given in constProp use it. If not use the composition
+    // to set tho
+    if (this->constProps_.rho0() == -1)
+    {
+        const label idGas = this->composition().idGas();
+        const label idLiquid = this->composition().idLiquid();
+        const label idSolid = this->composition().idSolid();
+
+        const scalarField& Ygas = this->composition().Y0(idGas);
+        const scalarField& Yliq = this->composition().Y0(idLiquid);
+        const scalarField& Ysol = this->composition().Y0(idSolid);
+
+        const scalar p0 =
+            this->composition().thermo().thermo().p()[parcel.cell()];
+        const scalar T0 = this->constProps_.T0();
+
+        parcel.rho() = this->composition().rho(Ygas, Yliq, Ysol, T0, p0);
+    }
 }
 
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C b/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C
index 0035a7c29a..b2b68dae3f 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C
@@ -174,6 +174,21 @@ void Foam::ReactingMultiphaseCloud<CloudType>::setParcelThermoProperties
     parcel.YGas() = this->composition().Y0(idGas);
     parcel.YLiquid() = this->composition().Y0(idLiquid);
     parcel.YSolid() = this->composition().Y0(idSolid);
+
+
+    // If rho0 was given in constProp use it. If not use the composition
+    // to set tho
+    if (constProps_.rho0() == -1)
+    {
+        const scalarField& Ygas = this->composition().Y0(idGas);
+        const scalarField& Yliq = this->composition().Y0(idLiquid);
+        const scalarField& Ysol = this->composition().Y0(idSolid);
+        const scalar p0 =
+            this->composition().thermo().thermo().p()[parcel.cell()];
+        const scalar T0 = constProps_.T0();
+
+        parcel.rho() = this->composition().rho(Ygas, Yliq, Ysol, T0, p0);
+    }
 }
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
index 439c509d1f..b2e3758c3d 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
@@ -67,7 +67,7 @@ inline Foam::KinematicParcel<ParcelType>::constantProperties::constantProperties
     dict_(parentDict.subOrEmptyDict("constantProperties")),
     parcelTypeId_(dict_, "parcelTypeId", -1),
     rhoMin_(dict_, "rhoMin", 1e-15),
-    rho0_(dict_, "rho0"),
+    rho0_(dict_, "rho0", -1),
     minParcelMass_(dict_, "minParcelMass", 1e-15)
 {}
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingHeterogeneousParcel/ReactingHeterogeneousParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingHeterogeneousParcel/ReactingHeterogeneousParcel.C
index acd0723a1b..20886adfb5 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingHeterogeneousParcel/ReactingHeterogeneousParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingHeterogeneousParcel/ReactingHeterogeneousParcel.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2019 OpenCFD Ltd.
+    Copyright (C) 2018-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -80,6 +80,30 @@ Foam::scalar Foam::ReactingHeterogeneousParcel<ParcelType>::LEff
 
 // * * * * * * * * * * *  Protected Member Functions * * * * * * * * * * * * //
 
+
+template<class ParcelType>
+template<class TrackCloudType>
+Foam::scalar Foam::ReactingHeterogeneousParcel<ParcelType>::updatedDeltaVolume
+(
+    TrackCloudType& cloud,
+    const scalarField& dMass,
+    const scalar p,
+    const scalar T
+)
+{
+    const auto& composition = cloud.composition();
+
+    scalarField dVolSolid(dMass.size(), Zero);
+    forAll(dVolSolid, i)
+    {
+        dVolSolid[i] =
+            dMass[i]/composition.solids().properties()[i].rho();
+    }
+
+    return sum(dVolSolid);
+}
+
+
 template<class ParcelType>
 template<class TrackCloudType>
 void Foam::ReactingHeterogeneousParcel<ParcelType>::calc
@@ -191,16 +215,52 @@ void Foam::ReactingHeterogeneousParcel<ParcelType>::calc
     (void)this->updateMassFraction(mass0, dMassSolid, this->Y_);
 
 
-    // Update particle density or diameter
-    if (cloud.constProps().constantVolume())
+    if
+    (
+        cloud.constProps().volUpdateType()
+     == constantProperties::volumeUpadteType::mUndefined
+    )
     {
-        this->rho_ = mass1/this->volume();
+        if (cloud.constProps().constantVolume())
+        {
+            this->rho_ = mass1/this->volume();
+        }
+        else
+        {
+            this->d_ = cbrt(mass1/this->rho_*6/pi);
+        }
     }
     else
     {
-        this->d_ = cbrt(mass1/this->rho_*6.0/pi);
-    }
+        switch (cloud.constProps().volUpdateType())
+        {
+            case constantProperties::volumeUpadteType::mConstRho :
+            {
+                this->d_ = cbrt(mass1/this->rho_*6/pi);
+                break;
+            }
+            case constantProperties::volumeUpadteType::mConstVol :
+            {
+                this->rho_ = mass1/this->volume();
+                break;
+            }
+            case constantProperties::volumeUpadteType::mUpdateRhoAndVol :
+            {
+                scalar deltaVol =
+                    updatedDeltaVolume
+                    (
+                        cloud,
+                        dMassSolid,
+                        pc,
+                        T0
+                    );
 
+                this->rho_ = mass1/(this->volume() + deltaVol);
+                this->d_ = cbrt(mass1/this->rho_*6.0/pi);
+                break;
+            }
+        }
+    }
     // Correct surface values due to emitted species
     this->correctSurfaceValues(cloud, td, Ts, Cs, rhos, mus, Prs, kappas);
     Res = this->Re(rhos, U0, td.Uc(), this->d_, mus);
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingHeterogeneousParcel/ReactingHeterogeneousParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingHeterogeneousParcel/ReactingHeterogeneousParcel.H
index 82f55f5ce4..120cfcd4ab 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingHeterogeneousParcel/ReactingHeterogeneousParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingHeterogeneousParcel/ReactingHeterogeneousParcel.H
@@ -168,6 +168,16 @@ protected:
 
     // Protected Member Functions
 
+        //- Return change of volume due to mass exchange
+        template<class TrackCloudType>
+        scalar updatedDeltaVolume
+        (
+            TrackCloudType& cloud,
+            const scalarField& dMass,
+            const scalar p,
+            const scalar T
+        );
+
 
         //- Calculate surface reactions
         template<class TrackCloudType>
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
index 7c50967678..60084c87ef 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -139,6 +139,48 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::updateMassFractions
 }
 
 
+template<class ParcelType>
+template<class TrackCloudType>
+Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::updatedDeltaVolume
+(
+    TrackCloudType& cloud,
+    const scalarField& dMassGas,
+    const scalarField& dMassLiquid,
+    const scalarField& dMassSolid,
+    const label idG,
+    const label idL,
+    const label idS,
+    const scalar p,
+    const scalar T
+)
+{
+    const auto& props = cloud.composition().phaseProps()[idG];
+    const auto& thermo = cloud.composition().thermo();
+
+    scalarField dVolGas(dMassGas.size(), Zero);
+    forAll(dMassGas, i)
+    {
+        label cid = props.carrierIds()[i];
+        dVolGas[i] = -dMassGas[i]/thermo.carrier().rho(cid, p, T);
+    }
+
+    scalarField dVolLiquid(dMassLiquid.size(), Zero);
+    forAll(dMassLiquid, i)
+    {
+        dVolLiquid[i] =
+            -dMassLiquid[i]/thermo.liquids().properties()[i].rho(p, T);
+    }
+
+    scalarField dVolSolid(dMassSolid.size(), Zero);
+    forAll(dMassSolid, i)
+    {
+        dVolSolid[i] = -dMassSolid[i]/thermo.solids().properties()[i].rho();
+    }
+
+    return (sum(dVolGas) + sum(dVolLiquid) + sum(dMassSolid));
+}
+
+
 // * * * * * * * * * * *  Protected Member Functions * * * * * * * * * * * * //
 
 template<class ParcelType>
@@ -188,6 +230,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     const vector& U0 = this->U_;
     const scalar T0 = this->T_;
     const scalar mass0 = this->mass();
+    const scalar rho0 = this->rho_;
+
 
     const scalar pc = td.pc();
 
@@ -256,9 +300,11 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         d0,
         T0,
         mass0,
+        rho0,
         idL,
         YMix[LIQ],
         YLiquid_,
+        YMix[SLD]*YSolid_,
         dMassPC,
         Sh,
         Ne,
@@ -382,16 +428,57 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
 
     (void)updateMassFractions(mass0, dMassGas, dMassLiquid, dMassSolid);
 
-     // Update particle density or diameter
-    if (cloud.constProps().constantVolume())
+    if
+    (
+        cloud.constProps().volUpdateType()
+     == constantProperties::volumeUpadteType::mUndefined
+    )
     {
-        this->rho_ = mass1/this->volume();
+        if (cloud.constProps().constantVolume())
+        {
+            this->rho_ = mass1/this->volume();
+        }
+        else
+        {
+            this->d_ = cbrt(mass1/this->rho_*6/pi);
+        }
     }
     else
     {
-        this->d_ = cbrt(mass1/this->rho_*6.0/pi);
-    }
+        switch (cloud.constProps().volUpdateType())
+        {
+            case constantProperties::volumeUpadteType::mConstRho :
+            {
+                this->d_ = cbrt(mass1/this->rho_*6/pi);
+                break;
+            }
+            case constantProperties::volumeUpadteType::mConstVol :
+            {
+                this->rho_ = mass1/this->volume();
+                break;
+            }
+            case constantProperties::volumeUpadteType::mUpdateRhoAndVol :
+            {
+                scalar deltaVol =
+                    updatedDeltaVolume
+                    (
+                        cloud,
+                        dMassGas,
+                        dMassLiquid,
+                        dMassSolid,
+                        idG,
+                        idL,
+                        idS,
+                        pc,
+                        T0
+                    );
 
+                this->rho_ = mass1/(this->volume() + deltaVol);
+                this->d_ = cbrt(mass1/this->rho_*6/pi);
+                break;
+            }
+        }
+    }
     // Correct surface values due to emitted species
     this->correctSurfaceValues(cloud, td, Ts, Cs, rhos, mus, Prs, kappas);
     Res = this->Re(rhos, U0, td.Uc(), this->d_, mus);
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
index 1595d9de27..757dc1bc5d 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
@@ -209,6 +209,22 @@ protected:
 
     // Protected Member Functions
 
+
+        //- Return change of volume due to mass exchange
+        template<class TrackCloudType>
+        scalar updatedDeltaVolume
+        (
+            TrackCloudType& cloud,
+            const scalarField& dMassGas,
+            const scalarField& dMassLiquid,
+            const scalarField& dMassSolid,
+            const label idG,
+            const label idL,
+            const label idS,
+            const scalar p,
+            const scalar T
+        );
+
         //- Calculate Devolatilisation
         template<class TrackCloudType>
         void calcDevolatilisation
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
index 3a07a45f4b..8434cf11cc 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -35,6 +36,29 @@ using namespace Foam::constant::mathematical;
 
 // * * * * * * * * * * *  Protected Member Functions * * * * * * * * * * * * //
 
+template<class ParcelType>
+template<class TrackCloudType>
+Foam::scalar Foam::ReactingParcel<ParcelType>::updatedDeltaVolume
+(
+    TrackCloudType& cloud,
+    const scalarField& dMass,
+    const scalar p,
+    const scalar T
+)
+{
+    const auto& composition = cloud.composition();
+
+    scalarField dVolLiquid(dMass.size(), Zero);
+    forAll(dVolLiquid, i)
+    {
+        dVolLiquid[i] =
+            dMass[i]/composition.liquids().properties()[i].rho(p, T);
+    }
+
+    return sum(dVolLiquid);
+}
+
+
 template<class ParcelType>
 template<class TrackCloudType>
 void Foam::ReactingParcel<ParcelType>::calcPhaseChange
@@ -49,9 +73,11 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
     const scalar d,
     const scalar T,
     const scalar mass,
+    const scalar rho,
     const label idPhase,
     const scalar YPhase,
     const scalarField& Y,
+    const scalarField& Ysol,
     scalarField& dMassPC,
     scalar& Sh,
     scalar& N,
@@ -62,7 +88,8 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
     const auto& composition = cloud.composition();
     auto& phaseChange = cloud.phaseChange();
 
-    if (!phaseChange.active() || (YPhase < SMALL))
+    // Some models allow evaporation and condensation
+    if (!phaseChange.active())
     {
         return;
     }
@@ -70,7 +97,6 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
     scalarField X(composition.liquids().X(Y));
 
     scalar Tvap = phaseChange.Tvap(X);
-
     if (T < Tvap)
     {
         return;
@@ -89,16 +115,27 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
         Pr,
         d,
         nus,
+        rho,
         Tdash,
         Tsdash,
         td.pc(),
         td.Tc(),
-        X,
+        X,                // components molar fractions purely in the liquid
+        Ysol*mass,        // total solid mass
+        YPhase*Y*mass,    // total liquid mass
         dMassPC
     );
 
     // Limit phase change mass by availability of each specie
-    dMassPC = min(mass*YPhase*Y, dMassPC);
+    forAll (Y, i)
+    {
+        // evaporation
+        if (dMassPC[i] > 0)
+        {
+            dMassPC[i] = min(mass*YPhase*Y[i], dMassPC[i]);
+        }
+        // condensation Do something?
+    }
 
     const scalar dMassTot = sum(dMassPC);
 
@@ -394,6 +431,7 @@ void Foam::ReactingParcel<ParcelType>::calc
     const vector& U0 = this->U_;
     const scalar T0 = this->T_;
     const scalar mass0 = this->mass();
+    const scalar rho0 = this->rho_;
 
 
     // Calc surface values
@@ -455,9 +493,11 @@ void Foam::ReactingParcel<ParcelType>::calc
         d0,
         T0,
         mass0,
+        rho0,
         0,
         1.0,
         Y_,
+        scalarField(0),
         dMassPC,
         Sh,
         Ne,
@@ -474,14 +514,52 @@ void Foam::ReactingParcel<ParcelType>::calc
 
     this->Cp_ = composition.Cp(0, Y_, td.pc(), T0);
 
-    // Update particle density or diameter
-    if (cloud.constProps().constantVolume())
+    if
+    (
+        cloud.constProps().volUpdateType()
+     == constantProperties::volumeUpadteType::mUndefined
+    )
     {
-        this->rho_ = mass1/this->volume();
+        // Update particle density or diameter
+        if (cloud.constProps().constantVolume())
+        {
+            this->rho_ = mass1/this->volume();
+        }
+        else
+        {
+            this->d_ = cbrt(mass1/this->rho_*6/pi);
+        }
     }
     else
     {
-        this->d_ = cbrt(mass1/this->rho_*6.0/pi);
+        switch (cloud.constProps().volUpdateType())
+        {
+            case constantProperties::volumeUpadteType::mConstRho :
+            {
+                this->d_ = cbrt(mass1/this->rho_*6/pi);
+                break;
+            }
+            case constantProperties::volumeUpadteType::mConstVol :
+            {
+                this->rho_ = mass1/this->volume();
+                break;
+            }
+            case constantProperties::volumeUpadteType::mUpdateRhoAndVol :
+            {
+                scalar deltaVol =
+                    updatedDeltaVolume
+                    (
+                        cloud,
+                        dMass,
+                        td.pc(),
+                        T0
+                    );
+                this->rho_ = mass1/(this->volume() - deltaVol);
+                this->d_ = cbrt(mass1/this->rho_*6/pi);
+                break;
+            }
+
+        }
     }
 
     // Remove the particle when mass falls below minimum threshold
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
index 3ece75e178..74a3245cb5 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
@@ -79,19 +79,39 @@ public:
     static const std::size_t sizeofFields;
 
 
+
     //- Class to hold reacting parcel constant properties
     class constantProperties
     :
         public ParcelType::constantProperties
     {
-        // Private data
 
+    public:
+
+            //- Type of activity coefficient models
+            enum volumeUpadteType
+            {
+                mConstRho,
+                mConstVol,
+                mUpdateRhoAndVol,
+                mUndefined
+            };
+
+    private:
+
+        // Private data
             //- Minimum pressure [Pa]
             demandDrivenEntry<scalar> pMin_;
 
-            //- Constant volume flag - e.g. during mass transfer
+            //- Method to update particle rho and diameter
+            //- 0: constant rho
+            //- 1: constant volume
+            //- 2: recalculation of rho and volume based on the lost mass
             demandDrivenEntry<bool> constantVolume_;
 
+            //- Method to update vol and rho
+            demandDrivenEntry<label> volUpdateType_;
+
 
     public:
 
@@ -114,6 +134,9 @@ public:
 
             //- Return const access to the constant volume flag
             inline bool constantVolume() const;
+
+            //- Return const access to the constant volume flag
+            inline label volUpdateType() const;
     };
 
 
@@ -181,6 +204,17 @@ protected:
 
     // Protected Member Functions
 
+            //- Return change of volume due to mass exchange
+            template<class TrackCloudType>
+            scalar updatedDeltaVolume
+            (
+                TrackCloudType& cloud,
+                const scalarField& dMass,
+                const scalar p,
+                const scalar T
+            );
+
+
         //- Calculate Phase change
         template<class TrackCloudType>
         void calcPhaseChange
@@ -195,9 +229,11 @@ protected:
             const scalar d,            // diameter
             const scalar T,            // temperature
             const scalar mass,         // mass
+            const scalar rho,          // density
             const label idPhase,       // id of phase involved in phase change
             const scalar YPhase,       // total mass fraction
-            const scalarField& YComponents, // component mass fractions
+            const scalarField& YLiq,   // liquid component  mass fractions
+            const scalarField& YSol,   // solid component mass fractions
             scalarField& dMassPC,      // mass transfer - local to parcel
             scalar& Sh,                // explicit parcel enthalpy source
             scalar& N,                 // flux of species emitted from parcel
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H
index 01244daa85..033b34cd10 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H
@@ -34,7 +34,8 @@ Foam::ReactingParcel<ParcelType>::constantProperties::constantProperties()
 :
     ParcelType::constantProperties(),
     pMin_(this->dict_, 0.0),
-    constantVolume_(this->dict_, false)
+    constantVolume_(this->dict_, false),
+    volUpdateType_(this->dict_, mUndefined)
 {}
 
 
@@ -46,7 +47,8 @@ inline Foam::ReactingParcel<ParcelType>::constantProperties::constantProperties
 :
     ParcelType::constantProperties(cp),
     pMin_(cp.pMin_),
-    constantVolume_(cp.constantVolume_)
+    constantVolume_(cp.constantVolume_),
+    volUpdateType_(cp.volUpdateType_)
 {}
 
 
@@ -58,8 +60,54 @@ inline Foam::ReactingParcel<ParcelType>::constantProperties::constantProperties
 :
     ParcelType::constantProperties(parentDict),
     pMin_(this->dict_, "pMin", 1000.0),
-    constantVolume_(this->dict_, word("constantVolume"))
-{}
+    constantVolume_(this->dict_, "constantVolume", false),
+    volUpdateType_(this->dict_, "volumeUpdateMethod")
+{
+    // If constantVolume found use it
+    if (this->dict_.found("constantVolume"))
+    {
+        volUpdateType_.setValue(mUndefined);
+    }
+    else
+    {
+        if (!this->dict_.found("volumeUpdateMethod"))
+        {
+            FatalErrorInFunction
+                << "Either 'constantVolume' or 'volumeUpdateMethod' " << nl
+                << " must be provided. " << nl
+                << " The new keyword is 'volumeUpdateMethod'. " << nl
+                << " Available methods are  : " << nl
+                << " constantRho, constantVolume or updateRhoAndVol. " << nl
+                << " 'constantVolume' is either true/false " << nl
+                << nl << exit(FatalError);
+        }
+        const word volumeUpdateMethod
+        (
+            this->dict_.getWord("volumeUpdateMethod")
+        );
+
+        if (volumeUpdateMethod == "constantRho")
+        {
+            volUpdateType_.setValue(mConstRho);
+        }
+        else if (volumeUpdateMethod == "constantVolume")
+        {
+            volUpdateType_.setValue(mConstVol);
+        }
+        else if (volumeUpdateMethod == "updateRhoAndVol")
+        {
+            volUpdateType_.setValue(mUpdateRhoAndVol);
+        }
+        else
+        {
+            FatalErrorInFunction
+                << "The method provided is not correct " << nl
+                << " Available methods are  : " << nl
+                << " constantRho, constantVolume or updateRhoAndVol. " << nl
+                << nl << exit(FatalError);
+        }
+    }
+}
 
 
 template<class ParcelType>
@@ -155,6 +203,14 @@ Foam::ReactingParcel<ParcelType>::constantProperties::constantVolume() const
 }
 
 
+template<class ParcelType>
+inline Foam::label
+Foam::ReactingParcel<ParcelType>::constantProperties::volUpdateType() const
+{
+    return volUpdateType_.value();
+}
+
+
 // * * * * * * * * * * ThermoParcel Member Functions * * * * * * * * * * * * //
 
 template<class ParcelType>
diff --git a/src/lagrangian/intermediate/parcels/include/makeReactingParcelPhaseChangeModels.H b/src/lagrangian/intermediate/parcels/include/makeReactingParcelPhaseChangeModels.H
index 0e885edb39..2c11a1b699 100644
--- a/src/lagrangian/intermediate/parcels/include/makeReactingParcelPhaseChangeModels.H
+++ b/src/lagrangian/intermediate/parcels/include/makeReactingParcelPhaseChangeModels.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -33,6 +34,7 @@ License
 #include "NoPhaseChange.H"
 #include "LiquidEvaporation.H"
 #include "LiquidEvaporationBoil.H"
+#include "LiquidEvapFuchsKnudsen.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -41,7 +43,8 @@ License
     makePhaseChangeModel(CloudType);                                           \
     makePhaseChangeModelType(NoPhaseChange, CloudType);                        \
     makePhaseChangeModelType(LiquidEvaporation, CloudType);                    \
-    makePhaseChangeModelType(LiquidEvaporationBoil, CloudType);
+    makePhaseChangeModelType(LiquidEvaporationBoil, CloudType);                \
+    makePhaseChangeModelType(LiquidEvapFuchsKnudsen, CloudType);
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.C b/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.C
index fba181f044..e2b88b4f6f 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.C
+++ b/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -261,11 +262,14 @@ Foam::tmp<Foam::scalarField> Foam::CompositionModel<CloudType>::X
             }
             break;
         }
-        default:
+        case phaseProperties::SOLID:
         {
-            FatalErrorInFunction
-                << "Only possible to convert gas and liquid mass fractions"
-                << abort(FatalError);
+            forAll(Y, i)
+            {
+                X[i] = Y[i]/thermo_.solids().properties()[i].W();
+                WInv += X[i];
+            }
+            break;
         }
     }
 
@@ -535,6 +539,67 @@ Foam::scalar Foam::CompositionModel<CloudType>::L
 }
 
 
+template<class CloudType>
+Foam::scalar Foam::CompositionModel<CloudType>::rho
+(
+    const scalarField& Ygas,
+    const scalarField& Yliq,
+    const scalarField& Ysol,
+    const scalar T,
+    const scalar p
+) const
+{
+    const scalarField& YMix = this->YMixture0();
+
+    const auto& carrier = this->carrier();
+    const auto& thermo = this->thermo();
+
+    scalarField Xgas(Ygas.size(), 0);
+    scalar WInv = 0;
+    forAll (Ygas, i)
+    {
+        label cid = phaseProps_[idGas()].carrierIds()[i];
+        Xgas[i] = YMix[idGas()]*Ygas[i]/carrier.W(cid);
+        WInv += Xgas[i];
+    }
+
+    scalarField Xliq(Yliq.size(), 0);
+    forAll (Yliq, i)
+    {
+        Xliq[i] = YMix[idLiquid()]*Yliq[i]/thermo.liquids().properties()[i].W();
+        WInv += Xliq[i];
+    }
+
+    scalarField Xsol(Ysol.size(), 0);
+    forAll (Ysol, i)
+    {
+        Xsol[i] = YMix[idSolid()]*Ysol[i]/thermo.solids().properties()[i].W();
+        WInv += Xsol[i];
+    }
+
+    Xgas /= (WInv + ROOTVSMALL);
+    Xliq /= (WInv + ROOTVSMALL);
+    Xsol /= (WInv + ROOTVSMALL);
+
+
+    scalar rho = 0;
+    forAll (Xgas, i)
+    {
+        label cid = phaseProps_[idGas()].carrierIds()[i];
+        rho += Xgas[i]*carrier.rho(cid, p, T);
+    }
+    forAll (Xliq, i)
+    {
+        rho += Xliq[i]*thermo.liquids().properties()[i].rho(p, T);
+    }
+    forAll (Xsol, i)
+    {
+        rho += Xsol[i]*thermo.solids().properties()[i].rho();
+    }
+
+    return rho;
+
+}
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #include "CompositionModelNew.C"
diff --git a/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H b/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H
index 0da52a8000..a6d9c57f8c 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H
+++ b/src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -265,6 +266,16 @@ public:
                 const scalar p,
                 const scalar T
             ) const;
+
+            //- Return rho of the full composition
+            virtual scalar rho
+            (
+                const scalarField& Ygas,
+                const scalarField& Yliq,
+                const scalarField& Ysol,
+                const scalar T,
+                const scalar p
+            ) const;
 };
 
 
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvapFuchsKnudsen/LiquidEvapFuchsKnudsen.C b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvapFuchsKnudsen/LiquidEvapFuchsKnudsen.C
new file mode 100644
index 0000000000..7e567abd7e
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvapFuchsKnudsen/LiquidEvapFuchsKnudsen.C
@@ -0,0 +1,338 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "LiquidEvapFuchsKnudsen.H"
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::LiquidEvapFuchsKnudsen<CloudType>::calcXcSolution
+(
+    const scalar massliq,
+    const scalar masssol,
+    scalar& Xliq,
+    scalar& Xsol
+) const
+{
+    const scalar Yliq = massliq/(massliq + masssol);
+    const scalar Ysol = 1 - Yliq;
+    Xliq = Yliq/liquids_.properties()[liqToLiqMap_].W();
+    Xsol = Ysol/this->owner().thermo().solids().properties()[solToSolMap_].W();
+    Xliq /= (Xliq + Xsol);
+    Xsol = 1 - Xliq;
+}
+
+
+template<class CloudType>
+Foam::tmp<Foam::scalarField> Foam::LiquidEvapFuchsKnudsen<CloudType>::calcXc
+(
+    const label celli
+) const
+{
+    scalarField Xc(this->owner().thermo().carrier().Y().size());
+
+    forAll(Xc, i)
+    {
+        Xc[i] =
+            this->owner().thermo().carrier().Y()[i][celli]
+           /this->owner().thermo().carrier().W(i);
+    }
+
+    return Xc/sum(Xc);
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::LiquidEvapFuchsKnudsen<CloudType>::Sh
+(
+    const scalar Re,
+    const scalar Sc
+) const
+{
+    return cbrt(1 + Sc*Re)*max(1, pow(Re, 0.077));
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::LiquidEvapFuchsKnudsen<CloudType>::activityCoeff
+(
+    const scalar Xliq,
+    const scalar Xsol
+) const
+{
+    switch (method_)
+    {
+        case pUNIFAC:
+        {
+            FatalErrorInFunction
+                << "Activity coefficient UNIFAC is not implemented " << nl
+                << abort(FatalError);
+            break;
+        }
+        case pHoff:
+        {
+            const scalar ic = this->coeffDict().getScalar("ic");
+            return inv((1 + ic*Xsol/(Xliq + ROOTVSMALL)));
+            break;
+        }
+    }
+    return -1;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::LiquidEvapFuchsKnudsen<CloudType>::LiquidEvapFuchsKnudsen
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    PhaseChangeModel<CloudType>(dict, owner, typeName),
+    method_(pHoff),
+    gamma_(this->coeffDict().getScalar("gamma")),
+    alpha_(this->coeffDict().getScalar("alpham")),
+    liquids_(owner.thermo().liquids()),
+    solution_(this->coeffDict().lookup("solution")),
+    liqToCarrierMap_(-1),
+    liqToLiqMap_(-1),
+    solToSolMap_(-1)
+{
+    if (solution_.size() > 2)
+    {
+        FatalErrorInFunction
+            << "Solution is not well defined. It should be (liquid solid)"
+            << nl << endl;
+    }
+    else
+    {
+        Info<< "Participating liquid-solid species:" << endl;
+
+        Info<< "    " << solution_[0] << endl;
+        liqToCarrierMap_ =
+            owner.composition().carrierId(solution_[0]);
+
+        // Determine mapping between model active liquids and global liquids
+        const label idLiquid = owner.composition().idLiquid();
+        liqToLiqMap_ =
+            owner.composition().localId(idLiquid, solution_[0]);
+
+        // Mapping for the solid
+        const label idSolid = owner.composition().idSolid();
+
+        solToSolMap_ =
+            owner.composition().localId(idSolid, solution_[1]);
+
+        const word activityCoefficientype
+        (
+            this->coeffDict().getWord("activityCoefficient")
+        );
+
+        if (activityCoefficientype == "Hoff")
+        {
+            method_ = pHoff;
+        }
+        else if (activityCoefficientype == "UNIFAC")
+        {
+            method_ = pUNIFAC;
+        }
+        else
+        {
+            FatalErrorInFunction
+                << "activityCoefficient must be either 'Hoff' or 'UNIFAC'"
+                << nl << exit(FatalError);
+        }
+    }
+}
+
+
+template<class CloudType>
+Foam::LiquidEvapFuchsKnudsen<CloudType>::LiquidEvapFuchsKnudsen
+(
+    const LiquidEvapFuchsKnudsen<CloudType>& pcm
+)
+:
+    PhaseChangeModel<CloudType>(pcm),
+    method_(pcm.method_),
+    gamma_(pcm.gamma_),
+    alpha_(pcm.alpha_),
+    liquids_(pcm.owner().thermo().liquids()),
+    solution_(pcm.solution_),
+    liqToCarrierMap_(pcm.liqToCarrierMap_),
+    liqToLiqMap_(pcm.liqToLiqMap_),
+    solToSolMap_(pcm.solToSolMap_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::LiquidEvapFuchsKnudsen<CloudType>::~LiquidEvapFuchsKnudsen()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::LiquidEvapFuchsKnudsen<CloudType>::calculate
+(
+    const scalar dt,
+    const label celli,
+    const scalar Re,
+    const scalar Pr,
+    const scalar d,
+    const scalar nu,
+    const scalar rho,
+    const scalar T,
+    const scalar Ts,
+    const scalar pc,
+    const scalar Tc,
+    const scalarField& X,
+    const scalarField& solMass,
+    const scalarField& liqMass,
+    scalarField& dMassPC
+) const
+{
+    const scalar rhog = this->owner().thermo().thermo().rho()()[celli];
+
+    const label gid = liqToCarrierMap_;
+    const label lid = liqToLiqMap_;
+    const label sid = solToSolMap_;
+
+    const scalar W = liquids_.properties()[lid].W();
+
+    const scalar YeInf = this->owner().thermo().carrier().Y()[gid][celli];
+
+    scalar sigma = liquids_.properties()[lid].sigma(pc, Ts);
+
+    // Kelvin effect
+    const scalar Ke = exp(4*sigma*W/(RR*rho*d*T));
+
+    // vapour diffusivity [m2/s]
+    const scalar Dab = liquids_.properties()[lid].D(pc, Ts);
+
+    // saturation pressure for species i [pa]
+    const scalar pSat = liquids_.properties()[lid].pv(pc, T);
+
+    scalar Xliq(0), Xsol(0);
+    calcXcSolution(liqMass[lid], solMass[sid], Xliq, Xsol);
+
+    // Activity Coefficient (gammaE*Xe)
+    const scalar gamma = activityCoeff(Xliq, Xsol);
+
+    // water concentration at surface
+    const scalar Rliq = RR/W;
+    const scalar YeSurf = max(gamma*Ke*pSat/(Rliq*T*rhog), 0);
+
+    const scalar Kn = 2*gamma_/d;
+    const scalar Cm =
+        (1+Kn)/(1+ (4/(3*alpha_) + 0.377)*Kn + sqr(Kn)*4/(3*alpha_));
+
+    // Schmidt number
+    const scalar Sc = nu/(Dab + ROOTVSMALL);
+
+    // Sherwood number
+    const scalar Sherwood = Sh(Re, Sc);
+
+    // mass flux  [kg/s]
+    const scalar Ni = (rhog*Sherwood*Dab*Cm/d)*log((1 - YeInf)/(1 - YeSurf));
+
+    // mass transfer [kg]
+    dMassPC[lid] += Ni*dt;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::LiquidEvapFuchsKnudsen<CloudType>::dh
+(
+    const label idc,
+    const label idl,
+    const scalar p,
+    const scalar T
+) const
+{
+    scalar dh = 0;
+
+    typedef PhaseChangeModel<CloudType> parent;
+    switch (parent::enthalpyTransfer_)
+    {
+        case (parent::etLatentHeat):
+        {
+            dh = liquids_.properties()[idl].hl(p, T);
+            break;
+        }
+        case (parent::etEnthalpyDifference):
+        {
+            scalar hc = this->owner().composition().carrier().Ha(idc, p, T);
+            scalar hp = liquids_.properties()[idl].h(p, T);
+
+            dh = hc - hp;
+            break;
+        }
+        default:
+        {
+            FatalErrorInFunction
+                << "Unknown enthalpyTransfer type" << abort(FatalError);
+        }
+    }
+
+    return dh;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::LiquidEvapFuchsKnudsen<CloudType>::Tvap
+(
+    const scalarField& X
+) const
+{
+    return Zero;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::LiquidEvapFuchsKnudsen<CloudType>::TMax
+(
+    const scalar p,
+    const scalarField& X
+) const
+{
+    // If liquid is present calculates pvInter
+    if (sum(X) > SMALL)
+    {
+        return liquids_.pvInvert(p, X);
+    }
+    else
+    {
+        return GREAT;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvapFuchsKnudsen/LiquidEvapFuchsKnudsen.H b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvapFuchsKnudsen/LiquidEvapFuchsKnudsen.H
new file mode 100644
index 0000000000..96bb45ee4d
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvapFuchsKnudsen/LiquidEvapFuchsKnudsen.H
@@ -0,0 +1,214 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::LiquidEvapFuchsKnudsen
+
+Group
+    grpLagrangianIntermediatePhaseChangeSubModels
+
+Description
+    Liquid evaporation/condensation model for solution of liquid and solid.
+
+    This model takes into account the Fuchs-Knudsen number correction, the
+    modified Raoult's law is used to obtain the concenration of the evapora
+    ble component on the surface and the activity coefficient is used.
+    The correction Kelvin effect is used.
+
+    Reference:
+    \verbatim
+        Xiaole Chen, Yu Feng, Wenqi Zhon, Clement Kleinstreuer.
+        Numerical investigation of the interaction, transport and deposition
+        of multicomponent droplets in a a simple mouth-throat model.
+        Journal of Aerosol Science, 105(2017), 108-127.
+        DOI:10.1016/j.jaerosci.2016.12.001
+    \endverbatim
+
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef LiquidEvapFuchsKnudsen_H
+#define LiquidEvapFuchsKnudsen_H
+
+#include "PhaseChangeModel.H"
+#include "liquidMixtureProperties.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                     Class LiquidEvapFuchsKnudsen Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class LiquidEvapFuchsKnudsen
+:
+    public PhaseChangeModel<CloudType>
+{
+public:
+
+    // Public Enumerations
+
+        //- Type of activity coefficient models
+        enum activityCoeffMethodType
+        {
+            pUNIFAC,
+            pHoff
+        };
+
+
+protected:
+
+    // Protected Data
+
+        //- Method used
+        activityCoeffMethodType method_;
+
+        //- Mean gas free path
+        scalar gamma_;
+
+        //- The mass thermal accomodation
+        scalar alpha_;
+
+        //- Global liquid properties data
+        const liquidMixtureProperties& liquids_;
+
+        //- List of active liquid names i.e (liquidName solidName)
+        List<word> solution_;
+
+        //- Mapping between liquid and carrier species
+        label liqToCarrierMap_;
+
+        //- Mapping between local and global liquid species
+        label liqToLiqMap_;
+
+        //- Mapping between local and global solid species
+        label solToSolMap_;
+
+
+    // Protected Member Functions
+
+        //- Sherwood number as a function of Reynolds and Schmidt numbers
+        scalar Sh(const scalar Re, const scalar Sc) const;
+
+        //- Calculate the carrier phase component volume fractions at celli
+        tmp<scalarField> calcXc(const label celli) const;
+
+        //- Calculate volumetric fractions of components in the solution
+        void calcXcSolution
+        (
+            const scalar massliq,
+            const scalar masssol,
+            scalar& Xliq,
+            scalar& Xsol
+        ) const;
+
+        //- Return activity coefficient
+        scalar activityCoeff(const scalar Xliq, const scalar Ysol) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("liquidEvapFuchsKnudsen");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        LiquidEvapFuchsKnudsen(const dictionary& dict, CloudType& cloud);
+
+        //- Construct copy
+        LiquidEvapFuchsKnudsen(const LiquidEvapFuchsKnudsen<CloudType>& pcm);
+
+        //- Construct and return a clone
+        virtual autoPtr<PhaseChangeModel<CloudType>> clone() const
+        {
+            return autoPtr<PhaseChangeModel<CloudType>>
+            (
+                new LiquidEvapFuchsKnudsen<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~LiquidEvapFuchsKnudsen();
+
+
+    // Member Functions
+
+        //- Update model
+        virtual void calculate
+        (
+            const scalar dt,
+            const label celli,
+            const scalar Re,
+            const scalar Pr,
+            const scalar d,
+            const scalar nu,
+            const scalar rho,
+            const scalar T,
+            const scalar Ts,
+            const scalar pc,
+            const scalar Tc,
+            const scalarField& X,
+            const scalarField& Xsol,
+            const scalarField& liqMass,
+            scalarField& dMassPC
+        ) const;
+
+        //- Return the enthalpy per unit mass
+        virtual scalar dh
+        (
+            const label idc,
+            const label idl,
+            const scalar p,
+            const scalar T
+        ) const;
+
+        //- Return vapourisation temperature
+        virtual scalar Tvap(const scalarField& X) const;
+
+        //- Return maximum/limiting temperature
+        virtual scalar TMax(const scalar p, const scalarField& X) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "LiquidEvapFuchsKnudsen.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
index 67e8d9136d..61c0491bb1 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -139,11 +140,14 @@ void Foam::LiquidEvaporation<CloudType>::calculate
     const scalar Pr,
     const scalar d,
     const scalar nu,
+    const scalar rho,
     const scalar T,
     const scalar Ts,
     const scalar pc,
     const scalar Tc,
     const scalarField& X,
+    const scalarField& solMass,
+    const scalarField& liqMass,
     scalarField& dMassPC
 ) const
 {
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.H b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.H
index 6c994d34bc..b038e10240 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.H
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -119,11 +120,14 @@ public:
             const scalar Pr,
             const scalar d,
             const scalar nu,
+            const scalar rho,
             const scalar T,
             const scalar Ts,
             const scalar pc,
             const scalar Tc,
             const scalarField& X,
+            const scalarField& solMass,
+            const scalarField& liqMass,
             scalarField& dMassPC
         ) const;
 
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.C b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.C
index 9351adc974..911f4b919b 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.C
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -139,11 +140,14 @@ void Foam::LiquidEvaporationBoil<CloudType>::calculate
     const scalar Pr,
     const scalar d,
     const scalar nu,
+    const scalar rho,
     const scalar T,
     const scalar Ts,
     const scalar pc,
     const scalar Tc,
     const scalarField& X,
+    const scalarField& solMass,
+    const scalarField& liqMass,
     scalarField& dMassPC
 ) const
 {
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.H b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.H
index d49eb55774..49e6e879e8 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.H
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -129,11 +130,14 @@ public:
             const scalar Pr,
             const scalar d,
             const scalar nu,
+            const scalar rho,
             const scalar T,
             const scalar Ts,
             const scalar pc,
             const scalar Tc,
             const scalarField& X,
+            const scalarField& solMass,
+            const scalarField& liqMass,
             scalarField& dMassPC
         ) const;
 
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/NoPhaseChange/NoPhaseChange.C b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/NoPhaseChange/NoPhaseChange.C
index 8c7bba25ca..6b736e125e 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/NoPhaseChange/NoPhaseChange.C
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/NoPhaseChange/NoPhaseChange.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -75,11 +76,14 @@ void Foam::NoPhaseChange<CloudType>::calculate
     const scalar Pr,
     const scalar d,
     const scalar nu,
+    const scalar rho,
     const scalar T,
     const scalar Ts,
     const scalar pc,
     const scalar Tc,
     const scalarField& X,
+    const scalarField& solMass,
+    const scalarField& liqMass,
     scalarField& dMassPC
 ) const
 {
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/NoPhaseChange/NoPhaseChange.H b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/NoPhaseChange/NoPhaseChange.H
index 395414b598..577a204264 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/NoPhaseChange/NoPhaseChange.H
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/NoPhaseChange/NoPhaseChange.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -94,11 +95,14 @@ public:
             const scalar Pr,
             const scalar d,
             const scalar nu,
+            const scalar rho,
             const scalar T,
             const scalar Ts,
             const scalar pc,
             const scalar Tc,
             const scalarField& X,
+            const scalarField& solMass,
+            const scalarField& liqMass,
             scalarField& dMassPC
         ) const;
 };
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.H b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.H
index e31c90821c..70a40083df 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.H
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/PhaseChangeModel/PhaseChangeModel.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -166,11 +167,14 @@ public:
             const scalar Pr,
             const scalar d,
             const scalar nu,
+            const scalar rho,
             const scalar T,
             const scalar Ts,
             const scalar pc,
             const scalar Tc,
             const scalarField& X,
+            const scalarField& solMass,
+            const scalarField& liqMass,
             scalarField& dMassPC
         ) const = 0;
 
diff --git a/src/thermophysicalModels/thermophysicalProperties/liquidProperties/liquidMixtureProperties/liquidMixtureProperties.C b/src/thermophysicalModels/thermophysicalProperties/liquidProperties/liquidMixtureProperties/liquidMixtureProperties.C
index 1dd8b8d577..acea47583f 100644
--- a/src/thermophysicalModels/thermophysicalProperties/liquidProperties/liquidMixtureProperties/liquidMixtureProperties.C
+++ b/src/thermophysicalModels/thermophysicalProperties/liquidProperties/liquidMixtureProperties/liquidMixtureProperties.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -110,7 +111,7 @@ Foam::scalar Foam::liquidMixtureProperties::Tc(const scalarField& X) const
         vTc += x1*properties_[i].Tc();
     }
 
-    return vTc/vc;
+    return vTc/ (vc + ROOTVSMALL);
 }
 
 
@@ -258,7 +259,7 @@ Foam::scalarField Foam::liquidMixtureProperties::Y(const scalarField& X) const
         sumY += Y[i];
     }
 
-    Y /= sumY;
+    Y /= (sumY + ROOTVSMALL);
 
     return Y;
 }
@@ -275,7 +276,7 @@ Foam::scalarField Foam::liquidMixtureProperties::X(const scalarField& Y) const
         sumX += X[i];
     }
 
-    X /= sumX;
+    X /= (sumX + ROOTVSMALL);
 
     return X;
 }
@@ -307,7 +308,7 @@ Foam::scalar Foam::liquidMixtureProperties::rho
         }
     }
 
-    return sumY/v;
+    return sumY/(v + ROOTVSMALL);
 }
 
 
@@ -333,7 +334,7 @@ Foam::scalar Foam::liquidMixtureProperties::pv
         }
     }
 
-    return pv/sumY;
+    return pv/(sumY + ROOTVSMALL);
 }
 
 
@@ -359,7 +360,7 @@ Foam::scalar Foam::liquidMixtureProperties::hl
         }
     }
 
-    return hl/sumY;
+    return hl/(sumY + ROOTVSMALL);
 }
 
 
@@ -385,7 +386,7 @@ Foam::scalar Foam::liquidMixtureProperties::Cp
         }
     }
 
-    return Cp/sumY;
+    return Cp/(sumY + ROOTVSMALL);
 }
 
 
@@ -411,7 +412,7 @@ Foam::scalar Foam::liquidMixtureProperties::sigma
         XsSum += Xs[i];
     }
 
-    Xs /= XsSum;
+    Xs /= (XsSum + ROOTVSMALL);
 
     forAll(properties_, i)
     {
@@ -468,7 +469,7 @@ Foam::scalar Foam::liquidMixtureProperties::kappa
         pSum += phii[i];
     }
 
-    phii /= pSum;
+    phii /= (pSum + ROOTVSMALL);
 
     scalar K = 0.0;
 
@@ -513,7 +514,7 @@ Foam::scalar Foam::liquidMixtureProperties::D
         }
     }
 
-    return 1.0/Dinv;
+    return 1.0/(Dinv + ROOTVSMALL);
 }
 
 
diff --git a/tutorials/lagrangian/reactingHeterogenousParcelFoam/rectangularDuct/constant/reactingCloud1Properties b/tutorials/lagrangian/reactingHeterogenousParcelFoam/rectangularDuct/constant/reactingCloud1Properties
index cf72e76116..a45d80fa85 100644
--- a/tutorials/lagrangian/reactingHeterogenousParcelFoam/rectangularDuct/constant/reactingCloud1Properties
+++ b/tutorials/lagrangian/reactingHeterogenousParcelFoam/rectangularDuct/constant/reactingCloud1Properties
@@ -64,7 +64,7 @@ constantProperties
     Cp0             850; //Initial particle Cp (overwritten by composition)
 
     hRetentionCoeff 0;
-    constantVolume  true;
+    volumeUpdateMethod  constantVolume;
 }
 
 
diff --git a/tutorials/lagrangian/reactingParcelFoam/parcelInBox/constant/reactingCloud1Properties b/tutorials/lagrangian/reactingParcelFoam/parcelInBox/constant/reactingCloud1Properties
index 2f5a91d527..e4b9a76324 100644
--- a/tutorials/lagrangian/reactingParcelFoam/parcelInBox/constant/reactingCloud1Properties
+++ b/tutorials/lagrangian/reactingParcelFoam/parcelInBox/constant/reactingCloud1Properties
@@ -60,7 +60,7 @@ constantProperties
     T0              350;
     Cp0             4100;
 
-    constantVolume  false;
+    volumeUpdateMethod  constantRho;
 }
 
 
diff --git a/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/H2O b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/H2O
new file mode 100644
index 0000000000..ea2cb03310
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/H2O
@@ -0,0 +1,42 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1806                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      H2O;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0.01;
+
+boundaryField
+{
+     walls
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 0.01;
+        value           uniform 0.01;
+    }
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 0.01;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/N2 b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/N2
new file mode 100644
index 0000000000..7282998db7
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/N2
@@ -0,0 +1,42 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1806                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      N2;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0.8;
+
+boundaryField
+{
+    walls
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        value           uniform 0.8;
+        inletValue      uniform 0.8;
+    }
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 0.8;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/O2 b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/O2
new file mode 100644
index 0000000000..f6ff93daa1
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/O2
@@ -0,0 +1,42 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1806                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      O2;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0.19;
+
+boundaryField
+{
+    walls
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 0.19;
+        value           uniform 0.2;
+    }
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 0.19;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/T b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/T
new file mode 100644
index 0000000000..f72426fdc0
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/T
@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1806                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 300.0;
+
+boundaryField
+{
+    walls
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 300.0;
+    }
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 300.0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/U b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/U
new file mode 100644
index 0000000000..a010d9d094
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/U
@@ -0,0 +1,42 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus.master.develop                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    location    "0";
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0.1 0 0);
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+    outlet
+    {
+        type            pressureInletOutletVelocity;
+        value           $internalField;
+    }
+    walls
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/alphat b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/alphat
new file mode 100644
index 0000000000..8c3ac9e91e
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/alphat
@@ -0,0 +1,43 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1806                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alphat;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    outlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    walls
+    {
+        type            compressible::alphatWallFunction;
+        Prt             0.85;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/k b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/k
new file mode 100644
index 0000000000..4189f4324b
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/k
@@ -0,0 +1,43 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1806                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 3.75e-9;
+
+boundaryField
+{
+    inlet
+    {
+        type            turbulentIntensityKineticEnergyInlet;
+        intensity       0.16;
+        value           uniform 3.75e-9;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 3.75e-9;
+    }
+    walls
+    {
+        type            kqRWallFunction;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/nut b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/nut
new file mode 100644
index 0000000000..9ca9214167
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/nut
@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1806                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      nut;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    outlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    walls
+    {
+        type            nutkWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/omega b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/omega
new file mode 100644
index 0000000000..8c018fb0e5
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/omega
@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1806                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      omega;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 -1 0 0 0 0];
+
+internalField   uniform 4.5e-3;
+
+boundaryField
+{
+    inlet
+    {
+        type            turbulentMixingLengthFrequencyInlet;
+        mixingLength    0.007;
+        k               k;
+        value           uniform 4.5e-3;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 4.5e-3;
+    }
+    walls
+    {
+        type            omegaWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           $internalField;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/p b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/p
new file mode 100644
index 0000000000..092220adeb
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/p
@@ -0,0 +1,42 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1806                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 100000;
+
+boundaryField
+{
+    inlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+    outlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+    walls
+    {
+        type            calculated;
+        value           $internalField;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/p_rgh b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/p_rgh
new file mode 100644
index 0000000000..2b8a79c028
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/0.orig/p_rgh
@@ -0,0 +1,40 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1806                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 1e5;
+
+boundaryField
+{
+    inlet
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+        type            totalPressure;
+        p0              $internalField;
+    }
+    walls
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/Allclean b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/Allclean
new file mode 100755
index 0000000000..fb1f384730
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/Allclean
@@ -0,0 +1,8 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/CleanFunctions      # Tutorial clean functions
+#------------------------------------------------------------------------------
+
+cleanCase0
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/Allrun b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/Allrun
new file mode 100755
index 0000000000..4749cc3693
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/Allrun
@@ -0,0 +1,16 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+restore0Dir
+
+runApplication blockMesh
+
+runApplication potentialFoam
+# Remove incompatible (volumetric) flux field
+\rm -f 0/phi 2>/dev/null
+
+runApplication $(getApplication)
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/constant/chemistryProperties b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/constant/chemistryProperties
new file mode 100644
index 0000000000..e4a6386238
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/constant/chemistryProperties
@@ -0,0 +1,28 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1806                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      chemistryProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+chemistryType
+{
+    solver            noChemistrySolver;
+}
+
+chemistry       off;
+
+initialChemicalTimeStep 1e-07;  // NOT USED
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/constant/combustionProperties b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/constant/combustionProperties
new file mode 100644
index 0000000000..f49b11d2e3
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/constant/combustionProperties
@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1806                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      combustionProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+combustionModel none;
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/constant/foam.dat b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/constant/foam.dat
new file mode 100644
index 0000000000..db40f0010c
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/constant/foam.dat
@@ -0,0 +1,82 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1806                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      foam.dat;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+O2
+{
+    specie
+    {
+        molWeight       31.9988;
+    }
+    thermodynamics
+    {
+        Tlow            200;
+        Thigh           5000;
+        Tcommon         1000;
+        highCpCoeffs    ( 3.69758 0.00061352 -1.25884e-07 1.77528e-11 -1.13644e-15 -1233.93 3.18917 );
+        lowCpCoeffs     ( 3.21294 0.00112749 -5.75615e-07 1.31388e-09 -8.76855e-13 -1005.25 6.03474 );
+    }
+    transport
+    {
+        As              1.67212e-06;
+        Ts              170.672;
+    }
+}
+
+H2O
+{
+    specie
+    {
+        molWeight       18.0153;
+    }
+    thermodynamics
+    {
+        Tlow            200;
+        Thigh           5000;
+        Tcommon         1000;
+        highCpCoeffs    ( 2.67215 0.00305629 -8.73026e-07 1.201e-10 -6.39162e-15 -29899.2 6.86282 );
+        lowCpCoeffs     ( 3.38684 0.00347498 -6.3547e-06 6.96858e-09 -2.50659e-12 -30208.1 2.59023 );
+    }
+    transport
+    {
+        As              1.67212e-06;
+        Ts              170.672;
+    }
+}
+
+N2
+{
+    specie
+    {
+        molWeight       28.0134;
+    }
+    thermodynamics
+    {
+        Tlow            200;
+        Thigh           5000;
+        Tcommon         1000;
+        highCpCoeffs    ( 2.92664 0.00148798 -5.68476e-07 1.0097e-10 -6.75335e-15 -922.798 5.98053 );
+        lowCpCoeffs     ( 3.29868 0.00140824 -3.96322e-06 5.64152e-09 -2.44486e-12 -1020.9 3.95037 );
+    }
+    transport
+    {
+        As              1.67212e-06;
+        Ts              170.672;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/constant/foam.inp b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/constant/foam.inp
new file mode 100644
index 0000000000..2189d7c94b
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/constant/foam.inp
@@ -0,0 +1,10 @@
+species
+(
+    O2
+    N2
+    H2O
+)
+;
+
+reactions
+{}
diff --git a/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/constant/g b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/constant/g
new file mode 100644
index 0000000000..037697dd98
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/constant/g
@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1806                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       uniformDimensionedVectorField;
+    location    "constant";
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           (0 -9.81 0);
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/constant/radiationProperties b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/constant/radiationProperties
new file mode 100644
index 0000000000..cb46544fd0
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/constant/radiationProperties
@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1806                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      radiationProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+radiation       off;
+
+radiationModel  none;
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/constant/reactingCloud1Properties b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/constant/reactingCloud1Properties
new file mode 100644
index 0000000000..03c36f64f8
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/constant/reactingCloud1Properties
@@ -0,0 +1,175 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1806                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      reactingCloud1Properties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solution
+{
+    active          yes;
+    coupled         true;
+    transient       yes;
+
+    maxCo           0.3;
+    cellValueSourceCorrection off;
+
+    sourceTerms
+    {
+        resetOnStartup  false;
+        schemes
+        {
+            rho             explicit 1;
+            U               explicit 1;
+            Yi              explicit 1;
+            h               explicit 1;
+            radiation       explicit 1;
+        }
+    }
+
+    interpolationSchemes
+    {
+        rho             cell;
+        U               cellPoint;
+        thermo:mu       cell;
+        T               cell;
+        Cp              cell;
+        kappa           cell;
+        p               cell;
+    }
+
+    integrationSchemes
+    {
+        U               Euler;
+        T               analytical;
+    }
+}
+
+
+constantProperties
+{
+    //Initial particle density, used if defined, if not  by composition
+    //rho0            5100;
+    T0              310; //Initial particle temperature
+    Cp0             850; //Initial particle Cp (overwritten by composition)
+
+    hRetentionCoeff 0;
+    volumeUpdateMethod  updateRhoAndVol; //constantRho, constantVolume,
+}
+
+
+subModels
+{
+    particleForces
+    {
+        //sphereDrag;
+        //gravity;
+    }
+
+    injectionModels
+    {
+        model1
+        {
+            type            patchInjection;
+            patch           inlet;
+            parcelBasisType mass;
+            U0              (0.4 0 0);
+            massTotal        1.38;
+            parcelsPerSecond 123;
+            SOI              0;
+            duration         1;
+            flowRateProfile  constant 1;
+
+            sizeDistribution
+            {
+                type        fixedValue;
+                fixedValueDistribution
+                {
+                    value       0.011;
+                }
+            }
+        }
+    }
+
+    dispersionModel none;
+
+    patchInteractionModel standardWallInteraction;
+
+    heatTransferModel RanzMarshall;
+
+    compositionModel singleMixtureFraction;
+
+    phaseChangeModel liquidEvapFuchsKnudsen;
+
+    liquidEvapFuchsKnudsenCoeffs
+    {
+        gamma               6.8e-8;     // Mean gas free path
+        alpham              1;          // The mass thermal accomodation
+        solution            (H2O NaCl); // Solution (liquid solid)
+
+        activityCoefficient Hoff;
+        ic                  1.85;
+        enthalpyTransfer    enthalpyDifference;
+    }
+
+    devolatilisationModel none;
+
+    surfaceReactionModel none;
+
+    stochasticCollisionModel none;
+
+    surfaceFilmModel none;
+
+    radiation       off;
+
+    standardWallInteractionCoeffs
+    {
+        type            rebound;
+    }
+
+    RanzMarshallCoeffs
+    {
+        BirdCorrection  off;
+    }
+
+    heterogeneousReactingModel  none;
+
+    singleMixtureFractionCoeffs
+    {
+        phases
+        (
+            gas
+            {
+            }
+            liquid
+            {
+                H2O  1;
+            }
+            solid
+            {
+                NaCl 1;
+            }
+        );
+        YGasTot0        0;
+        YLiquidTot0     0.0;
+        YSolidTot0      1;
+    }
+}
+
+
+cloudFunctions
+{
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/constant/thermophysicalProperties b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/constant/thermophysicalProperties
new file mode 100644
index 0000000000..52e0d23543
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/constant/thermophysicalProperties
@@ -0,0 +1,58 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1806                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         reactingMixture;
+    transport       sutherland;
+    thermo          janaf;
+    energy          sensibleEnthalpy;
+    equationOfState perfectGas;
+    specie          specie;
+}
+
+chemistryReader foamChemistryReader;
+
+foamChemistryFile "<constant>/foam.inp";
+
+foamChemistryThermoFile "<constant>/foam.dat";
+
+inertSpecie   N2;
+
+liquids
+{
+    H2O;
+}
+
+solids
+{
+    NaCl
+    {
+        defaultCoeffs   no;
+        rho             2165;
+        Cp              850;
+        kappa           0.04;
+        Hf              0;
+        emissivity      1.0;
+        W               58;
+    }
+}
+
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/constant/turbulenceProperties b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/constant/turbulenceProperties
new file mode 100644
index 0000000000..3b314c3574
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/constant/turbulenceProperties
@@ -0,0 +1,29 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1806                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  laminar;
+
+RAS
+{
+    RASModel            kOmegaSST; // kEpsilon;
+
+    turbulence          on;
+
+    printCoeffs         on;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/system/blockMeshDict b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/system/blockMeshDict
new file mode 100644
index 0000000000..87cbf30079
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/system/blockMeshDict
@@ -0,0 +1,75 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1806                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version         2.0;
+    format          ascii;
+    class           dictionary;
+    location    "system";
+    object          blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scale   1;
+
+vertices
+(
+    (0 0 0)
+    (90 0 0)
+    (90 10 0)
+    (0 10 0)
+    (0 0 10)
+    (90 0 10)
+    (90 10 10)
+    (0 10 10)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (100 40 10) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    inlet
+    {
+        type patch;
+        faces
+        (
+            (0 4 7 3)
+        );
+    }
+
+    outlet
+    {
+        type patch;
+        faces
+        (
+            (2 6 5 1)
+        );
+    }
+
+    walls
+    {
+        type wall;
+        faces
+        (
+            (3 7 6 2)
+            (1 5 4 0)
+            (0 3 2 1)
+            (4 5 6 7)
+        );
+    }
+);
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/system/controlDict b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/system/controlDict
new file mode 100644
index 0000000000..4d5d49b286
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/system/controlDict
@@ -0,0 +1,59 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1806                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     reactingParcelFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         15;
+
+deltaT          8.13e-3;
+
+writeControl    adjustableRunTime;
+
+writeInterval   0.5;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  10;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+adjustTimeStep  yes;
+
+maxCo           1.0;
+
+maxDeltaT       1;
+
+functions
+{
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/system/decomposeParDict b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/system/decomposeParDict
new file mode 100644
index 0000000000..6cad762db7
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/system/decomposeParDict
@@ -0,0 +1,30 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1806                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 4;
+
+method          scotch;
+
+coeffs
+{
+    n           (2 2 1);
+}
+
+distributed     no;
+
+roots           ( );
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/system/fvSchemes b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/system/fvSchemes
new file mode 100644
index 0000000000..16f8cdf93a
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/system/fvSchemes
@@ -0,0 +1,63 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1806                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+    div(phi,U)      Gauss upwind;
+    div(phid,p)     Gauss upwind;
+    div(phi,K)      Gauss linear;
+    div(phi,h)       Gauss upwind;
+    div(phi,k)       Gauss upwind;
+    div(phi,epsilon)  Gauss upwind;
+    div(phi,omega)  Gauss upwind;
+    div(phi,Yi_h)   Gauss upwind;
+    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear uncorrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         uncorrected;
+}
+
+wallDist
+{
+    method meshWave;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/system/fvSolution b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/system/fvSolution
new file mode 100644
index 0000000000..09949a4de8
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/rectangularChannel/system/fvSolution
@@ -0,0 +1,112 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1806                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    rho
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-05;
+        relTol          0.1;
+    }
+
+    rhoFinal
+    {
+        $rho;
+        tolerance       1e-05;
+        relTol          0;
+    }
+
+    "(U|k|epsilon)"
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-06;
+        relTol          0.1;
+    }
+
+    "(U|k|epsilon|omega)Final"
+    {
+        $U;
+        tolerance       1e-06;
+        relTol          0;
+    }
+
+    p_rgh
+    {
+        solver          GAMG;
+        tolerance       0;
+        relTol          0.1;
+        smoother        GaussSeidel;
+    }
+
+    p_rghFinal
+    {
+        $p_rgh;
+        tolerance       1e-06;
+        relTol          0;
+    }
+
+    "(Yi|O2|N2|H2O)"
+    {
+        solver          PBiCGStab;
+        preconditioner  DILU;
+        tolerance       1e-6;
+        relTol          0;
+    }
+
+    h
+    {
+        $Yi;
+        relTol          0.1;
+    }
+
+    hFinal
+    {
+        $Yi;
+    }
+}
+
+potentialFlow
+{
+    nNonOrthogonalCorrectors 5;
+}
+
+
+PIMPLE
+{
+    transonic       no;
+    nOuterCorrectors 1;
+    nCorrectors     2;
+    nNonOrthogonalCorrectors 0;
+    momentumPredictor yes;
+}
+
+relaxationFactors
+{
+    fields
+    {
+    }
+    equations
+    {
+        ".*"            1;
+    }
+}
+
+
+// ************************************************************************* //