港湾空港技術研究所報告第 62 巻第 2 号（2023.6）

テイラー級数近似と重みつき最小二乗法を用いたMaterial Point Method の高度化
中村

圭太*・松村

聡*・水谷

要

崇亮**

旨

本稿は Computer Methods in Applied Mechanics and Engineering で出版された論文 “Taylor particle-

in-cell transfer and kernel correction for material point method” の日本語訳版である．
Material point method（MPM）は，工学分野において大きな変位や変形を伴う問題の解決に広く用いられてきた．MPM のオリジナルの定式化は，FLIP（fluid-implicit-particle）法を採用しており，PIC（particle-in-cell）法に比べ，計算時の散逸は少ないが安定性も低い．Affine PIC（APIC）法はアフィン速度を導入し，PIC 法における角運動量の散逸を克服しながら安定したシミュレーションを実現した．本論文では，APIC 法の一種である一次のテイラー級数近似に基づくアフィン速度と PIC 法を組み合わせた Taylor–PIC（TPIC）法を示す．TPIC 法は極めてシンプルで（粒子速度勾配のみ必要），散逸の少なさや安定性といったオリジナルの APIC 法の主要な利点を受け継いでいる．APIC 法とは異なり TPIC 法は角運動量を厳密には保存しないが，粒子と背面格子の間の輸送において速度勾配が保存される．この速度勾配は，その反対称成分に角運動の情報を含んでいるため，TPIC 法は APIC 法と同様に角度運動を適切に記述することができる．
また，MPM は境界条件が明示的に課される場合，すなわち背面格子の境界上の速度（または運動量）が設定される場合，境界付近で応力振動を引き起こすことがある．アフィン型の輸送を用いる場
合，この不安定性は悪化し，シミュレーションが容易に破綻する．そこで，境界付近の粒子にカーネル関数を補正するカーネル補正法を提案する．提案した TPIC 法とカーネル補正法は，線形弾性体，vonMises，Newtonian fluid，Drucker–Prager の材料モデルを用いて 5 種類のシミュレーションで検証した．シミュレーションでは，小さな変形でも応力振動による解析誤差が観察されたが，カーネル補正法を適用することでこれらの振動を除去することに成功し，結果は解析解と一致した．また，解析結
果より TPIC 法の精度とロバスト性が確認できた．
キーワード：Material point method，Particle-in-cell 法，Fluid-implicit-particle 法，Affine particle-in-cell法，カーネル補正法，境界条件***

地盤研究領域 基礎工研究グループ 主任研究官
地盤研究領域 基礎工研究グループ グループ長〒239-0826 横須賀市長瀬 3-1-1 国立研究開発法人海上・港湾・航空技術研究所 港湾空港技術研究所
電話：046-844-5057 Fax：046-844-0618 e-mail：nakamura-ke@p.mpat.go.jp-3-

REPORT OF THE PORT AND
AIRPORT RESEARCH INSTITUTE
Vol.62, No.2 (Jun. 2023)

Advanced Material Point Method Using
Taylor Series Approximation and Weighted Least Squares
Keita NAKAMURA*
Satoshi MATSUMURA*
Takaaki MIZUTANI**
Synopsis
This paper is a Japanese translation of the article “Taylor particle-in-cell transfer and kernel
correction for material point method” published in Computer Methods in Applied Mechanics.
The material point method (MPM) has been extensively used to solve problems involving large
displacements and deformations. The standard MPM formulation, which adopts fluid-implicitparticle (FLIP) transfer, is less dissipative but also less stable than particle-in-cell (PIC) transfer.
The affine PIC (APIC) transfer introduces the affine velocity and was developed to realize stable
simulations while overcoming dissipations of the angular momentum in PIC. This paper presents
Taylor–PIC (TPIC) transfer, which is a type of APIC transfer that combines the affine velocity
based on the first-order Taylor series approximation and PIC transfer. TPIC is simple (onlythe particle velocity gradient is required) and inherits the key advantages of the original APIC,
such as less dissipation and stability. Although TPIC does not conserve angular momentum, incontrast to APIC, the velocity gradient is preserved in the transfer between particles and the grid.
This velocity gradient contains the angular information in its skew-symmetric component, whichallows TPIC to adequately describe angular motion, similar to APIC.
Furthermore, the MPM can cause stress oscillations near boundaries when boundary conditions are explicitly imposed, i.e., when the boundary grid velocity (or momentum) is set to zero.
When affine-type transfers are used, this instability is exacerbated, and simulations can easily
fail. Therefore, we propose a kernel correction method for particles near boundaries. The proposed TPIC transfer and kernel correction are validated through five types of simulations, eachusing a different material—linear elastic, von Mises, Newtonian fluid, and Drucker–Prager. In
the simulations, incorrect results due to stress oscillations are observed even for small deformations. Applying the corrected kernels successfully removes these spurious oscillations, and theresults are consistent with the analytical solutions. Moreover, the numerical results confirm theaccuracy and robustness of the proposed TPIC transfer.
Key Words：Material point method， Particle-in-cell， Fluid-implicit-particle， Affine particle-in-cell，
Kernel correction，Boundary condition***

Senior Researcher, Foundations Group, Geotechnical Engineering Department
Head of Foundations Group, Geotechnical Engineering Department
3-1-1 Nagase, Yokosuka, 239-0826 Japan Port and Airport Research Institute
Phone: +81-46-844-5057 Fax: +81-46-844-0618 e-mail: nakamura-ke@p.mpat.go.jp

-4-

目

次

1. はじめに ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅
2. Material point method (MPM) 定式化 ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅

2.1
2.2
2.3
2.4
2.5
2.6
2.7
2.8

6

支配方程式 ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅

7

粒子から背面格子への輸送 ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅

7

他の粒子の物理量の更新 ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅

8

空間離散化 ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅

7

節点速度の更新 ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅

7

時間ステップの算出方法 ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅

8

時間離散化 ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅

7

背面格子から粒子への輸送 ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅

893. 輸送法の違いの解析的検証 ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 10
4. カーネル補正法 ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 11
5. 解析例 ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 13

5.1
5.2
5.3
5.4
5.5

線形弾性体を用いた単純せん断解析 ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 13
支持力解析 ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 16
円盤の自由回転解析 ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 18
ダムブレイク解析 ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 18
砂柱崩壊解析 ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 20

6. おわりに ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 21
参考文献 ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 26

5

中村

圭太・松村

1. はじめにFLIP 法を採用している 1, 5) ．しかし，FLIP 方は計算が不安定になりやすいため， PIC 法と FLIP 法の両者の析に広く利用されている．これらの手法の大きな利点

手法を補完係数を用いて取り入れた手法 14) もよく用は，有限要素法（FEM）で発生する要素の潰れを回避いられる．背面格子から粒子への輸送に加え，粒子かしながら，非常に大きな変形をシミュレートできるこら背面格子への輸送の手法も研究されている．速度勾とである．MPM（Material Point Method）1, 2) は，PIC（Particle In Cell）法
法

崇亮

PIC 法は散逸性が高いため，オリジナルの MPM では

メッシュフリー法は，様々な工学分野の数値計算や解

4, 5)

聡・水谷

3)

配を利用した手法 15) は，広義にはアフィン変換と呼ばと FLIP（Fluid-implicit-Particle）

れる手法である．このアフィン変換を角運動量が保存を拡張したメッシュフリー法として広く用いらされるように再定義し， PIC 法と組み合わせた手法はれてきた．MPM では，連続体は背面格子（backgroundAffine-PIC（APIC）16) と呼ばれる．本手法で明らかに

grid）を移動する物質点（material points）に離散化されなった重要な点は，アフィン変換によって得られた速度

る（以後，断りなく物質点を単に粒子と呼ぶ場合もあを用いれば，計算の安定性を失うことなく PIC 法の散逸

る）
．この粒子は，速度や応力などの物理量を有してお的な性質を解決することができるという点である．本手り，それらの情報は一度背面格子に輸送される．次に，法は，様々なシミュレーションに応用されている 16, 17) ．

輸送された背面格子上の速度・力から運動方程式を解

また，このアフィン変換のアルゴリズムを一般化したき，節点における速度を更新する．この更新された節点

moving least squares–MPM（MLS-MPM）18) も提案され速度に基づいて，FEM と同様の手続きで背面格子の変ている．APIC 法は角運動量の保存は満たすが，運動エ形を解き，物質点の変形量（ひずみ等）や速度・位置がネルギーの散逸性は解決されていないため，polynomial

更新される．最終的に，時間ステップの終了時には変形

PIC method 法 19) という手法も提案されている．

したメッシュは元の状態に戻すため，FEM で発生する

APIC 法にあるように，粒子から背面格子への輸送，ような要素の潰れは生じない．

背面格子から粒子への輸送の組み合わせは，シミュレー

MPM によるシミュレーションの精度は用いるカーネ

ションの精度と安定性に大きく影響する．オリジナルのル関数（ここでは，距離によって値が変わる重み関数

APIC 法は広く研究されているが 16, 20) ，同じくアフィの意味で用いる）に強く依存するため，多くの研究がン速度場を形成する速度勾配を用いた手法 15) は十分な

カーネル関数による補間手法に焦点を当ててきた．例え

研究がなされていない．したがって，本研究の主な目的

ば，オリジナルの MPM1) は，カーネル関数として通常の一つは速度勾配を用いたアフィン変換と PIC 法の組

FEM で用いられる線形関数を採用しているため，そのみ合わせた手法の有効性を評価することである．このア1 回微分は背面格子の要素を跨ぐ際に不連続となり結果

フィン変換は本質的には粒子速度の一次テイラー級数近として計算が不安定になる．この問題を解決するため似であるため，オリジナルの APIC 法と区別するためにに，GIMP（generalized interpolation material point）法が本稿では Taylor–PIC（TPIC）法と呼ぶ．また，本研究提案された．GIMP 法の拡張として，CPDI（convected

では，FLIP，PIC，APIC，TPIC 法の違いを，散逸（＝particle domain interpolation）と呼ばれる新たな補完手法も提案された ．近年では，total Lagrangian MPM
6)

