流体シミュレーションを実装してみたいその２（シミュレーション編）

はじめに

こんにちは。株式会社スパーククリエイティブCTOの広本です。
普段は SPARK GEARの開発をしたり、UEでエンジンの改造をしたり、Unityでレンダラーや揺れもの物理を作ったりしています。

前回はレンダリング周りの技術について解説しましたが、今回はその裏側で動いている心臓部、流体シミュレーションの実装についてを雑に解説していきます。

単に教科書通りのナビエ・ストークス方程式を解くだけでは、ゲームや映像で使えるような「派手な爆発」や「カッコいい煙」にはなりません。
今回は、基本実装からスタートし、マコーマック法や渦度閉じ込めといった技術をどのように組み込んで理想の挙動に近づけていったのか、そのアルゴリズムを解説します。

数式がたくさん出てきますが、私は中学数学までしか分からないので何が書いてあるかは全くわかりません。
しかし、分からないことは分かる人に教えてもらうという事で何とかなるものです。

はじめに
1. ナビエ・ストークス方程式
2. マコーマック法 (MacCormack Method)
3. 渦度閉じ込め（Vorticity Confinement）
まとめ

1. ナビエ・ストークス方程式

まずは基礎となる流体シミュレーションの実装です。
今回は非圧縮性流体を対象とし、以下のナビエ・ストークス方程式をグリッドベースで解きます。

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

$\displaystyle \nabla \cdot \mathbf{u} = 0$

実装の手順は、Jos Stam氏の "Stable Fluids" に基づく標準的なものです。

Advection (移流): 速度場に従って物理量を運ぶ
Diffuse (拡散): 粘性項
Add Force (外力): 浮力や重力を加算
Projection (圧力計算): 発散なしの状態へ補正

これをセミラグランジュ法で実装した結果がこちらです。

動いてはいますが、「煙っぽくない」のが一目瞭然です。
全体的にボヤけており、煙特有の細かい乱流や巻き込みが消失しています。
これは、セミラグランジュ法における線形補間がローパスフィルタとして作用してしまう事が原因です。
解像度を上げても、この「ヌルヌル感」は容易には消えません。

2. マコーマック法 (MacCormack Method)

ここで、流体シミュレーションにおける「移流」の精度を劇的に向上させるテクニック、マコーマック法について解説します。

「作った流体がなんだかボヤける」「解像度を上げてもヌルヌル感が消えない」といった悩みを持つ方、その原因は数値散逸にあります。
そして、その特効薬こそがマコーマック法です。

兄弟的な手法に BFECC もありますが マコーマック法 の方が実装しやすいためこちらを選択します。

1. なぜ流体はボヤけるのか？

流体シミュレーションの基本となるセミラグランジュ法は、非常に安定しており、計算が発散しにくいという素晴らしい特徴があります。

その処理はシンプルで、「ある位置 $x$ に到達する粒子は、 $\Delta t$ 秒前には $x - \mathbf{u}\Delta t$ にいたはずだ」と考え、その位置の値を線形補間で取得するだけです。

$\phi^{n+1}(\mathbf{x}) = \phi^n(\mathbf{x} - \mathbf{u}(\mathbf{x})\Delta t)$

しかし、この「線形補間」が原因です。
線形補間は、隣接するセルの値を混ぜ合わせる処理なので、信号処理的に見ると「ローパスフィルタ（ぼかし）」をかけているのと全く同じことになります。

シミュレーションではこれを毎フレーム、何千回と繰り返します。
その結果、少しずつ、しかし確実にディテールが削ぎ落とされ、全体が平均化されてボヤけてしまうのです。
これが数値散逸です。

2. マコーマック法の基本概念

マコーマック法を一言で言えば、「行って戻って、誤差を測って補正する」手法です。

通常の移流（行き）で発生した誤差を推定するために、逆方向の移流（帰り）を行い、「元の場所に正しく戻れたか？」を確認します。
戻れなかった分のズレこそが、移流によって生じた誤差であると仮定して補正を加えます。

アルゴリズムの3ステップ

手順は以下のようになります。

Forward Advection（予備段階） 通常通り、現在の速度場で移流させます。
Backward Advection（逆移流） 移動先から、逆ベクトルで移流させて戻ってきます。
理想的な世界では、行って帰ってきたら元の値に戻るはずですが、数値散逸によってズレが生じています。
Error Correction（誤差補正） 「元の値」と「行って帰ってきた値」の差が誤差です。これを半分の重みで補正に使います。

これだけで、1次精度だったセミラグランジュ法が、2次精度相当へと引き上げられます。

