SPH法による流体シミュレーションについて

前置き

こんにちは。

今日は、マーチンキューブについて記事を書きたかったんですが、そもそもマーチンキューブって何の役に立つのかを知らないと興味を持てないと思うので、まず先にマーチンキューブがよく活躍する流体シミュレーションについて紹介したいと思います。

流体シミュレーションとは何かというと、水とかの動きをシミュレーションする事です。ゲームに流体シミュレーションを活用すると、リアルな水の表現が可能になるかもしれません。

流体シミュレーションには色々な手法がありますが、私が知ってるのはSPH法です。なのでこの記事ではSPH法について紹介します。

SPHとはSmoothed Particle Hydrodynamicsの略です。直訳すると「平滑化された粒子の流体力学」です。名前の通り、SPH法では水を沢山の粒子（粒々）の集まりとして表現します。というのも、実際に水は水分子の粒々の集まりだからです。それぞれの粒々の動きを一つづつ物理シミュレーションする事で全体として水みたいな動きになります。

「沢山の粒々を物理シミュレーションするだけなら、SPHとか使わなくてもUnityで大量のスフィアコライダ突っ込めばいいんじゃね？」

はい、今読者の心の声が聞こえてきましたね。
たしかにぶっちゃけUnityで大量のスフィアコライダをぶち込むだけでも割とそれっぽくシミュレートできるかもしれません。

ただし、Unityの物理演算では水の粘性とかが考慮されないので、よりリアルな水表現を追い求めるならSPHを使うといいかもしれません。
ゲームに使うだけならそれっぽく見えればそれでいいかもしれませんが、流体力学を研究してる人とかは、当然一番リアルに水を再現できる手法を使いたいでしょうね。

SPHを使うとこんな水表現ができるらしいです↓

ちなみになんで私がSPHについて知ってるかと言うと、私は学生の頃の卒論でSPHとマーチンキューブで水を描画する研究やってました。十数年前の話ですが。

SPH法

というわけでSPHをやっていきましょう。

沢山の粒子で水を表現するという事なので、早速つぶつぶをばら撒いてみました。

さて、それぞれの粒子をどうやって動かせばいいんですか？
そのためには、流体を表現するナビエストークス方程式というのを使います。
SPH向けのナビエストークス方程式はこれです↓

「うわ…わけわかんねえ数式出てきた…やばい…こわい…ブラウザ閉じよ…」
まあ大丈夫です。冷静に眺めるとこの方程式の言わんとする事はこういう事です。

つまり、”ある粒子の加速度は、その粒子の位置の圧力項と粘性項と外力（重力とか）を足し合わせたものですよ”と言っているだけにすぎない。

簡単でしょ？

「でもなんか三角形のわけわからん記号∇とか出てるやん、こんなん学校で習ってねえぞ」
∇ってのはナブラって言って、なんかベクトルとかスカラー場に付けると勾配（傾き）を意味する？らしいです。偏微分とかして求めるらしいですが私もよく知りません。その内勉強するかもしれません。まあ今はぶっちゃけどうでもよくないっすか。どうせ実装する時はコードにナブラとか出てこないんで。

離散化

私もナブラとか言われてもよく分からんのと同じように、パソコンだってナブラとか言われてもよく分かりません。
パソコン「数式とかじゃなくてもっとソースコードで説明してくれないと分かんない」
パソコンが喋った！

パソコンは数式とか分からないので、パソコンでもわかるように近似的に計算で求められる形に離散化してあげる必要があります。

ところで、一つの水粒子の場所の密度とか圧力って、そこの周りにどんだけ他の粒子が集まってるかで決まりそうですよね。
つまり、近くの粒子同士が影響し合って、かつ粒子同士の距離が近いほど大きな影響を与えることになります。
その事を数式として表現するとこうなります。

Ｗなんちゃらは”重み関数”と言って、距離が離れるほど影響が小さくなる感じの関数が入ります。距離がhを超えると影響ゼロになるようになってます。

とにかくこれを使って密度やら圧力やらを計算できます。

密度を求める

密度を求めるには、密度ρを上の(2)の方程式のφのとこに突っ込めばＯＫです。

「Ｗpoly6って何やねん！？」
これはね、SPHの研究してる人が考えた、「これ使うと密度がそれっぽく求められるよ」っていう重み関数です。具体的には、

こんな感じです。これは３次元の時に上手く動くヤツで、2Dでシミュレートする時はまた別の重み関数を使います。
otherwiseってなんやねん？と言うと、つまり粒子間の距離がhより小さい場合だけ計算すればよくて、hより大きい時は０、つまり無視してくださいという意味です。

数式だとサッパリ意味不明ですが、実際にC#の疑似コードで書くと大して難しくないです。

//粒子
class Particle
{
    public Vector3 position; //座標
    public Vector3 velocity; //速度
    public Vectore force;    //力
    public float   density;  //密度
    public float   pressure; //圧力
}

float h = 0.012; //影響半径
float particleMass = 0.0002f; //粒子の質量
float densityCoef = particleMass * 315f / (64f * Mathf.PI * Mathf.Pow(h,9)); //密度計算で使うヤツ

