Commit fd094124 authored by andy's avatar andy
Browse files

ENH: Updated cloud coupling terms

parent 76676b9d
......@@ -518,7 +518,7 @@ inline Foam::scalar Foam::KinematicParcel<ParcelType>::Re
const scalar muc
) const
{
return rhoc*mag(U - Uc_)*d/muc;
return rhoc*mag(U - Uc_)*d/(muc + ROOTVSMALL);
}
......
......@@ -63,7 +63,7 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::CpEff
template<class ParcelType>
template<class TrackData>
Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::HEff
Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::HsEff
(
TrackData& td,
const scalar p,
......@@ -74,9 +74,9 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::HEff
) const
{
return
this->Y_[GAS]*td.cloud().composition().H(idG, YGas_, p, T)
+ this->Y_[LIQ]*td.cloud().composition().H(idL, YLiquid_, p, T)
+ this->Y_[SLD]*td.cloud().composition().H(idS, YSolid_, p, T);
this->Y_[GAS]*td.cloud().composition().Hs(idG, YGas_, p, T)
+ this->Y_[LIQ]*td.cloud().composition().Hs(idL, YLiquid_, p, T)
+ this->Y_[SLD]*td.cloud().composition().Hs(idS, YSolid_, p, T);
}
......@@ -326,7 +326,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
updateMassFractions(mass0, dMassGas, dMassLiquid, dMassSolid);
// Heat transfer
// ~~~~~~~~~~~~~
......@@ -383,25 +382,37 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
// Transfer mass lost from particle to carrier mass source
forAll(YGas_, i)
{
scalar dm = np0*dMassGas[i];
label gid = composition.localToGlobalCarrierId(GAS, i);
td.cloud().rhoTrans(gid)[cellI] += np0*dMassGas[i];
scalar hs = composition.carrier().Hs(gid, 0.5*(T0 + T1));
td.cloud().rhoTrans(gid)[cellI] += dm;
td.cloud().hsTrans()[cellI] += dm*hs;
}
forAll(YLiquid_, i)
{
scalar dm = np0*dMassLiquid[i];
label gid = composition.localToGlobalCarrierId(LIQ, i);
td.cloud().rhoTrans(gid)[cellI] += np0*dMassLiquid[i];
scalar hs = composition.carrier().Hs(gid, 0.5*(T0 + T1));
td.cloud().rhoTrans(gid)[cellI] += dm;
td.cloud().hsTrans()[cellI] += dm*hs;
}
/*
// No mapping between solid components and carrier phase
forAll(YSolid_, i)
{
scalar dm = np0*dMassSolid[i];
label gid = composition.localToGlobalCarrierId(SLD, i);
td.cloud().rhoTrans(gid)[cellI] += np0*dMassSolid[i];
scalar hs = composition.carrier().Hs(gid, 0.5*(T0 + T1));
td.cloud().rhoTrans(gid)[cellI] += dm;
td.cloud().hsTrans()[cellI] += dm*hs;
}
*/
forAll(dMassSRCarrier, i)
{
td.cloud().rhoTrans(i)[cellI] += np0*dMassSRCarrier[i];
scalar dm = np0*dMassSRCarrier[i];
scalar hs = composition.carrier().Hs(i, 0.5*(T0 + T1));
td.cloud().rhoTrans(i)[cellI] += dm;
td.cloud().hsTrans()[cellI] += dm*hs;
}
// Update momentum transfer
......@@ -421,36 +432,38 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
// Remove the particle when mass falls below minimum threshold
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (mass1 < td.cloud().constProps().minParticleMass())
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] += np0*mass1*YMix[GAS]*YGas_[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] +=
np0*mass1*YMix[LIQ]*YLiquid_[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] +=
np0*mass1*YMix[SLD]*YSolid_[i];
td.cloud().rhoTrans(gid)[cellI] += dm*YMix[SLD]*YSolid_[i];
}
*/
td.cloud().UTrans()[cellI] += np0*mass1*U1;
td.cloud().UTrans()[cellI] += dm*U1;
td.cloud().hsTrans()[cellI] += dm*HsEff(td, pc, T1, idG, idL, idS);
// enthalpy transfer accounted for via change in mass fractions
td.cloud().addToMassPhaseChange(dm);
}
}
......@@ -531,28 +544,35 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
Sh -= dMassTot*td.cloud().constProps().LDevol()/dt;
// Molar average molecular weight of carrier mix
const scalar Wc = this->rhoc_*specie::RR*this->Tc_/this->pc_;
// Update molar emissions
forAll(dMassDV, i)
if (td.cloud().heatTransfer().BirdCorrection())
{
// Molar average molecular weight of carrier mix
const scalar Wc =
max(SMALL, this->rhoc_*specie::RR*this->Tc_/this->pc_);
// Note: hardcoded gaseous diffusivities for now
// TODO: add to carrier thermo
const scalar beta = sqr(cbrt(15.0) + cbrt(15.0));
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
const scalar Dab =
3.6059e-3*(pow(1.8*Ts, 1.75))*sqrt(1.0/W + 1.0/Wc)/(this->pc_*beta);
N += Ni;
NCpW += Ni*Cp*W;
Cs[id] += Ni*d/(2.0*Dab);
}
forAll(dMassDV, i)
{
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
const scalar Dab =
3.6059e-3*(pow(1.8*Ts, 1.75))
*sqrt(1.0/W + 1.0/Wc)
/(this->pc_*beta);
N += Ni;
NCpW += Ni*Cp*W;
Cs[id] += Ni*d/(2.0*Dab);
}
}
}
......
......@@ -134,9 +134,9 @@ private:
const label idS
) const;
//- Return the mixture effective enthalpy
//- Return the mixture effective sensible enthalpy
template<class TrackData>
scalar HEff
scalar HsEff
(
TrackData& td,
const scalar p,
......
......@@ -213,9 +213,17 @@ void Foam::ReactingParcel<ParcelType>::correctSurfaceValues
sumYiCbrtW += Ys[i]*cbrtW;
}
Cps = max(Cps, ROOTVSMALL);
rhos *= pc_/(specie::RR*T);
rhos = max(rhos, ROOTVSMALL);
mus /= sumYiSqrtW;
mus = max(mus, ROOTVSMALL);
kappas /= sumYiCbrtW;
kappas = max(kappas, ROOTVSMALL);
Prs = Cps*mus/kappas;
}
......@@ -335,7 +343,9 @@ void Foam::ReactingParcel<ParcelType>::calc
Res = this->Re(U0, d0, rhos, mus);
// Update particle component mass and mass fractions
scalar mass1 = updateMassFraction(mass0, dMassPC, Y_);
scalarField dMass(dMassPC);
scalar mass1 = updateMassFraction(mass0, dMass, Y_);
// Heat transfer
......@@ -390,11 +400,15 @@ void Foam::ReactingParcel<ParcelType>::calc
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (td.cloud().solution().coupled())
{
// Transfer mass lost from particle to carrier mass source
forAll(dMassPC, i)
// Transfer mass lost to carrier mass and enthalpy sources
forAll(dMass, i)
{
scalar dm = np0*dMass[i];
label gid = composition.localToGlobalCarrierId(0, i);
td.cloud().rhoTrans(gid)[cellI] += np0*dMassPC[i];
scalar hs = composition.carrier().Hs(gid, 0.5*(T0 + T1));
td.cloud().rhoTrans(gid)[cellI] += dm;
td.cloud().hsTrans()[cellI] += dm*hs;
}
// Update momentum transfer
......@@ -413,21 +427,27 @@ void Foam::ReactingParcel<ParcelType>::calc
// Remove the particle when mass falls below minimum threshold
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (mass1 < td.cloud().constProps().minParticleMass())
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);
td.cloud().rhoTrans(gid)[cellI] += np0*mass1*Y_[i];
scalar hs = composition.carrier().Hs(gid, T1);
td.cloud().rhoTrans(gid)[cellI] += dmi;
td.cloud().hsTrans()[cellI] += dmi*hs;
}
td.cloud().UTrans()[cellI] += np0*mass1*U1;
td.cloud().UTrans()[cellI] += dm*U1;
// enthalpy transfer accounted for via change in mass fractions
td.cloud().addToMassPhaseChange(dm);
}
}
......@@ -514,33 +534,44 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
// Add to cumulative phase change mass
td.cloud().addToMassPhaseChange(this->nParticle_*dMassTot);
// Average molecular weight of carrier mix - assumes perfect gas
const scalar Wc = this->rhoc_*specie::RR*this->Tc_/this->pc_;
forAll(YComponents, i)
forAll(dMassPC, 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 =
composition.liquids().properties()[idl].D(pc_, Ts, Wc);
// Update molar emissions
if (td.cloud().heatTransfer().BirdCorrection())
{
// Average molecular weight of carrier mix - assumes perfect gas
const scalar Wc = this->rhoc_*specie::RR*this->Tc_/this->pc_;
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)
N += Ni;
forAll(dMassPC, i)
{
const label idc = composition.localToGlobalCarrierId(idPhase, i);
const label idl = composition.globalIds(idPhase)[i];
// Sum of Ni*Cpi*Wi of emission species
NCpW += Ni*Cp*W;
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);
// Concentrations of emission species
Cs[idc] += Ni*d/(2.0*Dab);
const scalar Dab =
composition.liquids().properties()[idl].D(pc_, Ts, Wc);
// Molar flux of species coming from the particle (kmol/m^2/s)
N += Ni;
// Sum of Ni*Cpi*Wi of emission species
NCpW += Ni*Cp*W;
// Concentrations of emission species
Cs[idc] += Ni*d/(2.0*Dab);
}
}
}
......
......@@ -153,7 +153,10 @@ void Foam::ThermoParcel<ParcelType>::calcSurfaceValues
mus = td.muInterp().interpolate(this->position(), tetIs)/TRatio;
Pr = td.cloud().constProps().Pr();
Pr = max(ROOTVSMALL, Pr);
kappas = Cpc_*mus/Pr;
kappas = max(ROOTVSMALL, kappas);
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment