Stable Fluids

Jos StamさんのStable Fluids (1999)という論文についてのまとめ

シミュレーションの流れは以下の画像で、速度場wをさまざまなステップを経て更新してくことが、シミュレーションの実現に直結することがわかる

コード上では、vel_stepという関数にて、流体シミュレーションの本質である速度場を更新していく全体的なステップを記述している。

                    /* 速度更新ステップ (全体的なステップ) */
/* N: グリッドの一辺 
 * (u, v): 次の時間ステップの速度x, y成分, 
 * (u0, v0): 前の時間ステップの速度x, y成分 
*/
void Simulation::vel_step(
    int N, std::vector<float>& u, std::vector<float>& v, 
        std::vector<float>& u0, std::vector<float>& v0, float visc, float dt)
{
    /* add_forceで更新されたu0, v0を次の時間ステップのu, vに反映*/
	add_source(N, u, u0, dt);
	add_source(N, v, v0, dt);

    advect(N, 1, u, u0, u0, v0, dt); 
	advect(N, 2, v, v0, u0, v0, dt);

	std::swap(u0, u); 
	diffuse(N, 1, u, u0, visc, dt);

	std::swap(v0, v);
	diffuse(N, 2, v, v0, visc, dt);

	project(N, u, v, u0, v0);

	std::swap(u0, u);
	std::swap(v0, v);

	project(N, u, v, u0, v0);
}

ステップ1: 外力項の扱い(add force)

最初に Δt の間に外力は大きく変化しないと仮定し、外力の効果を加える
$\mathbf{w}_1(x) = \mathbf{w}_0(x) + \Delta t \mathbf{f}$

マウスクリックされたセルに対して外力適応

            /* ステップ1 外力項の加算（クリックしたセルに対して外力を適応) */
/* X, Y -> クリックした座標、
*  N -> サイズをスケールで割ったもの、
*  u, v -> マウスの動く方向に力の大きさ(force)をかけたもの
*/
void Simulation::add_force(int X, int Y, int N, float u, float v)
{
	if (X <= 0 || X > N || Y <= 0 || Y > N) {
		throw std::out_of_range("Index out of range.");
	}
    /* クリックしたセルにおいて、速度フィールドの直前の時間ステップの値に外力u, vを代入 */
	x_prev[IX(Y, X)] = u;
	y_prev[IX(Y, X)] = v;
    /* 次の時間ステップでsource項として機能 */
}

概要

主に、ユーザのインタラクション（マウスクリック・ドラッグ）による「局所的なまたは一点」に対して速度場(u, v)に外力を加えるための関数であり、ここでは特定の１セルに対して速度に該当する成分(u, v)を"前ステップの速度場"を表すx_prevやy_prevに代入している。x_prevやy_prevはadd_source関数を通じてx, yに反映される「入力バッファ」のような役割を持っている。

引数の説明

• x, y : クリックした座標
• N : サイズをスケールで割ったもの（グリッドの一辺）
• u, v : マウスの動く方向に力の大きさ（force）を掛けたもの。速度

関数の中身

• if文: グリッドの外は範囲外であると指定
• マウスクリックがあったセルに対してのみ外力を設定
• x_prev と y_prev は次の時間ステップでの速度場の更新に使用

全てのセルに外力の影響を適応

/* ステップ1: 外力項の加算
    (全てのセルに適応: フィールド全体に対して、外部からの影響を時間ステップに基づいて加算) */
void Simulation::add_source(int N, std::vector<float>& x,
                                std::vector<float>& s, float dt)
{
	for (int i = 0; i < size; i++) {
		x[i] += dt * s[i];
	}
}

概要

xが示す物理量（速度場、密度場など）に、sが表す外部から供給されるソース（源項）を時間ステップdtに基づいて「全体的」に加算する。sはadd_force関数やstampによって設定済みのx_prev(またはy_prev)配列が格納されている。つまり、add_source関数ではsを積分として反映しているというイメージになる。

引数の説明

• N : サイズをスケールで割ったもの（グリッドの一辺）
• x : 位置
• s : x_prevなどの前の時間ステップの速度
• dt : 微小時間Δt

関数の中身

