ENH: Updated cloud source term computations

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

@@ -104,13 +104,10 @@ void Foam::KinematicParcel<ParcelType>::calc
     // Define local properties at beginning of time step
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     const scalar np0 = nParticle_;
-    const scalar d0 = d_;
-    const vector U0 = U_;
-    const scalar rho0 = rho_;
     const scalar mass0 = mass();
 
     // Reynolds number
-    const scalar Re = this->Re(U0, d0, rhoc_, muc_);
+    const scalar Re = this->Re(U_, d_, rhoc_, muc_);
 
 
     // Sources
@@ -119,6 +116,9 @@ void Foam::KinematicParcel<ParcelType>::calc
     // Explicit momentum source for particle
     vector Su = vector::zero;
 
+    // Linearised momentum source coefficient
+    scalar Spu = 0.0;
+
     // Momentum transfer from the particle to the carrier phase
     vector dUTrans = vector::zero;
 
@@ -127,23 +127,7 @@ void Foam::KinematicParcel<ParcelType>::calc
     // ~~~~~~
 
     // Calculate new particle velocity
-    scalar Spu = 0.0;
+    this->U_ = calcVelocity(td, dt, cellI, Re, muc_, mass0, Su, dUTrans, Spu);
-    vector U1 =
-        calcVelocity
-        (
-            td,
-            dt,
-            cellI,
-            Re,
-            muc_,
-            d0,
-            U0,
-            rho0,
-            mass0,
-            Su,
-            dUTrans,
-            Spu
-        );
 
 
     // Accumulate carrier phase source terms
@@ -156,11 +140,6 @@ void Foam::KinematicParcel<ParcelType>::calc
         // Update momentum transfer coefficient
         td.cloud().UCoeff()[cellI] += np0*Spu;
     }
-
-
-    // Set new particle properties
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    U_ = U1;
 }
 
 
@@ -173,9 +152,6 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
     const label cellI,
     const scalar Re,
     const scalar mu,
-    const scalar d,
-    const vector& U,
-    const scalar rho,
     const scalar mass,
     const vector& Su,
     vector& dUTrans,
@@ -205,7 +181,7 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
     Spu = dt*Feff.Sp();
 
     IntegrationScheme<vector>::integrationResult Ures =
-        td.cloud().UIntegrator().integrate(U, dt, abp, bp);
+        td.cloud().UIntegrator().integrate(U_, dt, abp, bp);
 
     vector Unew = Ures.value();
 

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

@@ -297,9 +297,6 @@ protected:
             const label cellI,         // owner cell
             const scalar Re,           // Reynolds number
             const scalar mu,           // local carrier viscosity
-            const scalar d,            // diameter
-            const vector& U,           // velocity
-            const scalar rho,          // density
             const scalar mass,         // mass
             const vector& Su,          // explicit particle momentum source
             vector& dUTrans,           // momentum transfer to carrier

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

@@ -172,12 +172,11 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
 
     // Define local properties at beginning of timestep
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
     const scalar np0 = this->nParticle_;
     const scalar d0 = this->d_;
     const vector& U0 = this->U_;
-    const scalar rho0 = this->rho_;
     const scalar T0 = this->T_;
-    const scalar Cp0 = this->Cp_;
     const scalar mass0 = this->mass();
 
     const scalar pc = this->pc_;
@@ -189,11 +188,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
 
 
     // Calc surface values
-    // ~~~~~~~~~~~~~~~~~~~
     scalar Ts, rhos, mus, Prs, kappas;
     this->calcSurfaceValues(td, cellI, T0, Ts, rhos, mus, Prs, kappas);
-
-    // Reynolds number
     scalar Res = this->Re(U0, d0, rhos, mus);
 
 
@@ -203,16 +199,25 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     // Explicit momentum source for particle
     vector Su = vector::zero;
 
+    // Linearised momentum source coefficient
+    scalar Spu = 0.0;
+
     // Momentum transfer from the particle to the carrier phase
     vector dUTrans = vector::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;
 
 
