From c5c16150b7319adbf7265959bdf56080348dbee2 Mon Sep 17 00:00:00 2001
From: andy <a.heather@opencfd.co.uk>
Date: Mon, 16 Mar 2009 17:43:12 +0000
Subject: [PATCH] multiple updates

---
 .../Templates/KinematicCloud/KinematicCloud.C |  23 +--
 .../Templates/KinematicCloud/KinematicCloud.H |   3 -
 .../Templates/ReactingCloud/ReactingCloud.C   |   2 +-
 .../ReactingMultiphaseCloud.C                 |   2 +-
 .../Templates/ThermoCloud/ThermoCloud.C       |   2 +-
 .../KinematicParcel/KinematicParcel.C         |  53 +++---
 .../KinematicParcel/KinematicParcel.H         |   3 +-
 .../ReactingMultiphaseParcel.C                | 164 ++++++++++--------
 .../ReactingMultiphaseParcel.H                |  16 +-
 .../Templates/ReactingParcel/ReactingParcel.C |  83 ++++-----
 .../Templates/ReactingParcel/ReactingParcel.H |   2 +-
 .../Templates/ThermoParcel/ThermoParcel.C     |  59 ++++---
 .../LiquidEvaporation/LiquidEvaporation.C     |   3 -
 .../ConstantRateDevolatilisation.C            |   4 +-
 .../ConstantRateDevolatilisation.H            |   4 +-
 .../DevolatilisationModel.H                   |   4 +-
 .../NoDevolatilisation/NoDevolatilisation.C   |   4 +-
 .../NoDevolatilisation/NoDevolatilisation.H   |   4 +-
 .../SingleKineticRateDevolatilisation.C       |   4 +-
 .../SingleKineticRateDevolatilisation.H       |   4 +-
 .../NoSurfaceReaction/NoSurfaceReaction.C     |  33 ++--
 .../NoSurfaceReaction/NoSurfaceReaction.H     |  14 +-
 .../SurfaceReactionModel.H                    |  15 +-
 23 files changed, 258 insertions(+), 247 deletions(-)

diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
index 6df424c523a..f60f69261ce 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
@@ -232,7 +232,7 @@ void Foam::KinematicCloud<ParcelType>::evolve()
 
     if (debug)
     {
-        this->dumpParticlePositions();
+        this->writePositions();
     }
 
     if (coupled_)
@@ -261,25 +261,4 @@ void Foam::KinematicCloud<ParcelType>::info() const
 }
 
 
-template<class ParcelType>
-void Foam::KinematicCloud<ParcelType>::dumpParticlePositions() const
-{
-    OFstream pObj
-    (
-        this->db().time().path()/"parcelPositions_"
-      + this->name() + "_"
-      + name(this->injection().nInjections()) + ".obj"
-    );
-
-    forAllConstIter(typename KinematicCloud<ParcelType>, *this, iter)
-    {
-        const ParcelType& p = iter();
-        pObj<< "v " << p.position().x() << " " << p.position().y() << " "
-            << p.position().z() << nl;
-    }
-
-    pObj.flush();
-}
-
-
 // ************************************************************************* //
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
index 77a61b02d29..7f39808f83b 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
@@ -338,9 +338,6 @@ public:
             //- Print cloud information
             void info() const;
 
-            //- Dump particle positions to .obj file
-            void dumpParticlePositions() const;
-
 
             // Fields
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
index 15a088c000b..b14fcc17cd8 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
@@ -204,7 +204,7 @@ void Foam::ReactingCloud<ParcelType>::evolve()
 
     if (debug)
     {
-        this->dumpParticlePositions();
+        this->writePositions();
     }
 
     if (this->coupled())
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C b/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C
index 708a3f72239..29936f328ce 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C
@@ -181,7 +181,7 @@ void Foam::ReactingMultiphaseCloud<ParcelType>::evolve()
 
     if (debug)
     {
-        this->dumpParticlePositions();
+        this->writePositions();
     }
 
     if (this->coupled())
diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
index b54624bb812..40731ec2545 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
@@ -203,7 +203,7 @@ void Foam::ThermoCloud<ParcelType>::evolve()
 
     if (debug)
     {
-        this->dumpParticlePositions();
+        this->writePositions();
     }
 
     if (this->coupled())
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index c9df5a57768..637bf15fa4c 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -27,8 +27,6 @@ License
 #include "KinematicParcel.H"
 #include "dimensionedConstants.H"
 
