地盤材料の分岐挙動における塑性流動則の関連性の影響と分岐解析手法に関する研究
建設設計研究室
指導教官

1.

はじめに

地盤実験で地盤材料の破壊過程において材料の変形や局所化は分岐に起因し，分岐の発生は材料の構成式の塑性流動則による影響が大きい．一般に既存の分岐の研究では，
関連型の塑性流動則を用いているが，地盤材料のマクロ的力学挙動の精度に限界があり，より実モデルに近い非関連型を用いる必要性が問われている. だが非関連型を用いた
場合，どの程度分岐挙動に影響を及ぼすかは明確ではない．
そこで塑性流動則の違いが分岐挙動に対してどのような影響を及ぼすのかを調べることを本研究の目的とした．しかし，分岐解析手法が確立されていない接線剛性行列が非対称になりうることがあるため，従来の完全系分岐解析を行うことが出来ず，代替的手法として供試体に初期不整を与えて分岐挙動を近似的に評価する擬似分岐解析を用いた．
この解析結果より，荷重-変位関係や変形形態および局所的な視点で物理的な性質を比較した．また，最大荷重と初期不整の関係より完全系分岐荷重を比較した．その結果，塑
性流動則の相違は分岐挙動に対し影響を及ぼし，分岐挙動
は非関連型のほうが早い変形段階で分岐を生じ，非関連型のほうが局所化の進展が顕著で 分岐も早い段階で生じた．

2.

2.1 有限変形・弾塑性境界値問題
有限変形を考慮した一般的な弾塑性境界値問題の定式化についてまとめる．本研究では，運動力学的に厳密とされる変形勾配の乗算型分解 F = F e · F p と超弾性構成式に基づく定式化を採用する．滑らかな境界 ∂B を有する閉集合の領域 B ⊂ Rndim (ndim = 1, 2, 3) を考える．B は物体ないの参照配置を表し， ndim は領域の次元， R は実数の空間である．このとき，空間表示 (現配置参照型) の弱形式のつり合い方程式は，
dv
gϕt (τ ; η) =
τ : (η ⊗ ∇xt )
J
B


dv
T · η ds
− ρ0 B · η
−
J
B
∂t B
∀ η ∈ ν ϕt(2)

ここで，∂ϕ B は変位境界，Wは一階微分までが L (B)に属する Sobolev 空間である (p > 2)．つり合い式 (1) は1,p+ gϕt (τ t ; η) = 0∀ η ∈ ν ϕt(3)

左辺第 2 項の gϕt (τ t ; η) = 0 は不つり合い力 (残差力) に相当し，gϕt = 0 となるつり合い状態の解を求めることとなる．ここで，L∆u τ t は Kirchhoﬀ 応力の ∆u に関する
Lie 微分である．増分型の弾塑性構成式をL∆u τ t = cept : [∆u ⊗ ∇xt ](4)として，式 (3) に代入し，有限要素離散化を行い，増分ス
テップ毎に反復的に解く．ここで，cept は増分型の 4 階のconsistent 接線係数テンソルである．

2.2 弾塑性有限ひずみ理論と応力評価アルゴリズム(1) 超弾性構成式と塑性発展方程式
弾塑性材料の自由エネルギー関数 Φ は弾性ポテンシャル
W(be ) と塑性ひずみ硬化ポテンシャル H(ξ) を用いて，(5)として与えられるものとする．ここで，be = F e · F e,T は弾性の左 Cauchy-Green 変形テンソル，ξ は塑性ひずみ硬化に関するパラメータである．このとき超弾性構成式が，σ=22 ∂Ψ(be , ξ) e
2 ∂W(be ) e
·
b
=
2
·b
J
∂be
J ∂be

(6)

また，塑性ひずみ硬化則が，

q=−

∂H(ξ)
∂Ψ(be , ξ)
=−
∂ξ
∂ξ

(7)