物理情報の損失）の観点から解析的に分析も行った．

APIC 法では，B スプライン関数をカーネル関数とし

7)と呼ばれる手法も開発され，CPDI のような複雑な計算て用いており，それらの B スプライン関数は再帰的なをする必要なく問題の解決を試みている．また，B スプ関数として用いる方法 7, 9–13) との違いは，関数のピーク

ライン関数に基づく補間法も広く研究されている 8–13) ．

値が格子点上にくるようにしている 8) ．関数のピーク値

MPM の計算精度に関するもう一つの重要な点は，粒が格子上にくることで計算上扱いやすいが，境界付近で子と背面格子との間で行う情報の転送手法である．背面

応力振動が発生し不安定な挙動となることがわかってい格子から粒子への転送には，大きく分けて 2 つのアプ

る 21) ．APIC や TPIC などの手法をした場合，この不安ローチが存在する．

定性が悪化し即座に計算不能となる．この問題を解決するために，仮想粒子を境界外に配置する手法が提案され

1. グリッド速度から粒子速度を直接補間する PIC 法ている 21) ．本研究では，より直接的にカーネル関数を

2. グリッド加速度を用いて粒子速度を前進オイラー法

補正する手法を提案する．

により更新する FLIP 法

な お， 本 稿 の 内 容 は Computer Methods in Applied-6-テイラー級数近似と重みつき最小二乗法を用いた Material Point Method の高度化ゆえに，局所領域 Ω𝑝 で式 (5) は以下のように表される．

Mechanics and Engineering で出版された論文 “Taylor∑ 𝜼𝑖 ⋅ ∑ [∫ 𝜙𝑖 𝜙𝑗 𝜌 d𝑉] 𝒗𝑗̇

particle-in-cell transfer and kernel correction for material
point method” の日本語訳版である．

𝑖𝑖

2.1 支配方程式(1)

コーシー応力，⬚̇ は物質時間微分を表す．境界条件は以𝝈 ⋅ 𝒏 = 𝒕 ̄ on 𝜕Ω𝒕 ,𝒗 = 𝒗̄ on 𝜕Ω𝒗 ,𝑝1)+ ∫ 𝜼 ⋅ 𝜌𝒃 d𝑉 + ∫ 𝜼 ⋅ 𝒕 ̄ d𝐴,𝜕Ω𝒕

𝑚𝑖 = ∑ ∫ 𝜙𝑖 𝜌 d𝑉 = ∑ 𝑚𝑝 𝜙𝑖𝑝 .
𝑝

ここで，𝒇𝑖𝑛 = 𝒇𝑖

+ ∑ ∫ 𝜼 ⋅ 𝜌𝒃 d𝑉 + ∑ ∫𝑝𝜕Ω𝒕𝑝

𝜼 ⋅ 𝒕 ̄ d𝐴.

𝜼 = ∑ 𝜙𝑖 (𝒙)𝜼𝑖 .
𝑖

𝒗𝑖𝑛+1 − 𝒗𝑖𝑛= 𝒇𝑖𝑛 ,
Δ𝑡

(14)

+ 𝒇𝑖ext,𝑛 である．上付き 𝑛 は時間ス

2.4 粒子から背面格子への輸送

MPM では，粒子速度は背面格子に輸送される．この(5)過程では，オリジナルの MPM 1) のように，粒子質量でにおけるそれぞれの値 𝒗𝑖 ，𝜼𝑖 から以下のように近似で
𝑖(13)

テップを示し，Δ𝑡 = 𝑡𝑛+1 − 𝑡𝑛 である．

速度 𝒗 と試行関数 𝜼 は，形状関数 𝜙𝑖 (𝒙) を用いて，節点𝒗 = ∑ 𝜙𝑖 (𝒙)𝒗𝑖 ,

𝑚𝑛𝑖

int,𝑛

Ω𝑝

𝑝

すれば式 (9) は以下のように時間離散化できる．

うに離散化する．

∑ ∫ 𝜼 ⋅ 𝜌𝒗̇ d𝑉 = − ∑ ∫ ∇𝜼 ∶ 𝝈 d𝑉

Ω𝑝

集中質量行列 (13) を用いて，前進オイラー法を適用

式 (4) より，連続体を局所領域 Ω𝑝 を用いて以下のよ

きる．

(12)

2.3 時間離散化

2.2 空間離散化

Ω𝑝𝑝．集中質量行列は対角行列で，各行の値を足し合わせ(4)

質量保存則は自動的に満たされる 1) ．

𝑝

𝜙𝑖 𝒕 ̄ d𝐴] = ∑ 𝜙𝑖𝑝 𝑚𝑝 𝒃𝑝 + 𝒇𝑖̄ .

たものを対角上に配置することで得られる．

ここで，𝜼 はベクトル値の試行関数を示す．MPM では，

𝑝

𝜕Ω𝒕𝑝

(11)

量行列の代わりに簡便な集中質量行列がよく用いられる∫ 𝜼 ⋅ 𝜌𝒗̇ d𝑉 = − ∫ ∇𝜼 ∶ 𝝈 d𝑉

Ω𝑝

Ω𝑝

𝑝

た，𝑚𝑖𝑗 は整合質量行列と呼ばれる．MPM では整合質

れる．

𝑝

Ω𝑝

である．式 (11)–(12) では 𝜙𝑖𝑝 = 𝜙𝑖 (𝒙𝑝 ) としている．まよって，運動量保存式の弱形式は以下のように表さ

Ω

(10)

𝑝

Ω𝑝

𝒇𝑖ext = ∑ [∫ 𝜙𝑖 𝜌𝒃 d𝑉 + ∫

(3)

(8)

𝑚𝑖𝑗 = ∑ ∫ 𝜙𝑖 𝜙𝑗 𝜌 d𝑉 = ∑ 𝑚𝑝 𝜙𝑖𝑝 𝜙𝑗𝑝 ,𝑝(2)

𝜙𝑖 𝒕 ̄ d𝐴] .

(9)𝒇𝑖int = − ∑ ∫ ∇𝜙𝑖 ⋅ 𝝈 d𝑉 = − ∑ 𝑉𝑝 ∇𝜙𝑖𝑝 ⋅ 𝝈𝑝 ,𝜕Ω の単位法線ベクトルを表す．

Ω

𝜕Ω𝒕𝑝

𝑗

𝑝

ここで， ⬚̄ は境界値における値を表し， 𝒏 は物体表面ΩΩ𝑝∑ 𝑚𝑖𝑗 𝒗𝑗̇ = 𝒇𝑖int + 𝒇𝑖ext ,ここで，

ここで，𝜌 は密度，𝒃 は単位質量あたりの体積力，𝝈 は下の通りである．

Ω𝑝

𝜼𝑖 の任意性を考慮すれば以下を得る．

領域 Ω において，運動量保存則を Eulerian 形式で表𝜌𝒗̇ = ∇ ⋅ 𝝈Ω𝑝= ∑ 𝜼𝑖 ⋅ [− ∫ ∇𝜙𝑖 ⋅ 𝝈 d𝑉 + ∫ 𝜙𝑖 𝜌𝒃 d𝑉 + ∫2. Material point method (MPM) 定式化

すと次のようになる．

𝑗

重み付けした Shepard’s partition of unity 法 22) が用いら
れる．

𝒗𝑖𝑛 =

1
∑ 𝜙𝑛 𝑚 𝒗𝑛 .
𝑚𝑛𝑖 𝑝 𝑖𝑝 𝑝 𝑝(15)

APIC 法 16) では，アフィン速度 𝒗𝑝𝑛̃ を用いて以下のよ(6)

うに背面格子へ輸送される．

(7)

𝒗𝑖𝑛 =

-7-

1
∑ 𝜙𝑛 𝑚 𝒗𝑛̃ (𝒙𝑛 ).
𝑚𝑛𝑖 𝑝 𝑖𝑝 𝑝 𝑝 𝑖

(16)

中村

アフィン速度 𝒗𝑝𝑛̃ (𝒙𝑖𝑛 ) は以下のように表される．

圭太・松村

(17)

𝑛
𝑫𝑝𝑛 = ∑ 𝜙𝑖𝑝(𝒙𝑖𝑛 − 𝒙𝑝𝑛 ) ⊗ (𝒙𝑖𝑛 − 𝒙𝑝𝑛 ).

(18)

𝑖

APIC 法 16) が提案されている．APIC 法は，アフィン速度（節 2.4）を導入することで，PIC 法の散逸性を解消することに成功している．APIC 法では，新たな粒子の物理量として 𝑩𝑝𝑛 を導入しており，以下のように算出さである．APIC 法で必要となる新たな行列 𝑩𝑝𝑛 は角運動れる．

ている 16) ．

本研究で提案する TPIC 法は式 16 のようなアフィン

本手法で提案する TPIC 法も APIC と同様に，PIC 法

速度を用いるが，その速度は 1 次のテイラー級数近似を（式 (23)）による粒子速度の更新を採用している．TPIC

用いて以下のように表される．

では，粒子の速度勾配が必要となるが（節 2.4）
，速度勾配はどの手法においても変形勾配テンソル 𝑭 を算出す

(19)

ここで，∇𝒗𝑝𝑛 は粒子の速度勾配を示す．本手法は，すで15)る際に用いるため，実質的に APIC 法のような追加の物
理量を計算する必要はない．

であるが，PIC 法との組み合わせについては議論されていない．本研究では，式 (19) と2.7 他の粒子の物理量の更新

PIC 法を併用（すなわち TPIC）することによって，速他の粒子の物理量に関しては，輸送の手法の違いに関度勾配が保存されることを節 3 で明らかにしている．まわらず同様の手続きで行える．速度勾配テンソルは以下

た，TPIC の有用性については節 5 で示す．

のように得られる．

𝑛∇𝒗𝑝𝑛+1 = ∑ 𝒗𝑖𝑛+1 ⊗ ∇𝜙𝑖𝑝.

2.5 節点速度の更新𝑖