-#include "fvcGrad.H"
-
 // * * * * * * * * * * *  Protected Member Functions * * * * * * * * * * * * //
 
 template<class ParcelType>
@@ -66,30 +64,43 @@ void Foam::KinematicParcel<ParcelType>::calc
     const label cellI
 )
 {
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // 1. Calculate velocity - update U, dUTrans
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    scalar Cud = 0.0;
+    // Define local properties at beginning of time step
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    const scalar mass0 = mass();
+
+
+    // Initialise transfer terms
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Momentum
     vector dUTrans = vector::zero;
-    const vector U1 = calcVelocity(td, dt, vector::zero, mass(), Cud, dUTrans);
 
 
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // 2. Accumulate carrier phase source terms
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Motion
+    // ~~~~~~
+
+    // No additional forces
+    vector Fx = vector::zero;
+
+    // Calculate new particle velocity
+    scalar Cud = 0.0;
+    vector U1 = calcVelocity(td, dt, cellI, Fx, mass0, Cud, dUTrans);
+
+
+    // Accumulate carrier phase source terms
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     if (td.cloud().coupled())
     {
         // Update momentum transfer
         td.cloud().UTrans()[cellI] += nParticle_*dUTrans;
 
         // Coefficient to be applied in carrier phase momentum coupling
-        td.cloud().UCoeff()[cellI] += nParticle_*mass()*Cud;
+        td.cloud().UCoeff()[cellI] += nParticle_*mass0*Cud;
     }
 
 
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // 3. Set new particle properties
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Set new particle properties
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
     U_ = U1;
 }
 
