OpenFOAM 中的喷雾传热传质

从数学公式到C++代码全面解析喷雾传热传质过程原理

传热传质计算公式

首先是能量传递,气相通过温差热传导传热给液相,液相接受能量,使自己升温并蒸发。蒸发出来的热量重新回到气相。这里有一个潜热的关系,即蒸发出来的那部分液相的能量+潜热=蒸发出来的气相所包含的能量。所以说,由喷雾引起的气相能量源项是:

1VcpNp(TpTc)πDpλmNum1Vcpm˙pNph(Tp)\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)

其中,h(Tp)h(T_p) 是蒸发出来的气相处于液滴温度下的焓。如果气相求解的是显焓的输运方程,则这里是显焓。如果是绝对焓的输运方程,则这里是绝对焓。

另一方面,对于液相,
文献中液滴温度计算公式:

m˙p=dmpdt=πDpDShρmln(1+BM)\dot m_p = \frac{dm_p}{dt}=-\pi D_p\mathcal{D}Sh \rho_m \ln(1+B_M)

mpdhpdt=m˙phv(Tp)+πDκNu(TTp)fm_p\frac{dh_p}{dt}= \dot m_p h_v(T_p)+\pi D\kappa Nu(T-T_p)f

f=zez1,z=cp,vm˙pπDκNuf=\frac{z}{e^z-1}, z=-\frac{c_{p,v}\dot m_p}{\pi D\kappa Nu}

dTpdt=TTpτhf1cl,phv(Tp)τe\frac{dT_p}{dt}=\frac{T-T_p}{\tau_h}f-\frac{1}{c_{l,p}}\frac{h_v(T_p)}{\tau_e}

τe=mpπDDShρvln(1+B)\tau_e=\frac{m_p}{\pi D\mathcal{D}Sh\rho_v\ln(1+B)}

τh=mpCl,pπDκNu\tau_h=\frac{m_pC_{l,p}}{\pi D\kappa Nu}

代码解析

首先从求解器代码中可以获知气相能量输运方程中的源项是 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.Hreturn hsTrans_();
其中 hsTrans_ 的类型是 autoPtr<volScalarField::Internal>
ReactingParcelcalc 函数中计算它:

hsTrans()[cellI] += dm*hs 正数!    
hsTrans()[cellI] += np0*dhsTrans 负数!

dm 是蒸发质量,hs 是蒸发出来的气相的显焓(按照液滴温度计算的),dhsTranscalcHeatTransfer中计算。
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;

// Calc heat transfer coefficient
scalar htc = cloud.heatTransfer().htc(d, Re, Pr, kappa, NCpW);

// Calculate the integration coefficients
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_;

// Integrate to find the new parcel temperature
const scalar deltaT = cloud.TIntegrator().delta(T_, dt, acp + ancp, bcp);
const scalar deltaTncp = ancp*dt;
const scalar deltaTcp = deltaT - deltaTncp;

// Calculate the new temperature and the enthalpy transfer terms
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 函数中的数学原理:

As=πddAs=\pi dd

htc=NuκDphtc=Nu\frac{\kappa}{D_p}

bcp=htcAsmCp=NuκπDpmCpbcp=\frac{htc*As}{mC_p}=\frac{Nu\kappa\pi D_p}{mC_p}

acp=bcpTcell=NuκπDpTcellmCpacp=bcp*Tcell=\frac{Nu\kappa\pi D_p T_{cell}}{mC_p}

ancp=ShmCp=m˙phvmCpancp=\frac{Sh}{mC_p}=\frac{\dot m_p h_v}{mC_p}

delta(phi,dt,Alpha,Beta)=(AlphaBetaphi)dtEffdelta(phi, dt, Alpha, Beta)=(Alpha-Beta*phi)*dtEff

Alpha=acp+ancp=NuκπDTcmCp+ShmCpAlpha=acp+ancp=\frac{N_u\kappa\pi DTc}{mC_p}+\frac{Sh}{mC_p}

Beta=bcp=NuκπDmCpBeta=bcp=\frac{N_u\kappa\pi D}{mC_p}

phi=Tpphi=T_p

deltaT=(AlphaBetaTp)dtEff=(NuκπDpTcmCp+m˙phvdtNuκDpTpmCp)dtEffdeltaT=(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

dtEff=(1exp(Betadt))/BetadtEff=(1 - exp(- Beta*dt))/Beta

deltaTcp=deltaTdeltaTncp=deltaTancpdt=?deltaTcp=deltaT-deltaTncp=deltaT-ancp*dt=?

dhsTrans=mCpdeltaTcp=ancpdtmCpmCpdeltaTdhsTrans = -m*C_p*deltaTcp=ancp*dt*m*C_p -m*Cp_*deltaT

dhsTrans=m˙phvdtmCpdeltaTdhsTrans =\dot m_p h_v dt-m C_p deltaT

Sph=NuκπDpdtSph=Nu\kappa\pi D_pdt

上述公式中的CpC_p 是液滴的比热容,DpD_p 是液滴直径,D\mathcal{D}是vapour diffusivity,hvh_v是汽化潜热,TpT_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);