MPM では，運動方程式は背面格子で解かれるため，
節点速度は式 (14) を用いて以下のように更新できる．

𝒗𝑖𝑛+1 = 𝒗𝑖𝑛 + Δ𝑡𝒂𝑛𝑖 ,(20)

𝒂𝑛𝑖 = 𝒇𝑖𝑛 /𝑚𝑛𝑖 .

(21)

以下のように更新できる．

𝑭𝑝𝑛+1 = (𝑰 + Δ𝑡∇𝒗𝑝𝑛+1 ) ⋅ 𝑭𝑝𝑛 .

オリジナルの MPM では，粒子速度は FLIP 法更新された応力はさらに以下のように修正される．

𝝈𝑛+1= 𝝈𝑛𝑝 + 𝑻 ∶ Δ𝒅𝑝 + (Δ𝝎𝑝 ⋅ 𝝈𝑛𝑝 − 𝝈𝑛𝑝 ⋅ Δ𝝎𝑝 ), (28)𝑝

ここで，𝑻 は用いた構成則で計算される接線剛性テンソをルである．𝒅 は変形速度テンソル，𝝎 はスピンテンソル用いて以下のように更新される．

𝑛 𝑛
𝒗𝑝FLIP,𝑛+1 = 𝒗𝑝𝑛 + Δ𝑡 ∑ 𝜙𝑖𝑝
𝒂𝑖 ,𝑖

なお，PIC 法 3) では

𝑛 𝑛+1
𝒗𝑝PIC,𝑛+1 = ∑ 𝜙𝑖𝑝
𝒗𝑖 .
𝑖

(22)

を示し，それぞれ以下のように表される．

1
T[∇𝒗𝑝𝑛+1 + (∇𝒗𝑝𝑛+1 ) ] Δ𝑡2
1
TΔ𝝎𝑝 = [∇𝒗𝑝𝑛+1 − (∇𝒗𝑝𝑛+1 ) ] Δ𝑡
2
Δ𝒅𝑝 =

(23)

のように更新されるが， FLIP 法は PIC 法の散逸性を

(29)
(30)

最後に，粒子位置は以下のように更新される．
𝑛 𝑛+1
𝒙𝑝𝑛+1 = 𝒙𝑝𝑛 + Δ𝑡 ∑ 𝜙𝑖𝑝𝒗𝑖 .

低減することができるため，オリジナルの MPM では𝑖

FLIP 法を採用している．FLIP 法は散逸性が改善されているもの，数値解析の安定性はやや劣っている．そのため，両者の速度を係数 𝜉 によって補間する手法が提案さ(31)

なお，PIC 法では式 (31) は以下のように単純化できる．

𝒙𝑝𝑛+1 = 𝒙𝑝𝑛 + Δ𝑡𝒗𝑝𝑛+1 .

れている 14) ．

𝒗𝑝Blend,𝑛+1 = 𝜉𝒗𝑝FLIP,𝑛+1 + (1 − 𝜉)𝒗𝑝PIC,𝑛+1 ,(27)

客観的な応力速度に Jaumann 速度を選択した場合は，2.6 背面格子から粒子への輸送

4, 5)

(26)

この速度勾配テンソルを用いれば，変形勾配テンソルは

ここで，𝒂𝑛𝑖 は節点加速度を示し，以下のように定義さ
れる．

(25)

この行列 𝑩𝑝 は，角運動量が保存されるように定義されの際に更新される．

に提案された手法

𝑛 𝑛+1
𝑩𝑝𝑛+1 = ∑ 𝜙𝑖𝑝
𝒗𝑖 ⊗ (𝒙𝑖𝑛 − 𝒙𝑝𝑛 ).
𝑖の情報を有しており，これは背面格子から粒子への輸送

𝒗𝑝𝑛̃ (𝒙𝒊𝒏 ) = 𝒗𝑝𝑛 + ∇𝒗𝑝𝑛 ⋅ (𝒙𝑖𝑛 − 𝒙𝑝𝑛 ),崇亮

その後，PIC 法の問題点を改善する新たな手法として

𝒗𝑝𝑛̃ (𝒙𝒊𝒏 ) = 𝒗𝑝𝑛 + 𝑩𝑝𝑛 ⋅ (𝑫𝑝𝑛 )−1 ⋅ (𝒙𝑖𝑛 − 𝒙𝑝𝑛 ),

ここで，

聡・水谷

以上の輸送法の違いを図-1 にまとめる．

(24)
-8-

(32)

テイラー級数近似と重みつき最小二乗法を用いた Material Point Method の高度化

𝑛𝑚𝑛𝑖 = ∑ 𝑚𝑝 𝜙𝑖𝑝

1. 節点質量(式 (13))𝑝

2. 節点速度𝑛
𝑚𝑝 𝒗𝑝𝑛∑𝑝 𝜙𝑖𝑝

FLIP: 𝒗𝑖𝑛 =𝑚𝑛𝑖(式 (15))𝑛
𝑚𝑝 (𝒗𝑝𝑛 + 𝑩𝑝𝑛 ⋅ (𝑫𝑝𝑛 )−1 ⋅ (𝒙𝑖𝑛 − 𝒙𝑝𝑛 ))∑𝑝 𝜙𝑖𝑝

APIC: 𝒗𝑖𝑛 =𝑚𝑛𝑖(式 (16)，(17))

𝑛𝑚𝑝 (𝒗𝑝𝑛 + ∇𝒗𝑝𝑛 ⋅ (𝒙𝑖𝑛 − 𝒙𝑝𝑛 ))∑𝑝 𝜙𝑖𝑝TPIC: 𝒗𝑖𝑛 =𝑚𝑛𝑖(式 (16)，(19))

𝑛𝑛𝒇𝑖𝑛 = − ∑ 𝑉𝑝𝑛 ∇𝜙𝑖𝑝
⋅ 𝝈𝑛𝑝 + ∑ 𝜙𝑖𝑝
𝑚𝑝 𝒃𝑝 + 𝒇𝑖̄

3. 節点力

4. 節点速度の更新

5. 粒子速度の更新

𝑝𝑝𝒗𝑖𝑛+1 = 𝒗𝑖𝑛 + Δ𝑡𝒂𝑛𝑖 = 𝒗𝑖𝑛 + Δ𝑡

𝒇𝑖𝑛
𝑚𝑛𝑖

(式 (20)，(21))

𝑛 𝑛
FLIP: 𝒗𝑝𝑛+1 = 𝒗𝑝𝑛 + Δ𝑡 ∑ 𝜙𝑖𝑝
𝒂𝑖

APIC & TPIC:
5*. 追加の物理量の更新（APIC）(式 (11)，(12))(式 (22))

𝑖
𝑛 𝑛+1
𝑛+1
𝒗𝑝 = ∑ 𝜙𝑖𝑝 𝒗𝑖
𝑖(式 (23))

𝑛 𝑛+1
𝑩𝑝𝑛+1 = ∑ 𝜙𝑖𝑝
𝒗𝑖 ⊗ (𝒙𝑖𝑛 − 𝒙𝑝𝑛 )(式 (25))𝑖

𝑛 𝑛+1
𝒙𝑝𝑛+1 = 𝒙𝑝𝑛 + Δ𝑡 ∑ 𝜙𝑖𝑝
𝒗𝑖

6. 粒子位置の更新(式 (31))𝑖

図-1

FLIP・APIC・TPIC 法のまとめ

Remark 1. MLS–MPM 18) では，速度勾配 ∇𝒗𝑝𝑛+1 はア

用することが多い．弾性体における音速は以下のようにフィン速度行列 𝑩𝑝𝑛+1 ⋅ (𝑫𝑝𝑛 )−1 と同等となる．

計算される．

𝑛̂∇𝒗𝑝𝑛+1 = ∑ 𝒗𝑖𝑛+1 ⊗ ∇𝜙𝑖𝑝
𝑖𝑛 𝑛+1
= [∑ 𝜙𝑖𝑝𝒗𝑖 ⊗ (𝒙𝑖 − 𝒙𝑝𝑛 )] ⋅ (𝑫𝑝𝑛 )−1𝑖
𝑛+1= 𝑩𝑝 ⋅ (𝑫𝑝𝑛 )−1 ,

𝑐=

√
4
√
√𝐾 + 3𝐺
√

𝜌

,

(34)

ここで，𝐾 は体積弾性係数，𝐺 はせん断弾性係数，𝜌 は密度である．ゆえに，臨界時間ステップは以下のように(33)𝑛̂ は MLS 形状関数である 18, 23) ．そのため，
ここで， 𝜙𝑖𝑝

計算できる．

TPIC と APIC は MLS 形状関数を用いている場合は
MLS–MPM に帰着する．

Δ𝑡cr =ℎ,max(𝑐𝑝 + ‖𝒗𝑝 ‖)𝑝(35)

ここで，ℎ は背面格子の間隔，𝑐𝑝 は音速，𝒗𝑝 は粒子速

2.8 時間ステップの算出方法

度である．

陽的に時間積分する際は，時間ステップを決定すると
きに Courant–Friedrichs–Lewy 条件（CFL 条件）

24)

非圧縮ニュートン流体を用いるときは以下の式を用

を適

-9-

中村

いる 25) ．

ℎℎ2Δ𝑡cr = min (, ),max(𝑐𝑝 + ‖𝒗𝑝 ‖) 𝜈
𝑝

圭太・松村

聡・水谷

崇亮

示すように節点加速度を考慮していることから PIC 法に比べて散逸が少ない．

APIC 法では，粒子速度は以下のように背面格子に輸(36)

ここで，𝜈 = 𝜇/𝜌（𝜇 は粘性係数）は動粘性係数を示す．

送される．

テップは以下のように表される．

ゆえに，粒子速度は式 (23) を用いて以下のように保存

𝒗𝑖𝑛 = 𝒗𝑝𝑛 + 𝑩𝑝𝑛 ⋅ (𝑫𝑝𝑛 )−1 ⋅ (𝒙𝑖𝑛 − 𝒙𝑝𝑛 ).

クーラン（Courant number）𝐶r ≤ 1 を用いれば時間ス

