实现双密度松弛算法(double density relaxation algorithm)

点击查看原文

不知何故，我一直对流体模拟着迷。它们因物理、编程和创造力的相遇而诞生。它们的特征非常令人着迷：在受到的外力非常大时才会分散成小的液滴，否则会表现为一个整体。飞溅到四周后又再次恢复平静，等待下一次受力。

在此之前我已经多次尝试模拟流体，但从未真正成功过。基础概念看似简单，实现起来却难以置信的复杂。

我在互联网中鲜为人知的角落探索，我偶然发现了一篇名为《基于粒子的粘弹性流体模拟》（Particle-based Viscoelastic Fluid Simulation）的论文，作者是Simon Clavet，Philippe Beaudoin和Pierre Poulin。我花了很长一段时间才弄明白这篇论文（至少是其中一部分），但我终于使我的模拟能够运行了。我想和你分享我的实现（其中的要点）。如果你真的对完整的实现感兴趣，你可以去查看原文中标题背景图的源代码。

基础理论

如果我们更仔细地观察流体的行为，就会发现一些重要的特征：

同种流体之间相互吸引。这使得流体成为连贯的整体并产生表面张力效应。如果没有其他力量，流体将形成一个完美的球体。
流体有一个最大密度。流体从高密度区域流向低密度区域（这也被称为压力渐变）。
流体通常（大多数时候）不可被压缩：在某一点挤压它们会使它们流动到其他地方。

如果我们想象一种由大量粒子（分子）组成的流体，基于一些规则，粒子间表现为相互吸引或相互排斥。我们将要看到双密度松弛算法是如何实现这些规则的。

算法

首先，所有的粒子都受到一堆力的的作用。这些力会加速或减速它们。根据它们的速度将所有粒子移动到新的位置（即使这会将它们移动到其他粒子中或容器外）。

之后计算每个粒子的双密度。该密度基于一定距离（相互作用半径）内的相邻粒子的数量。

第一密度被用于计算相互作用半径内相邻粒子整体的吸引和排斥作用。第二密度（接近密度）仅用于推动彼此太靠近的粒子。这种吸引和排斥被称为松弛作用。

具体来说，该算法需要我们在对所有粒子的遍历中进行三个操作。该过程由下列三个步骤组成：

步骤 1 ：

将力作用于粒子并根据它们的速度更新位置
将粒子的位置存储在空间散列映射中，以便于在后面的步骤中进行查询

步骤2：

找到相关（足够近）的相邻粒子
计算粒子所处位置的压力和接近压力
应用松弛作用

步骤3：

将粒子限制在一定范围内
计算它们的新速度

性能评估

在继续深入前，我们要考虑一下性能问题：由于我们将要遍历1000-2000个粒子，使用类型化数组 (Typed Array) 是绝对有必要的。类型化数组只能存储预定义类型的数字，因此浏览器可以极大地优化对它们的访问。

我们将跟踪每个粒子的以下属性：

const state = {
    x: new Float32Array(PARTICLE_COUNT), // x location
    y: new Float32Array(PARTICLE_COUNT), // y location
    oldX: new Float32Array(PARTICLE_COUNT), // previous x location
    oldY: new Float32Array(PARTICLE_COUNT), // previous y location
    vx: new Float32Array(PARTICLE_COUNT), // horizontal velocity
    vy: new Float32Array(PARTICLE_COUNT), // vertical velocity
    p: new Float32Array(PARTICLE_COUNT), // pressure
    pNear: new Float32Array(PARTICLE_COUNT), // pressure near
    g: new Float32Array(PARTICLE_COUNT), // 'nearness' to neighbour
    mesh: [] // Three.js mesh for rendering
};

Mesh 属性必须是一个常规数组，因为它将存储Three.js Mesh对象，这些对象不是数字（当然，使用Three.js实现此算法并不是必需的）。

步骤 1

施力并根据速度更新粒子

我们算法的第一步是根据作用在粒子上的力更新粒子的速度（vx 和 vy），再根据粒子的新速度来刷新粒子的位置。不过在此之前千万不要忘记存储初始位置（oldX 和oldY），我们在后面将会用到这些值。

// Pass 1

for (let i = 0; i < PARTICLE_COUNT; i++) {

    // Update old position
    state.oldX[i] = state.x[i];
    state.oldY[i] = state.y[i];
    applyGlobalForces(i, dt);

    // Update positions
    state.x[i] += state.vx[i] * dt;
    state.y[i] += state.vy[i] * dt;

    // Update hashmap
    const gridX = (state.x[i] / canvasRect.w + 0.5) * GRID_CELLS;
    const gridY = (state.y[i] / canvasRect.h + 0.5) * GRID_CELLS;
    hashMap.add(gridX, gridY, i);

}