• for文 : グリットの数(size = N * N)だけ回す
• x[i] += dt * s[i]; : 時間 * 速度によって位置を求めている

ステップ2: 移流項、対流項の扱い（advect）

1. 移流項
   • 流体の移動による物体や量が運ばれる現象
   • 非線形方程式であり、直接解くのが難しい
   • ナビエストークス方程式の以下の項に値する
     $(u \cdot \nabla ) u$
   • この論文では特性曲線法を使用
   • 一般的には有限差分を使用
   
   
2. 特性曲線法
   • 変微分方程式の数値解法の一つで、オイラー座標系で実装
   • 一つ前の時間ステップの粒子の動きを追跡することで、移流の効果を計算
   • 安定性が高い

                /* ステップ2: 移流項、対流項 (Advection Term) */
/* N: グリッドの一辺の長さ, b: 境界条件指定パラメータ, 
 * d: 現在の移流後の値(次のステップでの密度や速度場), 
 * d0: 一つ前の時間ステップでの密度や速度場,
 * (u, v): xy成分の速度, dt: 時間ステップの大きさ
*/
void Simulation::advect(
    int N, int b, std::vector<float>& d, std::vector<float>& d0,
        std::vector<float>& u, std::vector<float>& v, float dt)
{
    int i, j, i0, j0, i1, j1;
        float x, y, s0, t0, s1, t1, dt0;
        dt0 = dt * N;  // 時間ステップとグリッドサイズに基づくスケーリング係数
        for (i = 1; i <= N; i++) {
            for (j = 1; j <= N; j++) {  // グリッド内の全てのセルに対して行う
                x = i - dt0 * u[IX(i, j)];  // x方向の移流後の位置を逆にたどる
                y = j - dt0 * v[IX(i, j)];  // y方向の移流後の位置を逆にたどる
                
                // xとyの範囲を0.5からN+0.5にクリップして境界外に出ないようにする
                if (x < 0.5) x = 0.5;
                if (x > N + 0.5) x = N + 0.5;
                i0 = (int)x;  // 移流元の整数部分のインデックス
                i1 = i0 + 1;  // 隣接するインデックス（補間に使用）
                
                if (y < 0.5) y = 0.5;
                if (y > N + 0.5) y = N + 0.5;
                j0 = (int)y;  // 移流元の整数部分のインデックス
                j1 = j0 + 1;  // 隣接するインデックス（補間に使用）

                s1 = x - i0;  // x方向の補間比率
                s0 = 1 - s1;  // x方向の逆補間比率
                t1 = y - j0;  // y方向の補間比率
                t0 = 1 - t1;  // y方向の逆補間比率

                // 移流後の値を計算し、補間によって求める
                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)]);
            }
        }
        set_bnd(N, b, d);  // 境界条件を設定
}

関数の中身

1. スケーリング係数 dt0

dt0 = dt * N;  // 時間ステップとグリッドサイズに基づくスケーリング係数

流体の移流を安定かつ正確にシミュレーションするためのスケーリング係数。dt だけではなく N も掛ける理由は、速度場 u と v のスケールがグリッドサイズに依存しているため。

• 移流項では、流体の粒子がどれだけの距離を移動するかを計算するため、速度uやvの影響を適切にスケーリングする必要がある
• グリッドが大きい（セル数が多い）場合、各セルのサイズが小さくなる。これにより、速度の影響を大きくした方が正確な移動距離を反映でできる

つまり、dt0 を導入することで、速度場に基づく移流がセルサイズに応じて適切にスケーリングされ、より安定で正確なシミュレーションが可能になる。

2. x, y について

x = i - dt0 * u[IX(i, j)];
y = j - dt0 * v[IX(i, j)];  

• i と j は現在のグリッド位置
• x と y一つ前の時間ステップにおける位置
• u[IX(i, j)]とv[IX(i, j)]は、それぞれ現在位置(i, j)でのx方向とy方向の速度

このコードは、現在位置(i, j)にある流体がどこから来たのかを求めている。速度ベクトル(u, v)にdt0を掛けて逆向きにすることで、一つ前の時間ステップにおける位置(x, y)を計算している。これは「逆追跡」と呼ばれる手法で、次の時間ステップの値を計算する際に使われる。

