有限変形理論に基づく土構造物の動的挙動解析
地盤防災研究室 山田健士
指導教官 大塚 悟

3. Newmark-β法
時間積分には Newmark-β法を採用した．この手法
では時刻 t + ∆t における速度，加速度を次のように仮定する．
 γ γγ k∆ u ( ) +  1 −  ut − ∆ t  1 −ut +∆t =
 ut
β ∆t
 β
 2β ut +∆t =(3.1)

 1
11kut − ∆u ( ) −− 1  ut
2β ∆tβ ∆t
 2β時刻 t + ∆t における運動方程式を数値解析するが，方程式は変位とともに変化するために，運動方程式は収束計算が必要になる．収束過程において k 回目と k-1k回目の差分量 δ u ( ) (修正子)を用いると運動方程式は次のように得られる．

M(t +∆t )ut(+∆) t + C(t +∆t)ut(+∆) t + K (t +∆t)δ u ( )k −1k −1k

= ft +∆t − Qt(+∆t )k −1kk −1k(3.2)

ここに M ， C は質量，減衰マトリクス，f ，Q は節点力，内力ベクトルである．質量マトリクスにはconsistent 質量マトリクスを用いている．式(3.1)より運 1γk −1k −1k −1 kM (t + ∆t ) +C (t + ∆t ) + K (t + ∆t )  δ u ( )2∆∆ttββ
(3.3)= f t + ∆t − Qt(+ ∆t ) − M (t + ∆t )ut(+ ∆t ) − C (t + ∆t )ut(+ ∆t )
k −1

k −1

k −1

k −1k −1

