zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

updates to the parcel coupling routines

Browse Source
This commit is contained in:
andy
2009-06-05 16:41:48 +01:00
parent 8dcb8422f2
commit c9f759aab8
8 changed files with 283 additions and 195 deletions
Download Patch File Download Diff File Expand all files Collapse all files

41
src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
View File

@ -103,6 +103,13 @@ void Foam::KinematicParcel<ParcelType>::calc
const scalar mass0 = mass();
// Explicit momentum source for particle
vector Su = vector::zero;
// Momentum transfer from the particle to the carrier phase
vector dUTrans = vector::zero;
// Motion
// ~~~~~~
@ -110,7 +117,7 @@ void Foam::KinematicParcel<ParcelType>::calc
vector Fx = vector::zero;
// Calculate new particle velocity
vector U1 = calcVelocity(td, dt, cellI, d0, U0, rho0, mass0, Fx);
vector U1 = calcVelocity(td, dt, cellI, d0, U0, rho0, mass0, Su, dUTrans);
// Accumulate carrier phase source terms
@ -118,7 +125,7 @@ void Foam::KinematicParcel<ParcelType>::calc
if (td.cloud().coupled())
{
// Update momentum transfer
td.cloud().UTrans()[cellI] += np0*mass0*(U0 - U1);
td.cloud().UTrans()[cellI] += np0*dUTrans;
}
@ -139,29 +146,41 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
const vector& U,
const scalar rho,
const scalar mass,
const vector& Fx
const vector& Su,
vector& dUTrans
) const
{
const polyMesh& mesh = this->cloud().pMesh();
// Return linearised term from drag model
scalar Cud = td.cloud().drag().Cu(U - Uc_, d, rhoc_, rho, muc_);
// Momentum transfer coefficient
const scalar ptc =
td.cloud().drag().ptc(U - Uc_, d, rhoc_, muc_) + ROOTVSMALL;
// Calculate particle forces
vector Ftot = td.cloud().forces().calc(cellI, dt, rhoc_, rho, Uc_, U);
// Momentum source due to particle forces
const vector PFCoupled =
mass*td.cloud().forces().calcCoupled(cellI, dt, rhoc_, rho, Uc_, U);
const vector PFNonCoupled =
mass*td.cloud().forces().calcNonCoupled(cellI, dt, rhoc_, rho, Uc_, U);
// New particle velocity
//~~~~~~~~~~~~~~~~~~~~~~
// Update velocity - treat as 3-D
const vector ap = Uc_ + (Ftot + Fx)/(Cud + VSMALL);
const scalar bp = Cud;
const scalar As = this->areaS(d);
const vector ap = Uc_ + (PFCoupled + PFNonCoupled + Su)/(ptc*As);
const scalar bp = 6.0*ptc/(rho*d);
vector Unew = td.cloud().UIntegrator().integrate(U, dt, ap, bp).value();
IntegrationScheme<vector>::integrationResult Ures =
td.cloud().UIntegrator().integrate(U, dt, ap, bp);
// Apply correction to velocity for reduced-D cases
vector Unew = Ures.value();
dUTrans += dt*(ptc*As*(Ures.average() - Uc_) - PFCoupled);
// Apply correction to velocity and dUTrans for reduced-D cases
meshTools::constrainDirection(mesh, mesh.solutionD(), Unew);
meshTools::constrainDirection(mesh, mesh.solutionD(), dUTrans);
return Unew;
}

3
src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
View File

@ -243,7 +243,8 @@ protected:
const vector& U, // velocity
const scalar rho, // density
const scalar mass, // mass
const vector& Fx // additional forces
const vector& Su, // explicit particle momentum source
vector& dUTrans // momentum transfer to carrier
) const;

241
src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
View File