と与えられる．ここで，σ は Cauchy 応力，q は塑性ひず
み硬化に関する応力のオーダーを持つ内部変数である．さらに，本研究では，弾塑性応答が最大塑性散逸の原理に従い，散逸不等式より，流れ則および発展則には関連型を用いるとそれぞれ，(1)と与えられる．ここで，ϕ(B) は変形写像，J = det F は
変形勾配 F = ϕ ⊗ ∇X のヤコビアン，τ は Kirchhoﬀ 応
力，ρ0 は質量密度，B は物体力，T は ∂t B 上で与えられる表面力，η は斉次境界条件を満たす可容変位変分の試験関数であり，その空間 νϕ は，νϕB = η : ϕ(B) → Rndim η ∈ [W 1,p (B)]ndim ,η(ϕ(X)) = 0 for X ∈ ∂ϕ B

変位に関して非線形であるため，有限の大きさを持つ増分
変位 ∆u を用いて線形化すると，dvL∆u τ t + (∆u ⊗∇xt ) · τ t : (η ⊗∇xt )JB

Ψ(be , ξ) = W(be ) + H(ξ)

有限変形弾塑性分岐解析の定式化= 0,

藪内 めぐみ鳥居 邦夫

1∂g(σ, q)− (Lv be ) · be,−1 = γ̇
2
∂σ

(8)

∂g(σ, q)ξ˙ = γ̇∂q(9)と導かれる．ここで，g(σ, q) は塑性ポテンシャル関数，γは塑性乗数，Lv be は be に関する Lie 微分である．これに，
f (σ, q) として Kuhn-Tucker 形式の負荷・除荷条件

γ̇ ≥ 0,f (σ, q) ≤ 0,γ̇f (σ, q) = 0(10)と塑性適応条件˙γ̇ f(σ,q) = 0pが課される．

(11)

(2) 降伏関数と塑性ポテンシャル関数
本研究では，静水圧依存型の弾塑性モデルである Drucker–
Pragar の降伏条件を用いて応力評価を行った．超弾性方程式 (5) について材料の等方性を仮定すると，弾性ポテン
シャル W は弾性の主ストレッチ λeA (A = 1, 2, 3) の関数であり，対数ひずみ εeA = log(λeA ) を用いた Henky モデル
W(εe1 , εe2 , εe3 ) =H0

H 0 / W0 = 2.00

2
1  e
λ ε1 + εe2 + εe32+µ (εe1 )2 + (εe2 )2 + (εe3 )2

y
x

W0

(12)

を用いる．ここで，λ > 0, µ > 0 は Lamé 定数である．これにより，Cauchy 応力の主応力が，

σA =

2
∂W
γ e
ε1 + εe2 + εe3 + µεeA ,
=
e
∂εA
2

(13)

また，Drucker–Prager の降伏関数は，

 

1 
f (σ, ξ) = √ ||dev[σ]|| + gf (σ̄) − σ̄y0 − q(ξ) ≤02(14)を用いる．ここで，σ̄ = 13 tr[σ] および dev[σ] = σ−σ̄1 はそれぞれ Kirchhoﬀ 応力 σ の等方成分，偏差成分，||dev[σ]|| =dev[σ] : dev[σ] は偏差応力のノルム，σ̂y0 は初期降伏応力である．また，g(σ̄) は静水圧の塑性への影響を表す関数
であり，gf (σ̄) = βf σ̄(15)という形式を用いることとすると，定数 βf は主応力空間における Drucker–Prager の降伏面による円錐の頂点に対応する．これらのパラメータ σ̄y0 , βf については Mohr–
Coulomb の粘着力 c̄ および摩擦角 θ̄ により，6c̄ cos θσ̄y0 = √,3(3 + A sin θ̄)

6 sin θ̄1
√ βf = √23(3 + A sin θ̄)(16)と決定できる (−1 ≤ A ≤ 1)．ここで，A は定数であり，
A = −1, 1 のときにはそれぞれ，Drucker–Prager 降伏面
は 3 次元応力空間において Mohr–Coulomb 降伏面の外側，内側の頂点を通る包絡線となる．また，塑性ポテンシャル
関数は，
1 dev[σ] + gg (σ̄) − [σ̄y0 − q(ξ)] ≤ 0 (17)g(σ, ξ) = √2

図–1 平面ひずみの供試体の有限要素解析モデル

