流體模擬：基礎

目錄

• 流體模擬
  • 1.矢量微積分
    • 1.1梯度（Grant）
    • 1.2散度（Divergence）
    • 1.3旋度(Curl)
    • 1.4 拉普拉斯算子（Laplacian）
  • 2.Navier-Stokes偏微分方程
    • 2.1符號標記
    • 2.2動量方程（The Momentum Equation）
    • 2.3 拉格朗日描述與歐拉描述
    • 2.4不可壓縮性
    • 2.5 丟棄粘度項
    • 2.6 邊界條件
  • 3.N-S方程的分步求解
    • 3.1 時間步長
    • 3.2 網格結構
  • 4.對流算法
    • 4.1半拉格朗日對流算法（Semi-Lagrangian Advection）
    • 4.2 邊界情況
    • 4.3時間步長大小
    • 4.4數值耗散
• 參考資料

流體模擬

1.矢量微積分

1.1梯度（Grant）

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

\[\nabla f(x, y)=\left(\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y}\right) \tag{1.1}\label{1.1} \]

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

\[\nabla f=\frac{\partial f}{\partial \vec{x}} \tag{1.2}\label{1.2} \]

有限微分形式(中心差分法)如下

\[\frac{f_{i+1, j}-f_{i-1, j}}{2 \delta x}, \frac{f_{i, j+1}-f_{i, j-1}}{2 \delta y} \tag{1.3}\label{1.3} \]

1.2散度（Divergence）

散度算子僅僅應用於向量場，它衡量在某一點處相應的矢量聚集或者發散程度，測量方向為徑向，結果為標量。同時散度有重要的物理意義，它是“密度”離開給定空間區域的速率。

2維、3維形式的散度算子如下所示

\[\begin{gathered} \nabla \cdot \vec{u}=\nabla \cdot(u, v)=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y} \\ \nabla \cdot \vec{u}=\nabla \cdot(u, v, w)=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z} \end{gathered} \tag{1.4}\label{1.4} \]

散度符號\(\nabla \cdot \vec{u}\)可以理解為梯度$\nabla \(與矢量\)\vec{u}$的點乘

\[\nabla \cdot \vec{u}=\left(\frac{\partial}{\partial x}, \frac{\partial}{\partial y}, \frac{\partial}{\partial z}\right) \cdot(u, v, w)=\frac{\partial}{\partial x} u+\frac{\partial}{\partial y} v+\frac{\partial}{\partial z} w \tag{1.5}\label{1.5} \]

有限微分形式(中心差分法)：

\[\frac{u_{i+1, j}-u_{j-1, j}}{2 \delta x}+\frac{v_{i, j+1}-v_{i, j-1}}{2 \delta y} \tag{1.6}\label{1.6} \]

若矢量場散度為0，則稱該矢量場無散度。在后面的Navier-Stokes方程中，他被應用於流速，他測量的是經過一小塊流體表面的速率變化，讓流體的的散度永遠是零（即不可壓縮狀態）。

1.3旋度(Curl)

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

\[\nabla \times \vec{u}=\nabla \times(u, v, w)=\left(\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}\right) \tag{1.7}\label{1.7} \]

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

\[\nabla \times \vec{u}=\nabla \times(u, v)=\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y} \tag{1.8}\label{1.8} \]

同樣地,旋度符號\(\nabla \times \vec{u}\)我們可以理解為梯度\(\nabla\)與矢量場\(\vec{u}\)的叉乘:

\[\nabla \times \vec{u}=\left(\frac{\partial}{\partial x}, \frac{\partial}{\partial y}, \frac{\partial}{\partial z}\right) \times(u, v, w) \tag{1.9}\label{1.9} \]

若矢量場旋度為0，則稱該矢量場無旋度。若流體的旋度為0，即流體是不可旋轉的狀態。

1.4 拉普拉斯算子（Laplacian）

拉普拉斯算子定義為梯度的散度，符號表示為\(\nabla \cdot \nabla\),顯然\(\nabla \cdot\)是散度,而后面的\(\nabla\)則為梯度，故拉普拉斯算子即梯度的散度，是一個二階微分算子。2維、3維形式分別如下：

\[\begin{gathered} \nabla \cdot \nabla f=\frac{\partial^{2} f}{\partial x^{2}}+\frac{\partial^{2} f}{\partial y^{2}} \\ \nabla \cdot \nabla f=\frac{\partial^{2} f}{\partial x^{2}}+\frac{\partial^{2} f}{\partial y^{2}}+\frac{\partial^{2} f}{\partial z^{2}} \end{gathered} \tag{1.10}\label{1.10} \]

