「平滑化」粒子法流体力学（SPH法）の基礎知識

この2部構成のブログシリーズでは、粒子法流体力学（SPH: Smoothed Particle Hydrodynamics）の基本を理解し、より伝統的な有限体積（FV）数値計算法と比較した利点と欠点を説明し、nanoFluidX（現Altair CFD）でのSPHの実践について少し述べます。

煩雑な文献の入れ子参照を避けるため、この入門解説の作成に使用された、詳細な洞察を与えてくれる文献のいくつかを、興味のある読者に紹介したいと思います。始めに紹介するのは、Price [1] [2]によって書かれた2つの論文です。この論文では、SPH方程式を導き出し、この手法の興味深い点を非常にわかりやすく説明しています。Adamiの論文[3]は、nanoFluidXの特定のSPHプログラミング実装に興味がある人には良い参考文献です。乱流の挙動についてもっと知りたい人は、Adamiの別の論文[4]を読むことをお勧めします。その他、この手法の創始者であるLucy [5]やGingold and Monaghan [6]の代表的な著作も良い読み物です。

SPH法とは

SPH法の考え方は非常にシンプルで、特性の担い手である離散粒子群から連続場をどのように再構築するかというところから始まります。この問題は3次元であり、問題の解は時間に独立でなければなりません。よく考えてみると、これは3次元の補間問題であり、散在するデータを補間して、空間と時間のどの点においても連続場を再構成できるようにする問題であると言えます。

ここで、同じ問題を2次元で扱ってみると、3次元と同じように理解できるでしょう。下の図1は、SPHが2次元でどのように機能するかを示しています。各粒子をクローネッカーデルタ関数（鋭い信号関数）として表現するのではなく、ガウス関数に似たものを導入して粒子表現を平滑化しているので、「平滑化」粒子法流体力学（“smoothed“ particle hydrodynamics）という名前になっています。各粒子を平滑化する長さを平滑化長、ガウス関数に似た関数をカーネル関数と呼びます。平滑化された粒子を積み重ねると、バルクでは一定の連続した関数が再構成されることがわかりますが、これはまさに私たちが最初に目指したことです。

図1：異なるガウス様曲線（カーネル関数）を用いた2次元の連続関数の再構成

ここでは、2次元の表現にこだわって、いくつかの点を指摘することにします。まず、カーネル関数とスムージングの長さを変えると、再構成されたフィールドの結果に大きな影響を与えることがあることに注意してください。次に、軸に沿った粒子分布が再構成の結果に大きな影響を与えることがすぐに結論付けられます。3つ目は、図では分かりにくいのですが、粒子分布が特性勾配に影響を与えるということです。以上のことから、SPHでは、粒子分布とカーネル関数の特性が、正確な場の再構成を得るために重要であることがわかります。

3Dでも原理は変わりません。本当に変わるのは、3次元のSPH粒子を2次元の平面/紙で便利に表現できるようになることだけです。

ここで、それぞれの粒子に対して、質量、圧力、速度、密度、体積といった実際の特性を割り当てていきます。また、平滑化長半径で定義される体積の積分の合計が1になるようなカーネル関数にしてみましょう。もう一つ、カーネル関数が満たすべき条件は、平滑化長の半径で囲まれた隣接する粒子が複数あることです。これは、カーネル関数の体積積分と圧力の局所的な値を掛け合わせることで、圧力などのあらゆる特性を再構成できることを意味します。これを式にすると次のようになります。

粒子法流体力学（SPH法）の基礎知識

ここで、p̃は再構成された圧力、rは位置ベクトル、r¹ は平滑化長半径内の点の位置ベクトル、hは平滑化長、Wはカーネル関数です。これをSPH法で離散化すると

粒子法流体力学（SPH法）の基礎知識

ここで、インデックスはいわゆるオーナー粒子（場を再構成している粒子）、インデックスはいわゆる近傍粒子を表します。離散化された式では、体積積分の近似であることを考慮し、近傍粒子の体積である , を含んでいることに注意してください。このことは、カーネル和の単位が1/体積であることも意味しています。以上のような原理で、空間と時間の任意の点における速度、密度、温度など、あらゆる物性を再構成することができます。

結論から言うと、SPHは特定のカーネル関数を使った変数の3次元補間に基づく手法です。隣接する粒子がオーナー粒子に与える影響を計量し、これを領域全体で繰り返し、各粒子はある時点でオーナー粒子となります。

SPH法ではどのように物事が動いているか？

以上、SPHの基本原理を説明しましたが、勾配がどのように形成されるのか、言い換えれば、何が物事を動かすのかについては何も述べませんでした。これは実は簡単なことで、どんな変数の勾配も通常の再構成と同じように和をとることで表現できることがわかります。なお、カーネル関数はほとんどの場合多項式であり、その微分をとることは解析的にできるので、コーディングには非常に都合がよいのです。圧力勾配を離散形式で表すと次のようになります。

