Fragile Points Method（FPM）の基礎理論

概要

　亀裂進展解析では、XFEM（拡張有限要素法）やペリダイナミックス、フェーズフィールド法などが挙げられます。近年、定式化されたFragile Point Method（FPM）は、メッシュレス解析手法として提案された手法です。FPMは、点間の結合を切断するだけで亀裂を表現できるため、リメッシュや形状関数の拡張が不要という大きな利点があります。また、FPMの方程式は弱形式で定式化されているため、自由境界条件の取り扱いが簡単になります。

　計算工学会でFPMの発表があり興味を持ちました。

　今回の記事は、亀裂進展解析までは実装できていないため、FPMの基礎的な理論について紹介したいと思います。

　以下の文献を参考にしています。

支配方程式

　線形弾性体の支配方程式は以下で表される。

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

$\boldsymbol{\sigma}$は応力テンソル、${\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*}

FPMの定式化

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

\begin{align*}
\int_{\Gamma} \mathrm{d}\Gamma v_i t_i(u) - \int_{\Omega} \mathrm{d}\Omega \varepsilon_{ij}(v)\sigma_{ij}(u) 
    + \int_{\Omega} \mathrm{d}\Omega v_i b_i  = 0
\end{align*}

第2項に注目し、部分積分を2回実行すると

\begin{align*}
\int_{\Omega_e} \mathrm{d}\Omega \varepsilon_{ij}(v)\sigma_{ij}(u) &= \int_{\Omega_e} \mathrm{d}\Omega \varepsilon_{ij}(v) D_{ijkl}\varepsilon_{kl}(u) \\
&= \frac{1}{2}\int_{\Gamma_e} \mathrm{d}\Gamma \varepsilon_{ij}(v) D_{ijkl} (u_k n_l + u_l n_k) - \frac{1}{2}\int_{\Omega_e} \mathrm{d}\Omega  \left(u_k \frac{\partial}{\partial x_l} + u_l \frac{\partial}{\partial x_k} \right)D_{ijkl}\varepsilon_{ij}(v) \\
&= \int_{\Gamma_e} \mathrm{d}\Gamma (u_k t_k^*(v) - u_k t_k(v)) + \int_{\Omega_e} \mathrm{d}\Omega \varepsilon_{ij}(v) D_{ijkl}\varepsilon_{kl}(u)
\end{align*}

と書ける。2行目の第1項を${\bf t^*}$（数値フラックス項）として定義し、3行目で再度部分積分を実行し生まれた項を${\bf t}$とした。また、$\Omega_e,\Gamma_e$ の下付き添字は要素$e$の領域を意味する。

　ジャンプ演算子と平均演算子を以下の2重括弧で定義する。

\begin{align*}
[\![ a ]\!] = a^{(e)} - a^{(k)} \ \ , \ \ \{\![ a ]\!\} = \frac{1}{2}(a^{(e)} + a^{(k)})
\end{align*}

　以下の図に示すように、

部分領域 $e$に隣接する部分領域$k$との境界における積分は、

\begin{align*}
\int_{\Gamma_e} \mathrm{d}\Gamma v_i \sigma_{ij}n_j =& \int_{e} \mathrm{d}\Gamma v_i^{(e)} \sigma_{ij}^{(e)}n_j^{(e)} + \int_{e} \mathrm{d}\Gamma v_i^{(k)} \sigma_{ij}^{(k)}n_j^{(k)}  \\
=& \int_{e} \mathrm{d}\Gamma v_i^{(e)} \sigma_{ij}^{(e)}n_j^{(e)} - \int_{e} \mathrm{d}\Gamma v_i^{(k)} \sigma_{ij}^{(k)}n_j^{(e)} \\
=&\int_{e} \mathrm{d}\Gamma [\![v_i]\!] \{\![\sigma_{ij}n_j]\!\} + \int_{e} \mathrm{d}\Gamma \{\![v_i ]\!\}  [\![ \sigma_{ij}n_j]\!]
\end{align*}