偏微分方程\(\nabla \cdot \nabla f=0\)被稱為拉普拉斯方程;若右邊為某個非0常數，即\(\nabla \cdot \nabla f=q\),如熱傳導方程，我們稱之為泊松方程。更一般地，如果梯度再乘上一個標量a（如\(\frac {1}{\rho}\))，即\(\nabla \cdot (a\nabla f)=q\)，我們依舊稱之為泊松問題。

有限微分形式(中心差分法)（2維形式）：

\[\frac{f_{i+1, i}-2 f_{i, j}+f_{i-1, j}}{(\delta x)^{2}}+\frac{f_{i, j+1}-2 f_{i, i}+f_{i, j-1}}{(\delta y)^{2}} \tag{1.11}\label{1.11} \]

如果網格單元是正方形（即\(\delta x=\delta y\)）,拉普拉斯算子可簡化為：

\[\nabla^{2} p=\frac{f_{i+1, j}+f_{i-1, j}+f_{i, j+1}+f_{i, j-1}-4 f_{i, j}}{(\delta x)^{2}} \tag{1.12}\label{1.12} \]

2.Navier-Stokes偏微分方程

流體模擬器的構建主要圍繞著名的不可壓縮Navier−StokesNavier−Stokes方程展開，它是一個流體力學領域的偏微分方程，方程形式如下：

\[\begin{eqnarray*} \frac{\partial \vec{u}}{\partial t}+\vec{u} \cdot \nabla \vec{u}+\frac{1}{\rho} \nabla p=\vec{g}+\nu \nabla \cdot \nabla \vec{u} \tag{2.1} \label{2.1} \\ \nabla \cdot \vec{u}=0 \tag{2.2} \label{2.2} \end{eqnarray*} \]

現在我們將其剖析成一個個比較簡單的部分來理解。

2.1符號標記

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

希臘字符\(\rho\)是流體的密度，對於水，該值大約為\(1000kg=m/s^3\)，而空氣則大約為\(1.3kg=m/s^3\)

字符\(p\)代表壓力，流體對任何物體施加的單位面積力

字符\(\vec{g}\)則是我們熟悉的重力加速度，通常取\((0,-9.81,0)m/s^3\)。我們約定y軸向上，而x軸和z軸在水平面上。另外補充一點，我們把其他的一些類似的力都累加\(\vec{g}\)上，也就是我們統一\(\vec{g}\)表示所有類似力之和，這類力我們稱之為體積力（稱之為體積力是因為它們的力是作用到整個流體而不只是流體的表面）；

希臘字符\(\vec{v}\)流體的運動粘度，它衡量流體的黏滯性。糖蜜之類的流體有非常高的粘度，而像酒精之類的流體有很低的粘度;

2.2動量方程（The Momentum Equation）

偏微分方程\(\eqref{2.1}\)我們稱之為動量方程，它本質上就是我們熟悉的牛頓定律\(\vec{F}=m\vec{a}\)的形式，描述了施加在流體上的力是如何影響流體的運動。

假設我們用粒子系統來模擬流體，每個粒子代表流體的一小滴，每個粒子有各自的質量m、體積V和速度\(\vec{u}\)。為了讓整個粒子系統運作起來，我們必須弄清楚每個粒子所承受的力的作用。牛頓第二定律告訴我們：\(\vec{F}=m\vec{a}\)，而根據加速度定義，我們有

\[\vec{a}=\frac{D\vec{u}}{Dt} \tag{2.3}\label{2.3} \]

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

\[m\frac{D\vec{u}}{Dt}=\vec{F} \tag{2.4}\label{2.4} \]

流體粒子承受的力有：

• 一個最簡單的力就是重力：\(m\vec{g}\)

• 其他的流體質點（或其他流體粒子）也會對當前的流體粒子產生作用力:

  流體內部的相互作用力之一便是壓力，較大壓力的區域會向較低壓力區域產生作用力。值得注意的是，我們只關注施加在粒子上的壓力的凈合力。例如,若施加在粒子上壓力在每個方向上都相等，那么它的壓力的合力便為0。

  我們用壓力的負梯度（取負是因為方向是由壓力大的區域指向壓力小的區域）來衡量在當前流體粒子處壓力的不平衡性，即取\(-\nabla p\)。那么流體粒子所承受的壓力就是對\(-\nabla p\)在整個流體粒子的體積上進行積分,為了簡化，我們簡單地將\(V\)與\(-\nabla p\)相乘，故粒子壓力部分為\(-\nabla p\)。