Δ𝑡 = 𝐶r Δ𝑡cr .

される．

𝑛
𝑛
𝒗𝑝𝑛+1 = 𝒗𝑝𝑛 ∑ 𝜙𝑖𝑝+ 𝑩𝑝𝑛 ⋅ (𝑫𝑝𝑛 )−1 ⋅ ∑ 𝜙𝑖𝑝(𝒙𝑖𝑛 − 𝒙𝑝𝑛 )(37)= 𝒗𝑝𝑛 ,

3. 輸送法の違いの解析的検証
ここでは，これまでに紹介した輸送法（FLIP， PIC，𝑖𝑖𝑛
= 𝑩𝑝 ⋅ (𝑫𝑝𝑛 )−1 ,

(38)

𝑖

∑ 𝒙𝑖 ⊗ ∇𝜙𝑖𝑝 = 𝑰,𝑖

ここで，𝑰 は 2 階の恒等テンソルである．

APIC 法とは対照的に，TPIC 法では角運動量が保存されない代わりに速度勾配が保存される．式 (19) より，(40)

節点速度は以下のように表される．

𝒗𝑖𝑛 = 𝒗𝑝𝑛 + ∇𝒗𝑝𝑛 ⋅ (𝒙𝑖𝑛 − 𝒙𝑝𝑛 ).

(41)

速度勾配 ∇𝒗𝑝 は以下のように保存される．

𝑛𝑛𝒗𝑝𝑛+1 = 𝒗𝑝𝑛 ∑ 𝜙𝑖𝑝+ ∇𝒗𝑝𝑛 ⋅ ∑ 𝜙𝑖𝑝(𝒙𝑖𝑛 − 𝒙𝑝𝑛 ) = 𝒗𝑝𝑛 ,格子に以下のように輸送される．
𝑛
𝜙𝑖𝑝
𝑚𝑝

= 𝒗𝑝𝑛 .

𝑖𝑖

る総和記号 ∑𝑝 は消える．外力がない場合を考えれば，

𝒗𝑖𝑛+1 = 𝒗𝑖𝑛 と置くことができ，𝒗𝑖𝑛+1 = 𝒗𝑝𝑛 が導かれる．

ここで， FLIP 法では 𝒂𝑛𝑖

れているため，TPIC 法も角運動をよく表現することができる（節 5.4）
．以上の結果を概略図として示すと図-2のようになる．

MLS–MPM 18) は，∇𝒗𝑝𝑛 = 𝑩𝑝𝑛 ⋅ (𝑫𝑝𝑛 )−1 が成り立つこ(43)

式 (38) を利用した．一方，粒子の速度勾配 ∇𝒗𝑝 は以下とから，角速度と速度勾配の両者を保存する．しかし，

MLS–MPM の形状関数はオリジナルの形状とは異なっのように散逸してしまう．

𝑖

(50)

分（skew-symmetric part）に角運動に関する情報が含ま= 0 を用いて， PIC 法では𝑛∇𝒗𝑝 = 𝒗𝑝𝑛 ⊗ ∑ ∇𝜙𝑖𝑝
= 0,

𝑖

(49)

ここで，式 (23) を用いた．速度勾配は，その反対称部

ゆえに，FLIP 法（式 (22)）と PIC 法（式 (23)）の両者

𝒗𝑝FLIP,𝑛+1 = 𝒗𝑝PIC,𝑛+1 = 𝒗𝑝𝑛 .

𝑖𝑛𝑛∇𝒗𝑝𝑛+1 = 𝒗𝑝𝑛 ∑ ∇𝜙𝑖𝑝
+ ∇𝒗𝑝𝑛 ⋅ ∑(𝒙𝑖𝑛 − 𝒙𝑝𝑛 ) ⊗ ∇𝜙𝑖𝑝
= ∇𝒗𝑝𝑛 ,(42)

一つの粒子しか存在しないので，式 (13) と (15) におけにおいて，粒子速度 𝒗𝑝𝑛 は保存される．

(48)

ゆえに，𝒗𝑖𝑛+1 = 𝒗𝑖𝑛 を仮定すれば，粒子速度 𝒗𝑝 と粒子在すると仮定すると，式 (15) を用いて粒子速度は背面

𝒗𝑖𝑛 =

(47)

ていないが，APIC 法では角運動量は保存される 16) ．

ある速度 𝒗𝑝 と速度勾配 ∇𝒗𝑝 を持った粒子が一つ存𝑛𝜙𝑖𝑝
𝑚𝑝 𝒗𝑝𝑛𝑖

ここで，式 (40) と (41) を用いた．速度勾配は保存され(39)

これらの空間勾配は以下のように表される．

∑ ∇𝜙𝑖𝑝 = 0,

(46)

𝑛
𝑛
∇𝒗𝑝𝑛+1 = 𝒗𝑝𝑛 ⊗ ∑ ∇𝜙𝑖𝑝
+ 𝑩𝑝𝑛 ⋅ (𝑫𝑝𝑛 )−1 ⋅ ∑(𝒙𝑖𝑛 − 𝒙𝑝𝑛 ) ⊗ ∇𝜙𝑖𝑝と linear field reproduction の両者を満たすものとする．

∑ 𝜙𝑖𝑝 𝒙𝑖 = 𝒙𝑝 .

𝑖

いては完全には保存されない．

この際，カーネル関数は以下のような partition of unity𝑖𝑖

ここで，式 (38) と (39) を用いた．一方，速度勾配につ

APIC，TPIC）の違いについて解析的観点から考察する．
∑ 𝜙𝑖𝑝 = 1,(45)た形の関数となってしまう（節 4）．
本節で示した速度勾配の保存は，あくまで粒子が一つ(44)しかない場合について考察し，輸送法の特徴について述

ここで，式 (40) を用いた．FLIP 法と PIC 法の両者におべたものであるため，複数の粒子が存在する際は当然成いて粒子の速度勾配は散逸したが，FLIP 法は式 (22) にり立たない．しかし，速度が空間に線形に分布している- 10 -

テイラー級数近似と重みつき最小二乗法を用いた Material Point Method の高度化

< TPIC >< FLIP/PIC >particle-to-gridtransferparticle-to-gridtransfer

Lost
grid-to-particle
transfer
Particles

grid-to-particle
transfer
Particles

Grid

図-2

Grid

FLIP・PIC・TPIC 輸送の概略図

場合（速度勾配が一定）について考えれば，TPIC 法は再付近で満足されない．この問題は 8) によって修正法が提びこの保存性を持つことになる（1 次のテイラー級数近案されている（図 4a）．しかし，この方法は partition of

似を用いているため）．以上のような特徴により，TPIC

unity（式 (38)）満足されても，linear field reproduction

は散逸性が少なく，安定した計算を可能にしていると考（式 (39)）は満足されない．

えられる．

式 (38) と (39) の両者を満足するために，本研究ではカーネル関数 𝜙𝑖 (𝒙) を以下のように補正する．

𝜙𝑖̂ (𝒙) = 𝐶(𝒙, 𝒙𝑖 )𝜙𝑖 (𝒙𝑝 ).

4. カーネル補正法(53)

補 正 関 数 𝐶 は 重 み つ き 最 小 二 乗 法 （weighted least

MPM は用いているカーネル関数によっては，境界条件として節点速度（あるいは運動量）を設定すると応力squares（WLS））を用いて次のように構築する．

振動が生じることが知られている．APIC 法や TPIC 法はこれらの影響を強く受け，計算が不安定化し計算不能となる．これまで，仮想粒子を用いてこの応力振動を抑𝐶(𝒙, 𝒙𝑖 ) = P(𝒙 − 𝒙𝑝 )M(𝒙𝑝 )−1 P⊤ (𝒙𝑖 − 𝒙𝑝 ),(54)

M(𝒙𝑝 ) = ∑ 𝜙𝑖𝑝 P⊤ (𝒙𝑖 − 𝒙𝑝 )P(𝒙𝑖 − 𝒙𝑝 ).

(55)

ここで，行列 M は以下のように表される．

制する手法が提案されている 21) ．本研究では，この問題𝑖に対してカーネル関数を直接補正し， partition of unityと linear field reproduction の両者を満足させることで解

線形の多項式基底 P は 3 次元では以下となる．

P(𝒙) = [1, 𝑥, 𝑦, 𝑧],

決を図る．
この 2 つの条件は，図 3 に示すように，境界付近でな

𝑚 = 4.

(56)

次に，行列 M−1 を次のようなブロック表示に書き換ければ通常は満足される．しかし，境界付近では満足されず，この不整合さが境界付近での応力振動を生じさせえる．

る．この問題はあらゆるカーネル関数に起こりうるが，本研究では以下の 2 次のスプライン関数を用いる．

𝐶𝑪2

M−1 = [ 1

𝑪2⊤
],
𝑪3

(57)

ここで，𝐶1 はスカラー，𝑪2 は長さ 𝑑（𝑑 は空間の次元

3120 ≤ |𝜉| <⎧−‖𝜉‖ + 4
2
⎪
1 2 3
9 1
3
𝑤(𝜉) =(51)‖𝜉‖ − ‖𝜉‖ +≤ |𝜉| < ,
2
8 2
2
⎨2
⎪
3
≤ |𝜉|
⎩0
2
𝑥𝑝 − 𝑥𝑖𝑦𝑝 − 𝑦𝑖𝑧𝑝 − 𝑧𝑖
𝜙𝑖𝑝 = 𝜙𝑖 (𝒙𝑝 ) = 𝑤 (
)𝑤(
)𝑤(
),
ℎ
ℎ
ℎ
(52)

を表す）のベクトル，𝑪3 は 𝑑 × 𝑑 の行列である．最終的に，補正されたカーネル関数 (53) は以下のように単純化できる．

̂ = 𝜙𝑖̂ (𝒙𝑝 ) = 𝐶(𝒙𝑝 , 𝒙𝑖 )𝜙𝑖 (𝒙𝑝 ) = (𝐶1 +𝑪2 ⋅(𝒙𝑖 −𝒙𝑝 ))𝜙𝑖𝑝 .
𝜙𝑖𝑝(58)
以上の定式化において，物質点 𝒙𝑝 の位置は空間に固ここで，ℎ は格子間隔を示す．