@ -174,14 +174,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::cellValueSourceCorrection
const scalar cpc = td.cpInterp().psi()[cellI];
this->cpc_ = (massCell*cpc + addedMass*cpEff)/massCellNew;
const scalar fCarrier =
(1.0 - td.cloud().constProps().hRetentionCoeff())
/td.cloud().constProps().hRetentionCoeff();
const scalar dh =
td.cloud().hsTrans()[cellI] - fCarrier*td.cloud().hcTrans()[cellI];
this->Tc_ += dh/(this->cpc_*massCellNew);
this->Tc_ += td.cloud().hsTrans()[cellI]/(this->cpc_*massCellNew);
}
@ -211,34 +204,17 @@ 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;
// Explicit momentum source for particle
vector Su = vector::zero;
// Momentum transfer from the particle to the carrier phase
vector dUTrans = vector::zero;
// Phase change in liquid phase
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Explicit enthalpy source for particle
scalar Sh = 0.0;
// Mass transfer due to phase change
scalarField dMassPC(YLiquid_.size(), 0.0);
// Return enthalpy source and calc mass transfer due to phase change
scalar ShPC =
calcPhaseChange
(
td,
dt,
cellI,
d0,
T0,
U0,
mass0,
idL,
YMix[LIQ],
YLiquid_,
dMassPC
);
// Sensible enthalpy transfer from the particle to the carrier phase
scalar dhsTrans = 0.0;
// Devolatilisation
@ -247,21 +223,22 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
// Mass transfer due to devolatilisation
scalarField dMassDV(YGas_.size(), 0.0);
// Return enthalpy source and calc mass transfer due to devolatilisation
scalar ShDV =
calcDevolatilisation
(
td,
dt,
T0,
mass0,
this->mass0_,
idG,
YMix[GAS],
YGas_,
canCombust_,
dMassDV
);
// Calc mass and enthalpy transfer due to devolatilisation
calcDevolatilisation
(
td,
dt,
T0,
mass0,
this->mass0_,
idG,
YMix[GAS],
YGas_,
canCombust_,
dMassDV,
Sh,
dhsTrans
);
// Surface reactions
@ -273,30 +250,53 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
scalarField dMassSRSolid(YSolid_.size(), 0.0);
scalarField dMassSRCarrier(td.cloud().carrierSpecies().size(), 0.0);
// Return enthalpy source and calc mass transfer(s) due to surface reaction
scalar hReaction =
calcSurfaceReactions
(
td,
dt,
cellI,
d0,
T0,
mass0,
canCombust_,
dMassDV, // assuming d(mass) due to phase change is non-volatile
YMix,
YGas_,
YLiquid_,
YSolid_,
dMassSRGas,
dMassSRLiquid,
dMassSRSolid,
dMassSRCarrier
);
// Clac mass and enthalpy transfer due to surface reactions
calcSurfaceReactions
(
td,
dt,
cellI,
d0,
T0,
mass0,
canCombust_,
dMassDV, // assuming d(mass) due to phase change is non-volatile
YMix,
YGas_,
YLiquid_,
YSolid_,
dMassSRGas,
dMassSRLiquid,
dMassSRSolid,
dMassSRCarrier,
Sh,
dhsTrans
);
// Heat of reaction retained by particle
const scalar ShSR = td.constProps().hRetentionCoeff()*hReaction;
// Phase change in liquid phase
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Mass transfer due to phase change
scalarField dMassPC(YLiquid_.size(), 0.0);
// Calc mass and enthalpy transfer due to phase change
calcPhaseChange
(
td,
dt,
cellI,
d0,
T0,
U0,
mass0,
idL,
YMix[LIQ],
YLiquid_,
dMassPC,
Sh,
dhsTrans
);
// Update component mass fractions
@ -313,27 +313,18 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
// 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;
scalar T1 =
calcHeatTransfer(td, dt, cellI, d0, U0, rho0, T0, cp0, Sh, dhsTrans);
// Motion
// ~~~~~~
// No additional forces
vector Fx = vector::zero;
// Calculate new particle velocity
vector U1 = calcVelocity(td, dt, cellI, d0, U0, rho0, mass0, Fx);
vector U1 = calcVelocity(td, dt, cellI, d0, U0, rho0, mass0, Su, dUTrans);
dUTrans += 0.5*(mass0 - mass1)*(U0 + U1);
// Accumulate carrier phase source terms
@ -344,35 +335,48 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
// Transfer mass lost from particle to carrier mass source
forAll(YGas_, i)
{
label id = td.cloud().composition().localToGlobalCarrierId(GAS, i);
td.cloud().rhoTrans(id)[cellI] += np0*dMassGas[i];
label gid = td.cloud().composition().localToGlobalCarrierId(GAS, i);
td.cloud().rhoTrans(gid)[cellI] += np0*dMassGas[i];
td.cloud().hcTrans()[cellI] +=
np0
*dMassGas[i]
*td.cloud().composition().carrierSpecies()[gid].H(T0);
}
forAll(YLiquid_, i)
{
label id = td.cloud().composition().localToGlobalCarrierId(LIQ, i);
td.cloud().rhoTrans(id)[cellI] += np0*dMassLiquid[i];
label gid = td.cloud().composition().localToGlobalCarrierId(LIQ, i);
td.cloud().rhoTrans(gid)[cellI] += np0*dMassLiquid[i];
td.cloud().hcTrans()[cellI] +=
np0
*dMassLiquid[i]
*td.cloud().composition().carrierSpecies()[gid].H(T0);
}
/*
// No mapping between solid components and carrier phase
forAll(YSolid_, i)
{
label id = td.cloud().composition().localToGlobalCarrierId(SLD, i);
td.cloud().rhoTrans(id)[cellI] += np0*dMassSolid[i];
label gid = td.cloud().composition().localToGlobalCarrierId(SLD, i);
td.cloud().rhoTrans(gid)[cellI] += np0*dMassSolid[i];
td.cloud().hcTrans()[cellI] +=
np0
*dMassSolid[i]
*td.cloud().composition().carrierSpecies()[gid].H(T0);
}
*/
forAll(dMassSRCarrier, i)
{
td.cloud().rhoTrans(i)[cellI] += np0*dMassSRCarrier[i];
td.cloud().hcTrans()[cellI] +=
np0
*dMassSRCarrier[i]
*td.cloud().composition().carrierSpecies()[i].H(T0);
}
// Update momentum transfer
td.cloud().UTrans()[cellI] += np0*(mass0*U0 - mass1*U1);
td.cloud().UTrans()[cellI] += np0*dUTrans;
// Update sensible enthalpy transfer
td.cloud().hsTrans()[cellI] += np0*(mass0*H0 - mass1*H1);
// Update chemical enthalpy transfer
td.cloud().hcTrans()[cellI] -= np0*ShSR;
td.cloud().hsTrans()[cellI] += np0*dhsTrans;
}
@ -388,29 +392,30 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
// Absorb parcel into carrier phase
forAll(YGas_, i)
{
label id =
label gid =
td.cloud().composition().localToGlobalCarrierId(GAS, i);
td.cloud().rhoTrans(id)[cellI] += np0*mass1*YMix[GAS]*YGas_[i];
td.cloud().rhoTrans(gid)[cellI] += np0*mass1*YMix[GAS]*YGas_[i];
}
forAll(YLiquid_, i)
{
label id =
label gid =
td.cloud().composition().localToGlobalCarrierId(LIQ, i);
td.cloud().rhoTrans(id)[cellI] +=
td.cloud().rhoTrans(gid)[cellI] +=
np0*mass1*YMix[LIQ]*YLiquid_[i];
}
/*
// No mapping between solid components and carrier phase
forAll(YSolid_, i)
{
label id =
label gid =
td.cloud().composition().localToGlobalCarrierId(SLD, i);
td.cloud().rhoTrans(id)[cellI] +=
td.cloud().rhoTrans(gid)[cellI] +=
np0*mass1*YMix[SLD]*YSolid_[i];
}
*/
td.cloud().UTrans()[cellI] += np0*mass1*U1;
td.cloud().hsTrans()[cellI] += np0*mass1*H1;
td.cloud().hsTrans()[cellI] +=
np0*mass1*HEff(td, pc, T1, idG, idL, idS);
}
}
@ -420,9 +425,9 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
else
{
this->U_ = U1;
this->T_ = T1;
this->cp_ = cpEff(td, pc, T1, idG, idL, idS);
this->T_ = T1;
this->U_ = U1;
// Update particle density or diameter
if (td.constProps().constantVolume())
@ -439,7 +444,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
template<class ParcelType>
template<class TrackData>
Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
(
TrackData& td,
const scalar dt,
@ -450,7 +455,9 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
const scalar YVolatileTot,
const scalarField& YVolatile,
bool& canCombust,
scalarField& dMassDV
scalarField& dMassDV,
scalar& Sh,
scalar& dhsTrans
) const
{
// Check that model is active, and that the parcel temperature is
@ -462,7 +469,7 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
|| T < td.constProps().Tbp()
)
{
return 0.0;
return;
}
// Total mass of volatiles evolved
@ -482,13 +489,13 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
td.cloud().addToMassDevolatilisation(this->nParticle_*dMassTot);
return -dMassTot*td.constProps().LDevol();
Sh -= dMassTot*td.constProps().LDevol()/dt;
}
template<class ParcelType>
template<class TrackData>
Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
void Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
(
TrackData& td,
const scalar dt,
@ -505,13 +512,15 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
scalarField& dMassSRGas,
scalarField& dMassSRLiquid,
scalarField& dMassSRSolid,
scalarField& dMassSRCarrier
scalarField& dMassSRCarrier,
scalar& Sh,
scalar& dhsTrans
) const
{
// Check that model is active
if (!td.cloud().surfaceReaction().active() || !canCombust)
{
return 0.0;
return;
}
// Update surface reactions
@ -542,7 +551,13 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::calcSurfaceReactions
*(sum(dMassSRGas) + sum(dMassSRLiquid) + sum(dMassSRSolid))
);
return hReaction;
const scalar xsi = min(T/5000.0, 1.0);
const scalar hRetentionCoeffMod =
(1.0 - xsi*xsi)*td.constProps().hRetentionCoeff();
Sh += hRetentionCoeffMod *hReaction/dt;
dhsTrans += (1.0 - hRetentionCoeffMod)*hReaction;
}

