Stam J. Real-Time Fluid Dynamics for Games[J]. Journal of Graphics Tools, 2003, 6.

感謝朱老師的指導..本文對論文的實現程式碼進行分析，並試圖補全論文中說得不夠詳細的部分，側重於對solver的分析，opengl的部分比較簡單，附件為已添加註釋的程式碼。

我使用canvas api實現了一個基於web的版本，限於瀏覽器效能，在幀率上不如原文。可以在此連結中瀏覽效果：http://pigcage.tech/MyProjects/stablefluid/

 

一、系統的建立

1、網格系統

這裡使用的是歐氏方法，程式碼中建立了一個N*N的網格（其中N=64，越大的N值意味著越高的效能開銷）。在網格的四周有一層邊界（0與N+1），亦即實際網格的大小為size=(N+2)*(N+2)。對於網格中的每一個單元格（cell），假定N*N網格的邊長為1，那麼每個單元格的邊長為h=1/N。

基於此網格系統，程式使用一些全域性變數來儲存流體的基本屬性，見附件中demo.c的註釋如圖：

時間步長為dt，一般來說，應滿足dt < h * |u|，以保證模擬的精度。

每個單元格對應有密度（density）與速度（velocity）兩個屬性。其中速度為向量，此處為二維模擬，故分為u、v兩個方向表示，此處u指上下方，v指左右方。對於這些屬性，分別使用形如dens、dens_prev形式的兩個陣列， 表示網格該屬性的當前狀態和前一狀態。使用SWAP(x0,x)來交替這兩個陣列以達到狀態的更新，其中x0為上一狀態。這些陣列均為一維陣列，程式中使用形如dens[IX(i , j)]的形式，表示第j行i列的單元格屬性。

2、準備知識

1.為了便於求解，使用到了Splitting方法，大體來說是這樣的，感覺比較顯然就不具體描述了吧..

2.N-S方程

3.順便貼一下物質導數D

二、密度步驟

由於這裡是處理煙霧，密度的高低等價於煙霧的濃度。使用密度資料即可繪製對應的煙霧檢視。對密度的處理分為上圖幾個步驟，在每一幀中，依次新增力、處理擴散和移動。

1、力和源的新增

這一步比較簡單，使用add_source函式。

void add_source ( int N, float * x, float * s, float dt )
{
    int i, size=(N+2)*(N+2);
    for ( i=0 ; i<size ; i++ ) x[i] += dt*s[i];
}

使用者輸入的源（source）或力（force）被儲存在 s[ ] 中，那麼把這些狀態直接加到新一幀的對應屬性 x[ ] 即可。至於乘以dt也就是因為使用了splitting的方法，並不是衝量..下同。

2、擴散（diffusion）

上面已經處理了外部輸入源的情況。對於接下來的擴散步驟（煙霧的自發擴散）。根據N-S方程，我們知道擴散項為 k · ∇^2 u ，其中k為擴散係數，u為密度場。於是我們對某個區域性的區域 x 而言，有 x - x0 = k · ∇^2 u 。這裡是對x、y兩個方向求偏導，也就相當於只考慮每個cell與其四領域的互動了。

設每一個幀的時間dt，擴散係數diff，對於上式我們有：

x[IX(i,j)] - x0[IX(i,j)] = a*(x[IX(i-1,j)]+x[IX(i+1,j)]+x[IX(i,j-1)]+x[IX(i,j+1)]-4*x[IX(i,j)]);

其中a=dt*diff*N*N，由於這裡和後面projection用到了相同的方法，我想把對這條式子右值的解釋留到最後。對上式進行化簡整理，可以得到：

x[IX(i,j)] = (x0[IX(i,j)] + a*(x[IX(i-1,j)]+x[IX(i+1,j)]+x[IX(i,j-1)]+x[IX(i,j+1)]))/(1+4*a);

令t=a/(1+4*a)，陣列x0為已知的上一狀態。我們可以把這個等式寫成線性方程組，如下圖：

（勘誤：上圖-1應為-t/a）

對於最左矩陣的第IX(i,j)行而言，對應的四個t的位置分別為第IX(i-1,j)，IX(i,j-1)，IX(i,j+1)，IX(i+1,j)列，而-t/a的位置為第IX(i,j)列。這是一個稀疏矩陣，求解這個線性方程組用到了高斯-賽德爾迭代。關於這個方法的本身可參考https://www.cfd-online.com/Wiki/Gauss-Seidel_method。程式碼實現如下：