定されているため，手法は WLS となる．一方，物質点 𝒙𝑝 を変数として扱う場合は，手法は移動最小二乗

図-3 に示すように，partition of unity（式 (38)）は境界- 11 -

中村

圭太・松村

聡・水谷

崇亮

Domain near boundary
Boundary of problem domain
Nodal point
Nodal point outside boundary
Focused material point

Local domain

Interior domain
Inner boundary of two domains
図-3

粒子が境界近傍に存在するときの概略図（2 次の B-スプライン関数）とから，要求されるサンプルポイントの数 𝑚 = 𝑑 + 1 と法（moving least squares（MLS））となる．この違いは，
カーネル関数の空間勾配に影響を与える

26)．

同じか常に多くなる．

Remark 2. 補正されたカーネル関数 𝜙𝑖̂ (𝒙) を用いれば，

最 小 二 乗 法 を 用 い た ア プ ロ ー チ は， reproducing

任意の連続関数 𝑢(𝒙) はその節点での値 𝑢𝑖 = 𝑢(𝒙𝑖 ) を用

kernel method 28) ，smooth particle hydrodynamics 29–31) ，いて以下のように近似できる．

𝑢(𝒙) = ∑ 𝜙𝑖̂ (𝒙)𝑢𝑖 .
𝑖

MPM 18, 23, 32) など，これまでにも粒子法の分野で広く用いられてきた．本研究ではこれを MPM の境界面に(59)おける安定化に応用している．粒子が境界面から遠い
場合は 𝐶1 = 1 と 𝑪2 = 0 が満たされ，粒子が境界面に𝑢(𝒙) = P(𝒙 − 𝒙𝑝 ) と 多 項 式 関 数 で 表 さ れ る 場 合，式 (53)–(55) と式 (59) を用いれば以下を得る．

P(𝒙 − 𝒙𝑝 ) = P(𝒙 − 𝒙𝑝 )M(𝒙𝑝 )−1 ∑ 𝜙𝑖𝑝 P⊤ (𝒙𝑖 − 𝒙𝑝 )P(𝒙𝑖 − 𝒙𝑝 )= P(𝒙 − 𝒙𝑝 ).

𝑖

(60)

ゆえに，1 次の多項式基底 P (式 (56)) を用いれば，parti-

近い場合は 𝐶1 ≠ 1 と 𝑪2 ≠ 0 が得られる 28) ．ゆえに，
式 (58) は以下のように粒子が境界面に近いときのみ用いれば良い．

𝜙Φ𝑖 (𝒙𝑝 ) = { 𝑖𝑝̂
𝜙𝑖𝑝

if 𝒙𝑝 ∈ Ωcotherwise (if 𝒙𝑝 ∈ Ω̂ c ).

(61)

式 (53) の空間勾配は以下のように表される 23) ．

tion of unity（式 (38)）と linear field reproduction（式 (39)）̂ = ∇𝜙𝑖̂ (𝒙)||∇𝜙𝑖𝑝= (𝑪2 + 𝑪3 ⋅ (𝒙𝑖 − 𝒙𝑝 ))𝜙𝑖𝑝 . (62)𝒙=𝒙の両者が満たされることがわかる．

𝑝

Remark 3. WLS，あるいは MLS のスキームを粒子法でこの補正されたカーネル関数の空間勾配は境界面付近用いている場合，式 (54) におけるモーメントマトリク

であっても式 (40) と (41) を満たす．しかし，図-4b に

ス M はサンプル点（通常は粒子）が少ないなどの悪条示すように，その形状はオリジナルのカーネル関数の空間勾配 ∇𝜙𝑖𝑝 とは異なってしまう．節 5.2 で，この僅か件により逆行列を求められない場合がある．これを回避するために，モーメントマトリクスの対角化がよく用いな差による数値計算結果の違いを示している．カーネルられる 27) ．しかし，本研究で提案している WLS を用い関数の空間勾配の形状を維持するには，図-4c に示すよ

たカーネル補正法ではモーメントマトリクス M の逆行

うに WLS の代わりに MLS を用いれば良い．しかし，

列は常に求めることができる．これは，本手法では節点

MLS を用いた際のカーネル関数の空間勾配は式 (54) のM(𝒙)−1 P(𝒙𝑖 − 𝒙) の項を 𝒙 について微分する必要があるを WLS のサンプル点として利用しているためである．
特に，線形の多項式（式 (56)）を用いている場合，物質点を取り囲む節点の数は少なくとも 2𝑑（粒子が存在しことから，式の導出が極めて複雑化してしまう 26, 28) ．

ている要素の節点数で 𝑑 は空間の次元を表す）となるこ

最も単純な解決策は，境界から遠い場合にはオリジ

ナルのカーネル関数の空間勾配 ∇𝜙𝑖𝑝 を用いることで- 12 -

テイラー級数近似と重みつき最小二乗法を用いた Material Point Method の高度化

ある．

∇𝜙∇Φ𝑖 (𝒙𝑝 ) = { 𝑖𝑝̂
∇𝜙𝑖𝑝

if 𝒙𝑝 ∈ Ωcotherwise (if 𝒙𝑝 ∈ Ω̂ c ).

きな移動を含んだシミュレーションを実施して検証している．B-スプライン関数をカーネル関数として用いている場合は，カーネル補正法を適用していない関数につい(63)ては，既往の研究でよく用いられる図-4a で示した関数

ここで，この手法の整合性を保つためには，式 (63) が連

8)えれば，任意の 𝑥𝑝 ∈ 𝜒𝑗 について ∇𝜙𝑖 (𝑥𝑝 ) = ∇𝜙𝑖̂ (𝑥𝑝 ) が続関数でなければならない．これは，1 次元について考を用いている．

5.1 線形弾性体を用いた単純せん断解析満たされる必要がることを意味している（𝜒𝑗 は B-スプ

最初の解析例では，線形弾性体を微小変形の範囲内でライン関数のノットを示し，このノットは図-3 の inner

単純せん断し，カーネル補正法の妥当性検証を行った．

ン関数では 𝑪3 =ないプレートに上下の両面で接している．下のプレー

図-6 に示すように，サンプルは 1 m × 1 m で全く滑らboundary 𝜕Ωc 上にある）．これは， 2 次の B-スプライ
𝑪3 =

3

ℎ2

4

ℎ2

𝑰， 3 次の B-スプライン関数ではトは空間に固定し，上のプレートを水平方向（𝑥 軸方𝑰 を用いれば，簡単に導くことができる（簡略化された 𝑪3 は既往の研究

16)向）に速度 1 mm/s で変位が 0.1 mm に達するまで動かで示されている)．式 (52)

した．これらの境界条件は背面格子に与えており，直接で示されているように，2 次元以降のカーネル関数は 1

節点速度 𝒗𝑖𝑛+1 を修正することで考慮している．比較と次元でのカーネル関数値の乗算で算出されることから，して，FEM による解析結果（解析解として扱う）を図-6

多次元においてもこの関数の連続性は保持される．最終的に，式 (61)，(63) を用いれば補正されたカーネル関数に示す．

は図-4d のようになる．

（𝐸 = 103 MPa）
， ポ ア ソ ン 比 （𝜈 = 0.333）， 密 度

線形弾性体の材料パラメーターはヤング率

Remark 4. ∑𝑖 𝜙𝑖𝑝 = 1 が満足されているとき，すなわ（𝜌 = 103 kg/m3 ）である．シミュレーションは，FLIP 法

ち，粒子が境界面から遠い場合は，式 (62) で示した WLS（オリジナル MPM で用いられている 1) )，APIC 法 16) ，を用いたカーネル補正法は，MLS–MPM の定式化 18) とTPIC 法，MLS–MPM 18, 23) を用いてそれぞれ実施した．

一致する 23) ．そのため，簡単のために WLS を全領域に

カーネル補正法の効果を検証するために， MLS-MPM

渡って適応している手法（図-4b）は MLS–MPM と呼以外の手法（節 4 に示したように，MLS-MPM はそも

ぶことにする．

そも補正する必要がない）では，カーネル補正法を適用

式 (61) と (63) に基けば，WLS を用いたカーネル関数

した関数とそうでない関数の両者についてシミュレーの補正は粒子が境界に近い時のみ適用すれば十分であるションを実施した．本シミュレーションでは，2 次の B-

ため，計算コストは低いと言える．なお，2 次の B-スプ

スプライン関数をカーネル関数として用いている．全てのケースについて，クラーン数は 𝐶r = 1（式 (35)）とライン関数では，MLS（図-4c）と WLS（図-4d）とで関数の形が同じであるが，これは一般的には成り立たな

して時間ステップを計算している．

い．例えば，3 次の B-スプライン関数（図-5）を見れば，

図-7 は von Mises 応力の分布を示している．全てのMLS（図-5c）のカーネル関数の空間勾配は境界面側かケースにおいて，カーネル補正法を用いていない場合

ら最初に増加して減少しているが，WLS では（図-5d）は，上下の境界付近で応力振動が見られている．これはは単調に減少する．

既往の研究 21) でも同様に確認されている．図に示すように，この応力振動により von Mises 応力の分布も正し

5. 解析例

く得られていないことがわかる．この結果は，メッシュを細かくしたり，時間ステップを短くするのでは改善さ本研究で提案した TPIC 法とカーネル補正法を 5 種類

れないことを別途確認している．一方，補正されたカーの数値シミュレーションを通して妥当性を検証する．最ネル関数（図-4d）を用いたケースでは，応力振動が解初の 2 つの解析例では，カーネル補正法をそれぞれ弾性

消され，FEM の結果と同等の結果が得られていること

体，および弾塑性体の物質で検証する．これらのシミュがわかる．この結果は，以下に示すような相対誤差を用

レーションでは，変形は微小になるようにしている．こいた結果（図-8）を見ても同様に確認できる．