3. 「オーバーシュート」問題

「これで完璧だ！」と思い実装すると、画面上に謎のノイズやアーティファクトが発生することがあります。
これは、マコーマック法の補正によって「本来ありえない値」が生まれてしまう現象、オーバーシュート（およびアンダーシュート）が原因です。

補正項が大きすぎる場合、計算結果が、周囲のセルの最大値よりも大きくなったり、最小値よりも小さくなったりしてしまいます。
物理的に、移流だけで値が増減することはありません（質量保存）。

この「はみ出し」が蓄積すると、シミュレーションが不安定になり、最終的に崩壊します。

4. 「リミッター」の実装

この暴走を防ぐために必須となるのが、リミッター処理です。
「補正後の値は、移流元の近傍セルの最大値・最小値の範囲内に収まっていなければならない」という制約を強制的にかけます。

具体的には、Forward移流でサンプリングした位置（ $\mathbf{x} - \mathbf{u}\Delta t$ ）の周囲 $2 \times 2$ （または $3 \times 3$ ）のセルの最小値（Min）と最大値（Max）を求め、最終的な値をこの範囲に clamp します。

$\phi^{final} = \text{clamp}(\phi^{corrected}, \phi_{min}, \phi_{max})$

これにより、2次精度のシャープさを保ちつつ、絶対的な安定性を保証することができます。

5. HLSLによる実装

というわけで、Compute Shader (HLSL) での実装例です。

Texture2D<float4> _SourceTex; // 移流対象（密度や速度）
Texture2D<float4> _VelocityTex; // 速度場
RWTexture2D<float4> _DestTex; // 書き込み先
SamplerState _LinearClamp;
float _DeltaTime;
float2 _TexelSize; // 1.0 / TextureResolution

[numthreads(8, 8, 1)]
void MacCormackAdvect(uint3 id : SV_DispatchThreadID)
{
    float2 uv = (id.xy + 0.5) * _TexelSize;

    // 1. 速度ベクトルの取得
    float2 velocity = _VelocityTex.SampleLevel(_LinearClamp, uv, 0).xy;
    float2 displacement = velocity * _DeltaTime;

    // --- マコーマック法 Step 1: Forward ---
    // 過去の位置 (p - dt * v) を計算
    float2 fwdUV = uv - displacement * _TexelSize; // UV空間での移動量に変換
    float4 phi_n_plus_1 = _SourceTex.SampleLevel(_LinearClamp, fwdUV, 0);

    // --- マコーマック法 Step 2: Backward ---
    // さらに過去から現在へ戻る (p + dt * v)
    // ※fwdUVを起点に逆ベクトルで移動
    float2 bwdUV = fwdUV + displacement * _TexelSize; 
    float4 phi_n_bwd = _SourceTex.SampleLevel(_LinearClamp, bwdUV, 0);

    // --- マコーマック法 Step 3: Correction ---
    float4 phi_n = _SourceTex.SampleLevel(_LinearClamp, uv, 0); // 現在の位置の元の値
    float4 result = phi_n_plus_1 + 0.5 * (phi_n - phi_n_bwd);

    // --- Limiter (Clamp) 処理 ---
    // サンプリング点(fwdUV)の周囲4テクセルを取得し、Min/Maxを求める
    float2 st = fwdUV / _TexelSize - 0.5;
    float2 iuv = floor(st);
    float2 fuv = frac(st);

    // 4近傍のフェッチ (Point Sampling的に取得)
    float4 tl = _SourceTex[int2(iuv) + int2(0, 1)];
    float4 tr = _SourceTex[int2(iuv) + int2(1, 1)];
    float4 bl = _SourceTex[int2(iuv) + int2(0, 0)];
    float4 br = _SourceTex[int2(iuv) + int2(1, 0)];

    // 最小値と最大値を計算
    float4 minVal = min(min(tl, tr), min(bl, br));
    float4 maxVal = max(max(tl, tr), max(bl, br));

    // 結果を範囲内に収める
    result = clamp(result, minVal, maxVal);

    _DestTex[id.xy] = result;
}

マコーマック法を導入することで、以下の効果が得られます。

ディテールの維持: 煙の細かい渦や切れ込みが長時間消えずに残る。  
乱流の保存: 速度場自体にも適用することで、エネルギー減衰が緩やかになり、動きがダイナミックになる。

計算コストは、テクスチャ参照回数が1回から3回（リミッター含めるとそれ以上）に増えるため重くなりますが、解像度を2倍にするよりはずっと安上がりで効果的です。