と計算できる。2行目で${\bf n}^{(e)} = -{\bf n}^{(k)}$を用いた。3行目の2項目は、要素間の流れる量は等しくなるのでゼロになる。

　数値フラックス項は、以下のように定義される。

\begin{align*}
t_i^*(v) = \sigma_{ij}^*n_j = \frac{\eta}{h_e} n_i [\![ v_j]\!] n_j = \frac{\eta}{h_e} [\![ v_i]\!]
\end{align*}

$h_e$は境界の長さ、$\eta$はペナルティーパラメータで、今回の記事ではヤング率と同じ値にする。

　全体の境界$\Gamma$からノイマン・ディリクレ境界条件で与えられる境界$\Gamma_t,\Gamma_u$ を除いた内部境界を境界を以下で定義する。

\begin{align*}
\Gamma_h = \Gamma -\Gamma_u - \Gamma_t
\end{align*}

最終的に、FPMの方程式は以下となる。

\begin{align*}
\sum_e \int_{\Omega_e} \mathrm{d}\Omega \boldsymbol{\varepsilon}(v)^T\boldsymbol{\sigma}(u) &-\sum_{e\in\Gamma_h} \int_{e} \mathrm{d}\Gamma  [\![{\bf v}]\!]^T \{\![{\bf n}_e\boldsymbol{\sigma}(u)]\!\}-\sum_{e\in\Gamma_h} \int_{e} \mathrm{d}\Gamma \{\![{\bf n}_e\boldsymbol{\sigma}(v)]\!\}^T [\![{\bf u}]\!]   \\
&+\sum_{e\in\Gamma_h} \frac{\eta}{h_e}\int_{e} \mathrm{d}\Gamma  [\![{\bf v}]\!] [\![{\bf u}]\!] -\sum_{e\in\Gamma_t} \int_{e} \mathrm{d}\Gamma {\bf v}^T \bar{{\bf t}}(u)- \sum_e \int_{\Omega_e} \mathrm{d}\Omega {\bf v}^T {\bf b}  = {\bf 0}
\end{align*}

Generalized Finite Difference

　2次元問題を考える。部分領域を代表する粒子の座標を$(x^0,y^0)$とし、変位を$(u_x^0,u_y^0)$とする。部分領域内の任意の位置$(x,y)$の変位を、計算点$(x^0,y^0)$まわりでテイラー展開により以下のように近似する。

\begin{align*}
{\bf u}^h = 
\begin{pmatrix}
u_x^0 \\
u_y^0 \\
\end{pmatrix}
+ 
\begin{pmatrix}
{\bf h} & 0 \\
0 & {\bf h} \\
\end{pmatrix}
\begin{pmatrix}
{\bf a}_x \\
{\bf a}_y \\
\end{pmatrix}
\end{align*}

${\bf h}$と${\bf a}_x,{\bf a}_y$は以下で表せる。

• 1次のテイラー展開の場合

\begin{align*}
{\bf h}(x,y) = 
\begin{pmatrix}
x-x^0 & y-y^0) 
\end{pmatrix}
\end{align*}

\begin{align*}
{\bf a}_x &= 
\begin{pmatrix}
\frac{\partial u_x}{\partial x} & \frac{\partial u_x}{\partial y} 
\end{pmatrix}\bigg|_{(x^0,y^0)}
\\
{\bf a}_y &= 
\begin{pmatrix}
\frac{\partial u_y}{\partial x} & \frac{\partial u_y}{\partial y} 
\end{pmatrix}\bigg|_{(x^0,y^0)}
\end{align*}

• 2次のテイラー展開の場合

\begin{align*}
{\bf h}(x,y) = 
\begin{pmatrix}
x-x^0 & y-y^0 & \frac{1}{2}(x-x^0)^2 & \frac{1}{2}(y-y^0)^2 & (x-x^0)(y-y^0) 
\end{pmatrix}
\end{align*}

