From 25575df60ee54320c36f0a938df52e8f5220eeaf Mon Sep 17 00:00:00 2001
From: sergio <sergio>
Date: Fri, 4 Jan 2019 17:08:53 -0800
Subject: [PATCH] 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
---
 .../reactingParcelFoam/createClouds.H         |   2 +-
 .../reactingHeterogenousParcelFoam/Make/files |   3 +
 .../Make/options                              |  48 ++
 .../reactingHeterogenousParcelFoam.C          |  43 ++
 .../reactingParcelFoam/reactingParcelFoam.C   |  13 +-
 .../COxidationDiffusionLimitedRate.C          |  12 +-
 .../COxidationDiffusionLimitedRate.H          |   8 +-
 .../COxidationHurtMitchell.C                  |  11 +-
 .../COxidationHurtMitchell.H                  |  10 +-
 .../COxidationIntrinsicRate.C                 |  12 +-
 .../COxidationIntrinsicRate.H                 |  10 +-
 .../COxidationKineticDiffusionLimitedRate.C   |  12 +-
 .../COxidationKineticDiffusionLimitedRate.H   |   6 +-
 .../COxidationMurphyShaddix.C                 |  11 +-
 .../COxidationMurphyShaddix.H                 |   6 +-
 src/lagrangian/intermediate/Make/files        |   7 +
 .../Templates/KinematicCloud/KinematicCloud.C |  30 +-
 .../Templates/KinematicCloud/KinematicCloud.H |   3 +
 .../KinematicCloud/KinematicCloudI.H          |  17 +-
 .../Templates/ReactingCloud/ReactingCloud.C   |   4 +-
 .../ReactingHeterogeneousCloud.C              | 248 ++++++++++
 .../ReactingHeterogeneousCloud.H              | 272 +++++++++++
 .../ReactingHeterogeneousCloudI.H             |  60 +++
 .../reactingHeterogeneousCloud.C              |  36 ++
 .../reactingHeterogeneousCloud.H              |  76 ++++
 .../basicHeterogeneousReactingCloud.H         |  68 +++
 .../ReactingHeterogeneousParcel.C             | 397 ++++++++++++++++
 .../ReactingHeterogeneousParcel.H             | 426 ++++++++++++++++++
 .../ReactingHeterogeneousParcelI.H            | 189 ++++++++
 .../ReactingHeterogeneousParcelIO.C           | 354 +++++++++++++++
 .../ReactingMultiphaseParcel.C                |   8 +-
 .../ReactingMultiphaseParcel.H                |  15 +-
 .../ReactingMultiphaseParcelIO.C              |  27 +-
 .../basicHeterogeneousReactingParcel.H        |  74 +++
 .../defineBasicHeterogeneousReactingParcel.C  |  38 ++
 ...asicHeterogeneousReactingParcelSubmodels.C |  69 +++
 ...eactingParcelHeterogeneousReactingModels.H |  47 ++
 .../HeterogeneousReactingModel.C              | 107 +++++
 .../HeterogeneousReactingModel.H              | 218 +++++++++
 .../HeterogeneousReactingModelNew.C           |  61 +++
 .../MUCSheterogeneousRate.C                   | 239 ++++++++++
 .../MUCSheterogeneousRate.H                   | 198 ++++++++
 .../NoheterogeneousReacting.C                 |  92 ++++
 .../NoheterogeneousReacting.H                 | 131 ++++++
 .../NoSurfaceReaction/NoSurfaceReaction.C     |  11 +-
 .../NoSurfaceReaction/NoSurfaceReaction.H     |   6 +-
 .../SurfaceReactionModel.C                    |  13 +-
 .../SurfaceReactionModel.H                    |   8 +-
 src/lagrangian/turbulence/Make/files          |   3 +
 .../makeBasicHeterogeneousParcelSubmodels.C   |  40 ++
 .../solidProperties/C/C.C                     |   4 +-
 .../solidProperties/CaCO3/CaCO3.C             |   4 +-
 .../solidProperties/ash/ash.C                 |   4 +-
 .../solidProperties/solidProperties.C         |  16 +-
 .../solidProperties/solidProperties.H         |  13 +-
 .../solidProperties/solidPropertiesI.H        |   7 +-
 .../reactingHeterogenousParcelFoam/0/N2       |  42 ++
 .../reactingHeterogenousParcelFoam/0/O2       |  42 ++
 .../reactingHeterogenousParcelFoam/0/T        |  41 ++
 .../reactingHeterogenousParcelFoam/0/U        |  43 ++
 .../reactingHeterogenousParcelFoam/0/alphat   |  43 ++
 .../reactingHeterogenousParcelFoam/0/k        |  43 ++
 .../reactingHeterogenousParcelFoam/0/nut      |  46 ++
 .../reactingHeterogenousParcelFoam/0/omega    |  48 ++
 .../reactingHeterogenousParcelFoam/0/p        |  43 ++
 .../reactingHeterogenousParcelFoam/0/p_rgh    |  40 ++
 .../constant/chemistryProperties              |  28 ++
 .../constant/combustionProperties             |  21 +
 .../constant/foam.dat                         |  82 ++++
 .../constant/foam.inp                         |   9 +
 .../reactingHeterogenousParcelFoam/constant/g |  22 +
 .../constant/particleTrackDict                |  28 ++
 .../constant/radiationProperties              |  22 +
 .../constant/reactingCloud1Properties         | 184 ++++++++
 .../constant/thermophysicalProperties         |  73 +++
 .../constant/turbulenceProperties             |  29 ++
 .../system/blockMeshDict                      |  75 +++
 .../system/controlDict                        |  81 ++++
 .../system/decomposeParDict                   |  30 ++
 .../system/fvSchemes                          |  63 +++
 .../system/fvSolution                         | 111 +++++
 81 files changed, 5002 insertions(+), 134 deletions(-)
 create mode 100644 applications/solvers/lagrangian/reactingParcelFoam/reactingHeterogenousParcelFoam/Make/files
 create mode 100644 applications/solvers/lagrangian/reactingParcelFoam/reactingHeterogenousParcelFoam/Make/options
 create mode 100644 applications/solvers/lagrangian/reactingParcelFoam/reactingHeterogenousParcelFoam/reactingHeterogenousParcelFoam.C
 create mode 100644 src/lagrangian/intermediate/clouds/Templates/ReactingHeterogeneousCloud/ReactingHeterogeneousCloud.C
 create mode 100644 src/lagrangian/intermediate/clouds/Templates/ReactingHeterogeneousCloud/ReactingHeterogeneousCloud.H
 create mode 100644 src/lagrangian/intermediate/clouds/Templates/ReactingHeterogeneousCloud/ReactingHeterogeneousCloudI.H
 create mode 100644 src/lagrangian/intermediate/clouds/baseClasses/reactingHeterogeneousCloud/reactingHeterogeneousCloud.C
 create mode 100644 src/lagrangian/intermediate/clouds/baseClasses/reactingHeterogeneousCloud/reactingHeterogeneousCloud.H
 create mode 100644 src/lagrangian/intermediate/clouds/derived/basicHeterogeneousReactingCloud/basicHeterogeneousReactingCloud.H
 create mode 100644 src/lagrangian/intermediate/parcels/Templates/ReactingHeterogeneousParcel/ReactingHeterogeneousParcel.C
 create mode 100644 src/lagrangian/intermediate/parcels/Templates/ReactingHeterogeneousParcel/ReactingHeterogeneousParcel.H
 create mode 100644 src/lagrangian/intermediate/parcels/Templates/ReactingHeterogeneousParcel/ReactingHeterogeneousParcelI.H
 create mode 100644 src/lagrangian/intermediate/parcels/Templates/ReactingHeterogeneousParcel/ReactingHeterogeneousParcelIO.C
 create mode 100644 src/lagrangian/intermediate/parcels/derived/basicHeterogeneousReactingParcel/basicHeterogeneousReactingParcel.H
 create mode 100644 src/lagrangian/intermediate/parcels/derived/basicHeterogeneousReactingParcel/defineBasicHeterogeneousReactingParcel.C
 create mode 100644 src/lagrangian/intermediate/parcels/derived/basicHeterogeneousReactingParcel/makeBasicHeterogeneousReactingParcelSubmodels.C
 create mode 100644 src/lagrangian/intermediate/parcels/include/makeHeterogeneousReactingParcelHeterogeneousReactingModels.H
 create mode 100644 src/lagrangian/intermediate/submodels/HeterogeneousReactingModel/HeterogeneousReactingModel/HeterogeneousReactingModel.C
 create mode 100644 src/lagrangian/intermediate/submodels/HeterogeneousReactingModel/HeterogeneousReactingModel/HeterogeneousReactingModel.H
 create mode 100644 src/lagrangian/intermediate/submodels/HeterogeneousReactingModel/HeterogeneousReactingModel/HeterogeneousReactingModelNew.C
 create mode 100644 src/lagrangian/intermediate/submodels/HeterogeneousReactingModel/MUCSheterogeneousRate/MUCSheterogeneousRate.C
 create mode 100644 src/lagrangian/intermediate/submodels/HeterogeneousReactingModel/MUCSheterogeneousRate/MUCSheterogeneousRate.H
 create mode 100644 src/lagrangian/intermediate/submodels/HeterogeneousReactingModel/NoheterogeneousReacting/NoheterogeneousReacting.C
 create mode 100644 src/lagrangian/intermediate/submodels/HeterogeneousReactingModel/NoheterogeneousReacting/NoheterogeneousReacting.H
 create mode 100644 src/lagrangian/turbulence/parcels/derived/basicHeterogeneousReactingParcel/makeBasicHeterogeneousParcelSubmodels.C
 create mode 100644 tutorials/lagrangian/reactingHeterogenousParcelFoam/0/N2
 create mode 100644 tutorials/lagrangian/reactingHeterogenousParcelFoam/0/O2
 create mode 100644 tutorials/lagrangian/reactingHeterogenousParcelFoam/0/T
 create mode 100644 tutorials/lagrangian/reactingHeterogenousParcelFoam/0/U
 create mode 100644 tutorials/lagrangian/reactingHeterogenousParcelFoam/0/alphat
 create mode 100644 tutorials/lagrangian/reactingHeterogenousParcelFoam/0/k
 create mode 100644 tutorials/lagrangian/reactingHeterogenousParcelFoam/0/nut
 create mode 100644 tutorials/lagrangian/reactingHeterogenousParcelFoam/0/omega
 create mode 100644 tutorials/lagrangian/reactingHeterogenousParcelFoam/0/p
 create mode 100644 tutorials/lagrangian/reactingHeterogenousParcelFoam/0/p_rgh
 create mode 100644 tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/chemistryProperties
 create mode 100644 tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/combustionProperties
 create mode 100644 tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/foam.dat
 create mode 100644 tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/foam.inp
 create mode 100644 tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/g
 create mode 100644 tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/particleTrackDict
 create mode 100644 tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/radiationProperties
 create mode 100644 tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/reactingCloud1Properties
 create mode 100644 tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/thermophysicalProperties
 create mode 100644 tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/turbulenceProperties
 create mode 100644 tutorials/lagrangian/reactingHeterogenousParcelFoam/system/blockMeshDict
 create mode 100644 tutorials/lagrangian/reactingHeterogenousParcelFoam/system/controlDict
 create mode 100644 tutorials/lagrangian/reactingHeterogenousParcelFoam/system/decomposeParDict
 create mode 100644 tutorials/lagrangian/reactingHeterogenousParcelFoam/system/fvSchemes
 create mode 100644 tutorials/lagrangian/reactingHeterogenousParcelFoam/system/fvSolution

