从数学公式到 C++ 代码分析喷雾粒子速度计算，以及喷雾引起的气相速度源项

本文下标 p 表示 parcel，g 表示 gas，无下标的也是表示环境气体。

Niklas Nordin 博士论文中的公式：

为了计算粒子的速度，构造速度方程：

$m_p \dfrac{d \mathbf{u}_p}{d t}=\mathbf{F_p}$

$\mathbf{F_p}=-\dfrac{\pi D^{2}}{8} \rho C_D\left|\mathbf{u}_{p}-\mathbf{u}\right|\left(\mathbf{u}_{p}-\mathbf{u}\right)$

其中：

$C_D=\left\{\begin{array}{ll}{\dfrac{24}{Re_{p}}\left(1+\dfrac{1}{6} R e_{p}^{2 / 3}\right)} & {Re_p<1000} \\ {0.424} & {Re_p>1000}\end{array}\right.$

$Re_{p}=\dfrac{\rho\left|\mathbf{u}_{p}-\mathbf{u}\right| D_p}{\mu_g}$

可以推导出：

$a_p=\dfrac{d \mathbf{u}_{p}}{d t}=-\dfrac{\mathbf{u}_{p}-\mathbf{u}}{\tau_{\mathbf{u}}}$

$\tau_{\mathrm{u}} =\dfrac{8 m_{p}}{\pi \rho C_D D_p^{2}\left|\mathbf{u}_{p}-\mathbf{u}\right|} =\dfrac{4}{3} \dfrac{\rho_{p} D_p}{\rho C_{p}\left|\mathbf{u}_{p}-\mathbf{u}\right|}$

其中$a_p$ 是粒子的加速度，$D_p$ 是粒子的直径。

可以变形为：

$\begin{array}{ll} \dfrac{d \mathbf{u}_{p}}{d t}\\ =-\dfrac{\mathbf{u}_p-\mathbf{u}}{\tau_u}\\ =-\dfrac{\mathbf{u}_p-\mathbf{u}}{\dfrac{4}{3} \dfrac{\rho_{p} D_p}{\rho C_{p}\left|\mathbf{u}_{p}-\mathbf{u}\right|}}\\ =-\dfrac{3(\mathbf{u}_p-\mathbf{u})\rho C_D \left|\mathbf{u}_{p}-\mathbf{u}\right|}{4 \rho_p D_p}\\ =-\dfrac{3(\mathbf{u}_p-\mathbf{u})\mu_g C_D}{4\rho_p D^2}\dfrac{\rho\left|\mathbf{u}_{p}-\mathbf{u}\right| D_p}{\mu_g}\\ =-\dfrac{3}{4}\dfrac{\mu_g C_D Re_p}{\rho_p D_p^2}(\mathbf{u}_{p}-\mathbf{u}) \end{array}$

还原到牛顿第二定律中：

$F_p = m_p \dfrac{d \mathbf{u}_p}{d t}$

$=m_p \dfrac{3}{4}\dfrac{\mu_g C_D Re_p}{\rho_p D_p^2}(\mathbf{u}-\mathbf{u}_{p})$

$=\dfrac{3m_p\mu_gC_DRe}{4\rho_p D_p^2}(\mathbf{u}-\mathbf{u}_{p})$

以下代码中的 mass 经常是指当前 parcel 中包含的一个 particle 的质量，需要乘以粒子数才等于当前 parcel 的总质量。
ReactingParcel::calc 中调用：

// 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;

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

KinematicParcel::calcVelocity

template<class ParcelType>
template<class TrackCloudType>
const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
(
    TrackCloudType& cloud,
    trackingData& td,
    const scalar dt,
    const scalar Re,
    const scalar mu,
    const scalar mass,
    const vector& Su,
    vector& dUTrans,
    scalar& Spu
) const
{
    const typename TrackCloudType::parcelType& p =
        static_cast<const typename TrackCloudType::parcelType&>(*this);
    typename TrackCloudType::parcelType::trackingData& ttd =
        static_cast<typename TrackCloudType::parcelType::trackingData&>(td);


    const typename TrackCloudType::forceType& forces = cloud.forces();


    // Momentum source due to particle forces
    const forceSuSp Fcp = forces.calcCoupled(p, ttd, dt, mass, Re, mu);
    const forceSuSp Fncp = forces.calcNonCoupled(p, ttd, dt, mass, Re, mu);
    const scalar massEff = forces.massEff(p, ttd, mass);
    

    // Shortcut splitting assuming no implicit non-coupled force ...
    // Calculate the integration coefficients
    const vector acp = (Fcp.Sp()*td.Uc() + Fcp.Su())/massEff;
    const vector ancp = (Fncp.Su() + Su)/massEff;
    const scalar bcp = Fcp.Sp()/massEff;


    // Integrate to find the new parcel velocity
    const vector deltaU = cloud.UIntegrator().delta(U_, dt, acp + ancp, bcp);
    const vector deltaUncp = ancp*dt;
    const vector deltaUcp = deltaU - deltaUncp;


    // Calculate the new velocity and the momentum transfer terms
    vector Unew = U_ + deltaU;


    dUTrans -= massEff*deltaUcp;


    Spu = dt*Fcp.Sp();


    // Apply correction to velocity and dUTrans for reduced-D cases
    const polyMesh& mesh = cloud.pMesh();
    meshTools::constrainDirection(mesh, mesh.solutionD(), Unew);
    meshTools::constrainDirection(mesh, mesh.solutionD(), dUTrans);


    return Unew;
}

forces.calcCoupled 和 forces.calcNonCoupled 都定义于 ParticleForce 中，但是返回值是 0，在其派生类中重新定义真正的函数，当我们选用 SphereDragForce 时，返回值是$\dfrac{3m_p\mu C_DRe}{4\rho D_p^2}$。

然后再回到 KinematicParcel::calcVelocity：

代码中的变量	数学公式
Fcp	$\dfrac{3m_p\mu C_D Re}{4\rho_p D_p^2}$
Fncp	$0$
massEff	$m_p$
acp	$\mathbf{u}_{c}$ Fcp$/m_p=\dfrac{3 \mathbf{u}_{c} \mu C_D Re}{4\rho_p D_p^2}$
ancp	0
bcp	Fcp$/m_p$
deltaU	$\dfrac{\text{Fcp}}{m_p}(\mathbf{u_c}-\mathbf{u_p})dt=\dfrac{3\mu C_D Re}{4\rho_p D_p^2}(\mathbf{u_c}-\mathbf{u_p})dt$
deltaUncp	$0$
deltaUcp	deltaU
dUTrans	$-m_d*$deltaU
Spu	dt*Fcp

其中$\mathbf{u}_{c}$ 是粒子当前所在位置的插值气相速度。
可以发现这里的 deltaU 与上文提到的文献中的公式是一样的。这就是粒子速度方程了。

接下来讨论粒子对于气相速度的影响。首先看一下速度输运方程：
首先看一下速度输运方程：

tmp<fvVectorMatrix> tUEqn
(
    fvm::ddt(rho, U) + fvm::div(phi, U)
  + MRF.DDt(rho, U)
  + turbulence->divDevRhoReff(U)
==
    rho()*g
  + parcels.SU(U)
  + fvOptions(rho, U)
);

if (pimple.momentumPredictor())
{
    solve(UEqn == -fvc::grad(p));

    fvOptions.correct(U);
    K = 0.5*magSqr(U);
}

KinematicCloud::SU

template<class CloudType>
inline Foam::tmp<Foam::fvVectorMatrix>
Foam::KinematicCloud<CloudType>::SU(volVectorField& U) const
{
    if (debug)
    {
        Info<< "UTrans min/max = " << min(UTrans()).value() << ", "
            << max(UTrans()).value() << nl
            << "UCoeff min/max = " << min(UCoeff()).value() << ", "
            << max(UCoeff()).value() << endl;
    }


    if (solution_.coupled())
    {
        if (solution_.semiImplicit("U"))
        {
            const volScalarField::Internal
                Vdt(mesh_.V()*this->db().time().deltaT());


            return UTrans()/Vdt - fvm::Sp(UCoeff()/Vdt, U) + UCoeff()/Vdt*U;
        }
        else
        {
            tmp<fvVectorMatrix> tfvm(new fvVectorMatrix(U, dimForce));
            fvVectorMatrix& fvm = tfvm.ref();


            fvm.source() = -UTrans()/(this->db().time().deltaT());


            return tfvm;
        }
    }


    return tmp<fvVectorMatrix>(new fvVectorMatrix(U, dimForce));
}

ReactingParcel::calc 中：

// 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;
// np0 是当前 parcel 中的粒子数

再看一下文献：

$\frac{\partial \overline{\rho} \tilde{u}_{i}}{\partial t}+\frac{\partial \overline{\rho} \tilde{u}_{i} \tilde{u}_{j}}{\partial x_{j}}=-\frac{\partial \overline{p}}{\partial x_{i}}+\frac{\partial \overline{\tau}_{i j}}{\partial x_{j}}-\frac{\partial \overline{\rho} \Gamma_{i j}}{\partial x_{j}}+\overline{\dot{S}}_i$

$\overline{\dot{S}}_{i}=\frac{1}{V_{f}} \sum_{m=1}^{n p} \left(n_{d, m} \dot{m}_{d, m} u_{d, i, m}-n_{d, m} F_{d, i, m}\right)$

可以发现公式第一部分是液滴蒸发变成气相带来的速度，第二部分是液滴和气相之间的相互作用力。这两部分分别对应代码中的 dm*U0 和 np0*dUTrans。 其中的 dUTrans 在 KinematicParcel::calcVelocity 中计算：$dUTrans=-m_pdeltaU$。然后文献中的$F=m_d\dfrac{deltaU}{dt}$。这是所有 parcel 的合力，即$F_i=\sum_{m=1}^{n p} F_{d,i,m}$，这里的$i$ 表示力的三个方向，$m$ 表示对所有 parcel 循环，np 是当前网格 parcel 的个数。$n_{d,m}$ 表示第 m 个 parcel 中 particle 的个数，$\cot m_{d,m}$ 表示第 m 个 parcel 的蒸发速率，$u_{d,m}$ 表示第 m 个 parcel 的速度。

事实上，代码中：

UTrans = dm*U0 + np0*dUTrans
//dm是蒸发出来的总质量（已经乘过np0了），np0是当前parcel所包含的particle个数
SU.source() = - UTrans/deltaT = - (dm*U0 + np0*dUTrans)/deltaT 
= - (dm*U0 - np0*m*deltaU)/deltaT
= - (dm/deltaT*U0 - np0*m*deltaU/deltaT)

SU.source()$=-(n_d\dot m u_d - F)$