\begin{align*}
{\bf a}_x &= 
\begin{pmatrix}
\frac{\partial u_x}{\partial x} & \frac{\partial u_x}{\partial y} & \frac{\partial^2 u_x}{\partial x^2} & \frac{\partial^2 u_x}{\partial y^2} & \frac{\partial^2 u_x}{\partial x\partial y}
\end{pmatrix}\bigg|_{(x^0,y^0)}
\\
{\bf a}_y &= 
\begin{pmatrix}
\frac{\partial u_y}{\partial x} & \frac{\partial u_y}{\partial y} & \frac{\partial^2 u_y}{\partial x^2} & \frac{\partial^2 u_y}{\partial y^2} & \frac{\partial^2 u_y}{\partial x\partial y}
\end{pmatrix}\bigg|_{(x^0,y^0)}
\end{align*}

${\bf h}$は既知であるが、${\bf a}_x,{\bf a}_y$は未知量である。これらは最小二乗法により求める。以下、要素に隣接する$m$個の要素の代表粒子を用いて、行列とベクトルを定義する。

\begin{align*}
{\bf A}^T = 
\begin{pmatrix}
{\bf h}(x_1) & {\bf h}(x_2) & \cdots & {\bf h}(x_m)
\end{pmatrix}
\end{align*}

\begin{align*}
{\bf u}_{xm} =
\begin{pmatrix}
u_x^1 - u_x^0  \\
u_x^2 - u_x^0  \\
\ldots \\
u_x^m - u_x^0 
\end{pmatrix}
\ \ , \ \  
{\bf u}_{ym} =
\begin{pmatrix}
u_y^1 - u_y^0  \\
u_y^2 - u_y^0  \\
\ldots \\
u_y^m - u_y^0 
\end{pmatrix}
\end{align*}

目的関数$J_x,J_y$を以下で定義する。

\begin{align*}
J_x =& ({\bf A} {\bf a}_x  - {\bf u}_{xm})^T ({\bf A} {\bf a}_x  - {\bf u}_{xm}) \\
J_y =& ({\bf A} {\bf a}_y  - {\bf u}_{ym})^T ({\bf A} {\bf a}_y  - {\bf u}_{ym})
\end{align*}

従って、

\begin{align*}
\frac{\partial J_x}{\partial {\bf a}_x} = 0  \ \ , \ \ 
\frac{\partial J_y}{\partial {\bf a}_y} = 0  
\end{align*}

より

\begin{align*}
{\bf a}_x = {\bf C} {\bf u}_{xm} \ \ , \ \  
{\bf a}_y = {\bf C} {\bf u}_{ym} 
\end{align*}

が得られる。また${\bf C}$は以下である。

\begin{align*}
{\bf C} = ({\bf A}^T {\bf A})^{-1} {\bf A}^T
\end{align*}

　最終的に変位場の近似式は以下で表される。

\begin{align*}
\mathbf{u}^h(x, y) = \mathbf{N}(x, y) \mathbf{u}^E
\end{align*}

ここで、$\mathbf{N}$は形状関数（$2 \times 2(m+1)$ 行列）であり、$\mathbf{u}^E$は変位ベクトルで以下のように定義される。

\begin{align*}
\mathbf{u}^E = \begin{pmatrix}
u_{x,1} & u_{y,1} & u_{x,2} & u_{y,2} & \cdots & u_{x,m} & u_{y,m} & u_{x,0} & u_{y,0}
\end{pmatrix}^T
\end{align*}

ひずみ場の近似式は以下で表される。

\begin{align*}
\boldsymbol{\varepsilon}^h(x, y) = \mathbf{B}(x, y) \mathbf{u}^E
\end{align*}

\begin{align*}
\boldsymbol{\varepsilon}^h = \begin{pmatrix}
\varepsilon_{xx} \\
\varepsilon_{yy} \\
\gamma_{xy}
\end{pmatrix}
= \begin{pmatrix}
\frac{\partial u_x}{\partial x} \\
\frac{\partial u_y}{\partial y} \\
\frac{\partial u_x}{\partial y} + \frac{\partial u_y}{\partial x}
\end{pmatrix}
\end{align*}