diff --git a/applications/solvers/lagrangian/reactingParcelFoam/createClouds.H b/applications/solvers/lagrangian/reactingParcelFoam/createClouds.H
index 954b74e069..ceaca29f59 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/createClouds.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/createClouds.H
@@ -1,5 +1,5 @@
 Info<< "\nConstructing reacting cloud" << endl;
-basicReactingMultiphaseCloud parcels
+basicReactingTypeCloud parcels
 (
     "reactingCloud1",
     rho,
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/reactingHeterogenousParcelFoam/Make/files b/applications/solvers/lagrangian/reactingParcelFoam/reactingHeterogenousParcelFoam/Make/files
new file mode 100644
index 0000000000..c140bdda79
--- /dev/null
+++ b/applications/solvers/lagrangian/reactingParcelFoam/reactingHeterogenousParcelFoam/Make/files
@@ -0,0 +1,3 @@
+reactingHeterogenousParcelFoam.C
+
+EXE = $(FOAM_APPBIN)/reactingHeterogenousParcelFoam
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/reactingHeterogenousParcelFoam/Make/options b/applications/solvers/lagrangian/reactingParcelFoam/reactingHeterogenousParcelFoam/Make/options
new file mode 100644
index 0000000000..841ef8c9c7
--- /dev/null
+++ b/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
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/reactingHeterogenousParcelFoam/reactingHeterogenousParcelFoam.C b/applications/solvers/lagrangian/reactingParcelFoam/reactingHeterogenousParcelFoam/reactingHeterogenousParcelFoam.C
new file mode 100644
index 0000000000..ee4e7ed4ca
--- /dev/null
+++ b/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"
+
+// ************************************************************************* //
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
index 92623f9ce3..f6609a1583 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
+++ b/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)
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationDiffusionLimitedRate/COxidationDiffusionLimitedRate.C b/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationDiffusionLimitedRate/COxidationDiffusionLimitedRate.C
index 5a2d0c0688..c51e973eba 100644
--- a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationDiffusionLimitedRate/COxidationDiffusionLimitedRate.C
+++ b/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,
diff --git a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationDiffusionLimitedRate/COxidationDiffusionLimitedRate.H b/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationDiffusionLimitedRate/COxidationDiffusionLimitedRate.H
index 7fa288b87b..3329e9c29e 100644
--- a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationDiffusionLimitedRate/COxidationDiffusionLimitedRate.H
+++ b/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,
diff --git a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationHurtMitchell/COxidationHurtMitchell.C b/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationHurtMitchell/COxidationHurtMitchell.C
index d4250fe65e..2dbc1b482e 100644
--- a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationHurtMitchell/COxidationHurtMitchell.C
+++ b/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,
diff --git a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationHurtMitchell/COxidationHurtMitchell.H b/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationHurtMitchell/COxidationHurtMitchell.H
index 15dbfd9d5b..a8a45f231c 100644
--- a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationHurtMitchell/COxidationHurtMitchell.H
+++ b/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,
diff --git a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationIntrinsicRate/COxidationIntrinsicRate.C b/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationIntrinsicRate/COxidationIntrinsicRate.C
index 4a6cad7d77..8ceb9441f1 100644
--- a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationIntrinsicRate/COxidationIntrinsicRate.C
+++ b/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,
diff --git a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationIntrinsicRate/COxidationIntrinsicRate.H b/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationIntrinsicRate/COxidationIntrinsicRate.H
index fff61aa539..84aa7b7818 100644
--- a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationIntrinsicRate/COxidationIntrinsicRate.H
+++ b/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,
diff --git a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationKineticDiffusionLimitedRate/COxidationKineticDiffusionLimitedRate.C b/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationKineticDiffusionLimitedRate/COxidationKineticDiffusionLimitedRate.C
index 73f1000e5d..d5f66d8754 100644
--- a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationKineticDiffusionLimitedRate/COxidationKineticDiffusionLimitedRate.C
+++ b/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,
diff --git a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationKineticDiffusionLimitedRate/COxidationKineticDiffusionLimitedRate.H b/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationKineticDiffusionLimitedRate/COxidationKineticDiffusionLimitedRate.H
index 485a6244ed..6e7ac7ee0c 100644
--- a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationKineticDiffusionLimitedRate/COxidationKineticDiffusionLimitedRate.H
+++ b/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,
diff --git a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationMurphyShaddix/COxidationMurphyShaddix.C b/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationMurphyShaddix/COxidationMurphyShaddix.C
index 5889f416ab..299b0ccae2 100644
--- a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationMurphyShaddix/COxidationMurphyShaddix.C
+++ b/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,
diff --git a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationMurphyShaddix/COxidationMurphyShaddix.H b/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationMurphyShaddix/COxidationMurphyShaddix.H
index 09fee068be..403a03b1b1 100644
--- a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationMurphyShaddix/COxidationMurphyShaddix.H
+++ b/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,
diff --git a/src/lagrangian/intermediate/Make/files b/src/lagrangian/intermediate/Make/files
index 8789d3c0bd..4d6e12f6dc 100644
--- a/src/lagrangian/intermediate/Make/files
+++ b/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
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
index 3c36b4b3e4..f5e5d96739 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
+++ b/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();
-    reduce(linearMomentum, sumOp<vector>());
+    const vector linearMomentum =
+        returnReduce(linearMomentumOfSystem(), sumOp<vector>());
 
