平成 22 年度 建設工学専攻修士論文用要旨
平成 23 年 3 月

MPS 法に基づく土構造物
づく土構造物の土構造物の弾塑性変形解析に
弾塑性変形解析に関する研究する研究

1. はじめに

環境防災研究室

八木 啓介

指導教官

大塚 悟，磯部 公一

第二項も同様に発散モデルを用いて次のように表される．

本研究の目的は，地震および水害による地盤構造物の崩壊メカニズムを明らかにし，限界状態・残留変位量をdiv(2 µ [ε ]ij ) = 2 µ適切に解析できる手法を構築し，地盤構造物の性能規定[ε ]ij rij .......................................... (6)
2dwij
0 ∑
n j ≠ i rij rij0

ここで，d は次元数，n0 は粒子数密度，w は重み関数で化の一助とすることにある．
そこで本研究では，大変形問題や複雑界面追跡に優れた解析手法として粒子法のひとつである MPS 法に着目ある．重み関数は，粒子間距離の関数であり，以下の式で表される．

し，地盤変形解析への適用性を検討するために第一段階として，Dracker-Prager 型降伏関数を用いた弾塑性動的解析コードを構築し，その妥当性について検討した．

( )

w rij

2. MPS 法に基づく弾塑性動的解析
2.1 弾性モデル1)

MPS 法の弾性動的解析手法に関しては越塚ら によっ
=


( 0 < r < r ) ...................................... (7)re−1rijije(r ≤ r )

0eij

( )

ni0 = ∑ w0 rij0 ...................................................................... (8)j ≠iて定式化されている.本研究では,FEM で一般的に用いre は影響半径という計算パラメータでありこれにより，られる弾塑性解析手法を導入し,MPS 法における弾塑性着目粒子 i から影響半径内に位置する近傍粒子を決定し，動的解析コードを構築した．以下にその計算方法につい重み関数を計算する（図-1）．本解析では，re=2.9l0 と設て示す．弾塑性体の粒子 i における運動程式および粒子定した．粒子数密度は，初期の粒子配置における着目粒

間 ij の応力の式は以下のように書ける．

子 i の全ての近傍粒子との重み関数の総和である．

Drɺ ∂σρi i = i + K i = divσ i + K i .........................................(1)Dt∂x○ :着目粒子 i[σ ]ij = λ (tr [ε ]ij )I + 2µ [ε]ij ................................................(2)rew (|rij|)○ :近傍粒子 jre :影響半径

ここで， rɺ ：速度ベクトル，ρ ：密度，σσ ：応力テンソ
ル，εε：ひずみテンソル，K ：外力ベクトル，λ，µ：ラl0 :初期の粒子間配置距離l0メの定数である．指標 i は粒子 i において，指標 ij は粒子 i と粒子 j の粒子間 ij において，という意味である．
粒子間の全ひずみテンソルは粒子座標位置から計算され

図-1 粒子間相互作用

る粒子間距離ベクトル rij と剛体回転成分を除去した粒

粒子間にせん断力が作用するとトルクが発生するので子間変位ベクトル uij のテンソル積から計算する．なお，角運動量保存則が保存されない．角運動量を保存させる

全ひずみテンソルは対称テンソルで定義する．

ためには，このトルクを打ち消すように付加的にトルク[ε ] ij =1 (u ij ⊗ rij + rij ⊗ u ij )2
rij rij0

uijx rijx
rij rij0= 1  u xr y + u yr x ij ij  ij ij
 2  rij rij0......(3)x yy x1  uij rij + uij rij  0εxx
2  rij rij  = ε yy ε uijy rijy
 xy ij0rij rij
 ij

粒子 i における応力の発散は MPS 法の加重平均，発散モデルを用いて以下の式で考える．(2)式第一項は以下の式から求められる.

λ tr [ε ]ij = λ(d∑ (ε xx + ε yy + ε zz )ij wij = pi ...........................(4)n0 j ≠i)div λ (tr [ε ]ij ) =(p + p )d∑ i j rij wij ...................................(5)n 0 j ≠ i rij rij0を与えなければならない．

1 ∂ω I i  i  = − ∑ Tij ............................................................... (9)2 j ≠i ∂t 

Tij = −rij × Fij ......................................................................... (10)s ∂rɺij 
mi 2d t ij
Fij = mi =wij ............................................. (11)00 ∂t  s ρi n rij