れは，変形が大きくなるとカーネル補正法なしでは即座

Error = √にシミュレーションが破綻してしまうためである．次の3 つの解析例では，TPIC 法の妥当性について物質の大- 13 -(𝝈 − 𝝈∗ ) ∶ (𝝈 − 𝝈∗ )× 100 %,𝝈∗ ∶ 𝝈 ∗

(64)

中村

圭太・松村

聡・水谷

1
r (x)

(x)

1

崇亮

0.5

0

1
0

0

1

2

3
x

4

5

6

0

1

2

3
x

4

5

6

0

1

2

3
x

4

5

6

1

2

3
x

4

5

6

1

2

3
x

4

5

6

(a) Steffen et al.8)

1
r (x)

(x)

1

0.5

0

1
0

0

1

2

3
x

4

5

6

(b) MLS–MPM 18, 23)

1
r (x)

(x)

1

0.5

0

1
0

0

1

2

3
x

4

5

6

0

(c) MLS を用いたカーネル補正．

1
r (x)

(x)

1

0.5

0

1
0

0

1

2

3
x

4

5

6

0

(d) WLS を用いたカーネル補正．

図-42 次の B-スプライン関数- 14 -

テイラー級数近似と重みつき最小二乗法を用いた Material Point Method の高度化

1
r (x)

(x)

1

0.5

0

1
0

0

1

2

3
x

4

5

6

0

1

2

3
x

4

5

6

0

1

2

3
x

4

5

6

1

2

3
x

4

5

6

1

2

3
x

4

5

6

(a) Steffen et al.8)

1
r (x)

(x)

1

0.5

0

1
0

0

1

2

3
x

4

5

6

(b) MLS–MPM 18, 23)

1
∇φ(x)

φ(x)

1

0.5

0
−1

0

0

1

2

3
x

4

5

6

0

(c) MLS を用いたカーネル補正．

1
∇φ(x)

φ(x)

1

0.5

0
−1

0

0

1

2

3
x

4

5

6

0

(d) WLS を用いたカーネル補正．

図-53 次の B-スプライン関数- 15 -

中村

圭太・松村

聡・水谷

崇亮

1 mm/s1m

No-slip boundary

Grid spacing: 0.02 m
Number of particles per cell: 4 points

No-slip boundary

1m

図-6

0

20
40
60
80
von Mises stress (kPa)100

単純せん断解析における幾何学条件と境界条件（左）．FEM 解析による von Mises 応力分布（右）ここで， 𝝈∗ は FEM 解析で得られた応力テンソルである． MLS–MPM は， partition of unity と linear field（𝜈 = 0.333），密度（𝜌 = 103 kg/m3 ），粘着力（𝑐 = 10 kPa）

reproduction（式 (38)–(41)) の両者をもともと満たして単純せん断の解析例 (節 5.1) と同様に，シミュレー

である．

いるため，正しい解析結果が得られている．

ションは FLIP 法， APIC 法， TPIC 法（カーネル補正法の有無），MLS-MPM でそれぞれ実施した．カーネル

5.2 支持力解析

関数は 2 次の B-スプライン関数を用いた．クーラン数（式 (35)）は 𝐶r = 1 として時間ステップを計算した．

前節におけるシミュレーションは単純な弾性体を扱っていたため，次に弾塑性体を用いてカーネル補正法の妥

図-10 は基礎の支持力，図-11 は偏差ひずみの分布を

当性検証を行う．シミュレーションは基礎の支持力問題示している．カーネル補正法を適用しない場合は，全てを採用した．von Mises モデルを用いている場合，基礎のケースにおいて Prandtl の解析解に一致せず，地盤内の支持力は Prandtl によって以下のような解析解が導かれている 33) ．

𝑞u = (2 + 𝜋)𝑐𝐵,のすべり線もぼやけてしまっている．これらの結果は，メッシュの細かさや時間ステップの短さを変えても改(65)

ここで，𝑐 は土の粘着力，𝐵 は基礎の幅を示す．本支持力

善しないことを別途確認している．さらに，APIC 法や線対称であるため，片側半分の領域のみシミュレーショ
ンを実施する．基礎は速度 4 mm/s 一定で鉛直下向き方

たところ，時間ステップを短くすることで APIC 法でもvon Mises モデルの降伏関数は以下のように表される．

図-10 に示すように，カーネル補正法はこれらの問

TPIC 法では，シミュレーションは不安定となってしま

問題の幾何学条件と境界条件は図-9 に示す．本解析はい，APIC 法では解析が破綻した．追加で検証を実施し解析を実施することができたが，その時の結果は TPIC

向に変位が 0.4 mm に達するまで移動させた．

ℱ(𝝈) = √3𝐽2 (𝝈) − 𝜎Y = 0,法と同様に正しい解は得られなかった．
題を大きく改善し， Prandtl の解析解も得られている．

(66)

ここで，𝐽2 は偏差応力の第二普遍量で，𝜎Y > 0 は降伏

さらに，偏差ひずみの分布（図-11）では，滑り線がぼやけることなくくっきりと現れている． MLS-MPM で応力を示す材料パラメーターである．平面ひずみ条件では，降伏応力 𝜎Y はモール・クーロンの破壊規準で用いられる粘着力 𝑐 を用いて以下のように表される．

𝜎Y = √3𝑐.

は，荷重–変位関係は解析解が得られている．しかし，滑り線はカーネル補正法を適用した FLIP 法，APIC 法，

TPIC 法と比べてややぼやける結果となった．これは，
MLS-MPM ではカーネル関数の形状（図-4b）がオリ(67)

ジナルのものと異なっていることが原因として考えら塑性ひずみの発展則は関連流れ則とした．材料パラ
メーターは，ヤング率（𝐸 = 103 MPa）
，ポアソン比れる．

- 16 -

テイラー級数近似と重みつき最小二乗法を用いた Material Point Method の高度化

FLIP

APIC

TPIC

FLIP
(with kernel correction)

APIC(with kernel correction)

TPIC(with kernel correction)

0

100

図-7

20
40
60
80
von Mises stress (kPa)

単純せん断解析における von Mises 応力分布

FLIP

APIC

TPIC

FLIP
(with kernel correction)

APIC(with kernel correction)

TPIC(with kernel correction)

0

100

図-8

MLS-MPM

20

40
60
80
Relative error (%)

単純せん断解析における応力テンソルの相対誤差- 17 -

MLS-MPM

中村

圭太・松村

聡・水谷

0.5 m

崇亮

いない．

続いて，∑𝑝 𝒙𝑝 × 𝑚𝑝 𝒗𝑝 によって計算される角運動量

4 mm/sについて検証を行った．図-15 では，角速度はその解析

解 𝐼𝜔0（𝐼 は円盤の慣性モーメント）で正規化している．

TPIC では，厳密には角運動量が満たされないが，本解析結果では顕著な散逸は見られなかった．図-16 は，運Number of particles per cell: 4 points

動エネルギー ∑𝑝 𝑚𝑝 𝒗𝑝 ⋅ 𝒗𝑝 /2 をその理論値 𝐼𝜔02 /2 で正5m

Grid spacing: 0.02 m

規化した値を図示している．角運動量とは異なり，運動
エネルギーは全てのケースにおいて徐々に減少し散逸が生じる結果となった．

5.4 ダムブレイク解析
本節では，TPIC 法の有効性についてダムブレイクの5m

図-9

シミュレーションを通して検証する．基礎の支持力問題（節 5.2）で示したように，APIC 法，TPIC 法はカーネ

支持力解析における幾何学条件と境界条件

ル補正法なしでは，たとえ微小変形の範囲内でも計算が不安定化してしまう．そのため，ダムブレイクのような粒子の移動が大きいシミュレーションでは解析が即座に5.3 円盤の自由回転解析破綻する．よって，本解析例では全てのケースにおいて

本節では， TPIC 法の妥当性について回転する円盤

カーネル補正法を適用した．解析の幾何学条件と境界条の解析をもとに検証する．図-12 に示すように，円盤は半径 0.5 m の線形弾性体（𝐸 = 20 MPa， 𝜈 = 0.3，

件は図-17 に示す．全ての境界において摩擦は全く生じ𝜌 = 103 kg/m3 ）である．見かけ上の壁面（境界条件は課ないものとしている．
流体には，非圧縮のニュートン流体を用いた．この場さない）も配置することで，カーネル補正法の妥当性検合，応力ひずみ関係は以下のように与えられる．

2
𝝈 = −𝑝𝑰 − 𝜇tr(𝒅)𝑰 + 2𝜇𝒅,3証も行った．円盤は各粒子に初期速度を与えることで角
速度 𝜔0 = 2𝜋 rad/s で回転させた．粒子の初速度 𝒗0 は

𝒗0 = 𝝎0 × 𝒓0 .

(69)

ここで，𝑰 は 2 階の恒等テンソル，𝜇 は粘性係数，𝒅 は

以下のように計算できる．

変形速度テンソル，𝑝 は流体の圧力を示す．流体の圧力(68)

ここで，𝒓0 は円盤の中心を原点とした際の粒子の初期は以下に示すような Morris の式 34) を用いて計算した．

𝑝 = 𝑐2 (𝜌 − 𝜌ref ),位置ベクトルである．

(70)

ここで，𝑐 は音速，𝜌ref は基準密度を示す．本解析では，

これまでの解析例と同様に，シミュレーションは FLIP

𝜇 = 1.01 × 10−3 Pa s，𝑐 = 50 m/s，𝜌ref = 103 kg/m3 と法，APIC 法，TPIC 法（カーネル補正法の有無），MLS-

MPM でそれぞれ実施した．本解析では，より影響範囲

した．音速は，実際の音速を用いる必要はなく，流体のが広い 3 次の B-スプライン関数を用いている．クーラ
ン数（式 (35)）は 𝐶r = 0.8 とした．

非圧縮性が保たれるように決定する 35) ．
シミュレーションはカーネル補正法を用いた FLIP

図-13 は粒子速度のノルムを示している．カーネル補法，APIC 法，TPIC 法で実施した．カーネル関数には，正法を適用しない APIC 法と TPIC 法では，境界付近

