updates to the reacting classes

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

@@ -87,146 +87,102 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Define local properties at beginning of timestep
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    const vector U0 = this->U_;
     const scalar mass0 = this->mass();
-    const scalar cp0 = this->cp_;
     const scalar np0 = this->nParticle_;
-    scalar T0 = this->T_;
+    const scalar T0 = this->T_;
     const scalar pc = this->pc_;
-    scalarField& Y = this->Y_;
+    scalarField& YMix = this->Y_;
-    scalar mass1 = mass0;
+    const label idL = td.cloud().composition().idLiquid();
 
-    label idGas = td.cloud().composition().idGas();
-    label idLiquid = td.cloud().composition().idLiquid();
-    label idSolid = td.cloud().composition().idSolid();
-
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Initialise transfer terms
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Momentum transfer from the particle to the carrier phase
-    vector dUTrans = vector::zero;
-
-    // Enthalpy transfer from the particle to the carrier phase
-    scalar dhTrans = 0.0;
 
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Calculate phase change in liquid phase
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Mass transfer from particle to carrier phase
-    // - components exist in particle already
+    scalarList dMassPC(td.cloud().gases().size(), 0.0);
-    scalarList dMassMT(td.cloud().gases().size(), 0.0);
+    scalar shPC = 
+        calcPhaseChange(td, dt, cellI, T0, idL, YMix[idL], YLiquid_, dMassPC);
 
-    // Mass transfer due to surface reactions from particle to carrier phase
+    // Update particle component mass fractions
-    // - components do not necessarily exist in particle already
+    updateMassFraction(mass0, dMassPC, YLiquid_);
-    scalarList dMassSR(td.cloud().gases().size(), 0.0);
-
-
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Calculate velocity - update U
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    scalar Cud = 0.0;
-    const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
-
-
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Calculate heat transfer - update T
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    scalar htc = 0.0;
-    scalar T1 = calcHeatTransfer(td, dt, cellI, htc, dhTrans);
-
-
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Calculate phase change - update mass, Y, cp, T, dhTrans
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    const scalar dMassPC = calcPhaseChange(td, dt, cellI, T1, Y[idLiquid], dMassMT);
-
-    // Update particle mass
-    mass1 -= dMassPC;
-
-    // Update particle liquid component mass fractions
-    this->updateMassFraction(mass0, dMassMT, YLiquid_);
-
-    // New specific heat capacity of mixture
-    scalar cp1 = cpEff(td, pc, T1);
-
-    // Correct temperature due to evaporated components
-    // TODO: use hl function in liquidMixture???
-//    scalar dhPC = -dMassPCTot*td.cloud().composition().L(0, Y_, pc, T0);
-    scalar Lvap = td.cloud().composition().L(idLiquid, this->YLiquid_, pc, T0);
-    T1 -= Lvap*dMassPC/(0.5*(mass0 + mass1)*cp1);
-
-    // Correct dhTrans to account for the change in enthalpy due to the
-    // liquid phase change
-    dhTrans +=
-        dMassPC
-       *td.cloud().composition().H(idLiquid, YLiquid_, pc, 0.5*(T0 + T1));
-
-//????????????????????????????????????????????????????????????????????????????????
-    // Store temperature for the start of the next process
-    T0 = T1;
-//????????????????????????????????????????????????????????????????????????????????
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Calculate Devolatilisation
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
-    const scalar dMassD = calcDevolatilisation(td, dt, T0, dMassMT);
+    // Mass transfer from particle to carrier phase
-
+    // - components exist in particle already
-    // Update particle mass
+    scalarList dMassDV(td.cloud().gases().size(), 0.0);
-    mass1 -= dMassD;
+    scalar shDV = calcDevolatilisation(td, dt, T0, mass0, idGas, YMix, dMassDV);
-
-    // New specific heat capacity of mixture
-    cp1 = cpEff(td, pc, T1);
-
-    // Update gas and solid component mass fractions
-//    updateMassFraction(mass0, mass1, dMassMT, YGas_);
-//    updateMassFraction(mass0, mass1, dMassMT, YSolid_);
-
-    // Correct particle temperature to account for latent heat of
-    // devolatilisation
-    T1 -=
-        td.constProps().Ldevol()
-       *sum(dMassMT)
-       /(0.5*(mass0 + mass1)*cp1);
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // 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);
-    calcSurfaceReactions(td, dt, cellI, T0, T1, dMassMT, dMassSR, dhRet);
-
-    // New total mass
-    mass1 -= sum(dMassSR);
-
 
     // Enthalpy retention divided between particle and carrier phase by the
     // fraction fh and (1 - fh)