12
src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H
View File

@ -223,7 +223,7 @@ protected:
//- Calculate Devolatilisation
template<class TrackData>
scalar calcDevolatilisation
void calcDevolatilisation
(
TrackData& td,
const scalar dt, // timestep
@ -234,12 +234,14 @@ protected:
const scalar YVolatileTot, // total volatile mass fraction
const scalarField& YVolatile, // volatile component mass fractions
bool& canCombust, // 'can combust' flag
scalarField& dMassDV // mass transfer - local to particle
scalarField& dMassDV, // mass transfer - local to particle
scalar& Sh, // explicit particle enthalpy source
scalar& dhsTrans // sensible enthalpy transfer to carrier
) const;
//- Calculate surface reactions
template<class TrackData>
scalar calcSurfaceReactions
void calcSurfaceReactions
(
TrackData& td,
const scalar dt, // timestep
@ -256,7 +258,9 @@ protected:
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
scalarField& dMassSRCarrier, // carrier phase mass transfer
scalar& Sh, // explicit particle enthalpy source
scalar& dhsTrans // sensible enthalpy transfer to carrier
) const;

95
src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
View File

@ -138,10 +138,17 @@ 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;
// Explicit momentum source for particle
vector Su = vector::zero;
// Momentum transfer from the particle to the carrier phase
vector dUTrans = vector::zero;
// Explicit enthalpy source for particle
scalar Sh = 0.0;
// Sensible enthalpy transfer from the particle to the carrier phase
scalar dhsTrans = 0.0;
// Phase change
@ -150,9 +157,23 @@ void Foam::ReactingParcel<ParcelType>::calc
// Mass transfer due to phase change
scalarField dMassPC(Y_.size(), 0.0);
// Return enthalpy source and calc mass transfer due to phase change
scalar ShPC =
calcPhaseChange(td, dt, cellI, d0, T0, U0, mass0, 0, 1.0, Y_, dMassPC);
// Calc mass and enthalpy transfer due to phase change
calcPhaseChange
(
td,
dt,
cellI,
d0,
T0,
U0,
mass0,
0,
1.0,
Y_,
dMassPC,
Sh,
dhsTrans
);
// Update particle component mass and mass fractions
scalar mass1 = updateMassFraction(mass0, dMassPC, Y_);
@ -162,12 +183,8 @@ void Foam::ReactingParcel<ParcelType>::calc
// ~~~~~~~~~~~~~
// Calculate new particle temperature
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;
scalar T1 =
calcHeatTransfer(td, dt, cellI, d0, U0, rho0, T0, cp0, Sh, dhsTrans);
// Motion
@ -177,8 +194,9 @@ void Foam::ReactingParcel<ParcelType>::calc
vector Fx = vector::zero;
// Calculate new particle velocity
vector U1 = calcVelocity(td, dt, cellI, d0, U0, rho0, mass0, Fx);
vector U1 = calcVelocity(td, dt, cellI, d0, U0, rho0, mass0, Su, dUTrans);
dUTrans += 0.5*(mass0 - mass1)*(U0 + U1);
// Accumulate carrier phase source terms
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -187,15 +205,19 @@ void Foam::ReactingParcel<ParcelType>::calc
// Transfer mass lost from particle to carrier mass source
forAll(dMassPC, i)
{
label id = td.cloud().composition().localToGlobalCarrierId(0, i);
td.cloud().rhoTrans(id)[cellI] += np0*dMassPC[i];
label gid = td.cloud().composition().localToGlobalCarrierId(0, i);
td.cloud().rhoTrans(gid)[cellI] += np0*dMassPC[i];
td.cloud().hcTrans()[cellI] +=
np0
*dMassPC[i]
*td.cloud().composition().carrierSpecies()[gid].H(T0);
}
// Update momentum transfer
td.cloud().UTrans()[cellI] += np0*(mass0*U0 - mass1*U1);
td.cloud().UTrans()[cellI] += np0*dUTrans;
// Update sensible enthalpy transfer
td.cloud().hsTrans()[cellI] += np0*(mass0*H0 - mass1*H1);
td.cloud().hsTrans()[cellI] += np0*dhsTrans;
}
@ -210,12 +232,13 @@ void Foam::ReactingParcel<ParcelType>::calc
// Absorb parcel into carrier phase
forAll(Y_, i)
{
label id =
label gid =
td.cloud().composition().localToGlobalCarrierId(0, i);
td.cloud().rhoTrans(id)[cellI] += np0*mass1*Y_[i];
td.cloud().rhoTrans(gid)[cellI] += np0*mass1*Y_[i];
}
td.cloud().UTrans()[cellI] += np0*mass1*U1;
td.cloud().hsTrans()[cellI] += np0*mass1*H1;
td.cloud().hcTrans()[cellI] +=
np0*mass1*td.cloud().composition().H(0, Y_, pc_, T1);
}
}
@ -225,9 +248,9 @@ void Foam::ReactingParcel<ParcelType>::calc
else
{
this->U_ = U1;
this->T_ = T1;
this->cp_ = td.cloud().composition().cp(0, Y_, pc_, T1);
this->T_ = T1;
this->U_ = U1;
// Update particle density or diameter
if (td.constProps().constantVolume())
@ -244,7 +267,7 @@ void Foam::ReactingParcel<ParcelType>::calc
template<class ParcelType>
template<class TrackData>
Foam::scalar Foam::ReactingParcel<ParcelType>::calcPhaseChange
void Foam::ReactingParcel<ParcelType>::calcPhaseChange
(
TrackData& td,
const scalar dt,
@ -256,7 +279,9 @@ Foam::scalar Foam::ReactingParcel<ParcelType>::calcPhaseChange
const label idPhase,
const scalar YPhase,
const scalarField& YComponents,
scalarField& dMassPC
scalarField& dMassPC,
scalar& Sh,
scalar& dhsTrans
)
{
if
@ -266,7 +291,7 @@ Foam::scalar Foam::ReactingParcel<ParcelType>::calcPhaseChange
|| YPhase < SMALL
)
{
return 0.0;
return;
}
// Calculate mass transfer due to phase change
@ -275,7 +300,7 @@ Foam::scalar Foam::ReactingParcel<ParcelType>::calcPhaseChange
dt,
cellI,
d,
min(T, td.constProps().Tbp()), // Limit to boiling temperature
T,
pc_,
this->Tc_,
this->muc_/(this->rhoc_ + ROOTVSMALL),
@ -291,10 +316,20 @@ Foam::scalar Foam::ReactingParcel<ParcelType>::calcPhaseChange
// Add to cumulative phase change mass
td.cloud().addToMassPhaseChange(this->nParticle_*dMassTot);
// Effective latent heat of vaporisation
scalar LEff = td.cloud().composition().L(idPhase, YComponents, pc_, T);
// Enthalphy transfer to carrier phase
forAll(YComponents, i)
{
label gid;
return -dMassTot*LEff;
gid = td.cloud().composition().localToGlobalCarrierId(idPhase, i);
const scalar hv = td.cloud().composition().carrierSpecies()[gid].H(T);
gid = td.cloud().composition().globalIds(idPhase)[i];
const scalar hl =
td.cloud().composition().liquids().properties()[gid].h(pc_, T);
Sh += dMassPC[i]*(hl - hv)/dt;
}
}

