粒子有限要素法（PFEM : Particle Finite Element Method）の紹介

Posted at 2026-03-31

概要

　土木の分野において、粒子法の研究が盛んに行われております。特に土石流や津波、液状化現象などのシミュレーションで活用されています。

　今回の記事では、粒子有限要素法（PFEM: Particle Finite Element Method）について紹介します。PFEMは弱形式で定式化されているため自由境界条件が自然に導入され、リメッシングによってメッシュねじれの問題に対応できる手法です。

　本記事は以下の論文の内容を参考にしています。

　Matlabですが、SPFEMのプログラムが公開されています。

支配方程式と弱形式

　支配方程式は以下で表される。

\begin{align*}
\rho \boldsymbol{a}=\nabla \cdot \boldsymbol{\sigma} + {\bf b} 
\end{align*}

$\boldsymbol{\sigma}$は応力テンソル、${\bf a}$は加速度（変位の2階時間微分）、${\bf b}$は物体力である。${\bf u}$を変位として、ノイマン・ディリクレ境界条件は以下で与えられる。

\begin{align*}
 \begin{cases}
  \boldsymbol{\sigma} \cdot {\bf n} = \bar {\bf t}  & \mbox{on } \Gamma_t \\ 
    {\bf u} = \bar {\bf u}  & \mbox{on } \Gamma_u 
\end{cases}
\end{align*}

応力$\boldsymbol{\sigma}$とひずみ$\boldsymbol{\varepsilon}$の関係は、線形弾性体の構成則により

\begin{align*}
\sigma_{ij} = D_{ijkl} \varepsilon_{kl}
\end{align*}

応力テンソルと外向きの単位ベクトル${\bf n}$を用いて、表面力${\bf t}$は以下で表せる。

\begin{align*}
t_i(u) = \sigma_{ij}(u)n_j
\end{align*}

支配方程式に試験関数${\bf v}$を乗じて積分を行うと、以下の弱形式が得られる。

\begin{align*}
\int_{\Omega} \mathrm{d}\Omega \boldsymbol{v}^T \boldsymbol{a} =\int_{\Gamma} \mathrm{d}\Gamma \boldsymbol{v}^T \boldsymbol{t} - \int_{\Omega} \mathrm{d}\Omega \boldsymbol{\varepsilon}^T(v)\boldsymbol{\sigma}(u) 
    + \int_{\Omega} \mathrm{d}\Omega \boldsymbol{v}^T \boldsymbol{b}  \end{align*}

Smoothed Particle Finite Element Method

　FEMでは、ひずみ・応力は要素内の積分点で評価される。しかし、PFEMはリメッシュを行うため、要素の情報が初期化され、ひずみなどの履歴変数の情報が失われてしまう。そのため、SPFEMでは平滑化手法を導入し、積分点ではなく粒子でひずみを評価することで、リメッシング後も情報を保持する。

三角形メッシュにおける形状関数とその微分

　節点$1, 2, 3$の座標を$(x_1, y_1)$, $(x_2, y_2)$, $(x_3, y_3)$とすると、要素面積$A_e$は外積を用いて以下のように計算される。

\begin{align*}
A_e = \frac{1}{2} \left| (x_2 - x_1)(y_3 - y_1) - (x_3 - x_1)(y_2 - y_1) \right|
\end{align*}

　三角形要素における線形形状関数$N_i$ ($i=1,2,3$)は以下で定義され

\begin{align*}
N_i(x,y) = \frac{1}{2A_e}(a_i + b_i x + c_i y)
\end{align*}

各節点における係数$a_i$, $b_i$, $c_i$は以下のように定義される。

\begin{align*}
a_1 &= x_2 y_3 - x_3 y_2  \ \ , \ \ 
b_1 = y_2 - y_3 \ \ , \ \ 
c_1 = x_3 - x_2　\\
a_2 &= x_3 y_1 - x_1 y_3 \ \ , \ \ 
b_2 = y_3 - y_1 \ \ , \ \ 
c_2 = x_1 - x_3　\\
a_3 &= x_1 y_2 - x_2 y_1 \ \ , \ \ 
b_3 = y_1 - y_2 \ \ , \ \ 
c_3 = x_2 - x_1
\end{align*}

形状関数は以下の性質（Partition of Unity）を満たす。

\begin{align*}
\sum_{i=1}^{3} N_i = 1
\end{align*}

形状関数の空間微分は以下となり、

\begin{align*}
\frac{\partial N_i}{\partial x} = \frac{b_i}{2A_e}, \quad \frac{\partial N_i}{\partial y} = \frac{c_i}{2A_e}
\end{align*}