• 其他的流體相互作用力則是由流體的黏性產生的：

  我們直觀地把這種力理解為盡可能使得粒子以周圍區域的平均速度移動的力，也就是使得粒子的速度與周圍區域粒子速度的差距最小化

  拉普拉斯算子是衡量一個量與之周圍區域該量平均值之差的算符，因而\(\nabla \cdot \nabla \vec{u}\)是當前粒子速度矢量與周圍區域平均速度矢量之差。

  為了計算粘滯力，我們同樣對\(\nabla \cdot \nabla \vec{u}\)整個粒子體積\(V\)上進行積分，與前面類似，我們簡單取\(V\nabla \cdot \nabla \vec{u}\)。除此之外，我們還引進一個稱為動力粘度系數的物理量，符號為\(\mu\)。因而粘滯力為\(V\mu\nabla \cdot \nabla \vec{u}\)

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

\[m \frac{D \vec{u}}{D t}=\vec{F}=m \vec{g}-V \nabla p+V \mu \nabla \cdot \nabla \vec{u} \tag{2.5}\label{2.5} \]

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

\[\rho \frac{D \vec{u}}{D t}=\rho \vec{g}- \nabla p+ \mu \nabla \cdot \nabla \vec{u} \tag{2.5}\label{2.6} \]

兩邊同時除以\(\rho\).同時定義運動粘度為\(v=\mu/\rho\),變成了

\[\frac{D \vec{u}}{D t}+ \frac{1}{\rho} \nabla p= \vec{g}+ v \nabla \cdot \nabla \vec{u} \tag{2.7}\label{2.7} \]

2.3 拉格朗日描述與歐拉描述

現在我們要把物質導數（Material derivative）\(\frac{D \vec{u}}{D t}\)弄清楚，我們需要了解兩種描述方法：拉格朗日描述和歐拉描述。

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

• 拉格朗日描述（ Lagrangian viewpoint）

這種描述方法把物體看成是由類似於粒子系統的形式組成，固體或流體的每個點看作一個獨立的粒子，粒子有各自相應的位置\(\vec{x}\)和速度\(\vec{u}\)。我們可以把粒子理解為組成物體的分子。對於我們通常采用拉格朗日描述法進行建模模擬，即用一系列離散的粒子集來構建，粒子之間通過網格相聯系。

• 歐拉描述（Eulerian viewpoint）

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

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

把拉格朗日描述法和歐拉描述法聯系起來的關鍵點就是物質導數。

從拉格朗日描述法出發：假設有一群粒子，每個粒子都有各自的位置\(\vec{x}\)和速度\(\vec{v}\),記\(q\)為通用的物理量（如密度、速度和溫度等），每個粒子有其對應的\(q\)值，方程\(q(t,\vec{x})\)描述在時間點\(t\)而位置\(\vec{x}\)的粒子對應的物理量值\(q\)。我們取對時間\(t\)的導數（注意用到了求導鏈式法則，以及\(\frac{\partial q}{\partial \vec{x}}=\nabla q\)和\(\vec{u}=\frac{d\vec{x}}{dt}\)):

\[\frac{d}{d t} q(t, \vec{x})=\frac{\partial q}{\partial t}+\nabla q \cdot \frac{d \vec{x}}{d t}=\frac{\partial q}{\partial t}+\nabla q \cdot \vec{u} \equiv \frac{D q}{D t} \tag{2.8}\label{2.8} \]

這就是物質導數。將其代入式子我們就得到了流體動量方程。物質導數針對的是流體質點在這里就是流體粒子）而不是空間的固定點，上式子寫得完整些就是：

\[\frac{D q}{D t}=\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.9}\label{2.9} \]

對於給定的速度場\(\vec{u}\),流體的物理性質如何在這個速度場\(\vec{u}\)下的變化計算我們稱之為對流（advection）。一個最簡單的對流方程，就是其物理量的物質導數為0，如下所示：

\[\frac{D q}{D t}=0 \Longrightarrow \frac{\partial q}{\partial t}+\vec{u} \cdot \nabla q=0 \tag{2.10}\label{2.10} \]

上式的意義即在拉格朗日視角觀察下，每個流體粒子的物理量保持不變。

當物質導數作用於向量時，如：RGB,速度場時，一個簡單方法就是分別對待每個分量。公式如下：

\[\frac{D \vec{u}}{D t}=\left[\begin{array}{c} D u / D t \\ D v / D t \\ D w / D t \end{array}\right]=\left[\begin{array}{c} \partial u / \partial t+\vec{u} \cdot \nabla u \\ \partial v / \partial t+\vec{u} \cdot \nabla v \\ \partial w / \partial t+\vec{u} \cdot \nabla w \end{array}\right]=\frac{\partial \vec{u}}{\partial t}+\vec{u} \cdot \nabla \vec{u} \tag{2.11}\label{2.11} \]

如果你想知道偏導數的具體細節：

