OpenFOAM 中的喷雾速度计算

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

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

Niklas Nordin 博士论文中的公式:

为了计算粒子的速度,构造速度方程:

mpdupdt=Fpm_p \dfrac{d \mathbf{u}_p}{d t}=\mathbf{F_p}

Fp=πD28ρCDupu(upu)\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)

其中:

CD={24Rep(1+16Rep2/3)Rep<10000.424Rep>1000C_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.

Rep=ρupudpμgRe_{p}=\dfrac{\rho\left|\mathbf{u}_{p}-\mathbf{u}\right| d_p}{\mu_g}

可以推导出:

ap=dupdt=upuτua_p=\dfrac{d \mathbf{u}_{p}}{d t}=-\dfrac{\mathbf{u}_{p}-\mathbf{u}}{\tau_{\mathbf{u}}}

τu=8mpπρCDdp2upu=43ρpdpρCpupu\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|}

其中apa_p 是粒子的加速度,dpd_p 是粒子的直径。

可以变形为:

dupdt=upuτu=upu43ρpdpρCpupu=3(upu)ρCDupu4ρpdp=3(upu)μgCD4ρpD2ρupudpμg=34μgCDRepρpdp2(upu)\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}

还原到牛顿第二定律中:

Fp=mpdupdtF_p = m_p \dfrac{d \mathbf{u}_p}{d t}

=mp34μgCDRepρpdp2(uup)=m_p \dfrac{3}{4}\dfrac{\mu_g C_D Re_p}{\rho_p d_p^2}(\mathbf{u}-\mathbf{u}_{p})

=3mpμgCDRe4ρpdp2(uup)=\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.calcCoupledforces.calcNonCoupled 都定义于 ParticleForce 中,但是返回值是 0,在其派生类中重新定义真正的函数,当我们选用 SphereDragForce 时,返回值是3mpμCDRe4ρdp2\dfrac{3m_p\mu C_DRe}{4\rho d_p^2}

然后再回到 KinematicParcel::calcVelocity

代码中的变量数学公式
Fcp3mpμCDRe4ρpdp2\dfrac{3m_p\mu C_DRe}{4\rho_p d_p^2}
Fncp00
massEffmpm_p
acpuc\mathbf{u}_{c} Fcp/mp=3ucμCDRe4ρpdp2/m_p=\dfrac{3 \mathbf{u}_{c} \mu C_D Re}{4\rho_p d_p^2}
ancp0
bcpFcp/mp/m_p
deltaUFcpmp(ucup)dt=3μCDRe4ρpdp2(ucup)dt\dfrac{\text{Fcp}}{m_p}(\mathbf{u_c}-\mathbf{u_p})dt=\dfrac{3\mu C_DRe}{4\rho_p d_p^2}(\mathbf{u_c}-\mathbf{u_p})dt
deltaUncp00
deltaUcpdeltaU
dUTransmd-m_d*deltaU
Spudt*Fcp

其中uc\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 中的粒子数

再看一下文献:

ρu~it+ρu~iu~jxj=pxi+τijxjρΓijxj+S˙i\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

S˙i=1Vfm=1np(nd,mm˙d,mud,i,mnd,mFd,i,m)\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*U0np0*dUTrans。 其中的 dUTransKinematicParcel::calcVelocity 中计算:dUTrans=mpdeltaUdUTrans=-m_pdeltaU。然后文献中的F=mddeltaUdtF=m_d\dfrac{deltaU}{dt}。这是所有 parcel 的合力,即Fi=m=1npFd,i,mF_i=\sum_{m=1}^{n p} F_{d,i,m},这里的ii 表示力的三个方向,mm 表示对所有 parcel 循环,np 是当前网格 parcel 的个数。nd,mn_{d,m} 表示第 m 个 parcel 中 particle 的个数,cotmd,m\cot m_{d,m} 表示第 m 个 parcel 的蒸发速率,ud,mu_{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()=(ndm˙udF)=-(n_d\dot m u_d - F)

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