OpenFOAM 中的喷雾子模型

从数学公式到 C++ 代码

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

dispersion

传统 RANS 框架下的公式

粒子速度方程:

ddtxp=upddtup=CDτpRep24(ugup)=CDτpRep24vrel\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}

其中的未知量CDC_D,vrel\mathbf{v}_{\mathbf{r e l}},τp\tau_{p} 计算如下:

CD=24Rep(1+16Rep2/3)Rep<1000CD=0.424Rep1000\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}

Rep=ρgvreldp/μgR e_{p}=\rho_{g}\left|\mathbf{v}_{\mathrm{rel}}\right| d_{p} / \mu_{g}

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

τp=ρldp2/18ρgvg\tau_{p}=\rho_{l} d_{p}^{2} / 18 \rho_{g} v_{g}

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

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

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

G(up,i)=1σ2πexp(up,i2σ)2G\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}

tturb=min(kε,Cpsk3/2ε1u~+uup)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)

其中Cps=0.16432C_{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,即u~+u\tilde{\mathbf{u}}+\mathbf{u}^{\prime} 赋值给了 td.Uc(),即当前 parcel 所在网格的气相速度,但是并没有真正的修改气相速度。
这里修改了传进去 UTurb_,但是在其它代码中并没有用到。
此函数最关键的用处就是更新了 td.Uc(),用于后续的 KinematicParcel::calcVelocity,见这篇本站博文

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

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

其中的 UTurb 就是上边公式中的u\mathbf{u}^{\prime}Uc 是公式中的u~\tilde{\mathbf{u}}U 是公式中的upu_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.

tturb=Cturb2ΔuSgSudt_{turb}=C_{turb} \frac{2 \Delta}{\left|\mathbf{u}_{S g S}-\mathbf{u}_{d}\right|}

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

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

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

G(usto,i)=12σsgs2πexp(usto,i22σsgs2)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 的返回值翻译成公式就是3mpμgCDRe4ρpD2\dfrac{3m_p\mu_g C_DRe}{4\rho_p D^2}
我们有:

Fp=mpdupdt=3mpμgCDRe4ρpD2(uup)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 返回值即(uup)(\mathbf{u}-\mathbf{u}_{p}) 前边的部分。
该返回值用于 KinematicParcel::calcVelocity 计算新的粒子速度,详细过程见这篇本站博文

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