平滑化粒子流体力学

平滑化粒子流体力学（SPH ）は、固体力学や流体の流れなどの連続体媒体の力学をシミュレートするために使用される計算手法です。これは、1977年にギンゴールドとモナハン^[2]およびルーシー^{[3]によって、当初は天体物理学の問題のために開発されました。これは、}天体物理学、弾道学、火山学、海洋学を含む多くの研究分野で使用されています。これはメッシュフリーの ラグランジアン法（座標が流体とともに移動する）であり、この方法の解像度は、密度などの変数に関して簡単に調整できます。

• SPH は構造上、メッシュフリー法であるため、自由表面の流れや大きな境界変位などの複雑な境界ダイナミクスが支配的な問題をシミュレートするのに最適です。
• メッシュがないため、多くのコアを持つアーキテクチャであっても、モデルの実装と並列化が大幅に簡素化されます。^{[4] [5]}

• SPH は、モデリング物理学で説明されているように、さまざまな分野に簡単に拡張でき、他のモデルとハイブリッド化できます。
• 弱圧縮性 SPH のセクションで説明したように、この方法は優れた保存特性を備えています。
• 関心のある測定基準が流体の密度（例えば、密度変動の確率密度関数）に関連している場合、SPHシミュレーションの粒子数あたりの計算コストは​​、グリッドベースのシミュレーションのセル数あたりのコストよりも大幅に低くなります。^[6] これは、SPHでは解像度が問題の場所に置かれるためです。

• ^{SPHでは、入口と出口[7]}や壁^[8]などの境界条件を設定することは、グリッドベースの方法よりも困難です。実際、「境界条件の扱いは、確かにSPH法の最も難しい技術的ポイントの1つです」と言われています。^[9] この課題は、SPHでは境界近くの粒子が時間とともに変化するためです。^[10] それにもかかわらず、SPHの壁境界条件は利用可能です。^{[8] [10] [11]}
• 関心のある測定基準が（直接）密度に関連していない場合（例えば、運動エネルギースペクトル）、粒子数あたりのSPHシミュレーションの計算コストは​​、セル数あたりのグリッドベースシミュレーションのコストよりも大幅に大きくなります。^[6] したがって、並列処理の高速化の問題を無視すると、一定密度の流れ（例えば、外部空気力学）のシミュレーションは、SPHよりもグリッドベースの方法の方が効率的です。

平滑化粒子流体力学は、流体の動きをモデル化するためにもますます使用されるようになっています。これは、従来のグリッドベースの手法に比べていくつかの利点があるためです。まず、粒子自体が質量を表すため、SPH では余分な計算を行わずに質量が保存されることが保証されます。次に、SPH では、線形方程式系を解くのではなく、隣接する粒子の重み付けされた寄与から圧力を計算します。最後に、流体の境界を追跡する必要があるグリッドベースの手法とは異なり、SPH では、粒子が密度の高い流体 (通常は水) を表し、空のスペースが軽い流体 (通常は空気) を表すため、2 相相互作用流体の自由表面を直接作成します。これらの理由から、SPH を使用して流体の動きをリアルタイムでシミュレートすることが可能です。ただし、グリッドベースの手法と SPH の手法のどちらも、メタボールやマーチング キューブ、ポイント スプラッティング、または「カーペット」視覚化などのポリゴン化手法を使用して、レンダリング可能な自由表面ジオメトリを生成する必要があります。ガス力学の場合、カーネル関数自体を使用してガス柱密度のレンダリングを生成する方が適切です (たとえば、SPLASH 視覚化パッケージで行われるように)。

グリッドベースの手法の欠点の 1 つは、同等の解像度のシミュレーションを生成するために多数のパーティクルが必要になることです。均一グリッドと SPH パーティクル手法の両方の一般的な実装では、レンダリングされない水ボリュームを埋めるために多数のボクセルまたはパーティクルが使用されます。ただし、高度なグリッドベースの手法、特にパーティクル メソッド (パーティクル レベル セットなど) と組み合わせた手法では、これらのシステムで非圧縮条件を適用するのが容易になるため、精度が大幅に向上します。流体シミュレーション用の SPH は、精度よりもインタラクティブ性が重視されるリアルタイム アニメーションやゲームでますます使用されています。

流体シミュレーションにおける SPH の最近の取り組みにより、パフォーマンス、精度、および適用領域が向上しました。

• B.ソレンサラー（2009年）は、より優れた非圧縮性制約を可能にする予測修正SPH（PCISPH）を開発しました^[12]
• M. Ihmsenらは2010年に、PCISPHの境界処理と適応型時間ステップを導入し、正確な剛体相互作用を実現した^[13]。
• K. Bodin et al., 2011は、標準的な状態圧力方程式を密度制約に置き換え、変分時間積分器を適用した^[14]
• R. Hoetzleinは2012年に、Fluids v.3で大規模なシーン向けの効率的なGPUベースのSPHを開発しました^[15]
• N. Akinciらは2012年に、流体力学に基づいた汎用的な境界処理と双方向SPH-剛体結合技術を導入した。このアプローチは、さまざまなタイプのSPHソルバーに適用できる^[16]
• M. Macklin et al., 2013は、より大きなタイムステップで、位置ベースダイナミクスフレームワーク内で非圧縮性の流れをシミュレートします^[17]
• N. Akinciらは2013年に、現実に観察される様々な興味深い物理的効果をシミュレートすることを可能にする多目的な表面張力と双方向の流体-固体接着技術を紹介した^[18]
• J. カイルとE. テレル（2013年）は、SPHをフルフィルム潤滑に適用した^[19]
• A. MahdaviとN. Talebbeydokhtiは2015年に、固体境界条件の実装と鋭角堰堤上の流れをシミュレートするためのハイブリッドアルゴリズムを提案した^[20]
• S. Tavakkolらは2016年に、粒子の水平方向と垂直方向のサイズを独立させ、曲面境界に沿って均一な質量分布を生成するcurvSPHを開発した^[21]
• W. KostorzとA. Esmail-Yakasは、2020年に、区分平面境界付近の正規化係数を評価するための一般的で効率的かつ単純な方法を提案した^[11]
• コラグロッシら（2019年）は、自由表面に近い円筒周りの流れを研究し、他の手法と比較した^[1]