\[\frac{D \vec{u}}{D t}=\left[\begin{array}{c} \frac{\partial u}{\partial t}+u \frac{\partial u}{\partial x}+v \frac{\partial u}{\partial y}+w \frac{\partial u}{\partial z} \\ \frac{\partial v}{\partial t}+u \frac{\partial v}{\partial x}+v \frac{\partial v}{\partial y}+w \frac{\partial v}{\partial z} \\ \frac{\partial w}{\partial t}+u \frac{\partial w}{\partial x}+v \frac{\partial w}{\partial y}+w \frac{\partial w}{\partial z} \end{array}\right] \tag{2.12}\label{2.12} \]

2.4不可壓縮性

關於流體的壓縮性的物理細節描述，可以看下3blue1brown的視頻散度與旋度：麥克斯韋方程組、流體等所用到的語言

我們需要知道通常情況下流體的體積變化非常小（除開一些極端的情況，而且這些極端情況我們日常生活中較少出現）。可壓縮流體的模擬涉及到非常復雜的情況，往往需要昂貴的計算資源開銷，為此在計算機流體模擬中我們通常把所有的流體當作是不可壓縮的，即它們的體積不會發生變化。

任取流體的一部分，設其體積為\(\Omega\)而其邊界閉合曲面為\(\partial\Omega\),我們可以通過圍繞邊界曲面\(\partial\Omega\)對流體速度\(\vec{u}\)在曲面法線方向上的分量進行積分來衡量這塊部分流體的體積變化速率：

\[\frac{d}{d t} \operatorname{Volume}(\Omega)=\iint_{\partial \Omega} \vec{u} \cdot n \tag{2.13}\label{2.13} \]

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

\[\iint_{\partial \Omega} \vec{u} \cdot n=0 \tag{2.14}\label{2.14} \]

由高斯散度定理（高斯公式），我們可以把上式轉化為體積分：

\[\iint_{\partial \Omega} \vec{u} \cdot n=\iiint_{\Omega} \nabla \cdot \vec{u}=0 \tag{2.15}\label{2.15} \]

上式對任意的\(\Omega\)成立，即無論\(\Omega\)取何值，積分值均為0。這種情況下只有令積分函數取0方可成立，即對0積分無論\(\Omega\)取和值結果均為0。所以有：

\[\nabla \cdot \vec{u}=0 \]

這就是Navier−StokesNavier−Stokes方程中的不可壓縮條件\(\eqref{2.2}\)。滿足不可壓縮條件的速度場被稱為是無散度的，即在該速度場下流體體積既不膨脹也不坍縮，而是保持在一個常量。模擬不可壓縮流體的關鍵部分就是使得流體的速度場保持無散度的狀態，這也是流體內部壓力的來源。

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

\[\nabla \cdot \frac{\partial \vec{u}}{\partial t}+\nabla \cdot(\vec{u} \cdot \nabla \vec{u})+\nabla \cdot \frac{1}{\rho} \nabla p=\nabla \cdot(\vec{g}+\nu \nabla \cdot \nabla \vec{u}) \tag{2.16}\label{2.16} \]

對於上式第一項，我們轉變一下求導次序：

\[\frac{\partial}{\partial t}\nabla \cdot\vec{u} \tag{2.17}\label{2.17} \]

如果滿足流體不可壓縮條件，那么上式取值0（無散度），我們調整下上式可得關於壓力的方程：

\[\nabla \cdot \frac{1}{\rho} \nabla p=\nabla \cdot(- \vec{u} \cdot \nabla \vec{u}+\vec{g}+\nu \nabla \cdot \nabla \vec{u}) \label{2.18} \tag{2.18} \]

2.5 丟棄粘度項

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

丟棄了粘度項的Navier−StokesNavier−Stokes方程被稱為歐拉方程，而這種理想的流體則是無粘度的。丟棄了粘度項的歐拉方程如下：

\[\begin{eqnarray*} \frac{D \vec{u}}{D t}+\frac{1}{\rho} \nabla p=\vec{g} \tag{2.19}\label{2.19} \\ \nabla \cdot \vec{u}=0 \tag{2.20}\label{2.20} \end{eqnarray*} \]

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

2.6 邊界條件

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

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

\[\vec{u} \cdot n=0 \tag{2.21}\label{2.21} \]

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

\[\vec{v} \cdot n=\vec{u_{solid}} \cdot n \tag{2.22}\label{2.22} \]

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

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

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

\[[p]=\lambda k \tag{2.23}\label{2.23} \]

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

\[p=\lambda k\tag{2.24}\label{2.24} \]

3.N-S方程的分步求解