上式の右辺を残差ベクトルと定義し，これが 0 となるまで繰り返し計算を行う．
4. 等価線形解析
本解析では次に示すせん断ひずみ γ を使用した．
21(4.1)ε x − ε y ) + γ xy2(2さらに，地震加速度の入力によって地盤内部に発生するせん断ひずみの最大値 γ max に対して，有効せん断ひずみを次式で仮定する．
γ eff = αγ max(4.2)
ここに，係数α は任意定数であるが，本解析では一般に用いられている α = 2 3 とした．せん断剛性係数 Gと減衰定数 h のひずみ依存性は図-1 中のγ eff を変数とした G − γ 関係，h − γ 関係によって表されるものとす
る．G0 は初期せん断弾性係数である．

γ=10.30.80.20.610.4 G =G0 1 + 50 ⋅ γ 0.60.2

h=0.25 ⋅ γ+ 0.02
0.005 + γ
0.1

00.00001 0.0001 0.001 0.01
shear strain γ

0.1

damping ratio h2. 有限変形理論
本解析では現在配置を常に基準配置と見なす
updated Lagrange 形式による定式化を行い，応力速度には客観性のある Truesdell 応力速度を採用した．
i(2.1)
σ = σ + σ trl − l ⋅ σ − σ ⋅ lT
なお，構成式は次を用いた．
(2.2)
σ = Ce : de
ここに，C は線形弾性マトリクスである．

動方程式は次式となる．

modulus ratio G/G0

1. はじめに新潟県中越地震では道路や宅地等の盛土構造物に多くの被害が見られた．阪神大震災以降，土構造物に対しても直下型地震のレベルⅡ地震動に対する設計法が確立されてきたが，残留変形量を精度良く求めるためには未だ多くの課題が残されている．本研究では地震時に変形する物体に作用する慣性力の影響を考慮することを目的に有限変形理論に基づく動的解析プログ
ラムを作成し，その妥当性について検討する．新潟県中越地震では余震による土構造物の被害の拡大が生じていることから，構造物の剛性低下の視点から変形量の進展について考察する．

01

図-1 せん断剛性と減衰定数のひずみ依存曲線
5. 1 次元解析による性能評価
まず，静的解析を行った．解析モデルは，図-2 に示す 1 次元要素とし，圧縮荷重を載荷した．土質定数は表-1 に示す．
図-3 に鉛直応力-変位関係の本解析による数値解と理論解を示す．両者はともに線形な関係を示すとともに完全に一致している．
次に動的解析を実施した．解析モデル，土質定数は静的解析と同様のものを使用し，100kN/m2 の一定荷重を載荷した．Rayleigh Damping の係数α，βはα= 0.4， 6. せん断波速度の検討
30.
β= 0 とした．このときの上端節点における鉛直方向
表-2 地盤定数27.の変位時刻歴を図-4 に示す．減衰を考慮しない場合，
E
1 MPa
要素は約 1.3 (m)の振幅で振動する．一方，減衰を考慮
24.
ν
0.3
した場合，時間の経過とともに振幅は著しく抑制され，21.
最終的には静的問題の解に収束する．減衰を考慮する
ρ
2.0 t / m3
18.
ことにより加速度 u が減少し，復元力のみによって釣図-5 に示す 30 (m)×0.5 (m)合った状態となる．以上より，1 次元動的問題におけ30.
30 15.の単純なビルディングモデル
る本解析の妥当性を確認できた．
12.
を用いてせん断波の伝播を検9.
討した．用いた地盤定数は
σ 226.
表-2 に示す．平面ひずみ条件
にて下端 xy 固定，側 y 方固定
3.
u
9.0 (m)とした．せん断波の伝播状況l0に着目するため，減衰は考慮
l0.5
単位: (m)しないものとする．入力波は1
固有振動数2Hz の規則せん断
図-5 ビルディング
11.0 (m)波 400sin4 π t (gal)とした．図-5モデル及びせん断に時刻 t =1.0 (s)におけるせん
図-2 1 次元モデル
波伝播状況
断波の伝播状況を示す．せん表-1 土質定数
断ひずみは下端から上部へと進行し，約 14 (m)まで伝E
1 MPa
播していることが分かる．弾性せん断波速度は理論式によって次のように求められる．
ν0.3

ρ

19.8 t m3

Vs = G ρ = 13.87 ( m s )両者を比較すると，本解析による結果は妥当性が得られる範囲であることを確認できる．

stress σ22 (kN/m2)

600
500

7. 斜面モデルによる動的挙動解析

4003007015P2200

P1数値解理論解
123displacement u (m)static0130

単位：( m )図-6 斜面モデル

dampingなし
damping-120254-0.510P415

図-3 垂直応力-変位関係-1.50P31010000displacement (m)(6.1)

30time (s)図-4 変位時刻歴における減衰効果

表-3 地盤定数

E5.0 MPaνρ

0.3

3.0 t / m3

始めに，線形弾性解析を行った．図-6 に斜面モデルを，表-3 に地盤定数を示す．Newmark-β法の係数γ，βはそれぞれγ= 0.5，β= 0.25 とする．入力加速度は400sin4 π t (gal)の規則波とする．平面ひずみ条件にて，下端 xy 固定，側方 x 固定の拘束条件とした．時間間隔はΔt =0.05 (s)とし，200 ステップ，計 10 (s)の解析を

20acceleration (m/s2)

行った．ここでは有限変形，微小変形についての解を比較するとともに，減衰の有無による動的挙動の変化について検討した．なお，減衰時の Rayleigh Dampingの 係数は 斜面 モデル の固 有値か らα = 0.0108 ，β= 0.0333 と決定した．
図-7，8 に水平変位最大時におけるせん断ひずみ分布を示す．非減衰の場合，斜面内のひずみは点々と分布し，せん断変形が数ヶ所で生じていることが分かる．

100-10-20000.0055

0.011

0.016520.022

input wave
46time (s)

FDID
810

図-11 斜面上端の加速度時刻歴 (非減衰)図-7 変位最大時のせん断ひずみ分布 (非減衰)00.00225

0.00450.00675

0.009

acceleration (m/s2)

840-4-802input wave
468time (s)

FDID
10

図-12 斜面上端の加速度時刻歴 (減衰)図-8 変位最大時のせん断ひずみ分布 (減衰)displacement (m)0.40.20-0.2-0.40246time (s)8FDID10図-9 斜面上端の水平変位時刻歴 (非減衰)displacement (m)0.40.20-0.2-0.40246time (s)

8FDID
10

図-10 斜面上端の変位時刻歴 (減衰)

一方，減衰を考慮した場合，斜面下部においてせん断ひずみが卓越し単純な変形形態を示している．なお，せん断ひずみの分布に有限変形と微小変形との明確な差異は現れなかった．
図-9，10 に減衰・非減衰時における斜面上端節点(P3)の変位時刻歴を示す．減衰の有無に関わらず有限変形(FD)と微小変形(ID)の解に差異が現れている．これは，せん断変形に伴う斜面要素のスピンの影響，変形した物体に作用する慣性力の影響であると考えられる．
図-11，12 に減衰・非減衰時における P3 の応答加速度時刻歴を示す．減衰を考慮しない場合，有限変形と微小変形で差異が発現した．有限変形理論では変形により要素の重心位置が刻々と変化するため，このような差異となって現れると考えられる．減衰を考慮すると，両者ともに加速度応答に増幅が見られず，入力加速度とほぼ同位相･同振幅となっている．
次に，等価線形解析を行った．解析モデル，拘束条
件，地盤定数は線形弾性解析の場合と同様のものを扱った．図-13 に斜面各点における水平変位時刻歴を示
す．高所に位置する点ほど著しい変位を示し，最も大きく変位した点は，斜面上端部であることが分かる．
さらに，図-14 に示した変位最大時の水平変位の分布を見ると，斜面上部は慣性力により高い変位を示し，そこを中心とした円弧状に水平変位が生じている．

acceleration (gal)

図-15 に最大変位時における斜面上部のせん断ひず より，地震の繰返しによる変形量，せん断ひずみの拡み分布を示す．せん断ひずみは斜面の底部に分布して 大を表現することができた．
いる．このとき，下部要素は斜面前方へと局所的なせ800
ん断変形を生じ，斜面中上部は剛体的に変位している．
0-4000-800
0-0.2-0.40246time (s)8P1P2P3P40.3

100.2

図-13 斜面上端の水平変位時刻歴

510152025time (s)図-16 入力加速度波形

displacement (m)displacement (m)

0.2

400

1波目
2波目

0.1
0-0.1-0.43-0.35-0.27-0.19-0.110(m)-0.2-0.30510152025time (s)図-17 1 波目，2 波目の斜面上端の水平変位時刻歴
図-14 変位最大時の水平変位分布00.00480.01320.0180.024

図-15 変位最大時のせん断ひずみ分布8. 結論
以下に，本研究で得られた知見を列挙する．
① 1 次元解析において，静的問題に対して本解析による数値解と理論解が完全に一致した．また，動的問題に対して減衰強制振動の収束値が静的問題の解に一致した．以上より 1 次元問題に対する本解析の妥当性が得られた．
② 斜面モデルによる線形弾性解析を実施し，有限変形理論と微小変形理論の解に差異があることを確認した．
③ 斜面モデルによる等価線形解析を実施し，実現象に近い応答変位が得られた．
④ 等価線形解析を応用し，本震・余震の繰返し載荷による動的挙動の変化について検討し，2 波目における水平変位，せん断ひずみの拡大を表現した．

最後に等価線形解析を応用し，本震･余震の繰返しによる動的挙動の変化について検討する．入力加速度は図-16 に示した，昨年 10 月 23 日に発生した新潟中越地震において長岡で観測された実波形データを用いる．最大加速度は 800 (gal)である．本震が伝播すると斜面内部のひずみ集中箇所では剛性が著しく低下する．
その剛性を余震時の初期せん断剛性として入力する． 【参考文献】
この解析では，斜面内の剛性低下による挙動変化に着 【1】 地盤工学会編：地震時の斜面の不安定化メカニズムと設計法に関する研究委員会 報告書，地盤工学会，
目しているため，本震と余震は同じ波形とした．
pp.117-144，1999
図-17 に 1 波目と 2 波目の斜面上端の応答水平変位時刻歴を示す．1 波目と 2 波目の水平変位に差異が現 【2】 社団法人日本塑性加工学会編：非線形有限要素法，
コロナ社，2002
れている．全体的に 2 波目の水平変位に拡大の傾向が見られる．最も大きな水平変位の差異は約 5 (cm)であ 【3】 Javier Bonet ， Richard D.Wood ： NONLINEARCONTINUUM MECHANICS FOR FINITE ELEMENTった．せん断ひずみ分布は斜面底部に集中し，本震･ANALYSIS，Cambridge University Press，1997
余震で変形形態に大きな差異はないものの，斜面内のせん断ひずみは．余震時では大きな値を示した．これ