diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
index ecbb9266511717bc7c8bcad067324b16b11c8319..efbdea63898eae4773007327e2a2f1267e424666 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
@@ -155,21 +155,7 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
false
),
mesh_,
- dimensionedVector("zero", dimensionSet(1, 1, -1, 0, 0), vector::zero)
- ),
- UCoeff_
- (
- IOobject
- (
- this->name() + "UCoeff",
- this->db().time().timeName(),
- this->db(),
- IOobject::NO_READ,
- IOobject::NO_WRITE,
- false
- ),
- mesh_,
- dimensionedScalar("zero", dimensionSet(1, 0, -1, 0, 0), 0.0)
+ dimensionedVector("zero", dimMass*dimVelocity, vector::zero)
)
{}
@@ -187,7 +173,6 @@ template<class ParcelType>
void Foam::KinematicCloud<ParcelType>::resetSourceTerms()
{
UTrans_.field() = vector::zero;
- UCoeff_.field() = 0.0;
}
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
index 785b9819199ca23a269206ce4f8ba90756551161..50280469e2f9f024789eff4e996784c039d42c3a 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
@@ -165,9 +165,6 @@ class KinematicCloud
//- Momentum
DimensionedField<vector, volMesh> UTrans_;
- //- Coefficient for carrier phase U equation
- DimensionedField<scalar, volMesh> UCoeff_;
-
// Private Member Functions
@@ -309,14 +306,8 @@ public:
//- Return reference to momentum source
inline DimensionedField<vector, volMesh>& UTrans();
- //- Coefficient for carrier phase U equation
- inline DimensionedField<scalar, volMesh>& UCoeff();
-
//- Return tmp momentum source term - fully explicit
- inline tmp<DimensionedField<vector, volMesh> > SU1() const;
-
- //- Return tmp momentum source term - semi-implicit
- inline tmp<fvVectorMatrix> SU2(volVectorField& U) const;
+ inline tmp<DimensionedField<vector, volMesh> > SU() const;
// Check
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
index dde35a4298fc380ddaaea4fc5abfe4c2063ecb1f..8c44470c44bb867583212eaadfea65aa51666e6b 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
@@ -203,25 +203,17 @@ Foam::KinematicCloud<ParcelType>::UTrans()
}
-template<class ParcelType>
-inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
-Foam::KinematicCloud<ParcelType>::UCoeff()
-{
- return UCoeff_;
-}
-
-
template<class ParcelType>
inline Foam::tmp<Foam::DimensionedField<Foam::vector, Foam::volMesh> >
-Foam::KinematicCloud<ParcelType>::SU1() const
+Foam::KinematicCloud<ParcelType>::SU() const
{
- tmp<DimensionedField<vector, volMesh> > tSU1
+ tmp<DimensionedField<vector, volMesh> > tSU
(
new DimensionedField<vector, volMesh>
(
IOobject
(
- this->name() + "SU1",
+ this->name() + "SU",
this->db().time().timeName(),
this->mesh(),
IOobject::NO_READ,
@@ -237,29 +229,10 @@ Foam::KinematicCloud<ParcelType>::SU1() const
)
);
- vectorField& SU1 = tSU1().field();
- SU1 = UTrans_/(mesh_.V()*this->db().time().deltaT());
-
- return tSU1;
-}
-
-
-template<class ParcelType>
-inline Foam::tmp<Foam::fvVectorMatrix>
-Foam::KinematicCloud<ParcelType>::SU2(volVectorField& U) const
-{
- if (debug)
- {
- Info<< "UTrans min/max = "
- << min(UTrans_) << ", " << max(UTrans_) << endl;
- Info<< "UCoeff min/max = "
- << min(UCoeff_) << ", " << max(UCoeff_) << endl;
- }
+ vectorField& SU = tSU().field();
+ SU = UTrans_/(mesh_.V()*this->db().time().deltaT());
- return
- UTrans_/(mesh_.V()*this->db().time().deltaT())
- - fvm::Sp(UCoeff_/mesh_.V(), U)
- + UCoeff_/mesh_.V()*U;
+ return tSU;
}
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H
index 6fe5bf68a65a9a7db9959fb02cc18f5fff2adb7a..0cc8bfc96b9812fe45c400c98e968e3272316a32 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H
@@ -180,15 +180,13 @@ public:
inline PtrList<DimensionedField<scalar, volMesh> >&
rhoTrans();
- //- Return tmp mass source for field i
- // Fully explicit
+ //- Return tmp mass source for field i - fully explicit
inline tmp<DimensionedField<scalar, volMesh> >
- Srho1(const label i) const;
+ Srho(const label i) const;
//- Return tmp total mass source for carrier phase
- // Fully explicit
- inline tmp<DimensionedField<scalar, volMesh> >
- Srho1() const;
+ // - fully explicit
+ inline tmp<DimensionedField<scalar, volMesh> > Srho() const;
// Check
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudI.H b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudI.H
index 1414b6b16d78e2c2fe141c5144f8836a9456acee..98d9eee8feff5ea659c9b51d58318a3cf950494b 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloudI.H
@@ -92,7 +92,7 @@ Foam::ReactingCloud<ParcelType>::rhoTrans()
template<class ParcelType>
inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
-Foam::ReactingCloud<ParcelType>::Srho1(const label i) const
+Foam::ReactingCloud<ParcelType>::Srho(const label i) const
{
return rhoTrans_[i]/(this->db().time().deltaT()*this->mesh().V());
}
@@ -100,7 +100,7 @@ Foam::ReactingCloud<ParcelType>::Srho1(const label i) const
template<class ParcelType>
inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
-Foam::ReactingCloud<ParcelType>::Srho1() const
+Foam::ReactingCloud<ParcelType>::Srho() const
{
tmp<DimensionedField<scalar, volMesh> > trhoTrans
(
diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
index 6872f2428eadeb99b6bbcbaa135624c1460312fc..62c926f1eaca456af40190185489245a4df6f4a7 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
@@ -115,20 +115,6 @@ Foam::ThermoCloud<ParcelType>::ThermoCloud
),
this->mesh(),
dimensionedScalar("zero", dimensionSet(1, 2, -2, 0, 0), 0.0)
- ),
- hCoeff_
- (
- IOobject
- (
- this->name() + "hCoeff",
- this->db().time().timeName(),
- this->db(),
- IOobject::NO_READ,
- IOobject::NO_WRITE,
- false
- ),
- this->mesh(),
- dimensionedScalar("zero", dimensionSet(1, 2, -3, -1, 0), 0.0)
)
{}
@@ -147,7 +133,6 @@ void Foam::ThermoCloud<ParcelType>::resetSourceTerms()
{
KinematicCloud<ParcelType>::resetSourceTerms();
hTrans_.field() = 0.0;
- hCoeff_.field() = 0.0;
}
diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H
index 9789bd0f5f08258613879bd186423def7c15ae20..7cdc7b12636a3c0bc54aa030878b3eb4e691a913 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H
@@ -102,9 +102,6 @@ class ThermoCloud
//- Enthalpy transfer
DimensionedField<scalar, volMesh> hTrans_;
- //- Coefficient for carrier phase h equation
- DimensionedField<scalar, volMesh> hCoeff_;
-
// Private Member Functions
@@ -179,14 +176,8 @@ public:
//- Return reference to enthalpy source
inline DimensionedField<scalar, volMesh>& hTrans();
- //- Coefficient for carrier phase h equation
- inline DimensionedField<scalar, volMesh>& hCoeff();
-
//- return tmp enthalpy source term - fully explicit
- inline tmp<DimensionedField<scalar, volMesh> > Sh1() const;
-
- //- Return tmp enthalpy source term - semi-implicit
- inline tmp<fvScalarMatrix> Sh2(volScalarField& h) const;
+ inline tmp<DimensionedField<scalar, volMesh> > Sh() const;
// Radiation - overrides thermoCloud virtual abstract members
diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
index 2f3df96f41b9db4483a9a8fa3b241df58e49f159..d85d043ec55216d5c024af352acddea9cc61d9ff 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H
@@ -83,25 +83,17 @@ Foam::ThermoCloud<ParcelType>::hTrans()
}
-template<class ParcelType>
-inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
-Foam::ThermoCloud<ParcelType>::hCoeff()
-{
- return hCoeff_;
-}
-
-
template<class ParcelType>
inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
-Foam::ThermoCloud<ParcelType>::Sh1() const
+Foam::ThermoCloud<ParcelType>::Sh() const
{
- tmp<DimensionedField<scalar, volMesh> > tSh1
+ tmp<DimensionedField<scalar, volMesh> > tSh
(
new DimensionedField<scalar, volMesh>
(
IOobject
(
- this->name() + "Sh1",
+ this->name() + "Sh",
this->db().time().timeName(),
this->mesh(),
IOobject::NO_READ,
@@ -118,30 +110,10 @@ Foam::ThermoCloud<ParcelType>::Sh1() const
)
);
- scalarField& Sh1 = tSh1().field();
- Sh1 = hTrans_/(this->mesh().V()*this->db().time().deltaT());
-
- return tSh1;
-}
-
-
-template<class ParcelType>
-inline Foam::tmp<Foam::fvScalarMatrix>
-Foam::ThermoCloud<ParcelType>::Sh2(volScalarField& h) const
-{
- const volScalarField cp = carrierThermo_.Cp();
-
- if (debug)
- {
- Info<< "hTrans min/max = "
- << min(hTrans_) << ", " << max(hTrans_) << endl;
- Info<< "hCoeff min/max = "
- << min(hCoeff_) << ", " << max(hCoeff_) << endl;
- }
+ scalarField& Sh = tSh().field();
+ Sh = hTrans_/(this->mesh().V()*this->db().time().deltaT());
- return hTrans_/(this->mesh().V()*this->db().time().deltaT())
- - fvm::Sp(hCoeff_/(cp*this->mesh().V()), h)
- + hCoeff_/(cp*this->mesh().V())*h;
+ return tSh;
}
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index ad564e6d89232c59e7b449e41d8533ebcc1f4d5f..69cc25869e5a29075863f01a5c37fe6f6544ad56 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -92,15 +92,6 @@ void Foam::KinematicParcel<ParcelType>::calc
const scalar rho0 = rho_;
const scalar mass0 = mass();
- const polyMesh& mesh = this->cloud().pMesh();
-
-
- // Initialise transfer terms
- // ~~~~~~~~~~~~~~~~~~~~~~~~~
-
- // Momentum
- vector dUTrans = vector::zero;
-
// Motion
// ~~~~~~
@@ -109,9 +100,7 @@ void Foam::KinematicParcel<ParcelType>::calc
vector Fx = vector::zero;
// Calculate new particle velocity
- scalar Cud = 0.0;
- vector U1 =
- calcVelocity(td, dt, cellI, d0, U0, rho0, mass0, Fx, Cud, dUTrans);
+ vector U1 = calcVelocity(td, dt, cellI, d0, U0, rho0, mass0, Fx);
// Accumulate carrier phase source terms
@@ -119,10 +108,7 @@ void Foam::KinematicParcel<ParcelType>::calc
if (td.cloud().coupled())
{
// Update momentum transfer
- td.cloud().UTrans()[cellI] += np0*dUTrans;
-
- // Coefficient to be applied in carrier phase momentum coupling
- td.cloud().UCoeff()[cellI] += np0*mass0*Cud;
+ td.cloud().UTrans()[cellI] += np0*mass0*(U0 - U1);
}
@@ -143,15 +129,13 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
const vector& U,
const scalar rho,
const scalar mass,
- const vector& Fx,
- scalar& Cud,
- vector& dUTrans
+ const vector& Fx
) const
{
const polyMesh& mesh = this->cloud().pMesh();
// Return linearised term from drag model
- Cud = td.cloud().drag().Cu(U - Uc_, d, rhoc_, rho, muc_);
+ scalar Cud = td.cloud().drag().Cu(U - Uc_, d, rhoc_, rho, muc_);
// Calculate particle forces
vector Ftot = td.cloud().forces().calc(cellI, dt, rhoc_, rho, Uc_, U);
@@ -169,10 +153,6 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
// Apply correction to velocity for reduced-D cases
meshTools::constrainDirection(mesh, mesh.solutionD(), Unew);
- // Calculate the momentum transfer to the continuous phase
- // - do not include gravity impulse
- 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 d5ccfd84b2f662c0c4534437070c1cbff2342742..7a57dd5fb017b90c0f63dcb943fa202d65b26c29 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
@@ -243,9 +243,7 @@ protected:
const vector& U, // velocity
const scalar rho, // density
const scalar mass, // mass
- const vector& Fx, // additional forces
- scalar& Cud, // linearised drag term coeff
- vector& dUTrans // momentum transfer to carrier phase
+ const vector& Fx // additional forces
) const;
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
index 7692c2011f3abc07f0b377f55f24041ce5683916..3d67afa89b7ca401d9c5b3aaf9740ca1cfd1643e 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
@@ -32,10 +32,10 @@ template<class ParcelType>
const Foam::label Foam::ReactingMultiphaseParcel<ParcelType>::GAS(0);
template<class ParcelType>
-const Foam::label Foam::ReactingMultiphaseParcel<ParcelType>::LIQUID(1);
+const Foam::label Foam::ReactingMultiphaseParcel<ParcelType>::LIQ(1);
template<class ParcelType>
-const Foam::label Foam::ReactingMultiphaseParcel<ParcelType>::SOLID(2);
+const Foam::label Foam::ReactingMultiphaseParcel<ParcelType>::SLD(2);
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
@@ -54,8 +54,8 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::cpEff
{
return
this->Y_[GAS]*td.cloud().composition().cp(idG, YGas_, p, T)
- + this->Y_[LIQUID]*td.cloud().composition().cp(idL, YLiquid_, p, T)
- + this->Y_[SOLID]*td.cloud().composition().cp(idS, YSolid_, p, T);
+ + this->Y_[LIQ]*td.cloud().composition().cp(idL, YLiquid_, p, T)
+ + this->Y_[SLD]*td.cloud().composition().cp(idS, YSolid_, p, T);
}
@@ -73,8 +73,27 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::HEff
{
return
this->Y_[GAS]*td.cloud().composition().H(idG, YGas_, p, T)
- + this->Y_[LIQUID]*td.cloud().composition().H(idL, YLiquid_, p, T)
- + this->Y_[SOLID]*td.cloud().composition().H(idS, YSolid_, p, T);
+ + this->Y_[LIQ]*td.cloud().composition().H(idL, YLiquid_, p, T)
+ + this->Y_[SLD]*td.cloud().composition().H(idS, YSolid_, p, T);
+}
+
+
+template<class ParcelType>
+template<class TrackData>
+Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::LEff
+(
+ TrackData& td,
+ const scalar p,
+ const scalar T,
+ const label idG,
+ const label idL,
+ const label idS
+) const
+{
+ return
+ this->Y_[GAS]*td.cloud().composition().L(idG, YGas_, p, T)
+ + this->Y_[LIQ]*td.cloud().composition().L(idL, YLiquid_, p, T)
+ + this->Y_[SLD]*td.cloud().composition().L(idS, YSolid_, p, T);
}
@@ -94,14 +113,14 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::updateMassFractions
scalar dMassSolidTot = sum(dMassSolid);
this->updateMassFraction(mass0*YMix[GAS], dMassGas, YGas_);
- this->updateMassFraction(mass0*YMix[LIQUID], dMassLiquid, YLiquid_);
- this->updateMassFraction(mass0*YMix[SOLID], dMassSolid, YSolid_);
+ this->updateMassFraction(mass0*YMix[LIQ], dMassLiquid, YLiquid_);
+ this->updateMassFraction(mass0*YMix[SLD], dMassSolid, YSolid_);
scalar massNew = mass0 - (dMassGasTot + dMassLiquidTot + dMassSolidTot);
YMix[GAS] = (mass0*YMix[GAS] - dMassGasTot)/massNew;
- YMix[LIQUID] = (mass0*YMix[LIQUID] - dMassLiquidTot)/massNew;
- YMix[SOLID] = 1.0 - YMix[GAS] - YMix[LIQUID];
+ YMix[LIQ] = (mass0*YMix[LIQ] - dMassLiquidTot)/massNew;
+ YMix[SLD] = 1.0 - YMix[GAS] - YMix[LIQ];
return massNew;
}
@@ -148,32 +167,18 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
const label idL = td.cloud().composition().idLiquid();
const label idS = td.cloud().composition().idSolid();
+ // Intial ethalpy state
+ scalar H0H = HEff(td, pc, T0, idG, idL, idS);
+ scalar H0L = LEff(td, pc, T0, idG, idL, idS);
+ scalar H0 = H0H - H0L;
- // Initialise transfer terms
- // ~~~~~~~~~~~~~~~~~~~~~~~~~
- // Momentum
- vector dUTrans = vector::zero;
-
- // Enthalpy
- scalar dhTrans = 0.0;
+ // Phase change in liquid phase
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Mass transfer due to phase change
scalarField dMassPC(YLiquid_.size(), 0.0);
- // Mass transfer due to devolatilisation
- scalarField dMassDV(YGas_.size(), 0.0);
-
- // Change in carrier phase composition due to surface reactions
- scalarField dMassSRGas(YGas_.size(), 0.0);
- scalarField dMassSRLiquid(YLiquid_.size(), 0.0);
- scalarField dMassSRSolid(YSolid_.size(), 0.0);
- scalarField dMassSRCarrier(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
@@ -185,7 +190,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
T0,
U0,
idL,
- YMix[LIQUID],
+ YMix[LIQ],
YLiquid_,
dMassPC
);
@@ -194,6 +199,9 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
// Devolatilisation
// ~~~~~~~~~~~~~~~~
+ // Mass transfer due to devolatilisation
+ scalarField dMassDV(YGas_.size(), 0.0);
+
// Return enthalpy source and calc mass transfer due to devolatilisation
scalar ShDV =
calcDevolatilisation
@@ -214,8 +222,14 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
// Surface reactions
// ~~~~~~~~~~~~~~~~~
+ // Change in carrier phase composition due to surface reactions
+ scalarField dMassSRGas(YGas_.size(), 0.0);
+ scalarField dMassSRLiquid(YLiquid_.size(), 0.0);
+ scalarField dMassSRSolid(YSolid_.size(), 0.0);
+ scalarField dMassSRCarrier(td.cloud().gases().size(), 0.0);
+
// Return enthalpy source and calc mass transfer(s) due to surface reaction
- scalar ShSR =
+ scalar HReaction =
calcSurfaceReactions
(
td,
@@ -233,34 +247,14 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
dMassSRGas,
dMassSRLiquid,
dMassSRSolid,
- dMassSRCarrier,
- dhTrans
+ dMassSRCarrier
);
+ // Heat of reaction split between component retained by particle
+ const scalar ShSR = td.constProps().hRetentionCoeff()*HReaction;
- // Heat transfer
- // ~~~~~~~~~~~~~
-
- // Total enthalpy source
- scalar Sh = ShPC + ShDV + ShSR;
-
- // Calculate new particle temperature
- scalar htc = 0.0;
- scalar T1 =
- calcHeatTransfer
- (
- td,
- dt,
- cellI,
- d0,
- U0,
- rho0,
- T0,
- cp0,
- Sh,
- htc,
- dhTrans
- );
+ // ...and component added to the carrier phase
+ const scalar ShSRc = (1.0 - td.constProps().hRetentionCoeff())*HReaction;
// Update component mass fractions
@@ -274,6 +268,22 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
updateMassFractions(mass0, dMassGas, dMassLiquid, dMassSolid);
+ // Heat transfer
+ // ~~~~~~~~~~~~~
+
+ // Total enthalpy source
+ scalar Sh = ShPC + ShDV + ShSR;
+
+ // Calculate new particle temperature
+ scalar T1 = calcHeatTransfer(td, dt, cellI, d0, U0, rho0, T0, cp0, Sh);
+
+ // Calculate new enthalpy state based on updated composition at new
+ // temperature
+ scalar H1H = HEff(td, pc, T1, idG, idL, idS);
+ scalar H1L = LEff(td, pc, T1, idG, idL, idS);
+ scalar H1 = H1H - H1L;
+
+
// Motion
// ~~~~~~
@@ -281,21 +291,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
vector Fx = vector::zero;
// Calculate new particle velocity
- scalar Cud = 0.0;
- vector U1 =
- calcVelocity
- (
- td,
- dt,
- cellI,
- d0,
- U0,
- rho0,
- 0.5*(mass0 + mass1),
- Fx,
- Cud,
- dUTrans
- );
+ vector U1 = calcVelocity(td, dt, cellI, d0, U0, rho0, mass0, Fx);
// Accumulate carrier phase source terms
@@ -311,13 +307,13 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
}
forAll(YLiquid_, i)
{
- label id = td.cloud().composition().localToGlobalGasId(LIQUID, i);
+ label id = td.cloud().composition().localToGlobalGasId(LIQ, i);
td.cloud().rhoTrans(id)[cellI] += np0*dMassLiquid[i];
}
// // No mapping between solid components and carrier phase
// forAll(YSolid_, i)
// {
-// label id = td.cloud().composition().localToGlobalGasId(SOLID, i);
+// label id = td.cloud().composition().localToGlobalGasId(SLD, i);
// td.cloud().rhoTrans(id)[cellI] += np0*dMassSolid[i];
// }
forAll(dMassSRCarrier, i)
@@ -326,16 +322,10 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
}
// Update momentum transfer
- td.cloud().UTrans()[cellI] += np0*dUTrans;
-
- // Coefficient to be applied in carrier phase momentum coupling
- td.cloud().UCoeff()[cellI] += np0*mass0*Cud;
+ td.cloud().UTrans()[cellI] += np0*(mass0*U0 - mass1*U1);
// Update enthalpy transfer
- td.cloud().hTrans()[cellI] += np0*dhTrans;
-
- // Coefficient to be applied in carrier phase enthalpy coupling
- td.cloud().hCoeff()[cellI] += np0*htc*this->areaS();
+ td.cloud().hTrans()[cellI] += np0*(mass0*H0 - (mass1*H1 + ShSRc));
}
@@ -357,29 +347,20 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
forAll(YLiquid_, i)
{
label id =
- td.cloud().composition().localToGlobalGasId(LIQUID, i);
+ td.cloud().composition().localToGlobalGasId(LIQ, i);
td.cloud().rhoTrans(id)[cellI] +=
- np0
- *mass1
- *YMix[LIQUID]
- *YLiquid_[i];
+ np0*mass1*YMix[LIQ]*YLiquid_[i];
}
// // No mapping between solid components and carrier phase
// forAll(YSolid_, i)
// {
// label id =
-// td.cloud().composition().localToGlobalGasId(SOLID, i);
+// td.cloud().composition().localToGlobalGasId(SLD, i);
// td.cloud().rhoTrans(id)[cellI] +=
-// np0
-// *mass1
-// *YMix[SOLID]
-// *YSolid_[i];
+// np0*mass1*YMix[SLD]*YSolid_[i];
// }
- td.cloud().hTrans()[cellI] +=
- np0
- *mass1
- *HEff(td, pc, T1, idG, idL, idS);
+ td.cloud().hTrans()[cellI] += np0*mass1*H1;
td.cloud().UTrans()[cellI] += np0*mass1*U1;
}
}
@@ -421,7 +402,7 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
const scalarField& YVolatile,
bool& canCombust,
scalarField& dMassDV
-)
+) const
{
// Check that model is active, and that the parcel temperature is
// within necessary limits for devolatilisation to occur
@@ -448,10 +429,7 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
);
// Volatile mass transfer - equal components of each volatile specie
- forAll(YVolatile, i)
- {
- dMassDV[i] = YVolatile[i]*dMassTot;
- }
+ dMassDV = YVolatile*dMassTot;
td.cloud().addToMassDevolatilisation(this->nParticle_*dMassTot);
@@ -478,9 +456,8 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
scalarField& dMassSRGas,
scalarField& dMassSRLiquid,
scalarField& dMassSRSolid,
- scalarField& dMassSRCarrier,
- scalar& dhTrans
-)
+ scalarField& dMassSRCarrier
+) const
{
// Check that model is active
if (!td.cloud().surfaceReaction().active() || !canCombust)
@@ -489,18 +466,15 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
}
// Update surface reactions
- const scalar pc = this->pc_;
- const scalar Tc = this->Tc_;
- const scalar rhoc = this->rhoc_;
const scalar HReaction = td.cloud().surfaceReaction().calculate
(
dt,
cellI,
d,
T,
- Tc,
- pc,
- rhoc,
+ this->Tc_,
+ this->pc_,
+ this->rhoc_,
mass,
YGas,
YLiquid,
@@ -519,27 +493,7 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
*(sum(dMassSRGas) + sum(dMassSRLiquid) + sum(dMassSRSolid))
);
- // Add enthalpy of consumed components to the carrier phase enthalpy
- // transfer
- const label idG = td.cloud().composition().idGas();
- const label idL = td.cloud().composition().idLiquid();
- const label idS = td.cloud().composition().idSolid();
- scalar dhGas =
- sum(dMassSRGas)*td.cloud().composition().H(idG, YGas, pc, T);
-
- scalar dhLiquid =
- sum(dMassSRLiquid)*td.cloud().composition().H(idL, YLiquid, pc, T);
-
- scalar dhSolid =
- sum(dMassSRSolid)*td.cloud().composition().H(idS, YSolid, pc, T);
-
- dhTrans += dhGas + dhLiquid + dhSolid;
-
- // Heat of reaction divided between particle and carrier phase by the
- // fraction fh and (1 - fh)
- dhTrans += (1.0 - td.constProps().hRetentionCoeff())*HReaction;
-
- return td.constProps().hRetentionCoeff()*HReaction;
+ return HReaction;
}
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
index f1bace3a8b973956d55a5b30d0b3173661e0415d..b5f45310cd4feb610d5ec9da1c4bc23649157484 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
@@ -71,11 +71,11 @@ public:
// IDs of phases in ReacingParcel phase list (Y)
static const label GAS;
- static const label LIQUID;
- static const label SOLID;
+ static const label LIQ;
+ static const label SLD;
- //- Class to hold reacting particle constant properties
+ //- Class to hold reacting multiphase particle constant properties
class constantProperties
:
public ReactingParcel<ParcelType>::constantProperties
@@ -177,6 +177,18 @@ private:
const label idS
) const;
+ //- Return the mixture effective latent heat
+ template<class TrackData>
+ scalar LEff
+ (
+ TrackData& td,
+ const scalar p,
+ const scalar T,
+ const label idG,
+ const label idL,
+ const label idS
+ ) const;
+
//- Update the mass fractions (Y, YGas, YLiquid, YSolid)
scalar updateMassFractions
(
@@ -223,7 +235,7 @@ protected:
const scalarField& YVolatile, // volatile component mass fractions
bool& canCombust, // 'can combust' flag
scalarField& dMassDV // mass transfer - local to particle
- );
+ ) const;
//- Calculate surface reactions
template<class TrackData>
@@ -237,16 +249,15 @@ protected:
const scalar mass, // mass
const bool canCombust, // 'can combust' flag
const scalarField& dMassVolatile, // mass transfer of volatiles
- const scalarField& YMix, // mixture mass fractions
+ const scalarField& YMix, // mixture mass fractions
const scalarField& YGas, // gas-phase mass fractions
const scalarField& YLiquid,// liquid-phase mass fractions
const scalarField& YSolid, // solid-phase mass fractions
scalarField& dMassSRGas, // gas-phase mass transfer - local
scalarField& dMassSRLiquid,// liquid-phase mass transfer - local
scalarField& dMassSRSolid, // solid-phase mass transfer - local
- scalarField& dMassSRCarrier,// carrier phase mass transfer
- scalar& dhTrans // enthalpy transfer to carrier phase
- );
+ scalarField& dMassSRCarrier// carrier phase mass transfer
+ ) const;
public:
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelIO.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelIO.C
index 64829732e8fee556a0c7b9ca0c5e260e83f08ab8..b37f7c715e27ea0272a038bef14676f10384e9cd 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelIO.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelIO.C
@@ -70,8 +70,8 @@ Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
// scale the mass fractions
const scalarField& YMix = this->Y_;
YGas_ /= YMix[GAS] + ROOTVSMALL;
- YLiquid_ /= YMix[LIQUID] + ROOTVSMALL;
- YSolid_ /= YMix[SOLID] + ROOTVSMALL;
+ YLiquid_ /= YMix[LIQ] + ROOTVSMALL;
+ YSolid_ /= YMix[SLD] + ROOTVSMALL;
}
// Check state of Istream
@@ -153,7 +153,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::readFields
forAllIter(typename Cloud<ParcelType>, c, iter)
{
ReactingMultiphaseParcel<ParcelType>& p = iter();
- p.YLiquid_[j] = YLiquid[i++]/(p.Y()[LIQUID] + ROOTVSMALL);
+ p.YLiquid_[j] = YLiquid[i++]/(p.Y()[LIQ] + ROOTVSMALL);
}
}
// Populate YSolid for each parcel
@@ -172,7 +172,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::readFields
forAllIter(typename Cloud<ParcelType>, c, iter)
{
ReactingMultiphaseParcel<ParcelType>& p = iter();
- p.YSolid_[j] = YSolid[i++]/(p.Y()[SOLID] + ROOTVSMALL);
+ p.YSolid_[j] = YSolid[i++]/(p.Y()[SLD] + ROOTVSMALL);
}
}
}
@@ -235,7 +235,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::writeFields
forAllConstIter(typename Cloud<ParcelType>, c, iter)
{
const ReactingMultiphaseParcel<ParcelType>& p0 = iter();
- YLiquid[i++] = p0.YLiquid()[j]*p0.Y()[LIQUID];
+ YLiquid[i++] = p0.YLiquid()[j]*p0.Y()[LIQ];
}
YLiquid.write();
@@ -259,7 +259,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::writeFields
forAllConstIter(typename Cloud<ParcelType>, c, iter)
{
const ReactingMultiphaseParcel<ParcelType>& p0 = iter();
- YSolid[i++] = p0.YSolid()[j]*p0.Y()[SOLID];
+ YSolid[i++] = p0.YSolid()[j]*p0.Y()[SLD];
}
YSolid.write();
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
index 4529d4e07798d4c55a886a9628be55e10541cd3e..1db979c6da654a75e2722d8a271a7abea384a2ca 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
@@ -102,23 +102,18 @@ void Foam::ReactingParcel<ParcelType>::calc
const scalar cp0 = this->cp_;
const scalar mass0 = this->mass();
+ // Intial ethalpy state
+ scalar H0H = td.cloud().composition().H(0, Y_, pc_, T0);
+ scalar H0L = td.cloud().composition().L(0, Y_, pc_, T0);
+ scalar H0 = H0H - H0L;
- // Intialise transfer terms
- // ~~~~~~~~~~~~~~~~~~~~~~~~
- // Momentum
- vector dUTrans = vector::zero;
-
- // Enthalpy
- scalar dhTrans = 0.0;
+ // Phase change
+ // ~~~~~~~~~~~~
// Mass transfer due to phase change
scalarField dMassPC(Y_.size(), 0.0);
-
- // Phase change
- // ~~~~~~~~~~~~
-
// Return enthalpy source and calc mass transfer due to phase change
scalar ShPC =
calcPhaseChange(td, dt, cellI, d0, T0, U0, 0, 1.0, Y_, dMassPC);
@@ -126,54 +121,30 @@ void Foam::ReactingParcel<ParcelType>::calc
// Update particle component mass fractions
updateMassFraction(mass0, dMassPC, Y_);
+ // Update mass
+ scalar mass1 = mass0 - sum(dMassPC);
+
// Heat transfer
// ~~~~~~~~~~~~~
// Calculate new particle temperature
- scalar htc = 0.0;
- scalar T1 =
- calcHeatTransfer
- (
- td,
- dt,
- cellI,
- d0,
- U0,
- rho0,
- T0,
- cp0,
- ShPC,
- htc,
- dhTrans
- );
+ scalar T1 = calcHeatTransfer(td, dt, cellI, d0, U0, rho0, T0, cp0, ShPC);
+
+ // Calculate new enthalpy state
+ scalar H1H = td.cloud().composition().H(0, Y_, pc_, T1);
+ scalar H1L = td.cloud().composition().L(0, Y_, pc_, T1);
+ scalar H1 = H1H - H1L;
// Motion
// ~~~~~~
- // Update mass
- scalar mass1 = mass0 - sum(dMassPC);
-
// No additional forces
vector Fx = vector::zero;
// Calculate new particle velocity
- scalar Cud = 0.0;
- vector U1 =
- calcVelocity
- (
- td,
- dt,
- cellI,
- d0,
- U0,
- rho0,
- 0.5*(mass0 + mass1),
- Fx,
- Cud,
- dUTrans
- );
+ vector U1 = calcVelocity(td, dt, cellI, d0, U0, rho0, mass0, Fx);
// Accumulate carrier phase source terms
@@ -188,16 +159,10 @@ void Foam::ReactingParcel<ParcelType>::calc
}
// Update momentum transfer
- td.cloud().UTrans()[cellI] += np0*dUTrans;
-
- // Coefficient to be applied in carrier phase momentum coupling
- td.cloud().UCoeff()[cellI] += np0*mass0*Cud;
+ td.cloud().UTrans()[cellI] += np0*(mass0*U0 - mass1*U1);
// Update enthalpy transfer
- td.cloud().hTrans()[cellI] += np0*dhTrans;
-
- // Coefficient to be applied in carrier phase enthalpy coupling
- td.cloud().hCoeff()[cellI] += np0*htc*this->areaS();
+ td.cloud().hTrans()[cellI] += np0*(mass0*H0 - mass1*H1);
}
@@ -216,8 +181,7 @@ void Foam::ReactingParcel<ParcelType>::calc
td.cloud().rhoTrans(id)[cellI] += np0*mass1*Y_[i];
}
td.cloud().UTrans()[cellI] += np0*mass1*U1;
- scalar HEff = td.cloud().composition().H(0, Y_, pc_, T1);
- td.cloud().hTrans()[cellI] += np0*mass1*HEff;
+ td.cloud().hTrans()[cellI] += np0*mass1*H1;
}
}
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
index 8b5a722d30ab7e6eef16e77357baee370aa79ecb..4d4d469d124b1ab7112185683e7b4d5e07c8cb4b 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
@@ -84,12 +84,8 @@ void Foam::ThermoParcel<ParcelType>::calc
const scalar cp0 = this->cp_;
const scalar mass0 = this->mass();
-
- // Initialise transfer terms
- // ~~~~~~~~~~~~~~~~~~~~~~~~~
-
- vector dUTrans = vector::zero;
- scalar dhTrans = 0.0;
+ // Initial enthalpy state
+ scalar H0 = cp0*T0;
// Heat transfer
@@ -99,22 +95,10 @@ void Foam::ThermoParcel<ParcelType>::calc
scalar Sh = 0.0;
// Calculate new particle velocity
- scalar htc = 0.0;
- scalar T1 =
- calcHeatTransfer
- (
- td,
- dt,
- cellI,
- d0,
- U0,
- rho0,
- T0,
- cp0,
- Sh,
- htc,
- dhTrans
- );
+ scalar T1 = calcHeatTransfer(td, dt, cellI, d0, U0, rho0, T0, cp0, Sh);
+
+ // Calculate new enthalpy state
+ scalar H1 = cp0 - T1;
// Motion
@@ -124,9 +108,7 @@ void Foam::ThermoParcel<ParcelType>::calc
vector Fx = vector::zero;
// Calculate new particle velocity
- scalar Cud = 0.0;
- vector U1 =
- calcVelocity(td, dt, cellI, d0, U0, rho0, mass0, Fx, Cud, dUTrans);
+ vector U1 = calcVelocity(td, dt, cellI, d0, U0, rho0, mass0, Fx);
// Accumulate carrier phase source terms
@@ -134,16 +116,10 @@ void Foam::ThermoParcel<ParcelType>::calc
if (td.cloud().coupled())
{
// Update momentum transfer
- td.cloud().UTrans()[cellI] += np0*dUTrans;
-
- // Coefficient to be applied in carrier phase momentum coupling
- td.cloud().UCoeff()[cellI] += np0*mass0*Cud;
+ td.cloud().UTrans()[cellI] += np0*mass0*(U0 - U1);
// Update enthalpy transfer
- td.cloud().hTrans()[cellI] += np0*dhTrans;
-
- // Coefficient to be applied in carrier phase enthalpy coupling
- td.cloud().hCoeff()[cellI] += np0*htc*this->areaS();
+ td.cloud().hTrans()[cellI] += np0*mass0*(H0 - H1);
}
// Set new particle properties
@@ -165,19 +141,16 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
const scalar rho,
const scalar T,
const scalar cp,
- const scalar Sh,
- scalar& htc,
- scalar& dhTrans
+ const scalar Sh
)
{
if (!td.cloud().heatTransfer().active())
{
- htc = 0.0;
return T;
}
// Calc heat transfer coefficient
- htc = td.cloud().heatTransfer().h
+ scalar htc = td.cloud().heatTransfer().h
(
d,
U - this->Uc_,
@@ -192,11 +165,8 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
// Determine new particle temperature
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- // Particle area
- scalar Ap = this->areaS(d);
-
// Determine ap and bp coefficients
- scalar ap = Tc_ + Sh/(htc*Ap + ROOTVSMALL);
+ scalar ap = Tc_ + Sh/(htc*this->areaS(d) + ROOTVSMALL);
scalar bp = 6.0*htc/(rho*d*cp);
if (td.cloud().radiation())
{
@@ -219,10 +189,6 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
IntegrationScheme<scalar>::integrationResult Tres =
td.cloud().TIntegrator().integrate(T, dt, ap, bp);
- // Enthalpy transfer
- // - Using average particle temperature
- dhTrans = dt*Ap*htc*(Tres.average() - Tc_);
-
return Tres.value();
}
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
index 18a36dc024fbfbc38437d581856f804b528f5469..aa2020aa10b94990ab4b4a555f2c1b86317a2a42 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
@@ -222,9 +222,7 @@ protected:
const scalar rho, // density
const scalar T, // temperature
const scalar cp, // specific heat capacity
- const scalar Sh, // additional enthalpy sources
- scalar& htc, // heat transfer coeff
- scalar& dhTrans // enthalpy transfer to carrier phase
+ const scalar Sh // additional enthalpy sources
);