// Calculate mass transfer due to phase change
phaseChange.calculate
(
dt,
this->cell(),
Re,
Pr,
d,
nus,
Tdash,
Tsdash,
td.pc(),
td.Tc(),
X,
dMassPC
);
//此函数的作用就是计算 scalarField& dMassPC,size=液相组分个数,其它全是输入参数

// Limit phase change mass by availability of each specie
dMassPC = min(mass*YPhase*Y, dMassPC);

const scalar dMassTot = sum(dMassPC);

// Add to cumulative phase change mass
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;
}


// Update molar emissions
if (cloud.heatTransfer().BirdCorrection())
{
// Average molecular weight of carrier mix - assumes perfect gas
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);

// 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[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是蒸发传质,正数。
}

这里的 dMassPC 是文献传质公式中的m˙pdt-\dot m_p dt,因为蒸发模型得到的 dMassPC 是正数,但是液滴传质是损失质量,所以传质公式中的m˙p\dot 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();


// Define local properties at beginning of time step
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

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();


// Calc surface values
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);


// Sources
// ~~~~~~~

// Explicit momentum source for particle
vector Su = Zero;

// Linearised momentum source coefficient
scalar Spu = 0.0;

// Momentum transfer from the particle to the carrier phase
vector dUTrans = Zero;

// Explicit enthalpy source for particle
scalar Sh = 0.0;

// Linearised enthalpy source coefficient
scalar Sph = 0.0;

// Sensible enthalpy transfer from the particle to the carrier phase
scalar dhsTrans = 0.0;


// 1. Compute models that contribute to mass transfer - U, T held constant
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

// Phase change
// ~~~~~~~~~~~~

// Mass transfer due to phase change
scalarField dMassPC(Y_.size(), 0.0);

// Molar flux of species emitted from the particle (kmol/m^2/s)
scalar Ne = 0.0;

// Sum Ni*Cpi*Wi of emission species
scalar NCpW = 0.0;

// Surface concentrations of emitted species
scalarField Cs(composition.carrier().species().size(), 0.0);

// Calc mass and enthalpy transfer due to phase change
calcPhaseChange
(
cloud,
td,
dt,
Res,
Prs,
Ts,
mus/rhos,
d0,
T0,
mass0,
0,
1.0,
Y_,
dMassPC,
Sh,
Ne,
NCpW,
Cs
);


// 2. Update the parcel properties due to change in mass
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

scalarField dMass(dMassPC);
scalar mass1 = updateMassFraction(mass0, dMass, Y_);

this->Cp_ = composition.Cp(0, Y_, td.pc(), T0);

// Update particle density or diameter
if (cloud.constProps().constantVolume())
{
this->rho_ = mass1/this->volume();
}
else
{
this->d_ = cbrt(mass1/this->rho_*6.0/pi);
}

// Remove the particle when mass falls below minimum threshold
if (np0*mass1 < cloud.constProps().minParcelMass())
{
td.keepParticle = false;

if (cloud.solution().coupled())
{
scalar dm = np0*mass0;

// Absorb parcel into carrier phase
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;
}

// Correct surface values due to emitted species
correctSurfaceValues(cloud, td, Ts, Cs, rhos, mus, Prs, kappas);
Res = this->Re(rhos, U0, td.Uc(), this->d(), mus);


// 3. Compute heat- and momentum transfers
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

// Heat transfer
// ~~~~~~~~~~~~~

// Calculate new particle temperature
this->T_ =
this->calcHeatTransfer
(
cloud,
td,
dt,
Res,
Prs,
kappas,
NCpW,
Sh,
dhsTrans,
Sph
);

this->Cp_ = composition.Cp(0, Y_, td.pc(), T0);


// Motion
// ~~~~~~

// Calculate new particle velocity
this->U_ =
this->calcVelocity(cloud, td, dt, Res, mus, mass1, Su, dUTrans, Spu);


// 4. Accumulate carrier phase source terms
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

if (cloud.solution().coupled())
{
// Transfer mass lost to carrier mass, momentum and enthalpy sources
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;
}

// Update momentum transfer
cloud.UTrans()[this->cell()] += np0*dUTrans;
cloud.UCoeff()[this->cell()] += np0*Spu;

// Update sensible enthalpy transfer
cloud.hsTrans()[this->cell()] += np0*dhsTrans;
cloud.hsCoeff()[this->cell()] += np0*Sph;

// Update radiation fields
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;
}
}
}
文章作者: Yan Zhang
文章链接: https://openfoam.top/sprayHeatMass/
版权声明: 本博客所有文章除特别声明外,均采用 CC BY-NC-SA 4.0 许可协议。转载请注明来自 OpenFOAM 成长之路
微信打赏给博主更多动力吧~