@@ -100,10 +111,11 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
 (
     TrackData& td,
     const scalar dt,
+    const label cellI,
     const vector& Fx,
     const scalar mass,
     scalar& Cud,
-    vector& dUTrans,
+    vector& dUTrans
 ) const
 {
     // Return linearised term from drag model
@@ -127,12 +139,13 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
     // Pressure gradient force
     if (td.cloud().forcePressureGradient())
     {
-        Ftot += rhoc_/rho_*(U_ & fvc::grad(Uc_))
+        const vector& d = this->mesh().deltaCoeffs()[cellI];
+        Ftot += rhoc_/rho_*(U_ & (d^Uc_));
     }
 
-    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Set new particle velocity
-    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // New particle velocity
+    //~~~~~~~~~~~~~~~~~~~~~~
 
     // Update velocity - treat as 3-D
     const vector ap = Uc_ + (Ftot + Fx)/(Cud + VSMALL);
@@ -142,8 +155,6 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
 
     // Calculate the momentum transfer to the continuous phase
     // - do not include gravity impulse
-
-    // TODO: USE AVERAGE PARTICLE MASS
     dUTrans = -mass*(Unew - U_ - dt*td.g());
 
     return Unew;
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
index b6df66d414c..35b46329dc6 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
@@ -84,7 +84,6 @@ public:
     //- Class to hold kinematic particle constant properties
     class constantProperties
     {
-
         // Private data
 
             //- Particle density [kg/m3] (constant)
@@ -120,7 +119,6 @@ public:
     :
         public Particle<ParcelType>::trackData
     {
-
         // Private data
 
             //- Reference to the cloud containing this particle
@@ -234,6 +232,7 @@ protected:
         (
             TrackData& td,
             const scalar dt,
+            const label cellI,
             const vector& Fx,
             const scalar mass,
             scalar& Cud,
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
index 816691d7700..8bcc94cef22 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
@@ -84,7 +84,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     const label cellI
 )
 {
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Define local properties at beginning of timestep
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     const scalar mass0 = this->mass();
@@ -92,97 +91,94 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     const scalar T0 = this->T_;
     const scalar pc = this->pc_;
     scalarField& YMix = this->Y_;
+    const label idG = td.cloud().composition().idGas();
     const label idL = td.cloud().composition().idLiquid();
 
 
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Calculate phase change in liquid phase
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Mass transfer from particle to carrier phase
+    // Initialise transfer terms
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Momentum
+    vector dUTrans = vector::zero;
+
+    // Enthalpy
+    scalar dhTrans = 0.0;
+
+    // Mass transfer due to phase change
     scalarList dMassPC(td.cloud().gases().size(), 0.0);
-    scalar shPC = 
+
+    // Mass transfer due to devolatilisation
+    scalarList dMassDV(td.cloud().gases().size(), 0.0);
+
+    // Change in carrier phase composition due to surface reactions
+    scalarList dMassSRc(td.cloud().gases().size(), 0.0);
+
+
+    // Phase change in liquid phase
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Return enthalpy source and calc mass transfer due to phase change
+    scalar ShPC =
         calcPhaseChange(td, dt, cellI, T0, idL, YMix[idL], YLiquid_, dMassPC);
 
     // Update particle component mass fractions
-    updateMassFraction(mass0, dMassPC, YLiquid_);
+    this->updateMassFraction(mass0, dMassPC, YLiquid_);
 
 
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Calculate Devolatilisation
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Mass transfer from particle to carrier phase
-    // - components exist in particle already
-    scalarList dMassDV(td.cloud().gases().size(), 0.0);
-    scalar shDV = calcDevolatilisation(td, dt, T0, mass0, idGas, YMix, dMassDV);
+    // Devolatilisation
+    // ~~~~~~~~~~~~~~~~
 
+    // Return enthalpy source and calc mass transfer due to devolatilisation
+    scalar ShDV = calcDevolatilisation(td, dt, T0, mass0, idG, YMix, dMassDV);
+
+
+    // Surface reactions
+    // ~~~~~~~~~~~~~~~~~
 
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Calculate surface reactions
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Mass transfer of volatile components from particle to carrier phase
     const scalarList dMassMT = dMassPC + dMassDV;
-    // Mass transfer due to surface reactions from particle to carrier phase
-    // - components do not necessarily exist in particle already
-    scalarList dMassSR(td.cloud().gases().size(), 0.0);
-    // Initialise enthalpy retention to zero
-    scalar dhRet = 0.0;
-    calcSurfaceReactions(td, dt, cellI, T0, dMassMT, dMassSR, dhRet);
-
-    // Enthalpy retention divided between particle and carrier phase by the
-    // fraction fh and (1 - fh)
-    scalar ShSR = td.constProps().fh()*dhRet;
-    dhTrans -= (1.0 - td.constProps().fh())*dhRet;
 
-    // Correct dhTrans to account for enthalpy of consumed solids
-    dhTrans +=
-        sum(dMassSR)*td.cloud().composition().H(idSolid, YSolid_, pc, T0);
-
-    // Correct dhTrans to account for enthalpy of evolved volatiles
-    dhTrans +=
-        sum(dMassMT)*td.cloud().composition().H(idGas, YGas_, pc, T0);
+    // Return enthalpy source and calc mass transfer(s) due to surface reaction
+    scalar ShSR =
+        calcSurfaceReactions(td, dt, cellI, T0, dMassMT, dMassSRc, dhTrans);
 
 
-    // ~~~~~~~~~~~~~~~~~~~~~~~
-    // Calculate heat transfer
-    // ~~~~~~~~~~~~~~~~~~~~~~~
-    scalar htc = 0.0;
+    // Heat transfer
+    // ~~~~~~~~~~~~~
 
     // Total enthalpy source
     scalar Sh = ShPC + ShDV + ShSR;
 
-    scalar T1 = calcHeatTransfer(td, dt, cellI, Sh, htc, shHT);
+    // Calculate new particle temperature
+    scalar htc = 0.0;
+    scalar T1 = calcHeatTransfer(td, dt, cellI, Sh, htc, dhTrans);
 
 
-    // ~~~~~~~~~~~~~~~~~~
-    // Calculate velocity
-    // ~~~~~~~~~~~~~~~~~~
-    // Update mass
-    scalar mass1 = mass0 - massPC - massD - massSR;
-    scalar Cud = 0.0;
-    vector dUTrans = vector::zero;
-    vector Fx = vector::zero;
-    vector U1 = calcVelocity(td, dt, Fx, 0.5*(mass0 + mass1), Cud, dUTrans);
+    // Motion
+    // ~~~~~~
 
+    // Update mass (not to include cMassSR)
+    scalar mass1 = mass0 - sum(dMassPC) - sum(dMassDV);
 
+    // No additional forces
+    vector Fx = vector::zero;
 
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Collect contributions to determine new particle thermo properties
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Update specific heat capacity
-    cp1 = cpEff(td, pc, T1);
+    // Calculate new particle velocity
+    scalar Cud = 0.0;
+    vector U1 =
+        calcVelocity(td, dt, cellI, Fx, 0.5*(mass0 + mass1), Cud, dUTrans);
 
 
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Accumulate carrier phase source terms
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
     if (td.cloud().coupled())
     {
         // Transfer mass lost from particle to carrier mass source
         forAll(dMassMT, i)
         {
             td.cloud().rhoTrans(i)[cellI] +=
-                np0*(dMassPC[i] + dMassDV[i] + dMassSR[i]);
+                np0*(dMassPC[i] + dMassDV[i] + dMassSRc[i]);
         }
 
         // Update momentum transfer
@@ -200,9 +196,9 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     }
 
 
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Remove the particle when mass falls below minimum threshold
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
     if (mass1 < td.constProps().minParticleMass())
     {
         td.keepParticle = false;
@@ -215,14 +211,16 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         td.cloud().hTrans()[cellI] += np0*mass1*HEff(td, pc, T1);
         td.cloud().UTrans()[cellI] += np0*mass1*U1;
     }
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
     // Set new particle properties
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
     else
     {
         this->U_ = U1;
         this->T_ = T1;
-        this->cp_ = cp1;
+        this->cp_ = cpEff(td, pc, T1);
 
         // Update particle density or diameter
         if (td.constProps().constantVolume())
@@ -246,8 +244,8 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
     const scalar T,
     const scalar mass,
     const label idVolatile,
-    scalarField_ YMixture,
-    scalarList& dMassMT
+    scalarField& YMixture,
+    scalarList& dMassDV
 )
 {
     // Check that model is active, and that the parcel temperature is
@@ -268,16 +266,16 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
         dt,
         this->mass0_,
         mass,
-        td.cloud().composition().YMixture0()[idVolatile],
-        YMix[0],
         T,
+        td.cloud().composition().YMixture0()[idVolatile],
+        YMixture[0],
         canCombust_
     );
 
     // Update (total) mass fractions
-    YMix[0] = (YMix[0]*mass - dMassTot)/(mass - dMassTot);
-    YMix[1] = YMix[1]*mass/(mass - dMassTot);
-    YMix[2] = 1.0 - YMix[0] - YMix[1];
+    YMixture[0] = (YMixture[0]*mass - dMassTot)/(mass - dMassTot);
+    YMixture[1] = YMixture[1]*mass/(mass - dMassTot);
+    YMixture[2] = 1.0 - YMixture[0] - YMixture[1];
 
     // Add to cummulative mass transfer
     forAll (YGas_, i)
@@ -286,26 +284,26 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
 
         // Volatiles mass transfer
         scalar volatileMass = YGas_[i]*dMassTot;
-        dMassMT[id] += volatileMass;
+        dMassDV[id] += volatileMass;
     }
 
     td.cloud().addToMassDevolatilisation(dMassTot);
 
-    return = td.constProps().Ldevol()*dMassTot;
+    return td.constProps().Ldevol()*dMassTot;
 }
 
 
 template<class ParcelType>
 template<class TrackData>
-void Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
+Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
 (
     TrackData& td,
     const scalar dt,
     const label cellI,
     const scalar T,
     const scalarList& dMassMT,
-    scalarList& dMassSR,
-    scalar& dhRet
+    scalarList& dMassSRc,
+    scalar& dhTrans
 )
 {
     // Check that model is active
@@ -316,7 +314,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
 
     // Update mass transfer(s)
     // - Also updates Y()'s
-    td.cloud().surfaceReaction().calculate
+    scalar HReaction = td.cloud().surfaceReaction().calculate
     (
         dt,
         cellI,
@@ -331,11 +329,27 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
         YLiquid_,
         YSolid_,
         this->Y_,
-        dMassSR,
-        dhRet
+        dMassSRc,
+        HReaction
     );
 
+    // Heat of reaction divided between particle and carrier phase by the
+    // fraction fh and (1 - fh)
+    dhTrans -= (1.0 - td.constProps().fh())*HReaction;
+
+/*
+    // Correct dhTrans to account for enthalpy of consumed solids
+    dhTrans +=
+        sum(dMassSR)*td.cloud().composition().H(idSolid, YSolid_, pc, T0);
+
+    // Correct dhTrans to account for enthalpy of evolved volatiles
+    dhTrans +=
+        sum(dMassMT)*td.cloud().composition().H(idGas, YGas_, pc, T0);
+
     // TODO: td.cloud().addToMassSurfaceReaction(sum(dMassSR));
+*/
+
+    return td.constProps().fh()*HReaction;
 }
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
index ca6315bde89..35ff4f296b6 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
@@ -185,23 +185,23 @@ protected:
             TrackData& td,
             const scalar dt,
             const scalar T,
-	    const scalar mass,
-	    const label idVolatile,
-            scalarList& dMassMT
+            const scalar mass,
+            const label idVolatile,
+            scalarField& YMixture,
+            scalarList& dMassDV
         );
 
         //- Calculate surface reactions
         template<class TrackData>
-        void calcSurfaceReactions
+        scalar calcSurfaceReactions
         (
             TrackData& td,
             const scalar dt,
             const label cellI,
-            const scalar T0,
-            const scalar T1,
+            const scalar T,
             const scalarList& dMassMT,
-            scalarList& dMassSR,
-            scalar& dhRet
+            scalarList& dMassSRc,
+            scalar& dhTrans
         );
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
index 1d18da31697..b751f5c9ebe 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
@@ -69,7 +69,6 @@ void Foam::ReactingParcel<ParcelType>::calc
     const label cellI
 )
 {
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Define local properties at beginning of time step
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     const scalar mass0 = this->mass();
@@ -77,51 +76,54 @@ void Foam::ReactingParcel<ParcelType>::calc
     const scalar T0 = this->T_;
 
 
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~
-    // 1. Calculate phase change
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Mass transfer from particle to carrier phase
+    // Intialise transfer terms
+    // ~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Momentum
+    vector dUTrans = vector::zero;
+
+    // Enthalpy
+    scalar dhTrans = 0.0;
+
+    // Mass transfer due to phase change
     scalarList dMassPC(td.cloud().gases().size(), 0.0);
+
+
+    // Phase change
+    // ~~~~~~~~~~~~
+
+    // Return enthalpy source and calc mass transfer due to phase change
     scalar ShPC = calcPhaseChange(td, dt, cellI, T0, 0, 1.0, Y_, dMassPC);
 
     // Update particle component mass fractions
     updateMassFraction(mass0, dMassPC, Y_);
 
 
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // 2. Calculate heat transfer
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Heat transfer
+    // ~~~~~~~~~~~~~
+
+    // Calculate new particle temperature
     scalar htc = 0.0;
-    scalar ShHT = 0.0;
-    scalar T1 = calcHeatTransfer(td, dt, cellI, ShPC, htc, ShHT);
+    scalar T1 = calcHeatTransfer(td, dt, cellI, ShPC, htc, dhTrans);
 
 
-    // ~~~~~~~~~~~~~~~~~~~~~
-    // 3. Calculate velocity
-    // ~~~~~~~~~~~~~~~~~~~~~
+    // Motion
+    // ~~~~~~
+
     // Update mass
     scalar mass1 = mass0 - sum(dMassPC);
-    scalar Cud = 0.0;
-    vector dUTrans = vector::zero;
-    vector Fx = vector::zero;
-    vector U1 = calcVelocity(td, dt, Fx, 0.5*(mass0 + mass1), Cud, dUTrans);
-
-
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // 4. Collect contributions to determine new particle thermo properties
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
 
-    // Total enthalpy transfer from the particle to the carrier phase
-    scalar dhTrans = ShHT + ShPC;
+    // No additional forces
+    vector Fx = vector::zero;
 
-    // New specific heat capacity of mixture - using new temperature
-    scalar cp1 = td.cloud().composition().cp(0, Y_, pc_, T1);
+    // Calculate new particle velocity
+    scalar Cud = 0.0;
+    vector U1 =
+        calcVelocity(td, dt, cellI, Fx, 0.5*(mass0 + mass1), Cud, dUTrans);
 
 
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // 5. Accumulate carrier phase source terms
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Accumulate carrier phase source terms
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     if (td.cloud().coupled())
     {
         // Transfer mass lost from particle to carrier mass source
@@ -137,16 +139,15 @@ void Foam::ReactingParcel<ParcelType>::calc
         td.cloud().UCoeff()[cellI] += np0*mass0*Cud;
 
         // Update enthalpy transfer
-        td.cloud().hTrans()[cellI] += np0*dhTrans;
+        td.cloud().hTrans()[cellI] += np0*(dhTrans + ShPC);
 
         // Coefficient to be applied in carrier phase enthalpy coupling
         td.cloud().hCoeff()[cellI] += np0*htc*this->areaS();
     }
 
 
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // 6a. Remove the particle when mass falls below minimum threshold
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Remove the particle when mass falls below minimum threshold
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     if (mass1 < td.constProps().minParticleMass())
     {
         td.keepParticle = false;
@@ -159,18 +160,20 @@ void Foam::ReactingParcel<ParcelType>::calc
                 td.cloud().rhoTrans(i)[cellI] += np0*mass1*Y_[i];
             }
             td.cloud().UTrans()[cellI] += np0*mass1*U1;
-            scalar HEff = td.cloud().composition().H(0, YComponents, pc_, T1);
+            scalar HEff = td.cloud().composition().H(0, Y_, pc_, T1);
             td.cloud().hTrans()[cellI] += np0*mass1*HEff;
         }
     }
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // 6b. Set new particle properties
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
+    // Set new particle properties
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
     else
     {
         this->U_ = U1;
         this->T_ = T1;
-        this->cp_ = cp1;
+        this->cp_ = td.cloud().composition().cp(0, Y_, pc_, T1);
 
         // Update particle density or diameter
         if (td.constProps().constantVolume())
@@ -195,7 +198,7 @@ Foam::scalar Foam::ReactingParcel<ParcelType>::calcPhaseChange
     const scalar T,
     const label idPhase,
     const scalar YPhase,
-    scalarField& YComponents,
+    const scalarField& YComponents,
     scalarList& dMass
 )
 {
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
index 3e4a2480061..cde36a297dc 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
@@ -190,7 +190,7 @@ protected:
             const scalar T,
             const label idPhase,
             const scalar YPhase,
-	    const scalarField& YComponents,
+            const scalarField& YComponents,
             scalarList& dMass
         );
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
index b1ff448686c..d8b999d0f58 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
@@ -53,7 +53,6 @@ void Foam::ThermoParcel<ParcelType>::calc
     const label cellI
 )
 {
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Define local properties at beginning of time step
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     const vector U0 = this->U_;
@@ -61,25 +60,37 @@ void Foam::ThermoParcel<ParcelType>::calc
     const scalar np0 = this->nParticle_;
 
 
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // 1. Calculate heat transfer - update T
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Initialise transfer terms
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    vector dUTrans = vector::zero;
+    scalar dhTrans = 0.0;
+
+
+    // Heat transfer
+    // ~~~~~~~~~~~~~
+
+    // No additional enthalpy sources
+    vector Sh = 0.0;
+
+    // Calculate new particle velocity
     scalar htc = 0.0;
-    scalar ShHT = 0.0;
-    const scalar T1 = calcHeatTransfer(td, dt, cellI, 0.0, htc, ShHT);
+    scalar T1 = calcHeatTransfer(td, dt, cellI, Sh, htc, dhTrans);
+
+
+    // Motion
+    // ~~~~~~
 
+    // No additional forces
+    vector Fx = vector::zero;
 
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // 2. Calculate velocity - update U
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Calculate new particle velocity
     scalar Cud = 0.0;
-    vector dUTrans = vector::zero;
-    const vector U1 = calcVelocity(td, dt, vector::zero, Cud, mass0, dUTrans);
+    vector U1 = calcVelocity(td, dt, cellI, Fx, Cud, mass0, dUTrans);
 
 
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // 3. Accumulate carrier phase source terms
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    //  Accumulate carrier phase source terms
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     if (td.cloud().coupled())
     {
         // Update momentum transfer
@@ -89,16 +100,15 @@ void Foam::ThermoParcel<ParcelType>::calc
         td.cloud().UCoeff()[cellI] += np0*mass0*Cud;
 
         // Update enthalpy transfer
-        td.cloud().hTrans()[cellI] += np0*ShHT;
+        td.cloud().hTrans()[cellI] += np0*dhTrans;
 
         // Coefficient to be applied in carrier phase enthalpy coupling
         td.cloud().hCoeff()[cellI] += np0*htc*this->areaS();
     }
 
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // 5. Set new particle properties
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    this->U() = U1;
+    // Set new particle properties
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    this->U_ = U1;
     T_ = T1;
 }
 