+    // 1. Compute models that contribute to mass transfer - U, T held constant
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
     // Phase change in liquid phase
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -310,28 +315,78 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     );
 
 
-    // Correct surface values due to emitted species
+    // 2. Update the parcel properties due to change in mass
-    this->correctSurfaceValues(td, cellI, Ts, Cs, rhos, mus, Prs, kappas);
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    Res = this->Re(U0, d0, rhos, mus);
-
-
-    // Update component mass fractions
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     scalarField dMassGas(dMassDV + dMassSRGas);
     scalarField dMassLiquid(dMassPC + dMassSRLiquid);
     scalarField dMassSolid(dMassSRSolid);
-
     scalar mass1 =
         updateMassFractions(mass0, dMassGas, dMassLiquid, dMassSolid);
 
+    this->Cp_ = CpEff(td, pc, T0, idG, idL, idS);
+
+    // Update particle density or diameter
+    if (td.cloud().constProps().constantVolume())
+    {
+        this->rho_ = mass1/this->volume();
+    }
+    else
+    {
+        this->d_ = cbrt(mass1/this->rho_*6.0/pi);
+    }
+
+    // Remove the particle when mass falls below minimum threshold
+    if (np0*mass1 < td.cloud().constProps().minParticleMass())
+    {
+        td.keepParticle = false;
+
+        if (td.cloud().solution().coupled())
+        {
+            scalar dm = np0*mass1;
+
+            // Absorb parcel into carrier phase
+            forAll(YGas_, i)
+            {
+                label gid = composition.localToGlobalCarrierId(GAS, i);
+                td.cloud().rhoTrans(gid)[cellI] += dm*YMix[GAS]*YGas_[i];
+            }
+            forAll(YLiquid_, i)
+            {
+                label gid = composition.localToGlobalCarrierId(LIQ, i);
+                td.cloud().rhoTrans(gid)[cellI] += dm*YMix[LIQ]*YLiquid_[i];
+            }
+/*
+            // No mapping between solid components and carrier phase
+            forAll(YSolid_, i)
+            {
+                label gid = composition.localToGlobalCarrierId(SLD, i);
+                td.cloud().rhoTrans(gid)[cellI] += dm*YMix[SLD]*YSolid_[i];
+            }
+*/
+            td.cloud().UTrans()[cellI] += dm*U0;
+
+            td.cloud().hsTrans()[cellI] += dm*HsEff(td, pc, T0, idG, idL, idS);
+
+            td.cloud().addToMassPhaseChange(dm);
+        }
+
+        return;
+    }
+
+    // Correct surface values due to emitted species
+    this->correctSurfaceValues(td, cellI, Ts, Cs, rhos, mus, Prs, kappas);
+    Res = this->Re(U0, this->d_, rhos, mus);
+
+
+    // 3. Compute heat- and momentum transfers
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     // Heat transfer
     // ~~~~~~~~~~~~~
 
     // Calculate new particle temperature
-    scalar Sph = 0.0;
+    this->T_ =
-    scalar T1 =
         this->calcHeatTransfer
         (
             td,
@@ -340,10 +395,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
             Res,
             Prs,
             kappas,
-            d0,
-            rho0,
-            T0,
-            Cp0,
             NCpW,
             Sh,
             dhsTrans,
@@ -351,35 +402,23 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         );
 
 
+    this->Cp_ = CpEff(td, pc, this->T_, idG, idL, idS);
+
+
     // Motion
     // ~~~~~~
 
     // Calculate new particle velocity
