【計算機圖形學】流體模拟渲染基礎

流體模拟渲染基礎

• 前言
• 矢量微積分
• Naiver-Stokes偏微分方程組
• N-S方程的分步求解
• 對流算法

前言

  本文主要參考文獻《FLUID SIMULATION SIGGRAPH 2007 Course Notes》，結合我的了解單純地講述一下流體渲染的一些基礎知識，本人水準有限，如有錯誤，歡迎指出。本文隻是單純針對流體模拟領域，可能一些地方不太嚴謹，但是對于虛拟模拟來說是可行的。即便如此，本文涉及到大量的數學方法。

一、矢量微積分

  高等數學中太多數讨論的是一維的微積分，而矢量微積分則是一維微積分的高維擴充。矢量微積分的三個基礎算子：梯度（符号為 ∇ ∇ ∇），散度（符号為 ∇ ⋅ ∇\cdot ∇⋅)，旋度（符号為 ∇ × ∇\times ∇×），在此基礎上流體力學中經常用到的還有拉普拉斯算子。

1、梯度（Gradient）

  梯度實際上就是矢量的空間偏導數，且結果依然是一個矢量， 2 2 2維的梯度如下：

(1.1) ∇ f ( x , y ) = ( ∂ f ∂ x , ∂ f ∂ y ) ∇f(x,y)=(\frac{\partial f}{\partial x},\frac{\partial f}{\partial y}) \tag {1.1} ∇f(x,y)=(∂x∂f​,∂y∂f​)(1.1)

  依此類推， 3 3 3維的梯度有如下形式：

(1.2) ∇ f ( x , y , z ) = ( ∂ f ∂ x , ∂ f ∂ y , ∂ f ∂ z ) ∇f(x,y,z)=(\frac{\partial f}{\partial x},\frac{\partial f}{\partial y},\frac{\partial f}{\partial z}) \tag {1.2} ∇f(x,y,z)=(∂x∂f​,∂y∂f​,∂z∂f​)(1.2)

  有時也會采用如下形式來表示梯度：

(1.3) ∇ f = ∂ f ∂ x ⃗ ∇f=\frac{\partial f}{\partial \vec x} \tag {1.3} ∇f=∂x

∂f​(1.3)

  梯度通常用來近似計算函數值（實際上就是一維形式的推廣)：

(1.4) f ( x ⃗ + Δ x ⃗ ) ≈ f ( x ⃗ ) + ∇ f ( x ⃗ ) ⋅ Δ x ⃗ f(\vec x+\Delta \vec x)\approx f(\vec x)+∇f(\vec x)\cdot \Delta \vec x \tag {1.4} f(x

+Δx

)≈f(x

)+∇f(x

)⋅Δx

(1.4)

  同樣的，多個函數的梯度就構成了一個矩陣：

(1.5) ∇ F ⃗ = ∇ ( f , g , h ) = ( ∂ f ∂ x ∂ f ∂ y ∂ f ∂ z ∂ g ∂ x ∂ g ∂ y ∂ g ∂ z ∂ h ∂ x ∂ h ∂ y ∂ h ∂ z ) = ( ∇ f ∇ g ∇ h ) ∇\vec F=∇(f,g,h)=\left( \begin{matrix} \frac{\partial f}{\partial x} & \frac{\partial f}{\partial y} & \frac{\partial f}{\partial z} \\ \frac{\partial g}{\partial x} & \frac{\partial g}{\partial y} & \frac{\partial g}{\partial z} \\ \frac{\partial h}{\partial x} & \frac{\partial h}{\partial y} & \frac{\partial h}{\partial z} \\ \end{matrix} \right) =\left( \begin{matrix}∇f\\ ∇g\\ ∇h\\ \end{matrix} \right) \tag {1.5} ∇F

=∇(f,g,h)=⎝⎜⎛​∂x∂f​∂x∂g​∂x∂h​​∂y∂f​∂y∂g​∂y∂h​​∂z∂f​∂z∂g​∂z∂h​​⎠⎟⎞​=⎝⎛​∇f∇g∇h​⎠⎞​(1.5)

2、散度（Divergence）

  散度算子僅僅應用于向量場，它衡量在某一點出相應的矢量聚集或者發散程度，測量方向為徑向，結果為标量。 2 2 2維、 3 3 3維形式的散度算子如下所示：

∇ ⋅ u ⃗ = ∇ ⋅ ( u , v ) = ∂ u ∂ x + ∂ v ∂ y ∇\cdot \vec u=∇\cdot (u,v)=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y} ∇⋅u

=∇⋅(u,v)=∂x∂u​+∂y∂v​

(1.6) ∇ ⋅ u ⃗ = ∇ ⋅ ( u , v , w ) = ∂ u ∂ x + ∂ v ∂ y + ∂ w ∂ z ∇\cdot \vec u=∇\cdot (u,v,w)=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z} \tag {1.6} ∇⋅u

=∇⋅(u,v,w)=∂x∂u​+∂y∂v​+∂z∂w​(1.6)

  輸入是矢量，而輸出為标量。類比梯度，散度符号 ∇ ⋅ u ⃗ ∇\cdot \vec u ∇⋅u

可以了解為梯度 ∇ ∇ ∇與矢量 u ⃗ \vec u u

的點乘：

(1.7) ∇ ⋅ u ⃗ = ( ∂ ∂ x , ∂ ∂ y , ∂ ∂ z ) ⋅ ( u , v , w ) = ∂ ∂ x u + ∂ ∂ y v + ∂ ∂ z w ∇\cdot \vec u=(\frac{\partial}{\partial x},\frac{\partial}{\partial y},\frac{\partial}{\partial z})\cdot (u,v,w)=\frac{\partial}{\partial x}u+\frac{\partial}{\partial y}v+\frac{\partial}{\partial z}w \tag {1.7} ∇⋅u

=(∂x∂​,∂y∂​,∂z∂​)⋅(u,v,w)=∂x∂​u+∂y∂​v+∂z∂​w(1.7)

  若矢量場散度為 0 0 0，則稱該矢量場無散度。

3、旋度（Curl）

  旋度衡量圍繞某一點的旋轉速度，測量方向為切向，三維形式的旋度是一個向量：

(1.8) ∇ × u ⃗ = ∇ × ( u , v , w ) = ( ∂ w ∂ y − ∂ v ∂ z , ∂ u ∂ z − ∂ w ∂ x , ∂ v ∂ x − ∂ u ∂ y ) ∇\times \vec u=∇\times (u,v,w) =(\frac{\partial w}{\partial y}-\frac{\partial v}{\partial z}, \frac{\partial u}{\partial z}-\frac{\partial w}{\partial x}, \frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}) \tag {1.8} ∇×u

=∇×(u,v,w)=(∂y∂w​−∂z∂v​,∂z∂u​−∂x∂w​,∂x∂v​−∂y∂u​)(1.8)

  倒推到 2 2 2維，我們取上式中的 w = 0 w=0 w=0，即矢量場為 ( u , v , 0 ) (u,v,0) (u,v,0)， 2 2 2維向量場的旋度是一個标量：

(1.9) ∇ × u ⃗ = ∇ × ( u , v ) = ∂ v ∂ x − ∂ u ∂ y ∇\times \vec u=∇\times (u,v)=\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y} \tag {1.9} ∇×u

=∇×(u,v)=∂x∂v​−∂y∂u​(1.9)

  同樣地，旋度符号 ∇ × u ⃗ ∇\times \vec u ∇×u

我們可以了解為梯度 ∇ ∇ ∇與矢量場 u ⃗ \vec u u

的叉乘：

(1.10) ∇ × u ⃗ = ( ∂ ∂ x , ∂ ∂ y , ∂ ∂ z ) × ( u , v , w ) ∇\times \vec u= (\frac{\partial }{\partial x}, \frac{\partial }{\partial y}, \frac{\partial }{\partial z})\times(u,v,w) \tag {1.10} ∇×u

=(∂x∂​,∂y∂​,∂z∂​)×(u,v,w)(1.10)

  若矢量場旋度為 0 0 0，則稱該矢量場無旋度。

4、拉普拉斯算子（Laplacian）

  拉普拉斯算子定義為梯度的散度，符号表示為 ∇ ⋅ ∇ ∇\cdot∇ ∇⋅∇，顯然 ∇ ⋅ ∇\cdot ∇⋅是散度，而後面的 ∇ ∇ ∇則為梯度，故拉普拉斯算子即梯度的散度，是一個二階微分算子。 2 2 2維、 3 3 3維形式分别如下：