mi l02 ............................................................................... (12)6
ここに，tsij：せん断力， I：慣性モーメント，Tij:粒子間Ii =で発生したトルク，m:粒子1個の質量である．トルクTijは粒子i,jの各々に1/2の割合で配分される．
以上を用いて，verlet法による時間積分を行い位置座標，速度，角速度および回転角を更新する．

平成 22 年度 建設工学専攻修士論文用要旨
平成 23 年 3 月

rik +1 = (2rik − rik −1 ) +Drɺ 2 ...................................................(13)∆tDt

よる土塊の崩落をシミュレーションした．
0.75 m0.75(r k +1 − rik ) .................................................................(14)rɺik +1 = i∆t∂ω ...............................................................(15)
=ω +
∆t
∂t

k +1i

=θ +ωωθ

ki

k +1i

ki

幅：0.75 m
高さ：0.8 m

0.8m
0.8

k +1i

【土粒子（青色）】

突出し長さ：0.25 m∆t ...............................................................(16)【壁粒子（緑色）】
完全拘束条件ここに k：解析ステップ， ∆t :計算時間間隔(=1.0×10-5sec)である．さらに，本解析では吉田ら2)の手法に基づい0.25 m0.25

図-2 解析粒子モデルて減衰項の導入を行った．

表-1 解析に用いた物性値

2.2 弾塑性モデル
本解析では，降伏関数に Dracker-Prager 型を用いる．
Dracker-Prager 型の降伏関数 f は応力の第一不変量 I1および偏差応力の第二不変量 J2 によって表される関数である．(17)式に引っ張りを正と定義した Dracker-Prager型の降伏関数 f を示す．

f (σ ) = J 2 + α I1 − κ ....................................................(17)
ここで，α,κ は土の強度定数 c,φ に関係付けられる係数

弾性係数 E [kN/m2]
ポアソン比 ν粘着力 c [kN/m2]
内部摩擦角 φ [deg]粒子の密度 ρ [kg/m3]
減衰比 h5000
0.35
5.05.016600.01

解析結果より得たモデルの経時変化を図-3に示す．各部分で土塊を維持したままの分裂や，衝突が再現できた．

である．塑性ひずみ増分は流れ則より(18)式から求められる．

{d ε } = d λ  ∂∂σg  .............................................................(18)
pdλ は一般的な FEM で用いられる(19)の式から与えられる．
Te ∂f   [ D ] {dε} ....................................................(19)σ∂d λ =  Te  ∂g  ∂f   [ D]   ∂σ  ∂σ [0.5 sec][1.0 sec][1.3 sec][1.6 sec]

これより塑性ひずみ増分および塑性ひずみが求められる．
弾性ひずみは全ひずみと塑性ひずみと差から求める．

[ε ] = [ε ] − [ε ] ...............................................................(20)
eij

pijij

弾塑性解析では，ひずみに弾性ひずみを用いて前に示した弾性解析の計算を行う．
3. 大変形問題に対する解析事例
MPS法は，解析対象に粒子を用いている．これより

図-3 土塊の崩落状況
4. まとめDEMとほぼ同様の計算方法から離散体同士の衝突のよMPS法の弾塑性動的解析コードを構築し，地盤材料にうな現象も計算することができる．土塊が分裂するよう対するシミュレーションすることでMPS法の地盤変形解な大変形問題では，初期の粒子配置において近傍粒子で析への適用性について検討した．MPS法では大変形や分なかった粒子，また近傍粒子であっても粒子間の応力状裂といった現象を扱えることから，これまで解析するこ態がDrackerPrager型降伏関数の尖角に達したとき粒子とのできなかった事象に対しても有効な可能性を秘めて間において連続体の関係から離散体へと切断処理を行う．

いる．今後も，検討および発展に期待できる．

切断処理は，粒子間の重み関数を0とすることで粒子間の＜参考文献＞相互作用をないものにする．これら離散体の関係にある

1) 越塚誠一：計算レクチャーシリーズ 5 粒子法 ， 丸粒子間が接近した際にDEMの計算手法に基づき，衝突応答の計算を行い，座標位置を修正する．
大変形の計算事例に用いた解析粒子モデルを図-2に，解析に用いた物性値を表-1に示す．本解析では，自重に善 ，2005．
2) 吉田ら：MPS 法を用いた地震応答解析のための基礎検討， 土木学会論文集 A，Vol. 66 .No.2，pp. 206-218，
2010．