-    scalar Spu = 0;
+    this->U_ =
-    vector U1 =
+        this->calcVelocity(td, dt, cellI, Res, mus, mass1, Su, dUTrans, Spu);
-        this->calcVelocity
-        (
-            td,
-            dt,
-            cellI,
-            Res,
-            mus,
-            d0,
-            U0,
-            rho0,
-            mass0,
-            Su,
-            dUTrans,
-            Spu
-        );
 
 
-    // Accumulate carrier phase source terms
+    // 4. Accumulate carrier phase source terms
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     if (td.cloud().solution().coupled())
     {
-        // Transfer mass lost from particle to carrier mass source
+        // Transfer mass lost to carrier mass, momentum and enthalpy sources
         forAll(YGas_, i)
         {
             scalar dm = np0*dMassGas[i];
@@ -421,76 +460,12 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
 
         // Update momentum transfer
         td.cloud().UTrans()[cellI] += np0*dUTrans;
-
-        // Update momentum transfer coefficient
         td.cloud().UCoeff()[cellI] += np0*Spu;
 
         // Update sensible enthalpy transfer
         td.cloud().hsTrans()[cellI] += np0*dhsTrans;
-
-        // Update sensible enthalpy coefficient
         td.cloud().hsCoeff()[cellI] += np0*Sph;
     }
-
-
-    // Remove the particle when mass falls below minimum threshold
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    if (np0*mass1 < td.cloud().constProps().minParticleMass())
-    {
-        td.keepParticle = false;
-
-        if (td.cloud().solution().coupled())
-        {
-            scalar dm = np0*mass1;
-
-            // Absorb parcel into carrier phase
-            forAll(YGas_, i)
-            {
-                label gid = composition.localToGlobalCarrierId(GAS, i);
-                td.cloud().rhoTrans(gid)[cellI] += dm*YMix[GAS]*YGas_[i];
-            }
-            forAll(YLiquid_, i)
-            {
-                label gid = composition.localToGlobalCarrierId(LIQ, i);
-                td.cloud().rhoTrans(gid)[cellI] += dm*YMix[LIQ]*YLiquid_[i];
-            }
-/*
-            // No mapping between solid components and carrier phase
-            forAll(YSolid_, i)
-            {
-                label gid = composition.localToGlobalCarrierId(SLD, i);
-                td.cloud().rhoTrans(gid)[cellI] += dm*YMix[SLD]*YSolid_[i];
-            }
-*/
-            td.cloud().UTrans()[cellI] += dm*U1;
-
-            td.cloud().hsTrans()[cellI] += dm*HsEff(td, pc, T1, idG, idL, idS);
-
-            td.cloud().addToMassPhaseChange(dm);
-        }
-    }
-
-
-    // Set new particle properties
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    else
-    {
-        this->Cp_ = CpEff(td, pc, T1, idG, idL, idS);
-        this->T_ = T1;
-        this->U_ = U1;
-
-        // Update particle density or diameter
-        if (td.cloud().constProps().constantVolume())
-        {
-            this->rho_ = mass1/this->volume();
-        }
-        else
-        {
-            this->d_ = cbrt(mass1/this->rho_*6.0/pi);
-        }
-    }
 }
 
 

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

@@ -267,21 +267,17 @@ void Foam::ReactingParcel<ParcelType>::calc
 
     // Define local properties at beginning of time step
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
     const scalar np0 = this->nParticle_;
     const scalar d0 = this->d_;
     const vector& U0 = this->U_;
-    const scalar rho0 = this->rho_;
     const scalar T0 = this->T_;
-    const scalar Cp0 = this->Cp_;
     const scalar mass0 = this->mass();
 
 
     // Calc surface values
-    // ~~~~~~~~~~~~~~~~~~~
     scalar Ts, rhos, mus, Prs, kappas;
     this->calcSurfaceValues(td, cellI, T0, Ts, rhos, mus, Prs, kappas);
-
-    // Reynolds number
     scalar Res = this->Re(U0, d0, rhos, mus);
 
 
@@ -291,16 +287,25 @@ void Foam::ReactingParcel<ParcelType>::calc
     // Explicit momentum source for particle
     vector Su = vector::zero;
 
+    // Linearised momentum source coefficient
+    scalar Spu = 0.0;
+
     // Momentum transfer from the particle to the carrier phase
     vector dUTrans = vector::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;
 
 
+    // 1. Compute models that contribute to mass transfer - U, T held constant
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
     // Phase change
     // ~~~~~~~~~~~~
 
@@ -338,96 +343,26 @@ void Foam::ReactingParcel<ParcelType>::calc
         Cs
     );
 
