【译】流体模拟
实现双密度松弛算法(double density relaxation algorithm)
点击查看原文javascript
不知何故，我一直对流体模拟着迷。它们因物理、编程和创造力的相遇而诞生。它们的特征很是使人着迷：在受到的外力很是大时才会分散成小的液滴，不然会表现为一个总体。飞溅到四周后又再次恢复平静，等待下一次受力。java
在此以前我已经屡次尝试模拟流体，但从未真正成功过。基础概念看似简单，实现起来却难以置信的复杂。git
我在互联网中不为人知的角落探索，我偶然发现了一篇名为《基于粒子的粘弹性流体模拟》（Particle-based Viscoelastic Fluid Simulation）的论文，做者是Simon Clavet，Philippe Beaudoin和Pierre Poulin。我花了很长一段时间才弄明白这篇论文（至少是其中一部分），但我终于使个人模拟可以运行了。我想和你分享个人实现（其中的要点）。若是你真的对完整的实现感兴趣，你能够去查看原文中标题背景图的源代码。github
基础理论
若是咱们更仔细地观察流体的行为，就会发现一些重要的特征：算法
- 同种流体之间相互吸引。这使得流体成为连贯的总体并产生表面张力效应。若是没有其余力量，流体将造成一个完美的球体。
- 流体有一个最大密度。流体从高密度区域流向低密度区域（这也被称为压力渐变）。
- 流体一般（大多数时候）不可被压缩：在某一点挤压它们会使它们流动到其余地方。
若是咱们想象一种由大量粒子（分子）组成的流体，基于一些规则，粒子间表现为相互吸引或相互排斥。咱们将要看到双密度松弛算法是如何实现这些规则的。编程
算法
首先，全部的粒子都受到一堆力的的做用。这些力会加速或减速它们。根据它们的速度将全部粒子移动到新的位置（即便这会将它们移动到其余粒子中或容器外）。canvas
以后计算每一个粒子的双密度。该密度基于必定距离（相互做用半径）内的相邻粒子的数量。数组
第一密度被用于计算相互做用半径内相邻粒子总体的吸引和排斥做用。第二密度（接近密度）仅用于推进彼此太靠近的粒子。这种吸引和排斥被称为松弛做用。浏览器
具体来讲，该算法须要咱们在对全部粒子的遍历中进行三个操做。该过程由下列三个步骤组成：性能优化
步骤 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 )将向容器外部(垂直于边界)移动。当速度被更新后，就会产生一个远离边界的净速度。
_
这就是流体模拟。
编码愉快！
- 1. 虚幻4流体模拟
- 2. Embergen 流体模拟工具
- 3. 流体模拟学习——基于物理的流体动画
- 4. 计算流体模拟理论2
- 5. 流体模拟边界问题
- 6. solidworks流体模拟分析概述
- 7. 模拟费用流
- 8. unity 模拟三体运动
- 9. 模拟体育竞技
- 10. JavaScript模拟自由落体
- 更多相关文章...
- • Eclipse 编译项目 - Eclipse 教程
- • XSL-FO 流 - XSL-FO 教程
- • 委托模式
- • ☆基于Java Instrument的Agent实现
-
每一个你不满意的现在，都有一个你没有努力的曾经。
- 1. 虚幻4流体模拟
- 2. Embergen 流体模拟工具
- 3. 流体模拟学习——基于物理的流体动画
- 4. 计算流体模拟理论2
- 5. 流体模拟边界问题
- 6. solidworks流体模拟分析概述
- 7. 模拟费用流
- 8. unity 模拟三体运动
- 9. 模拟体育竞技
- 10. JavaScript模拟自由落体