∇ ⋅ ∇ f = ∂ 2 f ∂ x 2 + ∂ 2 f ∂ y 2 ∇\cdot∇f=\frac{\partial^2f}{\partial x^2}+\frac{\partial^2f}{\partial y^2} ∇⋅∇f=∂x2∂2f​+∂y2∂2f​

(1.11) ∇ ⋅ ∇ f = ∂ 2 f ∂ x 2 + ∂ 2 f ∂ y 2 + ∂ 2 f ∂ z 2 ∇\cdot∇f=\frac{\partial^2f}{\partial x^2}+\frac{\partial^2f}{\partial y^2}+\frac{\partial^2f}{\partial z^2} \tag {1.11} ∇⋅∇f=∂x2∂2f​+∂y2∂2f​+∂z2∂2f​(1.11)

  簡言之，拉普拉斯算子定義如下：

(1.12) ∇ ⋅ ∇ f = Σ i = 1 n ∂ 2 f ∂ x i 2 ∇\cdot∇f=\Sigma_{i=1}^n\frac{\partial^2f}{\partial x_i^2} \tag {1.12} ∇⋅∇f=Σi=1n​∂xi2​∂2f​(1.12)

  偏微分方程 ∇ ⋅ ∇ f = 0 ∇\cdot ∇f=0 ∇⋅∇f=0被稱為拉普拉斯方程；而如果右邊為某個非 0 0 0常數，即 ∇ ⋅ ∇ f = q ∇\cdot ∇f=q ∇⋅∇f=q，我們稱之為泊松方程。更一般地，如果梯度再乘上一個标量 a a a（如 1 / ρ 1/\rho 1/ρ)，即 ∇ ⋅ ( a ∇ f ) = q ∇\cdot (a∇f)=q ∇⋅(a∇f)=q，我們依舊稱之為泊松問題。

二、 N a i v e r − S t o k e s Naiver-Stokes Naiver−Stokes偏微分方程組

  流體模拟器的建構主要圍繞著名的不可壓縮 N a v i e r − S t o k e s Navier-Stokes Navier−Stokes方程展開，它是一個流體力學領域的偏微分方程，方程形式如下：

(2.1) ∂ u ⃗ ∂ t + u ⃗ ⋅ ∇ u ⃗ + 1 ρ ∇ p = g ⃗ + ν ∇ ⋅ ∇ u ⃗ \frac{\partial \vec u}{\partial t}+\vec u\cdot ∇\vec u+\frac1\rho∇p=\vec g+\nu∇\cdot∇\vec u \tag {2.1} ∂t∂u

​+u

⋅∇u

+ρ1​∇p=g

​+ν∇⋅∇u

(2.1)

(2.2) ∇ ⋅ u ⃗ = 0 ∇\cdot\vec u=0 \tag {2.2} ∇⋅u

=0(2.2)

  這個方程組看起非常地複雜，接下來我們就把它剖析成一個個比較簡單的部分來了解。

1、符号标記

  我們有必要定義一些實體量的符号用以标記：

  符号 u ⃗ \vec u u

在流體力學中通常表示為流體的速度矢量，記 3 3 3維的速度矢量 u ⃗ = ( u , v , w ) \vec u=(u,v,w) u

=(u,v,w)；

  希臘字元 ρ \rho ρ是流體的密度，對于水，該值大約為 1000 k g / m 3 1000kg/m^3 1000kg/m3，而空氣則大約為 1.3 k g / m 3 1.3kg/m^3 1.3kg/m3；

  字元 p p p代表壓力，流體對任何物體施加的機關面積力；

  字元 g ⃗ \vec g g

​則是我們熟悉的重力加速度，通常取 ( 0 , − 9.81 , 0 ) m / s 2 (0,-9.81,0)m/s^2 (0,−9.81,0)m/s2。我們約定 y y y軸向上，而 x x x軸和 z z z軸在水準面上。另外補充一點，我們把其他的一些類似的力都累加到 g ⃗ \vec g g

​上，也就是我們統一用 g ⃗ \vec g g

​表示所有類似力之和，這類力我們稱之為體積力（稱之為體積力是因為它們的力是作用到整個流體而不隻是流體的表面）；

  希臘字元 ν \nu ν是流體的運動粘度，它衡量流體的黏滞性。糖蜜之類的流體有非常高的粘度，而像酒精之類的流體有很低的粘度；

  其它一些矢量微積分的符号算子前面已經提到過，不再贅述。

2、動量方程

  偏微分方程 ( 2.1 ) (2.1) (2.1)我們稱之為動量方程，它本質上就是我們熟悉的牛頓定律 F ⃗ = m a ⃗ \vec F=m\vec a F

=ma

的形式，描述了施加在流體上的力是如何影響流體的運動。

  假設我們用粒子系統來模拟流體，每個粒子代表流體的一小滴，每個粒子有各自的品質 m m m、體積 V V V和速度 u ⃗ \vec u u

。為了讓整個粒子系統運作起來，我們必須弄清楚每個粒子所承受的力的作用。牛頓第二定律告訴我們： F ⃗ = m a ⃗ \vec F=m\vec a F

=ma

，而根據加速度定義，我們有：

(2.3) a ⃗ = D u ⃗ D t \vec a=\frac{D\vec u}{Dt} \tag {2.3} a

=DtDu

​(2.3)

  符号 D D D是指物質導數，所謂物質導數，就是對流體質點求導數，而且是針對流體質點（在這裡就是流體粒子）而不是空間的固定點。因而牛頓第二定律就變成：

(2.4) m D u ⃗ D t = F ⃗ m\frac{D\vec u}{Dt}=\vec F \tag {2.4} mDtDu

​=F

(2.4)

  那麼流體粒子承受的力有哪些呢？一個最簡單的力就是重力： m g ⃗ m\vec g mg

​。而其他的流體質點（或其他流體粒子）也會對目前的流體粒子産生作用力。流體内部的互相作用力之一便是壓力，較大壓力的區域會向較低壓力區域産生作用力。值得注意的是，我們隻關注施加在粒子上的壓力的淨合力。例如，若施加在粒子上壓力在每個方向上都相等，那麼它的壓力的合力便為0。我們用壓力的負梯度（取負是因為方向是由壓力大的區域指向壓力小的區域）來衡量在目前流體粒子處壓力的不平衡性，即取 − ∇ p -∇p −∇p。那麼流體粒子所承受的壓力就是對 − ∇ p -∇p −∇p在整個流體粒子的體積上進行積分，為了簡化，我們簡單地将 V V V與 − ∇ p -∇p −∇p相乘，故粒子壓力部分為 − V ∇ p -V∇p −V∇p。

  其他的流體互相作用力則是由流體的黏性産生的，我們直覺地把這種力了解為盡可能使得粒子以周圍區域的平均速度移動的力，也就是使得粒子的速度與周圍區域粒子速度的差距最小化。拉普拉斯算子是衡量一個量與之周圍區域該量平均值之差的算符，因而 ∇ ⋅ ∇ u ⃗ ∇\cdot∇\vec u ∇⋅∇u

是目前粒子速度矢量與周圍區域平均速度矢量之差。為了計算粘滞力，我們同樣對 ∇ ⋅ ∇ u ⃗ ∇\cdot∇\vec u ∇⋅∇u

在整個粒子體積 V V V上進行積分，與前面類似，我們簡單取 V ∇ ⋅ ∇ u ⃗ V∇\cdot∇\vec u V∇⋅∇u

。除此之外，我們還引進一個稱為動力粘度系數的實體量，符号為 μ \mu μ。因而粘滞力為 V μ ∇ ⋅ ∇ u ⃗ V\mu∇\cdot∇\vec u Vμ∇⋅∇u

。

  把重力、壓力和粘滞力綜合一起，我們可得：

(2.5) m D u ⃗ D t = F ⃗ = m g ⃗ − V ∇ p + V μ ∇ ⋅ ∇ u ⃗ m\frac{D\vec u}{Dt}=\vec F=m\vec g-V∇p+V\mu∇\cdot∇\vec u \tag {2.5} mDtDu

​=F

=mg

​−V∇p+Vμ∇⋅∇u

(2.5)

  當粒子系統中的粒子數量趨于無窮大，而每個粒子大小趨于 0 0 0時，會産生一個問題：此時每個粒子的品質 m m m和體積 V V V變為 0 0 0，此時上式變得沒有意義。為此，我們把 ( 2.5 ) (2.5) (2.5)式調整一下，兩邊同除以體積 V V V，又因 ρ = m / V \rho=m/V ρ=m/V，故有：

(2.6) ρ D u ⃗ D t = ρ g ⃗ − ∇ p + μ ∇ ⋅ ∇ u ⃗ \rho\frac{D\vec u}{Dt}=\rho\vec g-∇p+\mu∇\cdot∇\vec u \tag {2.6} ρDtDu