3.

供試体の変形局所化挙動の擬似分岐解析手法

3.1 擬似分岐解析の概要
本研究では，供試体に初期不整を与えて分岐解析を行う
擬似分岐解析を行った．
この手法は供試体にあらかじめ初期不整を与えることにより，分岐解析で本来必要とされる手続きを経ずに近似的な分岐経路へと解を誘導し分岐挙動を擬似的に求める．そして，擬似分岐解析より得られた結果をもとに塑性流動則
による分岐に対する影響を考察する．
3.2 解析モデルと境界条件
本研究で用いた解析モデルを 図–1 に示す．平面ひずみ状態を仮定し，幅 W0 ，高さ H0 の矩形供試体の上端に一様な y 方向の変位が与えられるものとする．また，上下面の拘束条件については，剛体移動を止めるために下面中央部の１点のみを拘束して，それ以外は x 方向の変位を自由とした（端面摩擦なし）．また解析モデルの初期形状比は，標準的に用いられている H0 /W0 = 2.00 とする．
ここでの解析には，４節点アイソパラメトリック要素を用い，幅方向 48 要素，高さ方向 96 要素の 要素数 4,608(節点数 4,753) に均等分割したモデルを使用する．
3.3 材料定数と硬化関数
今回の解析で用いた硬化関数を示す．超弾性構成式 () 中の Lamé 定数 λ, µ は，ヤング率 E ，ポアソン比 ν とλ=νE,(1 + ν)(1 − 2ν)µ=E2(1 + ν)(19)

ここで，gg = βg σ̄ であり，βg はダイレンタンシー角 θ̃ に関して1
6 sin θ̃√ βg = √(18)
23(3 + A sin θ̃)と関係づけられる．ここで，ヤング率 E = 20 MPa，ポア
ソン比 ν = 0.3 とした．また，本研究ではひずみ軟化モデ
ルは用いないこととし，式 (),() 中の塑性ひずみ硬化関数
q(ξ) は，ξ n(20)q(ξ) = − τ̂y0 1 +
− τy0
ε̂y0

として表される．
実際の数値解析では，各荷重ステップの反復過程において数値積分点での応力評価を行うことになるが，これについては後退型 returen mapping 法を用いて厳密に降伏面上に位置する応力を求めるとともに，変形テンソルの弾性部
分 be をステップ毎に保存することで塑性変形履歴を更新する手続きを取る．

という n 乗硬化則の形式のひずみ硬化モデルを用いた．ただし実際の数値解析では硬化なしの完全塑性についてのみ検討を行った．ここで，τ̂y0 , ε̂y0 , n はそれぞれ降伏応力，
降伏ひずみ，塑性硬化則に関する材料パラメータであり，関連型・非関連ともに τ̂y0 = 100kPa, ε̂ y0 = 0.005, n =0.000, h = 0.00kPa とした．また，式 (??),() の定数は関連型が βf = 550kPa,, 非関連型は βg = 50kPa とした．

600

600

500

500

400

400

F / W0F / W0

a

300

δδ
δ
δ
δ

200
100

H0

0
0

5

300

100

10

15

(∆ H / H0 ) ×100

0
0

20

600

500

500

400

400

F / W0

F/ W0

図–2 初期不整の大きさ (反対称 1 次モードの例)

300

δ
δ
δ
δ

100
0
0

5

10

δ
δ
δ
δ

100

15

20

300
200

10

15

(b) 反対称 2 次

600

200

5

(∆ H / H0 ) ×100(a) 反対称 1 次

W0

δ
δ
δ
δ

200

0
0

20

5

10

15

(∆ H / H0 ) ×100(∆ H / H0 ) ×100(c) 対称 1 次(d) 対称 2 次

20

図–4 荷重-変位関係 ( 関連型 )
600

600

400

400

F/ W0

F/ W0

500

δ
δ
δ
δ

500

300

100

100

5

10

15

(∆ H / H0 ) ×100

0
0

20

(a) 反対称 1 次

δ=a× 100(%)
W0