这段代码中的 dt 参数是我们用秒做单位的时间步进值。为了达到 60FPS ，这个值应该约为 0.166 秒。我们使用单位时间的位移来衡量速度，因此这里通过加上速度与时间步进值 dt 的乘积来刷新粒子的位置。

applyGlobalForces 方法将外部的力施加到我们的粒子上，使粒子加速。

const applyGlobalForces = (i, dt) => {
    const force = GRAVITY;
    state.vx[i] += force[0] * dt;
    state.vy[i] += force[1] * dt;
};

加速度用单位时间增加的速度来表示。如果我们已经计算出了加速度，我们需要将它乘上时间步进值 dt 以求出增加的速度。因此速度的增量是 acceleration * dt 。

你可以使用下面的公式来计算一个力产生的加速度：acceleration = force / mass。假设所有粒子的质量都为 1 ，速度的增量则可以简单地表示为 force * dt 。

使用HashMap存储粒子的位置

在更新了粒子的位置之后，我们将新的位置存在一个HashMap中。

HashMap是一种对大型点集进行索引（和读取）的高效方式。我们将canvas划分成网格（一个数组），网格的每一个单元都是一个桶（也是一个数组）。在这个例子中使用的是 54*54 的网格。每一帧我们都会清除 HashMap 并将所有的粒子都放入桶中。

为了算出用来存放粒子的桶的索引，我们进行下面的计算：

const index = Math.round(cellX) + Math.round(cellY) * gridCellsInRow

在 HashMap 中使用计算得出的索引进行查询可以返回指定的桶。

该 HashMap 将左上角设置为坐标轴的原点，一格表示一个单位长度。本次模拟中，使用 canvas 的中心作为原点进行渲染，一个像素表示一个单位长度。这意味着我们在对 HashMap 进行累加或者查询前需要进行一次转换。具体到这个例子中，使用的是下面的方法：

const gridX = (state.x[i] / canvasRect.w + 0.5) * GRID_CELLS;
const gridY = (state.y[i] / canvasRect.h + 0.5) * GRID_CELLS;

复用表示桶的数组是一个重要的性能优化项。如果你以 175 次/秒的速度创建数组，这会是一项昂贵的操作。这就是为什么要用 splice 来清除数组的原因。

上述的 HashMap 实现放在了GitHub上，你可以点击查看。

步骤 2

计算密度以及邻近密度

现在我们有一打粒子，且能通过 HashMap 轻松地搜索到。它们的位置一直在被更新，但是目前为止我们仍没有考虑那些使它们表现得像流体的行为。这就是我们要在步骤 2 中完成的。

// Pass 2

for (let i = 0; i < PARTICLE_COUNT; i++) {

    const neighbours = getNeighboursWithGradients(i);
    updateDensities(i, neighbours);

    // perform double density relaxation
    relax(i, neighbours, dt);

}

找到相关（足够近）的近邻

为了计算双密度，首先我们必须找出相互作用半径内的近邻。我们将计算每个近邻的梯度(g)，表示它与粒子的距离。

在这里，梯度表示一个 0 到 1 之间的值。如果近邻有和粒子完全相同的位置 i，那么梯度值就是 1（非常接近）。当远离粒子时，值逐渐趋于零 0 直到距离等于相互作用半径时，最终等于 0（一点也不近）

该过程由 getNeighboursWithGradients 函数完成：

const getNeighboursWithGradients = i => {
  
    const gridX = (state.x[i] / canvasRect.w + 0.5) * GRID_CELLS;
    const gridY = (state.y[i] / canvasRect.h + 0.5) * GRID_CELLS;
    const radius = (INTERACTION_RADIUS / canvasRect.w) * GRID_CELLS;

    const results = hashMap.query(gridX, gridY, radius);
    const neighbours = [];

    for (let k = 0; k < results.length; k++) {

        const n = results[k];
        if (i === n) continue; // Skip itself

        const g = gradient(i, n);
        if (g === 0) continue

        state.g[n] = g; // Store the gradient
        neighbours.push(n); // Push the neighbour to neighbours

    }

    return neighbours;

};

我们在 HashMap 中查找所有在粒子某个半径范围内的桶。