形状関数$\mathbf{N}$とその微分の表現$\mathbf{B}$は、$\mathbf{C}$で表せる。

剛性行列

　FPMの弱形式を離散化すると、以下の連立一次方程式が得られる。

\begin{align*}
\mathbf{K}\mathbf{u} = \mathbf{Q}
\end{align*}

　ここで

• $\mathbf{K}$：全体剛性行列
• $\mathbf{u}$：変位ベクトル
• $\mathbf{Q}$：荷重ベクトル

剛性行列はバルク部分と境界部分に分けられる。

• バルク部分

\begin{align*}
{\bf K}_E = \int_{\Omega} \mathrm{d}\Omega {\bf B}^T {\bf D} {\bf B} 
\end{align*}

$\mathbf{D}$は弾性テンソルである。

• 境界部分

\begin{align*}
{\bf K}_h =&-\frac{1}{2} \int_{e} \mathrm{d}\Gamma {\bf N}_1^T{\bf n}_e  {\bf D} {\bf B}_1 
-\frac{1}{2} \int_{e} \mathrm{d}\Gamma {\bf B}_1^T{\bf D} {\bf n}_e^T  {\bf N}_1 
+\frac{\eta}{h_e} \int_{e} \mathrm{d}\Gamma {\bf N}_1^T {\bf N}_1 \\
&-\frac{1}{2} \int_{e} \mathrm{d}\Gamma {\bf N}_1^T{\bf n}_e  {\bf D} {\bf B}_2 
+\frac{1}{2} \int_{e} \mathrm{d}\Gamma {\bf B}_1^T{\bf D} {\bf n}_e^T  {\bf N}_2 
-\frac{\eta}{h_e} \int_{e} \mathrm{d}\Gamma {\bf N}_1^T {\bf N}_2 \\
&+\frac{1}{2} \int_{e} \mathrm{d}\Gamma {\bf N}_2^T{\bf n}_e  {\bf D} {\bf B}_1 
-\frac{1}{2} \int_{e} \mathrm{d}\Gamma {\bf B}_2^T{\bf D} {\bf n}_e^T  {\bf N}_1 
-\frac{\eta}{h_e} \int_{e} \mathrm{d}\Gamma {\bf N}_2^T {\bf N}_1 \\
&+\frac{1}{2} \int_{e} \mathrm{d}\Gamma {\bf N}_2^T{\bf n}_e  {\bf D} {\bf B}_2 
+\frac{1}{2} \int_{e} \mathrm{d}\Gamma {\bf B}_2^T{\bf D} {\bf n}_e^T  {\bf N}_2 
+\frac{\eta}{h_e} \int_{e} \mathrm{d}\Gamma {\bf N}_2^T {\bf N}_2 
\end{align*}

ここで、添字1と2は境界を共有する2つの部分領域を表す。また、荷重ベクトルは

\begin{align*}
{\bf Q} = \sum_{e\in\Gamma_t} \int_{e} \mathrm{d}\Gamma {\bf N}^T \bar{{\bf t}}+ \sum_e \int_{\Omega_e} \mathrm{d}\Omega {\bf N}^T {\bf b} 
\end{align*}

である。

数値積分

　${\bf K}_E,{\bf K}_h,{\bf Q}$を求めるには、数値積分を実行する必要がある。${\bf K}_h$における境界積分は、1点のガウス積分を実行した。

　境界$e$の数値積分は、$h_e$を境界の長さとして以下で計算する。

\begin{align*}
\int_{e} f(s) \, ds = \int_{-1}^{1} f(\xi) \, J(\xi) \, d\xi \approx f(\xi＝0) h_e 
 \end{align*}

$f(\xi=0)$は境界の中点の値に対応する。