本研究では初期不整 δ = 0.5%, 1%, 2%, 10% を与え，各々のケースについて関連型および非関連型の塑性流動則を用いて擬似分岐解析をおこなう．

4.

流れ則の相違や局所化におよぼす影響の検討

4.1 流れ則の相違による荷重−変位関係の相違について
前述の解析条件の下で得られた変位−荷重関係を初期不整のモードごとに示し，その相違を比較する．

20

500
400

F/ W0

F/ W0

(21)

15

(b) 反対称 2 次

400
300

δ
δ
δ
δ

200
100
0
0

10

600

500

3.4 擬似分岐解析
供試体の分岐モードでは sin θ,cos θ 波状拡散型 (数学的には調和型) のモードであることや，形状比 R =2.00 前後の供試体では多くの場合反対称 1 次モードの分岐が最小分
岐荷重を与えることが知られている．
そこで, 供試体に導入する初期不整のモードとして sin θ
形の不整を供試体側面に与え，対称 1 次，反対称 1 次，対称 2 次，反対称 2 次の不整形状を用いて検討を行う．
ここでは sin θ の半波長のモードを一次モードと称し，左右反対称性の有無によってそれぞれ対称モード反対称モードとする．図–3 に示すように，本研究では対称１次・対
称２次・反対称１次・反対称２次の四種類の変形モードを与える．また変形モードの大きさを初期不整といい，次式に大きさ δ の定義をしめす． W0 は供試体幅の長さ，a は与えたゆがみの長さである．(図–2 参照)

5(∆ H / H0 ) ×100600

図–3 変形形状 (変形モード) の種類

δ
δ
δ
δ

200

200

0
0

300

5

10

15

(∆ H / H0 ) ×100(c) 対称 1 次

300
200

δ
δ
δ
δ

100

20

0
0

5

10

15

20

(∆ H / H0 ) ×100(d) 対称 2 次

図–5 荷重-変位関係 ( 非関連型 )図–4 と図–5 より関連型・非関連型とも初期不整モード
によって，軸荷重に対する初期不整の大きさの影響の程度は異なるがモード毎のその大小については，関連型・非関連型を通じて同様の傾向がある．
図–6 をみると，全ケースを通じて関連型の場合よりも非関連型の場合のほうが軸荷重は小さい．軸荷重に及ぼす初期不整の大きさの影響の度合いは，初期不整モードによって様々である.
また，関連型では全ケースにおいて，変位に対する荷重の挙動が全体的に増加傾向にあるが，非関連型では一部それが見られない．これは材料モデルによる現象ではなく，数値解析的に非関連型のほうが変位の局所化の進展が大き
いため関連型のほうが荷重が有効に受け持っている領域が広くなり，非関連の方が早く荷重が低下すると予測される．

4.2 変形形状とせん断ひずみ分布
本節では，前節の解析結果をもとに発生するせん断ひず
みの分布や破壊形態について関連型と非関連型とひずみ局所化挙動の相違について比較を行った．図–7, 図–8 に供試

600

F / W0

400

初期形状

0
0

10

(∆ H / H0 ) ×100

関連型

非関連型

関連型

非関連型

関連型

非関連型

関連型

非関連型

初期不整 0.5%

200

20

初期形状

初期不整 1%

図–6 関連型と非関連型の挙動の比較 (反対称 1 次，初期不整 2%)

体の変形形状を示す．各図の右側に示してあるコンターレ
ベルは対数せん断ひずみの大きさである． 結果，関連型と
初期形状

初期形状

関連型

初期不整 0.5%

非関連型

初期不整 2%

初期形状

初期不整 10%

図–8 せん断ひずみの進展の比較 (例 対称 1 次)
初期形状

初期形状

関連型

非関連型

関連型

非関連型

初期不整 1%

初期不整 2%関連型

非関連型

初期不整 2%
図–9 せん断ひずみの分布 反対称 1 次 (初期不整 2%)

初期形状

関連型

初期不整 10%非関連型

図–7 せん断ひずみの進展の比較 (例 反対称 1 次)B*