@@ -118,7 +128,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
     if (!td.cloud().heatTransfer().active())
     {
         htc = 0.0;
-        dhTrans = 0.0;
+        ShHT = 0.0;
         return T_;
     }
 
@@ -134,7 +144,11 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
         this->muc_
     );
 
-    //- Assuming diameter = diameter at start of time step
+
+    // Set new particle temperature
+    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Assuming diameter = diameter at start of time step
     scalar Ap = this->areasS();
 
     // Determine ap and bp coefficients
@@ -161,7 +175,8 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
     IntegrationScheme<scalar>::integrationResult Tres =
         td.cloud().TIntegrator().integrate(T_, dt, ap, bp);
 
-    // Using average parcel temperature for enthalpy transfer calculation
+    // Enthalpy transfer
+    // - Using average particle temperature
     ShHT = dt*Ap*htc*(Tres.average() - Tc_);
 
     return Tres.value();
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
index fd33527bf53..584588d6958 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
@@ -199,9 +199,6 @@ void Foam::LiquidEvaporation<CloudType>::calculate
     scalarList& dMass
 ) const
 {
-    // initialise total mass transferred from the particle to carrier phase
-    scalar dMassTot = 0.0;
-
     // construct carrier phase species volume fractions for cell, cellI
     scalarField Xc = calcXc(cellI);
 
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/ConstantRateDevolatilisation/ConstantRateDevolatilisation.C b/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/ConstantRateDevolatilisation/ConstantRateDevolatilisation.C
index afe83f1baa8..77cea7e00b0 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/ConstantRateDevolatilisation/ConstantRateDevolatilisation.C
+++ b/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/ConstantRateDevolatilisation/ConstantRateDevolatilisation.C
@@ -66,9 +66,9 @@ Foam::scalar Foam::ConstantRateDevolatilisation<CloudType>::calculate
     const scalar dt,
     const scalar mass0,
     const scalar mass,
-    const scalar YVolatile0,
-    const scalarField& YVolatile,
     const scalar T,
+    const scalar YVolatile0,
+    const scalar YVolatile,
     bool& canCombust
 ) const
 {
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/ConstantRateDevolatilisation/ConstantRateDevolatilisation.H b/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/ConstantRateDevolatilisation/ConstantRateDevolatilisation.H
index 13abc8b2a96..a026bae9a94 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/ConstantRateDevolatilisation/ConstantRateDevolatilisation.H
+++ b/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/ConstantRateDevolatilisation/ConstantRateDevolatilisation.H
@@ -93,9 +93,9 @@ public:
             const scalar dt,
             const scalar mass0,
             const scalar mass,
-            const scalar YVolatile0,
-            const scalarField& YVolatile,
             const scalar T,
+            const scalar YVolatile0,
+            const scalar YVolatile,
             bool& canCombust
         ) const;
 };
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/DevolatilisationModel/DevolatilisationModel.H b/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/DevolatilisationModel/DevolatilisationModel.H
index 004482a040f..bb4e7ddc637 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/DevolatilisationModel/DevolatilisationModel.H
+++ b/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/DevolatilisationModel/DevolatilisationModel.H
@@ -135,9 +135,9 @@ public:
             const scalar dt,
             const scalar mass0,
             const scalar mass,