​=ρg

​−∇p+μ∇⋅∇u

(2.6)

  兩邊同除以 ρ \rho ρ，移項調整：

(2.7) D u ⃗ D t + 1 ρ ∇ p = g ⃗ + μ ρ ∇ ⋅ ∇ u ⃗ \frac{D\vec u}{Dt}+\frac1\rho∇p=\vec g+\frac\mu\rho∇\cdot∇\vec u \tag {2.7} DtDu

​+ρ1​∇p=g

​+ρμ​∇⋅∇u

(2.7)

  為了進一步簡化，定義運動粘度為 ν = μ / ρ \nu=\mu/\rho ν=μ/ρ，式 ( 2.7 ) (2.7) (2.7)變為：

(2.8) D u ⃗ D t + 1 ρ ∇ p = g ⃗ + ν ∇ ⋅ ∇ u ⃗ \frac{D\vec u}{Dt}+\frac1\rho∇p=\vec g+\nu∇\cdot∇\vec u \tag {2.8} DtDu

​+ρ1​∇p=g

​+ν∇⋅∇u

(2.8)

  我們已經快把動量方程推導出來，現在我們要把物質導數 D u ⃗ D t \frac{D\vec u}{Dt} DtDu

​弄清楚，為此，我們需要了解兩種描述方法：拉格朗日描述和歐拉描述。

3、拉格朗日描述與歐拉描述

  當我們嘗試研究流體或可變形固體的運動的時候，通常有兩種方法來描述：拉格朗日描述（ Lagrangian viewpoint）、歐拉描述（Eulerian viewpoint）。

  拉格朗日描述方法是我們比較熟悉的方法，這種描述方法把物體看成是由類似于粒子系統的形式組成，固體或流體的每個點看作一個獨立的粒子，粒子有各自相應的位置 x ⃗ \vec x x

和速度 u ⃗ \vec u u

。我們可以把粒子了解為組成物體的分子。對于我們通常采用拉格朗日描述法進行模組化模拟，即用一系列離散的粒子集來建構，粒子之間通過網格相聯系。

  歐拉描述方法則采用了完全不同的角度，它常被用于流體力學中。與拉格朗日描述追蹤每個物體粒子的方法不同，歐拉描述關注點是空間中的一個固定點，并考察在這個固定點上流體性質（如密度、速度、溫度等）是如何随着時間變化的。流體流動經過這個固定點可能會導緻這個固定點的實體性質發生一些變化（如一個溫度較高的流體粒子流經這個固定點，後面緊跟着一個溫度較低的流體粒子流過固定點，那麼這個固定點的溫度會降低，但是并沒有任何一個流體粒子的溫度發生了變化）。

  用天氣測量舉個簡單的例子：拉格朗日描述方法就是你乘坐在一個随風而飄的熱氣球上，測量周圍空氣的壓力、密度和渾濁度等天氣名額；而歐拉描述方法就是你固定在地面上，測量流過的空氣的天氣名額。

  歐拉描述法似乎看起來帶來了一些不必要的複雜度，但是目前大多數的流體模拟器都是基于歐拉描述法，這是因為歐拉描述法相比于拉格朗日描述法有一些不可比拟的優點：歐拉描述法能夠更加友善地計算一些實體量的空間導數（例如壓力梯度和粘度）；而如果用粒子方法的話（即拉格朗日描述法），那麼計算實體量相對于空間位置的變化是比較難的。

  把拉格朗日描述法和歐拉描述法聯系起來的關鍵點就是物質導數。首先從拉格朗日描述法出發，假設有一群粒子，每個粒子都有各自的位置 x ⃗ \vec x x

和速度 u ⃗ \vec u u

。記 q q q為通用的實體量（如密度、速度和溫度等），每個粒子有其對應的 q q q值。方程 q ( t , x ⃗ ) q(t,\vec x) q(t,x

)描述在時間點 t t t而位置為 x ⃗ \vec x x

的粒子對應的實體量值 q q q。則一個粒子的實體量 q q q随時間 t t t的變化率是多少？這是一個拉格朗日描述角度下的問題，我們取對時間 t t t的導數（注意用到了求導鍊式法則，以及 ∂ q ∂ x ⃗ = ∇ q \frac{\partial q}{\partial \vec x}=∇q ∂x

∂q​=∇q和 u ⃗ = d x ⃗ d t ） \vec u=\frac{d\vec x}{dt}） u

=dtdx

​）：

(2.9) d d t q ( t , x ⃗ ) = ∂ q ∂ t + ∇ q ⋅ d x ⃗ d t = ∂ q ∂ t + ∇ q ⋅ u ⃗ ≡ D q D t \frac d{dt}q(t,\vec x)=\frac{\partial q}{\partial t}+∇q\cdot\frac{d\vec x}{dt}=\frac{\partial q}{\partial t}+∇q\cdot\vec u\equiv\frac{Dq}{Dt} \tag {2.9} dtd​q(t,x

)=∂t∂q​+∇q⋅dtdx

​=∂t∂q​+∇q⋅u

≡DtDq​(2.9)

  這就是物質導數。把式 ( 2.9 ) (2.9) (2.9)代入式 ( 2.8 ) (2.8) (2.8)我們就得到了流體動量方程 ( 2.1 ) (2.1) (2.1)。物質導數針對的是流體質點（在這裡就是流體粒子）而不是空間的固定點。式 ( 2.9 ) (2.9) (2.9)寫完整一點就是：

(2.10) D q D t = ∂ q ∂ t + u ∂ q ∂ x + v ∂ q ∂ y + w ∂ q ∂ z \frac{Dq}{Dt}=\frac{\partial q}{\partial t}+u\frac{\partial q}{\partial x}+v\frac{\partial q}{\partial y}+w\frac{\partial q}{\partial z} \tag {2.10} DtDq​=∂t∂q​+u∂x∂q​+v∂y∂q​+w∂z∂q​(2.10)

  對于給定的速度場 u ⃗ \vec u u

， 流體的實體性質如何在這個速度場 u ⃗ \vec u u

下變化的計算我們稱之為對流（advection）。一個最簡單的對流方程，就是其實體量的物質導數為 0 0 0，如下所示：

(2.11) D q D t = 0    ⟹    ∂ q ∂ t + u ⃗ ⋅ ∇ q = 0 \frac{Dq}{Dt}=0\implies\frac{\partial q}{\partial t}+\vec u\cdot ∇q = 0 \tag {2.11} DtDq​=0⟹∂t∂q​+u

⋅∇q=0(2.11)

  公式 ( 2.11 ) (2.11) (2.11)的意義即在拉格朗日視角觀察下，每個流體粒子的實體量保持不變。

4、不可壓縮性

  關于流體的壓縮性在此不做過多的實體細節描述，隻需知道一點：通常情況下流體的體積變化非常小（除開一些極端的情況，而且這些極端情況我們日常生活中較少出現）。可壓縮流體的模拟涉及到非常複雜的情況，往往需要昂貴的計算資源開銷，為此在計算機流體模拟中我們通常把所有的流體當作是不可壓縮的，即它們的體積不會發生變化。

  任取流體的一部分，設其體積為 Ω \Omega Ω而其邊界閉合曲面為 ∂ Ω \partial\Omega ∂Ω，我們可以通過圍繞邊界曲面 ∂ Ω \partial\Omega ∂Ω對流體速度 u ⃗ \vec u u

在曲面法線方向上的分量進行積分來衡量這塊部分流體的體積變化速率：

(2.12) d d t V o l u m e ( Ω ) = ∫ ∫ ∂ Ω u ⃗ ⋅ n \frac d{dt}Volume(\Omega)=\int\int_{\partial\Omega}\vec u\cdot n \tag{2.12} dtd​Volume(Ω)=∫∫∂Ω​u

⋅n(2.12)

  對于不可壓縮的流體，其體積保持為某個常量，故其體積變化速率為 0 0 0：

(2.13) ∫ ∫ ∂ Ω u ⃗ ⋅ n = 0 \int\int_{\partial\Omega}\vec u\cdot n=0 \tag {2.13} ∫∫∂Ω​u

⋅n=0(2.13)

  由高斯散度定理，我們可以把式 ( 2.13 ) (2.13) (2.13)轉換為體積分：

(2.14) ∫ ∫ ∂ Ω u ⃗ ⋅ n = ∫ ∫ ∫ Ω ∇ ⋅ u ⃗ = 0 \int\int_{\partial\Omega}\vec u\cdot n=\int\int\int_\Omega∇\cdot \vec u=0 \tag{2.14} ∫∫∂Ω​u

⋅n=∫∫∫Ω​∇⋅u

