【譯】流體模擬
實現雙密度鬆弛算法(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模擬自由落體