6
src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H
View File

@ -190,7 +190,7 @@ protected:
//- Calculate Phase change
template<class TrackData>
scalar calcPhaseChange
void calcPhaseChange
(
TrackData& td,
const scalar dt, // timestep
@ -202,7 +202,9 @@ protected:
const label idPhase, // id of phase involved in phase change
const scalar YPhase, // total mass fraction
const scalarField& YComponents, // component mass fractions
scalarField& dMassPC // mass transfer - local to particle
scalarField& dMassPC, // mass transfer - local to particle
scalar& Sh, // explicit particle enthalpy source
scalar& dhsTrans // sensible enthalpy transfer to carrier
);
//- Update mass fraction

77
src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
View File

@ -96,31 +96,31 @@ void Foam::ThermoParcel<ParcelType>::calc
const scalar cp0 = this->cp_;
const scalar mass0 = this->mass();
// Initial enthalpy state
scalar H0 = cp0*T0;
// Explicit momentum source for particle
vector Su = vector::zero;
// Momentum transfer from the particle to the carrier phase
vector dUTrans = vector::zero;
// Explicit enthalpy source for particle
scalar Sh = 0.0;
// Sensible enthalpy transfer from the particle to the carrier phase
scalar dhsTrans = 0.0;
// Heat transfer
// ~~~~~~~~~~~~~
// No additional enthalpy sources
scalar Sh = 0.0;
// Calculate new particle velocity
scalar T1 = calcHeatTransfer(td, dt, cellI, d0, U0, rho0, T0, cp0, Sh);
// Calculate new enthalpy state
scalar H1 = cp0*T1;
scalar T1 =
calcHeatTransfer(td, dt, cellI, d0, U0, rho0, T0, cp0, Sh, dhsTrans);
// Motion
// ~~~~~~
// No additional forces
vector Fx = vector::zero;
// Calculate new particle velocity
vector U1 = calcVelocity(td, dt, cellI, d0, U0, rho0, mass0, Fx);
vector U1 = calcVelocity(td, dt, cellI, d0, U0, rho0, mass0, Su, dUTrans);
// Accumulate carrier phase source terms
@ -128,10 +128,10 @@ void Foam::ThermoParcel<ParcelType>::calc
if (td.cloud().coupled())
{
// Update momentum transfer
td.cloud().UTrans()[cellI] += np0*mass0*(U0 - U1);
td.cloud().UTrans()[cellI] += np0*dUTrans;
// Update sensible enthalpy transfer
td.cloud().hsTrans()[cellI] += np0*mass0*(H0 - H1);
td.cloud().hsTrans()[cellI] += np0*dhsTrans;
}
// Set new particle properties
@ -153,7 +153,8 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
const scalar rho,
const scalar T,
const scalar cp,
const scalar Sh
const scalar Sh,
scalar& dhsTrans
)
{
if (!td.cloud().heatTransfer().active())
@ -173,34 +174,44 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
this->muc_
);
const scalar As = this->areaS(d);
// Determine new particle temperature
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Determine ap and bp coefficients
scalar ap = Tc_ + Sh/(htc*this->areaS(d) + ROOTVSMALL);
scalar bp = 6.0*htc/(rho*d*cp);
if (td.cloud().radiation())
if (mag(htc) < ROOTVSMALL && !td.cloud().radiation())
{
// Carrier phase incident radiation field
// - The G field is not interpolated to the parcel position
// Instead, the cell centre value is applied directly
const scalarField& G = td.cloud().mesh().objectRegistry
::lookupObject<volScalarField>("G");
return T + dt*Sh/(this->mass()*cp);
}
// Helper variables
scalar ap;
scalar bp;
if(td.cloud().radiation())
{
const scalarField& G =
td.cloud().mesh().objectRegistry::lookupObject<volScalarField>("G");
const scalar sigma = radiation::sigmaSB.value();
const scalar epsilon = td.constProps().epsilon0();
const scalar D = epsilon*sigma*pow3(T)/(htc + ROOTVSMALL) + 1.0;
ap += 0.25*epsilon*G[cellI]/(htc + ROOTVSMALL);
ap /= D;
bp *= D;
ap =
(Sh/As + htc*Tc_ + epsilon*G[cellI]/4.0)
/(htc + epsilon*sigma*pow3(T));
bp =
6.0
*(Sh/As + htc*(Tc_ - T) + epsilon*(G[cellI]/4.0 - sigma*pow4(T)))
/(rho*d*cp*(ap - T));
}
else
{
ap = Tc_ + Sh/As/htc;
bp = 6.0*(Sh/As + htc*(Tc_ - T))/(rho*d*cp*(ap - T));
}
// Integrate to find the new parcel temperature
IntegrationScheme<scalar>::integrationResult Tres =
td.cloud().TIntegrator().integrate(T, dt, ap, bp);
dhsTrans += dt*htc*As*(Tres.average() - Tc_);
return Tres.value();
}

3
src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
View File

@ -222,7 +222,8 @@ protected:
const scalar rho, // density
const scalar T, // temperature
const scalar cp, // specific heat capacity
const scalar Sh // additional enthalpy sources
const scalar Sh, // explicit particle enthalpy source
scalar& dhsTrans // sensible enthalpy transfer to carrier
);