从数学公式到 C++ 代码

本文下标 p 表示 parcel，g 表示 gas。

dispersion

传统 RANS 框架下的公式

粒子速度方程：

$\begin{aligned} \frac{d}{dt} \mathbf{x}_{\mathbf{p}} &=\mathbf{u}_{\mathbf{p}} \\ \frac{d}{d t} \mathbf{u}_{\mathbf{p}} &=\frac{C_{D}}{\tau_{p}} \frac{R e_{p}}{24}\left(\mathbf{u}_{\mathbf{g}}-\mathbf{u}_{\mathbf{p}}\right)=\frac{C_{D}}{\tau_{p}} \frac{R e_{p}}{24} \mathbf{v}_{\mathbf{r e l}} \end{aligned}$

其中的未知量$C_D$,$\mathbf{v}_{\mathbf{r e l}}$,$\tau_{p}$ 计算如下：

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

$R e_{p}=\rho_{g}\left|\mathbf{v}_{\mathrm{rel}}\right| d_{p} / \mu_{g}$

$\mathbf{v}_{\mathrm{rel}}=\tilde{\mathbf{u}}+\mathbf{u}^{\prime}-\mathbf{u}_{\mathrm{p}}$

$\tau_{p}=\rho_{l} d_{p}^{2} / 18 \rho_{g} v_{g}$

$\tilde{\mathbf{u}}$ 是气相速度，$\mathbf{u}^{\prime}$ is stochastic velocity vector accounting for turbulence dispersion of parcels via interaction with surrounding gases.
$\mathbf{u}^{\prime}$ 是高斯分布，均值等于 0，方差等于$\sigma^2=2 k / 3$。

若随机变量X服从一个数学期望为 μ、方差为 σ^2 的正态分布，记为 N(μ，σ^2)。

In this way each component of$\mathbf{u}^{\prime}$ is chosen randomly from the Gaussian distribution

$G\left(\mathbf{u}_{\mathbf{p}, \mathbf{i}}^{\prime}\right)=\frac{1}{\sigma \sqrt{2 \pi}} \exp \left(-\frac{\mathbf{u}_{\mathbf{p}, \mathbf{i}}^{\prime}}{2 \sigma}\right)^{2}$

$t_{t u r b}=\min \left(\frac{k}{\varepsilon}, C_{p s} \frac{k^{3 / 2}}{\varepsilon} \frac{1}{\left|\tilde{\mathbf{u}}+\mathbf{u}^{\prime}-\mathbf{u}_{p}\right|}\right)$

其中$C_{ps}=0.16432$

OF 中的 stochstic turbulence dispersion 模型代码

template<class CloudType>
Foam::vector Foam::StochasticDispersionRAS<CloudType>::update
(
    const scalar dt,
    const label celli,
    const vector& U,
    const vector& Uc,
    vector& UTurb,
    scalar& tTurb
)
{
    Random& rnd = this->owner().rndGen();

    const scalar cps = 0.16432;

    const scalar k = this->kPtr_->primitiveField()[celli];
    const scalar epsilon =
        this->epsilonPtr_->primitiveField()[celli] + rootVSmall;

    const scalar UrelMag = mag(U - Uc - UTurb);

    const scalar tTurbLoc =
        min(k/epsilon, cps*pow(k, 1.5)/epsilon/(UrelMag + small));


    // Parcel is perturbed by the turbulence
    if (dt < tTurbLoc)
    {
        tTurb += dt;

        if (tTurb > tTurbLoc)
        {
            tTurb = 0;

            const scalar sigma = sqrt(2*k/3.0);

            // Calculate a random direction dir distributed uniformly
            // in spherical coordinates

            const scalar theta = rnd.scalar01()*twoPi;
            const scalar u = 2*rnd.scalar01() - 1;

            const scalar a = sqrt(1 - sqr(u));
            const vector dir(a*cos(theta), a*sin(theta), u);

            UTurb = sigma*mag(rnd.scalarNormal())*dir;
        }
    }
    else
    {
        tTurb = great;
        UTurb = Zero;
    }

    return Uc + UTurb;
}

这个函数在哪里被调用？答：