-    T1 += td.constProps().fh()*dhRet/(0.5*(mass0 + mass1)*cp0);
+    scalar ShSR = td.constProps().fh()*dhRet;
     dhTrans -= (1.0 - td.constProps().fh())*dhRet;
 
-    // Correct dhTrans to account for enthalpy of evolved volatiles
-    dhTrans +=
-        sum(dMassMT)
-       *td.cloud().composition().H(idGas, YGas_, pc, 0.5*(T0 + T1));
-
     // Correct dhTrans to account for enthalpy of consumed solids
     dhTrans +=
-        sum(dMassSR)
+        sum(dMassSR)*td.cloud().composition().H(idSolid, YSolid_, pc, T0);
-       *td.cloud().composition().H(idSolid, YSolid_, pc, 0.5*(T0 + T1));
+
+    // Correct dhTrans to account for enthalpy of evolved volatiles
+    dhTrans +=
+        sum(dMassMT)*td.cloud().composition().H(idGas, YGas_, pc, T0);
 
 
+    // ~~~~~~~~~~~~~~~~~~~~~~~
+    // Calculate heat transfer
+    // ~~~~~~~~~~~~~~~~~~~~~~~
+    scalar htc = 0.0;
+
+    // Total enthalpy source
+    scalar Sh = ShPC + ShDV + ShSR;
+
+    scalar T1 = calcHeatTransfer(td, dt, cellI, Sh, htc, shHT);
+
+
+    // ~~~~~~~~~~~~~~~~~~
+    // 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);
+
+
+
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Collect contributions to determine new particle thermo properties
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Update specific heat capacity
+    cp1 = cpEff(td, pc, T1);
+
+
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Accumulate carrier phase source terms
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     if (td.cloud().coupled())
     {
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-        // Accumulate carrier phase source terms
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
         // Transfer mass lost from particle to carrier mass source
         forAll(dMassMT, i)
         {
-            td.cloud().rhoTrans(i)[cellI] += np0*(dMassMT[i] + dMassSR[i]);
+            td.cloud().rhoTrans(i)[cellI] +=
+                np0*(dMassPC[i] + dMassDV[i] + dMassSR[i]);
         }
 
         // Update momentum transfer
@@ -266,7 +222,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     {
         this->U_ = U1;
         this->T_ = T1;
-        this->cp_ = cpEff(td, pc, T1);
+        this->cp_ = cp1;
 
         // Update particle density or diameter
         if (td.constProps().constantVolume())
@@ -288,6 +244,9 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
     TrackData& td,
     const scalar dt,
     const scalar T,
+    const scalar mass,
+    const label idVolatile,
+    scalarField_ YMixture,
     scalarList& dMassMT
 )
 {
@@ -296,37 +255,34 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
     if
     (
         !td.cloud().devolatilisation().active()
-     || this->T_<td.constProps().Tvap()
+     || T < td.constProps().Tvap()
-     || this->T_<td.constProps().Tbp()
+     || T < td.constProps().Tbp()
     )
     {
         return 0.0;
     }
 
     // Determine mass to add to carrier phase
-    const scalar mass = this->mass();
-    scalarField& Y = this->Y_;
     const scalar dMassTot = td.cloud().devolatilisation().calculate
     (
         dt,
         this->mass0_,
         mass,
-        td.cloud().composition().YMixture0(),
+        td.cloud().composition().YMixture0()[idVolatile],
-        Y,
+        YMix[0],
         T,
         canCombust_
     );
 
     // Update (total) mass fractions
-    Y[0] = (Y[0]*mass - dMassTot)/(mass - dMassTot);
+    YMix[0] = (YMix[0]*mass - dMassTot)/(mass - dMassTot);
-    Y[1] = Y[1]*mass/(mass - dMassTot);
+    YMix[1] = YMix[1]*mass/(mass - dMassTot);
-    Y[2] = 1.0 - Y[0] - Y[1];
+    YMix[2] = 1.0 - YMix[0] - YMix[1];
 
     // Add to cummulative mass transfer
-    label idGas = td.cloud().composition().idGas();
     forAll (YGas_, i)
     {
-        label id = td.cloud().composition().globalIds(idGas)[i];
+        label id = td.cloud().composition().globalIds(idVolatile)[i];
 
         // Volatiles mass transfer
         scalar volatileMass = YGas_[i]*dMassTot;
@@ -335,7 +291,7 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
 
     td.cloud().addToMassDevolatilisation(dMassTot);
 
-    return dMassTot;
+    return = td.constProps().Ldevol()*dMassTot;
 }
 
 
@@ -346,8 +302,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
     TrackData& td,
     const scalar dt,
     const label cellI,
-    const scalar T0,
+    const scalar T,
-    const scalar T1,
     const scalarList& dMassMT,
     scalarList& dMassSR,
     scalar& dhRet
@@ -366,8 +321,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
         dt,
         cellI,
         this->d_,
-        T0,
+        T,
-        T1,
         this->Tc_,
         this->pc_,
         this->rhoc_,

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

@@ -185,6 +185,8 @@ protected:
             TrackData& td,
             const scalar dt,
             const scalar T,
+	    const scalar mass,
+	    const label idVolatile,
             scalarList& dMassMT
         );
 

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

@@ -76,47 +76,41 @@ void Foam::ReactingParcel<ParcelType>::calc
     const scalar np0 = this->nParticle_;
     const scalar T0 = this->T_;
 
-    // ~~~~~~~~~~~~~~~~~~~~~
-    // 1. Calculate velocity
-    // ~~~~~~~~~~~~~~~~~~~~~
-
-    scalar Cud = 0.0;
-    vector dUTrans = vector::zero;
-    const vector U1 = calcVelocity(td, dt, vector::zero, mass0, Cud, dUTrans);
-
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~
-    // 2. Calculate phase change
+    // 1. Calculate phase change
     // ~~~~~~~~~~~~~~~~~~~~~~~~~
     // Mass transfer from particle to carrier phase
     scalarList dMassPC(td.cloud().gases().size(), 0.0);
-    const scalar dMassPCTot = calcPhaseChange(td, dt, cellI, T0, 1.0, dMassPC);
+    scalar ShPC = calcPhaseChange(td, dt, cellI, T0, 0, 1.0, Y_, dMassPC);
-
-    // Enthalpy change due to change in particle composition (sink)
-    scalar ShPC = -dMassPCTot*td.cloud().composition().L(0, Y_, pc, T0);
-
-    // Enthalpy change due to species released into the carrier (source)
-    scalar HEff = td.cloud().composition().H(0, Y_, pc, T0);
-    ShPC += dMassPCTot*HEff;
 
     // Update particle component mass fractions
     updateMassFraction(mass0, dMassPC, Y_);
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // 3. Calculate heat transfer
+    // 2. Calculate heat transfer
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
     scalar htc = 0.0;
     scalar ShHT = 0.0;
     scalar T1 = calcHeatTransfer(td, dt, cellI, ShPC, htc, ShHT);
 
 
+    // ~~~~~~~~~~~~~~~~~~~~~
+    // 3. Calculate velocity
+    // ~~~~~~~~~~~~~~~~~~~~~
+    // 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
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    // New mass
-    scalar mass1 = mass0 - dMassPCTot;
 
     // Total enthalpy transfer from the particle to the carrier phase
     scalar dhTrans = ShHT + ShPC;
@@ -165,6 +159,7 @@ 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);
             td.cloud().hTrans()[cellI] += np0*mass1*HEff;
         }
     }