我们对每个找到的近邻计算并存储梯度。由于存放近邻的数组也包含粒子本身，我们必须检查这个特例（if (i = n) continue）。丢弃梯度为 0 的近邻，因为它们不产生作用。将梯度存下来并返回近邻。

下面的函数用来计算梯度：

const gradient = (i, n) => {

    const particle = [state.x[i], state.y[i]]; // position of i
    const neighbour = [state.x[n], state.y[n]]; // position of n
  
    const lsq = lengthSq(subtract(particle, neighbour));
    if (lsq > INTERACTION_RADIUS_SQ) return 0;

    const distance = Math.sqrt(lsq)
    return 1 - distance / INTERACTION_RADIUS;

};

首先，我们计算粒子和近邻间距离的平方（lsq，长度的平方）。从中我们可以通过求 lsq 的平方根从而很容易地找到距离，但由于这是一个昂贵的操作，我们首先检查 lsq 是否小于交互半径的平方（因为距离也是平方）。如果 lsq 更大，近邻就落在交互半径之外，所以我们只返回 0。

然后我们计算粒子间实际的距离，我们会通过除以交互半径将其转化成梯度(g)。当位置相同时，会返回 0，当距离等于交互半径时会返回 1。你可能已经注意到这与我们想要的梯度值完全相反，这就是为什么我们要在返回之前计算 1 - distance / INTERACTION_RATIUS 将值反转。

设法记住梯度函数的结果可能是一个聪明的做法，因为每秒要计算太多太多次昂贵的开平方。

计算压力和邻近压力

下一步，计算在松弛函数中需要使用的压力(p)和邻近压力(pNear)。我们已经收集了所有的相关近邻和它们的梯度值。接下来我们要将他们的重量相加，以算出密度和邻近密度。

计算密度是这个算法至关重要的一个部分。我们将使用两者不同的核函数来分别计算密度（核函数是一种权重函数）：density = g * g and nearDensity = g * g * g。这意味着梯度值接近 0 的粒子几乎不会计入到这些密度中，而值接近 1 的粒子则会计入。对于邻近密度这种效果会更加明显，仅当距离非常近时才有效。

下面的 updatePressure 函数计算并存储了压力：

const updatePressure = (i, neighbours) => {

    let density = 0;
    let nearDensity = 0;

    for (let k = 0; k < neighbours.length; k++) {
        const g = state.g[neighbours[k]]; // Get g for this neighbour
        density += g * g;
        nearDensity += g * g * g;
    }

    state.p[i] = STIFFNESS * (density - REST_DENSITY);
    state.pNear[i] = STIFFNESS_NEAR * nearDensity;

};

这段代码中，STIFFNESS 和 STIFFNESS_NEAR 是两个常量，它们决定了压力效应对流体的影响有多显著。通过乘上 STIFFNESS_NEAR 和 nearDensity.pressureNear 计算出邻近压力恒为正，这意味着它总是施加一个排斥力。这强化了流体的不可压缩性。

对于常规密度，它的工作方式略有不同。在与 STIFFNESS 相乘之前，我们从密度中减去 REST_DENSITY。这意味着当密度高于静止密度时，压力将是正的，从而产生排斥力；低于静止密度时，压力将为负值，从而产生吸引力。实际上，这将导致粒子从高压区域移动到较低压力区域。此外，此部分负责了您可见的所有表面张力效果。

应用松弛作用

最后，来到了算法的核心部分——松弛作用。

const relax = (i, neighbours, dt) => {
  
    const pos = [state.x[i], state.y[i]];
  
    for (let k = 0; k < neighbours.length; k++) {
        const n = neighbours[k];
        const g = state.g[n];

        const nPos = [state.x[n], state.y[n]];
        const magnitude = state.p[i] * g + state.pNear[i] * g * g;

        const direction = unitApprox(subtract(nPos, pos))
        const force = multiplyScalar(direction, magnitude);

        const d = multiplyScalar(force, dt * dt);

        state.x[i] += d[0] * -.5;
        state.y[i] += d[1] * -.5;

        state.x[n] += d[0] * .5;
        state.y[n] += d[1] * .5;
    }

};

在这一步，我们遍历所有的近邻并基于之前计算出的压力 ( p and pNear ) 对其施加一个力。

为了算出这个力的大小，我们首先计算从当前粒子指向邻近粒子的单位向量。单位向量的长度总是 1 。如果你将它乘以 magnitude 值，就会得到一个从当前粒子指向邻近粒子的强度为 magnitude 的力。在我们的计算中，强度 ( magnitude) 再次使用加权压力计算：p * g + pNear * g * g。