-    // Correct surface values due to emitted species
-    correctSurfaceValues(td, cellI, Ts, Cs, rhos, mus, Prs, kappas);
-    Res = this->Re(U0, d0, rhos, mus);
 
-    // Update particle component mass and mass fractions
+    // 2. Update the parcel properties due to change in mass
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
     scalarField dMass(dMassPC);
-
     scalar mass1 = updateMassFraction(mass0, dMass, Y_);
 
+    this->Cp_ = composition.Cp(0, Y_, pc_, T0);
 
-    // Heat transfer
+    // Update particle density or diameter
-    // ~~~~~~~~~~~~~
+    if (td.cloud().constProps().constantVolume())
-
-    // Calculate new particle temperature
-    scalar Sph = 0.0;
-    scalar T1 =
-        this->calcHeatTransfer
-        (
-            td,
-            dt,
-            cellI,
-            Res,
-            Prs,
-            kappas,
-            d0,
-            rho0,
-            T0,
-            Cp0,
-            NCpW,
-            Sh,
-            dhsTrans,
-            Sph
-        );
-
-
-    // Motion
-    // ~~~~~~
-
-    // Calculate new particle velocity
-    scalar Spu = 0.0;
-    vector U1 =
-        this->calcVelocity
-        (
-            td,
-            dt,
-            cellI,
-            Res,
-            mus,
-            d0,
-            U0,
-            rho0,
-            mass0,
-            Su,
-            dUTrans,
-            Spu
-        );
-
-
-    // Accumulate carrier phase source terms
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    if (td.cloud().solution().coupled())
     {
-        // Transfer mass lost to carrier mass and enthalpy sources
+        this->rho_ = mass1/this->volume();
-        forAll(dMass, i)
+    }
-        {
+    else
-            scalar dm = np0*dMass[i];
+    {
-            label gid = composition.localToGlobalCarrierId(0, i);
+        this->d_ = cbrt(mass1/this->rho_*6.0/pi);
-            scalar hs = composition.carrier().Hs(gid, T0);
-
-            td.cloud().rhoTrans(gid)[cellI] += dm;
-            td.cloud().UTrans()[cellI] += dm*U0;
-            td.cloud().hsTrans()[cellI] += dm*hs;
-        }
-
-        // Update momentum transfer
-        td.cloud().UTrans()[cellI] += np0*dUTrans;
-
-        // Update momentum transfer coefficient
-        td.cloud().UCoeff()[cellI] += np0*Spu;
-
-        // Update sensible enthalpy transfer
-        td.cloud().hsTrans()[cellI] += np0*dhsTrans;
-
-        // Update sensible enthalpy coefficient
-        td.cloud().hsCoeff()[cellI] += np0*Sph;
     }
 
-
     // Remove the particle when mass falls below minimum threshold
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     if (np0*mass1 < td.cloud().constProps().minParticleMass())
     {
         td.keepParticle = false;
@@ -441,36 +376,81 @@ void Foam::ReactingParcel<ParcelType>::calc
             {
                 scalar dmi = dm*Y_[i];
                 label gid = composition.localToGlobalCarrierId(0, i);
-                scalar hs = composition.carrier().Hs(gid, T1);
+                scalar hs = composition.carrier().Hs(gid, T0);
 
                 td.cloud().rhoTrans(gid)[cellI] += dmi;
                 td.cloud().hsTrans()[cellI] += dmi*hs;
             }
-            td.cloud().UTrans()[cellI] += dm*U1;
+            td.cloud().UTrans()[cellI] += dm*U0;
 
             td.cloud().addToMassPhaseChange(dm);
         }
+
+        return;
     }
 
+    // Correct surface values due to emitted species
+    correctSurfaceValues(td, cellI, Ts, Cs, rhos, mus, Prs, kappas);
+    Res = this->Re(U0, this->d_, rhos, mus);
 