=0(2.14)

  式 ( 13 ) (13) (13)應該對任意的 Ω \Omega Ω成立，意即無論 Ω \Omega Ω取何值，積分值均為 0 0 0。這種情況下隻有令積分函數值取 0 0 0方可成立，即對 0 0 0積分無論 Ω \Omega Ω取何值結果均為 0 0 0。是以有：

(2.15) ∇ ⋅ u ⃗ = 0 ∇\cdot \vec u=0 \tag{2.15} ∇⋅u

=0(2.15)

  這就是 N a v i e r − S t o k e s Navier-Stokes Navier−Stokes方程中的不可壓縮條件 ( 2.2 ) (2.2) (2.2)。滿足不可壓縮條件的速度場被稱為是無散度的，即在該速度場下流體體積既不膨脹也不坍縮，而是保持在一個常量。模拟不可壓縮流體的關鍵部分就是使得流體的速度場保持無散度的狀态，這也是流體内部壓力的來源。

  為了把壓力與速度場的散度聯系起來，我們在動量方程 ( 2.1 ) (2.1) (2.1)兩邊同時取散度：

(2.16) ∇ ⋅ ∂ u ⃗ ∂ t + ∇ ⋅ ( u ⃗ ⋅ ∇ u ⃗ ) + ∇ ⋅ 1 ρ ∇ p = ∇ ⋅ ( g ⃗ + ν ∇ ⋅ ∇ u ⃗ ) ∇\cdot\frac{\partial \vec u}{\partial t}+∇\cdot(\vec u\cdot ∇\vec u)+∇\cdot\frac1\rho∇p=∇\cdot(\vec g+\nu∇\cdot∇\vec u) \tag {2.16} ∇⋅∂t∂u

​+∇⋅(u

⋅∇u

)+∇⋅ρ1​∇p=∇⋅(g

​+ν∇⋅∇u

)(2.16)

  對于上式 ( 2.16 ) (2.16) (2.16)第一項，我們轉變一下求導次序：

(2.17) ∂ ∂ t ∇ ⋅ u ⃗ \frac {\partial}{\partial t}∇\cdot\vec u \tag {2.17} ∂t∂​∇⋅u

(2.17)

  如果滿足流體不可壓縮條件，那麼式 ( 2.17 ) (2.17) (2.17)取值 0 0 0（因為無散度），然後我們調整一下式 ( 2.16 ) (2.16) (2.16)可得關于壓力的方程：

(2.18) ∇ ⋅ 1 ρ ∇ p = ∇ ⋅ ( − u ⃗ ⋅ ∇ u ⃗ + g ⃗ + ν ∇ ⋅ ∇ u ⃗ ) ∇\cdot\frac1\rho∇p=∇\cdot(-\vec u\cdot ∇\vec u+\vec g+\nu∇\cdot∇\vec u) \tag{2.18} ∇⋅ρ1​∇p=∇⋅(−u

⋅∇u

+g

​+ν∇⋅∇u

)(2.18)

5、丢棄粘度項

  在某些流體如蜂蜜、小水珠等的模拟中，粘滞力起着非常重要的作用。但是在大多數流體動畫模拟中，粘滞力的影響微乎其微，為此秉持着方程組越簡單越好的原則，我們常常丢棄粘度項。當然這也不可避免地帶來一些誤差，事實上，在計算流體力學中盡可能地減少丢棄粘度項帶來的誤差是一個非常大的挑戰。下面的叙述都是基于丢棄粘度項的前提。

  丢棄了粘度項的 N a v i e r − S t o k e s Navier-Stokes Navier−Stokes方程被稱為歐拉方程，而這種理想的流體則是無粘度的。丢棄了粘度項的歐拉方程如下：

(2.19) D u ⃗ D t + 1 ρ ∇ p = g ⃗ \frac{D\vec u}{Dt}+\frac1\rho∇p=\vec g \tag {2.19} DtDu

​+ρ1​∇p=g

​(2.19)

(2.20) ∇ ⋅ u ⃗ = 0 ∇\cdot\vec u=0 \tag{2.20} ∇⋅u

=0(2.20)

  大多數的流體模拟的計算方程都是歐拉方程。

6、邊界條件

  目前為止我們讨論的都是流體内部的情況，然而邊界部分也是流體模拟非常關鍵的部分。在流體模拟中我們僅僅關注兩種邊界條件：固體牆（solid walls）、自由面（free surfaces）。

  固體牆顧名思義就是流體與固體接觸的邊界，用速度來描述很簡單：流體既不會流進固體内部也不會從固體内部流出，是以流體在固體牆法線方向上的分量為 0 0 0：

(2.21) u ⃗ ⋅ n = 0 \vec u\cdot n=0 \tag {2.21} u

⋅n=0(2.21)

  當然，上述是固體自身不移動的情況下。通常來說，流體速度在法線方向上的分量與固體的移動速度在法線方向上的分量應該保持一緻：

(2.22) u ⃗ ⋅ n = u ⃗ s o l i d ⋅ n \vec u\cdot n=\vec u_{solid}\cdot n \tag{2.22} u

⋅n=u

solid​⋅n(2.22)

  上述的兩個公式都是僅對流體速度在法線方向上的分量做了限制，對于無粘度的流體，切線方向上的流體速度與固體的移動速度無必然的聯系。

  自由面是另外一個非常重要的邊界條件，它通常就是與另外一種流體相接壤的邊界部分。例如在模拟水花四濺時，水流表面不與固體接觸的都是自由面（如與空氣這種流體接觸）。因空氣密度遠小于水導緻空氣對水體的仿真影響非常小，為了簡化模拟，我們将空氣所占的空間設為某個固定大氣壓的區域，設為 0 0 0是最友善的方案，此時自由面就是壓強 p = 0 p=0 p=0的水體表面。

  在小規模的流體仿真中，自由面的表面張力占據着非常重要的地位。在微觀分子層面下，表面張力的存在是因為不同的分子互相吸引産生的力。從幾何的角度來解釋就是，表面張力就是促使流體的表面積盡可能小的一種力。實體學上，兩種不同的流體之間實際上存在着與表面平均曲率成正比的壓力驟變：

(2.23) [ p ] = λ k . [p]=\lambda k. \tag {2.23} [p]=λk.(2.23)

  公式 ( 2.23 ) (2.23) (2.23)中的 [ p ] [p] [p]記為壓力之差。 λ \lambda λ是表面張力系數，可以根據模拟的流體類型查找對應的張力系數（例如空氣與水在室溫下張力系數為 λ ≈ 0.073 N / m \lambda \approx 0.073N/m λ≈0.073N/m）。而 k k k就是平均曲率，機關為 m − 1 m^{-1} m−1。又因為我們常常設空氣的壓力為 0 0 0，是以水與空氣交界的自由面的壓力為：

(2.24) p = λ k p=\lambda k \tag {2.24} p=λk(2.24)

​

三、N-S方程的分步求解

  有了對以上對 N a v i e r − S t o k e s Navier-Stokes Navier−Stokes方程的理論支撐，接下來我們就要如何用計算機來對該組偏微分方程進行離散化求解。為了程式的松耦合性以及使計算盡可能地高效、簡單，在流體模型領域，我們将流體方程分成幾個獨立的步驟，然後按順序先後推進。對于不可壓縮的無粘度流體方程（即前面的歐拉方程 ( 2.19 ) (2.19) (2.19)和 ( 2.20 ) (2.20) (2.20)，我們将其離散化成對流項（advection）如公式 ( 3.1 ) (3.1) (3.1)、體積力項（body force）如公式 ( 3.2 ) (3.2) (3.2)、壓力/不可壓縮項如公式 ( 3.3 ) (3.3) (3.3)：

(3.1) D q D t = 0 \frac{Dq}{Dt}=0 \tag {3.1} DtDq​=0(3.1)

(3.2) ∂ u ⃗ ∂ t = g ⃗ \frac{\partial \vec u}{\partial t}=\vec g \tag {3.2} ∂t∂u

​=g

​(3.2)