非関連型とを比較すると，せん断ひずみのレベルはどの初期不整モードに形状においても関連型のほうが大きいが，ひずみの集中の度合いは非関連型のほうが局所化は顕著で
ある．

4.3

A*

関連型

除荷域の分布

本節では，実際に供試体のどの部分でどのような局所化が進行しているのかを明らかにするために局所的な視点での局所化の進展について検証する.
はじめに，例として反対称 1 次，初期不整 2% の場合を取り上げ，関連型および非関連型それぞれについて除荷・負荷の分布を 図–10 に示す.
同図はいずれもそれぞれ軸ひずみ 15%時のときの分布であり，黒い部分は除荷であり白い部分が負荷している部分となっている． 同図より関連型の方が，除荷と負荷の境界が明確であり，また非関連型に比べて幅が太いことがわ
かる．

D
C* *

非関連型

図–10 各要素の位置

4.4

局所的な挙動の分布

次にせん断ひずみ・偏差応力・等方応力や体積ひずみが供試体全体での変形である軸ひずみの進展に伴ってどのように変化していくかを調べ，その供試体内部の局所的な状態を調べる． 図–11 の応力経路より，各点における弾性・
塑性変形の履歴をみることが出来る．比較した結果，非関連型の方が中心付近 (点 B,C) で挙動が安定しており，偏差応力の減少が著しい．このため，非関連型の方が関連型よりも塑性変形の進展は著しいと考えられる．

200

300

300

200

100

100

p

p

q

q

200

100
100

200

0

0

200

p

200

点Ｂ

0

q
100

0

0

300

300

100

100

100

200

0

200

200

0

p

εv

0.6

0.1

点Ｃ

点Ｃ

0.6

300

100

100

εv

0.6

p

300

p

図–11 応力経路 (反対称 1 次，初期不整 2%)

100

100

0.1

εv

0.6

0.1

点Ａ

200

点Ｂ
300

q

q

εv

点Ｄ

関連型

非関連型

200

0.6

100

0.1

200

p

点Ｂ

εv

点Ｂ

p

200

100

0.1

点Ａ

関連型

q

0.6

εv

点Ｃ

200

100

0.1

200

p

p

0

300

100

100

0.5

εs

0
0

1

0.5

100

p

p

100

0
0

1

εs

点Ｃ

100

点 Ｄ

100

関連型

0
0.1

εv

εv

点Ｃ
q

q

100

0.1

0
0

点Ｄ

非関連型

200

0.2

0.6

0.6

0.5

εs

1

0
0

点Ｃ

0.5

εs

1

図–13 等方応力-体積ひずみ関係 (反対称 1 次，初期不整 2%)

点 Ｄ

非関連型

図–12 偏差応力-偏差ひずみ関係 (反対称 1 次，初期不整 2%)

図–12 で関連型と非関連型を比較すると，同じ偏差応力に対する偏差ひずみの増加は，中心周辺 (点 C,D) では関連型の方が大きい．また前述より，非関連型のほうが偏差
応力の減少が大きいことから，非関連型の方が関連型よりも強度が軟弱であることが考えられる．
図–13 で関連型と非関連型を比較すると，体積ひずみは非関連型のほうが小さく 1/10 程度の体積ひずみしか生じていない．これより体積変化率は関連型のほうが大きいことがわかる．

4.5

流れ則による局所化の様子の違い

関連型と非関連型のせん断ひずみの局所化挙動の相違を明確にするため，バンド状のひずみ集中域を横切った時のせん断ひずみ分布をそれぞれ示す.
図–14 に示すように供試体の対角寸法を D とし，供試
体左上方向への距離を a とする．対角線上には，図に示す
位置に A,B,C,D の点をとり，各点でのせん断ひずみを調べた. 図–14 の横軸には無次元化した距離 d = a/D をとり，縦軸にせん断ひずみ (対数ひずみの偏差不変量) をとる．結果，図–15 からわかるように，関連型・非関連型ともに対角寸法 D の 30 % 前後の領域にせん断ひずみが集中していることがわかるが，非関連型のほうがせん断ひずみの大きさ自体が小さいものの局所化の度合いが大きい．

D
a

d=

a
D

d :
a :
D:

図–14 せん断ひずみ分布の図示方法

4.6 擬似分岐解析による分岐荷重の予測
本節では前述した擬似分岐解析の結果をもとに各初期不整モードに対する分岐モードの発生荷重を線形近似が成り立つと仮定した上で外挿を行い, 完全分岐荷重を予測し，関連型と非関連型で比較する．方法は，初期不整における最大荷重をプロットし，そこから外挿により初期不整が 0.0% 時の分岐荷重を求める．分岐極大点が現れないケースについては土質試験による一般的方法と同様に軸ひずみが15%時のときの荷重を最大荷重と定義する．
また，佐藤5) により求められている関連型材料における分岐荷重と，ここで求める予測値を比較し，擬似分岐解析に基づいた予測値の精度についても検討した．各初期不整モードについて，初期不整の大きさを δ(%) と, 最大荷重
Fmax /W0 (kP a) との関係を関連型・非関連型それぞれの場合について図–16 に示す．横軸に初期不整を，縦軸に各初期不整ごとに最大荷重を示している．

小さくなる傾向がある．
2. 初期不整に誘発されて生じたせん断帯状のひずみ局
所化領域は，関連型の方がその幅が大きい．しかし非関連型のほうが局所化の進展が顕著である．
3. 擬似分岐解析より，分岐荷重の予測値は関連型のほうが 20%程度大きく, 非関連型のほうが低い荷重で分岐を生じる．

1

0.5

0
0

0.5

参考文献

1

図–15 関連型と非関連型のせん断ひずみの分布
500

500

.

400

400

0.9827

Pmax

Pmax

Pmax / W0 = −12.144δ + 398.33.

300

0.9992
Pmax /W
W0 = −6.272δ + 323.41

300

0

4

8

12

0

δ

関連型

4

δ

8

12

非関連型

(例) 反対称１次
図–16 擬似分岐解析による初期不整の予測

表–1 分岐荷重の予想

表–1 からわかるように，すべてのケースにおいて分岐荷重の予測値は，非関連型の方が関連型よりも小さい．これは非関連型の剛性行列は非対称になる場合があり，関連
型よりもゼロ固有値は早く発生する傾向があるため非関連型の方が最大分岐荷重が小さくなると考えられる．
また，関連型の場合について完全系の分岐解析により求められた分岐荷重との比較をすると 20% の誤差で比較的精度よく求められている．
関連型では線形近似式との適合がよいほど正確な分岐荷
重が擬似分岐解析により予測可能で関連型・非関連型ともに初期不整が小さい場合は最大荷重と初期不整との関係は線形関係にあるといえる．

5.

結論
1. 初期不整モードによらず，初期不整が大きいほど軸荷重の挙動は小さくなる. また，初期不整に対する軸荷重の感度は初期不整モードごとに異なり，全モー
ドを通じ，関連型よりも非関連型のほうが軸荷重が1) 山川 優樹, 池田 清宏, 須藤 良清, 寺井 信夫, 鳥居 邦夫: 砂の三軸試験供試体の分岐変形モードと寸法・形状効果, 土木学会論文集, No. 701/III-58, pp. 357–371 , 2002.
2) 山川 優樹, 寺田 賢二郎, 池田 清宏, 鳥居 邦夫: 圧縮場における弾塑性体の分岐解析とパスジャンプ挙動, 土木学会論文集, No. 701/III-58, pp. 73–86 , 2002.
3) 藤井 文夫, 野口 裕久: 固有値解析を必要としない大規模非線形構造系の分岐座屈モードの求め方, 構造工学論文集, Vol.
46A, pp. 241–250, 2000.
4) 野口 裕久, 久田 俊明: Scaled Corrector を用いた有限要素
分岐解析手法の開発, 日本機械学会論文集 (A 編), Vol. 58,
No. 555, pp. 2191–2198, 1992.
5) 佐藤啓介，山川 優樹, 池田 清宏, 鳥居 邦夫: 分岐解析によるひずみ局所か現象における供試体形状と境界条件の影響に関する検討, 応用力学論文集, Vol. 5, pp. 329–338, 2002.