template<class ParcelType>
template<class TrackCloudType>
void Foam::KinematicParcel<ParcelType>::calcDispersion
(
    TrackCloudType& cloud,
    trackingData& td,
    const scalar dt
)
{
    td.Uc() = cloud.dispersion().update
    (
        dt,
        this->cell(),
        U_,
        td.Uc(),
        UTurb_,
        tTurb_
    );
}

这个 calcDispersion 函数，把 dispersion 模型中计算出来的 Uc + UTurb，即$\tilde{\mathbf{u}}+\mathbf{u}^{\prime}$ 赋值给了 td.Uc()，即当前 parcel 所在位置的气相速度。
这里更新了传进去 UTurb_和tTurb_，这两个数都是只用于当前这个函数，更新它们是为了下一步使用。
此函数最关键的用处就是更新了 td.Uc()，用于后续的 KinematicParcel::calcVelocity，见这篇本站博文。

KinematicParcel::move 中调用了 p.calcDispersion(cloud, ttd, dt);。

另外还有 SprayParcel::calcBreakup， breakup 之后，child->calcDispersion(cloud, td, dt);。

其中的 UTurb 就是上边公式中的$\mathbf{u}^{\prime}$，Uc 是公式中的$\tilde{\mathbf{u}}$，U 是公式中的$u_p$。

---

针对 LES 的 dispersion 模型

Tsang C-W, Kuo C-W, Trujillo M, et al. Evaluation and validation of large-eddy simulation sub-grid spray dispersion models using high-fidelity volume-of-fluid simulation data and engine combustion network experimental data[J]. International Journal of Engine Research, 2019, 20(6):583–605.

$t_{turb}=C_{turb} \frac{2 \Delta}{\left|\mathbf{u}_{sgs}-\mathbf{u}_{d}\right|}$

$\mathbf{u}^{\prime}=\mathbf{u}_{sgs}+\mathbf{u}_{s t o}$

$\mathbf{u}_{s g s}=C_{s g s}(2 \tilde{\mathbf{u}}-3 \widetilde{\tilde{\mathbf{u}}}+\widetilde{\tilde{\mathbf{u}}})$

$\mathbf{u}_{sto}$ 是高斯分布，均值是0，方差是$\sigma_{s g s}^{2}=C_{sig} \frac{2}{3} k_{s g s}$，其中$C_{sig}$ 是模型参数。

$G\left(\mathbf{u}_{sto, i}\right)=\frac{1}{\sqrt{2 \sigma_{s g s}^{2} \pi}} \exp \left(-\frac{\mathbf{u}_{sto, i}^{2}}{2 \sigma_{s g s}^{2}}\right)$

---

drag

CdRe：

template<class CloudType>
Foam::scalar Foam::SphereDragForce<CloudType>::CdRe(const scalar Re)
{
    if (Re > 1000.0)
    {
        return 0.424*Re;
    }
    else
    {
        return 24.0*(1.0 + 1.0/6.0*pow(Re, 2.0/3.0));
    }
}

calcCoupled：

template<class CloudType>
Foam::forceSuSp Foam::SphereDragForce<CloudType>::calcCoupled
(
    const typename CloudType::parcelType& p,
    const typename CloudType::parcelType::trackingData& td,
    const scalar dt,
    const scalar mass,
    const scalar Re,
    const scalar muc
) const
{
    return forceSuSp(Zero, mass*0.75*muc*CdRe(Re)/(p.rho()*sqr(p.d())));
}

SphereDragForce::calcCoupled 的返回值翻译成公式就是$\dfrac{3m_p\mu_g C_DRe}{4\rho_p D^2}$。
我们有：

$F_p = m_p \dfrac{d \mathbf{u}_p}{d t} =\dfrac{3m_p\mu_gC_DRe}{4\rho_p D^2}(\mathbf{u}-\mathbf{u}_{p})$

SphereDragForce::calcCoupled 返回值即$(\mathbf{u}-\mathbf{u}_{p})$ 前边的部分。
该返回值用于 KinematicParcel::calcVelocity 计算新的粒子速度，详细过程见这篇本站博文。