ENH: updated parcel energy coupling + style updates

This commit is contained in:
andy
2011-03-23 16:07:59 +00:00
parent d94985a94b
commit 11556949d1
7 changed files with 89 additions and 56 deletions

View File

@ -157,12 +157,15 @@ Foam::scalar Foam::COxidationDiffusionLimitedRate<CloudType>::calculate
dMassSRCarrier[O2GlobalId_] -= dmO2;
dMassSRCarrier[CO2GlobalId_] += dmCO2;
const scalar HC = thermo.solids().properties()[CsLocalId_].H(T);
const scalar HsC = thermo.solids().properties()[CsLocalId_].Hs(T);
const scalar HCO2 = thermo.carrier().H(CO2GlobalId_, T);
const scalar HO2 = thermo.carrier().H(O2GlobalId_, T);
// carrier enthalpy transfer handled by change in composition
// const scalar HsO2 = thermo.carrier().Hs(O2GlobalId_, T);
// dhsTrans -= dmO2*HsO2;
// Heat of reaction [J]
return dmC*HC + dmO2*HO2 - dmCO2*HCO2;
return dmC*HsC - dmCO2*HCO2;
}

View File

@ -160,12 +160,15 @@ Foam::scalar Foam::COxidationKineticDiffusionLimitedRate<CloudType>::calculate
dMassSRCarrier[O2GlobalId_] -= dmO2;
dMassSRCarrier[CO2GlobalId_] += dmCO2;
const scalar HC = thermo.solids().properties()[CsLocalId_].H(T);
const scalar HsC = thermo.solids().properties()[CsLocalId_].Hs(T);
const scalar HCO2 = thermo.carrier().H(CO2GlobalId_, T);
const scalar HO2 = thermo.carrier().H(O2GlobalId_, T);
// carrier enthalpy transfer handled by change in composition
// const scalar HsO2 = thermo.carrier().Hs(O2GlobalId_, T);
// dhsTrans -= dmO2*HsO2;
// Heat of reaction [J]
return dmC*HC + dmO2*HO2 - dmCO2*HCO2;
return dmC*HsC - dmCO2*HCO2;
}

View File

@ -220,12 +220,15 @@ Foam::scalar Foam::COxidationMurphyShaddix<CloudType>::calculate
// Add to particle mass transfer
dMassSolid[CsLocalId_] += dOmega*WC_;
const scalar HC = thermo.solids().properties()[CsLocalId_].H(T);
const scalar HsC = thermo.solids().properties()[CsLocalId_].Hs(T);
const scalar HCO2 = thermo.carrier().H(CO2GlobalId_, T);
const scalar HO2 = thermo.carrier().H(O2GlobalId_, T);
// carrier enthalpy transfer handled by change in composition
// const scalar HsO2 = thermo.carrier().Hs(O2GlobalId_, T);
// dhsTrans -= dmO2*HsO2;
// Heat of reaction
return dOmega*(WC_*HC + WO2_*HO2 - (WC_ + WO2_)*HCO2);
return dOmega*(WC_*HsC - (WC_ + WO2_)*HCO2);
}

View File