2 次の B-スプライン関数と GIMP 関数 36) の 2 種類を用いた．クーラン数（式 (36)）には 𝐶r = 0.1 を用いた．

（見かけ上の壁面）で計算が不安定化し解析が破綻した．

FLIP 法では計算は可能であったが，不自然な速度の振図-18 に流体の圧力分布を示す． FLIP 法は計算が不動が見られた．一方，カーネル補正法を用いたケース
は，全てのケースにおいて解析結果の改善が見られた．

安定化しやすいため，シミュレーション初期の段階か確認できる．図-14 に示すように，FLIP 法は APIC 法，見られる．この圧力振動によって，流体の薄いシート状

これは，角速度の相対誤差 |𝜔 − 𝜔0 |/|𝜔0 | の分布を見てもら，B-スプライン，GIMP の両者において圧力の振動がの形状が破壊されてしまっている（𝑡 = 1.45 s での壁面

TPIC 法，MLS-MPM と比べて滑らかな結果が得られて- 18 -

テイラー級数近似と重みつき最小二乗法を用いた Material Point Method の高度化

Force (kN)

150

FLIP
APIC
TPIC
MLS-MPMFLIP*
APIC* (failed)
TPIC*
Prandtl solution

100

50

0

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

Displacement (mm)図-10

支持力解析における荷重–変位関係．アスタリスクはカーネル補正法を用いていないことを示す．APIC 法は計算が破綻したため示されていない．

Failed

FLIP

APIC

TPIC

FLIP
(with kernel correction)

APIC(with kernel correction)

TPIC(with kernel correction)

0

図-11

0.2e-3 0.4e-3 0.6e-4 0.8e-3 1e-3
Deviatoric strain

支持力解析における偏差ひずみ分布．

- 19 -

MLS-MPM

中村

圭太・松村

聡・水谷

崇亮

Pseudo wall
Grid spacing: 12.5 mm
Number of particles per cell: 4 points

Disk

1m

図-12

円盤の自由回転解析における幾何学条件と境界条件

上，𝑡 = 1.75 s での波の先頭付近）．対して，APIC 法と密度（𝜌 = 1.4 × 103 kg/m3 ）である．塑性流れ則には，TPIC 法ではこの圧力振動は大きく軽減されており圧力

以下に示すような非関連流れ則を用いた．

𝒢(𝝈) = √𝐽2 (𝝈) + 𝛽𝐼1 (𝝈),分布はより滑らかになっている．また，FLIP 法で破壊されてしまっていた流体の薄いシート状の形状も確認できる．なお，圧力振動を完全に抑制するには 𝑣–𝑝 法 35)

ここで，パラメーター 𝛽 はダイレイタンシー角 𝜓 を用などの他の手法の援用が必要となると考えられる．

いて以下のように表される．

𝛽=

本解析では，TPIC 法の解析結果は角運動を含めてほ
とんど APIC 法と変わらない結果となった．図-19 では，

TPIC 法が APIC 法と同等レベルで渦度を表現できていることがわかる．一方，FLIP 法では渦度も滑らかな結𝜓 = 0° とした．

.

(74)

壊シミュレーションはカーネル補正法を適用した FLIP

本節では，砂柱の崩壊解析を実施して，弾塑性モデル

法，APIC 法，TPIC 法で実施した．カーネル関数には，2を用いた大変形問題に対する TPIC のパフォーマンスの次の B-スプラインと GIMP 36) の 2 種類を用いた．カー

検証を行った．解析の幾何学条件と境界条件は図-20 にネル補正法を適用しなかった場合，APIC 法と TPIC 法示す．

では計算は即座に破綻してしまうが，その一例として計以下の Drucker–Prager モデルを砂の材料モデルとしℱ(𝝈) = √𝐽2 (𝝈) + 𝛼𝐼1 (𝝈) − 𝜅 = 0,

算破綻時の様子を図-21 に示す．クーラン数（式 (35)）
は 𝐶r = 1 を用いた．

図-22 は速度のノルムの分布を示している．FLIP 法で(71)

ここで， 𝐽2 は偏差応力の第二普遍量， 𝐼1 は応力の第一は応力振動に伴う粒子速度の不安定性が蓄積しており，普遍量，𝛼 と 𝜅 は材料パラメーターである．これらの材

不自然な挙動を示している 17) ．この不安定性は時間ステップを短くすることで抑制できる．例えば，FLIP–GIMPのケースではクーラン数をおよそ 𝐶r < 0.75 に設定す料パラメーターは，平面ひずみ条件を仮定すれば，モール・クーロンモデルの粘着力 𝑐 と内部摩擦角 𝜙 を用いて以下のように表される．

𝛼=2√9 + 12 tan 𝜓

ダムブレイクの解析例（節 5.4）と同様に，砂柱の崩5.5 砂柱崩壊解析

tan 𝜙

tan 𝜓

塑性ひずみの過剰な発生を防ぐために，本解析例では果が得られなかった．

て用いる．

(73)

√9 + 12 tan 𝜙
2

,

𝜅=

3𝑐

√9 + 12 tan 𝜙
2.

ることで改善された．一方，APIC 法や TPIC 法では不安定な挙動はほとんど改善された（APIC–GIMP では

𝑡 = 0.71 s 時にやや不安定な挙動が見られた）．

(72)

材料パラメータは，ヤング率（𝐸 = 10 MPa）
，ポアソン

比（𝜈 = 0.3）
，内部摩擦角（𝜙 = 35°）
，粘着力（𝑐 = 0 Pa），- 20 -

テイラー級数近似と重みつき最小二乗法を用いた Material Point Method の高度化

FailedFailed

FLIP

APIC

TPIC

FLIP
(with kernel correction)

APIC(with kernel correction)

TPIC(with kernel correction)

0

3.1

図-13

0.62 1.24 1.86 2.48
Norm of velocity (m/s)

MLS-MPM

円盤の回転解析における速度のノルムの分布

6. おわりに

また，大変形問題ではカーネル補正をしていないアフィン型の手法は即座に計算不能となった．このように，多本研究では，一次のテイラー級数近似に基づく速度近くの問題において，カーネル補正は適用することが望似と PIC 法を組み合わせた TPIC 法を提案した．APIC

ましい，あるいは必要不可欠と考えられる．数値計算の法で保存される角運動量は TPIC 法では保存されない

結果，補正されたカーネル関数は応力振動をうまく除が，粒子と背面格子間の輸送では速度勾配が保存され去し，結果は解析解とよく一致することが確認された．

る．したがって，TPIC 法は APIC 法と同様に角運動をAPIC 法と同様に，提案した TPIC 法は FLIP 法と PIC

適切に記述することができる．

法の両方の利点，すなわち安定性と非散逸性の両者を実

このようなアフィン型の手法は，境界付近の応力振動現している．

により不安定化し，シミュレーションが容易に破綻する

本研究では，角運動も含めて TPIC 法と APIC 法の間可能性がある．この応力振動を低減するためには，カーネル関数が partition of unity （∑𝑖 𝜙𝑖 (𝑥) = 1）と linearに大きな違いは見られなかったが，TPIC 法は APIC 法

field reproduction（∑𝑖 𝜙𝑖 (𝑥)𝑥𝑖 = 𝑥）の両方の条件を満たとは異なり運動量の保存性が厳密には満たされない．しす必要があることがわかった．そこで，本研究では WLS

スキームを用いたカーネル補正法を開発した．カーネル

る必要がある．

MPM はこれまで地盤工学分野において，地滑り 37) ，たがって，この違いの影響については今後さらに検討す補正は境界付近の粒子に対してのみ適用する必要がある

液状化 38) などの大変形問題において広く適用されてきため計算コストは少ない．

た．このような大規模なシミュレーションを行う際には提案した TPIC 法とカーネル補正法を 5 種類の数値シ解析手法のロバスト性や精度が大きく求められる．本研

ミュレーションで検証した．その結果，補正されていな究で提案した TPIC 法，カーネル補正法は，これら既往いカーネル関数を使用すると，小さな変形であっても応のシミュレーションにも容易に適用可能である．

力振動により誤った結果が得られることが確認された．

（2023 年 4 月 25 日受付）

- 21 -

中村

FLIP

聡・水谷

崇亮

FailedFailed

APIC

TPIC

0

FLIP
(with kernel correction)

圭太・松村

10

20
30
40
Relative error (%)

50

APIC
(with kernel correction)

TPIC(with kernel correction)

0

5

図-14

1

2
3
4
Relative error (%)

MLS-MPM

円盤の回転解析における角速度の相対誤差

Normalized angular momentum

1.2
FLIP
APIC
TPIC
MLS-MPM

1.1

1

0.9

0.8

0

1

2

3

4

5

Time (s)図-15円盤の回転解析における正規化された角運動量- 22 -

テイラー級数近似と重みつき最小二乗法を用いた Material Point Method の高度化

Normalized kinematic energy

1.2
FLIP
APIC
TPIC
MLS-MPM

1.1

1

0.9

0.8

0

1

2

3

4

5

Time (s)図-16円盤の回転解析における正規化された運動エネルギー

Grid spacing: 10 mm
Number of particles per cell: 4 points

Water

0.6 m

1.2 m

Free-slip boundary
3.22 m図-17

ダムブレイク解析における幾何学条件と境界条件- 23 -

中村

圭太・松村

聡・水谷

崇亮

Quadratic B-Spline(FLIP)

Quadratic B-Spline(APIC)

Quadratic B-Spline(TPIC)

GIMP(FLIP)

GIMP(APIC)GIMP(TPIC)

0

図-18

1

2
3
Pressure (kPa)

4

5

ダムブレイク解析における流体圧力分布

Quadratic B-Spline(FLIP)

Quadratic B-Spline(APIC)

-30

図-19

-14

2
18
Vorticity (1/s)

Quadratic B-Spline(TPIC)

34

50

ダムブレイク解析における渦度の分布（𝑡 = 3.4 s）- 24 -