由于力作用于粒子的加速度，而我们对实际的位移感兴趣，我们将力乘以 dt ，将加速度转化为速度的增加，然后再乘以 dt ，求出粒子的实际位移 (d)。

最后，我们从粒子 i 的位置减去这个位移的一半，然后将另一半加到相邻的粒子中。如果我们有一个正的压力，这意味着粒子被移动到远离彼此的位置。如果压力是负的，他们就会互相吸引。

我们将一半的位移应用于这两个粒子的原因是为了确保我们不违反牛顿第三定律：

对于每一个动作，都有一个相等和相反的作用。——牛顿

如果不能保持这一规律，将导致液体变得不平衡，使其自发地开始运动。

为了优化这个函数的性能，我使用了一个 unitApprox 函数，它是从 NickVogt 复制的。这是一个足够好的单位矢量近似，避免了昂贵的平方根计算。

步骤 3

将粒子包含在一定范围内

在这个过程中，所要做的就是确保我们的粒子包含在我们的容器中，为它们计算新的速度，并更新实际的 Three.js 网格位置，以反映模拟。

// Pass 3

for (let i = 0; i < PARTICLE_COUNT; i++) {

    // Constrain the particles to a container
    contain(i, dt);

    // Calculate new velocities
    calculateVelocity(i, dt);

    // Update
    state.mesh[i].position.set(state.x[i], state.y[i], 0);

}

包含粒子仅仅是将粒子的位置重置到容器的边界上的行为，如果粒子穿过容器边界的话。我们的容器有一个圆的形状，对此进行边界检查比较容易：

const contain = i => {

    const pos = [state.x[i], state.y[i]];

    if (lengthSq(pos) > canvasRect.radiusSq) {
    
        const unitPos = unit(pos)
        const newPos = multiplyScalar(clone(unitPos), canvasRect.radius)
    
        state.x[i] = newPos[0]
        state.y[i] = newPos[1]

    }

}

首先我们使用平方长度检查粒子是否通过边界，这样我们就不需要计算平方根了。如果它们确实越过了边界，我们取粒子位置的单位向量，并将其乘以包含圆的半径。这将有效地将其精确地放置在圆的边缘。

计算它们的新速度

接下来，我们为粒子计算新的速度：

const calculateVelocity = (i, dt) => {

    const pos = [state.x[i], state.y[i]];
    const old = [state.oldX[i], state.oldY[i]];

    const v = multiplyScalar(subtract(pos, old), 1 / dt);

    state.vx[i] = v[0];
    state.vy[i] = v[1];

};

使这个算法变得更快的一部分原因是，它不需要在每个单独的交互过程中跟踪粒子的速度。通过 1000-2000 个粒子之间的相互作用，你可以想象这将产生大量的相互作用。

相反，我们从当前位置减去我们在模拟步骤开始时存储的位置，以得到粒子移动的实际距离。我们把这个除以我们用来计算新速度的时间步长。它可能没有真实的流体准确，但它看起来很不错！

额外步骤：从容器边界释放粒子

如果你现在继续执行这些步骤，你会注意到一些奇怪的事情：粒子往往堆积在它们的容器的边界上。当我为此想破头之后，我认为我明白了：

一旦粒子在两个时间步长内保持在相同的位置(例如边界)，这将有效地将其速度重置为零。由于容器外没有粒子，因此没有压力将粒子推离边界，导致粒子卡在那里。虽然粒子最终将被容器内其他粒子的吸引而从边界中被拉出，但显而易见的，这种效果看起来并不是很好。

为了解决这个问题，我对包含函数做了一点修改：

const contain = (i, dt) => {

    const pos = [state.x[i], state.y[i]];
  
    if (lengthSq(pos) > canvasRect.radiusSq) {
    
        const unitPos = unit(pos)
        const newPos = multiplyScalar(clone(unitPos), canvasRect.radius)

        state.x[i] = newPos[0]
        state.y[i] = newPos[1]

        const antiStick = multiplyScalar(
            unitPos, 
            INTERACTION_RADIUS * dt
        )

        state.oldX[i] += antiStick[0]
        state.oldY[i] += antiStick[1]

    }

};

将位置约束到容器内之后，旧位置( oldX 和 oldy )将向容器外部(垂直于边界)移动。当速度被更新后，就会产生一个远离边界的净速度。

_

这就是流体模拟。

编码愉快！