平滑化粒子流体力学は、適応解像度、物理的に保存される量の数値的保存、および多くの桁にわたる現象をシミュレートする能力を備えているため、理論天体物理学の計算に最適です。^[22]

銀河形成、星形成、恒星衝突^[23] 、超新星^[24]、隕石衝突のシミュレーションは、この方法の天体物理学および宇宙論における幅広い用途の一部です。

SPHは、重力の影響を含む流体力学的流れをモデル化するために使用されます。放射伝達や磁場など、重要となる可能性のある他の天体物理学的プロセスを組み込むことは、天文学コミュニティで活発に研究されている分野であり、ある程度の成功を収めています。^{[25] [26]}

LiberskyとPetschek ^{[27] [28]} はSPHを固体力学に拡張しました。このアプリケーションにおけるSPHの主な利点は、グリッドベースの方法よりも大きな局所的な歪みを処理できることです。この機能は、金属成形、衝撃、亀裂成長、破壊、断片化など、固体力学の多くのアプリケーションで活用されています。

メッシュフリー法全般、特にSPHのもう一つの重要な利点は、この方法がメッシュフリーであるため、メッシュ依存性の問題が自然に回避されることです。特に、メッシュの位置合わせは亀裂を伴う問題に関係しており、カーネル関数の等方性サポートによりSPHでは回避されます。しかし、古典的なSPH定式化では、引張不安定性^[29] と一貫性の欠如^[30]の問題があります。 過去数年間、SPHソリューションの精度を向上させるためにさまざまな補正が導入され、LiuらによるRKPM ^[31]につながりました。Randles とLibersky ^[32] 、JohnsonとBeissel ^[33]は、 衝撃現象の研究で一貫性の問題を解決しようとしました。

Dykaら^{[34] [35]} とRandlesとLibersky ^[36]は SPHに応力点積分を導入し、Ted Belytschkoら^[37] は応力点技術が偽の特異モードによる不安定性を除去し、ラグランジュカーネルを使用することで引張不安定性を回避できることを示しました。SPH法の収束性を改善するための最近の研究は他にも数多く文献に見ることができます。

SPH の収束性と安定性に関する最近の理解の向上により、固体力学におけるより広範な応用が可能になりました。この方法のその他の応用例と開発例は次のとおりです。

• 金属成形シミュレーション^[38]
• ウィリアム・G・フーバーによる固体の衝撃破壊に対するSPHベースの方法SPAM（Smoothed Particle Applied Mechanics）^[39]
• 破壊と断片化のための修正SPH（SPH/MLSPH）。^[40]
• 固体中の衝撃波伝播のためのテイラーSPH（TSPH）。^[41]
• 一般化座標SPH（GSPH）は、粒子をデカルト座標系で不均一に割り当て、粒子が均一な間隔で整列する一般化座標系にマッピングして配置します。^[42]

平滑化粒子流体力学 (SPH) 法は、流体を粒子と呼ばれる個別の移動要素のセットに分割することによって機能します。粒子のラグランジュ特性により、速度を積分することで位置を設定できます。 ${\displaystyle i,j}$${\displaystyle \mathbf {r} _{i}}$${\displaystyle \mathbf {v} _{i}}$

${\displaystyle {\frac {\mathrm {d} {\boldsymbol {r}}_{i}}{\mathrm {d} t}}={\boldsymbol {v}}_{i}.}$

これらの粒子は、「平滑化長さ」と呼ばれる特性半径を持つカーネル関数を介して相互作用します。これは通常、方程式で として表されます。つまり、カーネルの範囲内にあるすべての粒子の関連プロパティを合計することで、任意の粒子の物理量を取得できます。カーネルは重み付け関数 として使用されます。これは 2 つの手順で理解できます。まず、任意のフィールドが との畳み込みとして記述されます。 ${\displaystyle h}$${\displaystyle W}$${\displaystyle A}$${\displaystyle W}$