@@ -198,36 +193,47 @@ Foam::scalar Foam::ReactingParcel<ParcelType>::calcPhaseChange
     const scalar dt,
     const label cellI,
     const scalar T,
-    const scalar YPhase, // TODO: NEEDED?????????????????????????????????????????
+    const label idPhase,
-    scalarList& dMassMT
+    const scalar YPhase,
+    scalarField& YComponents,
+    scalarList& dMass
 )
 {
     if
     (
         !td.cloud().phaseChange().active()
      || T < td.constProps().Tvap()
-     || YPhase > SMALL
+     || YPhase < SMALL
     )
     {
         return 0.0;
     }
 
-    scalar dMass = td.cloud().phaseChange().calculate
+    td.cloud().phaseChange().calculate
     (
         dt,
         cellI,
-        T,
+        min(T, td.constProps().Tbp()), // Limiting to boiling temperature
         pc_,
         this->d_,
         this->Tc_,
         this->muc_/this->rhoc_,
         this->U_ - this->Uc_,
-        dMassMT
+        dMass
     );
 
-    td.cloud().addToMassPhaseChange(dMass);
+    scalar dMassTot = sum(dMass);
 
-    return dMass;
+    // Add to cumulative phase change mass
+    td.cloud().addToMassPhaseChange(dMassTot);
+
+    // Effective latent heat of vaporisation
+    scalar LEff = td.cloud().composition().L(idPhase, YComponents, pc_, T);
+
+    // Effective heat of vaporised components
+    scalar HEff = td.cloud().composition().H(idPhase, YComponents, pc_, T);
+
+    return dMassTot*(HEff - LEff);
 }
 
 

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

