aokomoriuta's blog

青子守歌のブログ

MPS法における圧力ポアソン方程式の境界条件

MPS法やISPH法では圧力のポアソン方程式
\displaystyle \nabla^2 p = S(ここでSは生成項)
が出てくるのですが、色々なところの「MPS法実装してみました!」記事を見てみるとこのポアソン方程式の特に境界条件の扱い方が間違っているのが散見される(&よく「なんか上手く動かないんだけどどうしたらいい?」と聞かれる)ので、まとめておきます(別に新規性のあることじゃないし)。

まず、標準MPS法であれ高精度MPS法であれ、どんなラプラシアン・モデルを使っていても、基本的にラプラシアンは、自粒子と相手粒子との物理量の差の重み付け総和
\displaystyle \nabla^2 \phi_i = \sum_{j \ne i} \left( \phi_j - \phi_i \right) W\left( r_{ij} \right)
で表されるはずです(※このWはいわゆるMPS法の重み関数やSPH法の核関数とは限らない)。

で、これがよく見てみると \phi_iに関する線形和になっているので、このポアソン方程式は、 \phi_iを未知変数とする線形方程式A \mathbf{x} = \mathbf{b}として解けて、その時の係数行列Aのi行j列目成分a_{ij}
\displaystyle a_{ij} = \begin{cases}
\displaystyle W\left( r_{ij} \right) & \left(j \ne i \right)\\
\displaystyle -\sum_{j \ne i} W\left( r_{ij} \right) = -\sum_{j \ne i} a_{ij} & \left(j = i \right)\\
\end{cases}
と書けるところまでは、まぁ誰でも分かるでしょう。

というわけでこれを実装すれば話は済むんですが、さてここが問題。粒子の種類は以下のように色々ありますが、このiとjにはどれを含めてどれを含めるべきではないでしょうか?

  • 流体粒子(非圧縮性流れにおけるニュートン流体)
    • 自由表面の流体粒子(流体粒子のうち、自由表面であると判定された粒子)
    • 内部流体粒子(自由表面でない流体粒子)
  • 固定粒子(いわゆる圧力を計算する壁)
    • 自由表面の固定粒子(固定粒子のうち、自由表面であると判定された粒子)
    • 内部固定粒子(自由表面でない固定粒子)
  • ダミー粒子(いわゆる圧力を計算しない壁)
  • 無効粒子(計算に全く関与しない粒子、画面外に出た粒子とか。別名:幽霊粒子)

一番最後の無効粒子を計算に含める人は、まぁいないでしょう。
内部粒子は全て含める必要があるのも当然でしょう。

そうなると残るは自由表面粒子とダミー粒子の扱いです。
これらの扱いを間違えていることで圧力をきちんと計算できず、発散したり圧力の擾乱がひどくなる例があります。
分かりやすくするために鉛直1次元の簡単な例を用意してみました。
この例では、mainシートにある通り、初期粒子間距離0.5[m]で0.001[s]自由落下させた時の圧力を標準MPS法で計算します。

まず一番多い間違いは、ダミー粒子を計算に含めている場合です。
全部含めてしまった結果はp_allシートにあります。
ご覧のとおり、逆行列は存在せず計算できません。ただ、CG法のような反復解法を使うとなんとなーく計算それっぽい値が計算できたりしてしまいますが、負圧になったりする(でも、符号を逆転させるとなんとなくそれっぽい値に見えるので、そうしてしまってる実装があります)ので、やめましょう。
ダミー粒子はあくまでダミー(固定粒子の粒子数密度を計算するだめだけに存在する)粒子であり、圧力を含むいかなる物理量もそこには存在しませんので(※圧力が0ではありません)。
なので、ダミー粒子は無効粒子と同じ扱いで、ポアソン方程式の計算から除外する必要があります。
具体的には、そもそも方程式系に含めないのがいいです。ダミー粒子がm個あるなら方程式系の次数をm減らす感じで。具体的にはp_nodummyシートみたいな感じになります。