-            const scalarField& YMixture0,
-            const scalarField& YMixture,
             const scalar T,
+            const scalar YVolatile0,
+            const scalar YVolatile,
             bool& canCombust
         ) const = 0;
 };
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/NoDevolatilisation/NoDevolatilisation.C b/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/NoDevolatilisation/NoDevolatilisation.C
index 9d7ff147491..c0431a61f3b 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/NoDevolatilisation/NoDevolatilisation.C
+++ b/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/NoDevolatilisation/NoDevolatilisation.C
@@ -61,8 +61,8 @@ Foam::scalar Foam::NoDevolatilisation<CloudType>::calculate
     const scalar,
     const scalar,
     const scalar,
-    const scalarField&,
-    const scalarField&,
+    const scalar,
+    const scalar,
     const scalar,
     bool& canCombust
 ) const
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/NoDevolatilisation/NoDevolatilisation.H b/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/NoDevolatilisation/NoDevolatilisation.H
index 915547180a6..4dbea82607a 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/NoDevolatilisation/NoDevolatilisation.H
+++ b/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/NoDevolatilisation/NoDevolatilisation.H
@@ -76,8 +76,8 @@ public:
             const scalar,
             const scalar,
             const scalar,
-            const scalarField&,
-            const scalarField&,
+            const scalar,
+            const scalar,
             const scalar,
             bool&
         ) const;
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/SingleKineticRateDevolatilisation/SingleKineticRateDevolatilisation.C b/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/SingleKineticRateDevolatilisation/SingleKineticRateDevolatilisation.C
index 488e22978ef..2d09c134a97 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/SingleKineticRateDevolatilisation/SingleKineticRateDevolatilisation.C
+++ b/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/SingleKineticRateDevolatilisation/SingleKineticRateDevolatilisation.C
@@ -69,9 +69,9 @@ Foam::scalar Foam::SingleKineticRateDevolatilisation<CloudType>::calculate
     const scalar dt,
     const scalar mass0,
     const scalar mass,