了對以上對Navier−StokesNavier−Stokes方程的理論支撐，接下來我們就要如何用計算機來對該組偏微分方程進行離散化求解。為了程序的松耦合性以及使計算盡可能地高效、簡單，在流體模型領域，我們將流體方程分成幾個獨立的步驟，然后按順序先后推進。對於不可壓縮的無粘度流體方程（即前面的歐拉方程），我們將其離散化成對流項（advection）\(\eqref{3.1}\)，體積力項（body force）\(\eqref{3.2}\)，壓力/不可壓縮項 \(\eqref{3.3}\)：

\[\begin{eqnarray*} \frac{D q}{D t}=0 \tag{3.1} \label{3.1} \\ \frac{\partial \vec{u}}{\partial t}=\vec{g} \tag{3.2} \label{3.2}\\ \left\{\begin{array}{l} \frac{\partial \vec{u}}{\partial t}+\frac{1}{\rho} \nabla p=0 \\ \nabla \cdot \vec{u}=0 \end{array}\right. \tag{3.3} \label{3.3} \end{eqnarray*} \]

• 對於對流項：我們用了一個通用量的符號\(q\)是因為我們不僅僅要對流體的速度進行對流，還需要對其他物理量進行對流。我們記對流項公式\(\eqref{3.1}\)的對流計算算法為\(advect(\vec{u},\triangle t,q)\)，即對於給定的時間步長\(\triangle t\)和速度\(\vec{u}\)，對物理量\(q\)進行對流。
• 對於體積力項：我們采用簡單的前向歐拉法即可：\(\vec{u} \leftarrow \vec{u}+g \triangle t\)
• 對於壓力/不可壓縮項:我們用一個稱為\(project(\triangle t,\vec{u})\)的算法，通過\(project(\triangle t,\vec{u})\)計算出正確的壓力以確保速度場\(\vec{u}\)的無散度性質。歐拉方案不會着重研究具體粒子間的作用力，因而不會正向去求解\(\frac{1}{\rho} \nabla p\),它是利用流體不可壓縮的特性，將速度場\(\vec{u}\)投影到散度為0的空間上，間接地解算了壓力項。這種思想相當於，已知一個中間量\(\vec{u_{temp}}\),對這個中間量的唯一一個操作（如正向求解壓力\(\frac{1}{\rho} \nabla p\)）不可行，但是知道最終量\(\vec{u_{final}}\)符號的一個性質（散度為0），於是只要將\(\vec{u_{temp}}\)投影到符合散度為0的特性平面上，即可間接地還原正向求解壓力的操作，得到最終的速度場\(\vec{u_{temp}}\)。

對流項\(advect(\vec{u},\triangle t,q)\)的輸入速度場\(\vec{u}\)要確保為無散度的狀態，投影項\(project(\triangle t,\vec{u})\)確保流體體積保持不變，因而投影項輸出的速度場必然是無散度。所以我們只要確保投影項\(project(\triangle t,\vec{u})\)輸出的速度\(\vec{u}\)作為對流項\(advect(\vec{u},\triangle t,q)\)的輸入即可，這時我們的分步求解流體方程的優勢就體現出來了，其偽代碼如下所示:

算法1：Fluid Simulation(\(\vec{u_n},\triangle t\))

1. 初始化速度場\(\vec{u_n}\),使得\(\vec{u_n}\)無散度
2. 對於每個時間步\(n=0,1,2...\)
3. 決定一個合理的時間步長\(\triangle t=t_{n+1}-t_n\)
4. 對流項計算$ \vec {u_A}= advect(\vec{u_n},\triangle t,q)$
5. 體積力項計算\(\vec {u_B}=project(\triangle t,\vec{u})\)
6. 無散度計算\(\vec {u_{n+1}}=project(\triangle t,\vec{u_b})\)

3.1 時間步長

在流體模擬算法中，確定適當的時間步長是算法的第一步。因為計算流體模擬的最后結果是呈現在屏幕上的，所以\(\triangle t\)的選取與屏幕的刷新率有重要的關系。若選取的\(\triangle t\)有\(t_n+\triangle t>t_{frame}\),那么必須做一個截斷使\(\triangle t=t_{frame}-t_n\)。此外，流體模擬的三個步驟即對流項、體積力項、無散度投影項對時間步長\(\triangle t\)的要求不盡相同，要選擇一個滿足所有要求的最小時間步長能確保計算的收斂性。此外，一方面為了流體模擬的真實性，我們可能需要選取一個足夠小的時間步長來復現流體的高質量細節。另一方面，有時高性能的需求又使得我們不能選取太小的時間步長去渲染一幀。假設一幀至少要進行三個時間步的模擬，那么\(\triangle t\)應該至少設成幀間隔時間的三分之一。

3.2 網格結構

歐拉法的整個流程都是基於網格的，所以合理的網格結構是算法高效的關鍵點。Harlow和Welch提出了一種經典的MAC（marker and cell）網格結構，許多不可壓縮流體模擬的算法都在這個網格結構上呈現出了良好的效率。MAC網格是一種交叉排列的網格，不同類型的物理量被存儲於網格的不同位置。以二維的網格為例，如下圖所示，流體粒子的壓力數據存儲於網格的中心點\(P_{i,j}\)，而速度則沿着笛卡爾坐標被分成了兩部分。水平方向的\(u\)成分被存儲在了網格單元豎直邊的中心處，例如網格單元\((i,j)\)和\((i+1,j)\)之間的水平速度記為\(u_{i+1/2,j}\)。垂直方向的\(v\)成分則被存儲在了網格單元水平面的中心上。這樣的存儲方案十分有利於估算流體流進/流出某個網格單元的量。

擴展到三維的情況，MAC網格同樣是交錯排列的結構網格，如上圖所示。壓力數值存儲在立方體網格單元的中心，三個速度分量分別被記錄在立方體網格單元的三個表面的中心點上。在數值計算時，這樣的分配方式使得我們可以准確地采用中心差分法計算壓力梯度和速度的散度，同時克服了中心差分法的一個普遍的缺點。一維的情況為例，在網格頂點位置\(...,q_{i-1}，q_i,q_{i+1}...\)上估算量場\(q\)的導數，為了無偏（所謂無偏，就是不偏向左邊或者右邊）估計網格頂點i處的\(\frac{\partial q}{\partial x}\)，一種比較自然的方式就是采用一階中心差分法:

\[\left(\frac{\partial q}{\partial x}\right)_{i} \approx \frac{q_{i+1}-q_{i-1}}{2 \Delta x} \tag{3.4} \label{3.4} \]

上式\(\eqref{3.4}\)是無偏的，且精確度為\(O(\triangle x^2)\)。而前向歐拉差分法偏向右邊且精確度只有\(O(\triangle x)\)

\[\left(\frac{\partial q}{\partial x}\right)_{i} \approx \frac{q_{i+1}-q_{i}}{ \Delta x} \tag{3.5} \label{3.5} \]

然而，上式\(\eqref{3.4}\)存在着一個非常嚴重的問題：網格點\(i\)的估算導數完全忽略了\(q_i\)的值。數學上，只有常數函數的一階導數為零。但是上式遇到了鋸齒函數如\(q_i=(-1)^i\)時，它錯誤地將該類函數的導數估算為0，這種問題被稱為零空間問題（null-space problem）。

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

\[\left(\frac{\partial q}{\partial x}\right)_{i} \approx \frac{q_{i+1/2}-q_{i-1/2}}{ \Delta x} \tag{3.6} \label{3.6} \]

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

\[\begin{gathered} \vec{u}_{i, j}=\left(\frac{u_{i-1 / 2, j}+u_{i+1 / 2, j}}{2}, \frac{v_{i, j-1 / 2}+v_{i, j+1 / 2}}{2}\right) \\ \vec{u}_{i+1 / 2, j}=\left(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}\right) \\ \vec{u}_{i, j+1 / 2}=\left(\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}\right) \end{gathered} \tag{3.7} \label {3.7} \]