さて、しかしこの状態でもまだ逆行列は計算できません。
そう、自由表面粒子をそのまま計算に含めているからです。
これも、反復解法に突っ込むとなんとなーく計算すると値が出てくるので、あとは境界条件によって、最後に計算が終わった後に自由表面粒子だけ圧力を0にするという実装をよく見ます。
しかしこれは間違っています。これは粒子法に限らず「境界条件とは何か」が分かっていればダメなのが分かるはずなんですが、あまりその辺の基礎を勉強してこなかった方(自分の観測範囲だと、ゲーム屋さんやCG屋さん)は間違えやすいポイントのようです。
この自由表面粒子のように、求めたい未知変数が既知値であるような境界条件(ディレクレ境界条件)では、「計算してみると未知変数の値が指定の境界条件値に合った」ようにしなければいけません。自由表面粒子で言えば、最終的に未知変数(圧力)値が0になるようにするということです。そのために注意すべき点は3点あります。

  • まず生成項\mathbf{b}ですが、その物理的意味は知っての通り、直接的には粒子数密度の時間変化量のことで、遡れば、圧力勾配による粒子数密度の修正量に起因します。自由表面粒子はその判定方法の通り、粒子数密度が他に比べて極端に小さく粒子数密度の修正を諦めた粒子なのですから、当然、これに起因する生成項は強制的に0であるべきです。
  • そうなると、係数行列Aですが、
    • 自由表面粒子に対応する行では、生成項が0なので、非対角成分を0にします(対角成分は0以外なら何でもいいです)。そうすることで、未知変数(圧力)値は0以外の値を取り得ません。
    • それ以外(内部粒子)の行はそのままです。ここで「自由表面粒子の圧力は0なんだから、列から消してもいいじゃん」と単純に消してしまってる実装も見かけましたが、それは間違いです(後述します)。

この3点に気をつけて実装してみるとようやく正しい計算ができるようになりました

これで正しい計算はできるようになりましたが、自由表面粒子の行は非対角成分が0なので、逆行列は明らかであって、そもそも未知変数ではないので無駄に計算しています。
ではどうやったら省略できるか?
「計算しなくていいなら単純に消してしまえ」と、ダミー粒子・無効粒子のように扱うと結果が合わなくなります。この理由は簡単で、対角項の値が違ってくるからです。
確かに自由表面粒子では未知数が0なので、内部粒子の行において自由表面粒子から受ける影響量(自由表面粒子に対応する列W(r_{ij})との積)は0になり省略可能ですが、対角項の総和計算\displaystyle \sum W\left(r_{ij}\right)には含めておかなければならないからです。
なので、正しく計算を省略するには

  • 自由表面粒子の行は方程式系には含めない
  • 自由表面粒子の列も方程式系には含めないが、ただし内部粒子の行の対角項の計算する時には含める

というようにするのが正しい手順です。

以上、最終的にはp_lastシートのようになります。擬似コードで書くなら

i = 0;
foreach(particle)
{
	if((particle.type != ghost) && (particle.type != dummy)&& (particle.type != surface)) // 無効粒子とダミー粒子と自由表面粒子の行は完全に除外
	{
		A(i, i) = 0;
		foreach(neighbor)
		{
			if(neighbor.type != ghost) && (neighbor.type != dummy)) // 無効粒子とダミー粒子の列は完全に除外
			{
				r = distance(particle, neighbor)
				a_ij = W(r)
				if(neighbor.type != surface) // 自由表面粒子の列は除外
				{
					A(i, j) = a_ij;
				}
				A(i, i) -= a_ij; // ただし対角項には含める
			}
		}
		i++;
	}
}

となるはずです。

「実装してみたけどなんか圧力の計算がおかしいっぽい」という人はこれを疑ってみるのがよいと思います。