二維的NS-方程:
這個方程一定要拆分成部分才能解出來。
這里面我感覺只要把泊松方程解法搞定，基本快出山寫最簡單的 "計算流體" 完全沒問題
以下是做了個初始的source field,用python numpy 先快速擼了一遍算法。
並且重新實現3d版本在Houdini中，houdini有更好的可視化.
velocity advection 是RK2,
density advection 一階的,
然后發現density field確實跟不上velocity field
從零蛋編寫無粘流體最簡單的路線:
<0>為什么要使用MAC-GRID:
按照作者的說法，使用這種grid,最大的優勢是解決了中心差分法零空間問題。
速度場一定要存在face center.
所以在Houdini看vel field 是一定存在grid face center上的,而不是grid cell center!
而density,temperature....這類場是存在grid cell上的。
總結下一個inviscid fluid(歐拉無粘流體)實現過程:
<1>對流。
最重要的是理解材質導數.為什么這個玩意能等於0? 詳見歐拉與拉格朗日的觀點。
這個差分法是萬萬不可取，因為會出現null-space,試着想想如果有3個點，中間的點高一點，兩邊一樣高，那么用中心差分法算出來中間點的斜率將是0.
所以使用無條件穩定的 半拉格朗日,這個其實完全是解ODE問題了。而不是PDE了
比如要得到xp的值，只需用用 向前euler 向后追蹤法。當然這個方法比較廢物。作者推薦RK家族的算法比如:RK2:
def RK2(x, y, dt, ugrid, vgrid): nu = neibours_value(ugrid, x, y, "u") nv = neibours_value(vgrid, x, y, "v") xmid = x - 1.0 / 2.0 * dt * nu ymid = y - 1.0 / 2.0 * dt * nv umid = neibours_value(ugrid, xmid, ymid) vmid = neibours_value(vgrid, xmid, ymid) x -= dt * umid y -= dt * vmid return (x, y)
(你去看houdini上面的gas advect 還有RK5)
得到位置直接用bilerp()插值法求出xp的量。
def bilerp(f00, f10, f01, f11, tx, ty): """ FIGURE : first lerp in top x,then bottom x, then along y axis f00*----------.tx-----------*f10 | | | | | | | .ty | | | | | | | f01*----------.tx-----------*f11 """ return lerp(lerp(f00, f10, tx), lerp(f01, f11, tx), ty)
那么這個量就是作為下一個時間步的量。
測試這個最簡單的就是創建一個簡單的恆定區域網格速度，然后讓自己的初始的density是不是根據semi-lagrangian方法能advection.
<1.1>推出壓力方程:
梯度壓力方程:
離散，下面有具體的壓力梯度離散過程。
這個為壓力梯度方程。只要求出壓力p就可以得到無散速度場。
接下來不可壓縮流體條件:
使用中心差分法離散。這樣沒有null-space,因為速度場特殊的儲存方式
接下來推出怎么樣得到下一個時間步上的流體速度場是無散度,這里有個技巧，我們不能直接使用5.4式。但是:
首先把速度散度公式寫成下一步時間的離散形式:
把梯度壓力方程帶入可以得到:
觀察此方程，右邊就是不可壓縮 負的速度梯度, 這部分叫做rhs,右手方程,這個是已知的。
所以方程中只有壓力p未知。求出壓力p即可。
<2>求壓力右手方程。通常就叫做RHS,就是負的速度場的散度。-divergence(u)
這個是簡單的。
<3>求解壓力.
這個有好多方法求，由於是一個AX=b 的形式，大部分文章是以PCG/GAUSS-SEIDEI/GAUSS-SEIDEI SOR/JACOBI方式
由於在numpy方便，我直接把方程用切片法做了。
<4>把求的pressure field 帶入壓力梯度更新方程，求出無散度的速度。
Pressure gradient 壓力梯度方程離散:
這個方程我覺得才是套路中的套路。有了這個壓力，直接帶入這個公式 ，就可以求出無散的流體。