-    scalar linearKineticEnergy = linearKineticEnergyOfSystem();
-    reduce(linearKineticEnergy, sumOp<scalar>());
+    const scalar linearKineticEnergy =
+        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       = "
-        << returnReduce(this->size(), sumOp<label>()) << nl
+        << "    Current number of parcels       = " << nTotParcel << nl
         << "    Current mass in system          = "
         << returnReduce(massInSystem(), sumOp<scalar>()) << nl
-        << "    Linear momentum                 = "
-        << linearMomentum << nl
-        << "   |Linear momentum|                = "
-        << mag(linearMomentum) << nl
-        << "    Linear kinetic energy           = "
-        << linearKineticEnergy << nl;
+        << "    Linear momentum                 = " << linearMomentum << nl
+        << "   |Linear momentum|                = " << mag(linearMomentum) << nl
+        << "    Linear kinetic energy           = " << linearKineticEnergy << nl
+        << "    Average particle per parcel     = " << particlePerParcel << nl;
 
     injectors_.info(Info);
     this->surfaceFilm().info(Info);
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
index 99af097272..35af3acf59 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
+++ b/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;
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
index b10226c31c..095eabd9c6 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
+++ b/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)
     {
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
index 2cfd0a27ff..018a7ab7b7 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
+++ b/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)
 {}
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingHeterogeneousCloud/ReactingHeterogeneousCloud.C b/src/lagrangian/intermediate/clouds/Templates/ReactingHeterogeneousCloud/ReactingHeterogeneousCloud.C
new file mode 100644
index 0000000000..9069b759e0
--- /dev/null
+++ b/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);
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingHeterogeneousCloud/ReactingHeterogeneousCloud.H b/src/lagrangian/intermediate/clouds/Templates/ReactingHeterogeneousCloud/ReactingHeterogeneousCloud.H
new file mode 100644
index 0000000000..c0e61c2795
--- /dev/null
+++ b/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
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingHeterogeneousCloud/ReactingHeterogeneousCloudI.H b/src/lagrangian/intermediate/clouds/Templates/ReactingHeterogeneousCloud/ReactingHeterogeneousCloudI.H
new file mode 100644
index 0000000000..c9077e909b
--- /dev/null
+++ b/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();
+}
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/clouds/baseClasses/reactingHeterogeneousCloud/reactingHeterogeneousCloud.C b/src/lagrangian/intermediate/clouds/baseClasses/reactingHeterogeneousCloud/reactingHeterogeneousCloud.C
new file mode 100644
index 0000000000..b7ef77e4a3
--- /dev/null
+++ b/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);
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/clouds/baseClasses/reactingHeterogeneousCloud/reactingHeterogeneousCloud.H b/src/lagrangian/intermediate/clouds/baseClasses/reactingHeterogeneousCloud/reactingHeterogeneousCloud.H
new file mode 100644
index 0000000000..ec7d93013e
--- /dev/null
+++ b/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
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/clouds/derived/basicHeterogeneousReactingCloud/basicHeterogeneousReactingCloud.H b/src/lagrangian/intermediate/clouds/derived/basicHeterogeneousReactingCloud/basicHeterogeneousReactingCloud.H
new file mode 100644
index 0000000000..6c723e6893
--- /dev/null
+++ b/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
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingHeterogeneousParcel/ReactingHeterogeneousParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingHeterogeneousParcel/ReactingHeterogeneousParcel.C
new file mode 100644
index 0000000000..1f8c87e802
--- /dev/null
+++ b/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"
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingHeterogeneousParcel/ReactingHeterogeneousParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingHeterogeneousParcel/ReactingHeterogeneousParcel.H
new file mode 100644
index 0000000000..7e92e9d01d
--- /dev/null
+++ b/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
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingHeterogeneousParcel/ReactingHeterogeneousParcelI.H b/src/lagrangian/intermediate/parcels/Templates/ReactingHeterogeneousParcel/ReactingHeterogeneousParcelI.H
new file mode 100644
index 0000000000..4eeac61669
--- /dev/null
+++ b/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_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingHeterogeneousParcel/ReactingHeterogeneousParcelIO.C b/src/lagrangian/intermediate/parcels/Templates/ReactingHeterogeneousParcel/ReactingHeterogeneousParcelIO.C
new file mode 100644
index 0000000000..3465116bb8
--- /dev/null
+++ b/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;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
index 33f3412066..e25e197209 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
+++ b/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,
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
index d62d9d8a5e..99d6cda1ca 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
+++ b/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>
-class ReactingMultiphaseParcel;
+template<class ParcelType> class ReactingMultiphaseParcel;
 
 template<class ParcelType>