//粒子の密度計算
void CalcDensity(Particle[] particles)
{
    float h2 = h*h; //事前にhの二乗を計算しておく
    for(int i = 0; i< particles.length; i++) //一つづつ粒子の密度を計算
    {
        var nowParticle = particles[i]; //今回計算する粒子
        float sum = 0; //足し合わせる変数
        for(int j =0; j< particles.length; j++) //他の粒子全てについて
        {
            if(i==j){continue;} //自分自身だったらスキップ
            var nearParticle = particles[j];
            
            var diff = nearParticle.position - nowParticle.position; //粒子距離
            float r2 = Vector3.Dot(diff,diff); //粒子距離の２乗

            //粒子距離がhより小さい場合だけ計算する
            if( r2 < h2 )
            {
                float c = h2 - r2;
                sum += Mathf.Pow(c,3); //(h2-r2)の３乗
            }
        }

        nowParticle.density = sum * densityCoef; //密度が求まった
    }
}

//粒子

class Particle

{

public Vector3 position; //座標

public Vector3 velocity; //速度

public Vectore force; //力

public float density; //密度

public float pressure; //圧力

}

float h = 0.012; //影響半径

float particleMass = 0.0002f; //粒子の質量

float densityCoef = particleMass * 315f / (64f * Mathf.PI * Mathf.Pow(h,9)); //密度計算で使うヤツ

//粒子の密度計算

void CalcDensity(Particle[] particles)

{

float h2 = h*h; //事前にhの二乗を計算しておく

for(int i = 0; i< particles.length; i++) //一つづつ粒子の密度を計算

{

var nowParticle = particles[i]; //今回計算する粒子

float sum = 0; //足し合わせる変数

for(int j =0; j< particles.length; j++) //他の粒子全てについて

{

if(i==j){continue;} //自分自身だったらスキップ

var nearParticle = particles[j];

var diff = nearParticle.position - nowParticle.position; //粒子距離

float r2 = Vector3.Dot(diff,diff); //粒子距離の２乗

//粒子距離がhより小さい場合だけ計算する

if( r2 < h2 )

{

float c = h2 - r2;

sum += Mathf.Pow(c,3); //(h2-r2)の３乗

}

nowParticle.density = sum * densityCoef; //密度が求まった

}

簡単でしょ？
このコードは適当に書きましたが、まあこんな感じで密度が求まります。

圧力

圧力を求めるには、状態方程式というのを使います。

kは圧力を調整するための単なるパラメータです。ρ0は基本密度です。密度が基本密度に戻ろうとする圧力が流体にかかります。

C#の疑似コードで書くとこうです

float pressureStiffness = 200f; //圧力係数
float restDensity = 1000f; //静止密度

//粒子の圧力計算
void CalcPressure( Particle[] particles)
{
    for(int i = 0; i< particles.length; i++) //一つづつ粒子の圧力を計算
    {
        particles[i].pressure = pressureStiffness * ( particles[i].density - restDensity );
    }
}

float pressureStiffness = 200f; //圧力係数

float restDensity = 1000f; //静止密度

//粒子の圧力計算

void CalcPressure( Particle[] particles)

{

for(int i = 0; i< particles.length; i++) //一つづつ粒子の圧力を計算

{

particles[i].pressure = pressureStiffness * ( particles[i].density - restDensity );

}

圧力項

圧力が求まったので、圧力項も密度と同じノリで計算できます。
(1)のナビエストークス方程式の圧力項に(2)を適用するとこうなります。

当然こうでしょ？と思うんですが、実はこれだとダメだそうです。
なんか、これだと”作用反作用の法則”を満たさないから、”対称化”する必要があって、

こっちの式を使えだそうです。私も何のこっちゃ分からなくなってきたのですが、しょうがないので言われるがままに従います。

∇Wspikyとかいう重み関数（の勾配）は、

これを使います。今までのノリと同じように計算できると思うので、もうC#コードは載せません。（面倒くさくなってきた）
ノリさえ分かってれば、他の人が書いたSPHのコードとか見てもらえば大体意味が分かると思います。

粘性項

粘性項も圧力項と同じ感じで、これも対称化が必要だそうで、対称化した式はこうなります。

∇²Wviscとかいう重み関数（のラプラシアン）は、

こうです。

衝突判定

ここまででいい感じに圧力項とかが求まってきましたが、そういえば、水槽みたいな壁や床を用意してあげないと、このままでは粒子が重力で真っ逆さまに落ちていきますよね。

壁や床との衝突判定には、ペナルティ法という手法が利用できます。
ペナルティ法では壁を突き抜けそうな粒子に跳ね返す力を与えて壁を突き抜けないようにします。

数式だとこんな感じです。
例えば粒子のx,y,z座標がそれぞれ-1～1の範囲に収まるようにするには、左右、上下、奥手前の６回ペナルティ法の計算を行う事になるでしょう。

粒子の座標が求まった！

というわけで、やっと粒子の加速度が計算できます。
ある粒子の加速度aは、

a = 圧力項＋粘性項＋衝突判定の力＋重力

となります。
加速度さえあれば、速度も求められますね。
速度vは、

v = v + deltaTime * a;

速度があれば、座標もわかります。
座標posは、

pos = pos + deltaTime * v;

簡単ですね。
こんな風に普通に加速度から速度、座標を求める方法は、”前進オイラー法”という手法になります。
これは実は精度が低い方法です。普通のゲームの感じで1秒に60回の実行では発散（ぐっちゃぐちゃになる）してしまうかもしれません。もっとタイムステップを小さく刻む必要があるでしょう。

タイムステップ小さくするのも面倒だな～という場合は、前進オイラー法より安定度の高い手法を使う手があります。
例えばリープ・フロッグ法だとこうなります。

pos = pos + v * deltaTime + 0.5 * a * deltaTime * deltaTime;

v = v + 0.5 * ( prev_a + a ) * deltaTime; //prev_aは前回のフレームのaという意味