FEMで使われるひずみ-変位関係行列$\mathbf{B}$は以下で表現される。

\begin{align*}
\mathbf{B} = \begin{bmatrix}
\frac{\partial N_1}{\partial x} & 0 & \frac{\partial N_2}{\partial x} & 0 & \frac{\partial N_3}{\partial x} & 0 \\
0 & \frac{\partial N_1}{\partial y} & 0 & \frac{\partial N_2}{\partial y} & 0 & \frac{\partial N_3}{\partial y} \\
\frac{\partial N_1}{\partial y} & \frac{\partial N_1}{\partial x} & \frac{\partial N_2}{\partial y} & \frac{\partial N_2}{\partial x} & \frac{\partial N_3}{\partial y} & \frac{\partial N_3}{\partial x}
\end{bmatrix}
\end{align*}

平滑処理

　粒子$k$における平滑処理を実行したひずみ$\bar{\varepsilon}_k$は、以下で表現する。

\begin{align*}
\bar{\boldsymbol{\varepsilon}}_k = \int_{\Omega_k^s} \mathrm{d}\Omega \boldsymbol{\varepsilon}(x) \Phi_k(x)
\end{align*}

$\Phi_k$はPartition of Unityを満たし

\begin{align*}
\int_{\Omega_k^s} \mathrm{d}\Omega \Phi_k(x) =1
\end{align*}

具体的には、以下で定義される。

\begin{align*}
\Phi_k(x) =  \begin{cases}
 1/A_k^s  & x \in \Omega_k^s \\ 
 0 & x \notin \Omega_k^s
\end{cases}
\end{align*}

$\Omega_k^s$は以下の図のように、粒子$k$に隣接する要素の集合を意味する。

$A_k^s$はPartition of Unityを満たすように以下で定義される。

\begin{align*}
A_k^s = \int_{\Omega_k^s} \mathrm{d}\Omega =\frac{1}{3}\sum_{j=1}^{N_e^k} A_j^e
\end{align*}

また、重心と辺の中点を結んで得られる四角形の面積が、三角形の面積の$1/3$であることを利用している。

　最終的に粒子$k$におけるひずみ$\bar{\varepsilon}_k$は、以下となる。

\begin{align*}
\bar{\boldsymbol{\varepsilon}}_k = \frac{1}{A_k^s} \sum_{j=1}^{N_k^e} \frac{A_j^e}{3} \bar{\boldsymbol{\varepsilon}}_j
\end{align*}

SPFEMで使われる平滑化された行列$\bar{\mathbf{B}}$は、FEMで使われる$\mathbf{B}$を用いて以下のように表現できる。

\begin{align*}
\bar{{\bf B}}_k = \frac{1}{A_k^s} \sum_{j=1}^{N_k^e} \frac{A_j^e}{3}{\bf B}_j 
\end{align*}

　従って、弱形式のひずみエネルギー部分の項を離散化すると以下で表現できる。

\begin{align*}
\int_{\Omega} \mathrm{d}\Omega \mathbf{B}:\mathbf{\sigma} = \sum_{k=1}^{N_n} \sum_{i \in \text{adj}(k)} \bar{{\bf B}}_i^T \boldsymbol{\sigma}_k A_k^s
\end{align*}

$N_n$は全粒子数、$\text{adj}(k)$は粒子$k$に隣接する粒子の集合である。有限要素方程式の離散化形式は、以下で表現できる。

\begin{align*}
\sum_{k=1}^{N_n} \rho A_k^s \mathbf{a}_k = -\sum_{k=1}^{N_n} \sum_{i \in \text{adj}(k)} \bar{{\bf B}}_i^T \boldsymbol{\sigma}_k A_k^s + \int_{\Gamma} \mathrm{d}\Gamma \boldsymbol{N}^T \boldsymbol{t}
\end{align*}

　速度と変位は陽的時間積分により以下のように更新される。

\begin{align*}
\mathbf{v}_k &\leftarrow \mathbf{v}_k + (\mathbf{a}_k + \mathbf{g})\Delta t  \\ 
\mathbf{u}_k &\leftarrow \mathbf{u}_k + \mathbf{v}_k \Delta t
\end{align*}

$\mathbf{g}$は重力加速度である。物体力は弱形式の計算では考慮せず、速度更新時に直接加算する。

Alpha Shape アルゴリズム

　PFEMでは、メッシュの品質が低下したら、再度メッシュ分割（リメッシング）を実行する。しかし、凸包の性質により、ドローネ三角分割だけでは計算に不要な三角形も含まれてしまう。そこで、PFEMではAlpha Shapeアルゴリズムを用いて、余分な三角形を取り除く。

