DPM粒子受力分析

upc_ngh

在DPM中粒子所受力例如重力、虚拟质量力等都是通过forceSuSp中的vector Su以及scalar Sp来决定的，最后通过KinematicParcel.C来计算Fcp以及Fncp(这个地方是每个力单独耦合求解速度还是把所有的力都考虑进来一起耦合),不知道理解的对不对，现在我想要把每个粒子的单独受力比如重力、虚拟质量力等都单独写出来，不知道各位老师有没有什么办法去实现？Ps:我看另一个帖子中写到F=Su*(U-Ud)+Sp.这个是在哪里出现的我也没找到，Ud是啥，可以直接用吗？

这是单独一个力的以重力为例：

template<class CloudType>
Foam::forceSuSp Foam::GravityForce<CloudType>::calcNonCoupled
(
    const typename CloudType::parcelType& p,
    const scalar dt,
    const scalar mass,
    const scalar Re,
    const scalar muc
) const
{
    forceSuSp value(Zero, 0.0);

    value.Su() = mass*g_*(1.0 - p.rhoc()/p.rho());//这里求出来Su

    return value;
}

这是KinematicParcel.C出现通过力算速度的：

const forceType& forces = td.cloud().forces();
//这是最后粒子受力耦合求解速度
    // Momentum source due to particle forces
    const parcelType& p = static_cast<const parcelType&>(*this);
//这个地方是每个力单独耦合求解速度还是把所有的力都考虑进来一起耦合？？
    const forceSuSp Fcp = forces.calcCoupled(p, dt, mass, Re, mu);
    const forceSuSp Fncp = forces.calcNonCoupled(p, dt, mass, Re, mu);
    const forceSuSp Feff = Fcp + Fncp;
    const scalar massEff = forces.massEff(p, mass);
    // New particle velocity
    //~~~~~~~~~~~~~~~~~~~~~~
    // Update velocity - treat as 3-D
    const vector abp = (Feff.Sp()*Uc_ + (Feff.Su() + Su))/massEff;
    const scalar bp = Feff.Sp()/massEff;

    Spu = dt*Feff.Sp();

    IntegrationScheme<vector>::integrationResult Ures =
        td.cloud().UIntegrator().integrate(U_, dt, abp, bp);

    vector Unew = Ures.value();