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)Re<10000.424Re>1000C_D=\left\{\begin{array}{ll}{\dfrac{24}{Re_{p}}\left(1+\dfrac{1}{6} R e_{p}^{2 / 3}\right)} & {R e<1000} \\ {0.424} & {R e>1000}\end{array}\right.

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

可以推导出:

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

τu=8mpπρCDD2upu=43ρpDρCpupu\tau_{\mathrm{u}} =\dfrac{8 m_{p}}{\pi \rho C_D D^{2}\left|\mathbf{u}_{p}-\mathbf{u}\right|} =\dfrac{4}{3} \dfrac{\rho_{p} D}{\rho C_{p}\left|\mathbf{u}_{p}-\mathbf{u}\right|}

其中ada_d 是粒子的加速度,ρd\rho_d 是粒子的密度,DD 是粒子的直径。

可以变形为:

dupdt=upuτu=upu43ρpDρCpupu=3(upu)ρCDupu4ρpD=3(upu)μgCD4ρpD2ρupuDμg=34μgCDRepρpD2(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}{\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}\\ =-\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}{\mu_g}\\ = -\dfrac{3}{4}\dfrac{\mu_g C_D Re_p}{\rho_p D^2}(\mathbf{u}_{p}-\mathbf{u}) \end{array}

还原到牛顿第二定律中:

Fp=mpdupdt=mp34μgCDRepρpD2(uup)=3mpμgCDRe4ρpD2(uup)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^2}(\mathbf{u}-\mathbf{u}_{p})\\ =\dfrac{3m_p\mu_gC_DRe}{4\rho_p D^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ρD2\dfrac{3m_p\mu C_DRe}{4\rho D^2}

然后再回到 KinematicParcel::calcVelocity

Fcp=3mpμCDRe4ρD2Fncp=0massEff=mpacp=UcFcp/mpancp=0bcp=Fcp/mpdeltaU=Fcpmp(ucud)dt=Fcpmp(ucud)dt=3μCDRe4ρD2(ucud)dtdeltaUncp=0deltaUcp=deltaUdUTrans=mddeltaUSpu=dtFcp\begin{array}{ll} &Fcp=\dfrac{3m_p\mu C_DRe}{4\rho D^2}\\ &Fncp=0\\ &massEff=m_p\\ &acp=U_c Fcp/m_p\\ &ancp=0\\ &bcp=Fcp/m_p\\ &deltaU=\dfrac{Fcp}{m_p}(\mathbf{u}_c-\mathbf{u_d})dt=\dfrac{Fcp}{m_p}(\mathbf{u}_{c}-\mathbf{u_d})dt=\dfrac{3\mu C_DRe}{4\rho D^2}(\mathbf{u}_{c}-\mathbf{u_d})dt\\ &deltaUncp=0\\ &deltaUcp=deltaU\\ &dUTrans=-m_ddeltaU\\ &Spu=dt*Fcp \end{array}

可以发现这里的 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˙iS˙i=1Vfm=1np(nd,mm˙d,mud,i,mnd,mFd,i,m)\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*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)SU.source()=-(n_d\dot m u_d - F)

这里的负号可能是因为这篇博客里提到的。

add one line in stackedit.io!

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