[Paper]Stable Fluids

[Shimo888]

Jos Stam. Stable Fluids, 1999

• CG分野における流体シミュレーションの元祖的な論文
• 物理ベースのシミュレーションかつ高速な計算手法
• 見た目重視の計算手法なので精度は不十分で、数値的散逸は多い（実世界より流れの減衰が急）

[Shimo888]

2.1 Stable Navier-Stokes -基本方程式-

ナビエストークス方程式

$$\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu \nabla^2 \mathbf{u} + \mathbf{f} \cdots(1)$$ $$\nabla \cdot \mathbf{u} = 0 \cdots(2)$$

Helmholtz–Hodge 分解

任意のベクトル $\mathbf{w}$は、発散が0のベクトル $\mathbf{u}$と、あるスカラー量qの勾配の和として分解できる。

$$\mathbf{w} = \mathbf{u} + \nabla q \cdots(3)$$ $$\mathbf{\nabla} \cdot \mathbf{u} = 0 \cdots(4)$$

よって、上の式の発散をもとめると

$$\nabla \cdot \mathbf{w} = \nabla^2 q \cdots(5)$$

上のポアソン方程式を解けばqが求まることがわかる。

このとき、任意のベクトル $\mathbf{w}$に対しHelmholtz–Hodge 分解をしたときのuベクトルを取り出す演算子Pを以下のように定義することができる

$$\mathbf{u} = P\mathbf{w} = \mathbf{w} - \nabla q \cdots(6)$$

ナビエストークス方程式への適用

ナビエストークス方程式に上で定義した演算子Pを作用させると、

$$P (\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u}) = P(-\frac{1}{\rho}\nabla p + \nu \nabla^2 \mathbf{u} + \mathbf{f}) \cdots(7)$$ $$\frac{\partial {P \mathbf{u}}}{\partial t} + P(\mathbf{u} \cdot \nabla)\mathbf{u} = P(-\frac{1}{\rho}\nabla p + \nu \nabla^2 \mathbf{u} + \mathbf{f}) \cdots(8)$$

ここで、 $P\mathbf{u} = \mathbf{u}$と $P \nabla p = 0$ より、

$$\frac{\partial \mathbf{u}}{\partial t} = P\Big( - (\mathbf{u} \cdot \nabla) \mathbf{u} + \nu \nabla^2 \mathbf{u} + \mathbf{f} \Big)$$

Note

$P\mathbf{u} = \mathbf{u}$の説明
ここで(式2)と(式5)より

$$\mathbf{u_u} = P\mathbf{u} = \mathbf{u} - \nabla u_q$$ $$\nabla \cdot \mathbf{u} = \nabla^2 u_q = 0$$

ここでNeumann境界条件

$$\frac{\partial u_q}{\partial n} = 0 \quad \text{on } \partial D \qquad$$

より、Neumann 境界条件の下では調和関数の解は定数だけになるので、 $u_q$は定数となり

$$P\mathbf{u} = \mathbf{u}$$

となる

Note

$P \nabla p = 0$の説明
ヘルムホルツホッジ分解は発散0の成分とgradientに一意に分解するものであり、演算子Pは発散0の成分を抽出する演算子なので、演算子Pがgradientに適用した場合結果は0になるのは自明である。

[Shimo888]

2.2 Stable Navier-Stokes -Method Of Solution-

上のコメントで導出した式をもとにシミュレーションを行う

シミュレーションは上の画像のプロセスで行う。
上のコメントの式9を基本としており、最後に演算子Pを適用することで、次ステップの流速を求める。

1. $w_0$の定義

$$\mathbf{w}_0(\mathbf{x}) = \mathbf{u}(\mathbf{x}, t) \cdots (1)$$

2. 外力の計算

$$\mathbf{w}_1(\mathbf{x}) = \mathbf{w}_0(\mathbf{x}) + \Delta t \, \mathbf{f}(\mathbf{x}, t) \cdots(2)$$

3. 移流の計算

$$\mathbf{w}_2(\mathbf{x}) = \mathbf{w}_1(\mathbf{p}(\mathbf{x}, - \Delta t)) \cdots(3)$$

4. 粘性の計算

$$(\mathbf{I} - \nu\Delta t\nabla^2) \mathbf{w}_3(\mathbf{x}) = \mathbf{w}_2(\mathbf{x}) \cdots(4)$$

5. 演算子Pによる射影

$$\mathbf{w}_4 = \mathbf{w}_3 - \nabla q \cdots (5)$$