@@ -188,8 +188,10 @@ protected:
             const scalar dt,
             const label cellI,
             const scalar T,
+            const label idPhase,
             const scalar YPhase,
-            scalarList& dMassMT
+	    const scalarField& YComponents,
+            scalarList& dMass
         );
 
         //- Update mass fraction

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

@@ -61,33 +61,24 @@ void Foam::ThermoParcel<ParcelType>::calc
     const scalar np0 = this->nParticle_;
 
 
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // 1. Initialise transfer terms
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Momentum transfer from the particle to the carrier phase
-    vector dUTrans = vector::zero;
-
-    // Enthalpy transfer from the particle to the carrier phase
-    scalar dhTrans = 0.0;
-
-
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // 2. Calculate heat transfer - update T
+    // 1. Calculate heat transfer - update T
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     scalar htc = 0.0;
-    const scalar T1 = calcHeatTransfer(td, dt, cellI, 0.0, htc, dhTrans);
+    scalar ShHT = 0.0;
+    const scalar T1 = calcHeatTransfer(td, dt, cellI, 0.0, htc, ShHT);
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // 3. Calculate velocity - update U
+    // 2. Calculate velocity - update U
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     scalar Cud = 0.0;
+    vector dUTrans = vector::zero;
     const vector U1 = calcVelocity(td, dt, vector::zero, Cud, mass0, dUTrans);
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // 4. Accumulate carrier phase source terms
+    // 3. Accumulate carrier phase source terms
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     if (td.cloud().coupled())
     {
@@ -98,7 +89,7 @@ void Foam::ThermoParcel<ParcelType>::calc
         td.cloud().UCoeff()[cellI] += np0*mass0*Cud;
 
         // Update enthalpy transfer
-        td.cloud().hTrans()[cellI] += np0*dhTrans;
+        td.cloud().hTrans()[cellI] += np0*ShHT;
 
         // Coefficient to be applied in carrier phase enthalpy coupling
         td.cloud().hCoeff()[cellI] += np0*htc*this->areaS();
@@ -121,7 +112,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
     const label cellI,
     const scalar Sh,
     scalar& htc,
-    scalar& dhTrans
+    scalar& ShHT
 )
 {
     if (!td.cloud().heatTransfer().active())
@@ -171,7 +162,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
         td.cloud().TIntegrator().integrate(T_, dt, ap, bp);
 
     // Using average parcel temperature for enthalpy transfer calculation
-    dhTrans = dt*Ap*htc*(Tres.average() - Tc_);
+    ShHT = dt*Ap*htc*(Tres.average() - Tc_);
 
     return Tres.value();
 }

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

@@ -186,7 +186,7 @@ bool Foam::LiquidEvaporation<CloudType>::active() const
 
 
 template<class CloudType>
-Foam::scalar Foam::LiquidEvaporation<CloudType>::calculate
+void Foam::LiquidEvaporation<CloudType>::calculate
 (
     const scalar dt,
     const label cellI,
@@ -196,7 +196,7 @@ Foam::scalar Foam::LiquidEvaporation<CloudType>::calculate
     const scalar Tc,
     const scalar nuc,
     const vector& Ur,
-    scalarList& dMassMT
+    scalarList& dMass
 ) const
 {
     // initialise total mass transferred from the particle to carrier phase
@@ -242,12 +242,8 @@ Foam::scalar Foam::LiquidEvaporation<CloudType>::calculate
         scalar Ni = max(kc*(Cs - Cinf), 0.0);
 
         // mass transfer [kg]
-        scalar dm = Ni*A*liquids_->properties()[lid].W()*dt;
+        dMass[gid] += Ni*A*liquids_->properties()[lid].W()*dt;
-        dMassMT[gid] += dm;
-        dMassTot += dm;
     }
-
-    return dMassTot;
 }
 
 

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

@@ -104,7 +104,7 @@ public:
         bool active() const;
 
         //- Update model
