ＩＳＳＮ２４３６－０３５Ｘ

国立研究開発法人海上・港湾・航空技術研究所

港湾空港技術研究所
報告
REPORT OF
THE PORT AND AIRPORT RESEARCHINSTITUTE

VOL.60No.1

June 2021

NAGASE,YOKOSUKA,JAPANNATIONAL INSTITUTE OF MARITIME,PORT AND AVIATION TECHNOLOGY

総
1.

仁・川口

浩二・藤木

高二朗・久保田峻

峻博貴・田中 敦高二朗・久保田博貴・鶴田修己

プレート境界断層デコルマ帯におけるスロースリップ発生メカニズムに関する研究
杉山 友理・橘7.

匠夢・藤木

パラペット後退型護岸に働く波圧に関する検討鈴木

6.

仁・川口 浩二・岩本

仮設被覆工の耐波安定性について
鈴木

5.

善昭

富山湾・寄り回り波の力学機構
田村

4.

雅之・栗山

2019 年台風 15 号による横浜港に襲来したうねり性の波浪
田村

3.

次

波の遡上域における海浜地形変化に及ぼす潮汐変動の影響に関する検討
伴野

2.

目

伸也・森川

嘉之

MPM–剛体シミュレーションのための最小二乗法を用いた摩擦接触アルゴリズムの開発と開端杭の地
盤への貫入挙動への適用
中村

圭太・松村

聡・水谷

崇亮

港湾空港技術研究所報告第 60 巻第 1 号（2021.6）

MPM–剛体シミュレーションのための最小二乗法を用いた摩擦接触アルゴリズムの開発と開端杭の地盤への貫入挙動への適用
中村 圭太*・松村聡**・水谷

要

崇亮***

旨

本研究では，粒子と剛体の間の相互作用をシミュレートするために，particle-to-surface 接触に基づいた MPMによる新しい摩擦接触アルゴリズムを提案する。 格子ベースの multivelocity field technique は，MPM でよく用いられる手法であるが，このようなシミュレーションを行うには以下のようないくつかの欠点がある。(1) 早期の誤接触を防ぎ，また，非貫入条件を満たすためには特別な追加処理が必要である。(2) 剛体の接触検出は，剛体を構成する粒子の配置に依存するため，その配置によっては接触検出のパフォーマンスが低下する。これらの問題を解決するために，提案アルゴリズムでは，particle-to-surface 接触に基づくペナルティ法を用いる。しかし，
粒子が剛体表面に非常に接近することがあるため，物理量が粒子から剛体内にある格子節点に輸送されてしまう可能性がある。本研究では，重み付き最小二乗法を用いることで、形状関数を構築する際に partition-of-unity を損なうことなく，この問題を効果的に解決する。提案したアルゴリズムは，数値シミュレーションと解析解および FEM の結果を比較することで検証した。その結果，提案するアルゴリズムは，摩擦接触現象に対して妥当な解を得ることがわかった。
また，本解析手法の地盤工学問題への適用例として，開端杭の地盤への貫入問題を解いた。解析では，既往の実験の再現解析を行い比較・検証を行った。その結果，まだ改善点はあるものの，杭の支持力推定，及び開端杭の先端閉塞現象の解明に有効に活用できることが確認できた。

キーワード：MPM，摩擦接触，重み付き最小二乗法，ペナルティ法，開端杭，先端閉塞現象*

地盤研究領域研究官**

地盤研究領域

主任研究官***

地盤研究領域グループ長〒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.60, No.1 (June 2021)

Development of frictional contact algorithm
using weighted least squares for MPM–rigid body simulationand its application to behavior of open-ended pile driven into the ground

Keita NAKAMURA*
Satoshi MATSUMURA**
Takaaki MIZUTANI***

Synopsis

This study presents a new particle-to-surface frictional contact algorithm for the material point method(MPM) to simulate the interaction between material points and rigid bodies. Although a grid-based multivelocity field technique is a common method that works well with MPM, there are several disadvantages
in developing such simulations: (1) The technique needs additional treatment to prevent early contact andsatisfy the impenetrability condition. (2) The method performs poorly in the contact detection of rigid bodies, as it depends on the arrangement of their constituent material points. To overcome these problems, theproposed algorithm uses the penalty method based on particle-to-surface contact. However, the possibilityof physical quantities being transferred from the particles to grid nodes located within the rigid bodies arises
since a particle can be very close to the surface. The weighted least squares approximation could effectively
handle this problem, without forfeiting the partition-of-unity property while constructing the shape function, and is thus incorporated into the MPM framework in this study. The proposed algorithm is validatedby comparing numerical simulations with analytical solutions and FEM results. The numerical results showthat the proposed algorithm produces reasonable results for frictional contact phenomena.
As an application of the proposed algorithm to the geotechnical problem, behavior of open-ended piledriven into the ground is simulated. In the simulations, we compares the computed results with existing
experiments. As a consequence, the proposed algorithm can help to estimate the bearing capacity of piles
and elucidate the plugged effect.
キーワード：MPM，frictional contact，weighted least squares，penalty method，open-ended pile，plugged

effect******Researcher, Foundation Group, Geotechnical Engineering Department
Senior Researcher, Foundation Group, Geotechnical Engineering Department
Head, Foundation Group, Geotechnical Engineering Department
3-1-1 Nagase, Yokosuka, 239-0826 Japan
Phone: +81-46-844-5057

Fax: +81-46-844-0618-4-

e-mail: nakamura-ke@p.mpat.go.jp

目

次

1

はじめに

6

2

MPM（Material Point Method）とは

73

MPM における重み付き最小二乗近似

83.1

重み付き最小二乗近似を用いた局所近似（Local approximation） . . . . . . . . . . . . . . . . . . . . . . . . . .

8

3.2

Partition of unity を用いた全体近似（Global approximation） . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

MPM 離散化の概要

10

4.1

支配方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

4.2

空間離散化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

4.3

時間離散化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

接触アルゴリズム11

5.1

接触している粒子の形状関数とその勾配 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

5.2

垂直方向の接触力と剛体の速度 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

5.3

接線方向（摩擦に伴う）の接触力 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

摩擦による接触を考慮した WLS–MPM の計算手順

12

4

5

6
6.1

形状関数とその勾配 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

6.2

粒子から格子への物理量の輸送（Particle-to-grid transfer）

12. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6.3

節点速度の更新. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

6.4

格子から粒子への物理量の輸送（Grid-to-particle transfer） . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

6.5

時間ステップの決め方 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

6.6

MPM，WLS–MPM，および Hu et al. [2018] による手法の違い . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

数値シミュレーションによる解析手法の妥当性の確認

13

7
7.1

スライドボックス . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14

7.2

傾斜面を転がる円盤 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14

7.3

べた基礎の支持力問題 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

158

解析手法の適用例：開端杭の地盤への貫入

179

結論と今後の課題

19

謝辞

20

参考文献

20

-5-

中村

圭太・松村

聡・水谷 崇亮

また，計算量とコストを削減するために，local multimesh

contact method [Ma et al., 2010] も提案された。しかし，

1 はじめにLemiale et al. [2010] が指摘するように，このアルゴリズム

メッシュフリー法は，固体力学の数値計算や解析に広く[Bardenhagen et al., 2000] は，速度などの物理量が同じ格子利用されている。これらの手法の大きな利点は，有限要素

節点に輸送されたときに 2 つの物体が接触していると判定さ法（FEM）ではメッシュが歪むような極めて大きな変形をシれるため，物体が実際に接触するより前に接触を検出してしミュレートできることである。粒子法の代表的な手法の一つ

まう問題がある。これは，特に影響範囲が大きい高次の形状である MPM（Material Point Method; Sulsky et al. [1994,

関数を用いる場合に致命的となる。この問題を解決するため

1995]）は PIC 法（Particle-in-Cell; Harlow, 1964）を拡張すに，Lemiale et al. [2010] と Nairn [2013] は体積勾配と粒子る形で開発された。MPM では，連続体をバックグラウンドの位置を格子節点へ補完し，追加の接触判定基準を用いるこ

メッシュ（背面格子）内を移動する粒子で離散化する。速度やとで改善した。一方，Ma et al. [2014] は，ペナルティ関数を応力などの物理量は，粒子からバックグラウンドメッシュに形導入することで，接触検出の改善を図った。

状関数を用いて輸送される（particle-to-grid transfer）。そのこれらの手法は接触アルゴリズムの精度を向上させるが，後，バックグラウンドメッシュでは運動方程式が解かれ，格子

非貫入条件（impenetrability condition; 2 つの物体は互いに節点で更新された速度に基づいて FEM と同様にメッシュがオーバーラップしない）が破られる可能性がある。これは，粒変形する。最終的に，このメッシュの変形に基づいて各粒子の子の運動はバックグラウンドメッシュを用いて陰的に解かれ物理量が更新される（grid-to-particle transfer）。なお，時間るため，格子の節点で非貫入条件を満たしていても，粒子につステップの終了時には，変形したメッシュを元の状態に戻すた

いても同じように非貫入条件が満足されるとは限らないため

め，FEM のようにメッシュが潰れることなく，大変形までシである。Huang et al. [2011] はこの問題を解決するために，

ミュレート可能である。MPM シミュレーションの精度は，粒変形後の格子だけでなく時間ステップ開始時にリセットされ子–格子間の物理量の輸送（particle-to-grid/grid-to-particle

た格子についても非貫入条件を課した。また，本問題に対処