3. 境界の範囲内にクリッピング

if (x < 0.5) x = 0.5;
if (x > N + 0.5) x = N + 0.5

if (y < 0.5) y = 0.5;
if (y > N + 0.5) y = N + 0.5;

クリッピングとは、ある値が指定された範囲外に出てしまうのを防ぎ、範囲内に収める処理のこと。プログラムでは、変数が想定外の範囲に出ないように上限や下限を設定することで、計算結果が異常にならないよう使用される。

ここでは x, y それぞれの最小値が 0.5 で、最大値が N + 0.5 であると指定している。N はグリッドサイズであり、セル全体の集合すなわち、シミュレーション領域を表す。

4. 補間

「移流元の整数部分のインデックス」と「隣接するインデックス」

i0 = (int)x;  // 移流元の整数部分のインデックス
i1 = i0 + 1;  // 隣接するインデックス（補間に使用）

j0 = (int)y;  // 移流元の整数部分のインデックス
j1 = j0 + 1;  // 隣接するインデックス（補間に使用）

現在の位置から逆方向に辿った一つ前の時間ステップでの移流元の位置を、整数値で近似したセルのインデックスのことである。

移流元の位置 (x, y) は一般的に少数を含む座標となり、グリッド上のセルの中央や辺にピッタリ位置するわけではない。そのため、(x, y)を整数セルのインデックスに変換して扱う必要がある。

i0 = (int)x , j0 = (int)y は x, y の整数部分で、x, y の左側または上側に位置するセルを示す。

i1 = i0 + 1 , j1 = j0 + 1 は i0 , j0 に隣接するセルのインデックスで、x, yの右側または下側に位置するセルを指している。

補間比率
以下は、移流元の位置 (x, y) を整数インデックス (i0, j0) とその隣接セルに分解し、周囲のセルの値を用いて補間するための比率（補間比率）を計算している。

s1 = x - i0;  // x方向の補間比率
s0 = 1 - s1;  // x方向の逆補間比率
t1 = y - j0;  // y方向の補間比率
t0 = 1 - t1;  // y方向の逆補間比率

1. s1 = x - i0; // x方向の補間比率
   • s1 は x の小数部分に相当し、x の位置が i0 , i1 のどちらに近いかを示すための補間比率である
   • 例えば、 x = 3.7 の場合、 i0 = 3 であり、 s1 = 0.7 となる。この値は x が i1 = 4 に 70% の比率で近いことを意味する
   
   
2. s0 = 1 - s1; // x方向の逆補間比率
   • s0 は s1 の逆比率で、i0 にどれだけ近いかを示している
   • s1 が 0.7 ならば、s0 は 0.3 となり、x が i0 = 3 に 30% の比率で近いことを表す
   • この s0, s1 によって、 x が i0 と i1 の間でどの位置にあるかを示すことができ、周囲のセルの値を基に x 方向の線型補間が可能になる
   
   
3. t1 = y - j0; // y方向の補間比率
   • t0 は y の小数部分に相当し、y の位置が j0 , j1 のどちらに近いかを示すための補間比率である
   • 例えば、 y = 2.3 の場合、 j0 = 2 であり、 t1 = 0.3 となる。この値は y が j1 = 3 に 30% の比率で近いことを意味する
   
   
4. t0 = 1 - t1; // y方向の逆補間比率
   • t0 は t1 の逆比率で、j0 にどれだけ近いかを示している
   • t1 が 0.3 ならば、t0 は 0.7 となり、y が j0 = 2 に 70% の比率で近いことを表す
   • この t0, t1 によって、 y が j0 と j1 の間でどの位置にあるかを示すことができ、周囲のセルの値を基に y 方向の線型補間が可能になる

移流後の値を計算

// 移流後の値を計算し、周囲4つのセルからの補間によって求める
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)]);

1. 各要素の意味

   • d[IX(i, j)]
     • 現在の位置 (i, j) での値（例えば密度や速度）
     • IX(i, j) は (i, j) を１次元配列 d のインデックスに変換するための関数
   • d0[IX(i0, j0)], d0[IX(i0, j1)], d0[IX(i1, j0)], d0[IX(i1, j1)]
     • 移流元の周囲にある4つのセル (i0, j0), (i0, j1), (i1, j0), (i1, j1) のそれぞれの値を表している
   • s0, s1, t0, t1
     • 補間比率

