diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index b75b20f8456997de81b12c300c2917d569d1a350..9d1a918bbc44e6ee9163059e4c66a8fe176f449d 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/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;
- vector U1 =
- calcVelocity
- (
- td,
- dt,
- cellI,
- Re,
- muc_,
- d0,
- U0,
- rho0,
- mass0,
- Su,
- dUTrans,
- Spu
- );
+ this->U_ = calcVelocity(td, dt, cellI, Re, muc_, 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();
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
index a937ba3f7a369d32f70e13f05e6dd2e90fe5a3ce..19b953b7e78184004f4250a81449431dcbfd7018 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
+++ b/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
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
index 664fac9781053758c3c890505fb8919c21e8462f..ae95fa21ef2a2eb8c14adf5ff0f0dfd022a2d3e8 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
+++ b/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
- this->correctSurfaceValues(td, cellI, Ts, Cs, rhos, mus, Prs, kappas);
- Res = this->Re(U0, d0, rhos, mus);
-
-
- // Update component mass fractions
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ // 2. Update the parcel properties due to change in mass
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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;
- scalar T1 =
+ this->T_ =
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;
- vector U1 =
- this->calcVelocity
- (
- td,
- dt,
- cellI,
- Res,
- mus,
- d0,
- U0,
- rho0,
- mass0,
- Su,
- dUTrans,
- Spu
- );
+ this->U_ =
+ this->calcVelocity(td, dt, cellI, Res, mus, mass1, 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);
- }
- }
}
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
index 71330144640a7a641b30a1d2b32251f56afc8102..21e7768aa777951685425719c4872de4de739332 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
+++ b/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,22 +343,65 @@ 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
- scalarField dMass(dMassPC);
+ // 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);
+
+ // 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(Y_, i)
+ {
+ scalar dmi = dm*Y_[i];
+ label gid = composition.localToGlobalCarrierId(0, i);
+ scalar hs = composition.carrier().Hs(gid, T0);
+
+ td.cloud().rhoTrans(gid)[cellI] += dmi;
+ td.cloud().hsTrans()[cellI] += dmi*hs;
+ }
+ 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);
+
+
+ // 3. Compute heat- and momentum transfers
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Heat transfer
// ~~~~~~~~~~~~~
// Calculate new particle temperature
- scalar Sph = 0.0;
- scalar T1 =
+ this->T_ =
this->calcHeatTransfer
(
td,
@@ -362,45 +410,29 @@ void Foam::ReactingParcel<ParcelType>::calc
Res,
Prs,
kappas,
- d0,
- rho0,
- T0,
- Cp0,
NCpW,
Sh,
dhsTrans,
Sph
);
+ this->Cp_ = composition.Cp(0, Y_, pc_, T0);
+
// 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
- );
+ this->U_ =
+ this->calcVelocity(td, dt, cellI, Res, mus, mass1, Su, dUTrans, Spu);
- // Accumulate carrier phase source terms
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ // 4. Accumulate carrier phase source terms
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
if (td.cloud().solution().coupled())
{
- // Transfer mass lost to carrier mass and enthalpy sources
+ // Transfer mass lost to carrier mass, momentum and enthalpy sources
forAll(dMass, i)
{
scalar dm = np0*dMass[i];
@@ -414,64 +446,12 @@ void Foam::ReactingParcel<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(Y_, i)
- {
- scalar dmi = dm*Y_[i];
- label gid = composition.localToGlobalCarrierId(0, i);
- scalar hs = composition.carrier().Hs(gid, T1);
-
- td.cloud().rhoTrans(gid)[cellI] += dmi;
- td.cloud().hsTrans()[cellI] += dmi*hs;
- }
- td.cloud().UTrans()[cellI] += dm*U1;
-
- td.cloud().addToMassPhaseChange(dm);
- }
- }
-
-
- // Set new particle properties
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- else
- {
- this->Cp_ = composition.Cp(0, Y_, pc_, T1);
- 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);
- }
- }
}
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
index dbc63e7c2e09a6c02cdac58a9fdaed9cbb752dfc..f7d37f3ad1af996e11f93e3852f60f631b7fc71b 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
+++ b/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;
- scalar T1 =
+ this->T_ =
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;
- vector U1 =
- this->calcVelocity
- (
- td,
- dt,
- cellI,
- Re,
- mus,
- d0,
- U0,
- rho0,
- mass0,
- Su,
- dUTrans,
- Spu
- );
+ this->U_ =
+ this->calcVelocity(td, dt, cellI, Re, mus, 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);
- bp += 6.0*(epsilon*(Gc/4.0 - sigma*pow4(T)));
+ 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 /= 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;
}
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
index c8fe14ad1a6eaa337f7c64f62d8c030bd3a7faf2..3b494e746e7326a929809e3a18b5ea0840b4ea06 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
+++ b/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