transfer）に大きく依存するため，これまでに多くの改良手する簡単な方法の一つとして，ペナルティ法の適用がある。

法が提案されてきた。例えば，オリジナルの MPM（Sulsky

ペナルティ法は厳密には非貫入条件を満たさないため，制約

et al. [1994]）では，形状関数の微分が連続ではないため，粒としては弱いが，非貫入条件が満足されるまで制約を課し続子がバックグラウンドメッシュの要素間を移動する際に誤差けることができるため，陰的に粒子の移動が解かれる MPMが生じる（cell-crossing problem）。Bardenhagen and Koberとは相性が良い。Hamad et al. [2017] は，ペナルティ法を用[2004] は，GIMP (Generalized Interpolation Material Point)

いて，CPDI[Sadeghirad et al., 2011] における接触判定法を法を開発することでこれを解決した。Steffen et al. [2008]

改善した。ペナルティ法は簡単であるため広く用いられていは，Smooth Particle Hydrodynamics (SPH; [Monaghan and

るが，計算結果はペナルティ係数の選択に大きく依存する。

Lattanzio, 1985]) 等で広く使用されているカーネル関数であ

そのため，拡張ラグランジュ法（Simo and Laursen [1992],

る B-スプライン関数を用いた MPM シミュレーションをChen et al. [2017]）などの他のアプローチも提案されている。

行った。GIMP の拡張として，Sadeghirad et al. [2011] は

上記の MPM の接触アルゴリズムは，連続体は物質点でConvected particle domain interpolation（CPDI）を開発し構成されている。これに対して，2 つの物体のうち一方を点た。これらの手法は，粒子–格子間の物理量の輸送の精度をではなく面で構成する接触アルゴリズムも研究されており，向上させるが，角運動量は依然として保存されない [Jiang

これは，比較的軟らかい材料と剛体との相互作用をシミュ

et al., 2015]。この問題を解決するために，Klár et al. [2016]

レーションするのに適している。このようなアルゴリズムは，は affine transfer[Jiang et al., 2015] を適用し，砂のシミュ

MPM–FEM（Lian et al. [2011]; Chen et al. [2015]; Cheon

レーションを行なった。Hu et al. [2018] は，さらに，重みつき最小二乗近似（WLS; weighted least squares）を用いて，and Kim [2018]）や，MPM–DEM[Liu et al., 2018] の連成シ
ミュレーション等にも用いられる。Hu et al. [2018] は MPM

より一般的な形式の affine transfer アルゴリズムを開発した。

と剛体の連成シミュレーションを行なっている。

一方，Tran et al. [2019] は，移動最小二乗近似 (MLS; moving

本研究では，MPM と剛体の間の摩擦接触を考慮可能なシ

least squares) を適用し，要素内でにガウス求積を行う手法を

ンプルな手法を提案する。これまでのアルゴリズム（例えば，開発した。

Bardenhagen et al. [2000]）では，バックグラウンドメッシュ

摩擦接触現象は，様々な工学問題において必要性が生じるを用いて接触を検出していたが，提案するアルゴリズムではため，MPM においても，これまでにいくつかの摩擦接触ア

粒子と剛体表面の間で直接ペナルティ法を用いて接触検出を

ルゴリズムが提案されてきた。Bardenhagen et al. [2000] に行う（particle-to-surface contact algorithm）。そのため，影よって提案された格子ベースの multivelocity field technique

響範囲の大きい高次の形状関数であっても正確な接触検出が（節点において，2 物体の速度をそれぞれ分けて扱う）は，クー可能となる。しかし，この方法では，物理量が粒子から剛体

ロン摩擦モデルを用いて摩擦挙動を表現した先駆的な研究で内に位置する格子節点に輸送されてしまう可能性があり，形ある。Bardenhagen et al. [2000] はこのアルゴリズムに改良

状関数の Partition of Unity を維持するには（物理量を過不を加え，接触判定として接触面法線方向の圧縮力を考慮した。

足なく粒子から節点へ輸送するために必要），そのような無効-6-MPM–剛体シミュレーションのための最小二乗法を用いた摩擦接触アルゴリズムの開発と開端杭の地盤への貫入挙動への適用

Background mesh

Material points(particles)(i) Particle-to-grid transfer(ii) Deform grid(iii) Grid-to-particle transfer(iv) Reset grid

図 1: MPM の概略図

はテンソル積を表す（例えば，a ⊗ b = ai bj ）。

な節点を単純に除外するだけでは解決できない。そこで，本研究では，WLS（weighted least squares; 重み付き最小二乗法）近似 [Hu et al., 2018] を用いて有効な節点に対して形状

2 MPM（Material Point Method）とは関数の再構築を行う。最小二乗原理を用いた形状関数の構築
は，element-free Galerkin methods (EFG; Belytschko et al.

MPM は粒子法の一つであるが，FEM のようなメッシュ[1994]) や reproducing kernel particle method (RKPM; Liu（バックグラウンドメッシュ）も用いる。粒子とバックグラ

et al. [1995]) などの他のメッシュフリー法で広く用いられて

ウンドメッシュの両者間でうまく物理量を輸送しあうことでいる。提案手法は particle-to-surface アルゴリズムに基づい計算を行うが，バックグラウンドメッシュは運動方程式を解ているため，MPM–DEM や MPM–FEM の連成シミュレー

くための物理量の一時保存場所に過ぎないため，連続体は粒ションとも相性が良い。

子で離散化され，速度や応力などの物理量も粒子に保存され本論文は以下のように構成されている。第 2 では，MPMる。具体的な手順の概略図を図 1 に示す。図にあるように，手法の概略について，第 3 節では MPM における WLS 近

まず，速度や応力などの物理量が，粒子からバックグラウン

似の適用方法について，第 4 節では MPM の離散化方法に

ドメッシュに形状関数を用いて輸送される（particle-to-grid

ついて，第 5 節では新しい particle-to-surface アルゴリズム

transfer）。その後，バックグラウンドメッシュでは，運動方の詳細について，第 6 節では提案した摩擦接触解析法を用い

程式が解かれ，格子節点で更新された速度に基づいて FEM

た WLS–MPM の具体的な実装方法について説明する。第 7と同様にメッシュは変形する。最終的に，このメッシュの変節では，提案した定式化を 3 種類の数値シミュレーションに形に基づいて各粒子の物理量が更新される（grid-to-particle

よって検証する。第 8 節では，提案した摩擦接触解析法の工transfer）。なお，時間ステップの終了時には，変形したメッ

学分野への適用例について紹介する。最後に，第 9 節では，

シュを元の状態に戻すため，FEM のようにメッシュが潰れる

結論と今後の課題を述べる。なお，本論文の第 3–7 の内容にことなく，大変形までシミュレート可能である。
具体的なシミュレーションの例として，水の Dam break とついては，Nakamura et al. [2021] に掲載されている。
本論文で使用されている表記と記号は以下の通りである。

砂〜粘土の自由落下シミュレーションをそれぞれ図 2, 3 に示

太字はベクトルと行列を表し，記号「·」は隣接する添字のす。図にあるように，適切な構成則（応力–ひずみ関係）を与

1 つの縮約を表す（例えば，a · b = ai bi , c · a = cij aj ,えることで，様々な現象をシミュレートすることができる。

c · d = cij djk ）。 記号「:」は隣接する添字の 2 つの縮約を表

特に，通常の FEM ではメッシュの潰れによって計算が困難

し（例えば，A : B = Aij Bij , C : A = Cijkl Akl ），記号「⊗」

な条件であっても計算可能である点がこの粒子法の利点の一つである。一方で，決して万能な数値解析手法ではなく，微

-7-

中村

圭太・松村

聡・水谷 崇亮（1）WLS を用いた局所近似（local approximation）と，（2）

Partition-of-unity 法 Shepard [1968] に基づく局所近似を組み合わせた全体近似（global approximation）の 2 つの部分に分けられる。

3.1 重み付き最小二乗近似を用いた局所近似（Local approximation）
着目する領域 Ω において，座標 xp にある粒子からなる連続体があるとすると，WLS 近似では，関数 u の xp 付近の局所分布を次のように近似できる。

ũp (x) =m∑cj (xp )pj (x − xp ) ≡ c(xp ) · p(x − xp )(1)j

ここで，p は多項式の基底，c は係数ベクトルである。例えば，3 次元では，線形の多項式は次のように表される。

pT (x) = [1, x, y, z],

m=4

(2)

ベクトル c は，以下のような 2 乗誤差の和を最小化することで得られる。

Jp (c) =∑2

wi (xp ) [c(xp ) · p(xi − xp ) − ui ](3)i

図 2: Dam break シミュレーション

ここで，wi (xp ) = w(xp − xi ) は重み関数，ui はサンプル点

xi における関数 u の値である。重み関数は，xi から xp までの距離が大きくなるほど小さくなれば良く，B-spline 関数のようなよく知られた関数を適用することができる。ここで，(a) c = 0 kPa, ϕ = 30◦

wi (xp ) ̸= 0 を満たす局所領域 Ωp を定義する（一般に MPM
で は直方体である ）。この領域は，Element-free Galerkin（EFG; Belytschko et al. [1994]）法や SPH 法 [Monaghanand Lattanzio, 1985] における粒子の支配領域とみなすこと(b) c = 20 kPa, ϕ = 10◦もできる。図 4 に示すように，格子節点は領域に固定されたサンプルポイント xi の役割を果たす。
式 (3) の Jp を最小化することで c(xp ) が得られる。

