ENH: new ReactingHeterogeneousCloud

- number of particles per parcel info to kinematic cloud

- added turbulent dispersion to basicHeterogeneousReactingParcel

- corrected dhsTrans in MUCSheterogeneousRate::calculate

- added cloud macro system to reactingParcelFoam and fixed calculation
  of average particles per parcel

- added progress variable dimension to reacting model (nF)

- added ReactingHeterogeneous tutorial

applications/solvers/lagrangian/reactingParcelFoam/createClouds.H

@@ -1,5 +1,5 @@
 Info<< "\nConstructing reacting cloud" << endl;
-basicReactingMultiphaseCloud parcels
+basicReactingTypeCloud parcels
 (
     "reactingCloud1",
     rho,

applications/solvers/lagrangian/reactingParcelFoam/reactingHeterogenousParcelFoam/Make/files

@@ -0,0 +1,3 @@
+reactingHeterogenousParcelFoam.C
+
+EXE = $(FOAM_APPBIN)/reactingHeterogenousParcelFoam

applications/solvers/lagrangian/reactingParcelFoam/reactingHeterogenousParcelFoam/Make/options

@@ -0,0 +1,48 @@
+EXE_INC = \
+    -I. \
+    -I.. \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I${LIB_SRC}/sampling/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/ODE/lnInclude \
+    -I$(LIB_SRC)/combustionModels/lnInclude \
+    -I$(FOAM_SOLVERS)/combustion/reactingFoam
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lfvOptions \
+    -lsampling \
+    -lmeshTools \
+    -lturbulenceModels \
+    -lcompressibleTurbulenceModels \
+    -lspecie \
+    -lcompressibleTransportModels \
+    -lfluidThermophysicalModels \
+    -lreactionThermophysicalModels \
+    -lSLGThermo \
+    -lchemistryModel \
+    -lregionModels \
+    -lradiationModels \
+    -lsurfaceFilmModels \
+    -lsurfaceFilmDerivedFvPatchFields \
+    -llagrangian \
+    -llagrangianIntermediate \
+    -llagrangianTurbulence \
+    -lODE \
+    -lcombustionModels

applications/solvers/lagrangian/reactingParcelFoam/reactingHeterogenousParcelFoam/reactingHeterogenousParcelFoam.C

@@ -0,0 +1,43 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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/>.
+
+Application
+    reactingHeterogenousParcelFoam
+
+Group
+    grpLagrangianSolvers
+
+Description
+    Transient solver for the coupled transport of a single kinematic particle
+    cloud including the effect of the volume fraction of particles on the
+    continuous phase. Multi-Phase Particle In Cell (MPPIC) modeling is used to
+    represent collisions without resolving particle-particle interactions.
+
+\*---------------------------------------------------------------------------*/
+
+#define CLOUD_BASE_TYPE HeterogeneousReacting
+#define CLOUD_BASE_TYPE_NAME "HeterogeneousReacting"
+
+#include "reactingParcelFoam.C"
+
+// ************************************************************************* //

applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2017 OpenFOAM Foundation
@@ -37,7 +37,7 @@ Description
 
 #include "fvCFD.H"
 #include "turbulentFluidThermoModel.H"
-#include "basicReactingMultiphaseCloud.H"
+
 #include "surfaceFilmModel.H"
 #include "rhoReactionThermo.H"
 #include "CombustionModel.H"
@@ -48,6 +48,15 @@ Description
 #include "pressureControl.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
+#include "cloudMacros.H"
+
+#ifndef CLOUD_BASE_TYPE
+    #define CLOUD_BASE_TYPE ReactingMultiphase
+    #define CLOUD_BASE_TYPE_NAME "reacting"
+#endif
+
+#include CLOUD_INCLUDE_FILE(CLOUD_BASE_TYPE)
+#define basicReactingTypeCloud CLOUD_TYPE(CLOUD_BASE_TYPE)
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -91,20 +91,14 @@ Foam::COxidationDiffusionLimitedRate<CloudType>::COxidationDiffusionLimitedRate
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::COxidationDiffusionLimitedRate<CloudType>::
-~COxidationDiffusionLimitedRate()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class CloudType>
 Foam::scalar Foam::COxidationDiffusionLimitedRate<CloudType>::calculate
 (
     const scalar dt,
+    const scalar Re,
+    const scalar nu,
     const label celli,
     const scalar d,
     const scalar T,

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -50,7 +50,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// Forward class declarations
+// Forward declarations
 template<class CloudType>
 class COxidationDiffusionLimitedRate;
 
@@ -130,7 +130,7 @@ public:
 
 
     //- Destructor
-    virtual ~COxidationDiffusionLimitedRate();
+    virtual ~COxidationDiffusionLimitedRate() = default;
 
 
     // Member Functions
@@ -139,6 +139,8 @@ public:
         virtual scalar calculate
         (
             const scalar dt,
+            const scalar Re,
+            const scalar nu,
             const label celli,
             const scalar d,
             const scalar T,

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2012-2016 OpenFOAM Foundation
@@ -91,19 +91,14 @@ Foam::COxidationHurtMitchell<CloudType>::COxidationHurtMitchell
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::COxidationHurtMitchell<CloudType>::~COxidationHurtMitchell()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class CloudType>
 Foam::scalar Foam::COxidationHurtMitchell<CloudType>::calculate
 (
     const scalar dt,
+    const scalar Re,
+    const scalar nu,
     const label celli,
     const scalar d,
     const scalar T,

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2012-2016 OpenFOAM Foundation
@@ -73,9 +73,9 @@ class COxidationHurtMitchell
 :
     public SurfaceReactionModel<CloudType>
 {
-    // Private data
+    // Private Data
 
-        // Model constants
+        // Model Constants
 
             //- Stoichiometry of reaction
             const scalar Sb_;
@@ -143,7 +143,7 @@ public:
 
 
     //- Destructor
-    virtual ~COxidationHurtMitchell();
+    virtual ~COxidationHurtMitchell() = default;
 
 
     // Member Functions
@@ -152,6 +152,8 @@ public:
         virtual scalar calculate
         (
             const scalar dt,
+            const scalar Re,
+            const scalar nu,
             const label celli,
             const scalar d,
             const scalar T,

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2014-2016 OpenFOAM Foundation
@@ -103,20 +103,14 @@ Foam::COxidationIntrinsicRate<CloudType>::COxidationIntrinsicRate
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::COxidationIntrinsicRate<CloudType>::
-~COxidationIntrinsicRate()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class CloudType>
 Foam::scalar Foam::COxidationIntrinsicRate<CloudType>::calculate
 (
     const scalar dt,
+    const scalar Re,
+    const scalar nu,
     const label celli,
     const scalar d,
     const scalar T,

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2014-2016 OpenFOAM Foundation
@@ -63,9 +63,9 @@ class COxidationIntrinsicRate
 :
     public SurfaceReactionModel<CloudType>
 {
-    // Private data
+    // Private Data
 
-        // Model constants
+        // Model Constants
 
             //- Stoichiometry of reaction []
             const scalar Sb_;
@@ -148,7 +148,7 @@ public:
 
 
     //- Destructor
-    virtual ~COxidationIntrinsicRate();
+    virtual ~COxidationIntrinsicRate() = default;
 
 
     // Member Functions
@@ -157,6 +157,8 @@ public:
         virtual scalar calculate
         (
             const scalar dt,
+            const scalar Re,
+            const scalar nu,
             const label celli,
             const scalar d,
             const scalar T,

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -88,20 +88,14 @@ COxidationKineticDiffusionLimitedRate
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::COxidationKineticDiffusionLimitedRate<CloudType>::
-~COxidationKineticDiffusionLimitedRate()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class CloudType>
 Foam::scalar Foam::COxidationKineticDiffusionLimitedRate<CloudType>::calculate
 (
     const scalar dt,
+    const scalar Re,
+    const scalar nu,
     const label celli,
     const scalar d,
     const scalar T,

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -135,7 +135,7 @@ public:
 
 
     //- Destructor
-    virtual ~COxidationKineticDiffusionLimitedRate();
+    virtual ~COxidationKineticDiffusionLimitedRate() = default;
 
 
     // Member Functions
@@ -144,6 +144,8 @@ public:
         virtual scalar calculate
         (
             const scalar dt,
+            const scalar Re,
+            const scalar nu,
             const label celli,
             const scalar d,
             const scalar T,

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -103,19 +103,14 @@ Foam::COxidationMurphyShaddix<CloudType>::COxidationMurphyShaddix
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::COxidationMurphyShaddix<CloudType>::~COxidationMurphyShaddix()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class CloudType>
 Foam::scalar Foam::COxidationMurphyShaddix<CloudType>::calculate
 (
     const scalar dt,
+    const scalar Re,
+    const scalar nu,
     const label celli,
     const scalar d,
     const scalar T,

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -156,7 +156,7 @@ public:
 
 
     //- Destructor
-    virtual ~COxidationMurphyShaddix();
+    virtual ~COxidationMurphyShaddix() = default;
 
 
     // Member Functions
@@ -165,6 +165,8 @@ public:
         virtual scalar calculate
         (
             const scalar dt,
+            const scalar Re,
+            const scalar nu,
             const label celli,
             const scalar d,
             const scalar T,

src/lagrangian/intermediate/Make/files

@@ -12,6 +12,7 @@ $(BASECLOUDS)/kinematicCloud/kinematicCloud.C
 $(BASECLOUDS)/thermoCloud/thermoCloud.C
 $(BASECLOUDS)/reactingCloud/reactingCloud.C
 $(BASECLOUDS)/reactingMultiphaseCloud/reactingMultiphaseCloud.C
+$(BASECLOUDS)/reactingHeterogeneousCloud/reactingHeterogeneousCloud.C
 
 
 /* kinematic parcel sub-models */
@@ -44,6 +45,12 @@ $(REACTINGMPPARCEL)/defineBasicReactingMultiphaseParcel.C
 $(REACTINGMPPARCEL)/makeBasicReactingMultiphaseParcelSubmodels.C
 
 
+/* heterogeous reacting parcel sub-models */
+REACTINGHETERMPPARCEL=$(DERIVEDPARCELS)/basicHeterogeneousReactingParcel
+$(REACTINGHETERMPPARCEL)/defineBasicHeterogeneousReactingParcel.C
+$(REACTINGHETERMPPARCEL)/makeBasicHeterogeneousReactingParcelSubmodels.C
+
+
 /* kinematic MPPIC parcel sub-models */
 KINEMATICMPPICPARCEL=$(DERIVEDPARCELS)/basicKinematicMPPICParcel
 $(KINEMATICMPPICPARCEL)/defineBasicKinematicMPPICParcel.C

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

@@ -757,23 +757,29 @@ void Foam::KinematicCloud<CloudType>::autoMap(const mapPolyMesh& mapper)
 template<class CloudType>
 void Foam::KinematicCloud<CloudType>::info()
 {
-    vector linearMomentum = linearMomentumOfSystem();
+    const vector linearMomentum =
-    reduce(linearMomentum, sumOp<vector>());
+        returnReduce(linearMomentumOfSystem(), sumOp<vector>());
 
-    scalar linearKineticEnergy = linearKineticEnergyOfSystem();
+    const scalar linearKineticEnergy =
-    reduce(linearKineticEnergy, sumOp<scalar>());
+        returnReduce(linearKineticEnergyOfSystem(), sumOp<scalar>());
+
+    const label nTotParcel = returnReduce(this->size(), sumOp<label>());
+
+    const scalar particlePerParcel =
+    (
+        nTotParcel
+      ? (returnReduce(totalParticlePerParcel(), sumOp<scalar>()) / nTotParcel)
+      : 0
+    );
 
     Info<< "Cloud: " << this->name() << nl
-        << "    Current number of parcels       = "
+        << "    Current number of parcels       = " << nTotParcel << nl
-        << returnReduce(this->size(), sumOp<label>()) << nl
         << "    Current mass in system          = "
         << returnReduce(massInSystem(), sumOp<scalar>()) << nl
-        << "    Linear momentum                 = "
+        << "    Linear momentum                 = " << linearMomentum << nl
-        << linearMomentum << nl
+        << "   |Linear momentum|                = " << mag(linearMomentum) << nl
-        << "   |Linear momentum|                = "
+        << "    Linear kinetic energy           = " << linearKineticEnergy << nl
-        << mag(linearMomentum) << nl
+        << "    Average particle per parcel     = " << particlePerParcel << nl;
-        << "    Linear kinetic energy           = "
-        << linearKineticEnergy << nl;
 
     injectors_.info(Info);
     this->surfaceFilm().info(Info);

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

@@ -502,6 +502,9 @@ public:
             //- Total linear momentum of the system
             inline vector linearMomentumOfSystem() const;
 
+            //- Average particle per parcel
+            inline scalar totalParticlePerParcel() const;
+
             //- Total linear kinetic energy in the system
             inline scalar linearKineticEnergyOfSystem() const;
 

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

@@ -297,11 +297,26 @@ Foam::KinematicCloud<CloudType>::linearMomentumOfSystem() const
 }
 
 
+template<class CloudType>
+inline Foam::scalar
+Foam::KinematicCloud<CloudType>::totalParticlePerParcel() const
+{
+    scalar parPerParcel = 0;
+
+    for (const parcelType& p : *this)
+    {
+        parPerParcel += p.nParticle();
+    }
+
+    return parPerParcel;
+}
+
+
 template<class CloudType>
 inline Foam::scalar
 Foam::KinematicCloud<CloudType>::linearKineticEnergyOfSystem() const
 {
-    scalar linearKineticEnergy = 0.0;
+    scalar linearKineticEnergy = 0;
 
     for (const parcelType& p : *this)
     {

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2017 OpenFOAM Foundation
@@ -26,7 +26,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "ReactingCloud.H"
-
 #include "CompositionModel.H"
 #include "PhaseChangeModel.H"
 
@@ -199,7 +198,6 @@ Foam::ReactingCloud<CloudType>::ReactingCloud
     cloudCopyPtr_(nullptr),
     constProps_(),
     compositionModel_(c.compositionModel_->clone()),
-//    compositionModel_(nullptr),
     phaseChangeModel_(nullptr),
     rhoTrans_(0)
 {}

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

@@ -0,0 +1,248 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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 "ReactingHeterogeneousCloud.H"
+#include "HeterogeneousReactingModel.H"
+
+// * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::ReactingHeterogeneousCloud<CloudType>::setModels()
+{
+    heterogeneousReactionModel_.reset
+    (
+        HeterogeneousReactingModel<ReactingHeterogeneousCloud<CloudType>>::New
+        (
+            this->subModelProperties(),
+            *this
+        ).ptr()
+    );
+}
+
+
+template<class CloudType>
+void Foam::ReactingHeterogeneousCloud<CloudType>::cloudReset
+(
+    ReactingHeterogeneousCloud<CloudType>& c
+)
+{
+    CloudType::cloudReset(c);
+    heterogeneousReactionModel_.reset(c.heterogeneousReactionModel_.ptr());
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::ReactingHeterogeneousCloud<CloudType>::ReactingHeterogeneousCloud
+(
+    const word& cloudName,
+    const volScalarField& rho,
+    const volVectorField& U,
+    const dimensionedVector& g,
+    const SLGThermo& thermo,
+    bool readFields
+)
+:
+    CloudType(cloudName, rho, U, g, thermo, false),
+    reactingHeterogeneousCloud(),
+    cloudCopyPtr_(nullptr),
+    heterogeneousReactionModel_(nullptr)
+{
+    if (this->solution().active())
+    {
+        setModels();
+
+        if (readFields)
+        {
+            parcelType::readFields(*this, this->composition());
+            this->deleteLostParticles();
+        }
+    }
+}
+
+
+template<class CloudType>
+Foam::ReactingHeterogeneousCloud<CloudType>::ReactingHeterogeneousCloud
+(
+    ReactingHeterogeneousCloud<CloudType>& c,
+    const word& name
+)
+:
+    CloudType(c, name),
+    reactingHeterogeneousCloud(),
+    cloudCopyPtr_(nullptr),
+    heterogeneousReactionModel_(c.heterogeneousReactionModel_->clone())
+{}
+
+
+template<class CloudType>
+Foam::ReactingHeterogeneousCloud<CloudType>::ReactingHeterogeneousCloud
+(
+    const fvMesh& mesh,
+    const word& name,
+    const ReactingHeterogeneousCloud<CloudType>& c
+)
+:
+    CloudType(mesh, name, c),
+    reactingHeterogeneousCloud(),
+    cloudCopyPtr_(nullptr),
+    heterogeneousReactionModel_(c.heterogeneousReactionModel_->clone())
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::ReactingHeterogeneousCloud<CloudType>::setParcelThermoProperties
+(
+    parcelType& parcel,
+    const scalar lagrangianDt
+)
+{
+    CloudType::setParcelThermoProperties(parcel, lagrangianDt);
+    label idS = this->composition().idSolid();
+
+    // Set inital particle composition. Overwrite thermoProperties from
+    // ReactingCloud
+    parcel.Y() = this->composition().Y0(idS);
+
+    // Set initial progress to 0
+    parcel.F().setSize(heterogeneousReactionModel_->nF(), 0.0);
+
+    // Set the parcel to combust
+    parcel.canCombust() = 1;
+}
+
+
+template<class CloudType>
+void Foam::ReactingHeterogeneousCloud<CloudType>::checkParcelProperties
+(
+    parcelType& parcel,
+    const scalar lagrangianDt,
+    const bool fullyDescribed
+)
+{
+    CloudType::checkParcelProperties(parcel, lagrangianDt, false);
+
+    const label solId = this->composition().idSolid();
+    const label liqId = this->composition().idLiquid();
+    const label gasId = this->composition().idGas();
+
+    // Check YMixture is pure solid
+    if
+    (
+        this->composition().YMixture0()[solId] != 1.0
+     || this->composition().YMixture0()[liqId] != 0.0
+     || this->composition().YMixture0()[gasId] != 0.0
+    )
+    {
+        FatalErrorInFunction
+            << "The supplied composition must be : " << nl
+            << "    YGasTot0        0 : " << nl
+            << "    YLiquidTot0     0 : " << nl
+            << "    YSolidTot0      1 : " << nl
+            << "This Cloud only works with pure solid particles."
+            << abort(FatalError);
+    }
+    if (this->composition().liquids().size() > 0)
+    {
+        FatalErrorInFunction
+            << "The supplied composition has a liquid phase. " << nl
+            << "This Cloud only works with pure solid particles."
+            << abort(FatalError);
+    }
+}
+
+
+template<class CloudType>
+void Foam::ReactingHeterogeneousCloud<CloudType>::storeState()
+{
+    cloudCopyPtr_.reset
+    (
+        static_cast<ReactingHeterogeneousCloud<CloudType>*>
+        (
+            clone(this->name() + "Copy").ptr()
+        )
+    );
+}
+
+
+template<class CloudType>
+void Foam::ReactingHeterogeneousCloud<CloudType>::restoreState()
+{
+    cloudReset(cloudCopyPtr_());
+    cloudCopyPtr_.clear();
+}
+
+
+template<class CloudType>
+void Foam::ReactingHeterogeneousCloud<CloudType>::evolve()
+{
+    if (this->solution().canEvolve())
+    {
+        typename parcelType::trackingData td(*this);
+
+        this->solve(*this, td);
+    }
+}
+
+
+template<class CloudType>
+void Foam::ReactingHeterogeneousCloud<CloudType>::autoMap
+(
+    const mapPolyMesh& mapper
+)
+{
+    Cloud<parcelType>::autoMap(mapper);
+
+    this->updateMesh();
+}
+
+
+template<class CloudType>
+void Foam::ReactingHeterogeneousCloud<CloudType>::info()
+{
+    CloudType::info();
+    heterogeneousReactionModel_->info(Info);
+}
+
+
+template<class CloudType>
+void Foam::ReactingHeterogeneousCloud<CloudType>::writeFields() const
+{
+    CloudType::writeFields();
+}
+
+
+template<class CloudType>
+void Foam::ReactingHeterogeneousCloud<CloudType>::
+writeObjects(objectRegistry& obr) const
+{
+    CloudType::particleType::writeObjects(*this, this->composition(), obr);
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,272 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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::ReactingHeterogeneousCloud
+
+Group
+    grpLagrangianIntermediateClouds
+
+Description
+    Templated base class for reacting heterogeneous cloud
+
+    - Adds to reacting cloud:
+      - Heterogeneous reaction model
+
+SourceFiles
+    ReactingHeterogeneousCloudI.H
+    ReactingHeterogeneousCloud.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ReactingHeterogeneousCloud_H
+#define ReactingHeterogeneousCloud_H
+
+#include "reactingHeterogeneousCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+template<class CloudType>
+class HeterogeneousReactingModel;
+
+
+/*---------------------------------------------------------------------------*\
+                      Class ReactingHeterogeneousCloud Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class ReactingHeterogeneousCloud
+:
+    public CloudType,
+    public reactingHeterogeneousCloud
+{
+public:
+
+    // Public typedefs
+
+        //- Type of cloud this cloud was instantiated for
+        typedef CloudType cloudType;
+
+        //- Type of parcel the cloud was instantiated for
+        typedef typename CloudType::particleType parcelType;
+
+        //- Convenience typedef for this cloud type
+        typedef ReactingHeterogeneousCloud<CloudType>
+            reactingHeterogeneousCloudType;
+
+
+private:
+
+    // Private Data
+
+        //- Cloud copy pointer
+        autoPtr<ReactingHeterogeneousCloud<CloudType>> cloudCopyPtr_;
+
+
+    // Private Member Functions
+
+        //- No copy construct
+        ReactingHeterogeneousCloud(const ReactingHeterogeneousCloud&) = delete;
+
+        //- No copy assignment
+        void operator=(const ReactingHeterogeneousCloud&) = delete;
+
+
+protected:
+
+    // Protected Data
+
+        // References to the cloud sub-models
+
+            //- Heterogeneous reaction model
+            autoPtr
+            <
+                HeterogeneousReactingModel
+                <
+                    ReactingHeterogeneousCloud<CloudType>
+                >
+            > heterogeneousReactionModel_;
+
+
+    // Protected Member Functions
+
+        // Initialisation
+
+            //- Set cloud sub-models
+            void setModels();
+
+
+        // Cloud evolution functions
+
+            //- Reset state of cloud
+            void cloudReset(ReactingHeterogeneousCloud<CloudType>& c);
+
+
+public:
+
+    // Constructors
+
+        //- Construct given carrier gas fields
+        ReactingHeterogeneousCloud
+        (
+            const word& cloudName,
+            const volScalarField& rho,
+            const volVectorField& U,
+            const dimensionedVector& g,
+            const SLGThermo& thermo,
+            bool readFields = true
+        );
+
+        //- Copy constructor with new name
+        ReactingHeterogeneousCloud
+        (
+            ReactingHeterogeneousCloud<CloudType>& c, const word& name
+        );
+
+        //- Copy constructor with new name - creates bare cloud
+        ReactingHeterogeneousCloud
+        (
+            const fvMesh& mesh,
+            const word& name,
+            const ReactingHeterogeneousCloud<CloudType>& c
+        );
+
+        //- Construct and return clone based on (this) with new name
+        virtual autoPtr<Cloud<parcelType>> clone(const word& name)
+        {
+            return autoPtr<Cloud<parcelType>>
+            (
+                new ReactingHeterogeneousCloud(*this, name)
+            );
+        }
+
+        //- Construct and return bare clone based on (this) with new name
+        virtual autoPtr<Cloud<parcelType>> cloneBare(const word& name) const
+        {
+            return autoPtr<Cloud<parcelType>>
+            (
+                new ReactingHeterogeneousCloud(this->mesh(), name, *this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~ReactingHeterogeneousCloud() = default;
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return a reference to the cloud copy
+            inline const ReactingHeterogeneousCloud& cloudCopy() const;
+
+            //- Return progress variable dimension
+            inline label nF() const;
+
+
+        // Sub-models
+
+            //- Return reference to model
+            inline const HeterogeneousReactingModel
+            <
+                ReactingHeterogeneousCloud<CloudType>
+            >& heterogeneousReaction() const;
+
+
+            inline HeterogeneousReactingModel
+            <
+                ReactingHeterogeneousCloud<CloudType>
+            >& heterogeneousReaction();
+
+
+        // Cloud evolution functions
+
+            //- Set parcel thermo properties
+            void setParcelThermoProperties
+            (
+                parcelType& parcel,
+                const scalar lagrangianDt
+            );
+
+            //- Check parcel properties
+            void checkParcelProperties
+            (
+                parcelType& parcel,
+                const scalar lagrangianDt,
+                const bool fullyDescribed
+            );
+
+            //- Store the current cloud state
+            void storeState();
+
+            //- Reset the current cloud to the previously stored state
+            void restoreState();
+
+            //- Evolve the cloud
+            void evolve();
+
+
+        // Mapping
+
+            //- Remap the cells of particles corresponding to the
+            //  mesh topology change with a default tracking data object
+            virtual void autoMap(const mapPolyMesh&);
+
+
+        // I-O
+
+            //- Print cloud information
+            void info();
+
+            //- Write the field data for the cloud
+            virtual void writeFields() const;
+
+            //- Write particle fields as objects into the obr registry
+            virtual void writeObjects(objectRegistry& obr) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "ReactingHeterogeneousCloudI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "ReactingHeterogeneousCloud.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/clouds/Templates/ReactingHeterogeneousCloud/ReactingHeterogeneousCloudI.H

@@ -0,0 +1,60 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+inline const Foam::ReactingHeterogeneousCloud<CloudType>&
+Foam::ReactingHeterogeneousCloud<CloudType>::cloudCopy() const
+{
+    return *cloudCopyPtr_;
+}
+
+
+template<class CloudType>
+inline const
+Foam::HeterogeneousReactingModel<Foam::ReactingHeterogeneousCloud<CloudType>>&
+Foam::ReactingHeterogeneousCloud<CloudType>::heterogeneousReaction() const
+{
+    return *heterogeneousReactionModel_;
+}
+
+
+template<class CloudType>
+inline
+Foam::HeterogeneousReactingModel<Foam::ReactingHeterogeneousCloud<CloudType>>&
+Foam::ReactingHeterogeneousCloud<CloudType>::heterogeneousReaction()
+{
+    return *heterogeneousReactionModel_;
+}
+
+
+template<class CloudType>
+inline Foam::label Foam::ReactingHeterogeneousCloud<CloudType>::nF() const
+{
+    return *heterogeneousReactionModel_.nF();
+}
+
+// ************************************************************************* //

src/lagrangian/intermediate/clouds/baseClasses/reactingHeterogeneousCloud/reactingHeterogeneousCloud.C

@@ -0,0 +1,36 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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 "reactingHeterogeneousCloud.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(reactingHeterogeneousCloud, 0);
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/clouds/baseClasses/reactingHeterogeneousCloud/reactingHeterogeneousCloud.H

@@ -0,0 +1,76 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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::reactingHeterogeneousCloud
+
+Description
+    Virtual abstract base class for templated ReactingCloud
+
+SourceFiles
+    reactingHeterogeneousCloud.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef reactingHeterogeneousCloud_H
+#define reactingHeterogeneousCloud_H
+
+#include "typeInfo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                 Class reactingHeterogeneousCloud Declaration
+\*---------------------------------------------------------------------------*/
+
+class reactingHeterogeneousCloud
+{
+public:
+
+    //- Runtime type information
+    TypeName("reactingHeterogeneousCloud");
+
+
+    // Constructors
+
+        //- Null constructor
+        reactingHeterogeneousCloud() = default;
+
+
+    //- Destructor
+    virtual ~reactingHeterogeneousCloud() = default;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,68 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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::basicHeterogeneousReactingCloud
+
+Description
+    Cloud class to introduce heterogeneou reacting parcels
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef basicHeterogeneousReactingCloud_H
+#define basicHeterogeneousReactingCloud_H
+
+#include "Cloud.H"
+#include "KinematicCloud.H"
+#include "ThermoCloud.H"
+#include "ReactingCloud.H"
+#include "ReactingHeterogeneousCloud.H"
+#include "basicHeterogeneousReactingParcel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef ReactingHeterogeneousCloud
+    <
+        ReactingCloud
+        <
+            ThermoCloud
+            <
+                KinematicCloud
+                <
+                    Cloud
+                    <
+                        basicHeterogeneousReactingParcel
+                    >
+                >
+            >
+        >
+    > basicHeterogeneousReactingCloud;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,397 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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 "ReactingHeterogeneousParcel.H"
+#include "mathematicalConstants.H"
+
+using namespace Foam::constant::mathematical;
+
+
+// * * * * * * * * * * * *  Private Member Functions * * * * * * * * * * * * //
+
+template<class ParcelType>
+template<class TrackCloudType>
+Foam::scalar Foam::ReactingHeterogeneousParcel<ParcelType>::CpEff
+(
+    TrackCloudType& cloud,
+    trackingData& td,
+    const scalar p,
+    const scalar T,
+    const label idS
+) const
+{
+    return cloud.composition().Cp(idS, this->Y_, p, T);
+}
+
+
+template<class ParcelType>
+template<class TrackCloudType>
+Foam::scalar Foam::ReactingHeterogeneousParcel<ParcelType>::HsEff
+(
+    TrackCloudType& cloud,
+    trackingData& td,
+    const scalar p,
+    const scalar T,
+    const label idS
+) const
+{
+    return cloud.composition().Hs(idS, this->Y_, p, T);
+}
+
+
+template<class ParcelType>
+template<class TrackCloudType>
+Foam::scalar Foam::ReactingHeterogeneousParcel<ParcelType>::LEff
+(
+    TrackCloudType& cloud,
+    trackingData& td,
+    const scalar p,
+    const scalar T,
+    const label idS
+) const
+{
+    return cloud.composition().L(idS, this->Y_, p, T);
+}
+
+
+// * * * * * * * * * * *  Protected Member Functions * * * * * * * * * * * * //
+
+template<class ParcelType>
+template<class TrackCloudType>
+void Foam::ReactingHeterogeneousParcel<ParcelType>::calc
+(
+    TrackCloudType& cloud,
+    trackingData& td,
+    const scalar dt
+)
+{
+    typedef typename TrackCloudType::reactingCloudType reactingCloudType;
+
+    const CompositionModel<reactingCloudType>& composition =
+        cloud.composition();
+
+    // Define local properties at beginning of timestep
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    const scalar np0 = this->nParticle_;
+    const scalar d0 = this->d_;
+    const vector& U0 = this->U_;
+    const scalar T0 = this->T_;
+    const scalar mass0 = this->mass();
+
+    const scalar pc = td.pc();
+
+    const label idS = composition.idSolid();
+
+    // Calc surface values
+    scalar Ts, rhos, mus, Prs, kappas;
+    this->calcSurfaceValues(cloud, td, T0, Ts, rhos, mus, Prs, kappas);
+    scalar Res = this->Re(rhos, U0, td.Uc(), d0, mus);
+
+    // Sources
+    //~~~~~~~~
+
+    // Explicit momentum source for particle
+    vector Su = Zero;
+
+    // Linearised momentum source coefficient
+    scalar Spu = 0.0;
+
+    // Momentum transfer from the particle to the carrier phase
+    vector dUTrans = Zero;
+
+    // Explicit enthalpy source for particle
+    scalar Sh = 0.0;
+
+    // Linearised enthalpy source coefficient
+    scalar Sph = 0.0;
+
+    // Sensible enthalpy transfer from the particle to the carrier phase
+    scalar dhsTrans = 0.0;
+
+    // Molar flux of species emitted from the particle (kmol/m^2/s)
+    scalar Ne = 0.0;
+
+    // Surface concentrations of emitted species
+    scalarField Cs(composition.carrier().species().size(), 0.0);
+
+    // Sum Ni*Cpi*Wi of emission species
+    scalar NCpW = 0.0;
+
+
+    // Heterogeneous reactions
+    // ~~~~~~~~~~~~~~~~~
+
+    // Change in carrier phase composition due to surface reactions
+    scalarField dMassSRSolid(this->Y_.size(), 0.0);
+    scalarField dMassSRCarrier(composition.carrier().species().size(), 0.0);
+
+    // Calc mass and enthalpy transfer due to reactions
+    calcHeterogeneousReactions
+    (
+        cloud,
+        td,
+        dt,
+        Res,
+        mus/rhos,
+        d0,
+        T0,
+        mass0,
+        canCombust_,
+        Ne,
+        NCpW,
+        this->Y_,
+        F_,
+        dMassSRSolid,
+        dMassSRCarrier,
+        Sh,
+        dhsTrans
+    );
+
+    // 2. Update the parcel properties due to change in mass
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    scalarField dMassSolid(dMassSRSolid);
+
+    scalar mass1 = mass0 - sum(dMassSolid);
+
+    // Remove the particle when mass falls below minimum threshold
+    if (np0*mass1 < cloud.constProps().minParcelMass())
+    {
+        td.keepParticle = false;
+
+        return;
+    }
+
+    // Only solid is used to update mass fractions
+    (void)this->updateMassFraction(mass0, dMassSolid, this->Y_);
+
+
+    // Update particle density or diameter
+    if (cloud.constProps().constantVolume())
+    {
+        this->rho_ = mass1/this->volume();
+    }
+    else
+    {
+        this->d_ = cbrt(mass1/this->rho_*6.0/pi);
+    }
+
+    // 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);
+
+
+    // 3. Compute heat- and momentum transfers
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Heat transfer
+    // ~~~~~~~~~~~~~
+
+    // Calculate new particle temperature
+    this->T_ =
+        this->calcHeatTransfer
+        (
+            cloud,
+            td,
+            dt,
+            Res,
+            Prs,
+            kappas,
+            NCpW,
+            Sh,
+            dhsTrans,
+            Sph
+        );
+
+    //DebugVar(np0);
+
+    this->Cp_ = CpEff(cloud, td, pc, this->T_, idS);
+
+    // Motion
+    // ~~~~~~
+
+    // Calculate new particle velocity
+    this->U_ =
+        this->calcVelocity(cloud, td, dt, Res, mus, mass1, Su, dUTrans, Spu);
+
+
+    // 4. Accumulate carrier phase source terms
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    if (cloud.solution().coupled())
+    {
+        // No mapping between solid components and carrier phase
+        /*
+        forAll(this->Y_, i)
+        {
+            scalar dm = np0*dMassSolid[i];
+            label gid = composition.localToCarrierId(SLD, i);
+            scalar hs = composition.carrier().Hs(gid, pc, T0);
+            cloud.rhoTrans(gid)[this->cell()] += dm;
+            cloud.UTrans()[this->cell()] += dm*U0;
+            cloud.hsTrans()[this->cell()] += dm*hs;
+        }
+        */
+
+        forAll(dMassSRCarrier, i)
+        {
+            scalar dm = np0*dMassSRCarrier[i];
+            scalar hs = composition.carrier().Hs(i, pc, T0);
+            cloud.rhoTrans(i)[this->cell()] += dm;
+            cloud.UTrans()[this->cell()] += dm*U0;
+            cloud.hsTrans()[this->cell()] += dm*hs;
+        }
+
+        // Update momentum transfer
+        cloud.UTrans()[this->cell()] += np0*dUTrans;
+        cloud.UCoeff()[this->cell()] += np0*Spu;
+
+        // Update sensible enthalpy transfer
+        cloud.hsTrans()[this->cell()] += np0*dhsTrans;
+        cloud.hsCoeff()[this->cell()] += np0*Sph;
+
+        // Update radiation fields
+        if (cloud.radiation())
+        {
+            const scalar ap = this->areaP();
+            const scalar T4 = pow4(T0);
+            cloud.radAreaP()[this->cell()] += dt*np0*ap;
+            cloud.radT4()[this->cell()] += dt*np0*T4;
+            cloud.radAreaPT4()[this->cell()] += dt*np0*ap*T4;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+template<class ParcelType>
+template<class TrackCloudType>
+void Foam::ReactingHeterogeneousParcel<ParcelType>::calcHeterogeneousReactions
+(
+    TrackCloudType& cloud,
+    trackingData& td,
+    const scalar dt,
+    const scalar Re,
+    const scalar nu,
+    const scalar d,
+    const scalar T,
+    const scalar mass,
+    const label canCombust,
+    const scalar Ne,
+    scalar& NCpW,
+    const scalarField& YSolid,
+    scalarField& F,
+    scalarField& dMassSRSolid,
+    scalarField& dMassSRCarrier,
+    scalar& Sh,
+    scalar& dhsTrans
+) const
+{
+    // Check that model is active
+    if (!cloud.heterogeneousReaction().active())
+    {
+        return;
+    }
+
+    // Initialise demand-driven constants
+    (void)cloud.constProps().hRetentionCoeff();
+    (void)cloud.constProps().TMax();
+
+    // Check that model is active
+    if (canCombust != 1)
+    {
+        return;
+    }
+
+    // Update reactions
+    const scalar hReaction = cloud.heterogeneousReaction().calculate
+    (
+        dt,
+        Re,
+        nu,
+        this->cell(),
+        d,
+        T,
+        td.Tc(),
+        td.pc(),
+        td.rhoc(),
+        mass,
+        YSolid,
+        F,
+        Ne,
+        NCpW,
+        dMassSRSolid,
+        dMassSRCarrier
+    );
+
+    cloud.heterogeneousReaction().addToSurfaceReactionMass
+    (
+        this->nParticle_*sum(dMassSRSolid)
+    );
+
+    const scalar xsi = min(T/cloud.constProps().TMax(), 1.0);
+    const scalar coeff =
+        (1.0 - xsi*xsi)*cloud.constProps().hRetentionCoeff();
+
+    Sh += coeff*hReaction/dt;
+
+    dhsTrans += (1.0 - coeff)*hReaction;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class ParcelType>
+Foam::ReactingHeterogeneousParcel<ParcelType>::ReactingHeterogeneousParcel
+(
+    const ReactingHeterogeneousParcel<ParcelType>& p
+)
+:
+    ParcelType(p),
+    F_(p.F_),
+    canCombust_(p.canCombust_)
+{}
+
+
+template<class ParcelType>
+Foam::ReactingHeterogeneousParcel<ParcelType>::ReactingHeterogeneousParcel
+(
+    const ReactingHeterogeneousParcel<ParcelType>& p,
+    const polyMesh& mesh
+)
+:
+    ParcelType(p, mesh),
+    F_(p.F_),
+    canCombust_(p.canCombust_)
+{}
+
+
+// * * * * * * * * * * * * * * IOStream operators  * * * * * * * * * * * * * //
+
+#include "ReactingHeterogeneousParcelIO.C"
+
+// ************************************************************************* //

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

@@ -0,0 +1,426 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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::ReactingHeterogeneousParcel
+
+Group
+    grpLagrangianIntermediateParcels
+
+Description
+    Reacting heterogeneous Parcel
+
+SourceFiles
+    ReactingHeterogeneousParcelI.H
+    ReactingHeterogeneousParcel.C
+    ReactingHeterogeneousParcelIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ReactingHeterogeneousParcel_H
+#define ReactingHeterogeneousParcel_H
+
+#include "demandDrivenEntry.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<class ParcelType>
+class ReactingHeterogeneousParcel;
+
+template<class ParcelType>
+Ostream& operator<<
+(
+    Ostream&,
+    const ReactingHeterogeneousParcel<ParcelType>&
+);
+
+/*---------------------------------------------------------------------------*\
+                 Class ReactingHeterogeneousParcel Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class ParcelType>
+class ReactingHeterogeneousParcel
+:
+    public ParcelType
+{
+public:
+
+    //- Size in bytes of the fields
+    static const std::size_t sizeofFields;
+
+    //- Class to hold reacting particle constant properties
+    class constantProperties
+    :
+        public ParcelType::constantProperties
+    {
+        // Private data
+
+            //- Fraction of enthalpy retained by parcel due to surface
+            //  reactions
+            demandDrivenEntry<scalar> hRetentionCoeff_;
+
+    public:
+
+        // Constructors
+
+            //- Null constructor
+            constantProperties();
+
+            //- Copy constructor
+            constantProperties(const constantProperties& cp);
+
+            //- Construct from dictionary
+            constantProperties(const dictionary& parentDict);
+
+
+        // Access
+
+            //- Return const access to the fraction of enthalpy retained by
+            //  parcel due to surface reactions
+            inline scalar hRetentionCoeff() const;
+    };
+
+
+    //- Use base tracking data
+    typedef typename ParcelType::trackingData trackingData;
+
+
+private:
+
+    // Private Member Functions
+
+        //- Return the mixture effective specific heat capacity
+        template<class TrackCloudType>
+        scalar CpEff
+        (
+            TrackCloudType& cloud,
+            trackingData& td,
+            const scalar p,
+            const scalar T,
+            const label idS
+        ) const;
+
+        //- Return the mixture effective sensible enthalpy
+        template<class TrackCloudType>
+        scalar HsEff
+        (
+            TrackCloudType& cloud,
+            trackingData& td,
+            const scalar p,
+            const scalar T,
+            const label idS
+        ) const;
+
+        //- Return the mixture effective latent heat
+        template<class TrackCloudType>
+        scalar LEff
+        (
+            TrackCloudType& cloud,
+            trackingData& td,
+            const scalar p,
+            const scalar T,
+            const label idS
+        ) const;
+
+
+protected:
+
+    // Protected data
+
+        // Parcel properties
+
+            //- Progress variables for reactions
+            scalarField F_;
+
+            //- Flag to identify if the particle can devolatilise and combust
+            //  Combustion possible only after volatile content falls below
+            //  threshold value.  States include:
+            //  0 = can combust but can change
+            //  1 = can devolatilise, can combust
+            // -1 = cannot devolatilise or combust, and cannot change
+            label canCombust_;
+
+
+    // Protected Member Functions
+
+
+        //- Calculate surface reactions
+        template<class TrackCloudType>
+        void calcHeterogeneousReactions
+        (
+            TrackCloudType& cloud,
+            trackingData& td,
+            const scalar dt,           // timestep
+            const scalar Res,          // Re
+            const scalar nu,           // nu
+            const scalar d,            // diameter
+            const scalar T,            // temperature
+            const scalar mass,         // mass
+            const label canCombust,    // 'can combust' flag
+            const scalar N,            // flux of species emitted from particle
+            scalar& NCpW,
+            const scalarField& YSolid, // solid-phase mass fractions
+            scalarField& F,            // progress of each reaction
+            scalarField& dMassSRSolid, // solid-phase mass transfer - local
+            scalarField& dMassSRCarrier, // carrier phase mass transfer
+            scalar& Sh,                // explicit particle enthalpy source
+            scalar& dhsTrans           // sensible enthalpy transfer to carrier
+        ) const;
+
+
+public:
+
+    // Static data members
+
+        //- Runtime type information
+        TypeName("ReactingHeterogeneousParcel");
+
+        //- String representation of properties
+        AddToPropertyList
+        (
+            ParcelType,
+          + " nReactions(F1..FN)"
+        );
+
+        //- String representation of property types
+        AddToPropertyTypes
+        (
+            ParcelType,
+            " scalars "
+        );
+
+
+    // Constructors
+
+        //- Construct from mesh, position and topology
+        //  Other properties initialised as null
+        inline ReactingHeterogeneousParcel
+        (
+            const polyMesh& mesh,
+            const barycentric& coordinates,
+            const label celli,
+            const label tetFacei,
+            const label tetPti
+        );
+
+        //- Construct from a position and a cell, searching for the rest of the
+        //  required topology. Other properties are initialised as null.
+        inline ReactingHeterogeneousParcel
+        (
+            const polyMesh& mesh,
+            const vector& position,
+            const label celli
+        );
+
+        //- Construct from components
+        inline ReactingHeterogeneousParcel
+        (
+            const polyMesh& mesh,
+            const barycentric& coordinates,
+            const label celli,
+            const label tetFacei,
+            const label tetPti,
+            const label typeId,
+            const scalar nParticle0,
+            const scalar d0,
+            const scalar dTarget0,
+            const vector& U0,
+            const vector& f0,
+            const vector& angularMomentum0,
+            const vector& torque0,
+            const scalarField& Y,
+            const scalarField& F,
+            const constantProperties& constProps
+        );
+
+        //- Construct from Istream
+        ReactingHeterogeneousParcel
+        (
+            const polyMesh& mesh,
+            Istream& is,
+            bool readFields = true,
+            bool newFormat = true
+        );
+
+        //- Construct as a copy
+        ReactingHeterogeneousParcel(const ReactingHeterogeneousParcel& p);
+
+        //- Construct as a copy
+        ReactingHeterogeneousParcel
+        (
+            const ReactingHeterogeneousParcel& p,
+            const polyMesh& mesh
+        );
+
+        //- Construct and return a (basic particle) clone
+        virtual autoPtr<particle> clone() const
+        {
+            return autoPtr<particle>(new ReactingHeterogeneousParcel(*this));
+        }
+
+        //- Construct and return a (basic particle) clone
+        virtual autoPtr<particle> clone(const polyMesh& mesh) const
+        {
+            return autoPtr<particle>
+            (
+                new ReactingHeterogeneousParcel(*this, mesh)
+            );
+        }
+
+        //- Factory class to read-construct particles used for
+        //  parallel transfer
+        class iNew
+        {
+            const polyMesh& mesh_;
+
+        public:
+
+            iNew(const polyMesh& mesh)
+            :
+                mesh_(mesh)
+            {}
+
+            autoPtr<ReactingHeterogeneousParcel<ParcelType>> operator()
+            (
+                Istream& is
+            ) const
+            {
+                return autoPtr<ReactingHeterogeneousParcel<ParcelType>>
+                (
+                    new ReactingHeterogeneousParcel<ParcelType>
+                        (mesh_, is, true)
+                );
+            }
+        };
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return const access to F
+            inline const scalarField& F() const;
+
+            //- Return const access to the canCombust flag
+            inline label canCombust() const;
+
+
+        // Edit
+
+             //- Return access to F
+            inline scalarField& F();
+
+            //- Return access to the canCombust flag
+            inline label& canCombust();
+
+
+        // Main calculation loop
+
+
+            //- Update parcel properties over the time interval
+            template<class TrackCloudType>
+            void calc
+            (
+                TrackCloudType& cloud,
+                trackingData& td,
+                const scalar dt
+            );
+
+
+        // I-O
+
+            //- Read - composition supplied
+            template<class CloudType, class CompositionType>
+            static void readFields
+            (
+                CloudType& c,
+                const CompositionType& compModel
+            );
+
+            //- Read - no composition
+            template<class CloudType>
+            static void readFields(CloudType& c);
+
+            //- Write - composition supplied
+            template<class CloudType, class CompositionType>
+            static void writeFields
+            (
+                const CloudType& c,
+                const CompositionType& compModel
+            );
+
+            //- Read - no composition
+            template<class CloudType>
+            static void writeFields(const CloudType& c);
+
+            //- Write particle fields as objects into the obr registry
+            //  - no composition
+            template<class CloudType>
+            static void writeObjects
+            (
+                const CloudType& c,
+                objectRegistry& obr
+            );
+
+            //- Write particle fields as objects into the obr registry
+            template<class CloudType, class CompositionType>
+            static void writeObjects
+            (
+                const CloudType& c,
+                const CompositionType& compModel,
+                objectRegistry& obr
+            );
+
+
+    // Ostream Operator
+
+        friend Ostream& operator<< <ParcelType>
+        (
+            Ostream&,
+            const ReactingHeterogeneousParcel<ParcelType>&
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "ReactingHeterogeneousParcelI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "ReactingHeterogeneousParcel.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/parcels/Templates/ReactingHeterogeneousParcel/ReactingHeterogeneousParcelI.H

@@ -0,0 +1,189 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class ParcelType>
+inline Foam::ReactingHeterogeneousParcel<ParcelType>::constantProperties::
+constantProperties()
+:
+    ParcelType::constantProperties(),
+    hRetentionCoeff_(this->dict_, 0.0)
+{}
+
+
+template<class ParcelType>
+inline Foam::ReactingHeterogeneousParcel<ParcelType>::constantProperties::
+constantProperties
+(
+    const constantProperties& cp
+)
+:
+    ParcelType::constantProperties(cp),
+    hRetentionCoeff_(cp.hRetentionCoeff_)
+{}
+
+
+template<class ParcelType>
+inline Foam::ReactingHeterogeneousParcel<ParcelType>::constantProperties::
+constantProperties
+(
+    const dictionary& parentDict
+)
+:
+    ParcelType::constantProperties(parentDict),
+    hRetentionCoeff_(this->dict_, "hRetentionCoeff")
+{}
+
+
+template<class ParcelType>
+inline Foam::ReactingHeterogeneousParcel<ParcelType>::ReactingHeterogeneousParcel
+(
+    const polyMesh& mesh,
+    const barycentric& coordinates,
+    const label celli,
+    const label tetFacei,
+    const label tetPti
+)
+:
+    ParcelType(mesh, coordinates, celli, tetFacei, tetPti),
+    F_(0),
+    canCombust_(0)
+{}
+
+
+template<class ParcelType>
+inline Foam::ReactingHeterogeneousParcel<ParcelType>::ReactingHeterogeneousParcel
+(
+    const polyMesh& mesh,
+    const vector& position,
+    const label celli
+)
+:
+    ParcelType(mesh, position, celli),
+    F_(0),
+    canCombust_(0)
+{}
+
+
+template<class ParcelType>
+inline Foam::ReactingHeterogeneousParcel<ParcelType>::ReactingHeterogeneousParcel
+(
+    const polyMesh& mesh,
+    const barycentric& coordinates,
+    const label celli,
+    const label tetFacei,
+    const label tetPti,
+    const label typeId,
+    const scalar nParticle0,
+    const scalar d0,
+    const scalar dTarget0,
+    const vector& U0,
+    const vector& f0,
+    const vector& angularMomentum0,
+    const vector& torque0,
+    const scalarField& Y,
+    const scalarField& F,
+    const constantProperties& constProps
+)
+:
+    ParcelType
+    (
+        mesh,
+        coordinates,
+        celli,
+        tetFacei,
+        tetPti,
+        typeId,
+        nParticle0,
+        d0,
+        dTarget0,
+        U0,
+        f0,
+        angularMomentum0,
+        torque0,
+        Y,
+        constProps
+    ),
+    F_(F),
+    canCombust_(0)
+{}
+
+
+// * * * * * * * * * constantProperties Member Functions * * * * * * * * * * //
+
+
+template<class ParcelType>
+inline Foam::scalar
+Foam::ReactingHeterogeneousParcel<ParcelType>::constantProperties::
+hRetentionCoeff() const
+{
+    scalar value = hRetentionCoeff_.value();
+
+    if ((value < 0) || (value > 1))
+    {
+        FatalErrorInFunction
+            << "hRetentionCoeff must be in the range 0 to 1" << nl
+            << exit(FatalError) << endl;
+    }
+
+    return value;
+}
+
+
+// * * * * * * * * * * ThermoParcel Member Functions * * * * * * * * * * * * //
+
+template<class ParcelType>
+inline const Foam::scalarField& Foam::ReactingHeterogeneousParcel<ParcelType>::
+F() const
+{
+    return F_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalarField& Foam::ReactingHeterogeneousParcel<ParcelType>::
+F()
+{
+    return F_;
+}
+
+
+template<class ParcelType>
+inline Foam::label
+Foam::ReactingHeterogeneousParcel<ParcelType>::canCombust() const
+{
+    return canCombust_;
+}
+
+
+template<class ParcelType>
+inline Foam::label& Foam::ReactingHeterogeneousParcel<ParcelType>::canCombust()
+{
+    return canCombust_;
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/parcels/Templates/ReactingHeterogeneousParcel/ReactingHeterogeneousParcelIO.C

@@ -0,0 +1,354 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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 "ReactingHeterogeneousParcel.H"
+#include "IOstreams.H"
+#include "HeterogeneousReactingModel.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+template<class ParcelType>
+Foam::string Foam::ReactingHeterogeneousParcel<ParcelType>::propertyList_ =
+    Foam::ReactingHeterogeneousParcel<ParcelType>::propertyList();
+
+template<class ParcelType>
+Foam::string Foam::ReactingHeterogeneousParcel<ParcelType>::propertyTypes_ =
+    Foam::ReactingHeterogeneousParcel<ParcelType>::propertyTypes();
+
+template<class ParcelType>
+const std::size_t Foam::ReactingHeterogeneousParcel<ParcelType>::sizeofFields
+(
+    sizeof(scalar)
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class ParcelType>
+Foam::ReactingHeterogeneousParcel<ParcelType>::ReactingHeterogeneousParcel
+(
+    const polyMesh& mesh,
+    Istream& is,
+    bool readFields,
+    bool newFormat
+)
+:
+    ParcelType(mesh, is, readFields, newFormat),
+    F_(0),
+    canCombust_(1)
+{
+    if (readFields)
+    {
+        DynamicList<scalar> F;
+
+        is >> F;
+
+        F_.transfer(F);
+    }
+
+    is.check(FUNCTION_NAME);
+}
+
+
+template<class ParcelType>
+template<class CloudType>
+void Foam::ReactingHeterogeneousParcel<ParcelType>::readFields(CloudType& c)
+{
+    ParcelType::readFields(c);
+}
+
+
+template<class ParcelType>
+template<class CloudType, class CompositionType>
+void Foam::ReactingHeterogeneousParcel<ParcelType>::readFields
+(
+    CloudType& c,
+    const CompositionType& compModel
+)
+{
+    bool valid = c.size();
+
+    ParcelType::readFields(c);
+
+    IOField<scalar> mass0
+    (
+        c.fieldIOobject("mass0", IOobject::MUST_READ),
+        valid
+    );
+    c.checkFieldIOobject(c, mass0);
+
+    label i = 0;
+    for (ReactingHeterogeneousParcel<ParcelType>& p : c)
+    {
+        p.mass0() = mass0[i];
+        ++i;
+    }
+
+    const label idSolid = compModel.idSolid();
+    const wordList& solidNames = compModel.componentNames(idSolid);
+
+    // WIP until find a way to get info from Reacting Heterogeneous model
+    label nF(1);
+
+    // Set storage for each Y... for each parcel
+    for (ReactingHeterogeneousParcel<ParcelType>& p : c)
+    {
+        p.Y().setSize(solidNames.size(), Zero);
+        p.F_.setSize(nF, Zero);
+    }
+
+    for (label i = 0; i < nF; i++)
+    {
+        // Read F
+        IOField<scalar> F
+        (
+            c.fieldIOobject
+            (
+                "F" + name(i),
+                IOobject::MUST_READ
+            ),
+            valid
+        );
+
+        label j = 0;
+        for (ReactingHeterogeneousParcel<ParcelType>& p : c)
+        {
+            p.F_[i] = F[j];
+            ++j;
+        }
+    }
+
+
+    forAll(solidNames, j)
+    {
+        IOField<scalar> Y
+        (
+            c.fieldIOobject
+            (
+                "Y" + solidNames[j],
+                IOobject::MUST_READ
+            ),
+            valid
+        );
+
+        label i = 0;
+        for (ReactingHeterogeneousParcel<ParcelType>& p : c)
+        {
+            p.Y()[j] = Y[i];
+            ++i;
+        }
+    }
+}
+
+
+template<class ParcelType>
+template<class CloudType>
+void Foam::ReactingHeterogeneousParcel<ParcelType>::writeFields
+(
+    const CloudType& c
+)
+{
+    ParcelType::writeFields(c);
+}
+
+
+template<class ParcelType>
+template<class CloudType, class CompositionType>
+void Foam::ReactingHeterogeneousParcel<ParcelType>::writeFields
+(
+    const CloudType& c,
+    const CompositionType& compModel
+)
+{
+    // Skip Reacting layer. This class takes charge of
+    // writing Ysolids and F
+    ThermoParcel<KinematicParcel<particle>>::writeFields(c);
+
+    label np = c.size();
+
+    IOField<scalar> mass0(c.fieldIOobject("mass0", IOobject::NO_READ), np);
+
+    label nF = 0;
+    label i = 0;
+    for (const ReactingHeterogeneousParcel<ParcelType>& p : c)
+    {
+        mass0[i] = p.mass0_;
+        if (!i)
+        {
+            // Assume same size throughout
+            nF = p.F().size();
+        }
+        ++i;
+    }
+    mass0.write(np > 0);
+
+    for (label i = 0; i < nF; i++)
+    {
+        IOField<scalar> F
+        (
+            c.fieldIOobject
+            (
+                "F" + name(i),
+                IOobject::NO_READ
+            ),
+            np
+        );
+
+        for (const ReactingHeterogeneousParcel<ParcelType>& p0 : c)
+        {
+            F = p0.F()[i];
+        }
+
+        F.write(np > 0);
+    }
+
+    const label idSolid = compModel.idSolid();
+    const wordList& solidNames = compModel.componentNames(idSolid);
+
+    forAll(solidNames, j)
+    {
+        IOField<scalar> Y
+        (
+            c.fieldIOobject
+            (
+                "Y" + solidNames[j],
+                IOobject::NO_READ
+            ),
+            np
+        );
+
+        label i = 0;
+        for (const ReactingHeterogeneousParcel<ParcelType>& p0 : c)
+        {
+            Y[i] = p0.Y()[j];
+            ++i;
+        }
+
+        Y.write(np > 0);
+    }
+}
+
+
+template<class ParcelType>
+template<class CloudType>
+void Foam::ReactingHeterogeneousParcel<ParcelType>::writeObjects
+(
+    const CloudType& c,
+    objectRegistry& obr
+)
+{
+    ParcelType::writeObjects(c, obr);
+}
+
+
+template<class ParcelType>
+template<class CloudType, class CompositionType>
+void Foam::ReactingHeterogeneousParcel<ParcelType>::writeObjects
+(
+    const CloudType& c,
+    const CompositionType& compModel,
+    objectRegistry& obr
+)
+{
+    //ParcelType::writeObjects(c, obr);
+    // Skip Reacting layer
+    ThermoParcel<KinematicParcel<particle>>::writeObjects(c, obr);
+
+    label np = c.size();
+
+    // WIP
+    label nF = 0;
+    for (const ReactingHeterogeneousParcel<ParcelType>& p0 : c)
+    {
+        nF = p0.F().size();
+        break;
+    }
+
+    if (np > 0)
+    {
+        for (label i = 0; i < nF; i++)
+        {
+            const word fieldName = "F" + name(i);
+            IOField<scalar>& F
+            (
+                cloud::createIOField<scalar>(fieldName, np, obr)
+            );
+
+            label j = 0;
+            for (const ReactingHeterogeneousParcel<ParcelType>& p0 : c)
+            {
+                F[j] = p0.F()[i];
+                ++j;
+            }
+        }
+
+        const label idSolid = compModel.idSolid();
+        const wordList& solidNames = compModel.componentNames(idSolid);
+        forAll(solidNames, j)
+        {
+            const word fieldName = "Y" + solidNames[j];
+            IOField<scalar>& Y
+            (
+                cloud::createIOField<scalar>(fieldName, np, obr)
+            );
+
+            label i = 0;
+            for (const ReactingHeterogeneousParcel<ParcelType>& p0 : c)
+            {
+                Y[i] = p0.Y()[j];
+                ++i;
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+template<class ParcelType>
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const ReactingHeterogeneousParcel<ParcelType>& p
+)
+{
+    scalarField F(p.F());
+    if (os.format() == IOstream::ASCII)
+    {
+        os  << static_cast<const ParcelType&>(p)
+            << token::SPACE << F;
+    }
+    else
+    {
+        os  << static_cast<const ParcelType&>(p);
+        os  << F ;
+    }
+
+    os.check(FUNCTION_NAME);
+    return os;
+}
+
+
+// ************************************************************************* //

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2017 OpenFOAM Foundation
@@ -311,6 +311,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         td,
         dt,
         d0,
+        Res,
+        mus/rhos,
         T0,
         mass0,
         canCombust_,
@@ -607,6 +609,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
     trackingData& td,
     const scalar dt,
     const scalar d,
+    const scalar Re,
+    const scalar nu,
     const scalar T,
     const scalar mass,
     const label canCombust,
@@ -644,6 +648,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
     const scalar hReaction = cloud.surfaceReaction().calculate
     (
         dt,
+        Re,
+        nu,
         this->cell(),
         d,
         T,

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2016-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2017 OpenFOAM Foundation
@@ -52,18 +52,13 @@ SourceFiles
 namespace Foam
 {
 
-template<class ParcelType>
+template<class ParcelType> class ReactingMultiphaseParcel;
-class ReactingMultiphaseParcel;
 
 template<class ParcelType>
-Ostream& operator<<
+Ostream& operator<<(Ostream&, const ReactingMultiphaseParcel<ParcelType>&);
-(
-    Ostream&,
-    const ReactingMultiphaseParcel<ParcelType>&
-);
 
 /*---------------------------------------------------------------------------*\
-                 Class ReactingMultiphaseParcel Declaration
+                  Class ReactingMultiphaseParcel Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class ParcelType>
@@ -245,6 +240,8 @@ protected:
             trackingData& td,
             const scalar dt,           // timestep
             const scalar d,            // diameter
+            const scalar Re,           // Re
+            const scalar nu,           // nu
             const scalar T,            // temperature
             const scalar mass,         // mass
             const label canCombust,     // 'can combust' flag

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

@@ -138,7 +138,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::readFields
         label i = 0;
         for (ReactingMultiphaseParcel<ParcelType>& p : c)
         {
-            p.YGas_[j] = YGas[i++]/(p.Y()[GAS] + ROOTVSMALL);
+            p.YGas_[j] = YGas[i]/(p.Y()[GAS] + ROOTVSMALL);
+            ++i;
         }
     }
     // Populate YLiquid for each parcel
@@ -157,7 +158,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::readFields
         label i = 0;
         for (ReactingMultiphaseParcel<ParcelType>& p : c)
         {
-            p.YLiquid_[j] = YLiquid[i++]/(p.Y()[LIQ] + ROOTVSMALL);
+            p.YLiquid_[j] = YLiquid[i]/(p.Y()[LIQ] + ROOTVSMALL);
+            ++i;
         }
     }
     // Populate YSolid for each parcel
@@ -176,7 +178,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::readFields
         label i = 0;
         for (ReactingMultiphaseParcel<ParcelType>& p : c)
         {
-            p.YSolid_[j] = YSolid[i++]/(p.Y()[SLD] + ROOTVSMALL);
+            p.YSolid_[j] = YSolid[i]/(p.Y()[SLD] + ROOTVSMALL);
+            ++i;
         }
     }
 }
@@ -223,7 +226,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::writeFields
             label i = 0;
             for (const ReactingMultiphaseParcel<ParcelType>& p0 : c)
             {
-                YGas[i++] = p0.YGas()[j]*p0.Y()[GAS];
+                YGas[i] = p0.YGas()[j]*p0.Y()[GAS];
+                ++i;
             }
 
             YGas.write(np > 0);
@@ -246,7 +250,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::writeFields
             label i = 0;
             for (const ReactingMultiphaseParcel<ParcelType>& p0 : c)
             {
-                YLiquid[i++] = p0.YLiquid()[j]*p0.Y()[LIQ];
+                YLiquid[i] = p0.YLiquid()[j]*p0.Y()[LIQ];
+                ++i;
             }
 
             YLiquid.write(np > 0);
@@ -269,7 +274,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::writeFields
             label i = 0;
             for (const ReactingMultiphaseParcel<ParcelType>& p0 : c)
             {
-                YSolid[i++] = p0.YSolid()[j]*p0.Y()[SLD];
+                YSolid[i] = p0.YSolid()[j]*p0.Y()[SLD];
+                ++i;
             }
 
             YSolid.write(np > 0);
@@ -321,7 +327,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::writeObjects
             label i = 0;
             for (const ReactingMultiphaseParcel<ParcelType>& p0 : c)
             {
-                YGas[i++] = p0.YGas()[j]*p0.Y()[GAS];
+                YGas[i] = p0.YGas()[j]*p0.Y()[GAS];
+                ++i;
             }
         }
 
@@ -338,7 +345,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::writeObjects
             label i = 0;
             for (const ReactingMultiphaseParcel<ParcelType>& p0 : c)
             {
-                YLiquid[i++] = p0.YLiquid()[j]*p0.Y()[LIQ];
+                YLiquid[i] = p0.YLiquid()[j]*p0.Y()[LIQ];
+                ++i;
             }
         }
 
@@ -355,7 +363,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::writeObjects
             label i = 0;
             for (const ReactingMultiphaseParcel<ParcelType>& p0 : c)
             {
-                YSolid[i++] = p0.YSolid()[j]*p0.Y()[SLD];
+                YSolid[i] = p0.YSolid()[j]*p0.Y()[SLD];
+                ++i;
             }
         }
     }

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

@@ -0,0 +1,74 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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::basicHeterogeneousReactingParcel
+
+Description
+    Definition of reacting heterogeneous parcel
+
+SourceFiles
+    basicHeterogeneousReactingParcel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef basicHeterogeneousReactingParcel_H
+#define basicHeterogeneousReactingParcel_H
+
+#include "contiguous.H"
+#include "particle.H"
+#include "KinematicParcel.H"
+#include "ThermoParcel.H"
+#include "ReactingParcel.H"
+#include "ReactingHeterogeneousParcel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef ReactingHeterogeneousParcel
+    <
+        ReactingParcel
+        <
+            ThermoParcel
+            <
+                KinematicParcel
+                <
+                    particle
+                >
+            >
+        >
+    > basicHeterogeneousReactingParcel;
+
+    template<>
+    inline bool contiguous<basicHeterogeneousReactingParcel>()
+    {
+        return false;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/parcels/derived/basicHeterogeneousReactingParcel/defineBasicHeterogeneousReactingParcel.C

@@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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 "basicHeterogeneousReactingParcel.H"
+#include "Cloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTemplateTypeNameAndDebug(basicHeterogeneousReactingParcel, 0);
+    defineTemplateTypeNameAndDebug(Cloud<basicHeterogeneousReactingParcel>, 0);
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/parcels/derived/basicHeterogeneousReactingParcel/makeBasicHeterogeneousReactingParcelSubmodels.C

@@ -0,0 +1,69 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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 "basicHeterogeneousReactingCloud.H"
+
+#include "makeParcelCloudFunctionObjects.H"
+
+// Kinematic
+#include "makeThermoParcelForces.H" // thermo variant
+#include "makeParcelDispersionModels.H"
+#include "makeReactingParcelInjectionModels.H" // Reacting variant
+#include "makeParcelPatchInteractionModels.H"
+#include "makeParcelStochasticCollisionModels.H"
+#include "makeReactingParcelSurfaceFilmModels.H" // Reacting variant
+#include "makeHeterogeneousReactingParcelHeterogeneousReactingModels.H"
+
+// Thermodynamic
+#include "makeParcelHeatTransferModels.H"
+
+// Reacting
+#include "makeReactingMultiphaseParcelCompositionModels.H"
+#include "makeReactingParcelPhaseChangeModels.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makeParcelCloudFunctionObjects(basicHeterogeneousReactingCloud);
+
+// Kinematic sub-models
+makeThermoParcelForces(basicHeterogeneousReactingCloud);
+makeParcelDispersionModels(basicHeterogeneousReactingCloud);
+makeReactingParcelInjectionModels(basicHeterogeneousReactingCloud);
+makeParcelPatchInteractionModels(basicHeterogeneousReactingCloud);
+makeParcelStochasticCollisionModels(basicHeterogeneousReactingCloud);
+makeReactingParcelSurfaceFilmModels(basicHeterogeneousReactingCloud);
+
+// Thermo sub-models
+makeParcelHeatTransferModels(basicHeterogeneousReactingCloud);
+
+// Reacting sub-models
+makeReactingMultiphaseParcelCompositionModels(basicHeterogeneousReactingCloud);
+makeReactingParcelPhaseChangeModels(basicHeterogeneousReactingCloud);
+makeHeterogeneousReactingParcelHeterogeneousReactingModels
+(
+    basicHeterogeneousReactingCloud
+);
+
+// ************************************************************************* //

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

@@ -0,0 +1,47 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef makeHeterogeneousReactingParcelHeterogeneousReactingModels_H
+#define makeHeterogeneousReactingParcelHeterogeneousReactingModels_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "NoheterogeneousReacting.H"
+#include "MUCSheterogeneousRate.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeHeterogeneousReactingParcelHeterogeneousReactingModels(CloudType)  \
+                                                                               \
+    makeHeterogeneousReactingModel(CloudType);                                 \
+    makeHeterogeneousReactingModelType(NoheterogeneousReacting, CloudType);    \
+    makeHeterogeneousReactingModelType(MUCSheterogeneousRate, CloudType);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/HeterogeneousReactingModel/HeterogeneousReactingModel/HeterogeneousReactingModel.C

@@ -0,0 +1,107 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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 "HeterogeneousReactingModel.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::HeterogeneousReactingModel<CloudType>::HeterogeneousReactingModel
+(
+    CloudType& owner
+)
+:
+    CloudSubModelBase<CloudType>(owner),
+    dMass_(0.0),
+    nF_(0)
+{}
+
+
+template<class CloudType>
+Foam::HeterogeneousReactingModel<CloudType>::HeterogeneousReactingModel
+(
+    const dictionary& dict,
+    CloudType& owner,
+    const word& type
+)
+:
+    CloudSubModelBase<CloudType>(owner, dict, typeName, type),
+    dMass_(0.0),
+    nF_(this->coeffDict().template lookupOrDefault<label>("nF", 1))
+{}
+
+
+template<class CloudType>
+Foam::HeterogeneousReactingModel<CloudType>::HeterogeneousReactingModel
+(
+    const HeterogeneousReactingModel<CloudType>& srm
+)
+:
+    CloudSubModelBase<CloudType>(srm),
+    dMass_(srm.dMass_),
+    nF_(srm.nF_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::HeterogeneousReactingModel<CloudType>::addToSurfaceReactionMass
+(
+    const scalar dMass
+)
+{
+    dMass_ += dMass;
+}
+
+
+template<class CloudType>
+Foam::label Foam::HeterogeneousReactingModel<CloudType>::nF() const
+{
+    return nF_;
+}
+
+
+template<class CloudType>
+void Foam::HeterogeneousReactingModel<CloudType>::info(Ostream& os)
+{
+    const scalar mass0 = this->template getBaseProperty<scalar>("mass");
+    const scalar massTotal = mass0 + returnReduce(dMass_, sumOp<scalar>());
+
+    Info<< "    Mass transfer surface reaction  = " << massTotal << nl;
+
+    if (this->writeTime())
+    {
+        this->setBaseProperty("mass", massTotal);
+        dMass_ = 0.0;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "HeterogeneousReactingModelNew.C"
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/HeterogeneousReactingModel/HeterogeneousReactingModel/HeterogeneousReactingModel.H

@@ -0,0 +1,218 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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::HeterogeneousReactingModel
+
+Group
+    grpLagrangianIntermediateHeterogeneousReactionSubModels
+
+Description
+    Base class for heterogeneous reacting models
+
+SourceFiles
+    HeterogeneousReactingModel.C
+    HeterogeneousReactingModelNew.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef HeterogeneousReactingModel_H
+#define HeterogeneousReactingModel_H
+
+#include "IOdictionary.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+#include "CloudSubModelBase.H"
+#include "scalarField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class HeterogeneousReactingModel Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class HeterogeneousReactingModel
+:
+    public CloudSubModelBase<CloudType>
+{
+protected:
+
+    // Protected data
+
+        //- Net mass of lagrangian phase consumed
+        scalar dMass_;
+
+        //- Number of progress variables
+        label nF_;
+
+
+public:
+
+    //-Runtime type information
+    TypeName("heterogeneousReactionModel");
+
+
+    //- Declare runtime constructor selection table
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        HeterogeneousReactingModel,
+        dictionary,
+        (
+            const dictionary& dict,
+            CloudType& cloud
+        ),
+        (dict, cloud)
+    );
+
+
+    // Constructors
+
+        //- Construct null from owner
+        HeterogeneousReactingModel(CloudType& owner);
+
+        //- Construct from dictionary
+        HeterogeneousReactingModel
+        (
+            const dictionary& dict,
+            CloudType& cloud,
+            const word& type
+        );
+
+        //- Construct copy
+        HeterogeneousReactingModel
+        (
+            const HeterogeneousReactingModel<CloudType>& srm
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr
+        <
+            HeterogeneousReactingModel<CloudType>
+        > clone() const = 0;
+
+
+    //- Destructor
+    virtual ~HeterogeneousReactingModel() = default;
+
+
+    //- Selector
+    static autoPtr<HeterogeneousReactingModel<CloudType>> New
+    (
+        const dictionary& dict,
+        CloudType& cloud
+    );
+
+
+    // Member Functions
+
+        //- Update surface reactions
+        //  Returns the heat of reaction
+        virtual scalar calculate
+        (
+            const scalar dt,
+            const scalar Re,
+            const scalar nu,
+            const label celli,
+            const scalar d,
+            const scalar T,
+            const scalar Tc,
+            const scalar pc,
+            const scalar rhoc,
+            const scalar mass,
+            const scalarField& YSolid,
+            scalarField& F,
+            const scalar N,
+            scalar& NCpW,
+            scalarField& dMassSolid,
+            scalarField& dMassSRCarrier
+        ) const = 0;
+
+
+        //- Add to devolatilisation mass
+        void addToSurfaceReactionMass(const scalar dMass);
+
+        //- Number of reactions in the model
+        virtual label nReactions() const = 0;
+
+        //- Number of progress variable
+        virtual label nF() const;
+
+        //- Write injection info to stream
+        virtual void info(Ostream& os);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeHeterogeneousReactingModel(CloudType)                              \
+                                                                               \
+    typedef Foam::CloudType::reactingHeterogeneousCloudType                    \
+        reactingHeterogeneousCloudType;                                        \
+    defineNamedTemplateTypeNameAndDebug                                        \
+    (                                                                          \
+        Foam::HeterogeneousReactingModel<reactingHeterogeneousCloudType>,      \
+        0                                                                      \
+    );                                                                         \
+    namespace Foam                                                             \
+    {                                                                          \
+        defineTemplateRunTimeSelectionTable                                    \
+        (                                                                      \
+            HeterogeneousReactingModel<reactingHeterogeneousCloudType>,        \
+            dictionary                                                         \
+        );                                                                     \
+    }
+
+
+#define makeHeterogeneousReactingModelType(SS, CloudType)                      \
+                                                                               \
+    typedef Foam::CloudType::reactingHeterogeneousCloudType                    \
+        reactingHeterogeneousCloudType;                                        \
+    defineNamedTemplateTypeNameAndDebug                                        \
+        (Foam::SS<reactingHeterogeneousCloudType>, 0);                         \
+                                                                               \
+    Foam::HeterogeneousReactingModel<reactingHeterogeneousCloudType>::         \
+        adddictionaryConstructorToTable                                        \
+        <Foam::SS<reactingHeterogeneousCloudType>>                             \
+        add##SS##CloudType##reactingHeterogeneousCloudType##ConstructorToTable_;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "HeterogeneousReactingModel.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/HeterogeneousReactingModel/HeterogeneousReactingModel/HeterogeneousReactingModelNew.C

@@ -0,0 +1,61 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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 "HeterogeneousReactingModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::autoPtr<Foam::HeterogeneousReactingModel<CloudType>>
+Foam::HeterogeneousReactingModel<CloudType>::New
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+{
+    const word modelType(dict.get<word>("heterogeneousReactingModel"));
+
+    Info<< "Selecting surface reaction model " << modelType << endl;
+
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
+
+    if (!cstrIter.found())
+    {
+        FatalErrorInFunction
+            << "Unknown heterogeneousReactingModel type "
+            << modelType << nl << nl
+            << "Valid types :" << nl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<HeterogeneousReactingModel<CloudType>>
+    (
+        cstrIter()(dict, owner)
+    );
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/HeterogeneousReactingModel/MUCSheterogeneousRate/MUCSheterogeneousRate.C

@@ -0,0 +1,239 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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 "MUCSheterogeneousRate.H"
+#include "mathematicalConstants.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::MUCSheterogeneousRate<CloudType>::MUCSheterogeneousRate
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    HeterogeneousReactingModel<CloudType>(dict, owner, typeName),
+    D12_(this->coeffDict().getScalar("D12")),
+    epsilon_(this->coeffDict().getScalar("epsilon")),
+    gamma_(this->coeffDict().getScalar("gamma")),
+    sigma_(this->coeffDict().getScalar("sigma")),
+    E_(this->coeffDict().getScalar("E")),
+    A_(this->coeffDict().getScalar("A")),
+    Aeff_(this->coeffDict().getScalar("Aeff")),
+    Ea_(this->coeffDict().getScalar("Ea")),
+    nuFuel_(this->coeffDict().getScalar("nuFuel")),
+    nuOx_(this->coeffDict().getScalar("nuOx")),
+    nuProd_(this->coeffDict().getScalar("nuProd")),
+    O2GlobalId_(owner.composition().carrierId("O2")),
+    FuelLocalId_(-1),
+    ProdLocalId_(-1),
+    WO2_(0.0)
+{
+    // Determine Cs ids
+    label idSolid = owner.composition().idSolid();
+    FuelLocalId_ =
+        owner.composition().localId
+        (
+            idSolid,
+            this->coeffDict().getWord("fuel")
+        );
+
+    ProdLocalId_ =
+        owner.composition().localId
+        (
+            idSolid,
+            this->coeffDict().getWord("product")
+        );
+
+    // Set local copies of thermo properties
+    WO2_ = owner.thermo().carrier().W(O2GlobalId_);
+}
+
+
+template<class CloudType>
+Foam::MUCSheterogeneousRate<CloudType>::MUCSheterogeneousRate
+(
+    const MUCSheterogeneousRate<CloudType>& srm
+)
+:
+    HeterogeneousReactingModel<CloudType>(srm),
+    D12_(srm.D12_),
+    epsilon_(srm.epsilon_),
+    gamma_(srm.gamma_),
+    sigma_(srm.sigma_),
+    E_(srm.E_),
+    A_(srm.A_),
+    Aeff_(srm.Aeff_),
+    Ea_(srm.Ea_),
+    nuFuel_(srm.nuFuel_),
+    nuOx_(srm.nuOx_),
+    nuProd_(srm.nuProd_),
+    O2GlobalId_(srm.O2GlobalId_),
+    FuelLocalId_(srm.FuelLocalId_),
+    ProdLocalId_(srm.ProdLocalId_),
+    WO2_(srm.WO2_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::scalar Foam::MUCSheterogeneousRate<CloudType>::calculate
+(
+    const scalar dt,
+    const scalar Re,
+    const scalar nu,
+    const label celli,
+    const scalar d,
+    const scalar T,
+    const scalar Tc,
+    const scalar pc,
+    const scalar rhoc,
+    const scalar mass,
+    const scalarField& YSolid,
+    scalarField& F,
+    const scalar N,
+    scalar& NCpW,
+    scalarField& dMassSolid,
+    scalarField& dMassSRCarrier
+) const
+{
+    // Fraction of remaining combustible material
+    const scalar fComb = YSolid[FuelLocalId_];
+
+    // Surface combustion until combustible fraction is consumed
+    if (fComb < SMALL)
+    {
+        return 0.0;
+    }
+
+    const SLGThermo& thermo = this->owner().thermo();
+    const auto& composition = this->owner().composition();
+
+    const scalar WFuel = composition.solids().properties()[FuelLocalId_].W();
+    const scalar WProd = composition.solids().properties()[ProdLocalId_].W();
+
+    // O2 concentration [Kmol/m3]
+    const scalar Cb =
+        thermo.carrier().Y(O2GlobalId_)[celli]*rhoc/WO2_;
+
+    if (Cb < SMALL)
+    {
+        return 0.0;
+    }
+
+    // Reaction constant
+    const scalar k = A_*exp(-Ea_/(RR*T));
+
+    // Effective diffussivity
+    const scalar Deff = D12_*epsilon_/gamma_;
+
+     // Schmidt number
+    const scalar Sc = nu/(D12_ + ROOTVSMALL);
+
+    // Mass transfer coefficient [m/s]
+    const scalar alpha =
+        (2.0 + 0.6*Foam::sqrt(Re)*cbrt(Sc))*D12_/(d + ROOTVSMALL);
+
+    const scalar r = d/2;
+
+    const scalar f = F[FuelLocalId_];
+
+    const scalar rhof = composition.solids().properties()[FuelLocalId_].rho();
+
+    const scalar deltaRho0 = (nuOx_/nuFuel_)*rhof/WFuel;
+
+    // Progress variable rate
+    const scalar dfdt =
+        Aeff_*(Cb/deltaRho0)
+       /(
+           r/3/alpha
+         + sqr(r)*(1/cbrt(1-f)-1)/3/Deff
+         - (1/sqr(cbrt(1-f)))*r/k/sigma_/E_/3
+        );
+
+    // Update new progress variable
+    F[FuelLocalId_] += dfdt*dt;
+
+    // Interface radius
+    const scalar ri = r*cbrt(1-f);
+
+    // Interface radius rate
+    //const scalar dridt = -dfdt*(r/3)*pow(1-f, -2/3);
+    const scalar dridt = -dfdt*(pow3(r)/3)/sqr(ri);
+
+    // O2 flux [Kmol/sec]
+    const scalar q02 = deltaRho0*4*constant::mathematical::pi*sqr(ri)*dridt;
+
+    // Calculate the number of molar units reacted [Kmol]
+    const scalar dOmega = q02*dt;
+
+    // Heat of Reaction
+    const scalar Hc =
+        composition.solids().properties()[ProdLocalId_].Hf()
+      - composition.solids().properties()[FuelLocalId_].Hf();
+
+    //Stoichiometric mass ratio for fuel
+    const scalar sFuel = nuFuel_/(nuOx_);
+
+    //Stoichiometric mass ratio for product
+    const scalar sProd = nuProd_/(nuOx_);
+
+    // Add to carrier phase mass transfer [Kg]
+    dMassSRCarrier[O2GlobalId_] += dOmega*WO2_;
+
+    // Remove to particle mass transfer
+    dMassSolid[FuelLocalId_] -= dOmega*WFuel*sFuel;
+
+    // Add to particle product
+    dMassSolid[ProdLocalId_] += dOmega*WProd*sProd;
+
+    if (debug)
+    {
+        Pout<< "mass    = " << mass << nl
+            << "fComb   = " << fComb << nl
+            << "dfdt    = " << dfdt << nl
+            << "F       = " << F[FuelLocalId_] << nl
+            << "ri      = " << ri << nl
+            << "dridt   = " << dridt << nl
+            << "q02     = " << q02 << nl
+            << "dOmega  = " << dOmega << nl
+            << "Hr      = " << dOmega*WFuel*sFuel*Hc << endl;
+    }
+
+    // Heat of reaction [J]
+    return -dOmega*WFuel*sFuel*Hc;
+}
+
+
+template<class CloudType>
+Foam::label Foam::MUCSheterogeneousRate<CloudType>::nReactions() const
+{
+    return 1;
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/HeterogeneousReactingModel/MUCSheterogeneousRate/MUCSheterogeneousRate.H

@@ -0,0 +1,198 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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::MUCSheterogeneousRate
+
+Group
+    grpLagrangianIntermediateSurfaceReactionSubModels
+
+Description
+    Heteregeneous noncatalytic reaction MUCS approach.
+    Reference:
+        D. Papanastassiou and G. Bitsianes,
+        Modelling of Heterogeneous Gas-Solid Reactions,
+        Metallurgical Transsactions, 480. Volume 4. 1973
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef MUCSheterogeneousRate_H
+#define MUCSheterogeneousRate_H
+
+#include "HeterogeneousReactingModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declarations
+template<class CloudType>
+class MUCSheterogeneousRate;
+
+/*---------------------------------------------------------------------------*\
+                  Class MUCSheterogeneousRate Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class MUCSheterogeneousRate
+:
+    public HeterogeneousReactingModel<CloudType>
+{
+    // Private Data
+
+        // Model constants
+
+            //- Binary diffusivity [m2/s]
+            const scalar D12_;
+
+            //- Porosity of product layer
+            const scalar epsilon_;
+
+            //- Tortuosity
+            const scalar gamma_;
+
+            //- Spread factor in pellet
+            const scalar sigma_;
+
+            //- Effectiveness factor
+            const scalar E_;
+
+            //- Kinetic rate coefficient [m/s]
+            const scalar A_;
+
+            //- Effective areas [-]
+            const scalar Aeff_;
+
+            //- Kinetic rate activation energy [J/kmol]
+            const scalar Ea_;
+
+            //- Stoichiomatric solid reactant
+            const scalar nuFuel_;
+
+            //- Stoichiomatric oxydizer reactant
+            const scalar nuOx_;
+
+            //- Stoichiomatric solid product
+            const scalar nuProd_;
+
+
+        // Addressing
+
+            //- O2 position in global list
+            label O2GlobalId_;
+
+            //- Fuel solid local Id
+            label FuelLocalId_;
+
+            //- Product solid local Id
+            label ProdLocalId_;
+
+
+        // Local copies of thermo properties
+
+            //- Molecular weight of O2 [kg/kmol]
+            scalar WO2_;
+
+
+
+public:
+
+    //- Runtime type information
+    TypeName("MUCSheterogeneousRate");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        MUCSheterogeneousRate
+        (
+            const dictionary& dict,
+            CloudType& owner
+        );
+
+        //- Construct copy
+        MUCSheterogeneousRate
+        (
+            const MUCSheterogeneousRate<CloudType>& srm
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<HeterogeneousReactingModel<CloudType>> clone() const
+        {
+            return autoPtr<HeterogeneousReactingModel<CloudType>>
+            (
+                new MUCSheterogeneousRate<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~MUCSheterogeneousRate() = default;
+
+
+    // Member Functions
+
+        //- Update surface reactions
+        virtual scalar calculate
+        (
+            const scalar dt,
+            const scalar Re,
+            const scalar nu,
+            const label celli,
+            const scalar d,
+            const scalar T,
+            const scalar Tc,
+            const scalar pc,
+            const scalar rhoc,
+            const scalar mass,
+            const scalarField& YSolid,
+            scalarField& F,
+            const scalar N,
+            scalar& NCpW,
+            scalarField& dMassSolid,
+            scalarField& dMassSRCarrier
+        ) const;
+
+
+        //- Number of reactions in the model
+        virtual label nReactions() const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "MUCSheterogeneousRate.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/HeterogeneousReactingModel/NoheterogeneousReacting/NoheterogeneousReacting.C

@@ -0,0 +1,92 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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 "NoheterogeneousReacting.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::NoheterogeneousReacting<CloudType>::NoheterogeneousReacting
+(
+    const dictionary&,
+    CloudType& owner
+)
+:
+    HeterogeneousReactingModel<CloudType>(owner)
+{}
+
+
+template<class CloudType>
+Foam::NoheterogeneousReacting<CloudType>::NoheterogeneousReacting
+(
+    const NoheterogeneousReacting<CloudType>& srm
+)
+:
+    HeterogeneousReactingModel<CloudType>(srm.owner_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+bool Foam::NoheterogeneousReacting<CloudType>::active() const
+{
+    return false;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::NoheterogeneousReacting<CloudType>::calculate
+(
+    const scalar,
+    const scalar,
+    const scalar,
+    const label,
+    const scalar,
+    const scalar,
+    const scalar,
+    const scalar,
+    const scalar,
+    const scalar,
+    const scalarField&,
+    scalarField&,
+    const scalar,
+    scalar&,
+    scalarField&,
+    scalarField&
+) const
+{
+    return 0;
+}
+
+
+template<class CloudType>
+Foam::label Foam::NoheterogeneousReacting<CloudType>::nReactions() const
+{
+    return 0;
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/HeterogeneousReactingModel/NoheterogeneousReacting/NoheterogeneousReacting.H

@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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::NoheterogeneousReacting
+
+Group
+    grpLagrangianIntermediateSurfaceReactionSubModels
+
+Description
+    Dummy surface reaction model for 'none'
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef NoheterogeneousReacting_H
+#define NoheterogeneousReacting_H
+
+#include "HeterogeneousReactingModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class NoheterogeneousReacting Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class NoheterogeneousReacting
+:
+    public HeterogeneousReactingModel<CloudType>
+{
+public:
+
+    //- Runtime type information
+    TypeName("none");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        NoheterogeneousReacting(const dictionary& dict, CloudType& owner);
+
+        //- Construct copy
+        NoheterogeneousReacting
+        (
+            const NoheterogeneousReacting<CloudType>& srm
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<HeterogeneousReactingModel<CloudType>> clone() const
+        {
+            return autoPtr<HeterogeneousReactingModel<CloudType>>
+            (
+                new NoheterogeneousReacting<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~NoheterogeneousReacting() = default;
+
+
+    // Member Functions
+
+        //- Flag to indicate whether model activates devolatisation model
+        virtual bool active() const;
+
+        //- Update surface reactions
+        virtual scalar calculate
+        (
+            const scalar dt,
+            const scalar Re,
+            const scalar nu,
+            const label celli,
+            const scalar d,
+            const scalar T,
+            const scalar Tc,
+            const scalar pc,
+            const scalar rhoc,
+            const scalar mass,
+            const scalarField& YSolid,
+            scalarField& F,
+            const scalar N,
+            scalar& NCpW,
+            scalarField& dMassSolid,
+            scalarField& dMassSRCarrier
+        ) const;
+
+        //- Number of reactions in the model
+        virtual label nReactions() const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "NoheterogeneousReacting.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2017 OpenFOAM Foundation
@@ -50,13 +50,6 @@ Foam::NoSurfaceReaction<CloudType>::NoSurfaceReaction
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::NoSurfaceReaction<CloudType>::~NoSurfaceReaction()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -69,6 +62,8 @@ bool Foam::NoSurfaceReaction<CloudType>::active() const
 template<class CloudType>
 Foam::scalar Foam::NoSurfaceReaction<CloudType>::calculate
 (
+    const scalar,
+    const scalar,
     const scalar,
     const label,
     const scalar,

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -77,7 +77,7 @@ public:
 
 
     //- Destructor
-    virtual ~NoSurfaceReaction();
+    virtual ~NoSurfaceReaction() = default;
 
 
     // Member Functions
@@ -89,6 +89,8 @@ public:
         virtual scalar calculate
         (
             const scalar dt,
+            const scalar Re,
+            const scalar nu,
             const label celli,
             const scalar d,
             const scalar T,

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -36,7 +36,7 @@ Foam::SurfaceReactionModel<CloudType>::SurfaceReactionModel
 )
 :
     CloudSubModelBase<CloudType>(owner),
-    dMass_(0.0)
+    dMass_(0)
 {}
 
 
@@ -49,7 +49,7 @@ Foam::SurfaceReactionModel<CloudType>::SurfaceReactionModel
 )
 :
     CloudSubModelBase<CloudType>(owner, dict, typeName, type),
-    dMass_(0.0)
+    dMass_(0)
 {}
 
 
@@ -64,13 +64,6 @@ Foam::SurfaceReactionModel<CloudType>::SurfaceReactionModel
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::SurfaceReactionModel<CloudType>::~SurfaceReactionModel()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class CloudType>

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -63,7 +63,7 @@ class SurfaceReactionModel
 {
 protected:
 
-    // Protected data
+    // Protected Data
 
         //- Mass of lagrangian phase converted
         scalar dMass_;
@@ -110,7 +110,7 @@ public:
 
 
     //- Destructor
-    virtual ~SurfaceReactionModel();
+    virtual ~SurfaceReactionModel() = default;
 
 
     //- Selector
@@ -128,6 +128,8 @@ public:
         virtual scalar calculate
         (
             const scalar dt,
+            const scalar Re,
+            const scalar nu,
             const label celli,
             const scalar d,
             const scalar T,

src/lagrangian/turbulence/Make/files

@@ -19,4 +19,7 @@ $(REACTINGMPPARCEL)/makeBasicReactingMultiphaseParcelSubmodels.C
 KINEMATICMPPICPARCEL=$(DERIVEDPARCELS)/basicKinematicMPPICParcel
 $(KINEMATICMPPICPARCEL)/makeBasicKinematicMPPICParcelSubmodels.C
 
+HETEROGENEOUSREACTINGPARCEL=$(DERIVEDPARCELS)/basicHeterogeneousReactingParcel
+$(HETEROGENEOUSREACTINGPARCEL)/makeBasicHeterogeneousParcelSubmodels.C
+
 LIB = $(FOAM_LIBBIN)/liblagrangianTurbulence

src/lagrangian/turbulence/parcels/derived/basicHeterogeneousReactingParcel/makeBasicHeterogeneousParcelSubmodels.C

@@ -0,0 +1,40 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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 "basicHeterogeneousReactingCloud.H"
+
+#include "makeParcelTurbulenceDispersionModels.H"
+#include "makeThermoParcelTurbulenceForces.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makeThermoParcelTurbulenceForces(basicHeterogeneousReactingCloud);
+    makeParcelTurbulenceDispersionModels(basicHeterogeneousReactingCloud);
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/thermophysicalProperties/solidProperties/C/C.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2017 OpenFOAM Foundation
@@ -41,7 +41,7 @@ namespace Foam
 
 Foam::C::C()
 :
-    solidProperties(2010, 710, 0.04, 0.0, 1.0)
+    solidProperties(2010, 710, 0.04, 0.0, 1.0, 12.011)
 {
     if (debug)
     {

src/thermophysicalModels/thermophysicalProperties/solidProperties/CaCO3/CaCO3.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2017 OpenFOAM Foundation
@@ -41,7 +41,7 @@ namespace Foam
 
 Foam::CaCO3::CaCO3()
 :
-    solidProperties(2710, 850, 1.3, 0.0, 1.0)
+    solidProperties(2710, 850, 1.3, 0.0, 1.0, 100.086)
 {
     if (debug)
     {

src/thermophysicalModels/thermophysicalProperties/solidProperties/ash/ash.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2017 OpenFOAM Foundation
@@ -41,7 +41,7 @@ namespace Foam
 
 Foam::ash::ash()
 :
-    solidProperties(2010, 710, 0.04, 0.0, 1.0)
+    solidProperties(2010, 710, 0.04, 0.0, 1.0, 12.011)
 {
     if (debug)
     {

src/thermophysicalModels/thermophysicalProperties/solidProperties/solidProperties/solidProperties.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2017 OpenFOAM Foundation
@@ -36,6 +36,7 @@ namespace Foam
     defineRunTimeSelectionTable(solidProperties, dictionary);
 }
 
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::solidProperties::solidProperties
@@ -44,14 +45,16 @@ Foam::solidProperties::solidProperties
     scalar Cp,
     scalar kappa,
     scalar Hf,
-    scalar emissivity
+    scalar emissivity,
+    scalar W
 )
 :
     rho_(rho),
     Cp_(Cp),
     kappa_(kappa),
     Hf_(Hf),
-    emissivity_(emissivity)
+    emissivity_(emissivity),
+    W_(W)
 {}
 
 
@@ -61,7 +64,8 @@ Foam::solidProperties::solidProperties(const dictionary& dict)
     Cp_(dict.get<scalar>("Cp")),
     kappa_(dict.getCompat<scalar>("kappa", {{"K", 1612}})),
     Hf_(dict.get<scalar>("Hf")),
-    emissivity_(dict.get<scalar>("emissivity"))
+    emissivity_(dict.get<scalar>("emissivity")),
+    W_(dict.get<scalar>("W"))
 {}
 
 
@@ -74,6 +78,7 @@ void Foam::solidProperties::readIfPresent(const dictionary& dict)
     dict.readIfPresentCompat("kappa", {{"K", 1612}}, kappa_);
     dict.readIfPresent("Hf_", Hf_);
     dict.readIfPresent("emissivity", emissivity_);
+    dict.readIfPresent("W", W_);
 }
 
 
@@ -83,7 +88,8 @@ void Foam::solidProperties::writeData(Ostream& os) const
         << Cp_ << token::SPACE
         << kappa_ << token::SPACE
         << Hf_ << token::SPACE
-        << emissivity_;
+        << emissivity_ << token::SPACE
+        << W_;
 }
 
 

src/thermophysicalModels/thermophysicalProperties/solidProperties/solidProperties/solidProperties.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2017 OpenFOAM Foundation
@@ -53,7 +53,7 @@ namespace Foam
 
 class solidProperties
 {
-    // Private data
+    // Private Data
 
         //- Density [kg/m3]
         scalar rho_;
@@ -70,6 +70,9 @@ class solidProperties
         //- Emissivity
         scalar emissivity_;
 
+        //- Molar weight [Kg/Kmol]
+        scalar W_;
+
 
 public:
 
@@ -107,7 +110,8 @@ public:
             scalar Cp,
             scalar kappa,
             scalar Hf,
-            scalar emissivity
+            scalar emissivity,
+            scalar W
         );
 
         //- Construct from dictionary
@@ -155,6 +159,9 @@ public:
             //- Emissivity []
             inline scalar emissivity() const;
 
+            //- Molar weight [Kg/Kmol]
+            inline scalar W() const;
+
 
     // I-O
 

src/thermophysicalModels/thermophysicalProperties/solidProperties/solidProperties/solidPropertiesI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2017 OpenFOAM Foundation
@@ -64,5 +64,10 @@ inline Foam::scalar Foam::solidProperties::emissivity() const
     return emissivity_;
 }
 
+inline Foam::scalar Foam::solidProperties::W() const
+{
+    return W_;
+}
+
 
 // ************************************************************************* //

tutorials/lagrangian/reactingHeterogenousParcelFoam/0/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;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/reactingHeterogenousParcelFoam/0/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.2;
+
+boundaryField
+{
+    walls
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 0.2;
+        value           uniform 0.2;
+    }
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 0.2;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/reactingHeterogenousParcelFoam/0/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 900.0;
+
+boundaryField
+{
+    walls
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 900.0;
+    }
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 900.0;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/reactingHeterogenousParcelFoam/0/U

@@ -0,0 +1,43 @@
+/*--------------------------------*- 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 (2.0 0 0);
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform (2.0 0 0);
+    }
+    outlet
+    {
+        type            pressureInletOutletVelocity;
+        value            $internalField;
+    }
+    walls
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/reactingHeterogenousParcelFoam/0/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;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/reactingHeterogenousParcelFoam/0/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;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/reactingHeterogenousParcelFoam/0/nut

@@ -0,0 +1,46 @@
+/*--------------------------------*- 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;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/reactingHeterogenousParcelFoam/0/omega

@@ -0,0 +1,48 @@
+/*--------------------------------*- 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;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/reactingHeterogenousParcelFoam/0/p

@@ -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      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 1e5;
+
+boundaryField
+{
+
+    inlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+    outlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+    walls
+    {
+        type            calculated;
+        value           $internalField;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/reactingHeterogenousParcelFoam/0/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            fixedFluxPressure;
+    }
+    outlet
+    {
+        type            totalPressure;
+        p0              $internalField;
+    }
+    walls
+    {
+        type            fixedFluxPressure;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/reactingHeterogenousParcelFoam/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
+
+
+// ************************************************************************* //

tutorials/lagrangian/reactingHeterogenousParcelFoam/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;
+
+
+// ************************************************************************* //

tutorials/lagrangian/reactingHeterogenousParcelFoam/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;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/foam.inp

@@ -0,0 +1,9 @@
+species
+(
+    O2
+    N2
+)
+;
+
+reactions
+{}

tutorials/lagrangian/reactingHeterogenousParcelFoam/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);
+
+
+// ************************************************************************* //

tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/particleTrackDict

@@ -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      particleTrackDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+cloud           reactingCloud1Tracks;
+
+fields
+(
+    d
+    U
+    T
+);
+
+
+// ************************************************************************* //

tutorials/lagrangian/reactingHeterogenousParcelFoam/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;
+
+// ************************************************************************* //

tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/reactingCloud1Properties

@@ -0,0 +1,184 @@
+/*--------------------------------*- 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
+{
+    rho0            5100;//Particle density (overwritten by composition)
+    T0              303; //Initial particle temperature
+    Cp0             850; //Initial particle Cp (overwritten by composition)
+
+    hRetentionCoeff 0;
+    constantVolume  true;
+}
+
+
+subModels
+{
+    particleForces
+    {}
+
+    injectionModels
+    {
+        // Mass flow rate : massTotal/duration
+        // Volume flow rate : Mass flow rate/particleRho
+        // parcelsPerSecond : Volume flow rate/particleVol
+        model1
+        {
+            type            patchInjection;
+            patch           inlet;
+            parcelBasisType mass;
+            U0              (0.1 0 0);
+            massTotal        30;
+            parcelsPerSecond 8442;
+            SOI              0;
+            duration         1;
+            flowRateProfile  constant 1;
+
+            // As we want 1 particle per parcel, this avoid
+            // accumulated vol if nParticles is below 1
+            minParticlesPerParcel   0.7;
+
+            sizeDistribution
+            {
+                type        fixedValue;
+                fixedValueDistribution
+                {
+                    value       0.011;
+                }
+            }
+        }
+    }
+
+    dispersionModel gradientDispersionRAS;
+
+    patchInteractionModel standardWallInteraction;
+
+    heatTransferModel RanzMarshall;
+
+    compositionModel singleMixtureFraction;
+
+    phaseChangeModel none;
+
+    stochasticCollisionModel none;
+
+    surfaceFilmModel none;
+
+    radiation       off;
+
+    standardWallInteractionCoeffs
+    {
+        type            rebound;
+    }
+
+    RanzMarshallCoeffs
+    {
+        BirdCorrection  off;
+    }
+
+    heterogeneousReactingModel  MUCSheterogeneousRate;
+
+    MUCSheterogeneousRateCoeffs
+    {
+        D12         2.724e-4; //m2/s
+        epsilon     0.41;
+        gamma       3.07;
+        sigma       1;
+        E           1;
+        A           3.14e4;   // m/s
+        Aeff        0.7;
+        Ea          1.651e5;  // J/kmol
+        O2          O2;
+
+        // nuFuel*Fe3O4 + nuOx*O2 => nuProd*Fe2O3
+        nuFuel      2.0;
+        nuProd      3.0;
+        nuOx        0.5;
+        fuel        Fe3O4;
+        product     Fe2O3;
+    }
+
+    singleMixtureFractionCoeffs
+    {
+        phases
+        (
+            gas
+            {
+            }
+            liquid
+            {
+            }
+            solid
+            {
+                Fe3O4 1;
+                Fe2O3 0;
+            }
+        );
+        YGasTot0        0;
+        YLiquidTot0     0;
+        YSolidTot0      1;
+    }
+}
+
+
+cloudFunctions
+{
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/thermophysicalProperties

@@ -0,0 +1,73 @@
+/*--------------------------------*- 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
+{}
+
+solids
+{
+    Fe3O4
+    {
+        defaultCoeffs   no;
+        Fe3O4Coeffs
+        {
+            rho             5100;
+            Cp              850;
+            kappa           0.04;
+            Hf              0;
+            emissivity      1.0;
+            W               232;
+        }
+    }
+
+    Fe2O3
+    {
+        defaultCoeffs   no;
+        Fe2O3Coeffs
+        {
+            rho             5100;
+            Cp              710;
+            kappa           0.04;
+            Hf              525e3; //Heat of reaction HfProd - HfFuel
+            emissivity      1.0;
+            W               159;   //Kg/Kmole
+        }
+    }
+}
+
+
+
+// ************************************************************************* //

tutorials/lagrangian/reactingHeterogenousParcelFoam/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  RAS;
+
+RAS
+{
+    RASModel            kOmegaSST; // kEpsilon;
+
+    turbulence          on;
+
+    printCoeffs         on;
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/reactingHeterogenousParcelFoam/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;
+    note            "Created Wed Jul  1 19:20:21 2009. Blocks = 8, cells = 9340, vertices = 36";
+    class           dictionary;
+    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)
+        );
+    }
+);
+
+
+// ************************************************************************* //

tutorials/lagrangian/reactingHeterogenousParcelFoam/system/controlDict

@@ -0,0 +1,81 @@
+/*--------------------------------*- 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;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+DebugSwitches
+{
+    MUCSheterogeneousRate   0;
+}
+
+application     reactingHeterogenousParcelFoam;
+
+startFrom       latestTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         5;
+
+deltaT          0.01;
+
+writeControl    adjustableRunTime;
+
+writeInterval   0.5;
+
+purgeWrite      0;
+
+writeFormat     binary;
+
+writePrecision  10;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+adjustTimeStep  yes;
+
+maxCo           1;
+
+maxDeltaT       0.03;
+
+functions
+{
+//     surfaceFieldValue1
+//     {
+//         type            surfaceFieldValue;
+//         libs            ("libfieldFunctionObjects.so");
+//         enabled         yes;
+//         writeControl    writeTime;
+//         log             yes;
+//         writeFields     no;
+//         regionType      patch;
+//         name            outlet;
+//         operation       weightedAverage;
+//         weightField     phi;
+//         fields
+//         (
+//             H2O
+//             T
+//         );
+//     }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/reactingHeterogenousParcelFoam/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           ( );
+
+// ************************************************************************* //

tutorials/lagrangian/reactingHeterogenousParcelFoam/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;
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/reactingHeterogenousParcelFoam/system/fvSolution

@@ -0,0 +1,111 @@
+/*--------------------------------*- 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;
+    }
+}
+
+
+// ************************************************************************* //