粒子法流体力学（SPH法）の基礎知識

簡単でしょう？でも、圧力差は最初はどこから来るのでしょうか？その答えはとても直観的で、粒子の分布に関係しています。だから、その重要性を強調したのです。

すなわち、SPH法の最も一般的なアプローチは、弱圧縮法です。この近似法は、通常、非圧縮性の流体をシミュレーションしますが、人工的に元の密度値の1%程度圧縮させることを含みます。これにより、流体中の音速を効果的に減少させ、時間ステップサイズに関する数値条件を緩和することができます。

さて、場の再構成には粒子分布が重要であることは既に述べました。密度場を再構成する場合、3次元的に一様な粒子分布（一様密度）からスタートすることを想像してみましょう。ここで、粒子の一部をごくわずかに一方向に変位させると、局所積分の合計値も増加します。これにより、局所的に密度の高い領域が生じ、領域の体積が一定であると仮定すると、反対側に希薄化領域が生じます。そして、SPHのTait方程式（液体の密度を圧力に関連付けるために使用される方程式で、状態方程式の一種）としても知られている準非圧縮性状態方程式を通して、密度変化を圧力変化に結びつけることができます。

粒子法流体力学（SPH法）の基礎知識

ここで、数値開始圧力はいわゆるバックグラウンド圧力、剛性指数は和に基づく所有粒子の現在の密度であり、流体密度のデフォルト値であることの詳細な説明は、冒頭で指定した文献を参照してください。

そこで、実際にはまず密度場を再構成し、その密度揺らぎをもとに、上のTait方程式を使って各粒子の圧力を計算する、という流れになっています。そして、前述のように圧力場の勾配を求め、流れの駆動力とするのです。

最後に、少し前置きが長くなりましたが、実際の粒子加速度の形を示すと、これは単にナビエストークス方程式を陽に解いたものです。

粒子法流体力学（SPH法）の基礎知識

というのは、この形で書くと、単純にニュートンの第2法則になります。

粒子法流体力学（SPH法）の基礎知識

上記の加速度の形式は、先に示したオリジナルの方程式に対して右辺が若干異なっていますが、これは力の対称性を確保し、微分の精度を高めるために導入した細かい修正に過ぎないことに注意してほしいと思います。詳しくは参考文献を参照してください。

もちろん、この加速度を時間で積分して速度を求め、さらに粒子の位置をずらすことで、粒子（密度）分布が変化し、新たな状態での圧力更新と流れの駆動力となり、これが時間変化にそって繰り返されます。

次回は、このシリーズのパート2として、有限体積数値法と比較したSPHのメリットとデメリットを検証します。nanoFluidXの詳細については、お問い合わせいただくか、製品概要ページをご覧ください。

この寄稿は、FluiDyna GmbHのプロダクトマネージャであるMilos Stanic博士によって書かれました。nanoFluidXは、ドイツのFluiDyna GmbHによって開発された最先端の産業用SPHコードです。AltairはnFXソフトウェアの世界的な総代理店で、詳細はこちらでご覧いただけます。

続編：平滑化粒子流体力学（SPH）法と有限体積数値計算法の比較

参考文献

[1]	D. J. Price, “Smoothed particle hydrodynamics and magnetohydrodynamics,” Journal of Computational Physics , vol. 231, pp. 759-794, 2012.
[2]	D. J. Price, “Smoothed Particle Hydrodynamics: Things I wish my mother taught me,” 11 2011.
[3]	S. Adami, Modeling and Simulation of Multiphase Phenomena with Smoothed Particle Hydrodynamics, Lehrstuhl für Aerodynamik und Strömungsmechanik, Technische Universität München, 2014.
[4]	S. Adami, X. Hu and N. Adams, “Simulating three-dimensional turbulence with SPH,” in Proceedings of the Summer Program 2012, Stanford Univeristy, Paolo Alto, 2012.
[5]	L. Lucy, “A numerical approach to the testing of the fission hypothesis,” Astronomical Journal, vol. 82, pp. 1013-1024, #dec# 1977.
[6]	R. Gingold and J. Monaghan, “Smoothed particle hydrodynamics – Theory and application to non-spherical stars,”Monthly Notices of the Royal Astronomical Society, vol. 181, pp. 375-389, #nov# 1977.
[7]	A. Kajzer, J. Pozorski and K. Szewc, “Large-eddy simulations of 3D Taylor-Green vortex: comparison of Smoothed Particle Hydrodynamics, Lattice Boltzmann and Finite Volume methods,” Journal of Physics: Conference Series, vol. 530, p. 012019, 2014.