javascript - 我对这种流体动力学计算有什么误解？

Question

最后一个问题因为太短太含糊而被删除，所以这次我会尽量冗长。

在本文之后，我正在使用 Javascript 来实现一个基本的物理上准确的 SPH 。一开始我想要的只是让基础工作，只是为了通过第一个算法。

我从 (0.05, 0.05) 和 (0.06999, 0.05) 两个粒子的测试案例开始。他们的样本半径刚刚接触，我已经得到了奇怪的数字。

这就是我计算密度和压力的方式：

for (var i = 0; i < particles.length; i++)
{
    var pi = particles[i];

    pi.density = 0;

    for (var j = 0; j < pi.neighbors.length; j++)
    {
        var n  = pi.neighbors[j];
        var pj = n.p;

        var h  = (pi.sampleDist + pj.sampleDist) / 2;

        n.r    = subv(pi.pos, pj.pos); // subv() subtracts first vector from second
        n.dist = length(n.r);
        n.q    = n.dist/h;

        pi.density += pj.mass * W(n.q, h, 2);
    }

    pi.pressure = pi.stiffness * (Math.pow(pi.density / pi.restDensity, 7) - 1);
}

使用平滑内核

function W(q, h, d)
{
    var f = 0;

         if (0 <= q && q < 1) f = 2/3 - sqr(q) + 1/2*cube(q); // 2/3 - q² + 1/2q³
    else if (1 <= q && q < 2) f = 1/6 * cube(2 - q);          // 1/6(2-q)³

    return 3/Tau / Math.pow(h, d) * f;
}

这些是我正在使用的设置：

质量：0.001 kg
h：0.01 m
静止密度：1000 kg/m³
刚度：0.2

在只有两个具有这些设置的粒子的情况下，每个粒子的密度变为

pi.density = pj.mass * W(n.q, h, 2);= 0.001 * W(1.99899, 0.01, 2) = 0.001 * 0.00000079577 = 0.00000000079...

这个数字对我来说很有意义——周围几乎什么都没有，所以密度几乎为零。

但现在我要计算压力，我这样做是这样的：

这是我的实现：

var Fp = {x:0, y:0};

if (pi.density !== 0)
{
    var Vp = valueGradient(pi, 'pressure');

    Fp.x = -pi.mass/pi.density * Vp.x;
    Fp.y = -pi.mass/pi.density * Vp.y;
}

function valueGradient(pi, A) // A is a scalar
{
    var vA = {x: 0, y: 0};

    for (var j = 0; j < pi.neighbors.length; j++)
    {
        var n  = pi.neighbors[j];
        var pj = n.p;
        var h  = (pi.sampleDist + pj.sampleDist) / 2;

        var f = pj.mass / pj.density * pj[A];

        var vw = VW(n, h);

        vA.x += f * vw.x;
        vA.y += f * vw.y;
    }

    return vA;
}

function VW(n, h)
{
    var dw = 
          dW(n.q, h, 2)
        / (n.dist * h);

    return {
        x: dw * n.r.x,
        y: dw * n.r.y };
}


function dW(q, h, d) // derivative of kernel
{
    var f = 0;

         if (0 <= q && q < 1) f = -2*q + 3/2*sqr(q); // -2q + 3/2q²
    else if (1 <= q && q < 2) f = -1/2 * sqr(2 - q); // -1/2(2-q)²

    return 3/Tau / Math.pow(h, d) * f;
}

但是出了点问题，因为当粒子的采样半径刚刚接触时，排斥力是巨大的，当它们靠近时，它会变为零。我明白为什么会发生这种情况，因为根据我的设置，到处都有很多 0.000...，而且密度如此之小，除以所有这些 0.000... 最终会得到巨大的数字。所以为什么会有这些数字是有道理的，所以我怀疑我误解了一些基本的东西。

这感觉……倒退了。我犯了什么错误？我一直在看内核导数...

另一件事是，力并不是严格对称的，它们倾向于保持对角线，直到其中一个坐标的差值大约为零，然后它们突然翻转。

我觉得我在某个地方搞砸了微分方程......

可以单步执行此代码的实时版本位于www.bourt.com/particles。有问题的代码在simulation.js、physics.js和kernel.js中。

标签： javascriptmathsimulationfluid