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

传热传质计算公式

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

$\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(T_p)$ 是蒸发出来的气相处于液滴温度下的焓。如果气相求解的是显焓的输运方程，则这里是显焓。如果是绝对焓的输运方程，则这里是绝对焓。

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

$\dot m_p = \frac{dm_p}{dt}=-\pi D_p\mathcal{D}Sh \rho_m \ln(1+B_M)$

$m_p\frac{dh_p}{dt}= \dot m_p h_v(T_p)+\pi D\kappa Nu(T-T_p)f$

$f=\frac{z}{e^z-1}, z=-\frac{c_{p,v}\dot m_p}{\pi D\kappa Nu}$

$\frac{dT_p}{dt}=\frac{T-T_p}{\tau_h}f-\frac{1}{c_{l,p}}\frac{h_v(T_p)}{\tau_e}$

$\tau_e=\frac{m_p}{\pi D\mathcal{D}Sh\rho_v\ln(1+B)}$

$\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.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;

    // 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	$\pi D_p^2$
htc	$Nu\frac{\kappa}{D_p}$
bcp = htc*As/(m*Cp_)	$\frac{Nu\kappa\pi D_p}{m_p C_p}$
acp=bcp*Tc	$\frac{Nu\kappa\pi D_p T_c}{mC_p}$
Sh = -dMassPC[i]*dh/dt	$\dot m_p h_v$
ancp = Sh/(m*Cp_)	$\frac{\dot m_p h_v}{m_p C_p}$
delta(phi, dt, Alpha, Beta)=(Alpha-Beta*phi)*dtEff	-
phi	$T_p$
Alpha=acp+ancp	$\frac{Nu\kappa\pi D_p T_c}{m_p C_p} + \frac{\dot m_p h_v}{m_p C_p}$
Beta=bcp	$\frac{Nu\kappa\pi D_p}{m_p C_p}$
analytical:dtEff = (1 - exp(- Beta*dt))/Beta	$\dfrac{1-e^{-bcp*dt}}{bcp}$
Euler:dtEff = dt/(1 + Beta*dt)	$\dfrac{dt}{1+bcp*dt}$
deltaT	$dtEff \frac{ Nu\kappa\pi D_p (T_c-T_p) + \dot m_p h_v }{m_p C_p} $
Tnew = T_ + deltaT	-
deltaTncp = ancp*dt	$\frac{\dot m_p h_v dt}{m_p C_p}$
deltaTcp = deltaT - deltaTncp	$dtEff \frac{ Nu\kappa\pi D_p (T_c-T_p) + \dot m_p h_v }{m_p C_p} - \frac{\dot m_p h_v dt}{m_p C_p}$
dhsTrans = -m*Cp_*deltaTcp	$\dot m_p h_v dt - dtEff ( Nu\kappa\pi D_p T_c + \dot m_p h_v - Nu\kappa D_p T_p )$
Sph = dt*m*Cp_*bcp	$Nu\kappa\pi D_p dt$

上述公式中的$C_p$ 是液滴的比热容，$D_p$ 是液滴直径，$\mathcal{D}$是vapour diffusivity，$h_v$是汽化潜热，$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);

    // 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 是文献传质公式中的$-\dot m_p dt$，因为蒸发模型得到的 dMassPC 是正数，但是液滴传质是损失质量，所以传质公式中的$\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;
        }
    }
}