テイラー級数近似と重みつき最小二乗法を用いた Material Point Method の高度化

0.2 m

Grid spacing: 4 mm

0.3 m

Number of particles per cell: 4 points

Sand
Friction coefficient, 0.4

1.2 m

図-20

砂柱崩壊解析における幾何学条件と境界条件

0

図-21

0.24 0.48 0.72 0.96
Norm of velocity (m/s)

1.2

カーネル補正法を用いていない場合の不安定な挙動（𝑡 = 0.224 s）

Quadratic B-Spline(FLIP)

Quadratic B-Spline(APIC)

Quadratic B-Spline(TPIC)

GIMP(FLIP)

GIMP(APIC)GIMP(TPIC)

0図-220.24 0.48 0.72 0.96
Norm of velocity (m/s)

1.2

砂柱崩壊解析における速度のノルムの分布- 25 -

中村

圭太・松村

聡・水谷

参考文献

崇亮

using b-spline basis functions. International Journal for numerical methods in engineering, Vol. 113,

1) Deborah Sulsky, Zhen Chen, and Howard L Schreyer.

No. 3, pp. 411–431, 2018.

A particle method for history-dependent materials.

11) Y Bing, M Cortis, TJ Charlton, WM Coombs, andComputer methods in applied mechanics and engi-

CE Augarde. B-spline based boundary conditions inneering, Vol. 118, No. 1-2, pp. 179–196, 1994.

the material point method. Computers & Structures,2) Deborah Sulsky, Shi-Jian Zhou, and Howard L

Vol. 212, pp. 257–274, 2019.

Schreyer. Application of a particle-in-cell method to12) Yuya Yamaguchi, Shuji Moriguchi, and Kenjiro Ter-solid mechanics. Computer physics communications,ada. Extended b-spline-based implicit material point

Vol. 87, No. 1-2, pp. 236–252, 1995.

method. International Journal for Numerical Meth-

3) Francis H Harlow. The particle-in-cell computing

ods in Engineering, Vol. 122, No. 7, pp. 1746–1769,

method for fluid dynamics. Methods Comput. Phys.,2021.

Vol. 3, pp. 319–343, 1964.

13) Zheng Sun, Yong Gan, Jun Tao, Zhilong Huang, and

4) Jeremiah U Brackbill and Hans M Ruppel. Flip: A

Xiaomin Zhou. An improved quadrature scheme in

method for adaptively zoned, particle-in-cell calcula-b-spline material point method for large-deformation

tions of fluid flows in two dimensions. Journal of

problem analysis. Engineering Analysis with Bound-

Computational physics, Vol. 65, No. 2, pp. 314–343,

ary Elements, Vol. 138, pp. 301–318, 2022.

1986.

14) Yongning Zhu and Robert Bridson. Animating sand5) Jeremiah U Brackbill, Douglas B Kothe, and Hans M
Ruppel.

as a fluid. ACM Transactions on Graphics (TOG),

Flip: a low-dissipation, particle-in-cellVol. 24, No. 3, pp. 965–972, 2005.

method for fluid flow. Computer Physics Commu-

15) J. E. Guilkey P. C. Wallstedt. Improved velocity pro-nications, Vol. 48, No. 1, pp. 25–38, 1988.

jection for the material point method. Computer Mod-

6) Alireza Sadeghirad, Rebecca M Brannon, and Jeff

eling in Engineering & Sciences, Vol. 19, No. 3, pp.

Burghardt. A convected particle domain interpola-

223–232, 2007.

tion technique to extend applicability of the material16) Chenfanfu Jiang, Craig Schroeder, Andrew Selle,point method for problems involving massive defor-

Joseph Teran, and Alexey Stomakhin. The affine

mations. International Journal for numerical meth-particle-in-cell method. ACM Transactions on Graph-

ods in Engineering, Vol. 86, No. 12, pp. 1435–1456,

ics (TOG), Vol. 34, No. 4, pp. 1–10, 2015.

2011.

17) Gergely Klár, Theodore Gast, Andre Pradhana,7) Alban de Vaucorbeil, Vinh Phu Nguyen, and Christo-

Chuyuan Fu, Craig Schroeder, Chenfanfu Jiang, andpher R Hutchinson. A total-lagrangian material point

Joseph Teran.

Drucker-prager elastoplasticity formethod for solid mechanics problems involving large

sand animation.

ACM Transactions on Graphics

deformations.

(TOG), Vol. 35, No. 4, pp. 1–12, 2016.

Computer Methods in Applied Me-

chanics and Engineering, Vol. 360, p. 112783, 2020.

18) Yuanming Hu, Yu Fang, Ziheng Ge, Ziyin Qu, Yixin

8) Michael Steffen, Robert M Kirby, and Martin

Zhu, Andre Pradhana, and Chenfanfu Jiang. A mov-Berzins. Analysis and reduction of quadrature er-

ing least squares material point method with displace-rors in the material point method (mpm). Interna-

ment discontinuity and two-way rigid body coupling.

tional journal for numerical methods in engineering,

ACM Transactions on Graphics (TOG), Vol. 37,

Vol. 76, No. 6, pp. 922–948, 2008.

No. 4, pp. 1–14, 2018.

9) Roel Tielen, Elizaveta Wobbes, Matthias Möller, and19) Chuyuan Fu, Qi Guo, Theodore Gast, Chenfanfu

Lars Beuth. A high order material point method. Pro-

Jiang, and Joseph Teran. A polynomial particle-in-

cedia Engineering, Vol. 175, pp. 265–272, 2017.

cell method. ACM Transactions on Graphics (TOG),

10) Yong Gan, Zheng Sun, Zhen Chen, Xiong Zhang, andVol. 36, No. 6, pp. 1–12, 2017.

Yu Liu. Enhancement of the material point method

20) Chenfanfu Jiang, Craig Schroeder, and Joseph Teran.
- 26 -

テイラー級数近似と重みつき最小二乗法を用いた Material Point Method の高度化

An angular momentum conserving affine-particle-in-30) Javier Bonet and Sivakumar Kulasegaram. A simpli-

cell method. Journal of Computational Physics, Vol.

fied approach to enhance the performance of smooth

338, pp. 137–164, 2017.

particle hydrodynamics methods. Applied Mathemat-21) Stephan Schulz and Godehard Sutmann. A consistentics and Computation, Vol. 126, No. 2-3, pp. 133–155,

boundary method for the material point method-using

2002.

imge particles to reduce boundary artefacts. In PARTI-

31) Rade Vignjevic, Juan R Reveles, and James Camp-

CLES VI: proceedings of the VI International Confer-bell. Sph in a total lagrangian formalism. CMC-Techence on Particle-Based Methods: fundamentals andScience Press-, Vol. 4, No. 3, p. 181, 2006.

applications, pp. 522–533. CIMNE, 2019.

32) Quoc Anh Tran, Martin Berzins, and Wojciech T.

22) Donald Shepard. A two-dimensional interpolation
function for irregularly-spaced data.

Sołowski. An improved moving least squares method

In Proceed-

for the material point method. 2019.

ings of the 1968 23rd ACM national conference, pp.

33) Wai-Fah Chen. Limit analysis and soil plasticity. El-

517–524, 1968.

sevier, 2013.

23) Keita Nakamura, Satoshi Matsumura, and Takaaki34) Joseph P Morris, Patrick J Fox, and Yi Zhu. Modeling

Mizutani. Particle-to-surface frictional contact algo-

low reynolds number incompressible flows using sph.

rithm for material point method using weighted least

Journal of computational physics, Vol. 136, No. 1, pp.

squares. Computers and Geotechnics, Vol. 134, p.

214–226, 1997.

104069, 2021.

35) Zhen-Peng Chen, Xiong Zhang, Kam Yim Sze, Lei

24) Richard Courant, Kurt Friedrichs, and Hans Lewy.

Kan, and Xin-Ming Qiu. vp material point method forOn the partial difference equations of mathematical

weakly compressible problems. Computers & Fluids,

physics. IBM journal of Research and Development,

Vol. 176, pp. 170–181, 2018.

Vol. 11, No. 2, pp. 215–234, 1967.

36) Scott G Bardenhagen and Edward M Kober. The gen-25) Songdong Shao and Edmond YM Lo. Incompress-

eralized interpolation material point method. Com-

ible sph method for simulating newtonian and non-puter Modeling in Engineering and Sciences, Vol. 5,

newtonian flows with a free surface. Advances in wa-

No. 6, pp. 477–496, 2004.

ter resources, Vol. 26, No. 7, pp. 787–800, 2003.

37) Yuya Yamaguchi, Shinsuke Takase, Shuji Moriguchi,26) Xiaoying Zhuang, Hehua Zhu, and Charles Augarde.

and Kenjiro Terada.

Solid–liquid coupled mate-

An improved meshless shepard and least squares

rial point method for simulation of ground collapse

method possessing the delta property and requiring nowith fluidization. Computational Particle Mechanics,

singul r weight function. Computational Mechanics,

Vol. 7, pp. 209–223, 2020.

Vol. 53, No. 2, pp. 343–357, 2014.

38) Kewei Feng, Gang Wang, Duruo Huang, and Feng27) YY Lu, T Belytschko, and L12710120847 Gu. A new

Jin.

Material point method for large-deformation

implementation of the element free galerkin method.

modeling of coseismic landslide and liquefaction-

Computer methods in applied mechanics and engi-

induced dam failure. Soil Dynamics and Earthquake

neering, Vol. 113, No. 3-4, pp. 397–414, 1994.

Engineering, Vol. 150, p. 106907, 2021.

28) Wing Kam Liu, Sukky Jun, and Yi Fei Zhang. Reproducing kernel particle methods. International journal
for numerical methods in fluids, Vol. 20, No. 8-9, pp.
1081–1106, 1995.
29) Javier Bonet and Sivakumar Kulasegaram. Correction and stabilization of smooth particle hydrodynamics methods with applications in metal forming simulations. International journal for numerical methods
in engineering, Vol. 47, No. 6, pp. 1189–1214, 2000.
- 27 -