拉格朗日求解器一个bug

李东岳

在计算粒子的Eo的时候，OpenFOAM定义

template<class ParcelType>
inline Foam::scalar Foam::KinematicParcel<ParcelType>::Eo
(
    const vector& a,
    const scalar d,
    const scalar sigma
) const
{
    vector dir = U_/(mag(U_) + ROOTVSMALL);
    return mag(a & dir)*(rho_ - rhoc_)*sqr(d)/(sigma + ROOTVSMALL);
}

注意上面的(rho_ - rhoc_)，如果粒子密度小于连续相密度，在计算升力的时候：

scalar Eo = p.Eo(g, p.d(), sigma_);
scalar dH = p.d()*cbrt(1.0 + 0.163*pow(Eo, 0.757));

对于Eo的负值，pow(Eo, 0.757)为虚数

CFD中文网

拉格朗日求解器一个bug