🦔

Stable Fluids

2024/10/30に公開

物理シミュレーション

流体

tech

ステップ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方向の逆補間比率

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

各要素の意味
- 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
  - 補間比率
補間の計算手順
ここでは、まず 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);  // 境界条件を設定