(c) c = 100 kPa, ϕ = 0◦c(xp ) =図 3: Drucker–Prager モデルを用いた砂〜粘土の自由落下シ∑

ui wi (xp )qi (xp )(4)i

ミュレーション

ここで，qi (xp ) = M −1 (xp ) · p(xi − xp )(5)

また，行列 M （moment matrix）は，小変形理論で収まる範囲であれば，通常の FEM を用いた方が精度や計算速度の両者において有利である。

M (xp ) =∑w(xp − xi )p(xi − xp ) ⊗ p(xi − xp )(6)i

3 MPM における重み付き最小二乗近似である。式 (4) を式 (1) に代入すると，以下の式が得られる。

本 節 で は ，Hu et al. [2018] に よ っ て 提 案 さ れ た WLSũp (x) =∑ϕi (x)ui(7)（weighted least squares; 重み付き最小二乗法）を用いた

i

MPM について紹介する。この手法は Hu et al. [2018] は

ここで，ϕi (x) は以下のような WLS による形状関数である。

MLS–MPM（moving least squares material point method）

ϕi (x) = wi (xp )qi (xp ) · p(x − xp )と名付けていたが，多項式の係数が MLS のような任意の点ではなく，固定点によって決まることから，本研究では WLS–

(8)

この空間勾配 ϕi は，

MPM と呼ぶ。MPM における WLS 近似は，以下のように∇ϕi (x) =-8-∂p(x − xp )∂ϕi (x)= wi (xp )qi (xp ) ·
∂x
∂x

(9)

MPM–剛体シミュレーションのための最小二乗法を用いた摩擦接触アルゴリズムの開発と開端杭の地盤への貫入挙動への適用

problem domain
nodal pointmaterial point
focused material point

local domain

図 4: 解析領域 Ω，粒子（物質点）xp の局所領域 Ωp ，局所領域内の格子節点 xi 。WLS 近似は Ωp 内のサンプル点 xi に対して行わ

1

1

0.8

0.8

0.6

0.6

φ(x)

φ(x)

れる。

0.4
0.2
0

0.4
0.2

0

1

2

3

4

0

5

0

1

2

1

1

0.5

0.5

0
−0.5
−1

3

4

5

3

4

5

x

∇φ(x)

∇φ(x)

x

0
−0.5

0

1

2

3

4

−1

5

x

0

1

2
x

(a) Quadratic B-spline weight function(b) Cubic B-spline weight function

図 5: 1 次多項式と B-スプラインの重み関数を用いた WLS の形状関数とその勾配。

近似では，次のようなアプローチでこの問題に対処できる。

となる。
この WLS 近似の重要な点は，節点 xi を領域に固定された

例えば，1 次元の n 次の多項式の基底を考えた場合，粒子のサンプル点として局所的な近似に使用することである。これ

局所領域 Ωp は，WLS 近似のために少なくとも n + 1 個の節により，他のメッシュフリー法のような近傍粒子探索の必要点が必要である。より具体的に，粒子が境界上に位置すると性がなくなり，計算コストが大幅に削減される。さらに，係きに領域 Ωp 内の節点数が n + 1 であれば，n 次の多項式近似数 c が x ではなく xp に依存することで，式 (1) の p(x − xp )が関数 u を正確に補間することになるため，形状関数の値はのみが x に依存することになり，勾配 ∇ϕi (x) の表現が単純

境界上に位置している粒子に支配される。結果として，図 5化される。実際，1 次の多項式 (2) を式 (1) に適用すると，に示すように，境界上では Kronecker delta property が満足

式 (9) における ∇x p(x − xp ) は一定になる。これにより，することになる。このように，WLS–MPM は，他のメッシュ

重み関数の勾配 ∇x wi (xp ) を計算する必要がなくなる [Huフリー法（例えば Mukherjee and Mukherjee [1997]）のよう

et al., 2018]。

に境界で特別な処理を必要とせず，オリジナルの MPM と同

一般に WLS 近似は，Kronecker delta property u(xi ) ̸= ui

様に境界条件を扱うことができる。

を失う（節点での形状関数の値が 1 ではなくなる）ため，境界条件の扱いが困難となる。しかし，本定式化における WLS

-9-

中村

圭太・松村

聡・水谷 崇亮

4.2 空間離散化

3.2 Partition of unity を用いた全体近似（Global
approximation）

式 (16) から，連続体の領域は局所領域 Ωp に分割される。

∑∫

これで，Ωp における局所近似 ũp (x) が揃ったため，続いて，

u の局所近似をつなぎ合わせ，領域全体での近似を構築する。

Ωpp

局所近似は重なり合っており，一般的に互いに連続していなη · ρv̇ dV =

∑∫
p

+

いため，適切な手順をもってつなぎ合わせる必要がある。本∑∫p

研究では，重み付き平均により局所近似を全領域に拡張でき

ũ(x) =

ψp (x)ũp (x) =p

ψp (x)c(xp ) · p(x − xp ) (10)できる。

w(xp − x)ψp (x) = ∑q w(xq − x)(11)

ϕi (x)ηi(19)

したがって，局所領域 Ωp に対する式 (17) は次のようになる。

∑

(12)

ηi ·

[
∑ ∫

i

=

j

∑

[∫

MPM では，関数の再構成に up ではなく，局所近似 ũp (x) を]
ϕi ϕj ρ dV v̇jΩp∇ϕi · σ dV +ϕi t̄ dAϕi ρb dV +
∂Ωtp

Ωp

Ωp

]

∫

∫

ηi ·

i

(20)

用いることである。つまり，オリジナル MPM では，ũ(x) =∑ηi は任意の値を取り得るから次式を得る。
∑
mij v̇j = fiint + fiextp ψp (x)up を用いて，式 (12) を適用している。

(21)

j

4 MPM 離散化の概要

ここで，

本節では，MPM の離散化の概要を簡単に説明する。本手mij =法とオリジナルの MPM との違いは，WLS による形状関数∑∫pを定式化に適用したことのみである。

fiint =
fiext =

領域 Ω において，空間表示で力のつり合い式は次のように

ρv̇ = ∇ · σ + ρb

=

(13)

σ · n = t̄ on ∂Ωt

Ωp

∫

(23)

]

ϕi ρb dV +Ωp

Vp ∇ϕi (xp ) · σpp

ϕi t̄ dA∂Ωtp

ϕi (xp )mp bp + f¯i(24)と呼ばれる。MPM では，計算の単純化のため，整合質量マト
リクスの代わりに集中質量マトリクス（lamped mass matrix）(14)

on ∂Ωvを使用する。集中質量マトリクスは対角行列であり，対角線(15)上の各要素は整合質量マトリクスの対応する行の和である。

¯ は与えられた境界上での値を表す。n は，であり，記号上の □mi =境界表面上の単位法線ベクトルである。

∑∫
p

力のつり合い式の弱形式は以下となる。

∫
η · ρb dV +Ω∑

∇ϕi · σ dV =である。mij は整合質量マトリクス（consistent mass matrix）

˙ は物質時間微分を表す。境界条件は，
シー応力，記号上の □∇η : σ dV +

(22)

p

ここで，ρ は密度，b は単位質量あたりの体積力，σ はコー∫

mp ϕi (xp )ϕj (xp )

p

[
∑ ∫
∑

∑

Ωp

p

与えられる。

ϕi ϕj ρ dV =∑∫p

4.1 支配方程式Ω∑

η = η̃p (x) =WLS–MPM とオリジナルの MPM の大きな違いは，WLS–

η · ρv̇ dV =

(18)

i

w(xp − x)mpψp (x) = ∑q w(xq − x)mq∫ϕi (x)vi

i

として粒子の質量によって重み付けをした関数を用いている。

Ω∑

v = ṽp (x) =v = v̄

η · t̄ dA∂Ωtpp

ぞれ速度 v と試行関数 η を補間し，以下のように表すことが

例えば，Sulsky et al. [1994] と Hu et al. [2018] は，重み関数

∫

Ωp

p

ここで，

∑∫

WLS 近似により，節点速度 v と節点試行関数 η から，それ内挿関数 ψp を定義すると，全体近似は次のように表される。

∑

η · ρb dV +(17)

る Shepard 法 [Shepard, 1968] を用いる（図 6）。Shepard の∑∇η : σ dVΩpϕi ρ dV =Ωp∑

mp ϕi (xp )

(25)

p

4.3 時間離散化

η · t̄ dA∂Ωt

式. (21) に前進オイラー法を用いて時間離散化し，式 (25)(16)

の集中質量マトリクスを用いれば次式を得る。

ここで，η は任意のベクトル値の試行関数である。なお，MPM
では質量保存が自動的に満たされる [Sulsky et al., 1994]。

vi,n+1 = vi,n + ∆t
- 10 -

ext
int+ fi,n
fi,n
mi,n

(26)

MPM–剛体シミュレーションのための最小二乗法を用いた摩擦接触アルゴリズムの開発と開端杭の地盤への貫入挙動への適用

local approximation

global approximation

local domain
図 6: 局所近似から再構築した全体近似。

ここで，添え字 n は時間ステップを表し，∆t = tn+1 − tn とに表現するために，例外として接触した粒子を囲む格子節点なる。

については，剛体内にあっても有効化（active）する。

5 接触アルゴリズム

5.2 垂直方向の接触力と剛体の速度

本節では，MPM–剛体連成シミュレーションのためのペ法線方向の接触力は，粒子 p と剛体表面の間の距離に基づ

ナルティ法に基づく新しい接触アルゴリズムを提案する。ペいて計算される。距離 dp （図 8）を用いて，接触距離 d を次ナルティ法は厳密に非貫入条件（impenetrability condition;のように求めることができる。

2 つの物体は互いにオーバーラップしない）は満たさない

d = dp − d0 ,が，その簡便さから数値計算に広く応用されている。提案(27)

ここで，d0 は接触判定における閾値である（d = dp − d0 < 0するアルゴリズムの手順は以下の通りである。（1）最初に，の場合，接触が発生する）。この閾値は，2lp を粒子径とす

particle-to-surface 接触アルゴリズムに従って，法線方向のると，lp で与えることができる。1 次元では，例えば Np 個

接触力 f c,nor を粒子で計算し，その後，その力を格子節点にの粒子が h の間隔の格子の間に並んでいる場合，粒子径は輸送する。この方法は，純粋な格子ベースのアルゴリズム（例2lp = h/Np となる。厳密には，粒子サイズは現在の粒子サイ

えば Bardenhagen et al. [2000]) と比較して，特に高次の形