「流体がヌルいな」と思ったら、まずはこのマコーマック法を試してみてください。

だらだらと解説しましたが実装した結果になります。


▲ 適用前	▲ 適用後（微細な乱流が持続する）

最初の実装に比べると明らかに細かいディティールが保存されているのが見てとれます。

3. 渦度閉じ込め（Vorticity Confinement）

「マコーマック法」によって、移流ボケを大幅に改善することに成功しました。
しかし、流体シミュレーションにはもう一つ、問題が存在します。

それが「回転エネルギーの減衰」です。

ナビエ・ストークス方程式を離散化して解く過程で、本来保存されるべき小さな渦が、グリッドの粗さと数値粘性によって急速に消滅してしまいます。
結果として、流体はヌルヌルとした層流のような動きになりがちです。

今回は、この失われた乱流エネルギーを物理的な根拠に基づいて人工的に再注入するブースト技術、渦度閉じ込めについて、その数理的背景からコンピュートシェーダによる実装までを徹底解説します。

1. 渦度 (Vorticity) とは何か

まず、流体の「回転」を定量化する物理量である渦度 $\boldsymbol{\omega}$ を定義します。渦度は、速度ベクトル場 $\mathbf{u}$ の回転（Curl）として定義されます。

$\boldsymbol{\omega} = \nabla \times \mathbf{u}$

3次元デカルト座標系において、速度 $\mathbf{u} = (u, v, w)$ であるとき、渦度 $\boldsymbol{\omega} = (\omega_x, \omega_y, \omega_z)$ の各成分は以下の偏微分で表されます。

$\omega_x = \frac{\partial w}{\partial y} - \frac{\partial v}{\partial z}, \quad \omega_y = \frac{\partial u}{\partial z} - \frac{\partial w}{\partial x}, \quad \omega_z = \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y}$

2次元シミュレーションの場合、渦度はスカラー量（Z軸周りの回転成分 $\omega_z$ のみ）として扱われることが多いですが、本質的にはベクトル量らしいです。
良くわかりませんがComputeShaderで実装する分には特に気にしなくていい気がするので実装しましょう。

2. 閉じ込め力の方程式

渦度閉じ込めは、Steinhoff and Underhill (1994) によって提唱され、その後 Fedkiw et al. (2001) によってグラフィックス向けに定式化されたそうです。

その核心は、「失われた渦度を補うような外力 $\mathbf{f}_{conf}$ を速度場に加える」という事みたいです。

ステップ1: 渦の中心方向を探る

まず、渦度の大きさ（Magnitude） $|\boldsymbol{\omega}|$ の勾配を計算し、それを正規化します。

$\mathbf{N} = \frac{\nabla |\boldsymbol{\omega}|}{|\nabla |\boldsymbol{\omega}||}$

この正規化ベクトル $\mathbf{N}$ は、渦度の強さが変化する方向を向いています。一般に、渦の中心に向かうほど渦度は強くなるため、このベクトルは渦の中心（尾根）を指し示す法線ベクトルとして機能します。

ステップ2: 外力の生成

次に、このベクトル $\mathbf{N}$ と渦度 $\boldsymbol{\omega}$ の外積を取り、係数 $\epsilon$ を掛け合わせたものを「閉じ込め力」とします。

$\mathbf{f}_{conf} = \epsilon h (\mathbf{N} \times \boldsymbol{\omega})$

ここで $h$ はグリッド幅、 $\epsilon$ は強さを調整する係数です。
外積 $\mathbf{N} \times \boldsymbol{\omega}$ は、渦の中心方向と回転軸の両方に直交する方向、すなわち「回転の接線方向」を向きます。

つまり、この力 $\mathbf{f}_{conf}$ は、既存の回転をさらに加速させる方向に働き、数値計算によって減衰しようとする渦を維持・増幅させるのです。

3. HLSLによる実装

というわけでCompute Shaderに落とし込んでいきましょう。
処理は2パスで構成するのが効率的です。

Vorticity Kernel: 速度場から渦度を計算する  
Force Kernel: 渦度から閉じ込め力を計算し、速度場に加算する

Phase 1: 渦度の計算

まずは現在の速度場 _VelocityTex から渦度を計算します。
中心差分を用いて偏微分を近似します。

// 2次元シミュレーションの場合の例（結果はスカラー）
Texture2D<float2> _VelocityTex;
RWTexture2D<float> _VorticityTex;