　${\bf K}_E,{\bf Q}$の積分は、要素を三角形に分割して3点積分を実行した。三角形内の任意の座標$(x,y)$は面積座標$L_1,L_2,L_3$を用いて以下で表せる。

\begin{align*}
x &= L_1 x_1 + L_2 x_2 + L_3 x_3 \\
y &= L_1 y_1 + L_2 y_2 + L_3 y_3
\end{align*}

ここで

\begin{align*}
L_3 = 1-L_1-L_2
\end{align*}

関数$f(L_1,L_2)$の三角形上の数値積分は、積分点$(L_{i1},L_{i2})$と重み $w_i$ を用いて

\begin{align*}
\int_0^1 \mathrm{d}L_1 \int_0^{1-L_1} \mathrm{d}L_2 \mathrm{det} |J| f(L_1,L_2) \approx \sum_{i=1}^3 \mathrm{det} |J| w_if(L_{i1},L_{i2})
\end{align*}

と表せる。ヤコビアン$\mathrm{det}|J|$は三角形の面積の2倍である。

\begin{align*}
\mathrm{det} |J| = 2 \triangle
\end{align*}

　3点の積分点と重みの組み合わせは、以下である。

\begin{align*}
(L_{i1},L_{i2},w_i) = (1/2,1/2,1/6),\ \ (1/2,0,1/6), \ \ (0,1/2,1/6)
\end{align*}

簡単な検証

　
　プログラムの検証のため、正方形プレートに荷重をかけて理論解と一致するか検証した。材料定数および数値パラメータは以下の表に示す。

項目	記号	値	単位
ヤング率	$E$	210	$\mbox{GPa}$
ポアソン比	$\nu$	0.3	-
ペナルティパラメータ	$\eta$	$E$ (= 210 GPa)	$\mbox{GPa}$
分布荷重	-	1.0	$\mbox{GPa}$

　FPMの要素1つに対して、粒子1個が割り当てられる。要素はボロノイ図から生成しても良いが、簡単に検証するため四角形要素とした。四角形要素の場合、2次のテイラー展開における交差微分項$(x-x^0)(y-y^0)$が常にゼロとなってしまい、行列$\mathbf{C}$が特異になるので、1次のテイラー展開で計算を実施した。

正方形のプレートの引張

　以下の図は、青色の点は粒子を表し、バツ印はディリクレ境界条件（固定点）を表し、下部を完全固定（$u_x = u_y = 0$）、上部は$x$方向の変位を拘束（$u_x = 0$）した。オレンジ色の点はノイマン境界条件を表し上下方向に荷重をかける。

　このときの理論解は、$\sigma_{xx}=0.0,\sigma_{xy}=0.0,\sigma_{yy}=1.0\mbox{GPa}$である。以下の以下の図はFPMの計算結果である。理論解との誤差は最大で$1.0\times10^{-8}$程度であり、高精度に計算できていることが確認できる。また、ヒートマップは$1.0\mbox{GPa}$で正規化している。

正方形のプレートのせん断

　次にせん断方向に荷重をかけて計算を行った。

　このときの理論解は、$\sigma_{xx}=0.0,\sigma_{xy}=1.0\mbox{GPa},\sigma_{yy}=0.0$である。以下の図はFPMの計算結果である。理論解との誤差は最大で$1.0\times10^{-8}$程度であり、高精度に計算できていることが確認できる。また、ヒートマップは$1.0\mbox{GPa}$で正規化している。

　

まとめ

　本記事では、Fragile Points Method（FPM）の基礎理論と実装について解説した。FPMは、点間の結合を切断するだけで亀裂を表現できるため、リメッシュや形状関数の拡張が不要という大きな利点がある。土木分野におけるコンクリート構造物のひび割れ解析や岩盤の破壊シミュレーションへの応用が期待される。

　今後は、亀裂進展解析に向けたアルゴリズム作成を行う。ソースコードは githubの準備ができ次第公開する。