for ( k=0 ; k<20 ; k++ ) {
    for ( i=1 ; i<=N ; i++ ) {
        for ( j=1 ; j<=N ; j++ ) {
           x[IX(i,j)] = (x0[IX(i,j)] + a*(x[IX(i-1,j)]+x[IX(i+1,j)]+x[IX(i,j-1)]+x[IX(i,j+1)]))/(1+4*a);
        }
    }
    set_bnd ( N, b, x );
}

其中k為迭代的次數，此處設為20，較大的迭代次數可以保證好的精度，但效能的開銷也更大。set_bnd()函式的作用為設定邊界條件，後面再具體分析。

3、對流（advection）

由於速度場的存在，流體不僅有自發的擴散運動，也有因速度帶來的流動。按照通常的思維角度，也就是N-S方程的對流項 - (u · ∇) ρ，我們可以根據每個單元格當前的u、v分量（下圖a）計算出步長dt後的對應位置（下圖b）。但這個位置並不是總在單元格的中心點上（下圖b），這將導致計算出的結果難以使用。

因此，作者在此處採取了前向追溯的方式，使用半拉格朗日法（semi-lagrangian）通過當前單元格的速度來計算該處狀態在前一dt時的位置（下圖c）。這裡討論的是密度的對流，對於速度的對流也是一樣的道理，只是速度是個向量場，拆分成多個維度來計算。

實現的程式碼如下：

void advect ( int N, int b, float * d, float * d0, float * u, float * v, float dt )
{
    int i, j, i0, j0, i1, j1;
    float x, y, s0, t0, s1, t1, dt0;
    dt0 = dt*N;
    FOR_EACH_CELL
        //線性回溯上一時刻的位置(x,y)
        x = i-dt0*u[IX(i,j)]; y = j-dt0*v[IX(i,j)]; 

        if (x<0.5f) x=0.5f; if (x>N+0.5f) x=N+0.5f; i0=(int)x; i1=i0+1;
        if (y<0.5f) y=0.5f; if (y>N+0.5f) y=N+0.5f; j0=(int)y; j1=j0+1; 
        
        //四個鄰居到點(x,y)線性插值，求得點(x,y)的值
        s1 = x-i0; s0 = 1-s1; t1 = y-j0; t0 = 1-t1;
        d[IX(i,j)] = s0*(t0*d0[IX(i0,j0)]+t1*d0[IX(i0,j1)])+s1*(t0*d0[IX(i1,j0)]+t1*d0[IX(i1,j1)]);
    END_FOR
    set_bnd ( N, b, d );
}

這裡主要解釋一下d、d0的等量關係。

前一位置(x,y)可能位於網格的任意一點上，由於其不一定出於中心點，我們無法直接確定該點的密度。作者提出通過線性插值的方式，使用點(x,y)相鄰四個中心點的密度值，來算出點(x,y)處的密度。具體如下圖：

其中，紅點處為點(x,y)，相鄰的四個網格內的黑點為中心點，左上角網格為假定點(x,y)經過時間步長dt後到達的當前位置(i,j)，s0、s1、t0、t1為點(x,y)到四個相鄰中心點的橫縱距離。根據線性插值的規則不難得出(x,y)處的實際密度，而此處的密度，經過dt時長後到達點(i,j)處，故實際上就是當前中心點(i,j)的密度，如上述程式片段所示。對於後面速度步驟中，處理u、v分量，亦是相同的原理。

至此，對當前一幀密度狀態的處理完成，如下：

void dens_step ( int N, float * x, float * x0, float * u, float * v, float diff, float dt )
{
    add_source ( N, x, x0, dt );
    SWAP ( x0, x ); diffuse ( N, 0, x, x0, diff, dt );
    SWAP ( x0, x ); advect ( N, 0, x, x0, u, v, dt );
}

三、速度步驟

對速度的處理與對密度的處理大致類似，由於速度有u、v兩個方向的分量，因此對每一分量要經過同樣的處理。不同的是，對速度的處理需要經歷投影（project）步驟，得到一個不可壓縮場（質量守恆場），實際上就是處理N-S方程中的壓力項。