最后，在實現中下標索引一般沒有浮點數之說，前面直接采用\(i+1/2\)的記法是為了便於敘述。一般約定如下：

\[\begin{gathered} 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} \end{gathered} \tag{3.8} \label{3.8} \]

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

4.對流算法

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

\[\frac{Dq}{Dt}=0 \tag{4.1} \label{4.1} \]

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

\[q^{n+1}=advect(\vec{u},\triangle t,q^n) \tag{4.2} \label{4.2} \]

公式\(\eqref{4.2}\)中的各個符號含義：

• \(\vec{u}\):在MAC網格上的離散化的速度場；
• \(\triangle t\):時間步長
• \(q^n\):當前的物理量場\(q\)（如流體密度、速度、燃燒物濃度等）；
• \(q^{n-1}\):經過對流后得到的新的量場。

在這里要特別注意，輸入對流算法的速度\(\vec{u}\)必須是無散度的，否則模擬結果會出現一些奇怪的失真現象。

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

一維情況下，對流方程\(\eqref{4.1}\)寫出偏微分的形式如下：

\[\frac{\partial q}{\partial t}+u\frac{\partial q} {\partial x}=0 \tag{4.3} \label{4.3} \]

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

\[\frac{q_{i}^{n+1}-q_{i}^{n}}{\Delta t}+u_{i}^{n} \frac{q_{i+1}^{n}-q_{i-1}^{n}}{2 \Delta x}=0 \tag{4.4} \label {4.4} \]

