从数学公式到C++代码全面解析喷雾传热传质过程原理
传热传质计算公式首先是能量传递,气相通过温差热传导传热给液相,液相接受能量,使自己升温并蒸发。蒸发出来的热量重新回到气相。这里有一个潜热的关系,即蒸发出来的那部分液相的能量+潜热=蒸发出来的气相所包含的能量。所以说,由喷雾引起的气相能量源项是:
1 V c ∑ p N p ( T p − T c ) π D p λ m N u m − 1 V c ∑ p m ˙ p N p h ( T p ) \frac{ 1 }{ V_{c} } \sum\nolimits_p N_p( T_p - T_{c} )\pi {D_p}\lambda_m Nu_m - \frac{ 1 }{ V_{c} }\sum\nolimits_p \dot m_p N_p h(T_p) V c 1 ∑ p N p ( T p − T c ) π D p λ m N u m − V c 1 ∑ p m ˙ p N p h ( T p )
其中,h ( T p ) h(T_p) h ( T p ) 是蒸发出来的气相处于液滴温度下的焓。如果气相求解的是显焓的输运方程,则这里是显焓。如果是绝对焓的输运方程,则这里是绝对焓。
另一方面,对于液相, 文献中液滴温度计算公式:
m ˙ p = d m p d t = − π D p D S h ρ m ln ( 1 + B M ) \dot m_p = \frac{dm_p}{dt}=-\pi D_p\mathcal{D}Sh \rho_m \ln(1+B_M) m ˙ p = d t d m p = − π D p D S h ρ m ln ( 1 + B M )
m p d h p d t = m ˙ p h v ( T p ) + π D κ N u ( T − T p ) f m_p\frac{dh_p}{dt}= \dot m_p h_v(T_p)+\pi D\kappa Nu(T-T_p)f m p d t d h p = m ˙ p h v ( T p ) + π D κ N u ( T − T p ) f
f = z e z − 1 , z = − c p , v m ˙ p π D κ N u f=\frac{z}{e^z-1}, z=-\frac{c_{p,v}\dot m_p}{\pi D\kappa Nu} f = e z − 1 z , z = − π D κ N u c p , v m ˙ p
d T p d t = T − T p τ h f − 1 c l , p h v ( T p ) τ e \frac{dT_p}{dt}=\frac{T-T_p}{\tau_h}f-\frac{1}{c_{l,p}}\frac{h_v(T_p)}{\tau_e} d t d T p = τ h T − T p f − c l , p 1 τ e h v ( T p )
τ e = m p π D D S h ρ v ln ( 1 + B ) \tau_e=\frac{m_p}{\pi D\mathcal{D}Sh\rho_v\ln(1+B)} τ e = π D D S h ρ v ln ( 1 + B ) m p
τ h = m p C l , p π D κ N u \tau_h=\frac{m_pC_{l,p}}{\pi D\kappa Nu} τ h = π D κ N u m p C l , p
代码解析首先从求解器代码中可以获知气相能量输运方程中的源项是 parcels.Sh(he)
。
Sh
函数定义于 ThermoCloudI.H:
template <class CloudType >inline Foam::tmp<Foam::fvScalarMatrix>Foam::ThermoCloud<CloudType>::Sh (volScalarField& hs) const { if (this ->solution ().coupled ()) { if (this ->solution ().semiImplicit ("h" )) { const volScalarField Cp (thermo_.thermo().Cp()) ; const volScalarField::Internal Vdt (this ->mesh ().V ()*this ->db ().time ().deltaT ()); return hsTrans ()/Vdt - fvm::SuSp (hsCoeff ()/(Cp*Vdt), hs) + hsCoeff ()/(Cp*Vdt)*hs; } else { tmp<fvScalarMatrix> tfvm (new fvScalarMatrix (hs, dimEnergy/dimTime)); fvScalarMatrix& fvm = tfvm.ref (); fvm.source () = -hsTrans ()/(this ->db ().time ().deltaT ()); return tfvm; } } return tmp<fvScalarMatrix>(new fvScalarMatrix (hs, dimEnergy/dimTime)); }
其中的 hsTrans()
定义于 ThermoCloudI.H
:return hsTrans_();
。 其中 hsTrans_
的类型是 autoPtr<volScalarField::Internal>
。 在 ReactingParcel
的 calc
函数中计算它:
hsTrans ()[cellI] += dm*hs 正数!hsTrans ()[cellI] += np0*dhsTrans 负数!
dm
是蒸发质量,hs
是蒸发出来的气相的显焓(按照液滴温度计算的),dhsTrans
在 calcHeatTransfer
中计算。calcHeatTransfer
定义于 ThermalParcel
,调用于 ReactingParcel
中的 calc
函数,返回新的粒子温度。
template <class ParcelType >template <class TrackCloudType >Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer ( TrackCloudType& cloud, trackingData& td, const scalar dt, const scalar Re, const scalar Pr, const scalar kappa, const scalar NCpW, const scalar Sh, scalar& dhsTrans, scalar& Sph ) { if (!cloud.heatTransfer ().active ()) { return T_; } const scalar d = this ->d (); const scalar rho = this ->rho (); const scalar As = this ->areaS (d); const scalar V = this ->volume (d); const scalar m = rho*V; scalar htc = cloud.heatTransfer ().htc (d, Re, Pr, kappa, NCpW); const scalar bcp = htc*As/(m*Cp_); const scalar acp = bcp*td.Tc (); scalar ancp = Sh; if (cloud.radiation ()) { const tetIndices tetIs = this ->currentTetIndices (); const scalar Gc = td.GInterp ().interpolate (this ->coordinates (), tetIs); const scalar sigma = physicoChemical::sigma.value (); const scalar epsilon = cloud.constProps ().epsilon0 (); ancp += As*epsilon*(Gc/4.0 - sigma*pow4 (T_)); } ancp /= m*Cp_; const scalar deltaT = cloud.TIntegrator ().delta (T_, dt, acp + ancp, bcp); const scalar deltaTncp = ancp*dt; const scalar deltaTcp = deltaT - deltaTncp; scalar Tnew = T_ + deltaT; Tnew = min (max (Tnew, cloud.constProps ().TMin ()), cloud.constProps ().TMax ()); dhsTrans -= m*Cp_*deltaTcp; Sph = dt*m*Cp_*bcp; return Tnew; }
其中用到的其它函数定义:
template <class Type >inline Type Foam::integrationScheme::explicitDelta( const Type& phi, const scalar dtEff, const Type& Alpha, const scalar Beta ) { return (Alpha - Beta*phi)*dtEff; } template <class Type >inline Type Foam::integrationScheme::delta( const Type& phi, const scalar dt, const Type& Alpha, const scalar Beta ) const { return explicitDelta (phi, dtEff (dt, Beta), Alpha, Beta); } Foam::scalar Foam::integrationSchemes::Euler::dtEff ( const scalar dt, const scalar Beta ) const { return dt/(1 + Beta*dt); } Foam::scalar Foam::integrationSchemes::analytical::dtEff ( const scalar dt, const scalar Beta ) const { return mag (Beta*dt) > small ? (1 - exp (- Beta*dt))/Beta : dt; }
calcHeatTransfer
函数中的数学原理:
A s = π d d As=\pi dd A s = π d d
h t c = N u κ D p htc=Nu\frac{\kappa}{D_p} h t c = N u D p κ
b c p = h t c ∗ A s m C p = N u κ π D p m C p bcp=\frac{htc*As}{mC_p}=\frac{Nu\kappa\pi D_p}{mC_p} b c p = m C p h t c ∗ A s = m C p N u κ π D p
a c p = b c p ∗ T c e l l = N u κ π D p T c e l l m C p acp=bcp*Tcell=\frac{Nu\kappa\pi D_p T_{cell}}{mC_p} a c p = b c p ∗ T c e l l = m C p N u κ π D p T c e l l
a n c p = S h m C p = m ˙ p h v m C p ancp=\frac{Sh}{mC_p}=\frac{\dot m_p h_v}{mC_p} a n c p = m C p S h = m C p m ˙ p h v
d e l t a ( p h i , d t , A l p h a , B e t a ) = ( A l p h a − B e t a ∗ p h i ) ∗ d t E f f delta(phi, dt, Alpha, Beta)=(Alpha-Beta*phi)*dtEff d e l t a ( p h i , d t , A l p h a , B e t a ) = ( A l p h a − B e t a ∗ p h i ) ∗ d t E f f
A l p h a = a c p + a n c p = N u κ π D T c m C p + S h m C p Alpha=acp+ancp=\frac{N_u\kappa\pi DTc}{mC_p}+\frac{Sh}{mC_p} A l p h a = a c p + a n c p = m C p N u κ π D T c + m C p S h
B e t a = b c p = N u κ π D m C p Beta=bcp=\frac{N_u\kappa\pi D}{mC_p} B e t a = b c p = m C p N u κ π D
p h i = T p phi=T_p p h i = T p
d e l t a T = ( A l p h a − B e t a ∗ T p ) ∗ d t E f f = ( N u κ π D p T c m C p + m ˙ p h v d t − N u κ D p T p m C p ) d t E f f deltaT=(Alpha-Beta*T_p)*dtEff=(\frac{Nu\kappa\pi D_pT_c}{mC_p}+\frac{\dot m_ph_v}{dt}-\frac{Nu\kappa D_pT_p}{mC_p})dtEff d e l t a T = ( A l p h a − B e t a ∗ T p ) ∗ d t E f f = ( m C p N u κ π D p T c + d t m ˙ p h v − m C p N u κ D p T p ) d t E f f
d t E f f = ( 1 − e x p ( − B e t a ∗ d t ) ) / B e t a dtEff=(1 - exp(- Beta*dt))/Beta d t E f f = ( 1 − e x p ( − B e t a ∗ d t ) ) / B e t a
d e l t a T c p = d e l t a T − d e l t a T n c p = d e l t a T − a n c p ∗ d t = ? deltaTcp=deltaT-deltaTncp=deltaT-ancp*dt=? d e l t a T c p = d e l t a T − d e l t a T n c p = d e l t a T − a n c p ∗ d t = ?
d h s T r a n s = − m ∗ C p ∗ d e l t a T c p = a n c p ∗ d t ∗ m ∗ C p − m ∗ C p ∗ d e l t a T dhsTrans = -m*C_p*deltaTcp=ancp*dt*m*C_p -m*Cp_*deltaT d h s T r a n s = − m ∗ C p ∗ d e l t a T c p = a n c p ∗ d t ∗ m ∗ C p − m ∗ C p ∗ d e l t a T
d h s T r a n s = m ˙ p h v d t − m C p d e l t a T dhsTrans =\dot m_p h_v dt-m C_p deltaT d h s T r a n s = m ˙ p h v d t − m C p d e l t a T
S p h = N u κ π D p d t Sph=Nu\kappa\pi D_pdt S p h = N u κ π D p d t
上述公式中的C p C_p C p 是液滴的比热容,D p D_p D p 是液滴直径,D \mathcal{D} D 是vapour diffusivity,h v h_v h v 是汽化潜热,T p T_p T p 是液滴温度。
calcPhaseChange
定义于 ReactingParcel
,调用于 ReactingParcel
中的 calc
函数。
template <class ParcelType >template <class TrackCloudType >void Foam::ReactingParcel<ParcelType>::calcPhaseChange( TrackCloudType& cloud, trackingData& td, const scalar dt, const scalar Re, const scalar Pr, const scalar Ts, const scalar nus, const scalar d, const scalar T, const scalar mass, const label idPhase, const scalar YPhase, const scalarField& Y, scalarField& dMassPC, scalar& Sh, scalar& N, scalar& NCpW, scalarField& Cs ) { typedef typename TrackCloudType::reactingCloudType reactingCloudType; const CompositionModel<reactingCloudType>& composition = cloud.composition (); PhaseChangeModel<reactingCloudType>& phaseChange = cloud.phaseChange (); if (!phaseChange.active () || (YPhase < small)) { return ; } scalarField X (composition.liquids().X(Y)) ; scalar Tvap = phaseChange.Tvap (X); if (T < Tvap) { return ; } const scalar TMax = phaseChange.TMax (td.pc (), X); const scalar Tdash = min (T, TMax); const scalar Tsdash = min (Ts, TMax); scalarField hmm (dMassPC) ; phaseChange.calculate ( dt, this ->cell (), Re, Pr, d, nus, Tdash, Tsdash, td.pc (), td.Tc (), X, dMassPC ); dMassPC = min (mass*YPhase*Y, dMassPC); const scalar dMassTot = sum (dMassPC); phaseChange.addToPhaseChangeMass (this ->nParticle_*dMassTot); forAll(dMassPC, i) { const label cid = composition.localToCarrierId (idPhase, i); const scalar dh = phaseChange.dh (cid, i, td.pc (), Tdash); Sh -= dMassPC[i]*dh/dt; } if (cloud.heatTransfer ().BirdCorrection ()) { const scalar Wc = td.rhoc ()*RR*td.Tc ()/td.pc (); forAll(dMassPC, i) { const label cid = composition.localToCarrierId (idPhase, i); const scalar Cp = composition.carrier ().Cp (cid, td.pc (), Tsdash); const scalar W = composition.carrier ().W (cid); const scalar Ni = dMassPC[i]/(this ->areaS (d)*dt*W); const scalar Dab = composition.liquids ().properties ()[i].D (td.pc (), Tsdash, Wc); N += Ni; NCpW += Ni*Cp*W; Cs[cid] += Ni*d/(2.0 *Dab); } } }
另外还有一个 Sh
的临时 scalar
,申明并初始化为 0 于 ReactingParcel
中的 calc
函数。 注释说明它是:Explicit enthalpy source for particle。 随后通过 calcPhaseChange
计算 Sh
:
forAll(dMassPC, i) { const label cid = composition.localToCarrierId (idPhase, i); const scalar dh = phaseChange.dh (cid, i, td.pc (), Tdash); Sh -= dMassPC[i]*dh/dt; }
这里的 dMassPC
是文献传质公式中的− m ˙ p d t -\dot m_p dt − m ˙ p d t ,因为蒸发模型得到的 dMassPC
是正数,但是液滴传质是损失质量,所以传质公式中的m ˙ p \dot m_p m ˙ p 是负数。 最后在 calc
中通过 calcHeatTransfer
使用 Sh
。
前面已经多次提到这个 ReactingParcel
中的 calc
函数,这里给出它的代码:
template <class ParcelType >template <class TrackCloudType >void Foam::ReactingParcel<ParcelType>::calc( TrackCloudType& cloud, trackingData& td, const scalar dt ) { typedef typename TrackCloudType::reactingCloudType reactingCloudType; const CompositionModel<reactingCloudType>& composition = cloud.composition (); const scalar np0 = this ->nParticle_; const scalar d0 = this ->d_; const vector& U0 = this ->U_; const scalar T0 = this ->T_; const scalar mass0 = this ->mass (); scalar Ts, rhos, mus, Prs, kappas; this ->calcSurfaceValues (cloud, td, T0, Ts, rhos, mus, Prs, kappas); scalar Res = this ->Re (rhos, U0, td.Uc (), d0, mus); vector Su = Zero; scalar Spu = 0.0 ; vector dUTrans = Zero; scalar Sh = 0.0 ; scalar Sph = 0.0 ; scalar dhsTrans = 0.0 ; scalarField dMassPC (Y_.size(), 0.0 ) ; scalar Ne = 0.0 ; scalar NCpW = 0.0 ; scalarField Cs (composition.carrier().species().size(), 0.0 ) ; calcPhaseChange ( cloud, td, dt, Res, Prs, Ts, mus/rhos, d0, T0, mass0, 0 , 1.0 , Y_, dMassPC, Sh, Ne, NCpW, Cs ); scalarField dMass (dMassPC) ; scalar mass1 = updateMassFraction (mass0, dMass, Y_); this ->Cp_ = composition.Cp (0 , Y_, td.pc (), T0); if (cloud.constProps ().constantVolume ()) { this ->rho_ = mass1/this ->volume (); } else { this ->d_ = cbrt (mass1/this ->rho_*6.0 /pi); } if (np0*mass1 < cloud.constProps ().minParcelMass ()) { td.keepParticle = false ; if (cloud.solution ().coupled ()) { scalar dm = np0*mass0; forAll(Y_, i) { scalar dmi = dm*Y_[i]; label gid = composition.localToCarrierId (0 , i); scalar hs = composition.carrier ().Hs (gid, td.pc (), T0); cloud.rhoTrans (gid)[this ->cell ()] += dmi; cloud.hsTrans ()[this ->cell ()] += dmi*hs; } cloud.UTrans ()[this ->cell ()] += dm*U0; cloud.phaseChange ().addToPhaseChangeMass (np0*mass1); } return ; } correctSurfaceValues (cloud, td, Ts, Cs, rhos, mus, Prs, kappas); Res = this ->Re (rhos, U0, td.Uc (), this ->d (), mus); this ->T_ = this ->calcHeatTransfer ( cloud, td, dt, Res, Prs, kappas, NCpW, Sh, dhsTrans, Sph ); this ->Cp_ = composition.Cp (0 , Y_, td.pc (), T0); this ->U_ = this ->calcVelocity (cloud, td, dt, Res, mus, mass1, Su, dUTrans, Spu); if (cloud.solution ().coupled ()) { forAll(dMass, i) { scalar dm = np0*dMass[i]; label gid = composition.localToCarrierId (0 , i); scalar hs = composition.carrier ().Hs (gid, td.pc (), T0); cloud.rhoTrans (gid)[this ->cell ()] += dm; cloud.UTrans ()[this ->cell ()] += dm*U0; cloud.hsTrans ()[this ->cell ()] += dm*hs; } cloud.UTrans ()[this ->cell ()] += np0*dUTrans; cloud.UCoeff ()[this ->cell ()] += np0*Spu; cloud.hsTrans ()[this ->cell ()] += np0*dhsTrans; cloud.hsCoeff ()[this ->cell ()] += np0*Sph; if (cloud.radiation ()) { const scalar ap = this ->areaP (); const scalar T4 = pow4 (T0); cloud.radAreaP ()[this ->cell ()] += dt*np0*ap; cloud.radT4 ()[this ->cell ()] += dt*np0*T4; cloud.radAreaPT4 ()[this ->cell ()] += dt*np0*ap*T4; } } }