(3.3) { ∂ u ⃗ ∂ t + 1 ρ ∇ p = 0 ∇ ⋅ u ⃗ = 0 \begin{cases} \frac{\partial \vec u}{\partial t}+\frac{1}{\rho}∇p=0\\ ∇\cdot\vec u=0 \end{cases} \tag {3.3} {∂t∂u

​+ρ1​∇p=0∇⋅u

=0​(3.3)

  需要注意的是，在對流項公式 ( 3.1 ) (3.1) (3.1)中我們用了一個通用量的符号 q q q是因為我們不僅僅要對流體的速度進行對流，還需要對其他實體量進行對流。我們記對流項公式 ( 3.1 ) (3.1) (3.1)的對流計算算法為 a d v e c t ( u ⃗ , Δ t , q ) advect(\vec u, \Delta t, q) advect(u

,Δt,q)，即對于給定的時間步長 Δ t \Delta t Δt和速度場 u ⃗ \vec u u

，對實體量q進行對流。

  對于體積力項 ( 3.2 ) (3.2) (3.2)，我們采用簡單的前向歐拉法即可： u ⃗ ← u ⃗ + g Δ t \vec u \leftarrow \vec u + g\Delta t u

←u

+gΔt。

  對于壓力/不可壓縮項 ( 3.3 ) (3.3) (3.3)，我們用一個稱為 p r o j e c t ( Δ t , u ⃗ ) project(\Delta t, \vec u) project(Δt,u

)的算法，通過 p r o j e c t ( Δ t , u ⃗ ) project(\Delta t, \vec u) project(Δt,u

)計算出正确的壓力以確定速度場 u ⃗ \vec u u

的無散度性質。歐拉方案不會着重研究具體粒子間的作用力，因而不會正向去求解 1 ρ ∇ p \frac{1}{\rho}∇p ρ1​∇p，它是利用流體不可壓縮的特性，将速度場 u ⃗ \vec u u

投影到散度為 0 0 0的空間上，間接地解算了壓力項。這種思想相當于，已知一個中間量 u ⃗ t e m p \vec u_{temp} u

temp​，對這個中間量的唯一一個操作（如正向求解壓力 1 ρ ∇ p \frac{1}{\rho}∇p ρ1​∇p）不可行，但是直到最終量 u ⃗ f i a n l \vec u_{fianl} u

fianl​符号的一個性質（散度為 0 0 0），于是隻要将 u ⃗ t e m p \vec u_{temp} u

temp​投影到符合散度為 0 0 0的特性平面上，即可間接地還原正向求解壓力的操作，得到最終的速度場 u ⃗ t e m p \vec u_{temp} u

temp​。

  對流項 a d v e c t ( u ⃗ , Δ t , q ) advect(\vec u, \Delta t, q) advect(u

,Δt,q)的輸入速度場 u ⃗ \vec u u

要確定為無散度的狀态，投影項 p r o j e c t ( Δ t , u ⃗ ) project(\Delta t, \vec u) project(Δt,u

)確定了流體體積保持不變，因而投影項輸出的速度場必然是無散度的。是以我們隻要確定投影項 p r o j e c t ( Δ t , u ⃗ ) project(\Delta t, \vec u) project(Δt,u

)輸出的速度場 u ⃗ \vec u u

作為對流項 a d v e c t ( u ⃗ , Δ t , q ) advect(\vec u, \Delta t, q) advect(u

,Δt,q)的輸入即可，這時我們的分步求解流體方程的優勢就展現出來了，其僞代碼如下所示。

算法1 Fluid Simulation( u ⃗ n \vec u_n u

n​, Δ t \Delta t Δt):

1: 初始化速度場 u ⃗ n \vec u_n u

n​,使得 u ⃗ n \vec u_n u

n​無散度

2: 對于每個時間步 n = 0 , 1 , 2 , . . . n = 0,1,2,... n=0,1,2,...

3:   決定一個合理的時間步長 Δ t = t n + 1 − t n \Delta t = t_{n+1}-t_n Δt=tn+1​−tn​

4:   對流項計算 u ⃗ A = a d v e c t ( u ⃗ n , Δ t , q ⃗ ) \vec u_A=advect(\vec u_n,\Delta t,\vec q) u

A​=advect(u

n​,Δt,q

​)

5:   體積力項計算 u ⃗ B = u ⃗ A + Δ t g ⃗ \vec u_B=\vec u_A+\Delta t\vec g u

B​=u

A​+Δtg

​

6:   無散度投影 u ⃗ n + 1 = p r o j e c t ( Δ t , u ⃗ B ) \vec u_{n+1}=project(\Delta t,\vec u_B) u

n+1​=project(Δt,u

B​)

1、時間步長

  在流體模拟算法中，确定适當的時間步長是算法的第一步。因為計算流體模拟的最後結果是呈現在螢幕上的，是以 Δ t \Delta t Δt的選取與螢幕的重新整理率有重要的關系。若選取的 Δ t \Delta t Δt有 t n + Δ t > t f r a m e t_n+\Delta t > t_{frame} tn​+Δt>tframe​，那麼必須做一個截斷使 Δ t = t f r a m e − t n \Delta t=t_{frame}-t_n Δt=tframe​−tn​。此外，流體模拟的三個步驟即對流項、體積力項、無散度投影項對時間步長 Δ t \Delta t Δt的要求不盡相同，要選擇一個滿足所有要求的最小時間步長能確定計算的收斂性。此外，一方面為了流體模拟的真實性，我們可能需要選取一個足夠小的時間步長來複現流體的高品質細節。另一方面，有時高性能的需求又使得我們不能選取太小的時間步長去渲染一幀。假設一幀至少要進行三個時間步的模拟，那麼 Δ t \Delta t Δt應該至少設成幀間隔時間的三分之一。

2、網格結構

  歐拉法的整個流程都是基于網格的，是以合理的網格結構是算法高效的關鍵點。 H a r l o w Harlow Harlow和 W e l c h Welch Welch提出了一種經典的 M A C MAC MAC（marker and cell）網格結構，許多不可壓縮流體模拟的算法都在這個網格結構上呈現出了良好的效率。 M A C MAC MAC網格是一種交叉排列的網格，不同類型的實體量被存儲于網格的不同位置。以二維的網格為例，如圖3-1左圖所示，流體粒子的壓力資料存儲于網格的中心點 P i , j P_{i,j} Pi,j​，而速度則沿着笛卡爾坐标被分成了兩部分。水準方向的 u u u成分被存儲在了網格單元豎直邊的中心處，例如網格單元 ( i , j ) (i,j) (i,j)和 ( i + 1 , j ) (i+1,j) (i+1,j)之間的水準速度記為 u i + 1 / 2 , j u_{i+1/2,j} ui+1/2,j​。垂直方向的 v v v成分則被存儲在了網格單元水準面的中心上。這樣的存儲方案十分有利于估算流體流進/流出某個網格單元的量。

圖3-1 MAC網格,左圖二維,右圖三維   擴充到三維的情況，$MAC$網格同樣是交錯排列的結構網格，如圖3-1右圖所示。壓力數值存儲在立方體網格單元的中心，三個速度分量分别被記錄在立方體網格單元的三個表面的中心點上。在數值計算時，這樣的配置設定方式使得我們可以準确地采用中心差分法計算壓力梯度和速度的散度，同時克服了中心差分法的一個普遍的缺點。一維的情況為例，在網格頂點位置$...,q_{i-1},q_i,q_{i+1}...$上估算量場$q$的導數，為了無偏（所謂無偏，就是不偏向左邊或者右邊）估計網格頂點$i$處的$\frac{\partial q}{\partial x}$，一種比較自然的方式就是采用一階中心差分法： $$ (\frac{\partial q}{\partial x})_i\approx \frac{q_{i+1}-q_{i-1}}{2\Delta x} \tag {3.4} $$   公式$(3.4)$是無偏的，且精确度為$O(\Delta x^2)$。而前向歐拉差分法偏向右邊且精确度隻有$O(\Delta x)$： $$ (\frac{\partial q}{\partial x})_i\approx \frac{q_{i+1}-q_i}{\Delta x} \tag {3.5} $$   然而，公式$(3.4)$存在着一個非常嚴重的問題：網格點$i$的估算導數完全忽略了$q_i$的值。數學上，隻有常數函數的一階導數為零。但是公式$(3.4)$遇到了鋸齒函數如$q_i=(-1)^i$時，它錯誤地将該類函數的導數估算為$0$，這種問題被稱為零空間問題（null-space problem）。

  交叉錯排的 M A C MAC MAC網格完美地克服了中心差分法的零空間問題，同時也保持了它的無偏二階精度。在 M A C MAC MAC網格上運用中心差分法，網格點 i i i處的估算導數公式如下所示：

(3.6) ( ∂ q ∂ x ) i ≈ q i + 1 / 2 − q i − 1 / 2 Δ x (\frac{\partial q}{\partial x})_i\approx\frac{q_{i+1/2}-q_{i-1/2}}{\Delta x} \tag {3.6} (∂x∂q​)i​≈Δxqi+1/2​−qi−1/2​​(3.6)

   M A C MAC MAC網格确實給流體的壓力計算和不可壓縮性的處理帶來了很大的便利，但與此同時也帶來了一些其他方面的麻煩。如果我們要估算某個地方的速度向量，即便采樣點恰好在網格點上我們也要做一些插值才能擷取相應的速度向量。在網格點處，我們通常采用平均法，以二維為例：

(3.7) u ⃗ i , j = ( u i − 1 / 2 , j + u i + 1 / 2 , j 2 , v i , j − 1 / 2 + v i , j + 1 / 2 2 ) , u ⃗ i + 1 / 2 , j = ( u i + 1 / 2 , j , v i , j − 1 / 2 + v i , j + 1 / 2 + v i + 1 , j − 1 / 2 + v i + 1 , j + 1 / 2 4 ) , u ⃗ i , j + 1 / 2 = ( u i − 1 / 2 , j + u i + 1 / 2 , j + u i − 1 / 2 , j + 1 + u i + 1 / 2 , j + 1 4 , v i , j + 1 / 2 ) . \vec u_{i,j}=(\frac{u_{i-1/2,j}+u_{i+1/2,j}}{2},\frac{v_{i,j-1/2}+v_{i,j+1/2}}{2}),\\ \vec u_{i+1/2,j}=(u_{i+1/2,j},\frac{v_{i,j-1/2}+v_{i,j+1/2}+v_{i+1,j-1/2}+v_{i+1,j+1/2}}{4}),\\ \vec u_{i,j+1/2}=(\frac{u_{i-1/2,j}+u_{i+1/2,j}+u_{i-1/2,j+1}+u_{i+1/2,j+1}}{4},v_{i,j+1/2}).\tag {3.7} u

i,j​=(2ui−1/2,j​+ui+1/2,j​​,2vi,j−1/2​+vi,j+1/2​​),u

i+1/2,j​=(ui+1/2,j​,4vi,j−1/2​+vi,j+1/2​+vi+1,j−1/2​+vi+1,j+1/2​​),u

i,j+1/2​=(4ui−1/2,j​+ui+1/2,j​+ui−1/2,j+1​+ui+1/2,j+1​​,vi,j+1/2​).(3.7)

  最後，在實作中下标索引一般沒有浮點數之說，前面直接采用 i + 1 / 2 i+1/2 i+1/2的記法是為了便于叙述。一般約定如下：

(3.8) p ( i , j , k ) = p i , j , k , u ( i , j , k ) = u i − 1 / 2 , j , k , v ( i , j , k ) = v i , j − 1 / 2 , k , w ( i , j , k ) = w i , j , k − 1 / 2 . p(i,j,k)=p_{i,j,k},\\ u(i,j,k)=u_{i-1/2,j,k},\\ v(i,j,k)=v_{i,j-1/2,k},\\ w(i,j,k)=w_{i,j,k-1/2}. \tag{3.8} p(i,j,k)=pi,j,k​,u(i,j,k)=ui−1/2,j,k​,v(i,j,k)=vi,j−1/2,k​,w(i,j,k)=wi,j,k−1/2​.(3.8)

  因而對于 n x × n y × n z nx\times ny\times nz nx×ny×nz分辨率的網格，壓力數值存儲在 n x × n y × n z nx\times ny\times nz nx×ny×nz的數組中，速度的 u u u成分存儲在 ( n x + 1 ) × n y × n z (nx+1)\times ny\times nz (nx+1)×ny×nz數組中，速度的 v v v成分存儲在 n x × ( n y + 1 ) × n z nx\times (ny+1)\times nz nx×(ny+1)×nz數組中，速度的 w w w成分存儲在 n x × n y × ( n z + 1 ) nx\times ny\times (nz+1) nx×ny×(nz+1)數組中。

四、對流算法

  求解如下所示的對流方程是流體模拟的關鍵一步：

(4.1) D q D t = 0 \frac{Dq}{Dt}=0 \tag {4.1} DtDq​=0(4.1)

  我們把這個對流數值計算的算法記為：

(4.2) q n + 1 = a d v e c t ( u ⃗ , Δ t , q n ) q^{n+1}=advect(\vec u,\Delta t,q^n) \tag {4.2} qn+1=advect(u

,Δt,qn)(4.2)

  公式 ( 4.2 ) (4.2) (4.2)中的各個符号含義：

   u ⃗ \vec u u

：在 M A C MAC MAC網格上的離散化的速度場；

   Δ t \Delta t Δt：時間步長；

   q n q^n qn：目前的實體量場 q q q（如流體密度、速度、燃燒物濃度等）；

   q n + 1 q^{n+1} qn+1：經過對流後得到的新的量場。

  在這裡要特别注意，輸入對流算法的速度場 u ⃗ \vec u u

必須是無散度的，否則模拟結果會出現一些奇怪的失真現象。

1、半拉格朗日對流算法（Semi-Lagrangian Advection）

  一維情況下，對流方程 ( 4.1 ) (4.1) (4.1)寫成偏微分的形式如下：

(4.3) ∂ q ∂ t + u ∂ q ∂ x = 0 \frac{\partial q}{\partial t}+u\frac{\partial q}{\partial x}=0 \tag {4.3} ∂t∂q​+u∂x∂q​=0(4.3)

  分别采用前向歐拉差分法計算對時間的偏導和中心差分法計算對空間的偏導，我們有：

(4.4) q i n + 1 − q i n Δ t + u i n q i + 1 n − q i − 1 n 2 Δ x = 0 \frac{q^{n+1}_{i}-q^n_i}{\Delta t}+u^n_i\frac{q^n_{i+1}-q^n_{i-1}}{2\Delta x}=0 \tag {4.4} Δtqin+1​−qin​​+uin​2Δxqi+1n​−qi−1n​​=0(4.4)

  轉成以 q i n + 1 q^{n+1}_i qin+1​為計算目标的顯式公式，得：

(4.5) q i n + 1 = q i n − Δ t u i n q i + 1 n − q i − 1 n 2 Δ x q^{n+1}_i=q^n_i-\Delta t u^n_i\frac{q^n_{i+1}-q^n_{i-1}}{2\Delta x} \tag {4.5} qin+1​=qin​−Δtuin​2Δxqi+1n​−qi−1n​​(4.5)

  公式 ( 4.5 ) (4.5) (4.5)看起來沒什麼問題，但是卻存在非常嚴重的漏洞。首先，前向歐拉法被證明是無條件不穩定的空間離散方法：無論取多麼小Δ?，随着時間步的推進，累積誤差終将發散。即使使用更穩定的時間積分方法來取代前向歐拉方法，解決了時間上的PDE（Partial Differential Equation，偏微分方程）計算，空間上的PDE計算還是會帶來重大的麻煩。标準中心差分方法不可避免地會出現的零空間問題，具有高頻震蕩性質的速度場對空間的導數被錯誤地計算為 0 0 0或幾乎為 0 0 0，低離速度分量被分離出來，進而導緻模拟效果中出現許多奇怪的高頻擺動和震蕩。

  針對這些問題，研究者們提出了一個解然不同的、更加簡單和更具實體直覺意義的半拉格朗日法。之是以叫半拉格朗日法，是因為這種方法是以拉格朗日視角去解決歐拉視角的對流方程（“半”字的由來）。假設我們的目标是求解網格點 x ⃗ G \vec x_G x

G​的在第 n + 1 n+1 n+1個時間步時關于實體量 q q q的新值，記為 q G n + 1 q^{n+1}_G qGn+1​。在拉格朗日的視角下，我們可以尋找在第 n + 1 n+1 n+1時間步之前，是空間中的哪一個點上的流體粒子在速度場 u ⃗ \vec u u

的作用下“流向”了 x ⃗ G \vec x_G x

G​，我們記這個粒子在第 n n n個時間步時的網格位置為 x ⃗ P \vec x_P x

P​，則第 n + 1 n+1 n+1個時間步時 x ⃗ G \vec x_G x

G​的 q G n + 1 q^{n+1}_G qGn+1​即為第 n n n個時間步時 x ⃗ P \vec x_P x

P​的 q P n q^{n}_P qPn​。如下圖4-1為半拉格朗日對流法的示意圖。

圖4-1 半拉格朗日對流法

  半拉格朗日對流法的第一步就是要找出 x ⃗ P \vec x_P x

P​，為此我們根據 x ⃗ G \vec x_G x

G​做反向的追蹤。粒子位置對時間的導數就是速度場：

(4.6) d x ⃗ d t = u ⃗ ( x ⃗ ) \frac{d\vec x}{dt}=\vec u(\vec x) \tag {4.6} dtdx

​=u

(x

)(4.6)

  經過一個時間步長 Δ t \Delta t Δt之後，粒子由 x ⃗ P \vec x_P x

P​移動到 x ⃗ G \vec x_G x

G​。為了得到 x ⃗ P \vec x_P x

P​，最簡單的方法就是采用前向歐拉法進行倒推：

(4.7) x ⃗ P = x ⃗ G − Δ t u ⃗ ( x ⃗ G ) \vec x_P=\vec x_G-\Delta t\vec u(\vec x_G) \tag {4.7} x

P​=x

G​−Δtu

(x

G​)(4.7)

  然而前向歐拉法隻有一階的精度，若在不改變 Δ t \Delta t Δt的情況下提高精度，我們可以采用高階的龍格庫塔法（Runge-Kutta method）。采用二階的龍格庫塔法如下所示：

(4.7) x ⃗ m i d = x ⃗ G − 1 2 Δ t u ⃗ ( x ⃗ G ) , x ⃗ P = x ⃗ G − Δ t u ⃗ ( x ⃗ m i d ) . \vec x_{mid}=\vec x_G-\frac12\Delta t\vec u(\vec x_G),\\ \vec x_P=\vec x_G-\Delta t\vec u(\vec x_{mid}). \tag {4.7} x

mid​=x

G​−21​Δtu

(x

G​),x

P​=x

G​−Δtu

(x

mid​).(4.7)

  倒推得到 Δ t \Delta t Δt之前的網格位置 x ⃗ P \vec x_P x

P​一般不會恰好在網格頂點上，為此我們需要做些插值。三維模拟通常采用三線性插值，而二維的則采用雙線性插值。

(4.8) q G n + 1 = i n t e r p o l a t e ( q n , x ⃗ P ) q^{n+1}_G=interpolate(q_n,\vec x_P) \tag {4.8} qGn+1​=interpolate(qn​,x

P​)(4.8)

2、邊界情況

  若我們倒推得到的 x ⃗ P \vec x_P x

P​仍然在流體的内部，那麼做插值是完全沒問題的。但若 x ⃗ P \vec x_P x

P​在流體的邊界之外呢？這種情況的出現的原因通常有兩個：一個是 x ⃗ P \vec x_P x

P​确确實實在流體的外部且即将流入流體内部，另一個是由前向歐拉法或龍格庫塔法的數值計算方法帶來的誤差導緻。

  在一種情況下，我們應該知道當流體流入時其攜帶的實體量，此時我們将這個外部流入的實體量作為傳回值即可。例如，第 n n n個時間步時的外部流體以速度 U ⃗ \vec U U

和溫度 T T T在第 n + 1 n+1 n+1個時間步時注入流體内部 x ⃗ G \vec x_G x

G​的位置，那麼 T ⃗ G n + 1 \vec T^{n+1}_G T

Gn+1​的值就為 T T T。

  在第二種由誤差導緻的情況下，一個适當的政策就是根據邊界上的最近點外推出所求得實體量。在模拟某些流體時，外推變得很簡單。例如，在模拟煙霧時我們簡單地假設煙霧流體外部即空氣的速度風場為某個常數 U ⃗ \vec U U

（可能為 0 0 0），這樣邊界上的速度場都取 U ⃗ \vec U U

。但還有一些必須根據流體内部的已知量外推出未知量，這時情況就變得比較複雜了。具體如何外推将在後面介紹，目前我們隻需要知道大概的步驟：首先尋找邊界上的最近點，然後在最近點的領域内插值擷取相應的實體量場。

3、時間步長大小

  對任何一種數值計算方法的主要的考慮點就是它是否穩定。幸運的是，半拉格朗日對流法已經被證明是一種無條件穩定的算法：無論 Δ t \Delta t Δt取多大，它永遠不會出現數值爆炸的現象。因為每一個新值 q q q的确定，都是通過對舊值得插值，無論是線性插值、雙線性插值還是三線性插值， q q q的大小都是處于插值點之間，不會得到比原來插值點更大或者更小的值，因而 q q q是有上下界的。這使得我們可以盡情地根據所需的模拟品質和模拟效率去調整時間步長。

  但是在實踐中，時間步長的大小也不能選得太過極端，否則會産生一些奇觀的現象。Foster和Fekiw提出了一個對 Δ t \Delta t Δt的限制：流體粒子在 Δ t \Delta t Δt内的倒推軌迹最多經過某個常數個網格單元為宜，例如5個：

(4.9) Δ t ≤ 5 Δ x u m a x \Delta t \leq \frac{5\Delta x}{u_{max}} \tag {4.9} Δt≤umax​5Δx​(4.9)

  公式 ( 4.9 ) (4.9) (4.9)中， u m a x u_{max} umax​是速度場的最大值，我們可以簡單地取 存儲在網格中的最大速度值。一個更魯棒的方法考慮了體積力（如重力、浮力等）對最大速度的影響：

(4.10) u m a x = m a x ( ∣ u n ∣ ) + Δ t ∣ g ∣ u_{max}=max(|u^n|)+\Delta t|g| \tag {4.10} umax​=max(∣un∣)+Δt∣g∣(4.10)

  将不等式 ( 4.9 ) (4.9) (4.9)的最大值帶入公式 ( 4.10 ) (4.10) (4.10)，我們有：

(4.11) u m a x = m a x ( ∣ u n ∣ ) + 5 Δ x u m a x ∣ g ∣ u_{max}=max(|u^n|)+\frac{5\Delta x}{u_{max}}|g| \tag {4.11} umax​=max(∣un∣)+umax​5Δx​∣g∣(4.11)

  取一個簡單的速度上界（簡化了公式 ( 4.11 ) (4.11) (4.11)）， u m a x u_{max} umax​：

(4.12) u m a x = m a x ( ∣ u n ∣ ) + 5 Δ x g u_{max}=max(|u^n|)+\sqrt{5\Delta xg} \tag {4.12} umax​=max(∣un∣)+5Δxg

​(4.12)

  這樣確定了 u m a x u_{max} umax​始終為正，且避免公式 ( 4.9 ) (4.9) (4.9)的除 0 0 0錯誤。

  關于時間步長的讨論離不開 C F L CFL CFL（以Courant、Friedrichs、Lewy三人的名字命名）條件。 C F L CFL CFL條件是一個簡單而直覺的判斷計算是否收斂的必要條件。它的直覺實體解釋就是時間推進求解的速度必須大于實體擾動傳播的速度，隻有這樣才能将實體上所有的擾動俘獲到。滿足 C F L CFL CFL條件意味着當 Δ x \Delta x Δx和 Δ t \Delta t Δt趨于取極限 0 0 0時，數值計算所求的解就會收斂到原微分方程的解。

  對于半拉格朗日對流法，其滿足 C F L CFL CFL條件當且僅當在極限情況下，追蹤得到的粒子軌迹足夠逼近真實的軌迹。足夠逼近的意思是經過正确的網格插值能夠得到正确的依賴域（即差分格式的依賴域包含了原微分方程的依賴域），追蹤的軌迹就會收斂到正确真實的軌迹。

  因而，對于采用标準的顯式有限差分法的對流方程求解，為了保證收斂，我們要求 q n + 1 q^{n+1} qn+1的新值是由以目前網格點為中心、以 C Δ x C\Delta x CΔx（ C C C是一個小的整數常量）為半徑的鄰域範圍内插值得到：

(4.13) Δ t ≤ C Δ x ∣ u ⃗ ∣ \Delta t \leq C\frac{\Delta x}{|\vec u|} \tag {4.13} Δt≤C∣u

∣Δx​(4.13)

  公式 ( 4.13 ) (4.13) (4.13)中的 C C C被稱為 C F L CFL CFL數，因而不等式 ( 4.9 ) (4.9) (4.9)可以看成是公式 ( 4.13 ) (4.13) (4.13)取 C F L CFL CFL數為 5 5 5得到。

4、數值耗散

  對流算法在對流擷取新的實體量場 q i n + 1 q^{n+1}_i qin+1​時會進行一些插值操作，插值不可避免地會平滑實體量場，這帶來了一些數值耗散。一次兩次的數值耗散不會由太大的影響，但是在流體模拟中我們會在每個時間步都進行對流運算，反反複複的平滑操作将數值耗散不斷擴大，損失大量的流體細節。

  以一維的對流項計算為例，流體速度為常量 u > 0 u>0 u>0：

(4.14) ∂ q ∂ t + u ∂ q ∂ x = 0 \frac{\partial q}{\partial t}+u\frac{\partial q}{\partial x}=0 \tag {4.14} ∂t∂q​+u∂x∂q​=0(4.14)

  假設 Δ t &lt; Δ x u \Delta t &lt; \frac{\Delta x}{u} Δt<uΔx​，即單個時間步長内粒子追蹤軌迹長度小于單個網格單元的大小。我們的目标點是 x i x_i xi​，則倒推得到的粒子位置就落在了 [ x i − 1 , x i ] [x_{i-1},x_i] [xi−1​,xi​]上的 x i − Δ t u x_i-\Delta tu xi​−Δtu，然後進行線性插值得到 q i n + 1 q^{n+1}_i qin+1​：

(4.15) q n + 1 = Δ t u Δ x q i − 1 n + ( 1 − Δ t u Δ x ) q i n q^{n+1}=\frac{\Delta tu}{\Delta x}q^n_{i-1}+(1-\frac{\Delta tu}{\Delta x})q^n_i \tag {4.15} qn+1=ΔxΔtu​qi−1n​+(1−ΔxΔtu​)qin​(4.15)

  将公式 ( 4.15 ) (4.15) (4.15)整理一下，有：

(4.16) q i n + 1 = q i n − Δ t u q i n − q i − 1 n Δ x q^{n+1}_i=q^n_i-\Delta tu\frac{q^n_i-q^n_{i-1}}{\Delta x} \tag {4.16} qin+1​=qin​−ΔtuΔxqin​−qi−1n​​(4.16)

  公式 ( 4.16 ) (4.16) (4.16)實際上正好就是采用時間上的前向歐拉差分法和空間上的單向有限差分法的歐拉方案，把 q i n q^n_i qin​看成是 q n q^n qn關于 x i x_i xi​的函數，對 q i − 1 n q^n_{i-1} qi−1n​進行泰勒級數展開：

(4.17) q i − 1 n = q i n − ( ∂ q ∂ x ) i n Δ x + ( ∂ 2 q ∂ x 2 ) i n Δ x 2 2 + O ( Δ x 3 ) q^n_{i-1}=q^n_i-(\frac{\partial q}{\partial x})^n_i\Delta x+(\frac{\partial^2q}{\partial x^2})^n_i\frac{\Delta x^2}{2}+O(\Delta x^3) \tag {4.17} qi−1n​=qin​−(∂x∂q​)in​Δx+(∂x2∂2q​)in​2Δx2​+O(Δx3)(4.17)

  将公式 ( 4.17 ) (4.17) (4.17)代入公式 ( 4.16 ) (4.16) (4.16)，并做一些變量消去，可得：

(4.18) q i n + 1 = q i n − Δ t u ( ∂ q ∂ x ) i n + Δ t u Δ x ( ∂ 2 q ∂ x 2 ) i n + O ( Δ x 2 ) q^{n+1}_i=q^n_i-\Delta tu(\frac{\partial q}{\partial x})^n_i+\Delta tu\Delta x(\frac{\partial^2q}{\partial x^2})^n_i+O(\Delta x^2) \tag {4.18} qin+1​=qin​−Δtu(∂x∂q​)in​+ΔtuΔx(∂x2∂2q​)in​+O(Δx2)(4.18)

  在二階截斷誤差的情況下，結合公式 ( 4.18 ) (4.18) (4.18)和公式 ( 4.14 ) (4.14) (4.14)，有：

(4.19) ∂ q ∂ t + u ∂ q ∂ x = u Δ x ( ∂ 2 q ∂ x 2 ) \frac{\partial q}{\partial t}+u\frac{\partial q}{\partial x}=u\Delta x(\frac{\partial^2q}{\partial x^2}) \tag {4.19} ∂t∂q​+u∂x∂q​=uΔx(∂x2∂2q​)(4.19)

  右邊就是對流方程計算時引入的額外類似粘度乘上系數 u Δ x u\Delta x uΔx的項。這也就是說，當我們采用簡單的半拉格朗日法去求解無粘度的對流方程時，模拟的結果卻看起來我們像時在模拟有粘度的流體。這就是數值耗散! 當然，當 Δ x → 0 \Delta x\to 0 Δx→0時，這個數值耗散系數也會趨于 0 0 0，是以取時間步無窮小時能夠得到正确的模拟結果，但這需要耗費巨額的計算資源開銷。我們通常模拟的流體大多數都是無粘度的，是以如何減少這個數值耗散是個至關重要的難題。

  一個簡單有效的修複數值耗散的方法就是采用更加銳利的插值方法，進而盡可能地減少由插值帶來的數值耗散。在一維的情況時，我們采用三次插值（cubic interpolant）如下公式 ( 4.21 ) (4.21) (4.21)，而不是簡單的一次線性插值 ( 4.20 ) (4.20) (4.20)：

(4.20) q ≈ ( 1 − s ) x i + s x i + 1 q\approx(1-s)x_i+sx_{i+1} \tag {4.20} q≈(1−s)xi​+sxi+1​(4.20)

(4.21) q ≈ [ − 1 3 s + 1 2 s 2 − 1 6 s 3 ] q i − 1 + [ 1 − s 2 + 1 2 ( s 3 − s ) ] q i + [ s + 1 2 ( s 2 − s 3 ) ] q i + 1 + [ 1 6 ( s 3 − s ) ] q i + 2 q\approx[-\frac13s+\frac12s^2-\frac16s^3]q_{i-1}+[1-s^2+\frac12(s^3-s)]q_i\\ +[s+\frac12(s^2-s^3)]q_{i+1}+[\frac16(s^3-s)]q_{i+2} \tag {4.21} q≈[−31​s+21​s2−61​s3]qi−1​+[1−s2+21​(s3−s)]qi​+[s+21​(s2−s3)]qi+1​+[61​(s3−s)]qi+2​(4.21)

  擴充到二維或者三維就是雙三次插值（bicubic interpolation）或三三次插值（tricubic interpolation）。以二維情況為例，我們可以先沿着 x x x軸做第一遍的三次插值如公式 ( 4.22 ) (4.22) (4.22)，然後再沿着 y y y軸做第二遍插值如公式 ( 4.23 ) (4.23) (4.23)：

(4.22) q j − 1 = w − 1 ( s ) q i − 1 , j − 1 + w 0 ( s ) + q i , j − 1 + w 1 ( s ) q i + 1 , j − 1 + w 2 ( s ) q i + 2 , j − 1 , q j = w − 1 ( s ) q i − 1 , j + w 0 ( s ) + q i , j + w 1 ( s ) q i + 1 , j + w 2 ( s ) q i + 2 , j , q j + 1 = w − 1 ( s ) q i − 1 , j + 1 + w 0 ( s ) + q i , j + 1 + w 1 ( s ) q i + 1 , j + 1 + w 2 ( s ) q i + 2 , j + 1 , q j + 2 = w − 1 ( s ) q i − 1 , j + 2 + w 0 ( s ) + q i , j + 2 + w 1 ( s ) q i + 1 , j + 2 + w 2 ( s ) q i + 2 , j + 2 . q_{j-1}=w_{-1}(s)q_{i-1,j-1}+w_0(s)+q_{i,j-1}+w_1(s)q_{i+1,j-1}+w_2(s)q_{i+2,j-1},\\ q_{j}=w_{-1}(s)q_{i-1,j}+w_0(s)+q_{i,j}+w_1(s)q_{i+1,j}+w_2(s)q_{i+2,j},\\ q_{j+1}=w_{-1}(s)q_{i-1,j+1}+w_0(s)+q_{i,j+1}+w_1(s)q_{i+1,j+1}+w_2(s)q_{i+2,j+1},\\ q_{j+2}=w_{-1}(s)q_{i-1,j+2}+w_0(s)+q_{i,j+2}+w_1(s)q_{i+1,j+2}+w_2(s)q_{i+2,j+2}. \tag {4.22} qj−1​=w−1​(s)qi−1,j−1​+w0​(s)+qi,j−1​+w1​(s)qi+1,j−1​+w2​(s)qi+2,j−1​,qj​=w−1​(s)qi−1,j​+w0​(s)+qi,j​+w1​(s)qi+1,j​+w2​(s)qi+2,j​,qj+1​=w−1​(s)qi−1,j+1​+w0​(s)+qi,j+1​+w1​(s)qi+1,j+1​+w2​(s)qi+2,j+1​,qj+2​=w−1​(s)qi−1,j+2​+w0​(s)+qi,j+2​+w1​(s)qi+1,j+2​+w2​(s)qi+2,j+2​.(4.22)

(4.23) q = w − 1 ( t ) q j − 1 + w 0 ( t ) q j + w 1 ( t ) q j + 1 + w 2 ( t ) q j + 2 q=w_{-1}(t)q_{j-1}+w_0(t)q_j+w_1(t)q_{j+1}+w_2(t)q_{j+2} \tag {4.23} q=w−1​(t)qj−1​+w0​(t)qj​+w1​(t)qj+1​+w2​(t)qj+2​(4.23)

  當然也可以先沿着 y y y軸，然後再沿着 x x x軸做插值操作。