-Ostream& operator<<
-(
-    Ostream&,
-    const ReactingMultiphaseParcel<ParcelType>&
-);
+Ostream& operator<<(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
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelIO.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelIO.C
index 3701ed26b7..aa75ca3794 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelIO.C
+++ b/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;
             }
         }
     }
diff --git a/src/lagrangian/intermediate/parcels/derived/basicHeterogeneousReactingParcel/basicHeterogeneousReactingParcel.H b/src/lagrangian/intermediate/parcels/derived/basicHeterogeneousReactingParcel/basicHeterogeneousReactingParcel.H
new file mode 100644
index 0000000000..63b9156fdb
--- /dev/null
+++ b/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
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/parcels/derived/basicHeterogeneousReactingParcel/defineBasicHeterogeneousReactingParcel.C b/src/lagrangian/intermediate/parcels/derived/basicHeterogeneousReactingParcel/defineBasicHeterogeneousReactingParcel.C
new file mode 100644
index 0000000000..582f8290a7
--- /dev/null
+++ b/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);
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/parcels/derived/basicHeterogeneousReactingParcel/makeBasicHeterogeneousReactingParcelSubmodels.C b/src/lagrangian/intermediate/parcels/derived/basicHeterogeneousReactingParcel/makeBasicHeterogeneousReactingParcelSubmodels.C
new file mode 100644
index 0000000000..5214727a14
--- /dev/null
+++ b/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
+);
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/parcels/include/makeHeterogeneousReactingParcelHeterogeneousReactingModels.H b/src/lagrangian/intermediate/parcels/include/makeHeterogeneousReactingParcelHeterogeneousReactingModels.H
new file mode 100644
index 0000000000..a2b57b0bfc
--- /dev/null
+++ b/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
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/HeterogeneousReactingModel/HeterogeneousReactingModel/HeterogeneousReactingModel.C b/src/lagrangian/intermediate/submodels/HeterogeneousReactingModel/HeterogeneousReactingModel/HeterogeneousReactingModel.C
new file mode 100644
index 0000000000..1dc780b82a
--- /dev/null
+++ b/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"
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/HeterogeneousReactingModel/HeterogeneousReactingModel/HeterogeneousReactingModel.H b/src/lagrangian/intermediate/submodels/HeterogeneousReactingModel/HeterogeneousReactingModel/HeterogeneousReactingModel.H
new file mode 100644
index 0000000000..c90ff27dcc
--- /dev/null
+++ b/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
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/HeterogeneousReactingModel/HeterogeneousReactingModel/HeterogeneousReactingModelNew.C b/src/lagrangian/intermediate/submodels/HeterogeneousReactingModel/HeterogeneousReactingModel/HeterogeneousReactingModelNew.C
new file mode 100644
index 0000000000..9f326da207
--- /dev/null
+++ b/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)
+    );
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/HeterogeneousReactingModel/MUCSheterogeneousRate/MUCSheterogeneousRate.C b/src/lagrangian/intermediate/submodels/HeterogeneousReactingModel/MUCSheterogeneousRate/MUCSheterogeneousRate.C
new file mode 100644
index 0000000000..a8266d9d11
--- /dev/null
+++ b/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;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/HeterogeneousReactingModel/MUCSheterogeneousRate/MUCSheterogeneousRate.H b/src/lagrangian/intermediate/submodels/HeterogeneousReactingModel/MUCSheterogeneousRate/MUCSheterogeneousRate.H
new file mode 100644
index 0000000000..395169ba52
--- /dev/null
+++ b/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
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/HeterogeneousReactingModel/NoheterogeneousReacting/NoheterogeneousReacting.C b/src/lagrangian/intermediate/submodels/HeterogeneousReactingModel/NoheterogeneousReacting/NoheterogeneousReacting.C
new file mode 100644
index 0000000000..12d803b6c7
--- /dev/null
+++ b/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;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/HeterogeneousReactingModel/NoheterogeneousReacting/NoheterogeneousReacting.H b/src/lagrangian/intermediate/submodels/HeterogeneousReactingModel/NoheterogeneousReacting/NoheterogeneousReacting.H
new file mode 100644
index 0000000000..8c5a62465e
--- /dev/null
+++ b/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
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.C b/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.C
index 901e75b77f..e14119fc0e 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.C
+++ b/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,
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.H b/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.H
index 371ae9c464..5b70775df1 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.H
+++ b/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,
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/SurfaceReactionModel/SurfaceReactionModel.C b/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/SurfaceReactionModel/SurfaceReactionModel.C
index ae5ef36911..dc3a7b99bb 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/SurfaceReactionModel/SurfaceReactionModel.C
+++ b/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>
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/SurfaceReactionModel/SurfaceReactionModel.H b/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/SurfaceReactionModel/SurfaceReactionModel.H
index c261bea5ef..7f86a856aa 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/SurfaceReactionModel/SurfaceReactionModel.H
+++ b/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,
diff --git a/src/lagrangian/turbulence/Make/files b/src/lagrangian/turbulence/Make/files
index 730a4c9fff..a4277b430f 100644
--- a/src/lagrangian/turbulence/Make/files
+++ b/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
diff --git a/src/lagrangian/turbulence/parcels/derived/basicHeterogeneousReactingParcel/makeBasicHeterogeneousParcelSubmodels.C b/src/lagrangian/turbulence/parcels/derived/basicHeterogeneousReactingParcel/makeBasicHeterogeneousParcelSubmodels.C
new file mode 100644
index 0000000000..eeecfac240
--- /dev/null
+++ b/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);
+}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/thermophysicalProperties/solidProperties/C/C.C b/src/thermophysicalModels/thermophysicalProperties/solidProperties/C/C.C
index 9898750c20..d3b2529f5d 100644
--- a/src/thermophysicalModels/thermophysicalProperties/solidProperties/C/C.C
+++ b/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)
     {
diff --git a/src/thermophysicalModels/thermophysicalProperties/solidProperties/CaCO3/CaCO3.C b/src/thermophysicalModels/thermophysicalProperties/solidProperties/CaCO3/CaCO3.C
index b8b86dae87..8151d4e09a 100644
--- a/src/thermophysicalModels/thermophysicalProperties/solidProperties/CaCO3/CaCO3.C
+++ b/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)
     {
diff --git a/src/thermophysicalModels/thermophysicalProperties/solidProperties/ash/ash.C b/src/thermophysicalModels/thermophysicalProperties/solidProperties/ash/ash.C
index e74e01e37b..9a5f3b2b97 100644
--- a/src/thermophysicalModels/thermophysicalProperties/solidProperties/ash/ash.C
+++ b/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)
     {
diff --git a/src/thermophysicalModels/thermophysicalProperties/solidProperties/solidProperties/solidProperties.C b/src/thermophysicalModels/thermophysicalProperties/solidProperties/solidProperties/solidProperties.C
index 475d52d705..c8716eee38 100644
--- a/src/thermophysicalModels/thermophysicalProperties/solidProperties/solidProperties/solidProperties.C
+++ b/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_;
 }
 
 