@ -165,6 +165,11 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
const label cellI
)
{
typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
const CompositionModel<reactingCloudType>& composition =
td.cloud().composition();
// Define local properties at beginning of timestep
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
const scalar np0 = this->nParticle_;
@ -178,9 +183,9 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
const scalar pc = this->pc_;
const scalarField& YMix = this->Y_;
const label idG = td.cloud().composition().idGas();
const label idL = td.cloud().composition().idLiquid();
const label idS = td.cloud().composition().idSolid();
const label idG = composition.idGas();
const label idL = composition.idLiquid();
const label idS = composition.idSolid();
// Calc surface values
@ -221,7 +226,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
scalar NCpW = 0.0;
// Surface concentrations of emitted species
scalarField Cs(td.cloud().composition().carrier().species().size(), 0.0);
scalarField Cs(composition.carrier().species().size(), 0.0);
// Calc mass and enthalpy transfer due to phase change
this->calcPhaseChange
@ -271,10 +276,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
Cs
);
// Correct surface values due to emitted species
this->correctSurfaceValues(td, cellI, Ts, Cs, rhos, mus, Prs, kappas);
Res = this->Re(U0, d0, rhos, mus);
// Surface reactions
// ~~~~~~~~~~~~~~~~~
@ -283,13 +284,9 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
scalarField dMassSRGas(YGas_.size(), 0.0);
scalarField dMassSRLiquid(YLiquid_.size(), 0.0);
scalarField dMassSRSolid(YSolid_.size(), 0.0);
scalarField dMassSRCarrier
(
td.cloud().composition().carrier().species().size(),
0.0
);
scalarField dMassSRCarrier(composition.carrier().species().size(), 0.0);
// Clac mass and enthalpy transfer due to surface reactions
// Calc mass and enthalpy transfer due to surface reactions
calcSurfaceReactions
(
td,
@ -313,6 +310,11 @@ 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
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -324,6 +326,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
updateMassFractions(mass0, dMassGas, dMassLiquid, dMassSolid);
// Heat transfer
// ~~~~~~~~~~~~~
@ -380,25 +384,33 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
// Transfer mass lost from particle to carrier mass source
forAll(YGas_, i)
{
label gid = td.cloud().composition().localToGlobalCarrierId(GAS, i);
label gid = composition.localToGlobalCarrierId(GAS, i);
td.cloud().rhoTrans(gid)[cellI] += np0*dMassGas[i];
// td.cloud().hsTrans()[cellI] +=
// np0*dMassGas[i]*composition.carrier().Hs(gid, T0);
}
forAll(YLiquid_, i)
{
label gid = td.cloud().composition().localToGlobalCarrierId(LIQ, i);
label gid = composition.localToGlobalCarrierId(LIQ, i);
td.cloud().rhoTrans(gid)[cellI] += np0*dMassLiquid[i];
// td.cloud().hsTrans()[cellI] +=
// np0*dMassLiquid[i]*composition.carrier().Hs(gid, T0);
}
/*
// No mapping between solid components and carrier phase
forAll(YSolid_, i)
{
label gid = td.cloud().composition().localToGlobalCarrierId(SLD, i);
label gid = composition.localToGlobalCarrierId(SLD, i);
td.cloud().rhoTrans(gid)[cellI] += np0*dMassSolid[i];
// td.cloud().hsTrans()[cellI] +=
// np0*dMassSolid[i]*composition.carrier().Hs(gid, T0);
}
*/
forAll(dMassSRCarrier, i)
{
td.cloud().rhoTrans(i)[cellI] += np0*dMassSRCarrier[i];
// td.cloud().hsTrans()[cellI] +=
// np0*dMassSRCarrier[i]*composition.carrier().Hs(i, T0);
}
// Update momentum transfer
@ -427,14 +439,12 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
// Absorb parcel into carrier phase
forAll(YGas_, i)
{
label gid =
td.cloud().composition().localToGlobalCarrierId(GAS, i);
label gid = composition.localToGlobalCarrierId(GAS, i);
td.cloud().rhoTrans(gid)[cellI] += np0*mass1*YMix[GAS]*YGas_[i];
}
forAll(YLiquid_, i)
{
label gid =
td.cloud().composition().localToGlobalCarrierId(LIQ, i);
label gid = composition.localToGlobalCarrierId(LIQ, i);
td.cloud().rhoTrans(gid)[cellI] +=
np0*mass1*YMix[LIQ]*YLiquid_[i];
}
@ -442,8 +452,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
// No mapping between solid components and carrier phase
forAll(YSolid_, i)
{
label gid =
td.cloud().composition().localToGlobalCarrierId(SLD, i);
label gid = composition.localToGlobalCarrierId(SLD, i);
td.cloud().rhoTrans(gid)[cellI] +=
np0*mass1*YMix[SLD]*YSolid_[i];
}
@ -508,6 +517,11 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
return;
}
typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
const CompositionModel<reactingCloudType>& composition =
td.cloud().composition();
// Total mass of volatiles evolved
td.cloud().devolatilisation().calculate
(
@ -535,10 +549,9 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
// Note: hardcoded gaseous diffusivities for now
// TODO: add to carrier thermo
const scalar beta = sqr(cbrt(15.0) + cbrt(15.0));
const label id =
td.cloud().composition().localToGlobalCarrierId(GAS, i);
const scalar Cp = td.cloud().composition().carrier().Cp(id, Ts);
const scalar W = td.cloud().composition().carrier().W(id);
const label id = composition.localToGlobalCarrierId(GAS, i);
const scalar Cp = composition.carrier().Cp(id, Ts);
const scalar W = composition.carrier().W(id);
const scalar Ni = dMassDV[i]/(this->areaS(d)*dt*W);
// Dab calc'd using API vapour mass diffusivity function

View File

@ -252,6 +252,11 @@ void Foam::ReactingParcel<ParcelType>::calc
const label cellI
)
{
typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
const CompositionModel<reactingCloudType>& composition =
td.cloud().composition();
// Define local properties at beginning of time step
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
const scalar np0 = this->nParticle_;
@ -301,7 +306,7 @@ void Foam::ReactingParcel<ParcelType>::calc
scalar NCpW = 0.0;
// Surface concentrations of emitted species
scalarField Cs(td.cloud().composition().carrier().species().size(), 0.0);
scalarField Cs(composition.carrier().species().size(), 0.0);
// Calc mass and enthalpy transfer due to phase change
calcPhaseChange
@ -389,8 +394,10 @@ void Foam::ReactingParcel<ParcelType>::calc
// Transfer mass lost from particle to carrier mass source
forAll(dMassPC, i)
{
label gid = td.cloud().composition().localToGlobalCarrierId(0, i);
label gid = composition.localToGlobalCarrierId(0, i);
td.cloud().rhoTrans(gid)[cellI] += np0*dMassPC[i];
// td.cloud().hsTrans()[cellI] +=
// np0*dMassPC[i]*composition.carrier().Hs(gid, T0);
}
// Update momentum transfer
@ -418,13 +425,12 @@ void Foam::ReactingParcel<ParcelType>::calc
// Absorb parcel into carrier phase
forAll(Y_, i)
{
label gid =
td.cloud().composition().localToGlobalCarrierId(0, i);
label gid = composition.localToGlobalCarrierId(0, i);
td.cloud().rhoTrans(gid)[cellI] += np0*mass1*Y_[i];
}
td.cloud().UTrans()[cellI] += np0*mass1*U1;
td.cloud().hsTrans()[cellI] +=
np0*mass1*td.cloud().composition().H(0, Y_, pc_, T1);
np0*mass1*composition.H(0, Y_, pc_, T1);
}
}
@ -434,7 +440,7 @@ void Foam::ReactingParcel<ParcelType>::calc
else
{
this->Cp_ = td.cloud().composition().Cp(0, Y_, pc_, T1);
this->Cp_ = composition.Cp(0, Y_, pc_, T1);
this->T_ = T1;
this->U_ = U1;
@ -484,6 +490,11 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
return;
}
typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
const CompositionModel<reactingCloudType>& composition =
td.cloud().composition();
// Calculate mass transfer due to phase change
td.cloud().phaseChange().calculate
(
@ -511,19 +522,18 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
forAll(YComponents, i)
{
const label idc =
td.cloud().composition().localToGlobalCarrierId(idPhase, i);
const label idl = td.cloud().composition().globalIds(idPhase)[i];
const label idc = composition.localToGlobalCarrierId(idPhase, i);
const label idl = composition.globalIds(idPhase)[i];
const scalar dh = td.cloud().phaseChange().dh(idc, idl, pc_, T);
Sh -= dMassPC[i]*dh/dt;
// Update particle surface thermo properties
const scalar Dab =
td.cloud().composition().liquids().properties()[idl].D(pc_, Ts, Wc);
composition.liquids().properties()[idl].D(pc_, Ts, Wc);
const scalar Cp = td.cloud().composition().carrier().Cp(idc, Ts);
const scalar W = td.cloud().composition().carrier().W(idc);
const scalar Cp = composition.carrier().Cp(idc, Ts);
const scalar W = composition.carrier().W(idc);
const scalar Ni = dMassPC[i]/(this->areaS(d)*dt*W);
// Molar flux of species coming from the particle (kmol/m^2/s)

View File

@ -208,8 +208,8 @@ void Foam::ThermoParcel<ParcelType>::calc
// Sum Ni*Cpi*Wi of emission species
scalar NCpW = 0.0;
// Calculate new particle velocity
scalar Cuh = 0.0;
// Calculate new particle temperature
scalar Sph = 0.0;
scalar T1 =
this->calcHeatTransfer
(
@ -226,7 +226,7 @@ void Foam::ThermoParcel<ParcelType>::calc
NCpW,
Sh,
dhsTrans,
Cuh
Sph
);
@ -267,7 +267,7 @@ void Foam::ThermoParcel<ParcelType>::calc
td.cloud().hsTrans()[cellI] += np0*dhsTrans;
// Update sensible enthalpy coefficient
td.cloud().hsCoeff()[cellI] += np0*Cuh*this->areaS();
td.cloud().hsCoeff()[cellI] += np0*Sph;
}
// Set new particle properties
@ -294,7 +294,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
const scalar NCpW,
const scalar Sh,
scalar& dhsTrans,
scalar& Cuh
scalar& Sph
)
{
if (!td.cloud().heatTransfer().active())
@ -317,6 +317,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
htc = max(htc, ROOTVSMALL);
const scalar As = this->areaS(d);
scalar ap = Tc_ + Sh/As/htc;
scalar bp = 6.0*(Sh/As + htc*(Tc_ - T));
if (td.cloud().radiation())
@ -337,9 +338,9 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
scalar Tnew = max(Tres.value(), td.cloud().constProps().TMin());
dhsTrans += dt*htc*As*(0.5*(T + Tnew) - Tc_);
Sph = dt*htc*As;
Cuh = dt*bp;
dhsTrans += Sph*(0.5*(T + Tnew) - Tc_);
return Tnew;
}

View File

@ -243,7 +243,7 @@ protected:
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
scalar& Cuh // linearised heat transfer coefficient
scalar& Sph // linearised heat transfer coefficient
);