-        scalar calculate
+        void calculate
         (
             const scalar dt,
             const label cellI,
@@ -114,7 +114,7 @@ public:
             const scalar Tc,
             const scalar nuc,
             const vector& Ur,
-            scalarList& dMassMT
+            scalarList& dMass
         ) const;
 };
 

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

@@ -56,7 +56,7 @@ bool Foam::NoPhaseChange<CloudType>::active() const
 
 
 template<class CloudType>
-Foam::scalar Foam::NoPhaseChange<CloudType>::calculate
+void Foam::NoPhaseChange<CloudType>::calculate
 (
     const scalar,
     const label,
@@ -70,7 +70,6 @@ Foam::scalar Foam::NoPhaseChange<CloudType>::calculate
 ) const
 {
     // Nothing to do...
-    return 0.0;
 }
 
 

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

@@ -70,7 +70,7 @@ public:
         bool active() const;
 
         //- Update model
-        scalar calculate
+        void calculate
         (
             const scalar dt,
             const label cellI,
@@ -80,7 +80,7 @@ public:
             const scalar Tc,
             const scalar nuc,
             const vector& Ur,
-            scalarList& dMassMT
+            scalarList& dMass
         ) const;
 };
 

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

@@ -136,7 +136,7 @@ public:
         virtual bool active() const = 0;
 
         //- Update model
-        virtual scalar calculate
+        virtual void calculate
         (
             const scalar dt,
             const label cellI,
@@ -146,7 +146,7 @@ public:
             const scalar Tc,
             const scalar nuc,
             const vector& Ur,
-            scalarList& dMassMT
+            scalarList& dMass
         ) const = 0;
 };
 

src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/ConstantRateDevolatilisation/ConstantRateDevolatilisation.C

@@ -66,14 +66,14 @@ Foam::scalar Foam::ConstantRateDevolatilisation<CloudType>::calculate
     const scalar dt,
     const scalar mass0,
     const scalar mass,
-    const scalarField& YMixture0,
+    const scalar YVolatile0,
-    const scalarField& YMixture,
+    const scalarField& YVolatile,
     const scalar T,
     bool& canCombust
 ) const
 {
-    const scalar massVolatile0 = YMixture0[0]*mass0;
+    const scalar massVolatile0 = YVolatile0*mass0;
-    const scalar massVolatile  = YMixture[0]*mass;
+    const scalar massVolatile  = YVolatile*mass;
 
     if (massVolatile <= volatileResidualCoeff_*massVolatile0)
     {

src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/ConstantRateDevolatilisation/ConstantRateDevolatilisation.H

@@ -93,8 +93,8 @@ public:
             const scalar dt,
             const scalar mass0,
             const scalar mass,
-            const scalarField& YMixture0,
+            const scalar YVolatile0,
-            const scalarField& YMixture,
+            const scalarField& YVolatile,
             const scalar T,
             bool& canCombust
         ) const;

src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/SingleKineticRateDevolatilisation/SingleKineticRateDevolatilisation.C

@@ -29,7 +29,8 @@ License
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template <class CloudType>
-Foam::SingleKineticRateDevolatilisation<CloudType>::SingleKineticRateDevolatilisation
+Foam::SingleKineticRateDevolatilisation<CloudType>::
+SingleKineticRateDevolatilisation
 (
     const dictionary& dict,
     CloudType& owner
@@ -68,14 +69,14 @@ Foam::scalar Foam::SingleKineticRateDevolatilisation<CloudType>::calculate
     const scalar dt,
     const scalar mass0,
     const scalar mass,
-    const scalarField& YMixture0,
+    const scalar YVolatile0,
-    const scalarField& YMixture,
+    const scalarField& YVolatile,
     const scalar T,
     bool& canCombust
 ) const
 {
-    const scalar massVolatile0 = YMixture0[0]*mass;
+    const scalar massVolatile0 = YVolatile0*mass;
-    const scalar massVolatile = YMixture[0]*mass;
+    const scalar massVolatile = YVolatile*mass;
 
     if (massVolatile <= volatileResidualCoeff_*massVolatile0)
     {

src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/SingleKineticRateDevolatilisation/SingleKineticRateDevolatilisation.H

@@ -95,8 +95,8 @@ public:
             const scalar dt,
             const scalar mass0,
             const scalar mass,
-            const scalarField& YMixture0,
+            const scalar YVolatile0,
-            const scalarField& YMixture,
+            const scalarField& YVolatile,
             const scalar T,
             bool& canCombust
         ) const;