-    // Set new particle properties
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    else
+    // 3. Compute heat- and momentum transfers
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Heat transfer
+    // ~~~~~~~~~~~~~
+
+    // Calculate new particle temperature
+    this->T_ =
+        this->calcHeatTransfer
+        (
+            td,
+            dt,
+            cellI,
+            Res,
+            Prs,
+            kappas,
+            NCpW,
+            Sh,
+            dhsTrans,
+            Sph
+        );
+
+    this->Cp_ = composition.Cp(0, Y_, pc_, T0);
+
+
+    // Motion
+    // ~~~~~~
+
+    // Calculate new particle velocity
+    this->U_ =
+        this->calcVelocity(td, dt, cellI, Res, mus, mass1, Su, dUTrans, Spu);
+
+
+    // 4. Accumulate carrier phase source terms
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    if (td.cloud().solution().coupled())
     {
-        this->Cp_ = composition.Cp(0, Y_, pc_, T1);
+        // Transfer mass lost to carrier mass, momentum and enthalpy sources
-        this->T_ = T1;
+        forAll(dMass, i)
-        this->U_ = U1;
+        {
+            scalar dm = np0*dMass[i];
+            label gid = composition.localToGlobalCarrierId(0, i);
+            scalar hs = composition.carrier().Hs(gid, T0);
 
-        // Update particle density or diameter
+            td.cloud().rhoTrans(gid)[cellI] += dm;
-        if (td.cloud().constProps().constantVolume())
+            td.cloud().UTrans()[cellI] += dm*U0;
-        {
+            td.cloud().hsTrans()[cellI] += dm*hs;
-            this->rho_ = mass1/this->volume();
-        }
-        else
-        {
-            this->d_ = cbrt(mass1/this->rho_*6.0/pi);
         }
+
+        // Update momentum transfer
+        td.cloud().UTrans()[cellI] += np0*dUTrans;
+        td.cloud().UCoeff()[cellI] += np0*Spu;
+
+        // Update sensible enthalpy transfer
+        td.cloud().hsTrans()[cellI] += np0*dhsTrans;
+        td.cloud().hsCoeff()[cellI] += np0*Sph;
     }
 }
 

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

@@ -172,21 +172,16 @@ void Foam::ThermoParcel<ParcelType>::calc
     // Define local properties at beginning of time step
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     const scalar np0 = this->nParticle_;
-    const scalar d0 = this->d_;
-    const vector U0 = this->U_;
-    const scalar rho0 = this->rho_;
-    const scalar T0 = this->T_;
-    const scalar Cp0 = this->Cp_;
     const scalar mass0 = this->mass();
 
 
     // Calc surface values
     // ~~~~~~~~~~~~~~~~~~~
     scalar Ts, rhos, mus, Pr, kappas;
-    calcSurfaceValues(td, cellI, T0, Ts, rhos, mus, Pr, kappas);
+    calcSurfaceValues(td, cellI, this->T_, Ts, rhos, mus, Pr, kappas);
 
     // Reynolds number
-    scalar Re = this->Re(U0, d0, rhos, mus);
+    scalar Re = this->Re(this->U_, this->d_, rhos, mus);
 
 
     // Sources
@@ -195,12 +190,18 @@ void Foam::ThermoParcel<ParcelType>::calc
     // Explicit momentum source for particle
     vector Su = vector::zero;
 
+    // Linearised momentum source coefficient
+    scalar Spu = 0.0;
+
     // Momentum transfer from the particle to the carrier phase
     vector dUTrans = vector::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;
 
@@ -212,8 +213,7 @@ void Foam::ThermoParcel<ParcelType>::calc
     scalar NCpW = 0.0;
 
     // Calculate new particle temperature
-    scalar Sph = 0.0;
+    this->T_ =
-    scalar T1 =
         this->calcHeatTransfer
         (
             td,
@@ -222,10 +222,6 @@ void Foam::ThermoParcel<ParcelType>::calc
             Re,
             Pr,
             kappas,
-            d0,
-            rho0,
-            T0,
-            Cp0,
             NCpW,
             Sh,
             dhsTrans,
@@ -237,23 +233,8 @@ void Foam::ThermoParcel<ParcelType>::calc
     // ~~~~~~
 
     // Calculate new particle velocity
-    scalar Spu = 0.0;
+    this->U_ =
-    vector U1 =
+        this->calcVelocity(td, dt, cellI, Re, mus, mass0, Su, dUTrans, Spu);
-        this->calcVelocity
-        (
-            td,
-            dt,
-            cellI,
-            Re,
-            mus,
-            d0,
-            U0,
-            rho0,
-            mass0,
-            Su,
-            dUTrans,
-            Spu
-        );
 
 
     //  Accumulate carrier phase source terms
@@ -272,11 +253,6 @@ void Foam::ThermoParcel<ParcelType>::calc
         // Update sensible enthalpy coefficient
         td.cloud().hsCoeff()[cellI] += np0*Sph;
     }