6. 次ステップの流速

$$\mathbf{u}(\mathbf{x}, t + \Delta t) = \mathbf{w}_4(\mathbf{x}) \cdots (6)$$

[詳説] 外力の計算

$$\mathbf{w}_1(\mathbf{x}) = \mathbf{w}_0(\mathbf{x}) + \Delta t \, \mathbf{f}(\mathbf{x}, t) \cdots(2)$$

ただの離散化なのでスキップ。

[詳説] 移流の計算

$$\mathbf{w}_2(\mathbf{x}) = \mathbf{w}_1(\mathbf{p}(\mathbf{x}, - \Delta t)) \cdots(3)$$

先行の手法では有限差分的な手法を行なっていたため、時間ステップ幅が大きい場合などは計算の破綻が起こった。
この論文では、特性曲線法 (Method of Characteristics)というアプローチを採用することで上の問題を解決した。
この手法は「無条件安定性」であり、シミュレーションの破綻を防ぐことができる。
考え方はシンプルで、ある計算格子点 x の次の時刻での速度を求めたい場合、

1. その点 x に今いる流体粒子が、「Δt 時間前にどこから来たのか」を、現在の速度場を逆にたどって探す（バックトレース）。
2. 見つけ出した過去の位置での速度を、点 x の新しい速度として採用する

だけである。
$\mathbf{p}(\mathbf{x}, - \Delta t)$を $ \Delta t$秒前のxにある粒子の位置と定義して、上の式が成り立つ。
格子点以外の場所の値の補完処理だけなので、かなり簡単に実装可能である。
詳細な説明については論文の補足を参照。

[詳説] 粘性の計算

$$(\mathbf{I} - \nu\Delta t\nabla^2) \mathbf{w}_3(\mathbf{x}) = \mathbf{w}_2(\mathbf{x}) \cdots(4)$$

陽解法でとくと計算が安定しないので、陰解法で時間を進める。
陰解法で記述すると以下のようにかける( $\mathbf{w}_3$使うだけ)。

$$\mathbf{w}_3(\mathbf{x}) = \mathbf{w}_2(\mathbf{x}) + \nu\Delta t\nabla^2 \mathbf{w}_3(\mathbf{x})$$

この式を式変形すると、 $\nu\Delta t\nabla^2$は行列とみなせるので上の式4が求まる

[詳説] 演算子Pによる射影

$$\mathbf{w}_4 = \mathbf{w}_3 - \nabla q \cdots (5)$$

上のコメントの式5,6より、

$$\nabla^2 q = \nabla \cdot \mathbf{w}_3$$

とおくと、式5が成り立つのがわかる。
この式をちゃんと解くにポアソン方程式のSolverが必要になる。
ポアソン方程式ととく方法に緩和法と呼ばれる方法があるらしいが、収束が遅いのであまり使えないらしい。
計算の正確さが精度や見た目にモロに依存するので、論文ではポアソン方程式を空間的に離散化させて、粘性とPの射影のステップを大規模な疎行列で計算したらいい。 線形時間で動作するらしくO(N)らしい。
ただ、計算が複雑らしく周期境界条件のときのみ有効ではあるが、FFTを用いた簡単な実装(下のコメント)も提案している。

[Shimo888]

2.2 Stable Navier-Stokes -Periodic Boundaries and the FFT-

周期境界条件を保つ領域を考えるとアルゴリズムが非常にシンプルになるらしい。
周期性があるために速度場をフーリエ領域に変換することが可能になる（下で補足）。　
フーリエ領域における拡散演算子と射影演算子は、以下の式に変換できる。

• 拡散演算子

$$(\mathbf{I} - \nu\Delta t\nabla^2) \rightarrow (1 + \nu\Delta t k^2)$$

• 射影演算子

$$P\mathbf{w} \rightarrow \hat{P}\hat{\mathbf{w}}(\mathbf{k}) = \hat{\mathbf{w}}(\mathbf{k}) - \frac{1}{k^2}(\mathbf{k} \cdot \hat{\mathbf{w}}(\mathbf{k}))\mathbf{k}$$

計算フローは以下である。