2. 補間の計算手順
   ここでは、まず y 方向の補間を行い、その結果に対して x 方向の補間を行っている。
   それぞれの方向で内分点を求めている、と考えることもできる。

   1. y 方向の補間
   • t0 * d0[IX(i0, j0)] + t1 * d0[IX(i0, j1)] は i0 列における y 方向の補間
     • t0 は d0[IX(i0, j0)] の重み（上側セルに近いほど大きくなる）
     • t1 は d0[IX(i0, j1)] の重み（下側セルに近いほど大きくなる)
   • t0 * d0[IX(i1, j0)] + t1 * d0[IX(i1, j1)] は i1 列における y 方向の補間
   1. x 方向の補間
   • s0 は i0 側の重み（左側に近いほど大きくなる）
   • s1 は i1 側の重み（右側に近いほど大きくなる）

5. 境界条件の設定

以下で境界条件を設定する

set_bnd(N, b, d);  // 境界条件を設定

ステップ3: 粘性項の扱い（diffuse）

ナビエストークス方程式は以下のような形で表される。

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

粘性項は $\nu \nabla^2 \mathbf{u}$ の部分であり、流体が内部摩擦によって速度場がスムーズになっていく効果（速度分布が平滑化される拡散現象）を示していて、数値的にはこれは「拡散方程式」として扱われる。コードは以下。

                /* ステップ3: 粘性項の扱い（拡散方程式） */
// N: グリッドの一辺, b: 境界条件を指定するパラメータ,x: 計算後の値,
// x0: 直前の時間ステップの値, diff: 拡散係数(粘性係数), dt: 時間ステップ
void Simulation::diffuse(int N, int b, std::vector<float>& x,
                            std::vector<float>& x0, float diff, float dt)
{
	int i, j, k;
     // 粘性係数ν, Δt, 1 / Δx^2 = 1 / Δy^2 = 1/h^2をaでまとめたもの
	float a = dt * diff * N * N;
	for (k = 0; k < 25; 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);
			}
		}
        /* 境界条件の適応 */
        /*  b = 0: 境界で値をそのまま内部の値に等しくさせる。　
         *  b = 1, 2: 境界での速度成分を反転（壁での反射をシミュレート）*/
		set_bnd(N, b, x);
	}
}

概要

ここでは粘性項 $\nu \nabla^2 \mathbf{u}$を離散化することから始まる。∇すなわち2回微分をテイラー展開を用いながら、拡散方程式 $ \frac{\partial x}{\partial t} = \nu \nabla^2 x$ をn+1（陰的解法）で離散化すると求めるxと前の時間ステップのxについて以下の関係が成り立つ。

$ \chi_{i,j}^{n+1} = \left( \chi_{i,j}^{n} + a \left( \chi_{i-1,j}^{n+1} + \chi_{i+1,j}^{n+1} + \chi_{i,j-1}^{n+1} + \chi_{i,j+1}^{n+1} \right) \right) \cdot \frac{1}{1 + 4a} $

$ (a = \frac{\Delta t}{h^2} ν とおく) $

i, jは座標（グリットの位置）を表していている。式をみると、(i, j)のxを求めるためには、その上下左右の座標も必要になることがわかる。

引数の説明

• N : サイズをスケールで割ったもの（グリッドの一辺）
• b : 境界条件を指定するパラメータ
• x : 求めるの値（速度など）
• x0 : 前の時間ステップの値（速度）
• diff : 拡散係数（粘性係数νに相当）
• dt : 微小時間（時間ステップ）

関数の中身

• float a = dt * diff * N * N; : 式に出てくる定数をaとしてまとめる
• ガウス・ザイデルの反復 : 陰的解法を用いると、for文を回すだけでガウス・ザイデルが適応される。反復回数は25回
• 拡散方程式 : 上記の式で離散化した拡散方程式を用い、2重ループで全てのグリットに対して、xを求める
• set_bnd(N, b, x) : 境界条件を適応する