-
-    // Set new particle properties
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    this->U_ = U1;
-    T_ = T1;
 }
 
 
@@ -290,10 +266,6 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
     const scalar Re,
     const scalar Pr,
     const scalar kappa,
-    const scalar d,
-    const scalar rho,
-    const scalar T,
-    const scalar Cp,
     const scalar NCpW,
     const scalar Sh,
     scalar& dhsTrans,
@@ -302,9 +274,12 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
 {
     if (!td.cloud().heatTransfer().active())
     {
-        return T;
+        return T_;
     }
 
+    const scalar d = this->d();
+    const scalar rho = this->rho();
+
     // Calc heat transfer coefficient
     scalar htc = td.cloud().heatTransfer().htc(d, Re, Pr, kappa, NCpW);
 
@@ -313,7 +288,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
         return
             max
             (
-                T + dt*Sh/(this->volume(d)*rho*Cp),
+                T_ + dt*Sh/(this->volume(d)*rho*Cp_),
                 td.cloud().constProps().TMin()
             );
     }
@@ -322,7 +297,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
     const scalar As = this->areaS(d);
 
     scalar ap = Tc_ + Sh/As/htc;
-    scalar bp = 6.0*(Sh/As + htc*(Tc_ - T));
+    scalar bp = 6.0*(Sh/As + htc*(Tc_ - T_));
     if (td.cloud().radiation())
     {
         tetIndices tetIs = this->currentTetIndices();
@@ -330,20 +305,20 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
         const scalar sigma = physicoChemical::sigma.value();
         const scalar epsilon = td.cloud().constProps().epsilon0();
 
-        ap = (ap + epsilon*Gc/(4.0*htc))/(1.0 + epsilon*sigma*pow3(T)/htc);
+        ap = (ap + epsilon*Gc/(4.0*htc))/(1.0 + epsilon*sigma*pow3(T_)/htc);
-        bp += 6.0*(epsilon*(Gc/4.0 - sigma*pow4(T)));
+        bp += 6.0*(epsilon*(Gc/4.0 - sigma*pow4(T_)));
     }
-    bp /= rho*d*Cp*(ap - T) + ROOTVSMALL;
+    bp /= rho*d*Cp_*(ap - T_) + ROOTVSMALL;
 
     // Integrate to find the new parcel temperature
     IntegrationScheme<scalar>::integrationResult Tres =
-        td.cloud().TIntegrator().integrate(T, dt, ap*bp, bp);
+        td.cloud().TIntegrator().integrate(T_, dt, ap*bp, bp);
 
     scalar Tnew = max(Tres.value(), td.cloud().constProps().TMin());
 
     Sph = dt*htc*As;
 
-    dhsTrans += Sph*(0.5*(T + Tnew) - Tc_);
+    dhsTrans += Sph*(0.5*(T_ + Tnew) - Tc_);
 
     return Tnew;
 }

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

@@ -236,10 +236,6 @@ protected:
             const scalar Re,           // Reynolds number
             const scalar Pr,           // Prandtl number - surface
             const scalar kappa,        // Thermal conductivity - surface
-            const scalar d,            // diameter
-            const scalar rho,          // density
-            const scalar T,            // temperature
-            const scalar Cp,           // specific heat capacity
             const scalar NCpW,         // Sum of N*Cp*W of emission species
             const scalar Sh,           // explicit particle enthalpy source
             scalar& dhsTrans,          // sensible enthalpy transfer to carrier