ズでなければならず，また，粒子の形状も球形ではないが（例状関数を使用する場合には，より正確な接触検出が可能であ

えば，2 次元では長方形），計算のしやすから，ここでは初期る。しかし，この場合，粒子が剛体表面に極めて近くまで接の粒子径を使用し，球形の粒子を用いる。なお，効率的に接近するため，物理量が粒子から剛体内にある格子節点に輸送触している粒子を見つけるには，局所領域 Ωp 内に無効化さされる可能性がある。この問題に対しては，形状関数を構築れている節点が存在する粒子のみを，式 (27) を用いてチェッする際に，剛体内の節点を除外し，重み付き最小二乗法で形クすれば良い。

状関数を再構築することで対処する。局所領域 Ωp で WLS

法線方向の接触力は，貫入量 d に基づいて次のように計算近似を行えば，Partition of Unity は保持されるため，粒子–される。

格子間の物理量の伝達を適切に行うことができる。（2）法線

fpc,nor = −k nor dnp ,方向の接触力とは異なり，クーロン摩擦に基づく厳密な摩擦(28)

制約を満たすために，接線方向成分 f c,tan を格子節点で直接

ここで k nor は，材料の剛性に依存するペナルティ係数であ計算する。（3）式 (26) における節点速度の更新は，計算される。節 7.1 での数値シミュレーションでも示されているよう

c,nor
た接触力 fic = fi
+ fic,tan を用いてさらに修正する必要に，計算結果はペナルティ係数 k nor に大きく依存する。適切がある。以上の詳細な説明は以降の節で行う。

なペナルティ係数を得るには，ペナルティ法に基づく他の方法と同様に感度解析等が必要である。
c,nor

最後に，接触している粒子から，法線方向の接触力 fi

5.1 接触している粒子の形状関数とその勾配

rigidと剛体の速度 vi
が以下のように輸送される。

剛体に接触している粒子は，粒子の局所領域 Ωp と剛体が重

fic,nor =なるため，形状関数とその勾配を適切に処理する必要がある。

∑

ϕi (xp )fpc,nor(29)

ψi (xp )vprigid

(30)

p

図 7 に示すように，剛体内の格子節点は無効化（inactive），

virigid =すなわち形状関数（式 (8) と (9)，または式 (37) と (38)）の∑p

構築の際は除外する必要がある。ただし，粒子の移動を適切- 11 -

中村

圭太・松村聡・水谷 崇亮

active nodal point
inactive nodal point
contacted material pointmaterial point

local domain
rigid surface
図 7: 局所領域 Ωp における形状関数の構築の際に除外する格子節点。

particle

したがって，MPM と剛体の双方向の連成解析を行うためには，別の形式が適切である。例えば，最も簡単な形としては，

fistick = −k tan v̂itan ∆t が考えられる。ここで，k tan は接線方向のペナルティ係数を表す。より詳細な実装については，例
えば Wang and Chan [2014] が参考になる。

su
rigid

6 摩擦による接触を考慮した WLS–MPM

rface

の計算手順

図 8: Particle-to-surface 接触法の概略図。dp − d0 < 0 のと

この節では，摩擦による接触を考慮した WLS–MPM の詳

き接触と判定される。

細な計算手順を紹介する。

6.1 形状関数とその勾配
5.3 接線方向（摩擦に伴う）の接触力

重み関数 w については，本研究では 3 次の B-スプライン関数を使用する。

粒 子 と 剛 体 の 間 に 摩 擦 が な い 場 合 ，接 線 成 分 は ゼ ロ
c,tan

（fi

(

= 0）であるため，ノード i における接触力 f c は，wi (x) = N

fic,nor となる。一方，摩擦による接触を考えた場合は，接線方向の成分はクーロン摩擦に基づいて評価する。まず，粒子ように表される。

(31)

y − yih)(Nz − zih)(35)

0 ≤ |ξ| < 1
1 ≤ |ξ| < 2 .
2 ≤ |ξ|(36)

位置 xp,n における形状関数 (8) とその勾配 (9) は，それぞそのため，接線方向は以下となる。

れ次のように表される。

(32)

c,nor
/∥fic,nor ∥ である。ゆえに，接線方向の力
ここで ni = fi
fistick は以下のように計算される。

mi v̂itan
fistick = −
∆t

(
N


3|ξ|3 − 6|ξ|2 + 4
1
N (ξ) =
(2 − |ξ|)36

0の接線力を計算する。粒子と剛体表面の相対速度は，以下のv̂itan = v̂i − (v̂i · ni )ni)

ここで h は格子間隔で N は以下のように定義される。

の速度と剛体の速度が一致する粘着状態（stick condition）

v̂i = vi − virigidx − xihϕi (xp,n ) = wi (xp,n )qi (xp,n ) · p0(37)∇ϕi (xp,n ) = wi (xp,n )qi (xp,n ) · ∇p0(38)

ここで p0 = p(x − xp,n )|x=xp,n である。Hu et al. [2018] でも指摘されているように，WLS–MPM では，形状関数の空間(33)

勾配 ∇ϕi を計算する際に，重み関数の空間勾配 ∇wi (x) で計算する必要がないため，計算コストを削減することができる。

クーロン摩擦を考慮すると，次のような接線力（摩擦力）が導かれる。

() fistick
fic,tan = min µ∥fic,nor ∥, ∥fistick ∥∥fistick ∥6.2 粒子から格子への物理量の輸送（Particle-togrid transfer）(34)

なお，上の粘着条件の際の力は，剛体の運動が制御され，与え式 (26) を解くためには，質量マトリクス mi ，内力ベクトられていることを前提に計算されている（one-way coupling）。

ル fiint ，外力ベクトル fiext をそれぞれ式 (25)，(23)，(24)- 12 -

MPM–剛体シミュレーションのための最小二乗法を用いた摩擦接触アルゴリズムの開発と開端杭の地盤への貫入挙動への適用で計算する必要がある。式 (26) の残りの項 vi,n は全体近ここで，T は構成則の接線係数テンソルである。d と ω はそ

似（global approximation）によって評価することができる

れぞれ変形速度テンソルとスピンテンソルであり，次のよう（節 3.2 参照）。

に表される。