ステップ4: 投影（project）

ニュートン流体では非圧縮($\nabla \cdot \mathbf{u} = 0$）である必要があり、投影ステップでは速度場を発散のない（非圧縮性）場に変換するためのステップである。この操作はヘルムホルツ・ホッジ分解を基に行われている。

$\mathbf{u} = \mathbf{u}_{div\text{-}free} + \nabla p$

ここで、$\mathbf{u}_{div\text{-}free}$は、発散のない速度場、$p$はスカラー圧力場を示す。投影ステップでは、速度場の発散成分を除去することで非圧縮性を保証している。

/* ステップ4: 投影ステップ (Projection Step) */
/* N: グリッドの一辺、 (u, v): 次の時間ステップの速度x, y成分、 p: 圧力場、 div: 速度場の発散*/
void Simulation::project(int N, std::vector<float>& u, std::vector<float>& v,
                            std::vector<float>& p, std::vector<float>& div)
{
	int i, j, k;
	float h;
	h = 1.0 / N;    // グリッドの単位長さ
    
    // 発散場を計算
	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;    // 圧力場を初期化
		}
	}
	set_bnd(N, 0, div); set_bnd(N, 0, p);   // 発散場と圧力場に対して境界条件を適用
    
    // ポアソン方程式をガウス・ザイデル法で反復的に解く
	for (k = 0; k < 20; k++) {  // 反復回数20回
		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);   // 境界条件を適用
	}
    
    // 圧力場の勾配を引くことで速度場を非圧縮性にする。
	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;
		}
	}
	set_bnd(N, 1, u);   // 速度場 u の境界条件を適用
	set_bnd(N, 2, v);   // 速度場 v の境界条件を適用
}

概要

ヘルムホルツ・ホッジ分解を数値的に実装するために、

1. 現在の速度場から発散場を計算
2. 発散場に対するポアソン方程式を解いて圧力場を求める
3. 圧力場の勾配を速度場から引くことで発散のない速度場を得る

というステップを踏んでいる。

引数の説明

• N: サイズをスケールで割ったもの（グリッドの一辺)
• (u, v): 次の時間ステップの速度(x, y)成分
• p: 圧力場
• div: 速度場の発散

関数の中身

1. float h = 1.0 / N;: グリッドセルの単位長さを計算、これは後で速度さを長さで割る際に使用する
2. $\nabla \cdot \mathbf{v} = 0$（速度場の発散）を計算
   • 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)]);:
   • $\nabla \cdot \mathbf{v} = 0$（速度場の発散）を中心差分法で計算し、スケーリングを行う
   • 中央差分とは、前方差分のテイラー展開と後方差分のテイラー展開の差を取り、微分を求めると以下のようになるのが、この式になる理由である
     $f'(x) = \frac{f(x + h) - f(x - h)}{2h} + O(h^2)$
3. ポアソン方程式を解く
   発散しない速度場を求めるためには
   $\nabla \cdot \mathbf{u} = 0$
   を解く必要があるが、このままは扱いづらいので、ホルムヘルツ・ホッジ分解によって分解した
   $\mathbf{u} = \mathbf{u}_{div\text{-}free} + \nabla p$
   に対して右側から発散演算子∇との内積をとると、
   $\nabla \cdot \mathbf{u} - \nabla \cdot (\nabla p) = 0$
   となり、
   $\nabla^2 p = \nabla \cdot \mathbf{u} $
   というポアソン方程式を得られ、$\nabla^2 p$を解けば、発散しない速度場を求めることにつながることがわかる。これはステップ３の粘性項にて拡散方程式を取り扱った時と同じように離散化し、ガウス・ザイデル法による反復法を適用することで解くことができる。
   • 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;
4. 圧力場の勾配を引き、速度場の発散を除去する
   • x, y成分の速度から3で求めた圧力場の勾配を引くことで発散を除去する
   • x成分: u[IX(i, j)] -= 0.5 * (p[IX(i + 1, j)] - p[IX(i - 1, j)]) / h;
   • y成分: v[IX(i, j)] -= 0.5 * (p[IX(i, j + 1)] - p[IX(i, j - 1)]) / h;