はじめに
最近、地盤解析手法の1つとしてMaterial Point Method（以下、MPMとする）に注目しています。
MPMは格子法と粒子法のハイブリッド手法で、映画「アナと雪の女王」などで注目されています。
今回は、MPMの概要と進化の歴史について解説したいと思います。
なお、今回は全体像を理解することが主目的なので、計算式はなるべく省略しています。
MPM開発までの歴史
1964年にHarlowらがParticle-in-Cell（PIC）法と呼ばれる流体解析手法を提案しました。
このPIC法は、空間を動く粒子（Lagrangian描像）と固定された格子（Euler描像）を組み合わせて、プラズマや流体の物理方程式を解くハイブリッド型の数値解析手法です。
PIC法では以下の手順を繰返して計算します。
- Deposit: 粒子が持つ電荷や速度の情報を、まわりの背景格子（ノード）へ割り当てて密度や電流を計算します
- Field Solve: 格子上でポアソン方程式やマクスウェル方程式を解き、空間全体の電場や磁場を求めます
- Gather: 格子で求めた電磁場の強さを、各粒子の位置に割り戻します
- Push: ローレンツ力やニュートンの運動方程式を使い、各粒子の速度と位置を更新します
PIC法は、個々の粒子の動き（運動論的効果）をダイレクトに追跡できるため、従来の流体近似では扱えない複雑な非平衡プラズマ現象を高精度に再現できる一方、「粒子（連続体）」と「固定格子（メッシュ）」の間で情報を相互に変換するという基本構造に起因する、いくつかの欠点（限界）が存在します。
例えば、PIC法では、毎ステップごとに粒子の「すべての情報（位置・速度など）」を一度格子（グリッド）に集約し、格子上で新しい速度を計算して、それをまた粒子に戻しますが、情報がまるで平均化されて薄まるようにボヤけてしまうという欠点があります。
この欠点により、水のシミュレーションをしているはずなのに、計算誤差のせいで「ハチミツや泥のようにドロドロした、異常に粘り気の強い動き」になってしまい、水しぶきや激しい渦といったリアルな流体現象が完全に消し飛んでしまいました。
このような欠点を回避するために、1986年にBrackbillとRuppelがFLIP法という手法を提案しました。
PIC法では格子上で計算された「最新の速度」をそのまま粒子に上書きするのに対し、FLIP法では格子上で計算された「前ステップからの速度の変化量（差分：$\Delta v$）」だけを粒子に足し合わせるという違いがあります。
この「差分だけを適用する」というアプローチにより、FLIP法はエネルギーをほぼ完璧に保存できるようになり前述の欠点を回避することにしました。
FLIP法によって「ハイブリッド手法でもエネルギーを損なわずに計算できる」ことが証明されました。
1990年代に入り、このFLIP法の優れたエネルギー保存特性に目をつけた研究者たちが、「流体だけでなく、固体（ゴム、金属、プラスチックなど）の応力やひずみの計算にもこのサイクルが使えるのではないか？」と考えだしました。
その結果、FLIP法のアルゴリズムに固体の「構成式（物質の硬さや変形のルール）」を組み込む形で1994年にSulskyらの研究チームがMPMを提案しました。
MPMの概要
MPMは前述のPIC法をベースにLagrangian法（粒子法）とEuler法（格子法）を組み合わせたハイブリッド手法となります。
物質点はLagrangianに移動し、運動方程式の計算は背景Euler格子上で行うというのがMPMの基本方針となります。
MPMの計算の流れは以下の通りです。
- Particle to Grid: 粒子の質量・運動量・必要な力学量を格子節点へ写す
- Grid velocity computation: 節点質量と運動量から節点速度を求める
- Grid momentum update: 内力・外力・境界条件から節点加速度・速度を更新する
- Velocity gradient / deformation update: 格子速度場から粒子位置の速度勾配を計算し、変形勾配テンソル$\boldsymbol{F}$、応力、内部変数を更新する
- Grid to Particle: 格子速度または速度増分を粒子へ戻す
- Particle advection: 粒子座標を更新する
MPMの最大の利点は、有限要素法（FEM）では計算が破綻してしまうような「大いなる変形」や「破壊・飛散」を、固定グリッドを併用することで安定して計算できる点です。
その一方で初期MPMでは以下のような問題点がありました。
- cell-crossing error: 粒子が背景格子セルを横切る際の数値振動・応力ノイズ
- boundary representation: 物体境界が背景Euler格子と一致しないため、境界位置やNeumann境界条件の扱いが難しい
- particle distribution/quadrature error: 粒子が積分点を兼ねるため、粒子配置が偏ると数値積分精度が低下する
- contact treatment: 接触はMPMの長所でもあるが、摩擦・分離・複数物体を厳密に扱うと実装は複雑になる
特にcell-crossing errorは大変形を扱う上で重要な問題であるため、このcell-crossing errorを解決するための種々の方法が提案されています。
MPMの進化の歴史
cell-crossing errorを改善するための手法として、Generalized Interpolation Material Point（GIMP法）が提案されました。
粒子を点ではなく、粒子ごとに一定の支配領域（Particle dmain）を持つ「塊」として定義します。
標準MPMでは粒子位置一点で形状関数を評価するため、セル境界を越えたときに勾配が急変してcell-crossing noiseが出ます。
GIMPでは粒子に有限な領域を持たせ、その領域上で形状関数を平均化します。
これにより、MPMに比べて数値振動が圧倒的に少なくなり、エネルギー保存性が向上しました。
ただし、GIMP法では激しい回転やせん断（斜めの変形）」が起きると破綻するという問題がありました。
GIMP法の欠点を克服するため、Sadeghiradらは、粒子の領域を四角形に固定せず、粒子の4つの「角（コーナーホバー）」の位置を物体の変形に合わせてリアルタイムに追従・移動させる手法を開発しました。
この手法をConvected Particle Domain Interpolation（CPDI法）と呼びます。
これにより、粒子は平行四辺形や任意の多角形にねじれることができるようになり、激しいせん断や引き裂き、回転を伴う破壊シミュレーションでもエラーが起きなくなりました。
また、粒子の速度を近傍のセルに重みを付けて分配するというMPMの手法（PIC法のやり方）では、粒子→格子→粒子を繰り返すたびに高周波成分が平滑化され、強い数値散逸が入ります。
FLIP法の手法では、「速度」ではなく「格子速度の変化量」を戻して散逸を減らしました。
ただし今度はノイズが残りやすいという問題点があります。
この問題を解決するために、Affine Particle-in-Cell（APIC法）が開発されました。
APIC法では速度に加えてアフィン（回転・剪断）行列を粒子に持たせるという手法です。
これによって、単なる並進速度だけでは失われていた局所速度勾配や回転の情報をParticle to Grid/Grid to Particle間で保持しやすくなり、PICほど散逸せず、FLIPほどノイジーでもない転送が可能になります。
従来の粒子法では、粒子と格子（グリッド）の間で速度情報のやり取りを行う際に「回転のエネルギー（角運動量）」が失われがちでした。
APIC法では粒子が通常の速度情報に加え、周囲の変形や回転（アフィン変換）を表す速度勾配行列を保持するため、渦や回転の動きが綺麗に保存されます。
また、FLIP法で発生しやすかった、粒子が不自然にバラバラと飛び散るノイズ（Jittering）が発生しないというメリットがあります。
さて、このAPIC法を数学的に解明していくと、粒子の速度の行列成分は粒子位置における速度ベクトルの勾配を移動最小二乗法 (Moving Least Squares, MLS) で求めたものに一致するということがわかりました。
MLS とは、Least Squares（最小二乗法）の結果を局所的に得るため、適当な重み関数を用いて各データ点に重みを付けて最小二乗法による（主に多項式の）フィッティングを行う手法です。
APICが導入されたときは、「こういう情報を持たせれば良い転送になる」というアルゴリズム的な側面が強かったです。
一方、粒子周辺の速度場を最小二乗近似すれば、自然にこのaffine velocityが出てくるという見方ができます。
この試みは非常にうまくいき、結果としてMLS-MPMはMPMよりもシンプルで、計算が高速になり、APICと整合した手法になりました。
さて、MPM（PIC法に基づく）の問題点に対し、別の着目点から改良を試みた研究もあります。
おさらいですが、MPMの主な問題点は「粒子が背景格子セルを横切る際の数値振動・応力ノイズ」、「物体境界が背景Euler格子と一致しないため、境界位置やNeumann境界条件の扱いが難しい」でありました。
また、MPMはUpdated Lagrangianという手法を用いることが一般的です。
これは、現時点（または一歩前のステップ）の変形した構成配向を基準として、次のステップの増分を計算する方法です。
つまり、ステップごとに基準座標系が「更新」されます。
この手法は大変形を追跡する手法として便利な手法ではありますが、逐次更新で誤差が蓄積しやすいという欠点を抱えています。
以上の問題を解決するための手法として、A generalized particle in cell method（GPIC法）が提案されました。
GPIC法はこれまでの改良手法とは違いより格子法（有限要素法）に寄っているといえます。
従来のMPMでは背景格子と粒子の組合せで計算していましたが、GPIC法では背景格子と有限要素メッシュの組合せで計算します。
これにより、従来のMPMと比較し、大変形への適用性が下がりますが、Total Laglangian（初期状態の形状を基準座標系とする手法）有限要素法を用いることで誤差の蓄積を減らし、境界条件の扱いを簡単にしました。
文献等ではFEM-MPMハイブリッド手法と紹介されていたりしますが、粒子という考え方から離れてしまっているため、本当にMPMなの？という疑問が出てきます。
おわりに
固体力学ではFEMは微小変形の精度が高いが大変形に弱い、MPMは大変形に対応できるが微小変形の精度がFEMと比較し低いとされています。
ならば、FEMで計算し変形量等に応じてMPMに切り替えればいいのでは？と考えることもできるのではないかと思います。
このあたりは、現在調査中であるため、興味のある方はご連絡ください。
参考文献
- F. H. Harlow: The Particle-in-Cell Method for Numerical Solution of Problems in Fluid Dynamics, Proceedings of Symposium in Applied Mathematics, Vol.15, 1963.
- C. K. Birdsall and A. B. Langdon: Plasma Physics via Computer Simulation, CRC Press, 1985.
- D. Sulsky, Z. Chen and H. L. Schreyer: A particle method for history-dependent materials, Computer Methods in Applied Mechanics and Engineering, Vol.118, 1994.
- S. G. Bardenhagen and E. M. Kober: The Generalized Interpolation Material Point Method, Computer Modeling in Engineering & Sciences, Vol.5, 2004.
- A. Sadeghirad, R. M. Brannon and J. Burghardt: A convected particle domain interpolation technique to extend applicability of the material point method for problems involving massive deformations, International Journal for Numerical Methods in Engineering, Vol.86, 2011.
- C. Jiang, C. Schroeder, J. Teran, A. Stomakhin and A. Selle: The affine particle-in-cell method, ACM Transactions on Graphics, Vol.34, 2015.
- Y. Hu, Y. Fang, Z. Ge, Z. Qu, Y. Zhu, A. Pradhana and C. Jiang: A Moving Least Squares Material Point Method with Displacement Discontinuity and Two-Way Rigid Body Coupling, ACM Transactions on Graphics, Vol.37, 2018.
- V. P. Nguyen, A. de Vaucorbeil, S. Sinaie and G.-R. Liu, A generalized particle in cell method for explicit solid dynamics, Computer Methods in Applied Mechanics and Engineering, Vol.371, 2020.
- 解説: Jelly: https://blog.oimo.io/2022/05/11/mls-mpm-elastic/#fn:2, 2026年9月27日閲覧