図5.1: 左:オイラー的、右:ラグランジュ的
前章では、格子法による流体シミュレーションの作成方法について解説しました。本章では、もう一つの流体のシミュレーション方法である粒子法、特にSPH法を用いて流体の動きを表現していきます。多少噛み砕いて説明を行っているので、不十分な表現などありますがご了承ください。
流体の動きの観測方法として、オイラー的視点とラグランジュ的視点というものが存在します。オイラー的視点とは、流体に等間隔で観測点を固定し、その観測点での流体の動きを解析するものです。一方、ラグランジュ的視点とは、流体の流れに沿って動く観測点を浮かべ、その観測点での流体の動きを観測するものとなります(図5.1参照)。基本的に、オイラー的視点を用いた流体シミュレーション手法のことを格子法、ラグランジュ的視点を用いた流体シミュレーション手法のことを粒子法と呼びます。
図5.1: 左:オイラー的、右:ラグランジュ的
オイラー的視点とラグランジュ的視点では、微分の演算の仕方が異なります。はじめに、オイラー的視点で表された物理量*1を以下に示してみます。
[*1] 物理量とは、観測できる速度や質量などのことを指します。 端的には「単位が有るもの」と捉えて良いでしょう。
\phi = \phi (\overrightarrow{x}, t)
これは、時刻tで位置\overrightarrow{x}にある物理量\phiという意味になります。この物理量の時間微分は、
\frac{\partial \phi}{\partial t}
と表せます。もちろんこれは、物理量の位置が\overrightarrow{x}で固定されていますので、オイラー的視点での微分になります。
[*2] 流れに沿った観測点の移動のことを、移流と呼びます。
一方、ラグランジュ的視点では、観測点を流れに沿って移動*2させますので、観測点自体も時間の関数となっています。そのため、初期状態で\overrightarrow{x}_0にあった観測点は、時刻tで
\overrightarrow{x}(\overrightarrow{x}_0, t)
に存在します。 よって物理量の表記も
\phi = \phi (\overrightarrow{x}(\overrightarrow{x}_0, t), t)
となります。微分の定義に従って、現在の物理量と\Delta t秒後の物理量の変化量を見てみると
\displaystyle \lim_{\Delta t \to 0} \frac{\phi(\overrightarrow{x}(\overrightarrow{x}_0, t + \Delta t), t + \Delta t) - \phi(\overrightarrow{x}(\overrightarrow{x}_0, t), t)}{\Delta t}
= \sum_i \frac{\partial \phi}{\partial x_i} \frac{\partial x_i}{\partial t} + \frac{\partial \phi}{\partial t}
= \left( \left( \begin{matrix}u_1\\u_2\\u_3\end{matrix} \right)
\cdot
\left( \begin{matrix} \frac{\partial}{\partial x_1}\\\frac{\partial}{\partial x_2}\\\frac{\partial}{\partial x_3} \end{matrix} \right)
+ \frac{\partial}{\partial t}
\right) \phi\\
= (\frac{\partial}{\partial t} + \overrightarrow{u} \cdot {grad}) \phi
となります。これが、観測点の移動を考慮した物理量の時間微分となります。しかしながら、この表記を用いていては式が複雑になりますので、
\dfrac{D}{Dt} := \frac{\partial}{\partial t} + \overrightarrow{u} \cdot {grad}
という演算子を導入することで、短く表すことができます。これら、観測点の移動を考慮した一連の操作を、ラグランジュ微分と呼びます。一見複雑そうに見えますが、観測点が移動する粒子法では、ラグランジュ的視点で式を表した方が都合が良くなります。
流体は、流体の速度が音速よりも十分に小さい場合、体積の変化が起きないとみなすことができます。これは流体の非圧縮条件と呼ばれ、以下の数式で表されます。
\nabla \cdot \overrightarrow{u} = 0
これは、流体内で湧き出しや消失がないことを示しています。この式の導出には少し複雑な積分が入りますので、説明は割愛*3します。「流体は圧縮しない！」程度に捉えておいてください。
[*3] "Fluid Simulation for Computer Graphics - Robert Bridson" で詳しく解説されています。
粒子法では、流体を小さな粒子によって分割し、ラグランジュ的視点で流体の動きを観測します。この粒子が、前節の観測点にあたります。 一口に「粒子法」といっても、現在では多くの手法が提案されており、有名なものとして
などがあります。
はじめに、粒子法におけるナビエ・ストークス方程式(以下NS方程式)は、以下のように記述されます。
\dfrac{D \overrightarrow{u}}{Dt} = -\dfrac{1}{\rho}\nabla p + \nu \nabla \cdot \nabla \overrightarrow{u} + \overrightarrow{g}
\label{eq:navier}
前章の格子法で出てきたNS方程式とは少し形が異なりますね。移流項がまるまる抜けてしまっていますが、先程のオイラー微分とラグランジュ微分の関係を見てみると、うまくこの形に変形できることがわかります。粒子法では観測点を流れに沿って移動させますから、NS方程式計算時に移流項を考慮する必要がありません。移流の計算はNS方程式で算出した加速度をもとに粒子位置を直接更新することで済ませる事ができます。
現実の流体は分子の集まりですので、ある種のパーティクルシステムであると言うことができます。しかし、コンピュータで実際の分子の数の計算を行うのは不可能ですので、計算可能な大きさに調節してあげる必要があります。図5.2に示されているそれぞれの粒(*4)は、計算可能な大きさで分割した流体の一部分を表していています。これらの粒は、それぞれ質量m、位置ベクトル\overrightarrow{x}、速度ベクトル\overrightarrow{u}、体積Vを持つと考えることができます。
図5.2: 流体のパーティクル近似
これらそれぞれの粒について、外から受けた力\overrightarrow{f}を計算し、運動方程式m \overrightarrow{a} = \overrightarrow{f}を解くことで加速度が算出され、次のタイムステップでどのように移動するかを決めることができます。
[*4] 英語では'Blob'と呼ばれます
前述の通り、それぞれの粒子は周りから何らかの力を受けて動きますが、その「力」とは一体何でしょうか。簡単な例として、重力m \overrightarrow{g}があげられますが、それ以外に周りの粒子からも何らかの力を受けるはずです。これらの力について、以下に解説します。
流体粒子にかかる力の1つ目は、圧力です。流体は必ず圧力の高い方から低い方に向かって流れます。もし圧力がどの方向からも同じだけかかっていたとすると、力は打ち消されて動きが止まってしまいますから、圧力のバランスが不均一である場合を考えます。前章で述べられたように、圧力のスカラー場の勾配を取ることで、自分の粒子位置から見て最も圧力上昇率の高い方向を算出することができます。粒子が力を受ける方向は、圧力の高い方から低い方ですので、マイナスを取って-\nabla pとなります。また、粒子は体積を持っていますから、粒子にかかる圧力は、-\nabla pに粒子の体積をかけて算出します*5。最終的に、- V \nabla pという結果が導出されます。
[*5] 流体の非圧縮条件により、単に体積をかけるだけで粒子にかかる圧力の積分を表すことができます。
流体粒子にかかる力の２つ目は、粘性力です。粘性(ねばりけ)のある流体とは、はちみつや溶かしたチョコレートなどに代表される、変形しづらい流体のことを指します。粘性があるという言葉を粒子法の表現に当てはめてみると、粒子の速度は、周りの粒子速度の平均をとりやすいということになります。前章で述べられた通り、周囲の平均をとるという演算は、ラプラシアンを用いて行うことができます。
粘性の度合いを動粘性係数\muを用いて表すと、\mu \nabla \cdot \nabla \overrightarrow{u}と表す事ができます。
これらの力を運動方程式m \overrightarrow{a} = \overrightarrow{f}に当てはめて整理すると、
m \dfrac{D\overrightarrow{u}}{Dt} = - V \nabla p + V \mu \nabla \cdot \nabla \overrightarrow{u} + m\overrightarrow{g}
ここで、mは\rho Vであることから、変形して(Vが打ち消されます)
\rho \dfrac{D\overrightarrow{u}}{Dt} = - \nabla p + \mu \nabla \cdot \nabla \overrightarrow{u} + \rho \overrightarrow{g}
両辺\rhoで割り、
\dfrac{D\overrightarrow{u}}{Dt} = - \dfrac{1}{\rho}\nabla p + \dfrac{\mu}{\rho} \nabla \cdot \nabla \overrightarrow{u} + \overrightarrow{g}
最後に、粘性項の係数\dfrac{\mu}{\rho}に\nuを導入して、
\dfrac{D\overrightarrow{u}}{Dt} = - \dfrac{1}{\rho}\nabla p + \nu \nabla \cdot \nabla \overrightarrow{u} + \overrightarrow{g}
となり、はじめに挙げたNS方程式を導出することができました。
粒子法では、粒子自体が流体の観測点を表現しているので、移流項の計算は単に粒子位置を移動させるだけで完了します。実際の時間微分の計算では、無限に小さい時間を用いますが、コンピューターでの計算では無限を表現できないため、十分小さい時間\Delta tを用いて微分を表現します。これを差分と言い、\Delta tを小さくすればするほど、正確な計算を行うことができます。
加速度について、差分の表現を導入すると、
\overrightarrow{a} = \dfrac{D\overrightarrow{u}}{Dt} \equiv \frac{\Delta \overrightarrow{u}}{\Delta t}
となります。よって速度の増分\Delta \overrightarrow{u}は、
\Delta \overrightarrow{u} = \Delta t \overrightarrow{a}
となり、また、位置の増分についても同様に、
\overrightarrow{u} = \frac{\partial \overrightarrow{x}}{\partial t} \equiv \frac{\Delta \overrightarrow{x}}{\Delta t}
より、
\Delta \overrightarrow{x} = \Delta t \overrightarrow{u}
となります。
この結果を利用することで、次のフレームでの速度ベクトルと位置ベクトルを算出できます。現在のフレームでの粒子速度が\overrightarrow{u}_nであるとすると、次のフレームでの粒子速度は\overrightarrow{u}_{n+1}で、
\overrightarrow{u}_{n+1} = \overrightarrow{u}_n + \Delta \overrightarrow{u} = \overrightarrow{u}_n + \Delta t \overrightarrow{a}
と表せます。
現在のフレームでの粒子位置が\overrightarrow{x}_nであるとすると、次のフレームでの粒子位置は\overrightarrow{x}_{n+1}で、
\overrightarrow{x}_{n+1} = \overrightarrow{x}_n + \Delta \overrightarrow{x} = \overrightarrow{x}_n + \Delta t \overrightarrow{u}
と表せます。
この手法は、前進オイラー法と呼ばれます。これを毎フレーム繰り返すことで、各時刻での粒子の移動を表現することができます。
前節では、粒子法におけるNS方程式の導出方法について解説しました。もちろん、これらの微分方程式をコンピュータでそのまま解くことはできませんので、何らかの近似をしてあげる必要があります。その手法として、CG分野でよく用いられるSPH法について解説します。
SPH法は、本来宇宙物理学における天体同士の衝突シミュレーションに用いられていた手法ですが、1996年にDesbrunら*6によってCGにおける流体シミュレーションにも応用されました。また、並列化も容易で、現在のGPUでは大量の粒子の計算をリアルタイムに行うことが可能です。コンピュータシミュレーションでは、連続的な物理量を離散化して計算を行う必要がありますが、この離散化を、重み関数と呼ばれる関数を用いて行う手法をSPH法と呼びます。
[*6] Desbrun and Cani, Smoothed Particles: A new paradigm for animating highly deformable bodies, Eurographics Workshop on Computer Animation and Simulation (EGCAS), 1996.
SPH法では、粒子一つ一つが影響範囲を持っていて、他の粒子と距離が近いほどその粒子の影響が大きく受けるという動作をします。この影響範囲を図示すると図5.3のようになります。
図5.3: 2次元の重み関数
この関数を重み関数*7と呼びます。
[*7] 通常この関数はカーネル関数とも呼ばれますが、ComputeShaderにおけるカーネル関数と区別するためこの呼び方にしています。
SPH法における物理量を\phiとすると、重み関数を用いて以下のように離散化されます。
\phi(\overrightarrow{x}) = \sum_{j \in N}m_j\frac{\phi_j}{\rho_j}W(\overrightarrow{x_j} - \overrightarrow{x}, h)
N, m, \rho, hはそれぞれ、近傍粒子の集合、粒子の質量、粒子の密度、重み関数の影響半径です。また、関数Wが先ほど述べた重み関数になります。
さらに、この物理量には、勾配とラプラシアンなどの偏微分演算が適用でき、勾配は、
\nabla \phi(\overrightarrow{x}) = \sum_{j \in N}m_j\frac{\phi_j}{\rho_j} \nabla W(\overrightarrow{x_j} - \overrightarrow{x}, h)
ラプラシアンは、
\nabla^2 \phi(\overrightarrow{x}) = \sum_{j \in N}m_j\frac{\phi_j}{\rho_j} \nabla^2 W(\overrightarrow{x_j} - \overrightarrow{x}, h)
と表せます。式からわかるように、物理量の勾配及びラプラシアンは、重み関数に対してのみ適用されるイメージになります。重み関数Wは、求めたい物理量によって異なるものを使用しまが、この理由の説明については割愛*8します。
[*8] "CGのための物理シミュレーションの基礎 - 藤澤誠" で詳しく解説されています。
流体の粒子の密度は、先ほどの重み関数で離散化した物理量の式を利用して、
\rho(\overrightarrow{x}) = \sum_{j \in N}m_jW_{poly6}(\overrightarrow{x_j} - \overrightarrow{x}, h)
と与えられます。ここで、利用する重み関数Wは、以下で与えられます。
図5.4: Poly6重み関数
粘性項を離散化も密度の場合と同様重み関数を利用して、
f_{i}^{visc} = \mu\nabla^2\overrightarrow{u}_i = \mu \sum_{j \in N}m_j\frac{\overrightarrow{u}_j - \overrightarrow{u}_i}{\rho_j} \nabla^2 W_{visc}(\overrightarrow{x_j} - \overrightarrow{x}, h)
と表されます。ここで、重み関数のラプラシアン\nabla^2 W_{visc}は、以下で与えられます。
図5.5: Viscosity重み関数のラプラシアン
同様に、圧力項を離散化していきます。
f_{i}^{press} = - \frac{1}{\rho_i} \nabla p_i = - \frac{1}{\rho_i} \sum_{j \in N}m_j\frac{p_j - p_i}{2\rho_j} \nabla W_{spiky}(\overrightarrow{x_j} - \overrightarrow{x}, h)
ここで、重み関数の勾配W_{spiky}は以下で与えられます。
図5.6: Spiky重み関数の勾配
この時、粒子の圧力は事前に、Tait方程式と呼ばれる、
p = B\left\{\left(\frac{\rho}{\rho_0}\right)^\gamma - 1\right\}
で算出されています。 ここで、Bは気体定数です。非圧縮性を保証するためには、本来ポアソン方程式を解かなければならないのですが、リアルタイム計算には向きません。その代わりSPH法*9では、近似的に非圧縮性を確保する点で格子法よりも圧力項の計算が苦手であるといわれます。
[*9] Tait方程式を用いた圧力計算を行うSPH法を、特別にWCSPH法と呼びます。
サンプルはこちらのリポジトリ(https://github.com/IndieVisualLab/UnityGraphicsProgramming)のAssets/SPHFluid以下に掲載しています。今回の実装では、極力シンプルにSPHの手法を解説するために高速化や数値安定性は考慮していませんのでご了承ください。
シミュレーションに使用する諸々のパラメータの説明については、コード内コメントに記載しています。
リスト5.1: シミュレーションに使用するパラメータ(FluidBase.cs)
1: NumParticleEnum particleNum = NumParticleEnum.NUM_8K; // 粒子数 2: float smoothlen = 0.012f; // 粒子半径 3: float pressureStiffness = 200.0f; // 圧力項係数 4: float restDensity = 1000.0f; // 静止密度 5: float particleMass = 0.0002f; // 粒子質量 6: float viscosity = 0.1f; // 粘性係数 7: float maxAllowableTimestep = 0.005f; // 時間刻み幅 8: float wallStiffness = 3000.0f; // ペナルティ法の壁の力 9: int iterations = 4; // イテレーション回数 10: Vector2 gravity = new Vector2(0.0f, -0.5f); // 重力 11: Vector2 range = new Vector2(1, 1); // シミュレーション空間 12: bool simulate = true; // 実行 or 一時停止 13: 14: int numParticles; // パーティクルの個数 15: float timeStep; // 時間刻み幅 16: float densityCoef; // Poly6カーネルの密度係数 17: float gradPressureCoef; // Spikyカーネルの圧力係数 18: float lapViscosityCoef; // Laplacianカーネルの粘性係数
今回のデモシーンでは、コードに記載されているパラメータの初期化値とは異なる値をインスペクタで設定していますので注意してください。
重み関数の係数はシミュレーション中で変化しないため、初期化時にCPU側で計算しておきます。(ただし、実行途中でパラメータを編集する可能性も踏まえてUpdate関数内で更新しています)
今回、粒子ごとの質量はすべて一定にしているので、物理量の式内にある質量mはシグマの外に出て以下になります。
\phi(\overrightarrow{x}) = m \sum_{j \in N}\frac{\phi_j}{\rho_j}W(\overrightarrow{x_j} - \overrightarrow{x}, h)
そのため、係数計算の中に質量を含めてしまうことができます。
重み関数の種類で係数も変化してきますから、それぞれに関して係数を計算します。
リスト5.2: 重み関数の係数の事前計算(FluidBase.cs)
1: densityCoef = particleMass * 4f / (Mathf.PI * Mathf.Pow(smoothlen, 8)); 2: gradPressureCoef 3: = particleMass * -30.0f / (Mathf.PI * Mathf.Pow(smoothlen, 5)); 4: lapViscosityCoef 5: = particleMass * 20f / (3 * Mathf.PI * Mathf.Pow(smoothlen, 5));
最終的に、これらのCPU側で計算した係数(及び各種パラメータ)をGPU側の定数バッファに格納します。
リスト5.3: ComputeShaderの定数バッファに値を転送する(FluidBase.cs)
1: fluidCS.SetInt("_NumParticles", numParticles);
2: fluidCS.SetFloat("_TimeStep", timeStep);
3: fluidCS.SetFloat("_Smoothlen", smoothlen);
4: fluidCS.SetFloat("_PressureStiffness", pressureStiffness);
5: fluidCS.SetFloat("_RestDensity", restDensity);
6: fluidCS.SetFloat("_Viscosity", viscosity);
7: fluidCS.SetFloat("_DensityCoef", densityCoef);
8: fluidCS.SetFloat("_GradPressureCoef", gradPressureCoef);
9: fluidCS.SetFloat("_LapViscosityCoef", lapViscosityCoef);
10: fluidCS.SetFloat("_WallStiffness", wallStiffness);
11: fluidCS.SetVector("_Range", range);
12: fluidCS.SetVector("_Gravity", gravity);
リスト5.4: ComputeShaderの定数バッファ(SPH2D.compute)
1: int _NumParticles; // 粒子数 2: float _TimeStep; // 時間刻み幅(dt) 3: float _Smoothlen; // 粒子半径 4: float _PressureStiffness; // Beckerの係数 5: float _RestDensity; // 静止密度 6: float _DensityCoef; // 密度算出時の係数 7: float _GradPressureCoef; // 圧力算出時の係数 8: float _LapViscosityCoef; // 粘性算出時の係数 9: float _WallStiffness; // ペナルティ法の押し返す力 10: float _Viscosity; // 粘性係数 11: float2 _Gravity; // 重力 12: float2 _Range; // シミュレーション空間 13: 14: float3 _MousePos; // マウス位置 15: float _MouseRadius; // マウスインタラクションの半径 16: bool _MouseDown; // マウスが押されているか
リスト5.5: 密度の計算を行うカーネル関数(SPH2D.compute)
1: [numthreads(THREAD_SIZE_X, 1, 1)]
2: void DensityCS(uint3 DTid : SV_DispatchThreadID) {
3: uint P_ID = DTid.x; // 現在処理しているパーティクルID
4:
5: float h_sq = _Smoothlen * _Smoothlen;
6: float2 P_position = _ParticlesBufferRead[P_ID].position;
7:
8: // 近傍探索(O(n^2))
9: float density = 0;
10: for (uint N_ID = 0; N_ID < _NumParticles; N_ID++) {
11: if (N_ID == P_ID) continue; // 自身の参照回避
12:
13: float2 N_position = _ParticlesBufferRead[N_ID].position;
14:
15: float2 diff = N_position - P_position; // 粒子距離
16: float r_sq = dot(diff, diff); // 粒子距離の2乗
17:
18: // 半径内に収まっていない粒子は除外
19: if (r_sq < h_sq) {
20: // 計算には2乗しか含まれないのでルートをとる必要なし
21: density += CalculateDensity(r_sq);
22: }
23: }
24:
25: // 密度バッファを更新
26: _ParticlesDensityBufferWrite[P_ID].density = density;
27: }
本来であれば粒子を全数調査せず、適切な近傍探索アルゴリズムを用いて近傍粒子を探す必要がありますが、今回の実装では簡単のために全数調査を行っています(10行目のforループ)。また、自分と相手粒子との距離計算を行うため、11行目で自身の粒子同士で計算を行うのを回避しています。
重み関数の有効半径hによる場合分けは19行目のif文で実現します。密度の足し合わせ(シグマの計算)は、9行目で0で初期化しておいた変数に対してシグマ内部の計算結果を加算していくことで実現します。ここで、もう一度密度の計算式を示します。
\rho(\overrightarrow{x}) = \sum_{j \in N}m_jW_{poly6}(\overrightarrow{x_j} - \overrightarrow{x}, h)
密度の計算は上式のとおり、Poly6重み関数を用います。 Poly6重み関数はリスト5.6で計算します。
リスト5.6: 密度の計算(SPH2D.compute)
1: inline float CalculateDensity(float r_sq) {
2: const float h_sq = _Smoothlen * _Smoothlen;
3: return _DensityCoef * (h_sq - r_sq) * (h_sq - r_sq) * (h_sq - r_sq);
4: }
最終的にリスト5.5の25行目で書き込み用バッファに書き込みます。
リスト5.7: 粒子毎の圧力を計算する重み関数(SPH2D.compute)
1: [numthreads(THREAD_SIZE_X, 1, 1)]
2: void PressureCS(uint3 DTid : SV_DispatchThreadID) {
3: uint P_ID = DTid.x; // 現在処理しているパーティクルID
4:
5: float P_density = _ParticlesDensityBufferRead[P_ID].density;
6: float P_pressure = CalculatePressure(P_density);
7:
8: // 圧力バッファを更新
9: _ParticlesPressureBufferWrite[P_ID].pressure = P_pressure;
10: }
圧力項を解く前に、粒子単位の圧力を算出しておき、後の圧力項の計算コストを下げます。先程も述べましたが、圧力の計算では本来、以下の式のようなポアソン方程式と呼ばれる方程式を解く必要があります。
\nabla^2 p = \rho \frac{\nabla \overrightarrow{u}}{\Delta t}
しかし、コンピュータで正確にポアソン方程式を解く操作は非常に計算コストが高いため、以下のTait方程式を用いて近似的に求めます。
p = B\left\{\left(\frac{\rho}{\rho_0}\right)^\gamma - 1\right\}
リスト5.8: Tait方程式の実装(SPH2D.compute)
1: inline float CalculatePressure(float density) {
2: return _PressureStiffness * max(pow(density / _RestDensity, 7) - 1, 0);
3: }
リスト5.9: 圧力項・粘性項を計算するカーネル関数(SPH2D.compute)
1: [numthreads(THREAD_SIZE_X, 1, 1)]
2: void ForceCS(uint3 DTid : SV_DispatchThreadID) {
3: uint P_ID = DTid.x; // 現在処理しているパーティクルID
4:
5: float2 P_position = _ParticlesBufferRead[P_ID].position;
6: float2 P_velocity = _ParticlesBufferRead[P_ID].velocity;
7: float P_density = _ParticlesDensityBufferRead[P_ID].density;
8: float P_pressure = _ParticlesPressureBufferRead[P_ID].pressure;
9:
10: const float h_sq = _Smoothlen * _Smoothlen;
11:
12: // 近傍探索(O(n^2))
13: float2 press = float2(0, 0);
14: float2 visco = float2(0, 0);
15: for (uint N_ID = 0; N_ID < _NumParticles; N_ID++) {
16: if (N_ID == P_ID) continue; // 自身を対象とした場合スキップ
17:
18: float2 N_position = _ParticlesBufferRead[N_ID].position;
19:
20: float2 diff = N_position - P_position;
21: float r_sq = dot(diff, diff);
22:
23: // 半径内に収まっていない粒子は除外
24: if (r_sq < h_sq) {
25: float N_density
26: = _ParticlesDensityBufferRead[N_ID].density;
27: float N_pressure
28: = _ParticlesPressureBufferRead[N_ID].pressure;
29: float2 N_velocity
30: = _ParticlesBufferRead[N_ID].velocity;
31: float r = sqrt(r_sq);
32:
33: // 圧力項
34: press += CalculateGradPressure(...);
35:
36: // 粘性項
37: visco += CalculateLapVelocity(...);
38: }
39: }
40:
41: // 統合
42: float2 force = press + _Viscosity * visco;
43:
44: // 加速度バッファの更新
45: _ParticlesForceBufferWrite[P_ID].acceleration = force / P_density;
46: }
圧力項、粘性項の計算も、密度の計算方法と同様に行います。
初めに、以下の圧力項による力の計算を31行目にて行っています。
f_{i}^{press} = - \frac{1}{\rho_i} \nabla p_i = - \frac{1}{\rho_i} \sum_{j \in N}m_j\frac{p_j - p_i}{2\rho_j} \nabla W_{press}(\overrightarrow{x_j} - \overrightarrow{x}, h)
シグマの中身の計算は以下の関数で行われます。
リスト5.10: 圧力項の要素の計算(SPH2D.compute)
1: inline float2 CalculateGradPressure(...) {
2: const float h = _Smoothlen;
3: float avg_pressure = 0.5f * (N_pressure + P_pressure);
4: return _GradPressureCoef * avg_pressure / N_density
5: * pow(h - r, 2) / r * (diff);
6: }
次に、以下の粘性項による力の計算を34行目で行っています。
f_{i}^{visc} = \mu\nabla^2\overrightarrow{u}_i = \mu \sum_{j \in N}m_j\frac{\overrightarrow{u}_j - \overrightarrow{u}_i}{\rho_j} \nabla^2 W_{visc}(\overrightarrow{x_j} - \overrightarrow{x}, h)
シグマの中身の計算は以下の関数で行われます。
リスト5.11: 粘性項の要素の計算(SPH2D.compute)
1: inline float2 CalculateLapVelocity(...) {
2: const float h = _Smoothlen;
3: float2 vel_diff = (N_velocity - P_velocity);
4: return _LapViscosityCoef / N_density * (h - r) * vel_diff;
5: }
最後に、リスト5.9の39行目にて圧力項と粘性項で算出した力を足し合わせ、最終的な出力としてバッファに書き込んでいます。
リスト5.12: 衝突判定と位置更新を行うカーネル関数(SPH2D.compute)
1: [numthreads(THREAD_SIZE_X, 1, 1)]
2: void IntegrateCS(uint3 DTid : SV_DispatchThreadID) {
3: const unsigned int P_ID = DTid.x; // 現在処理しているパーティクルID
4:
5: // 更新前の位置と速度
6: float2 position = _ParticlesBufferRead[P_ID].position;
7: float2 velocity = _ParticlesBufferRead[P_ID].velocity;
8: float2 acceleration = _ParticlesForceBufferRead[P_ID].acceleration;
9:
10: // マウスインタラクション
11: if (distance(position, _MousePos.xy) < _MouseRadius && _MouseDown) {
12: float2 dir = position - _MousePos.xy;
13: float pushBack = _MouseRadius-length(dir);
14: acceleration += 100 * pushBack * normalize(dir);
15: }
16:
17: // 衝突判定を書くならここ -----
18:
19: // 壁境界(ペナルティ法)
20: float dist = dot(float3(position, 1), float3(1, 0, 0));
21: acceleration += min(dist, 0) * -_WallStiffness * float2(1, 0);
22:
23: dist = dot(float3(position, 1), float3(0, 1, 0));
24: acceleration += min(dist, 0) * -_WallStiffness * float2(0, 1);
25:
26: dist = dot(float3(position, 1), float3(-1, 0, _Range.x));
27: acceleration += min(dist, 0) * -_WallStiffness * float2(-1, 0);
28:
29: dist = dot(float3(position, 1), float3(0, -1, _Range.y));
30: acceleration += min(dist, 0) * -_WallStiffness * float2(0, -1);
31:
32: // 重力の加算
33: acceleration += _Gravity;
34:
35: // 前進オイラー法で次の粒子位置を更新
36: velocity += _TimeStep * acceleration;
37: position += _TimeStep * velocity;
38:
39: // パーティクルのバッファ更新
40: _ParticlesBufferWrite[P_ID].position = position;
41: _ParticlesBufferWrite[P_ID].velocity = velocity;
42: }
壁との衝突判定をペナルティ法を用いて行います(19-30行目)。ペナルティ法とは、境界位置からはみ出した分だけ強い力で押し返すという手法になります。
本来は壁との衝突判定の前に障害物との衝突判定も行うのですが、今回の実装ではマウスとのインタラクションを行うようにしています(213-218行目)。マウスが押されていれば、指定された力でマウス位置から遠ざかるような力を加えています。
33行目にて外力である重力を加算しています。重力の値をゼロにすると無重力状態になり、面白い視覚効果が得られます。また、位置の更新は前述の前進オイラー法で行い(36-37行目)、最終的な結果をバッファに書き込みます。
リスト5.13: シミュレーションの主要関数(FluidBase.cs)
1: private void RunFluidSolver() {
2:
3: int kernelID = -1;
4: int threadGroupsX = numParticles / THREAD_SIZE_X;
5:
6: // Density
7: kernelID = fluidCS.FindKernel("DensityCS");
8: fluidCS.SetBuffer(kernelID, "_ParticlesBufferRead", ...);
9: fluidCS.SetBuffer(kernelID, "_ParticlesDensityBufferWrite", ...);
10: fluidCS.Dispatch(kernelID, threadGroupsX, 1, 1);
11:
12: // Pressure
13: kernelID = fluidCS.FindKernel("PressureCS");
14: fluidCS.SetBuffer(kernelID, "_ParticlesDensityBufferRead", ...);
15: fluidCS.SetBuffer(kernelID, "_ParticlesPressureBufferWrite", ...);
16: fluidCS.Dispatch(kernelID, threadGroupsX, 1, 1);
17:
18: // Force
19: kernelID = fluidCS.FindKernel("ForceCS");
20: fluidCS.SetBuffer(kernelID, "_ParticlesBufferRead", ...);
21: fluidCS.SetBuffer(kernelID, "_ParticlesDensityBufferRead", ...);
22: fluidCS.SetBuffer(kernelID, "_ParticlesPressureBufferRead", ...);
23: fluidCS.SetBuffer(kernelID, "_ParticlesForceBufferWrite", ...);
24: fluidCS.Dispatch(kernelID, threadGroupsX, 1, 1);
25:
26: // Integrate
27: kernelID = fluidCS.FindKernel("IntegrateCS");
28: fluidCS.SetBuffer(kernelID, "_ParticlesBufferRead", ...);
29: fluidCS.SetBuffer(kernelID, "_ParticlesForceBufferRead", ...);
30: fluidCS.SetBuffer(kernelID, "_ParticlesBufferWrite", ...);
31: fluidCS.Dispatch(kernelID, threadGroupsX, 1, 1);
32:
33: SwapComputeBuffer(ref particlesBufferRead, ref particlesBufferWrite);
34: }
これまでに述べたComputeShaderのカーネル関数を、毎フレーム呼び出す部分です。それぞれのカーネル関数に対して適切なComputeBufferを与えてあげます。
ここで、タイムステップ幅\Delta tを小さくすればするほどシミュレーションの誤差が出にくくなることを思い出してみてください。60FPSで実行する場合、\Delta t = 1 / 60となりますが、これでは誤差が大きく出てしまい粒子が爆発してしまいます。さらに、\Delta t = 1 / 60より小さいタイムステップ幅をとると、1フレーム当たりの時間の進み方が実時間より遅くなり、スローモーションになってしまいます。これを回避するには、\Delta t = 1 / (60 \times {iterarion})として、メインルーチンを1フレームにつきiterarion回回します。
リスト5.14: 主要関数のイテレーション(FluidBase.cs)
1: // 計算精度を上げるために時間刻み幅を小さくして複数回イテレーションする
2: for (int i = 0; i<iterations; i++) {
3: RunFluidSolver();
4: }
こうすることで、小さいタイムステップ幅で実時間のシミュレーションを行うことができます。
通常のシングルアクセスのパーティクルシステムとは異なり、粒子同士が相互作用しますから、計算途中に他のデータが書き換わってしまっては困ります。これを回避するために、GPUで計算を行っている際に値を書き換えない読み込み用バッファと書き込み用バッファの2つを用意します。これらのバッファを毎フレーム入れ替えることで、競合なくデータを更新できます。
リスト5.15: バッファを入れ替える関数(FluidBase.cs)
1: void SwapComputeBuffer(ref ComputeBuffer ping, ref ComputeBuffer pong) {
2: ComputeBuffer temp = ping;
3: ping = pong;
4: pong = temp;
5: }
リスト5.16: パーティクルのレンダリング(FluidRenderer.cs)
1: void DrawParticle() {
2:
3: Material m = RenderParticleMat;
4:
5: var inverseViewMatrix = Camera.main.worldToCameraMatrix.inverse;
6:
7: m.SetPass(0);
8: m.SetMatrix("_InverseMatrix", inverseViewMatrix);
9: m.SetColor("_WaterColor", WaterColor);
10: m.SetBuffer("_ParticlesBuffer", solver.ParticlesBufferRead);
11: Graphics.DrawProcedural(MeshTopology.Points, solver.NumParticles);
12: }
10行目にて、流体粒子の位置計算結果を格納したバッファをマテリアルにセットし、シェーダーに転送します。11行目にて、パーティクルの個数分インスタンス描画をするよう命令しています。
リスト5.17: パーティクルのレンダリング(Particle.shader)
1: struct FluidParticle {
2: float2 position;
3: float2 velocity;
4: };
5:
6: StructuredBuffer<FluidParticle> _ParticlesBuffer;
7:
8: // --------------------------------------------------------------------
9: // Vertex Shader
10: // --------------------------------------------------------------------
11: v2g vert(uint id : SV_VertexID) {
12:
13: v2g o = (v2g)0;
14: o.pos = float3(_ParticlesBuffer[id].position.xy, 0);
15: o.color = float4(0, 0.1, 0.1, 1);
16: return o;
17: }
1-6行目にて、流体粒子の情報を受け取るための情報の定義を行います。この時、スクリプトからマテリアルに転送したバッファの構造体と定義を一致させる必要があります。位置データの受け取りは、14行目のようにid : SV_VertexIDでバッファの要素を参照することで行います。
あとは通常のパーティクルシステムと同様、図5.7のようにジオメトリシェーダーで計算結果の位置データを中心としたビルボード*10を作成し、粒子画像をアタッチしてレンダリングします。
図5.7: ビルボードの作成
[*10] 表が常に視点方向を向くPlaneのことを指します。
図5.8: レンダリング結果
動画はこちら(https://youtu.be/KJVu26zeK2w)に掲載しています。
本章では、SPH法を用いた流体シミュレーションの手法を示しました。SPH法を用いることで、流体の動きをパーティクルシステムのように汎用的に扱うことができるようになりました。
先述の通り、流体シミュレーションの手法はSPH法以外にもたくさんの種類があります。本章を通して、他の流体シミュレーション手法に加え、他の物理シミュレーション自体についても興味を持っていただき、表現の幅を広げていただければ幸いです。