-    const scalar YVolatile0,
-    const scalarField& YVolatile,
     const scalar T,
+    const scalar YVolatile0,
+    const scalar YVolatile,
     bool& canCombust
 ) const
 {
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/SingleKineticRateDevolatilisation/SingleKineticRateDevolatilisation.H b/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/SingleKineticRateDevolatilisation/SingleKineticRateDevolatilisation.H
index 62873a77d71..91bd5f06122 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/SingleKineticRateDevolatilisation/SingleKineticRateDevolatilisation.H
+++ b/src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/SingleKineticRateDevolatilisation/SingleKineticRateDevolatilisation.H
@@ -95,9 +95,9 @@ public:
             const scalar dt,
             const scalar mass0,
             const scalar mass,
-            const scalar YVolatile0,
-            const scalarField& YVolatile,
             const scalar T,
+            const scalar YVolatile0,
+            const scalar YVolatile,
             bool& canCombust
         ) const;
 };
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.C b/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.C
index 8556282b9df..4b92bb38dfb 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.C
+++ b/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.C
@@ -56,27 +56,26 @@ bool Foam::NoSurfaceReaction<CloudType>::active() const
 
 
 template<class CloudType>
-void Foam::NoSurfaceReaction<CloudType>::calculate
+Foam::scalar Foam::NoSurfaceReaction<CloudType>::calculate
 (
-    const scalar dt,
-    const label celli,
-    const scalar dp,
-    const scalar T0,
-    const scalar T1,
-    const scalar Tc,
-    const scalar pc,
-    const scalar rhoc,
-    const scalar massp,
-    const scalarList& dMassMT,
-    scalarField& YGas,
-    scalarField& YLiquid,
-    scalarField& YSolid,
-    scalarField& YMixture,
-    scalarList& dMassSR,
-    scalar& dhRet
+    const scalar,
+    const label,
+    const scalar,
+    const scalar,
+    const scalar,
+    const scalar,
+    const scalar,
+    const scalar,
+    const scalarList&,
+    scalarField&,
+    scalarField&,
+    scalarField&,
+    scalarField&,
+    scalarList&
 ) const
 {
     // do nothing
+    return 0.0;
 }
 
 
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.H b/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.H
index edb919a6214..75b729baa7f 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.H
+++ b/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/NoSurfaceReaction/NoSurfaceReaction.H
@@ -74,24 +74,22 @@ public:
         virtual bool active() const;
 
         //- Update surface reactions