${\displaystyle A({\boldsymbol {r}})=\int A\left({\boldsymbol {r^{\prime }}}\right)W(|{\boldsymbol {r}}-{\boldsymbol {r^{\prime }}}|,h)\,\mathrm {d} V\!\left({\boldsymbol {r'}}\right).}$

上記の近似を行う際の誤差は オーダーです。次に、積分は粒子上のリーマン和を使用して近似されます。 ${\displaystyle h^{2}}$

${\displaystyle A({\boldsymbol {r}})=\sum _{j}V_{j}A_{j}W(|{\boldsymbol {r}}-{\boldsymbol {r}}_{j}|,h),}$

ここで、 の合計にはシミュレーション内のすべての粒子が含まれます。は粒子の体積、は粒子の量の値、 は位置を表します。 たとえば、粒子の密度は次のように表すことができます。 ${\displaystyle j}$${\displaystyle V_{j}}$${\displaystyle j}$${\displaystyle A_{j}}$${\displaystyle A}$${\displaystyle j}$${\displaystyle {\boldsymbol {r}}}$${\displaystyle \rho _{i}}$${\displaystyle i}$

${\displaystyle \rho _{i}=\rho ({\boldsymbol {r}}_{i})=\sum _{j}m_{j}W_{ij},}$

ここで は粒子の質量と粒子の密度を表し、 は の短縮表記である。離散和で積分を近似する際に生じる誤差は、粒子の大きさ（つまりは空間次元）、および空間における粒子の配置に依存する。後者の影響はまだよくわかっていない。^[43]${\displaystyle m_{j}=\rho _{j}V_{j}}$${\displaystyle \rho _{j}}$${\displaystyle W_{ij}=W_{ji}}$${\displaystyle W(|{\boldsymbol {r}}_{i}-{\boldsymbol {r}}_{j}|,h)}$${\displaystyle h}$${\displaystyle V_{j}^{1/d}}$${\displaystyle d}$

一般的に使用されるカーネル関数には、ガウス関数、五次スプライン、ウェンドランドカーネルなどがあります。^[44]最後の2つのカーネルはコンパクトにサポートされており（有限距離離れた場所に小さな寄与があるガウスとは異なり）、サポートは に比例します。これには、遠くの粒子からの比較的小さな寄与を含めないことで計算労力を節約できるという利点があります。 ${\displaystyle C^{2}}$${\displaystyle h}$

スムージング長のサイズは空間的にも時間的にも固定できますが、これでは SPH の能力をフルに活用できません。各粒子に独自のスムージング長を割り当て、時間とともに変化させることで、シミュレーションの解像度をローカル条件に応じて自動的に適応させることができます。たとえば、多くの粒子が密集している非常に密度の高い領域では、スムージング長を比較的短くして、高い空間解像度を得ることができます。逆に、個々の粒子が離れており解像度が低い低密度領域では、スムージング長を増やして、関心領域の計算を最適化することができます。

質量が一定の粒​​子の場合、補間された密度を時間に関して 微分すると、${\displaystyle \rho _{i}}$

${\displaystyle {\frac {d\rho _{i}}{dt}}=\sum _{j}m_{j}\left({\boldsymbol {v}}_{i}-{\boldsymbol {v}}_{j}\right)\cdot \nabla W_{ij},}$

ここで はに関するの勾配です。この式をラグランジアン記述の連続方程式（物質導関数を使用）と比較すると、 ${\displaystyle \nabla W_{ij}=-\nabla W_{ji}}$${\displaystyle W_{ij}}$${\displaystyle {\boldsymbol {r}}_{i}}$

${\displaystyle {\frac {d\rho }{dt}}=-\rho \nabla \cdot {\boldsymbol {v}},}$

右辺が の近似であることは明らかです。したがって、離散発散演算子を次のように定義します。 ${\displaystyle -\rho \nabla \cdot \mathbf {v} }$

${\displaystyle \operatorname {D} _{i}\left\{{\boldsymbol {v}}_{j}\right\}=-{\frac {1}{\rho _{i}}}\sum _{j}m_{j}\left({\boldsymbol {v}}_{i}-{\boldsymbol {v}}_{j}\right)\cdot \nabla W_{ij}.}$

この演算子は、与えられた質量、位置、速度を持つ与えられた粒子セットに対して、粒子におけるの SPH 近似値を与えます。 ${\displaystyle \nabla \cdot \mathbf {v} }$${\displaystyle i}$${\displaystyle m_{j}}$${\displaystyle \left\{\mathbf {r} _{j}\right\}}$${\displaystyle \left\{\mathbf {v} _{j}\right\}}$

圧縮性非粘性流体に関するもう 1 つの重要な方程式は、運動量バランスに関する オイラー方程式です。

${\displaystyle {\frac {d{\boldsymbol {v}}}{dt}}=-{\frac {1}{\rho }}\nabla p+{\boldsymbol {g}}}$

連続性と同様に、離散勾配演算子を定義して次のように記述することが課題である。

${\displaystyle {\frac {d{\boldsymbol {v}}_{i}}{dt}}=-{\frac {1}{\rho }}\operatorname {\mathbf {G} } _{i}\left\{p_{j}\right\}+{\boldsymbol {g}}}$

一つの選択肢は

${\displaystyle \operatorname {\mathbf {G} } _{i}\left\{p_{j}\right\}=\rho _{i}\sum _{j}m_{j}\left({\frac {p_{i}}{\rho _{i}^{2}}}+{\frac {p_{j}}{\rho _{j}^{2}}}\right)\nabla W_{ij},}$

これは、上記の発散演算子と歪随伴であるという性質を持ち、

${\displaystyle \sum _{i}V_{i}{\boldsymbol {v}}_{i}\cdot \operatorname {\mathbf {G} } _{i}\left\{p_{j}\right\}=-\sum _{i}V_{i}p_{i}\operatorname {D} _{i}\left\{{\boldsymbol {v}}_{j}\right\},}$

これは連続体の離散バージョンである

${\displaystyle \int {\boldsymbol {v}}\cdot \operatorname {grad} p=-\int p\operatorname {div} \cdot {\boldsymbol {v}}.}$

この特性は優れた保全特性につながります。^[45]

また、この選択により対称発散演算子と反対称勾配が導かれることにも注意してください。オイラー方程式で圧力勾配を離散化する方法がいくつかありますが、上記の反対称形式が最もよく知られています。これは、線形運動量と角運動量の厳密な保存をサポートします。つまり、反対称性により、粒子ごとに作用する力は、有効方向の符号変更を含めて粒子ごとに作用する力に等しくなります。 ${\displaystyle i}$${\displaystyle j}$${\displaystyle j}$${\displaystyle i}$${\displaystyle \nabla W_{ij}=-\nabla W_{ji}}$

それにもかかわらず、数値的または物理的により優れた性能を発揮する可能性のある他の演算子が提案されている。例えば、これらの演算子の欠点の1つは、発散はゼロ次整合（つまり、定数ベクトル場に適用するとゼロになる）である一方、勾配は整合していないことが分かることである。この問題を回避するためにいくつかの手法が提案されており、再正規化演算子につながる（例えば^[46]を参照）。 ${\displaystyle \operatorname {D} }$${\displaystyle \operatorname {\mathbf {G} } }$

上記の SPH 支配方程式は、粒子系の ラグランジアンから始めて、最小作用原理から導くことができます。

${\displaystyle {\mathcal {L}}=\sum _{j}m_{j}\left({\tfrac {1}{2}}{\boldsymbol {v}}_{j}^{2}-e_{j}+{\boldsymbol {g}}\cdot {\boldsymbol {r}}_{j}\right)}$、

ここで、 は粒子固有の内部エネルギーです。変分力学のオイラー・ラグランジュ方程式は、各粒子について次のようになります。 ${\displaystyle e_{j}}$

${\displaystyle {\frac {\mathrm {d} }{\mathrm {d} t}}{\frac {\partial {\mathcal {L}}}{\partial {\boldsymbol {v}}_{i}}}={\frac {\partial {\mathcal {L}}}{\partial {\boldsymbol {r}}_{i}}}.}$

上記のラグランジアンに適用すると、次の運動量方程式が得られます。

${\displaystyle m_{i}{\frac {\mathrm {d} {\boldsymbol {v}}_{i}}{\mathrm {d} t}}=-\sum _{j}m_{j}{\frac {\partial e_{j}}{\partial {\boldsymbol {r}}_{i}}}+m_{i}{\boldsymbol {g}}=-\sum _{j}m_{j}{\frac {\partial e_{j}}{\partial \rho _{j}}}{\frac {\partial \rho _{j}}{\partial {\boldsymbol {r}}_{i}}}+m_{i}{\boldsymbol {g}}}$

ここで連鎖律が使われているのは、が に依存し、後者は粒子の位置に依存するからである。熱力学的性質を用いると、次のように書くことができる。 ${\displaystyle e_{j}}$${\displaystyle \rho _{j}}$${\displaystyle \mathrm {d} e=\left(p/\rho ^{2}\right)\mathrm {d} \rho }$

${\displaystyle m_{i}{\frac {\mathrm {d} {\boldsymbol {v}}_{i}}{\mathrm {d} t}}=-\sum _{j}m_{j}{\frac {p_{j}}{\rho _{j}^{2}}}{\frac {\partial \rho _{j}}{\partial {\boldsymbol {r}}_{i}}}+m_{i}{\boldsymbol {g}},}$

SPH密度補間を代入し、明示的に微分すると、 ${\displaystyle {\tfrac {\partial \rho _{j}}{\partial {\boldsymbol {r}}_{i}}}}$

${\displaystyle {\frac {\mathrm {d} {\boldsymbol {v}}_{i}}{\mathrm {d} t}}=-\sum _{j}m_{j}\left({\frac {p_{i}}{\rho _{i}^{2}}}+{\frac {p_{j}}{\rho _{j}^{2}}}\right)\nabla W_{ij}+{\boldsymbol {g}},}$

これはすでに述べたSPH運動量方程式であり、演算子を認識します。これにより、線形運動量が保存される理由が説明され、角運動量とエネルギーの保存も保存されます。^[47]${\displaystyle \operatorname {\mathbf {G} } }$

80 年代と 90 年代に行われた大型加速器内の点状粒子の数値積分に関する研究から、長期間にわたって正確な保存特性を持つ適切な時間積分器が開発されました。これらはシンプレクティック積分器と呼ばれます。SPH の文献で最もよく知られているのは、各粒子について次のように表されるリープフロッグ方式です。 ${\displaystyle i}$

${\displaystyle {\begin{aligned}{\boldsymbol {v}}_{i}^{n+1/2}&={\boldsymbol {v}}_{i}^{n}+{\boldsymbol {a}}_{i}^{n}{\frac {\Delta t}{2}},\\{\boldsymbol {r}}_{i}^{n+1}&={\boldsymbol {r}}_{i}^{n}+{\boldsymbol {v}}_{i}^{i+1/2}\Delta t,\\{\boldsymbol {v}}_{i}^{n+1}&={\boldsymbol {v}}_{i}^{n+1/2}+{\boldsymbol {a}}_{i}^{i+1}{\frac {\Delta t}{2}},\end{aligned}}}$

ここで、 は時間ステップ、上付き文字は時間の反復を表し、 は粒子の加速度であり、運動量方程式の右辺で与えられます。 ${\displaystyle \Delta t}$${\displaystyle {\boldsymbol {a}}_{i}}$

他にもシンプレクティック積分器は存在します（参考文献^[48]を参照）。多くの反復後に誤差が蓄積されるのを避けるために、高次の非シンプレクティック方式ではなく、シンプレクティック（低次でも）方式を使用することが推奨されます。

密度の積分については、あまり研究されていません (詳細については以下を参照してください)。

シンプレクティック スキームは保守的ですが明示的であるため、数値安定性には、Courant-Friedrichs-Lewy 条件 (以下を参照) に類似した安定性条件が必要です。

SPH畳み込みが境界の近く、つまりs · hよりも近くで実行される場合、積分サポートは切り捨てられます。実際、畳み込みが境界の影響を受ける場合、畳み込みは2つの積分に分割されます。

${\displaystyle A({\boldsymbol {r}})=\int _{\Omega ({\boldsymbol {r}})}A\left({\boldsymbol {r^{\prime }}}\right)W(|{\boldsymbol {r}}-{\boldsymbol {r^{\prime }}}|,h)d{\boldsymbol {r^{\prime }}}+\int _{B({\boldsymbol {r}})-\Omega ({\boldsymbol {r}})}A\left({\boldsymbol {r^{\prime }}}\right)W(|{\boldsymbol {r}}-{\boldsymbol {r^{\prime }}}|,h)d{\boldsymbol {r^{\prime }}},}$

ここで、 B( r )は半径s · hでrを中心とするコンパクトな支持球であり、Ω( r )は計算領域内のコンパクトな支持球の部分、Ω ∩ B( r )を表す。したがって、SPHにおける境界条件の課しは、右辺の第2積分の近似に完全に基づいている。もちろん、同じことが微分演算子の計算にも適用できる。

${\displaystyle \nabla A({\boldsymbol {r}})=\int _{\Omega ({\boldsymbol {r}})}A\left({\boldsymbol {r^{\prime }}}\right)\nabla W({\boldsymbol {r}}-{\boldsymbol {r^{\prime }}},h)d{\boldsymbol {r^{\prime }}}+\int _{B({\boldsymbol {r}})-\Omega ({\boldsymbol {r}})}A\left({\boldsymbol {r^{\prime }}}\right)\nabla W({\boldsymbol {r}}-{\boldsymbol {r^{\prime }}},h)d{\boldsymbol {r^{\prime }}}.}$

これまで、SPH で境界をモデル化するためのいくつかの手法が導入されてきました。

最も単純な境界モデルは積分を無視するものであり、

${\displaystyle \int _{B({\boldsymbol {r}})-\Omega ({\boldsymbol {r}})}A\left({\boldsymbol {r^{\prime }}}\right)\nabla W({\boldsymbol {r}}-{\boldsymbol {r^{\prime }}},h)d{\boldsymbol {r^{\prime }}}\simeq {\boldsymbol {0}},}$

バルク相互作用のみが考慮されるように、

${\displaystyle \nabla A_{i}=\sum _{j\in \Omega _{i}}V_{j}A_{j}\nabla W_{ij}.}$

これは、単相シミュレーションで自由表面を考慮する場合によく使われるアプローチです。^[49]

この境界条件の主な利点は、その明らかな単純さです。しかし、この境界技術を適用する際には、いくつかの一貫性の問題を考慮する必要があります。^[49]これは、実際にはその潜在的な応用に大きな制限となります。

おそらくSPHで境界条件を課す最も一般的な方法論、あるいは少なくとも最も伝統的な方法は、流体拡張技術です。この技術は、境界を横切るコンパクトサポートにいわゆるゴースト粒子を配置し、その場の値を都合よく課すことに基づいています。^[50]

この線に沿って、積分無視法は流体拡張の特殊なケースとして考えることができ、その場合、フィールドAは計算領域の外側では消滅します。

この方法論の主な利点は、境界寄与がバルク相互作用の一部として計算されるという単純さである。また、この方法論は文献で深く分析されている。^{[51] [50] [52]}

一方、切り捨てられた領域にゴースト粒子を配置することは簡単な作業ではなく、複雑な境界形状をモデル化することは面倒になります。空の領域にゴースト粒子を配置するための最も一般的な2つのアプローチは、ミラー粒子^[53]と固定粒子^{[50]です。}

最新の境界手法は境界積分法である。^[54]この方法では、空の体積積分が表面積分と繰り込みに置き換えられる。

${\displaystyle \nabla A_{i}={\frac {1}{\gamma _{i}}}\left(\sum _{j\in \Omega _{i}}V_{j}A_{j}\nabla W_{ij}+\sum _{j\in \partial \Omega _{i}}S_{j}A_{j}{\boldsymbol {n}}_{j}W_{ij}\right),}$

${\displaystyle \gamma _{i}=\sum _{j\in \Omega _{i}}V_{j}W_{ij},}$

ここでnj_は一般的なj番目の境界要素の法線である。表面項は半解析的表現を考慮して解くこともできる。^[54]

密度を決定する別の方法は、SPH平滑化演算子自体に基づいています。したがって、密度はSPH補間を利用して粒子分布から推定されます。カーネル切り捨てによる自由表面での望ましくない誤差を克服するために、密度の定式化を再び時間で積分することができます。 ^[54]

流体力学における弱圧縮性SPHは、圧縮性流体のナビエ・ストークス方程式またはオイラー方程式の離散化に基づいています。システムを閉じるために、圧力と密度をリンクする適切な状態方程式が利用されます。一般的に、SPHではいわゆるコール方程式^{[55] （誤って「}テイト方程式」 と呼ばれることもあります）が使用されます。 ${\displaystyle p}$${\displaystyle \rho }$

${\displaystyle p={\frac {\rho _{0}c^{2}}{\gamma }}\left(\left({\frac {\rho }{\rho _{0}}}\right)^{\gamma }-1\right)+p_{0},}$

ここで、は基準密度、は音速です。水の場合、が一般的に使用されます。負の圧力値を避けるために、 背景圧力が追加されます。${\displaystyle \rho _{0}}$${\displaystyle c}$${\displaystyle \gamma =7}$${\displaystyle p_{0}}$

水などの実際のほぼ非圧縮性の流体は、 のオーダーの非常に高い音速を特徴とします。したがって、圧力情報は実際の体積流に比べて速く伝わり、マッハ数は非常に小さくなります。運動量方程式から次の関係が導かれます。 ${\displaystyle 10^{3}\mathrm {m/s} }$${\displaystyle M}$

${\displaystyle {\frac {\delta \rho }{\rho _{0}}}\approx {\frac {|{\boldsymbol {v}}|^{2}}{c^{2}}}=M^{2}}$

ここで、 は密度の変化と速度ベクトルです。実際には、時間積分スキームで時間ステップが小さすぎることを避けるために、実際の値よりも小さい c の値が採用されます。一般的に、密度の変化が 1% 未満になるように数値的な音速が採用されます。これはいわゆる弱圧縮性仮定です。これはマッハ数が0.1 未満であることに対応し、次のことを意味します。 ${\displaystyle \rho }$${\displaystyle v}$

${\displaystyle c=10v_{\text{max}}}$

ここで、最大速度は、例えばトリチェリの法則や推測によって推定する必要がある。密度の変化は小さいため、線形状態方程式を採用することができる。^[56]${\displaystyle v_{\text{max}}}$

${\displaystyle p=c^{2}\left(\rho -\rho _{0}\right)}$

通常、弱圧縮スキームは、圧力場と密度場における高周波スプリアスノイズの影響を受けます。 ^[57] この現象は、音波の非線形相互作用と、スキームが時間的に明示的で空間的に中心にあるという事実によって引き起こされます。^[58]

長年にわたり、この問題を解決するためにいくつかの手法が提案されてきました。それらは 3 つのグループに分類できます。

1. 密度フィルタを採用した方式、
2. 連続方程式に拡散項を加えたモデル、
3. リーマンソルバーを使用して粒子の相互作用をモデル化するスキーム。

最初のグループの方式では、密度場に対して直接フィルターを適用して、不要な数値ノイズを除去します。最もよく使用されるフィルターは、MLS (移動最小二乗法) と Shepard フィルター^[57]で 、各時間ステップまたは n 時間ステップごとに適用できます。フィルター処理手順を頻繁に使用するほど、より規則的な密度場と圧力場が得られます。一方で、これは計算コストの増加につながります。長時間のシミュレーションでは、フィルター処理手順を使用すると、静水圧成分が乱れ、流体の全体体積と密度場の間に不整合が生じる可能性があります。さらに、動的自由表面境界条件の強制が保証されません。

密度と圧力の場を滑らかにする別の方法は、連続方程式（グループ2）内に拡散項を追加することです。

${\displaystyle {\displaystyle {\frac {d\rho _{i}}{dt}}=\sum _{j}m_{j}\left({\boldsymbol {v}}_{i}-{\boldsymbol {v}}_{j}\right)\cdot \nabla W_{ij}+{\mathcal {D}}_{i}(\rho ),}}$

このようなアプローチを採用した最初のスキームは、Ferrari ^[59] とMolteni ^{[56]によって記述され、拡散項は密度場のラプラシアンとしてモデル化されました。同様のアプローチは、FatehiとManzari [60]} でも使用されました。

Antuonoら ^[61]では、自由表面近くのいくつかの矛盾を除去するために、 Molteni ^[56] の拡散項の補正が提案された。この場合、採用された拡散項は密度場の高次微分演算子と同等である。 ^[62] この方式はδ-SPHと呼ばれ、拡散のないSPHのすべての保存特性（例えば、線形および角運動量、全エネルギー、 ^[63]を参照 ）と、密度場と圧力場の滑らかで規則的な表現を維持する。

3番目のグループには、リーマンソルバーで得られた数値フラックスを使用して粒子相互作用をモデル化するSPHスキームがあります。^{[64] [65] [66]}

リーマンソルバーに基づくSPH法では、粒子から粒子 を指す単位ベクトルに沿って粒子間リーマン問題が構築されます 。このリーマン問題では、初期の左および右の状態は それぞれ粒子および上にあります。および状態は ${\displaystyle \mathbf {e} _{ij}=-\mathbf {r} _{ij}/r_{ij}}$${\displaystyle i}$${\displaystyle j}$${\displaystyle i}$${\displaystyle j}$${\displaystyle L}$${\displaystyle R}$

$${\displaystyle {\begin{cases}(\rho _{L},U_{L},P_{L})=(\rho _{i},\mathbf {v} _{i}\cdot \mathbf {e} _{ij},P_{i})\\(\rho _{R},U_{R},P_{R})=(\rho _{j},\mathbf {v} _{j}\cdot \mathbf {e} _{ij},P_{j}).\end{cases}}}$$

リーマン問題の解は、不連続面から3つの波を発することになります。2つの波は衝撃波または希薄波で、最小または最大の波速で伝わります。中間の波は常に接触不連続面であり、およびで表される2つの中間状態を分離します。中間状態が および を満たすと仮定すると、滑らかな流れまたは中程度の衝撃のみの線形リーマンソルバーは次のように記述できます。 ${\displaystyle (\rho _{L}^{\ast },U_{L}^{\ast },P_{L}^{\ast })}$${\displaystyle (\rho _{R}^{\ast },U_{R}^{\ast },P_{R}^{\ast })}$${\displaystyle U_{L}^{\ast }=U_{R}^{\ast }=U^{\ast }}$${\displaystyle P_{L}^{\ast }=P_{R}^{\ast }=P^{\ast }}$

$${\displaystyle {\begin{cases}U^{\ast }={\overline {U}}+{\frac {1}{2}}{\frac {(P_{L}-P_{R})}{{\bar {\rho }}c_{0}}}\\P^{\ast }={\overline {P}}+{\frac {1}{2}}{\bar {\rho }}c_{0}{(U_{L}-U_{R})},\end{cases}}}$$

ここで、 および は粒子間平均である。リーマン問題の解、すなわちおよび を用いて、SPH法の離散化は次のように表される。 ${\displaystyle {\overline {U}}=(U_{L}+U_{R})/2}$${\displaystyle {\overline {P}}=(P_{L}+P_{R})/2}$${\displaystyle U^{\ast }}$${\displaystyle P^{\ast }}$

$${\displaystyle {\frac {d\rho _{i}}{dt}}=2\rho _{i}\sum _{j}{\frac {m_{j}}{\rho _{j}}}(\mathbf {v} _{i}-\mathbf {v} ^{\ast })\cdot \nabla _{i}W_{ij},}$$

$${\displaystyle {\frac {d\mathbf {v} _{i}}{dt}}=-2\sum _{j}m_{j}\left({\frac {P^{\ast }}{\rho _{i}\rho _{j}}}\right)\nabla _{i}W_{ij}.}$$

ここで 、これは粒子間の平均速度と圧力が、リーマン問題の解によって単純に置き換えられることを示しています。両者を比較すると、粒子間の平均からの中間速度と圧力が、それぞれ暗黙の散逸、つまり密度正規化と数値粘性に相当することがわかります。 ${\displaystyle \mathbf {v} ^{\ast }=U^{\ast }\mathbf {e} _{ij}+({\overline {\mathbf {v} }}_{ij}-{\overline {U}}\mathbf {e} _{ij})}$

上記の離散化は非常に散逸的であるため、中間圧力を制限することによって導入される暗黙の数値散逸を減らすためにリミッターを適用するのが簡単な修正である ^[67]

$${\displaystyle P^{\ast }={\overline {P}}+{\frac {1}{2}}\beta {\overline {\rho }}{(U_{L}-U_{R})},}$$

ここでリミッタは次のように定義される。

$${\displaystyle \beta =\min {\big (}\eta \max(U_{L}-U_{R},0),{\overline {c}}{\big )}.}$$

流体が膨張波の作用を受けているとき、すなわち では散逸が起こらないことを が保証し、流体が圧縮波の作用を受けているとき、すなわちではパラメータが散逸を調整するために使用されることに注意してください。数値実験により、 は一般に効果的であることがわかりました。また、中間速度によって導入される散逸は制限されないことにも注意してください。 ${\displaystyle \beta }$${\displaystyle U_{L}<U_{R}}$${\displaystyle \eta }$${\displaystyle U_{L}\geq U_{R}}$${\displaystyle \eta =3}$

一般に、流体力学の流れを記述するには、ナビエ・ストークス方程式の粘性をモデル化するために拡散過程を簡便に扱う必要があります。これにはラプラシアン微分演算子が関係するため、特別な考慮が必要です。直接計算では満足のいく結果が得られないため、拡散をモデル化するいくつかのアプローチが提案されています。

• 人工粘性

モナハンとギンゴールド ^[68]によって導入された 人工粘性は、高マッハ数流体の流れを扱うために使用されました。

${\displaystyle \Pi _{ij}={\begin{cases}{\dfrac {-\alpha {\bar {c}}_{ij}\phi _{ij}+\beta \phi _{ij}^{2}}{{\bar {\rho }}_{ij}}}&\quad {\boldsymbol {v}}_{ij}\cdot {\boldsymbol {r}}_{ij}<0\\0&\quad {\boldsymbol {v}}_{ij}\cdot {\boldsymbol {r}}_{ij}\geq 0\end{cases}}}$

ここで、は体積粘度を制御し、ノイマン・リヒトマイヤーの人工粘度と同様に機能します。は次のように定義されます。 ${\displaystyle \alpha }$${\displaystyle \beta }$${\displaystyle \phi _{ij}}$

${\displaystyle \phi _{ij}={\frac {h{\boldsymbol {v}}_{ij}\cdot {\boldsymbol {r}}_{ij}}{\Vert {\boldsymbol {r}}_{ij}\Vert ^{2}+\eta _{h}^{2}}},}$

ここで、η _{h は}hの小さな分数（例：0.01 h）であり、近距離で数値が無限大になるのを防ぎます。

人工粘性は、一般的な流れのシミュレーションの全体的な安定性を向上させることも示されています。したがって、次の形式で非粘性問題に適用されます。

${\displaystyle \Pi _{ij}=\alpha hc{\frac {{\boldsymbol {v}}_{ij}\cdot {\boldsymbol {r}}_{ij}}{\Vert {\boldsymbol {r}}_{ij}\Vert ^{2}+\eta _{h}^{2}}}.}$

このアプローチにより、非粘性シミュレーションを安定化できるだけでなく、物理的粘性をモデル化することも可能です。

${\displaystyle \alpha hc=2(n+2){\frac {\mu }{\rho }}}$

は上記の式に代入されます。ここで、はモデルの空間次元の数です。このアプローチでは、バルク粘性が導入されます。 ${\displaystyle n}$${\displaystyle \zeta ={\frac {5}{3}}\mu }$

• モリス

低レイノルズ数に対しては、モリス^[69]による粘性モデルが 提案された。

${\displaystyle [\nu \Delta {\boldsymbol {v}}]_{ij}=2\nu {\frac {m_{j}}{\rho _{j}}}\,{\frac {{\boldsymbol {r}}_{ij}\cdot \nabla w_{h,ij}}{\Vert {\boldsymbol {r}}_{ij}\Vert ^{2}+\eta _{h}^{2}}}\,{\boldsymbol {v}}_{ij}.}$

• ロシャオ

• 表面張力
• 熱伝達
• 乱気流

天体物理学では、純粋な流体力学に加えて自己重力をモデル化することがしばしば望まれます。SPHの粒子ベースの性質は、粒子ベースの重力ソルバー、例えばtree Gravity Code、^{[70] 、} particle mesh、particle-particle-meshと組み合わせるのに最適です。

固体力学の支配方程式を離散化するために、 まず、剛体回転を再現するための 補正行列^{[71] [72]が導入される。}${\displaystyle \mathbb {B} ^{0}}$

どこ

$${\displaystyle \nabla _{a}^{0}W_{a}={\frac {\partial W\left(|\mathbf {r} _{ab}^{0}|,h\right)}{\partial |\mathbf {r} _{ab}^{0}|}}\mathbf {e} _{ab}^{0}}$$

は、初期参照構成で評価されたカーネル関数の勾配を表します。下付き文字とが固体粒子を表すために使用され、平滑化の長さは流体方程式の離散化の場合と同じであることに注意してください。 ${\displaystyle a}$${\displaystyle b}$${\displaystyle h}$

初期構成を基準として、固体密度は次のように直接評価される。

ここで、 は変形テンソルのヤコビ行列式です。 ${\displaystyle J=\det(\mathbb {F} )}$${\displaystyle \mathbb {F} }$

運動量方程式を次の形で離散化することができる。

ここで、粒子間の平均第一ピオラ・キルヒホッフ応力は 次のように定義される。 ${\displaystyle {\tilde {\mathbb {P} }}}$

また、およびは、固体粒子に作用する流体圧力と粘性力にそれぞれ対応します。 ${\displaystyle \mathbf {f} _{a}^{F:p}}$${\displaystyle \mathbf {f} _{a}^{F:v}}$${\displaystyle a}$

流体構造連成では、周囲の固体構造が流体の移動境界として振る舞い、流体構造界面には滑りのない境界条件が課せられる。隣接する固体粒子の存在により流体粒子に作用する相互作用力とは次の式で得られる ^[73]。${\displaystyle \mathbf {f} _{i}^{S:p}}$${\displaystyle \mathbf {f} _{i}^{S:v}}$${\displaystyle i}$${\displaystyle a}$

そして

ここで、虚圧力と虚速度は次のように定義される。 ${\displaystyle p_{a}^{d}}$${\displaystyle \mathbf {v} _{a}^{d}}$

ここで、は固体構造の表面法線方向を表し、仮想粒子密度は状態方程式を通じて計算されます。 ${\displaystyle \mathbf {n} ^{S}}$${\displaystyle \rho _{a}^{d}}$

したがって、固体粒子に作用する相互作用力とは次の式で与えられる。 ${\displaystyle \mathbf {f} _{a}^{F:p}}$${\displaystyle \mathbf {f} _{a}^{F:v}}$${\displaystyle a}$

そして

カーネル関数の導関数の反対称特性により、相互作用する粒子の各ペアの運動量保存が保証されます。 ${\displaystyle i}$${\displaystyle a}$

粒状材料のシミュレーションに使用される離散要素法は、SPH に関連しています。

• Hoover, WG (2006); 滑らかな粒子応用力学: 最先端の技術、World Scientific。
• ステリングワーフ、RF;カリフォルニア州ウィンゲート。 「SPH によるインパクト モデリング」、Memorie della Societa Astronomia Italiana、Vol. 65、p. 1117年（1994年）。
• アマダ、T.井村正人；安室裕子;真鍋裕也;および千原和也 (2004)。 「GPU での粒子ベースの流体シミュレーション」、グラフィックス プロセッサ上の汎用コンピューティングに関する ACM ワークショップ議事録(2004 年 8 月、カリフォルニア州ロサンゼルス)。
• Desbrun, M.; および Cani, M.-P. (1996)。「Smoothed Particles: 高度に変形可能な物体をアニメーション化する新しいパラダイム」、Eurographics Workshop on Computer Animation and Simulation (1996 年 8 月、フランス、ポワティエ) の議事録。
• Hegeman, K.、Carr, NA、および Miller, GSP、「GPU での粒子ベースの流体シミュレーション」、Proceedings of International Conference on Computational Science (Reading、英国、2006 年 5 月)、Lecture Notes in Computer Science v. 3994/2006 (Springer-Verlag)。
• Kelager, M. (2006) 「平滑化粒子流体力学を使用したラグランジュ流体力学」、修士論文、コペンハーゲン大学。
• Kolb, A.、および Cuntz, N. (2005)、「GPU ベースの流体シミュレーションのための動的粒子結合」、第 18 回シミュレーション技術シンポジウムの議事録(2005) 722 ～ 727 ページ。
• Liu, GR; および Liu, MB; Smoothed Particle Hydrodynamics: a meshfree particles method、シンガポール: World Scientific (2003)。
• モナハン、ジョセフ・J. (1992)。「平滑化粒子流体力学」、天文学と天体物理学年次レビュー(1992)。30: 543–74。
• Muller, M.、Charypar, D.、および Gross, M.、「インタラクティブ アプリケーション向けの粒子ベースの流体シミュレーション」、Breen, D.、および Lin, M. (編)、Proceedings of Eurographics/SIGGRAPH Symposium on Computer Animation (2003)。
• Vesterlund, M.;平滑化粒子流体力学を使用した粘性流体のシミュレーションとレンダリング、修士論文、ウメオ大学、スウェーデン。
• Violeau, D.;流体力学とSPH法、オックスフォード大学出版局 (2012)。

• SPH を使用した初の大規模星形成シミュレーション
• SPHERIC (SPH 研究・エンジニアリング国際コミュニティ)
• ITVO は、数値シミュレーション アーカイブのデータベースを照会するために作成されたイタリア理論仮想観測所の Web サイトです。
• SPHC 画像ギャラリーには、SPH コード SPHC のさまざまなテスト ケース、実験検証、商用アプリケーションが示されています。
• ナビエ・ストークス方程式から始まるSPHモデルの導出

• AlgodooはSPHを使用した教育用の2Dシミュレーションフレームワークです。
• AQUAgpusphは研究者による研究者のための無料（GPLv3）SPHです。
• dive solutionsは、CFD目的の商用WebベースのSPHエンジニアリングソフトウェアです。
• DualSPHysics は、SPHysics をベースに GPU コンピューティングを使用する、ほぼオープンソースの SPH コードです。オープンソース コンポーネントは LGPL の下で利用できます。
• FLUIDS v.1 は、CPU および GPU 用の液体用の、C++ によるシンプルなオープン ソース (Zlib) リアルタイム 3D SPH 実装です。
• FluidixはOneZero Softwareから入手可能なGPUベースの粒子シミュレーションAPIです。
• GADGET [1]は宇宙論的N体/SPHシミュレーションのための無料で利用できる（ GPL ）コードである。
• GPUSPH 粘度付き SPH シミュレータ (GPLv3)
• Pasimodoは、SPHなどの粒子ベースのシミュレーション手法のためのプログラムパッケージです。
• LAMMPSは、SPHシミュレーションを実行できる大規模並列オープンソースの古典分子動力学コードです。
• Physics Abstraction Layerは、SPHサポートを備えたリアルタイム物理エンジンをサポートするオープンソースの抽象化システムです。
• PreonLabは、FIFTY2 Technologyが開発した暗黙のSPH法を実装した商用エンジニアリングソフトウェアです。
• Puntoは粒子シミュレーションのための無料で利用できる視覚化ツールです。
• pysph Python での平滑化粒子流体力学のためのオープンソースフレームワーク (新しい BSD ライセンス)
• Py-SPHViewer 平滑化粒子流体力学シミュレーション用のオープンソースPython可視化ツール。^[1]
• 映画業界向けの RealFlow Commercial SPH ソルバー。
• RheoCubeは、複雑な流体のレオロジーと安定性の研究と予測を目的としたLorenz Researchの商用SaaS製品です。
• SimPARTIXは、Fraunhofer IWMのSPHおよび離散要素法（DEM）シミュレーション用の商用シミュレーションパッケージです。
• SPHフロー
• スフェラ
• SPHinXsys は、オープンソースのマルチフィジックス、マルチ解像度の SPH ライブラリです。物理的に正確なシミュレーションのための C++ API を提供し、流体、固体、マルチボディダイナミクスなどを含む結合された産業用動的システムをモデル化することを目的としています。
• SPHysicsはFortranで書かれたオープンソースのSPH実装です。
• SPLASHはSPHシミュレーション用のオープンソース（GPL）可視化ツールです。
• SYMPLER: フライブルク大学のフリーウェアの SYMbolic ParticLE シミュレーター。
• Nauticle は、粒子ベースの数値法のための汎用計算ツールです。
• NDYNAMICS は、CENTROID LAB が開発した暗黙的 SPH に基づく商用流体シミュレーション ソフトウェアで、現在、内部/外部浸水/原子力/化学工学アプリケーションに使用されています。