使用 OpenFPM 的 SPH 密度求和结果的差异

问：

我正在使用 OpenFPM 框架进行平滑粒子流体动力学 （SPH） 模拟。我建立了一个 2D 无量纲化的 Poiseuille 流动模拟，并观察到密度求和结果中存在微小但一致的差异。

问题描述：

颗粒间距为 dp = 4，导致颗粒体积为 V = dp^2 = 16。理想情况下，SPH 总和应为 1/16 = 0.06250。 但是，我的结果略有偏差：

kernel_summation : 0.062503951537163940366248482405354

内核实现： 我正在使用五项样条内核：

inline double densitySummation(const double r_)
    const double ALPHA_2D = 7.0 / (478.0 * M_PI);
    double q = r_ / H;
    double result;
    if (q <= 1.0)
        result = std::pow((3.0 - q),5) - 6.0 * std::pow((2.0 - q),5) + 15.0 * std::pow((1.0 - q),5);
    else if (q <= 2.0)
        result = std::pow((3.0 - q),5) - 6.0 * std::pow((2.0 - q),5);
    else if (q <= 3.0)
        result = std::pow((3.0 - q),5);
    else
        result = 0.0;
    return ALPHA_2D * result / (H * H);

这是我用来验证密度求和的函数。我正在设置域并迭代它，并进行此求和，没有任何移位，但是求和一直存在错误。

inline void densitySummationNoCorrection(particles &vd, myCelllist &NN)
{
    // particle iterator
    auto part = vd.getDomainIterator();

    // For each particle ...
    while (part.isNext())
    {
        auto a = part.get();
        if (vd.template getProp<type>(a) == FLUID_PARTICLES)
        {
            Point<2, double> xa = vd.getPos(a);
            double massa = (vd.getProp<type>(a) == FLUID_PARTICLES) ? MassFluid : MassBound;
            auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.getPos(a)));

            double kernel
            double summationInLoop = 0.0;
            while (Np.isNext())
            {
                auto b = Np.get();
                double massb = (vd.getProp<type>(a) == FLUID_PARTICLES) ? MassFluid : MassBound;
                double rhob = vd.template getProp<rho>(a);

                // Get particle properties b
                Point<2, double> xb = vd.getPos(b);
                Point<2, double> dr = xa - xb;
                Point<2, double> DW;
                double r_ = sqrt(norm2(dr));
               

                double wab = Wab(r_);
                kernel += wab
                

                ++Np;
            }

            vd.template getProp<Volume>(a) = kernel;
            vd.template getProp<rho_densitySummation>(a) = massa*kernel;
           
        }
        ++part;
    }
}

常数：

constexpr double a2_ = 7.0 / 478.0 * oneOverPi * one_over_H * one_over_H;
constexpr double a2_grad = -5.0 * a2_ * one_over_H;

验证： 我使用另一个仿真套件验证了仿真参数，它们按预期工作。域设置和粒子间距显示正确，精确间距为 4.000。

我对这个微小的差异感到困惑。虽然它很小，但它是一致的，我很想了解它的起源。任何见解或建议将不胜感激。

先谢谢你！

最好 伊凡

如上所述