首先有以下概念：

至於為什麼是這樣來呢，以及project部分的詳細數學過程，我想留到最後解釋。這裡Δx = Δy =1/N = h。根據霍奇分解定理（Helmholtz-Hodge Decomposition），速度場w可以被分解為不可壓場u與梯度場∇p的和：

w = u + ∇p

u是不可壓場，有∇ · u = 0，對於等式(1)兩邊分別進行 ∇ · 操作，得到：

∇ · w = Δ p

離散化：

等式左側實際上就是散度場div[]，我們可以通過已知的u、v分量求得，程式碼如下。

for ( i=1 ; i<=N ; i++ ) {
    for ( j=1 ; j<=N ; j++ ) {
        div[IX(i,j)] = -0.5*h*(u[IX(i+1,j)]-u[IX(i-1,j)]+v[IX(i,j+1)]-v[IX(i,j-1)]);
        p[IX(i,j)] = 0;
    }
}

於是有：

類似前面在advection中所述，此處的式(1)同樣可以寫成線性方程組的形式，不再列出。因而這裡同樣可以用高斯-賽德爾迭代求解壓力場 p[] ，如下：

for ( k=0 ; k<20 ; k++ ) {
    for ( i=1 ; i<=N ; i++ ) {
        for ( j=1 ; j<=N ; j++ ) {
            p[IX(i,j)] = (div[IX(i,j)]+p[IX(i-1,j)]+p[IX(i+1,j)]+p[IX(i,j-1)]+p[IX(i,j+1)])/4;
        }
    }
    set_bnd ( N, 0, p );
}

解出p以後，就可以算出梯度場 ∇p ，最後用速度場 w - ∇p = u，得到不可壓場。

for ( i=1 ; i<=N ; i++ ) {
    for ( j=1 ; j<=N ; j++ ) {
        u[IX(i,j)] -= 0.5*(p[IX(i+1,j)]-p[IX(i-1,j)])/h;
        v[IX(i,j)] -= 0.5*(p[IX(i,j+1)]-p[IX(i,j-1)])/h;
    }
}

至此，project步驟完成。

四、邊界條件

所謂邊界，就是網格最外層的一圈（0和N+1）。在計算網格的某一cell的具體狀態時，我們通常要使用到與之相鄰的幾個cell的狀態，那麼對位於網格邊緣的cell而言，它們在一些方向上自然是沒有鄰居的。為了解決這個問題，我們在網格的外面額外增加一層邊界，並根據實際情況為其賦值。

先看實現：

void set_bnd ( int N, int b, float * x )
{
    int i;

    for ( i=1 ; i<=N ; i++ ) {
        x[IX(0  ,i)] = b==1 ? -x[IX(1,i)] : x[IX(1,i)];
        x[IX(N+1,i)] = b==1 ? -x[IX(N,i)] : x[IX(N,i)];
        x[IX(i,0  )] = b==2 ? -x[IX(i,1)] : x[IX(i,1)];
        x[IX(i,N+1)] = b==2 ? -x[IX(i,N)] : x[IX(i,N)];
    }
    x[IX(0  ,0  )] = 0.5f*(x[IX(1,0  )]+x[IX(0  ,1)]);
    x[IX(0  ,N+1)] = 0.5f*(x[IX(1,N+1)]+x[IX(0  ,N)]);
    x[IX(N+1,0  )] = 0.5f*(x[IX(N,0  )]+x[IX(N+1,1)]);
    x[IX(N+1,N+1)] = 0.5f*(x[IX(N,N+1)]+x[IX(N+1,N)]);
}

引數b等於1、2時分別處理u、v兩個方向的分量，也就是該速度方向到達對應的兩個邊界時發生碰撞，速度方向取反。當b為0時，處理的是密度值，密度不存在碰撞問題，故四個邊界均保留原值。

對於網格的四個頂點而言，這裡的狀態為其相鄰兩個單元格的平均值。

五、完整過程

最後放一下實際執行模擬的程式碼。這裡實際是先處理速度再處理密度，因為需要保證其不可壓縮性。

while ( simulating )
{
    get_from_UI ( dens_prev, u_prev, v_prev );
    vel_step ( N, u, v, u_prev, v_prev, visc, dt );
    dens_step ( N, dens, dens_prev, u, v, diff, dt );
    draw_dens ( N, dens );
}