轉成以\(q^{n+1}_i\)維計算目標的顯式公式，得：

\[q_{i}^{n+1}=q_{i}^{n}-\Delta t u_{i}^{n} \frac{q_{i+1}^{n}-q_{i-1}^{n}}{2 \Delta x} \tag{4.5} \label{4.5} \]

公式\(\eqref{4.1}\)有着非常嚴重的漏洞。

首先，前向歐拉法被證明是無條件不穩定的空間離散方法：無論取多么小Δ𝑡，隨着時間步的推進，累積誤差終將發散；然后標准中心差分方法不可避免地會出現的零空間問題，具有高頻震盪性質的速度場對空間的導數被錯誤地計算為0或幾乎為0，低離速度分量被分離出來，從而導致模擬效果中出現許多奇怪的高頻擺動和震盪。

針對這些問題，研究者們提出了一個解然不同的、更加簡單和更具物理直觀意義的半拉格朗日法。之所以叫半拉格朗日法，是因為這種方法是以拉格朗日視角去解決歐拉視角的對流方程（“半”字的由來）

假設我們的目標室求解網格\(\vec{x_g}\)的在第n+1個時間步時關於物理量\(q\)的新值，記為\(q^{n+1}_G\)。在拉格朗日的視角下，我們可以尋找在第\(n+1\)時間步之前，是空間中的哪一個點上的流體粒子在速度場\(\vec{u}\)的作用下“流向”了\(\vec{x_G}\),我們記這個粒子在第\(n\)個時間步時的網格位置為\(\vec{x_p}\),則第\(n+1\)個時間步時\(q^{n+1}_G\)即為第\(n\)個時間步時\(\vec{x_p}\)的\(q^n_p\)。如下圖為半拉格朗日對流法的示意圖

半拉格朗日對流法的第一步就是要找\(\vec{x_p}\)，為此我們根據\(\vec{x_G}\)做反向的追蹤。粒子位置對時間的導數就是速度場:

\[\frac{d\vec{x}}{d t}=\vec{u}(\vec{x}) \tag{4.6} \label{4.6} \]

經過一個時間步長\(\triangle t\)之后，粒子由\(\vec{x_p}\)移動到\(\vec{x_G}\)。為了得到\(\vec{x_p}\)，最簡單的方法就是采用前向歐拉法進行倒推：

\[\vec{x_p}=\vec{x_G}-\triangle t\vec{u}（\vec{x_G}）\tag{4.7} \label{4.7} \]

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

\[\begin{gathered} \vec{x}_{m i d}=\vec{x}_{G}-\frac{1}{2} \Delta t \vec{u}\left(\vec{x}_{G}\right) \\ \vec{x}_{P}=\vec{x}_{G}-\Delta t \vec{u}\left(\vec{x}_{m i d}\right) \end{gathered} \tag{4.7} \]

倒推得到\(\Delta t\)之前的網格位置\(\vec{x_P}\)一般不會恰好在網格頂點上，為此我們需要做些插值。三維模擬通常采用三線性插值，而二維的則采用雙線性插值。

\[q_{G}^{n+1}=\text { interpolate }\left(q_{n}, \vec{x}_{P}\right) \tag{4.8} \label{4.8} \]

4.2 邊界情況

若我們倒推得到\(\vec{x}_{P}\)仍然在流體的內部，那么做插值是完全沒問題的。但\(\vec{x}_{P}\)在流體的邊界之外呢？這種情況的出現的原因通常有兩個：一個\(\vec{x}_{P}\)確確實實在流體的外部且即將流入流體內部，另一個是由前向歐拉法或龍格庫塔法的數值計算方法帶來的誤差導致。

• 在一種情況下，我們應該知道當流體流入時其攜帶的物理量，此時我們將這個外部流入的物理量作為返回值即可。例如，第\(n\)個時間步時的外部流體以速度\(\vec{U}\)和溫度\(T\)在第\(n+1\)個時間步時注入流體內\(\vec{x}_{G}\)的位置，那\(\vec{T}_G^{n+1}\)的值就為\(T\)。
• 在第二種由誤差導致的情況下，一個適當的策略就是根據邊界上的最近點外推出所求得物理量。在模擬某些流體時，外推變得很簡單。例如，在模擬煙霧時我們簡單地假設煙霧流體外部即空氣的速度風場為某個常數\(\vec{U}\)（可能為0），這樣邊界上的速度場都取\(\vec{U}\)。但還有一些必須根據流體內部的已知量外推出未知量，這時情況就變得比較復雜了。具體如何外推將在后面介紹，目前我們只需要知道大概的步驟：首先尋找邊界上的最近點，然后在最近點的領域內插值獲取相應的物理量場。