　Alpha Shapeは、半径$\alpha h$の球の中に入らない三角形を削除するアルゴリズムである。$\alpha$は設定するパラメータであり、$h$は平均的な要素サイズであり、以下で計算される。

\begin{align*}
h = \frac{1}{N_e} \sum_{e=1}^{N_e} \min\{a_e,b_e,c_e\}
\end{align*}

$a_e,b_e,c_e$は、要素$e$の三角形の各辺の長さである。また、$N_e$は全要素数である。各要素$e$の外接円の半径$R_e$は以下で計算される。

\begin{align*}
R_e = \frac{a_eb_ec_e}{4A_e}
\end{align*}

以下の図のように、外接円の半径$R_e$が閾値$\alpha h$より小さい場合は要素を残し、大きい場合は削除する。

　下の図は、Alpha Shape アルゴリズムを適用前と適用後の結果である。余分な三角形が削除されていることが確認できる。

リメッシングを行わない場合のシミュレーション結果を以下に示す。メッシュねじれが発生し、計算が破綻していることが確認できる。

メッシュ品質評価

　PFEMでは、メッシュの品質が低下すると計算精度が悪化したり、数値的不安定性が生じる。そのため、定期的にメッシュ品質を評価し、必要に応じてリメッシングを実行する。

　面積$A$と辺$b,c$に挟まれる角度$\theta$ との関係は、以下で表現できる。

\begin{align*}
A = \frac{1}{2}b c \sin \theta
\end{align*}

従って、辺$b,c$に挟まれる角は、

\begin{align*}
\theta_a = \arcsin \left(\frac{2A_e}{b c}\right) 
\end{align*}

辺$a,c$に挟まれる角は、

\begin{align*}
\theta_b = \arcsin \left(\frac{2A_e}{a c}\right) 
\end{align*}

辺$a,b$に挟まれる角は、

\begin{align*}
\theta_c = \arcsin \left(\frac{2A_e}{a b}\right) 
\end{align*}

と表せる。全要素について最小内角$\theta_{\min}$

\begin{align*}
\theta_{\min} = \min_{e1=,\dots,N_e}\{\min(\theta_a,\theta_b,\theta_c)  \}
\end{align*}

を計算し、閾値$\theta_{threshold}$を下回った場合

\begin{align*}
\theta_{\min}<\theta_{threshold}
\end{align*}

にリメッシングを実行する。また、メッシュの品質をさらに向上させるため、境界でない粒子に対してラプラシアンスムージングを適用する。

\begin{align*}
\mathbf{x}_i^{\text{new}} &= 
\begin{cases}
\displaystyle\frac{1}{N_{\text{adj}}} \sum_{j \in \text{adj}(i)} \mathbf{x}_j & \text{if } d_{\text{min}}^i \leq \lambda h \\
\mathbf{x}_i & \text{otherwise}
\end{cases}
\\
d_{\text{min}}^i &= \min_{j \in \text{adj}(i)} \|\mathbf{x}_i - \mathbf{x}_j\|  \ \ , \ \ \lambda=0.2
\end{align*}

Drucker-Prager 非粘着性土壌斜面の崩壊

　構成則にDrucker-Prager モデルを採用し、非粘着性土壌斜面の崩壊シミュレーションを実施した。材料定数および数値パラメータは以下の表に示す。

項目	記号	値	単位
密度	$\rho$	$2142.0$	$\mbox{kg/}m^3$
ヤング率	$E$	$70.0$	$\mbox{MPa}$
ポアソン比	$\nu$	$0.3$	-
粘着力	$C$	$0.0$	$\mbox{Pa}$
内部摩擦角	$\phi$	$20.0,30.0,40.0$	-
ダイレイタンシー角	$\psi$	$0.0$	-
時間刻み幅	$\Delta t$	$10^{-4}$	$\mbox{s}$
初期重力計算時間	$\Delta t$	$15$	$\mbox{s}$
シミュレーション時間	$\Delta t$	$12$	$\mbox{s}$
粒子間隔	$h$	$0.5$	$\mbox{m}$

　1次の三角形要素を使用するため、アワーグラスモードが生じる。そのため、アワーグラス制御を導入する。
　
　解析モデルは、高さ$15\mbox{m}$の土台に、斜面角度$45°$、高さ$20\mbox{m}$の斜面を配置したモデルである。境界条件として、底面は完全固定し、左右端は$x$方向の速度をゼロとした。