-        virtual void calculate
+        virtual scalar calculate
         (
             const scalar dt,
-            const label celli,
-            const scalar dp,
-            const scalar T0,
-            const scalar T1,
+            const label cellI,
+            const scalar d,
+            const scalar T,
             const scalar Tc,
             const scalar pc,
             const scalar rhoc,
-            const scalar massp,
+            const scalar mass,
             const scalarList& dMassMT,
             scalarField& YGas,
             scalarField& YLiquid,
             scalarField& YSolid,
             scalarField& YMixture,
-            scalarList& dMassSR,
-            scalar& dhRet
+            scalarList& dMassSRc
         ) const;
 };
 
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/SurfaceReactionModel/SurfaceReactionModel.H b/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/SurfaceReactionModel/SurfaceReactionModel.H
index 6b4f5be31cb..5bc62e89eb0 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/SurfaceReactionModel/SurfaceReactionModel.H
+++ b/src/lagrangian/intermediate/submodels/ReactingMultiphase/SurfaceReactionModel/SurfaceReactionModel/SurfaceReactionModel.H
@@ -132,24 +132,23 @@ public:
         virtual bool active() const = 0;
 
         //- Update surface reactions
-        virtual void calculate
+        //  Returns the heat of reaction
+        virtual scalar calculate
         (
             const scalar dt,
-            const label celli,
-            const scalar dp,
-            const scalar T0,
-            const scalar T1,
+            const label cellI,
+            const scalar d,
+            const scalar T,
             const scalar Tc,
             const scalar pc,
             const scalar rhoc,
-            const scalar massp,
+            const scalar mass,
             const scalarList& dMassMT,
             scalarField& YGas,
             scalarField& YLiquid,
             scalarField& YSolid,
             scalarField& YMixture,
-            scalarList& dMassSR,
-            scalar& dhRet
+            scalarList& dMassSRc
         ) const = 0;
 };
 
-- 
GitLab