4.3時間步長大小

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

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

\[\Delta t \le \frac{5 \Delta x}{u_{max}} \tag{4.9} \label {4.9} \]

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

\[u_{\max }=\max \left(\left|u^{n}\right|\right)+\Delta t|g| \tag{4.10} \label {4.10} \]

將不等式(4.9)(4.9)的最大值帶入公式\(\eqref{4.10}\)，我們有：

\[u_{\max }=\max \left(\left|u^{n}\right|\right)+\frac{5 \Delta x}{u_{\max }}|g| \tag{4.11} \label {4.11} \]

取一個簡單的速度上界（簡化了公式\(\eqref{4.11}\)），\(u_{max}\):

\[u_{\max }=\max \left(\left|u^{n}\right|\right)+\sqrt{5 \Delta x g} \tag{4.12} \label {4.12} \]

這樣確保了\(u_{max}\)始終為正，且避免公式\(\eqref{4.9}\)的除0錯誤。

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

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

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

\[\Delta t \leq C \frac{\Delta x}{|\vec{u}|} \tag{4.13} \label {4.13} \]

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

4.4數值耗散

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

以一維的對流項計算為例，流體速度為常量\(U>0\):

\[\frac{\partial q}{\partial t}+u \frac{\partial q}{\partial x}=0 \tag{4.14} \label {4.14} \]

假設\(\Delta t < \frac{\Delta x}{u}\),即單個時間步長內粒子追蹤軌跡小於單個網格單元的大小。我們的目標點時\(x_i\),則倒推地道道阿粒子位置就落在了\({x_{i-1},x_i}\)上的\(x_i-\Delta t u\),r然后進行線性插值得到\(q^{n+1}_i\):

\[q^{n+1}=\frac{\Delta t u}{\Delta x} q_{i-1}^{n}+\left(1-\frac{\Delta t u}{\Delta x}\right) q_{i}^{n} \tag{4.15} \label {4.15} \]

將公式\(\eqref{4.15}\)整理一下，有：

\[q_{i}^{n+1}=q_{i}^{n}-\Delta t u \frac{q_{i}^{n}-q_{i-1}^{n}}{\Delta x} \tag{4.16} \label {4.16} \]

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

\[q_{i-1}^{n}=q_{i}^{n}-\left(\frac{\partial q}{\partial x}\right)_{i}^{n} \Delta x+\left(\frac{\partial^{2} q}{\partial x^{2}}\right)_{i}^{n} \frac{\Delta x^{2}}{2}+O\left(\Delta x^{3}\right) \tag{4.17} \label {4.17} \]

將公式\(\eqref{4.17}\)代入公式\(\eqref{4.16}\)，並做一些變量消去，可得：

\[q_{i}^{n+1}=q_{i}^{n}-\Delta t u\left(\frac{\partial q}{\partial x}\right)_{i}^{n}+\Delta t u \Delta x\left(\frac{\partial^{2} q}{\partial x^{2}}\right)_{i}^{n}+O\left(\Delta x^{2}\right) \tag{4.18} \label {4.18} \]

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

\[\frac{\partial q}{\partial t}+u \frac{\partial q}{\partial x}=u \Delta x\left(\frac{\partial^{2} q}{\partial x^{2}}\right) \tag{4.19} \label {4.19} \]

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

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

\[\begin{eqnarray*} q \approx(1-s) x_{i}+s x_{i+1} \tag{4.20} \label {4.20}\\ q \approx\left[-\frac{1}{3} s+\frac{1}{2} s^{2}-\frac{1}{6} s^{3}\right] q_{i-1}+\left[1-s^{2}+\frac{1}{2}\left(s^{3}-s\right)\right] q_{i} \tag{4.21} \label {4.21}\\ +\left[s+\frac{1}{2}\left(s^{2}-s^{3}\right)\right] q_{i+1}+\left[\frac{1}{6}\left(s^{3}-s\right)\right] q_{i+2} \end{eqnarray*} \]

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

\[\begin{eqnarray*} 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}, \tag{4.22} \label {4.22}\\ 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=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} \label {4.23} \end{eqnarray*} \]

當然也可以先沿着\(y\)軸，然后再沿着\(x\)軸做插值操作。

參考資料

有限差分法（FDM）：微分方程數值求解——有限差分法

文獻：《FLUID SIMULATION SIGGRAPH 2007 Course Notes》,《GPU Games》

博客：流體模擬Fluid Simulation：流體模擬基礎

視頻：散度與旋度：麥克斯韋方程組、流體等所用到的語言