[numthreads(8, 8, 1)]
void ComputeVorticity(uint3 id : SV_DispatchThreadID)
{
    uint2 up = id.xy + uint2(0, 1);
    uint2 down = id.xy - uint2(0, 1);
    uint2 right = id.xy + uint2(1, 0);
    uint2 left = id.xy - uint2(1, 0);

    // 速度のサンプリング
    float2 u_up = _VelocityTex[up];
    float2 u_down = _VelocityTex[down];
    float2 u_right = _VelocityTex[right];
    float2 u_left = _VelocityTex[left];

    // rot = dv/dx - du/dy
    // 中心差分: (f(x+1) - f(x-1)) / 2dx
    // 係数 0.5 は省略して後でスケーリングしても良いが、ここでは記述する
    float dv_dx = (u_right.y - u_left.y) * 0.5;
    float du_dy = (u_up.x - u_down.x) * 0.5;

    float vorticity = dv_dx - du_dy;

    _VorticityTex[id.xy] = vorticity;
}

Phase 2: 閉じ込め力の適用

次に、計算した渦度を用いて力を生成し、速度場を更新します。
ここが技術的なハイライトです。
ゼロ除算回避なども考慮します。

Texture2D<float> _VorticityTex;
RWTexture2D<float2> _VelocityTex;
float _Epsilon; // 強度係数 (例: 5.0 ~ 10.0)
float _DeltaTime;

[numthreads(8, 8, 1)]
void ApplyVorticityForce(uint3 id : SV_DispatchThreadID)
{
    // 近傍の渦度絶対値を取得 (|omega|)
    float v_c = abs(_VorticityTex[id.xy]);
    float v_r = abs(_VorticityTex[id.xy + uint2(1, 0)]);
    float v_l = abs(_VorticityTex[id.xy - uint2(1, 0)]);
    float v_u = abs(_VorticityTex[id.xy + uint2(0, 1)]);
    float v_d = abs(_VorticityTex[id.xy - uint2(0, 1)]);

    // 1. 渦度強度の勾配 (Gradient of Magnitude) を計算
    // N = (d|w|/dx, d|w|/dy, 0)
    float2 N;
    N.x = (v_r - v_l) * 0.5;
    N.y = (v_u - v_d) * 0.5;

    // 2. 正規化 (Normalization)
    // 極小値でのゼロ除算を防ぐための小さなイプシロンを加算
    float magSqr = dot(N, N);
    N = (magSqr > 0.00001) ? N * rsqrt(magSqr) : float2(0, 0);

    // 3. 外積による力の方向決定 (Cross Product)
    // 2Dの場合、N = (Nx, Ny, 0), omega = (0, 0, w) とみなす
    // Force = N x omega = (Ny * w, -Nx * w, 0)
    float vorticity = _VorticityTex[id.xy];
    float2 force = float2(N.y * vorticity, -N.x * vorticity);

    // 4. 速度場への加算
    // F = epsilon * Force
    // u_new = u + F * dt
    float2 currentVelocity = _VelocityTex[id.xy];
    _VelocityTex[id.xy] = currentVelocity + force * (_Epsilon * _DeltaTime);
}

4. 実装上の注意点とチューニング

係数 $\epsilon$ の魔力

変数 _Epsilon は魔法の数字です。

小さすぎる場合: 効果が薄く、すぐに拡散してしまいます。  
大きすぎる場合: エネルギーが過剰に注入され、シミュレーションが発散（爆発）します。

通常、グリッド解像度やタイムステップに依存しますが、0.1 から 10.0 程度の範囲で、見た目を確認しながら調整する必要があります。

境界条件

壁際などの境界付近では、勾配計算（微分）が不安定になりがちです。
境界条件を適切に設定し、壁の外側の値を0にする、あるいは壁の法線方向の微分を0にするなどの処理を入れておかないと、壁際から謎の無限エネルギーが発生することがあります。

実装結果

渦度閉じ込めは、物理的に「正しい」挙動ではありません。
エネルギー保存則を無視して、人工的に回転力を加えているからです。
しかし、グラフィックスにおいては「見た目のリッチさ」こそが正義です。


▲ 適用前（拡散して消える）	▲ 適用後（微細な乱流が持続する）

これだけで、単調だった流体が、生き物のようにうねるダイナミックな流体へと進化します。

まとめ

というわけで色々実装した結果です。

今回はシミュレーション部分を実装していきました。
ナビエストークスの実装は比較的によく見るのですが、マコーマック法や渦度閉じ込めになるとそんなに記事を見かけないので、今回はその辺を多めに解説してみました。

次回は燃焼シミュレーションに関しての記事になる予定です。

SPARK CREATIVE Tech Blog

https://spark-group.jp/