vi,n = ṽn (xi ) =∑ψp (xi )Cp,n · p(xi − xp,n )]1[
T∇vp,n+1 + (∇vp,n+1 ) ∆t
2
]
1[
T
∆ωp =
∇vp,n+1 − (∇vp,n+1 ) ∆t
2
∆dp =

(39)

p

ここで Cp,n は d × m の係数マトリクスである（d は解析の

(46)
(47)

最後に，粒子の位置は次のように更新できる。

次元）。
粒子が剛体と接触している場合，それに対応する法線方向
c,norの接触力 firigidと剛体の速度 vi

xp,n+1 = xp,n + ∆tvp,n+1(48)

も，それぞれ式 (29) と式 (30) を t = tn のときの値を用いて計算する必要がある。

6.5 時間ステップの決め方

6.3 節点速度の更新

陽解法による時間積分で安定した解を得るためには，CFL（Courant–Friedrichs–Levy condition）[Courant et al., 1967]

節点速度 vi は，式 (26) を用いて更新される。接触が検出を考慮する必要がある。その時間増分は次のように与えらされた場合，すなわち法線方向の接触力が節点に輸送されたれる。

場合，式 (26) で更新された節点速度に基づいて，クーロンモ
c,tan

デル（式 (34) 参照）で接線方向（摩擦）の接触力 fic,nor

算する必要がある。節点の接触力は fic = fi

∆tcr =を計+ fic,tan とvi,n+1 ← vi,n+1 + ∆t

fic
mi,n

(49)

p

ここで，h は格子間隔，cp と vp はそれぞれ粒子 p の音速として求めることができる。最後に，速度を次式のように更新する。

hmax(cp + ∥vp ∥)

速度を表す。音速は次のように計算できる。

√

(40)

c=

K + 43 Gρ(50)

ここで，K は体積弾性率，G はせん断弾性率，ρ は密度を表

6.4 格 子 か ら 粒 子 へ の 物 理 量 の 輸 送（Grid-toparticle transfer）

す。節 7 で示す数値シミュレーションでは ∆t = 0.5∆tcr を用いている。

式 (40) を用いて節点の速度を更新した後，粒子の物理量を6.6 MPM，WLS–MPM，および Hu et al. [2018] に更新する。まず，係数行列 C は次のように更新できる。

Cp,n+1 =∑vi,n+1 ⊗ wi (xp,n )qi (xp,n )よる手法の違い(41)
同様の定式化は，1 次の多項式と，粒子の質量で重み付けさiれた関数（式 (12) 参照）を用いて，Hu et al. [2018] によっ

これを用いて更新された粒子の速度は以下のように得られる。

て導かれている。したがって，粒子が境界にも剛体にも近接

vp,n+1 = ṽp,n+1 (xp,n )= Cp,n+1 · p0∑=
vi,n+1 ϕi (xp,n )

していない場合には，提案した定式化は Hu et al. [2018] のものに帰着される。本研究の貢献は，境界条件や剛体との摩擦接触を扱うためのアルゴリズム [Hu et al., 2018] の一般化

(42)

i

を含む。MPM [Sulsky et al., 1994] と WLS–MPM 定式化の比較を表 1 にまとめる。

また，その空間勾配は次式となる。

∇vp,n+1 = ∇ṽp,n+1 (xp,n )= Cp,n+1 · ∇p0∑=
vi,n+1 ⊗ ∇ϕi (xp,n )7 数値シミュレーションによる解析手法の妥当性の確認(43)i

本節では，提案した接触アルゴリズムを検証するために，更新された節点速度を用いれば，変形勾配テンソルは以下の3 種類の数値シミュレーションを行った。最初のシミュレー

ように更新される。

Fp,n+1 = (I + ∆t∇vp,n+1 ) · Fp,nションでは，ペナルティ係数 k nor と格子の解像度の影響を検証するために，スライドボックスを用いた。2 つ目のシミュ(44)

レーションでは，円盤が傾斜面上を転がる現象を再現する。
客観応力速度に Jaumann 応力速度を選択した場合，応力は

この解析例は，摩擦接触アルゴリズムを検証するために広次のように更新できる。

く用いられている [Bardenhagen et al., 2000, 2001, Huang

σp,n+1 = σp,n + T : ∆dp + (∆ωp · σp,n − σp,n · ∆ωp ) (45)et al., 2011, Hamad et al., 2017]。3 つ目は，提案したアルゴリズムでべた基礎の支持力問題を解き，有限要素シミュレー
ションと比較する。

- 13 -

中村

圭太・松村聡・水谷 崇亮

表 1: MPM と WLS–MPM における定式化の違い

MPM

WLS–MPM

Shape function, ϕip

wi (xp )wi (xp )qi (xp ) · p0

Shape function gradient, ∇ϕip

Particle position, xp,n+1∇wi (xp )∑ϕip mp∑p
Vp ∇ϕip · σp∑p
ϕip mp bp∑p
p ψp (xi )vp∑ fi
ϕip
vp,n + ∆t i
mi∑i vi,n+1 ⊗ ∇ϕip∑
xp,n + ∆t i vi,n+1 ϕip

Coefficient matrix for WLS, Cp,n+1–

Grid mass, mi
Grid internal force, fiint
Grid external force, fiext
Grid velocity, vi
Particle velocity, vp,n+1
Gradient of the particle velocity, ∇vp,n+1*

wi (xp )qi (xp ) · ∇p0∗∗∗∑

p ψp (xi )Cp · p(xi − xp )∑i vi,n+1 ϕip∗
xp,n + ∆tvp,n+1∑
i vi,n+1 ⊗ wi (xp,n )qi (xp,n )The equation is the same as that of the MPM

7.2 傾斜面を転がる円盤

7.1 スライドボックス
図 9 はスライドボックスシミュレーションのモデルを示し

図 12 にシミュレーションの形状を示す。斜面は水平方向ており，サンプルは 0.1 m × 0.1 m の正方形である。下の板に対して角度 θ で傾斜しており，その傾斜面は三角形の物体

は摩擦係数が µ の粗い板で，他の板は滑らかな板とした。シの 1 辺として表現される。円盤が剛体の場合，傾斜角 θ と摩ミュレーションでは，最初に上板が速度 vy で所定の変位が得

擦係数 µ の組み合わせにより，円盤が転がり落ちる場合と滑られるまで下降し，その後，下板が速度 vx で滑るようになっり落ちる場合の解析解を求めることができる。

ている。計算時に応力振動を抑えるため，上板の速度は 0.1に，0 から 0.5mm/s まで 0.1 秒間直線的に制御され，一定にx′c (t) ={x′0 + 12 gt2 (sin θ − µ cos θ), tan θ > 3µtan θ ≤ 3µx′0 + 13 gt2 sin θ,

保たれる。摩擦力 Fx は，下板を除くすべての板の表面が滑

ここで，x′c は x′ 軸に沿った重心の位置である。

秒で 0 から 0.5mm/s まで直線的に制御し，0.8 秒で一定にした後，0.1 秒で減少させた。同様に，下板は図 10 に示すよう(slip)(stick)(51)

らかであることから，垂直方向の力 Fy を用いて Fx = µFy

シミュレーションでは，半径 R = 0.05 m の円盤を，ヤングとして求めることができる。垂直方向の力は，サンプルが要率 E = 103 kPa，ポアソン比 ν = 0.3，密度 ρ = 103 kg/m3

素として振る舞うと仮定して計算した。

の線形弾性体として表現している。シミュレーションは，33

粒子には，ヤング率 E = 10 kPa，ポアソン比 nu = 0.3，種類の格子間隔 h = 0.4R, 0.2R, 0.1R で行った。また，円盤密度 103 kg/m3 の線形弾性体を使用した。また，下板の摩の解像度の影響を調べるために，格子の 1 要素あたりの粒子

擦係数は ¥ mu = 0.3 とした。ペナルティ係数は材料の剛性

数 Np を変化させた。傾斜角は，θ = 30◦ とした。摩擦係数

やグリッドの解像度に依存するため，3 つのペナルティ係数は µ = 0 (slip) と µ = 0.3 (stick) とした。ペナルティ係数はk nor = 1 × 103 , 2 × 103 , 3 × 103 kN/m と，2 種類の解像度
h = 0.02 m, h = 0.005 m で計算結果を比較した。初期状態

k nor = 2 × 103 kN/m である。
図 13 には，slip ケースの重心位置 x′c が表示されている。

において，1 要素あたりの粒子数 Np は 4 とした。

図 13a では，メッシュを粗くした場合（h = 0.4R）でも解析

図 11 は，底板に作用する垂直力と摩擦力の大きさを示し解と非常によく一致しているが，メッシュを細かくすることている。これらの力は節点接触力 fic の総和として計算した。

でさらに精度を上がることがわかる。さらに，1 要素あたり

図 11a より，k nor = 1 × 103 では，理論値よりも値が低く，の粒子数を増やすことでも計算精度は改善される（図 13b）。

k nor = 3 × 103 kN/m では大きくなっている。もっとも良い
結果は，k nor = 3 × 103 kN/m の時に得られた。また，図 11b

図 14 は，stick ケースの質量中心位置 x′c を示している。slipに示すように，細かいメッシュを用いることでも精度の向上

シュでは精度が低いが（図 14a），細かいメッシュでは解析を図ることができる。この場合，どのペナルティ係数の値に解とよく一致していることがわかる。また，図 14b に示すよ

対しても，解析解との良好な一致が見られる。このように，良うに，粒子数を増やすことでシミュレーションの精度が向上

い精度の結果を得るためには，十分に細かいメッシュを用いする。

ケースの場合とは異なり，stick ケースの場合は，粗いメッ

るか，感度解析によってペナルティ係数 k nor に適切な値を選

以上のように，本解析例では，提案された接触アルゴリズ

択する必要がある。今後，より適切な法線方向の接触力の算

ムが滑り接触を精度良くシミュレートできることがわかった。

出方法や，ペナルティ係数 k nor を合理的に決定する方法を検一方，摩擦条件を適用した場合，その結果は空間の解像度に討する必要がある。

大きく依存する。例えば，図 14a の精度の低い結果は，円盤- 14 -

MPM–剛体シミュレーションのための最小二乗法を用いた摩擦接触アルゴリズムの開発と開端杭の地盤への貫入挙動への適用

background grid

rigid body

図 9: スライドボックスシミュレーションの計算モデル。下板を除いたすべての板において摩擦はゼロとしている。図における粒子の配置は，粗いメッシュ（h = 0.02 m）の場合と同じである。

Velocity (mm/s)

10
|vy |
|vx |

7.3.1

地盤の自重を無視した場合の基礎の支持力

Prandtl により与えられた支持力の解析解は以下である。

5
qu = (2 + π)cB(52)

ここで，c は土の粘着力，B は基礎幅である。MPM およ0

び FEM によるシミュレーションでは。線形弾性体を用いた

0

0.5

1
Time t (s)1.52

von Mises モデルを適用している。von Mises モデルの降伏関数はF (σ) =

図 10: スライドボックスシミュレーションにおける上下の板の速度。

√

3J2 (σ) − σY = 0(53)

ここで J2 は偏差応力の第 2 不変量，σY > 0 は降伏応力を表す材料定数である。平面ひずみ条件では，降伏応力 σY は，Mohr–Coulomb の破壊基準における粘着力 c と次のように関の空間分解能が低いこと（つまり，円盤がでこぼこしている）

連づけられる。

で円盤が斜面上で跳ねてしまうのが理由である。スリップの

σY =

√

3c

(54)

場合には，ディスクが平面上で滑る（転がらない）ため，こ塑性ひずみの発展則には，関連流れ則を用いた。

のような挙動は起こらない。提案する接触アルゴリズムでは，粒子と剛体表面の間で直接接触判定を行うため，粒子の配置計算では，応力振動を抑えるために，図 16 に示すようや数が格子ベースの接触アルゴリズム [Bardenhagen et al.,に，0.1 秒で 0 から 0.5mm/s まで直線的に制御し，総変位が2000, 2001] よりも計算結果に影響を与える可能性がある。

0.6mm になるまで一定に保つようにした。この速度は，粒子の加速度が無視できるほど小さくなるように決定されているので，本研究では，支持力に対する載荷速度の影響は無視して

7.3 べた基礎の支持力問題いる。地盤の材料パラメータは，ヤング率 E = 106 kPa，ポ

本節では，得られた解を商用ソフトウェア PLAXIS を用

アソン比 ν = 0.3，密度 ρ = 103 kg/m3 ，粘着力 c = 10 kPa

いた有限要素シミュレーションの解と比較することで，提案とした。MPM，FEM ともに，基礎と地面の接触面で摩擦はしたアルゴリズムの検証を行った。まず，土の重さを無視し働かない。なお，自重がない場合，摩擦の有無によらず支持て von Mises モデルを適用した支持力問題を用いて，MPM

力は同じ値となる。ペナルティ係数は k nor = 2 × 106 kN/mと FEM のシミュレーションの両者を評価する。そこでは，とした。

この問題の解析解としてよく知られている Prandtl 解を用い

図 17 は，土の自重を無視した場合のべた基礎の支持力をる。さらに，土の自重を考慮し，摩擦がある基礎とない基礎示している。MPM と FEM の両者において，Prandtl 解よの両方で，MPM の解を FEM の解と比較する。この問題にりもやや硬い挙動を示している。この結果の精度を上げるたは理論値が存在しないため，FEM の結果を解析解とみなす。

めには，図 17 のようにより細かいメッシュを用いる必要が

図 15 にベタ基礎の支持力問題のモデルと境界条件を示す。

ある [Fatemizaddeh et al., 2015]。また，偏差ひずみの分布

本シミュレーションは左右対称であるため，基礎の片側半分（図 18）は，Prandtl 型の破壊メカニズムで予想される滑り線

のみをシミュレートしている。

と一致している。

- 15 -

圭太・松村聡・水谷 崇亮

Vertical force Fy (kN)

8
6
4
2
0

0

0.5

1

1.5

2
Frictional force Fx (kN)

Frictional force Fx (kN)

Vertical force Fy (kN)中村

8
k nor = 1 × 103
k nor = 2 × 103
k nor = 3 × 103
analytical solution

6
4
2
0

0

0.5

1
Time t (s)

1.5

2

8
6
4
2
0

0

0.5

1

1.5

2

1.5

2

8
k nor = 1 × 103
k nor = 2 × 103
k nor = 3 × 103
analytical solution

6
4
2
0

0

(a) Rough mesh (h = 0.02)

0.5

1Time t (s)(b) Fine mesh (h = 0.005)

図 11: 下板に作用する鉛直力と摩擦力の絶対値。

background grid

付けることができる材料パラメータで，平面ひずみの状態では，次のようになる。

α= √

tan ϕ
2

9 + 12 tan ϕ

,

κ= √

3c
9 + 12 tan2 ϕ(56)

塑性ひずみの発展則には，関連流れ則を用いた。
この場合も同様に，計算での応力振動を抑えるために，図 19に示すように，基礎の速度を 0.1 秒で 0 から 0.5mm/s まで直線的に制御し，総変位が 4mm になるまで一定に保つように

rigid body

した。地面の材料パラメータは，ヤング率 E = 106 kPa，ポアソン比 ν = 0.3，密度 ρ = 103 kg/m3 ，粘着力 c = 10 kPa，
内部摩擦角 ϕ = 30◦ とした。また，垂直応力 σy と水平応力

σx の初期値は以下のように設定した。
図 12: 傾斜面上の円盤の転がりシミュレーションにおける計σy = ρgH,

算モデル。図における粒子の配置は，h = 0.2R のメッシュの場合と同じである。

σ x = σy

ν
1−ν

(57)

ここで，H は地表面からの深さ，g は重力加速度である。FEM
では，フーチングと地面の接触を境界条件とし，水平方向の拘束は，基礎と地盤との間が滑らかな条件では自由，粗い条
件では固定としている。一方，MPM では，摩擦係数 µ = 0

7.3.2 地盤の自重を考慮した場合の基礎の支持力

を滑らかな状態で，µ = 0.6 を粗い状態で使用した。ペナル

MPM と FEM シミュレーションの両者が Prandtl 解と良ティ係数は k nor = 2 × 106 kN/m とした。

い一致を示したため（図 17），次に，提案された MPM の摩図 20 に，基礎の支持力を表す。基礎と地盤との間が，滑擦接触アルゴリズムを FEM シミュレーションと比較して検らかな状態でも，粗い状態でも，MPM の結果は FEM の結証する。計算モデルと境界条件は図 15 と同じだが，土の重量

果とよく一致している。図 21 に示された偏差ひずみ分布でを考慮して，構成則に Drucker–Prager モデルを適用した点は，滑らかな条件では Hill 型の破壊メカニズム，粗い条件でが異なる。Drucker–Prager モデルの降伏関数はF (σ) =√J2 (σ) + αI1 (σ) − κ = 0

は Prandtl 型の破壊メカニズムとなっており，Griffiths et al.

(55)[2002] による報告と同様に合理的な破壊メカニズムが観察された。

こ こ で I1 は 応 力 の 第 一 不 変 量 で あ る 。α と κ は Mohr-

Coulomb モデルにおける粘着力 c と内部摩擦角 ϕ に関係- 16 -

50

50

40

40
Position x0c (cm)Position x0c (cm)

MPM–剛体シミュレーションのための最小二乗法を用いた摩擦接触アルゴリズムの開発と開端杭の地盤への貫入挙動への適用

30

20

10

0

0

0.1

0.2

0.3

0

0.4

0

0.1

0.2

0.3

0.4

0.3

0.4

50
h = 0.4R, Np = 4
h = 0.2R, Np = 4
h = 0.1R, Np = 4
analytical solution

h = 0.4R, Np = 4
h = 0.4R, Np = 16
h = 0.4R, Np = 64
analytical solution

40
Position yc0 (cm)40
Position yc0 (cm)

20

10

50

30

20

10

0

30

30

20

10

0

0.1

0.2
Time t (s)

0.3

0

0.4

0

(a) 格子解像度の影響

0.1

0.2
Time t (s)(b) 1 要素あたりの粒子数の影響

図 13: 円盤の転がりシミュレーションにおける slip 条件での重心位置。

とから，解析では 5mm の格子幅とし，肉厚は 10mm とした。地盤材料モデルは，内部摩擦角 38◦ ，粘着力 c = 0kPa の

8 解析手法の適用例：開端杭の地盤へのDrucker–Prager モデルを用いた。流れ則は非関連流れ則で，ダイレイタンシー角は 0◦ である。解析における杭の貫入速貫入
本節では，提案した解析手法の適用例として，開端杭の地度は，計算時間との兼ね合いから 0.1m/s とした。これは実盤への貫入実験 [中村ら, 2019] の再現解析を行った。解析モ

験に比べてはるかに速いが，本解析で用いる Drucker–Prager

デルと境界条件を図 22 に示す。図にあるように，解析は軸モデルは時間依存性を有する構成則ではない。また，貫入速

対称問題とし，杭を剛体として扱うことで貫入問題を解いた。

度は一定としているため，加速度が計算結果に与える影響も

その際，杭の内側と外側の土（MPM の粒子）の領域を完全小さい。影響が考えられる点としては，静的解析とは異なり，に区別するためには，杭の肉厚は最低でも格子の 2 要素分を貫入に伴い MPM 粒子は鉛直下向きの速度が支配的となるた確保する必要がある。仮に 2 要素未満の幅で肉厚を設定した

め，杭と地盤との相対速度が小さくなり，結果として摩擦力場合は，杭の内側と外側の粒子が同一の格子節点に物理量をがやや小さく評価されてしまう可能性が考えられる。静的解輸送してしまう可能性があるため，特別な場合分けを行って析は陰解法が必要となるため，本研究の対象ではないが，貫別々の節点として処理されなければならない。本研究では，入速度の検証としてさらに速度が 10 倍速い 1m/s でも貫入をそのような複雑な取扱いはせずに，肉厚を 2 要素分取るこ行い，比較することとした。杭と地盤との摩擦係数には 0.53ととした。また，実験で用いた杭の肉厚は 2mm であるが，を用いた。
前節までで提案した解析手法では，法線方向の接触力 fpc,nor

その半分の 1mm を格子幅とすると計算時間が膨大であるこは，ペナルティ係数 k nor に基づいて式 (28) 計算を行ってい- 17 -

圭太・松村聡・水谷 崇亮

50

50

40

40
Position x0c (cm)Position x0c (cm)

中村

30

20

10

0

0

0.1

0.2

0.3

0.4

0

0.5

0

0.1

0.2

0.3

0.4

0.5

0.4

0.5

50
h = 0.4R, Np = 4
h = 0.2R, Np = 4
h = 0.1R, Np = 4
analytical solution

h = 0.4R, Np = 4
h = 0.4R, Np = 16
h = 0.4R, Np = 64
analytical solution

40
Position yc0 (cm)40
Position yc0 (cm)

20

10

50

30

20

10

0

30

30

20

10

0

0.1

0.2
0.3
Time t (s)

0.4

0

0.5

(a) 格子解像度の影響

0

0.1

0.2
0.3
Time t (s)(b) 1 要素あたりの粒子数の影響

図 14: 円盤の転がりシミュレーションにおける stick 条件での重心位置。

た。しかし，地盤に強制的に杭を貫入するような問題は，計ついては，著者らがいくつか行った数値計算では，2 次元 1 要

算条件が極めて厳しいため，適切なペナルティ係数を与える

素内に 4 粒子でおおよそ 0.6–0.9 程度が良好な解を与えるとのが困難である。そのため，本研究では，Wang and Chan

いう結果を得たが，今後より詳細な検証が必要である。なお，[2014] を参考にして，以下のような式をもって法線方向の接触力 fpnor を計算した。

本解析では ξ = 0.6 を用いた。

fpnor = −(1 − ξ)2mp
(dp − d0 )np(∆t)2

図 23 に，実験と解析で得られた変位–荷重関係を示す。解析では，Drucker–Prager モデルのような極めてシンプルなモ
デルを用いているものの，得られた荷重値は実験値と同等の(58)

オーダーを示しており，解析手法の有効性が認められる。実ここで，ξ(0 ≤ ξ < 1) は，接触力の大きさを決めるパラメー

験で得られた荷重は解析結果に比べてやや大きい値を示し，タである。上式では粒子の質量 mp を用いているが，MPM

また荷重の上昇率も大きいことから，実験ではより杭の先端の粒子はあくまで連続体の代表点であるため，粒子の質量を閉塞現象が顕著に現れていることがわかる。杭の先端閉塞現用いて接触力を計算するのは厳密には適切ではない。しかし，

象や支持力は，土のダイレイタンシーの影響を受けることか

1 要素内で分担する粒子数の数が大きくない場合は十分適用

ら [Peng et al., 2014]，解析のパフォーマンスをより向上させ可能であることが，著者らの検証でわかっている。ただし，初るには，ダイレイタンシーを適切に考慮したモデルを適用す

期粒子配置の際に，1 要素に対する粒子数を変更した場合は，る必要があると考えられる。なお，本解析ではダイレイタン

連続体としての質量の分布に変化がなくても，粒子の質量 mp

シー角を 0◦ としているが，通常の Drucker–Prager モデルのは変わるため，同じパラメータ ξ では同様の解が得られない

場合は 0◦ 以外の値にするとダイレイタンシーが生じ続ける計

可能性があることに注意が必要である。パラメータ ξ の値に算になってしまうため，限界状態モデルなどの適用が必要と- 18 -

MPM–剛体シミュレーションのための最小二乗法を用いた摩擦接触アルゴリズムの開発と開端杭の地盤への貫入挙動への適用

0.5 m

grid spacing: 0.05 m
number of particles per cell: 4 points

5m

footing (rigid body)5m

Velocity |vy | (mm/s)

図 15: 支持力問題での解析モデルと境界条件0.6値を示していることから，貫入速度の影響は小さいと考えられる。
図 24 と 25 に，それぞれ杭貫入時の平均応力と偏差応力分0.4布を示す。平均応力分布より，杭の貫入に伴い杭内部の平均応力が上昇し，先端閉塞現象が生じていることがわかる。ま0.2

た，杭の先端閉塞時に観測される杭内の地表面が低下する様

0子も確認できる。せん断応力分布を見ると，杭先端以深に楔

0

0.2

0.4

0.6 0.8
Time t (s)

1

1.2

1.4

形のせん断領域が生じていることがわかる。この領域はしば

図 16: 基礎の鉛直速度（土の重量を無視したケース）

象が生じる原因の一つと考えられている。今後，土の内部摩しば実験でも観測されており [菊池ら, 2008]，杭の先端閉塞現
擦角や粘着力，杭の肉厚や根入れ，杭径等を変化させてパラ
メトリックスタディを行うことで，杭の先端閉塞現象のメカ

60
Force (kN)

ニズム解明に役立つと考えられる。

40
MPM
MPM (fine mesh)
FEM
Prandtl solution

20
0

0

0.1

0.2
0.3
0.4
Displacement (mm)0.5

9 結論と今後の課題
本研究では、MPM のための particle-to-surface 接触アルゴリズムを提案した。精度の高い接触検出を実現するために，提案手法では，multivelocity field technique（節点での 2 物

0.6

体の速度をそれぞれ分けて扱う） の代わりに、粒子と剛体表面の間で直接ペナルティ法を用いた。物理量が粒子から剛体

図 17: 基礎の変位–支持力関係（土の重量を無視して，von内に位置する格子節点に輸送される可能性がある問題点は，Mises モデルを用いたケース）。MPM のメッシュが細かい

重み付き最小二乗法（WLS–MPM）を用いて解決した。このケースは，格子間隔が 0.025m である。

方法では，partition-of-unity を損なうことなく形状関数を構成できる。WLS-MPM は，境界におけるデルタ関数の特性を保持できるため，境界条件はオリジナルの MPM と同様と考なる。貫入速度の影響については，貫入速度 1m/s では荷重えることができる。

の変動値が大きくなる結果となったが，これは時間ステップ

提案したアルゴリズムは 3 つの数値シミュレーションで検をより細かく取れば緩和されると考えられる。値が変動する証した。その結果，すべてのシミュレーションにおいて，提ことを除けば，貫入速度 1m/s と 0.1m/s の荷重値は同程度の- 19 -

中村

圭太・松村聡・水谷 崇亮

Velocity |vy | (mm/s)

図 18: 基礎の支持力問題における偏差ひずみ分布（土の重量を無視して，von Mises モデルを用いたケース）

0.6

法の検証をより進めていくことで，杭の支持力推定に活用可能になると考える。また，連続体解析を行っている本手法では，これまで行われてきた DEM（discrete element method）

0.4

によるシミュレーションとは異なり，土の構成則を適用できるだけでなく，地盤内の応力状態を可視化することができる

0.2

ため，パラメトリックスタディを実施していくことで，杭の0

先端閉塞現象のメカニズム解明に役立てることができると考

0

2

4
6
Time t (s)

8

10

えられる。
（2021 年 6 月 4 日受付）

図 19: 基礎の鉛直速度（Drucker–Prager モデルのケース）謝辞案したアルゴリズムは解析および FEM による計算結果とよ

本研究は，日本鉄鋼連盟，及び JSPS 科研費（19K15093）

く一致した。しかし，スライドボックスシミュレーションでの助成を受けたものです。

は，所定の精度の結果を得るためには，ペナルティ法に基づ
く他のアルゴリズムと同様に，十分に細かいメッシュを用い

参考文献

るか，感度解析を用いてペナルティ係数の適切な値を選択する必要があることがわかった。法線方向の接触力の計算や，S. Bardenhagen, J. Brackbill, and D. Sulsky. The material-

ペナルティ係数については，より適切な手法を今後検討する

point method for granular materials. Computer methods

必要がある。円盤の転がりシミュレーションでは，本アルゴ

in applied mechanics and engineering, 187(3-4):529–541,
2000.

リズムが粒子と剛体表面の間で直接接触検出を行うため，格子ベースの接触アルゴリズムよりも，粒子の配置やその数がS. Bardenhagen, J. E. Guilkey, K. Roessig, J. Brackbill,

提案アルゴリズムの結果に大きく影響を与える可能性がある

W. Witzel, and J. Foster. Improved contact algorithm

ことわかった。べた基礎支持力問題のシミュレーションでは，

for the material point method and application to stresspropagation in granular material. 2001.

提案したアルゴリズムが，粗い基礎と滑らかな基礎の支持力

S. G. Bardenhagen and E. M. Kober. The generalized in-の違いや，地盤の破壊メカニズムをうまくシミュレートでき

terpolation material point method. Computer Modeling

ることを実証した。

in Engineering and Sciences, 5(6):477–496, 2004.

提案したアルゴリズムの適用例として，本研究では，開端杭の地盤への貫入シミュレーションを行った。解析では，

T. Belytschko, Y. Y. Lu, and L. Gu. Element-free galerkin

Drucker–Prager モデルのような極めてシンプルなモデルを

methods. International journal for numerical methods
in engineering, 37(2):229–256, 1994.

用いているものの，得られた荷重値は実験値と同等のオーダーを示しており，解析手法の有効性が認められた。現場のZ. Chen, X. Qiu, X. Zhang, and Y. Lian. Improved coupling

支持力推定は経験に基づいているため，近年顕著となっていof finite element method with material point method
based on a particle-to-surface contact algorithm. Com-る杭の大口径化には対応が困難なケースもある。今後，本手- 20 -

MPM–剛体シミュレーションのための最小二乗法を用いた摩擦接触アルゴリズムの開発と開端杭の地盤への貫入挙動への適用

600
Force (kN)

Force (kN)

600
400
200
0

MPM
FEM
0

1

2
3
Displacement (mm)

400
200
0

4

MPM
FEM
0

(a) smooth (MPM: µ = 0)1

2
3
Displacement (mm)4(b) rough (MPM: µ = 0.6)

図 20: WLS–MPM と FEM における基礎の支持力の比較(b) rough (µ = 0.6)(a) smooth (µ = 0)

図 21: 偏差ひずみ分布：（a）滑らかな基礎における Hill 型の破壊メカニズム（b）粗い基礎における Prandtl 型の破壊メカニズム

puter Methods in Applied Mechanics and Engineering,pacity of rough rigid strip footing on cohesive soil: prob-

293:1–19, 2015.

abilistic study. Journal of Geotechnical and Geoenvironmental Engineering, 128(9):743–755, 2002.

Z.-P. Chen, X. Zhang, X.-M. Qiu, and Y. Liu. A frictional
contact algorithm for implicit material point method.

F. Hamad, S. Giridharan, and C. Moormann. A penalty

Computer Methods in Applied Mechanics and Engineer-function method for modelling frictional contact in mpm.
Procedia Engineering, 175:116–123, 2017.

ing, 321:124–144, 2017.
Y.-J. Cheon and H.-G. Kim. An efficient contact algo-

Y. Hu, Y. Fang, Z. Ge, Z. Qu, Y. Zhu, A. Pradhana, andrithm for the interaction of material particles with finite

C. Jiang. A moving least squares material point method

elements. Computer Methods in Applied Mechanics andwith displacement discontinuity and two-way rigid bodyEngineering, 335:631–659, 2018.

coupling. ACM Transactions on Graphics (TOG), 37(4):
1–14, 2018.

R. Courant, K. Friedrichs, and H. Lewy. On the partial dif-

P. Huang, X. Zhang, S. Ma, and X. Huang. Contact algorithms for the material point method in impact and pen-ference equations of mathematical physics. IBM journalof Research and Development, 11(2):215–234, 1967.

etration simulation. International journal for numerical

S. Fatemizaddeh, F. Hamad, and C. Moormann. Effects

methods in engineering, 85(4):498–517, 2011.

of the mpm discretisations on soil-structure problems.
In PARTICLES IV: proceedings of the IV International

C. Jiang, C. Schroeder, A. Selle, J. Teran, and A. Stom-

Conference on Particle-Based Methods: fundamentals

akhin. The affine particle-in-cell method. ACM Transactions on Graphics (TOG), 34(4):1–10, 2015.

and applications, pages 635–645. CIMNE, 2015.

G. Klár, T. Gast, A. Pradhana, C. Fu, C. Schroeder,

D. Griffiths, G. A. Fenton, and N. Manoharan. Bearing ca-- 21 -

中村

圭太・松村聡・水谷 崇亮

0.01 m

0.075 maxisymmetricpile (rigid body)

grid spacing: 0.005 m
3m

number of particles per cell: 4 points

0.4 m
図 22: 杭貫入解析でのモデルと境界条件

Exp. STR02
Exp. STR04
Calc. 0.1m/s
Calc. 1m/s

Vertical load (kN)

60

40

C. Liu, Q. Sun, and G. G. Zhou. Coupling of material
point method and discrete element method for granular flows impacting simulations. International Journal
for Numerical Methods in Engineering, 115(2):172–188,
2018.
W. K. Liu, S. Jun, and Y. F. Zhang. Reproducing kernel

20

particle methods. International journal for numericalmethods in ﬂuids, 20(8-9):1081–1106, 1995.
0

0

500

1,000

1,500

J. Ma, D. Wang, and M. Randolph. A new contact al-

2,000

gorithm in the material point method for geotechnicalDisplacement (mm)simulations. International Journal for Numerical and

図 23: 杭の変位–荷重関係

Analytical Methods in Geomechanics, 38(11):1197–1210,
2014.
Z. Ma, X. Zhang, and P. Huang. An object-oriented mpm

C. Jiang, and J. Teran. Drucker-prager elastoplastic-

framework for simulation of large deformation and con-ity for sand animation. ACM Transactions on Graphics

tact of numerous grains. Computer Modeling in Engi-(TOG), 35(4):1–12, 2016.

neering and Sciences (CMES), 55(1):61, 2010.

V. Lemiale, J. Nairn, and A. Hurmane. Material point

J. J. Monaghan and J. C. Lattanzio. A refined particle
method for astrophysical problems. Astronomy and as-

method simulation of equal channel angular pressing involving large plastic strain and contact through sharp

trophysics, 149:135–143, 1985.

corners. 2010.

Y. Mukherjee and S. Mukherjee. On boundary conditions

Y. Lian, X. Zhang, and Y. Liu. Coupling of finite elementin the element-free galerkin method. Computational Me-

method with material point method by local multi-meshchanics, 19(4):264–270, 1997.

contact method. Computer Methods in Applied Mechan-

J. Nairn. Modeling imperfect interfaces in the material

ics and Engineering, 200(47-48):3482–3494, 2011.

- 22 -

MPM–剛体シミュレーションのための最小二乗法を用いた摩擦接触アルゴリズムの開発と開端杭の地盤への貫入挙動への適用

図 24: 杭貫入時の平均応力分布

D. Shepard. A two-dimensional interpolation function forpoint method using multimaterial methods. Comput.

irregularly-spaced data. In Proceedings of the 1968 23rdModel. Eng. Sci, 1(1):1–15, 2013.

ACM national conference, pages 517–524, 1968.

K. Nakamura, S. Matsumura, and T. Mizutani. Particle-

J. C. Simo and T. Laursen. An augmented lagrangian

to-surface frictional contact algorithm for material
point method using weighted least squares.
puters and Geotechnics,

134:104069,

treatment of contact problems involving friction. Com-

Com-

2021.

puters & Structures, 42(1):97–116, 1992.

doi:

M. Steffen, R. M. Kirby, and M. Berzins. Analysis andhttps://doi.org/10.1016/j.compgeo.2021.104069.

reduction of quadrature errors in the material point

S. Peng, C. W. W. Ng, and G. Zheng. The dilatant behaviour of sand–pile interface subjected to loading andmethod (mpm).

stress relief. Acta Geotechnica, 9(3):425–437, 2014.

methods in engineering, 76(6):922–948, 2008.

International journal for numericalA

D. Sulsky, Z. Chen, and H. L. Schreyer. A particle method

convected particle domain interpolation technique tofor history-dependent materials. Computer methods inextend applicability of the material point method for

applied mechanics and engineering, 118(1-2):179–196,

A. Sadeghirad, R. M. Brannon, and J. Burghardt.

1994.

problems involving massive deformations. International

D. Sulsky, S.-J. Zhou, and H. L. Schreyer. Application of

Journal for numerical methods in Engineering, 86(12):a particle-in-cell method to solid mechanics. Computer

1435–1456, 2011.

- 23 -

中村

圭太・松村聡・水谷 崇亮

図 25: 杭貫入時のせん断応力分布

physics communications, 87(1-2):236–252, 1995.
Q.-A. Tran, E. Wobbes, W. Solowski, M. Möller, andC. Vuik. Moving least squares reconstruction for b-spline
material point method. 2019.
J. Wang and D. Chan. Frictional contact algorithms in sph
for the simulation of soil–structure interaction. International Journal for Numerical and Analytical Methods inGeomechanics, 38(7):747–770, 2014.
中村圭太，元水佑介，松村聡，水谷崇亮，新谷聡，大下英治，
末政直晃. 開端杭の引抜き抵抗力に及ぼす施工過程と杭形状の影響に関する模型実験. 港湾空港技術研究所資料,

No.1364, 2019.
菊池喜昭，佐藤宇紘，森川嘉之. 開端杭貫入時の 杭先端周辺
地盤の変形挙動の観察. 港湾空港技術研究所資料, No.1177,

2008.

- 24 -

港湾空港技術研究所報告第 60 巻 第 1 号

２０２１．６

編集兼発行人

国立研究開発法人海上・港湾・航空技術研究所

発

港 湾 空 港 技 術 研 究 所

行

所

横 須 賀 市 長 瀬 ３ 丁 目 １ 番 １ 号
TEL. 046(844)5040 URL. http://www.pari.go.jp/

C （2021）by MPAT
Copyright ○All rights reserved. No part of this book must be reproduced by any means without the written
permission of the President of MPAT
この資料は、海上・港湾・航空技術研究所理事長の承認を得て刊行したものである。したがって、
本報告書の全部または一部の転載、複写は海上・港湾・航空技術研究所理事長の文書による承認を得ずしてこれを行ってはならない。

CONTENTS
1. Cyclic Beach Morphological Changes in the Swash Zone

due to Tidal Range Fluctuations
Masayuki BANNO・Yoshiaki KURIYAMA
2. Ocean swells induced by the Typhoon Faxai hit Port of Yokohama in 2019

Hitoshi TAMURA・Koji KAWAGUCHI・Takumu IWAMOTO・Takashi FUJIKI
3. Dynamics of ocean swells in Toyama Bay

Hitoshi TAMURA・Koji KAWAGUCHI・Takashi FUJIKI
4. Study on Stability of Temporary Armor Units against Waves
Kojiro SUZUKI・Hiroki KUBOTA・Tsutomu TANAKA
5. Study on Wave Pressure acting on Receding Parapet type Seawall
Kojiro SUZUKI・Hiroki KUBOTA・Naoki TSURUTA
6. Study of slow slip mechanism on plate boundary fault decollement zone
Yuri SUGIYAMA・Shinya TACHIBANA・Yoshiyuki MORIKAWA
7. Development of frictional contact algorithm using weighted least squares for MPM–rigid bodysimulation and its application to behavior of open-ended pile driven into the ground
Keita NAKAMURA・Satoshi MATSUMURA・Takaaki MIZUTANI