From fd094124d9f021d59b4a2394495fc5d653361f13 Mon Sep 17 00:00:00 2001 From: andy Date: Wed, 20 Apr 2011 09:46:43 +0100 Subject: [PATCH] ENH: Updated cloud coupling terms --- .../KinematicParcel/KinematicParcelI.H | 2 +- .../ReactingMultiphaseParcel.C | 84 ++++++++++++------- .../ReactingMultiphaseParcel.H | 4 +- .../Templates/ReactingParcel/ReactingParcel.C | 81 ++++++++++++------ .../Templates/ThermoParcel/ThermoParcel.C | 3 + 5 files changed, 114 insertions(+), 60 deletions(-) diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H index 754677803c..06e2b1e462 100644 --- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H +++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H @@ -518,7 +518,7 @@ inline Foam::scalar Foam::KinematicParcel::Re const scalar muc ) const { - return rhoc*mag(U - Uc_)*d/muc; + return rhoc*mag(U - Uc_)*d/(muc + ROOTVSMALL); } diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C index 0d605da15d..5972bd73ae 100644 --- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C +++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C @@ -63,7 +63,7 @@ Foam::scalar Foam::ReactingMultiphaseParcel::CpEff template template -Foam::scalar Foam::ReactingMultiphaseParcel::HEff +Foam::scalar Foam::ReactingMultiphaseParcel::HsEff ( TrackData& td, const scalar p, @@ -74,9 +74,9 @@ Foam::scalar Foam::ReactingMultiphaseParcel::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::calc updateMassFractions(mass0, dMassGas, dMassLiquid, dMassSolid); - // Heat transfer // ~~~~~~~~~~~~~ @@ -383,25 +382,37 @@ void Foam::ReactingMultiphaseParcel::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::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; - // enthalpy transfer accounted for via change in mass fractions + td.cloud().hsTrans()[cellI] += dm*HsEff(td, pc, T1, idG, idL, idS); + + td.cloud().addToMassPhaseChange(dm); } } @@ -531,28 +544,35 @@ void Foam::ReactingMultiphaseParcel::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); + 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); - N += Ni; - NCpW += Ni*Cp*W; - Cs[id] += Ni*d/(2.0*Dab); - } + // 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); + } + } } diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H index cbb6bb8ead..7e9e16bbc7 100644 --- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H +++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H @@ -134,9 +134,9 @@ private: const label idS ) const; - //- Return the mixture effective enthalpy + //- Return the mixture effective sensible enthalpy template - scalar HEff + scalar HsEff ( TrackData& td, const scalar p, diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C index b5a2706732..c330254867 100644 --- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C +++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C @@ -213,9 +213,17 @@ void Foam::ReactingParcel::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::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::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::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]; - } - td.cloud().UTrans()[cellI] += np0*mass1*U1; + scalar hs = composition.carrier().Hs(gid, T1); - // enthalpy transfer accounted for via change in mass fractions + td.cloud().rhoTrans(gid)[cellI] += dmi; + td.cloud().hsTrans()[cellI] += dmi*hs; + } + td.cloud().UTrans()[cellI] += dm*U1; + + td.cloud().addToMassPhaseChange(dm); } } @@ -514,33 +534,44 @@ void Foam::ReactingParcel::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); - 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); + // 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_; - // 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; + forAll(dMassPC, i) + { + const label idc = composition.localToGlobalCarrierId(idPhase, i); + const label idl = composition.globalIds(idPhase)[i]; - // Concentrations of emission species - Cs[idc] += Ni*d/(2.0*Dab); + 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); + + 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); + } } } diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C index 6397e29437..dbc63e7c2e 100644 --- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C +++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C @@ -153,7 +153,10 @@ void Foam::ThermoParcel::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); }