六、散度、梯度及其它

現在我們來把散度、梯度以及上面提及的使用高斯-賽德爾迭代求解的等式具體解釋一下。

首先說散度（divergence），某位置上的散度形容的是這個位置上向量場的發散程度。舉個可能不太恰當的例子，你周圍靜止地站著一些人，這時候你這個位置的散度div=0；由於你今天噴了香水，周圍的人紛紛向你湊過來聞，這時候你的散度就是負的了；然後你又放了一個屁，大家不堪其臭紛紛散開，於是你又成了一個正源，即div > 0。

回到流體中，對於一個可壓縮的流體來說，某點周圍的物質由於速度的存在，可能總體會朝這個點聚攏，這時候div<0，表現為該點附近的物質被壓縮了；類似地，當div>0時，物質紛紛離開這個點，看起來就是這裡的物質膨脹或者爆炸了。一般而言，我們研究的是不可壓流體，或者換句話說，我們希望這個流體看起來是質量守恆的（mass conserving），物質不會憑空壓縮或爆炸。也就是，我們需要保準散度為0，物質不能壓縮也不能爆炸。更準確地說，流體的體積應該是不變的，而控制這一過程的正是壓力p，也就是N-S方程的壓力項。

現在來看散度運算元 。div = ∇ · w ，其中∇ · 是散度運算元，w是一個梯度場，運算的結果是純量。我們討論二維的情況。令w=( u , v )， ∇ · w = ∂u/∂x +∂v/∂y 。中心差分法，對二維平面上的速度場u( i , j )、v( i , j ) 有：

∂u/∂x = ( u(i+1,j) - u(i-1,j) )/2Δx

同理有：

∂v/∂y = ( v(i,j+1) - v(i,j-1) )/2Δy

於是我們得到散度計算公式：

 ∇ · w = ∂u/∂x +∂v/∂y = ( u(i+1,j) - u(i-1,j) )/2Δx + ( v(i,j+1) - v(i,j-1) )/2Δy

已知Δx=Δy=h=1/N，得到：

散度div =  ( u(i+1,j) - u(i-1,j) + v(i,j+1) - v(i,j-1) ) / 2h

2. 接下來如何去除速度場的散度，使得 ∇ · w = 0，即實現質量守恆、不可壓性呢？

根據亥姆赫茲 - 霍奇分解定理，任一速度場 w（向量場），可以唯一分解為一個無散度場u與一個梯度場的和。即：

w = u +  ∇ p

所謂projection，也就是去掉∇ p，即取 w 在 u 上的投射。

那麼對於上式，u是無散場，∇ · u = 0，對等式兩邊∇ · 有：

∇ · w = Δ p

對於此式，∇ · w = ∂u/∂x +∂v/∂y ，此處已知，也就是我們前面求得散度場div。

對於等式右邊，Δ p = ∇ · ∇ p = ∂(∂p/∂x) / ∂x +∂(∂p/∂y) / ∂y ，不妨以x方向來舉例，設x方向上有三個離散的點 x(i-1) ， x(i) ， x(i+1) ，此外x(i-1) ， x(i)中點處有一點a，  x(i) ， x(i+1) 中點處有一點b。

中心差分法，對於a處我們有 f' (a) = [ f(i) - f(i-1) ]/h，同理對b有f' (b) = [ f(i+1) - f(i) ]/h

那麼，對於i處有f'' (i) = [ f' (b) - f' (a) ] / h = [ f(i+1) - 2f(i) + f(i-1) ] / h^2。令 ∇ · w = Δ p = div，於是：

於是得到p的迭代公式：

這裡注意，因為原程式中lin_solve()函式的寫法關係，上式div對應的變數係數只能為1，因此程式在計算散度div時實際計算的時d = - h^2 · div，再代入上式計算。於是有：

至此可以用高斯迭代計算出壓力場p。得到p後，求壓力梯度場 ∇ p = ( ∂p/∂x , ∂p/∂y )，最後速度場減去此梯度場，得到不可壓場。

 

七、附件

附件為添加了一定註釋的完整實現程式碼demo.c與solver.c，由於為邊分析邊添加註釋，難免錯漏。本文亦出於整理思路的需要寫成，錯漏之處還望指正。

附件：

solver

demo