$$\text{力の加算 (add force): } \mathbf{w}_1 = \mathbf{w}_0 + \Delta t >\mathbf{f}$$ $$\text{移流 (advect): } \mathbf{w}_2(\mathbf{x}) = \mathbf{w}_1(\mathbf{p}(\mathbf{x}, -\Delta t))$$ $$\text{フーリエ変換 (transform): } \hat{\mathbf{w}}_2 = \text{FFT}\{\mathbf{w}_2\}$$ $$\text{拡散 (diffuse): } \hat{\mathbf{w}}_3(\mathbf{k}) = \frac{\hat{\mathbf{w}}_2(\mathbf{k})}{1 + \nu\Delta t k^2}$$ $$\text{射影 (project): } \hat{\mathbf{w}}_4 = \hat{P} \hat{\mathbf{w}}_3$$ $$\text{逆フーリエ変換 (transform): } \mathbf{w}_4 = \text{FFT}^{-1}\{\hat{\mathbf{w}}_4\}$$

FFTの計算量はO(N Log N)なので、マルチグリッド法(計算量O(N))より計算量は大きいがめっちゃシンプルに実装できるらしい。
必要なのはFFTのライブラリくらい。

[Shimo888]

2.4 Stable Navier-Stokes -Moving Substances through the Fluid-

流体中を拡散して移動する物質(スカラー場)の式の説明。
イメージとしては、コーヒーに入れたミルクの模様やタバコから立ち上る煙などのこと。　
ここではそれらの物質を流体中を移動する任意のスカラー量 $a$とする。
$a$は、流体中を拡散する粒子の密度や、流体の温度、テクスチャ座標？（これはよくわからん）に対応する。
それらの式は以下の移流拡散方程式で記述できる。

$$\frac{\partial a}{\partial t} = -(\mathbf{u} \cdot \nabla)a + \kappa_a \nabla^2 a - \alpha_a a + S_a$$

このとき、$\kappa_{a}$は拡散系数、$a_{a}$は減衰率、$S_{a}$はソース項（生み出される量）である。
イメージ的には、

• 流体に流される時に働く移流項
• 粒子の拡散を表す拡散項
• 減衰項（なぜ自身の量に比例するのかは？、熱とかだとそうなる？）
• 粒子の蛇口であるソース項

からなる式で直感と一致する。

減衰項の計算

$$(1 + \Delta t \alpha_a) a(\mathbf{x}, t + \Delta t) = a(\mathbf{x}, t)$$

陰解法で記述すると、以下のようにかける。

$$a(\mathbf{x}, t + \Delta t) = a(\mathbf{x}, t) - \Delta t \alpha_a a(\mathbf{x}, t + \Delta t)$$

式変形すると上の式になる。

[Shimo888]

プログラムの説明

全体

while ( simulating ) {
    /* 表示とユーザーインタラクションの処理 */

    /* UIから力FとソースSsourceを取得 */
    Swap(U1, U0); 
    Swap(S1, S0);

    Vstep( U1, U0, visc, F, dt );
    Sstep( S1, S0, kS, aS, U1, Ssource, dt );
}

• U0は前ステップの流速
• U1は次ステップの流速
• S0は前ステップのスカラー場
• S1は次ステップのスカラー場
• dt は時間刻み幅
• viscは粘性系数
• Fは外力
• ks は拡散係数
• aSは減衰率
• Sourceはソース

速度のSolver

Vstep ( U1, U0, visc, F, dt )
    for(i=0; i<NDIM; i++)
        addForce ( U0[i], F[i], dt );

    for(i=0; i<NDIM; i++)
        Transport ( U1[i], U0[i], U0, dt );

    for(i=0; i<NDIM; i++)
        Diffuse ( U0[i], U1[i], visc, dt );
        
    Project ( U1, U0, dt );

スカラー場のSolver

Sstep ( S1, S0, k, a, U, source, dt )
    addForce ( S0, source, dt );
    Transport ( S1, S0, U, dt );
    Diffuse ( S0, S1, k, dt );
    Dissipate ( S1, S0, a, dt );

外力

数式から自明
入力に時間ステップを加算したものを足すだけ

減衰

数式(前コメントの式)から自明
1 + dt * a でわり加算するだけ。

移流

Transport ( S1, S0, U, dt )
    for each cell (i,j,k) do
        // セルの中心位置を計算
        X = O + (i+0.5, j+0.5, k+0.5) * D; 
        
        // 粒子を時間-dtだけ逆向きに追跡
        TraceParticle ( X, U, -dt, X0 );
        
        // 過去の位置X0で場の値を線形補間し、新しい値とする
        S1[i,j,k] = LinInterp ( X0, S0 );
    end

拡散と射影

周期境界ではないときは、そ行列のソルバーが必要
マルチグリッド法などが理論的には最適らしい。
todo: 補足に詳しい説明がされているので後でまとめる