diff --git a/src/thermophysicalModels/thermophysicalProperties/solidProperties/solidProperties/solidProperties.H b/src/thermophysicalModels/thermophysicalProperties/solidProperties/solidProperties/solidProperties.H
index 36adb9c6b6..7e8c204142 100644
--- a/src/thermophysicalModels/thermophysicalProperties/solidProperties/solidProperties/solidProperties.H
+++ b/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
 
diff --git a/src/thermophysicalModels/thermophysicalProperties/solidProperties/solidProperties/solidPropertiesI.H b/src/thermophysicalModels/thermophysicalProperties/solidProperties/solidProperties/solidPropertiesI.H
index f824fc6b55..7e033307f0 100644
--- a/src/thermophysicalModels/thermophysicalProperties/solidProperties/solidProperties/solidPropertiesI.H
+++ b/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_;
+}
+
 
 // ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingHeterogenousParcelFoam/0/N2 b/tutorials/lagrangian/reactingHeterogenousParcelFoam/0/N2
new file mode 100644
index 0000000000..7282998db7
--- /dev/null
+++ b/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;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingHeterogenousParcelFoam/0/O2 b/tutorials/lagrangian/reactingHeterogenousParcelFoam/0/O2
new file mode 100644
index 0000000000..f24eb76861
--- /dev/null
+++ b/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;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingHeterogenousParcelFoam/0/T b/tutorials/lagrangian/reactingHeterogenousParcelFoam/0/T
new file mode 100644
index 0000000000..97c0dcd9a7
--- /dev/null
+++ b/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;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingHeterogenousParcelFoam/0/U b/tutorials/lagrangian/reactingHeterogenousParcelFoam/0/U
new file mode 100644
index 0000000000..d5cc789130
--- /dev/null
+++ b/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);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingHeterogenousParcelFoam/0/alphat b/tutorials/lagrangian/reactingHeterogenousParcelFoam/0/alphat
new file mode 100644
index 0000000000..8c3ac9e91e
--- /dev/null
+++ b/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;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingHeterogenousParcelFoam/0/k b/tutorials/lagrangian/reactingHeterogenousParcelFoam/0/k
new file mode 100644
index 0000000000..4189f4324b
--- /dev/null
+++ b/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;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingHeterogenousParcelFoam/0/nut b/tutorials/lagrangian/reactingHeterogenousParcelFoam/0/nut
new file mode 100644
index 0000000000..af9bc49294
--- /dev/null
+++ b/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;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingHeterogenousParcelFoam/0/omega b/tutorials/lagrangian/reactingHeterogenousParcelFoam/0/omega
new file mode 100644
index 0000000000..1c72607137
--- /dev/null
+++ b/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;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingHeterogenousParcelFoam/0/p b/tutorials/lagrangian/reactingHeterogenousParcelFoam/0/p
new file mode 100644
index 0000000000..a924bbd16e
--- /dev/null
+++ b/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;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingHeterogenousParcelFoam/0/p_rgh b/tutorials/lagrangian/reactingHeterogenousParcelFoam/0/p_rgh
new file mode 100644
index 0000000000..19370a271e
--- /dev/null
+++ b/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;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/chemistryProperties b/tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/chemistryProperties
new file mode 100644
index 0000000000..e4a6386238
--- /dev/null
+++ b/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
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/combustionProperties b/tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/combustionProperties
new file mode 100644
index 0000000000..f49b11d2e3
--- /dev/null
+++ b/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;
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/foam.dat b/tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/foam.dat
new file mode 100644
index 0000000000..f3e7973bfc
--- /dev/null
+++ b/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;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/foam.inp b/tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/foam.inp
new file mode 100644
index 0000000000..197c73031a
--- /dev/null
+++ b/tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/foam.inp
@@ -0,0 +1,9 @@
+species
+(
+    O2
+    N2
+)
+;
+
+reactions
+{}
diff --git a/tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/g b/tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/g
new file mode 100644
index 0000000000..037697dd98
--- /dev/null
+++ b/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);
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/particleTrackDict b/tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/particleTrackDict
new file mode 100644
index 0000000000..b523c37e92
--- /dev/null
+++ b/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
+);
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/radiationProperties b/tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/radiationProperties
new file mode 100644
index 0000000000..cb46544fd0
--- /dev/null
+++ b/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;
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/reactingCloud1Properties b/tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/reactingCloud1Properties
new file mode 100644
index 0000000000..bb4eeb5d9b
--- /dev/null
+++ b/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
+{
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/thermophysicalProperties b/tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/thermophysicalProperties
new file mode 100644
index 0000000000..0870029d84
--- /dev/null
+++ b/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
+        }
+    }
+}
+
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/turbulenceProperties b/tutorials/lagrangian/reactingHeterogenousParcelFoam/constant/turbulenceProperties
new file mode 100644
index 0000000000..1433192686
--- /dev/null
+++ b/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;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingHeterogenousParcelFoam/system/blockMeshDict b/tutorials/lagrangian/reactingHeterogenousParcelFoam/system/blockMeshDict
new file mode 100644
index 0000000000..068a3e1e1a
--- /dev/null
+++ b/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)
+        );
+    }
+);
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingHeterogenousParcelFoam/system/controlDict b/tutorials/lagrangian/reactingHeterogenousParcelFoam/system/controlDict
new file mode 100644
index 0000000000..bf221166de
--- /dev/null
+++ b/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
+//         );
+//     }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingHeterogenousParcelFoam/system/decomposeParDict b/tutorials/lagrangian/reactingHeterogenousParcelFoam/system/decomposeParDict
new file mode 100644
index 0000000000..6cad762db7
--- /dev/null
+++ b/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           ( );
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingHeterogenousParcelFoam/system/fvSchemes b/tutorials/lagrangian/reactingHeterogenousParcelFoam/system/fvSchemes
new file mode 100644
index 0000000000..16f8cdf93a
--- /dev/null
+++ b/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;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingHeterogenousParcelFoam/system/fvSolution b/tutorials/lagrangian/reactingHeterogenousParcelFoam/system/fvSolution
new file mode 100644
index 0000000000..2949998aca
--- /dev/null
+++ b/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;
+    }
+}
+
+
+// ************************************************************************* //