応用力学論文集
13, pp.411-422
pp. 1–12（2010（2010年 8 月）
応用力学論文集 Vol.
Vol.13,年 8 月）

土木学会
土木学会

拡張下負荷面 Cam-clay モデルの有限変形理論に基づく定式化とリターンマッピングを用いた陰的応力更新法
Formulation and implicit stress-update algorithm of the extended subloading surface Cam-claymodel with kinematic hardening for ﬁnite strains山川優樹∗ , 山口洋介∗∗ , 橋口公一∗∗∗ , 池田清宏∗∗∗∗Yuki Yamakawa, Yosuke Yamaguchi, Koichi Hashiguchi, Kiyohiro Ikeda∗ 正会員 博（工） 東北大学准教授 大学院工学研究科土木工学専攻（〒980-8579 仙台市青葉区荒巻字青葉 6-6-06）∗∗ 学生会員 東北大学大学院工学研究科土木工学専攻（同上）∗∗∗ 正会員 工博・農博 第一工業大学教授 工学部社会環境工学科（〒899-4395 鹿児島県霧島市国分中央 1-10-2）∗∗∗∗ 正会員 Ph.D. 東北大学教授 大学院工学研究科土木工学専攻（〒980-8579 仙台市青葉区荒巻字青葉 6-6-06）

This paper presents a formulation and a stress-update algorithm for the extended subloading
surface Cam-clay model with kinematic hardening at ﬁnite strains. The model is re-formulatedon the basis of the framework of the multiplicative hyperelasto-plasticity. The hyperelastic law,yield function, and plastic evolution laws of the model are expressed in terms of tensorial variables deﬁned on the intermediate conﬁguration. The stress-update algorithm using the so-called
exponential map-based return mapping is developed and is implemented into a ﬁnite element
code with numerically-computed algorithmic tangent moduli. Numerical assessment conﬁrmed
accuracy and robustness of the proposed algorithm even for large incremental deformation.

Key Words : multiplicative ﬁnite-strain hyperelasto-plasticity, Cam-clay, extended subload-

ing surface model, anisotropy, kinematic hardening, exponential map-besed
return mapping

1.

は (1) が有する上掲の問題点を解決でき，応力計算の際にリターンマッピングを適用することにより，変形増分の大きさに依存しない高精度の解を得られるなどの数値解析上の利点を有している．

はじめに

弾塑性解析の分野では，Simo and Ortiz1) , Simo andTaylor2) , Ortiz and Simo3) らによる先駆的な研究により，数値解析を前提とした弾塑性理論の枠組みが体系化された．また，この枠組みのもとで，リターンマッピ
ング法を用いた高精度・高効率の陰的応力更新法が開発され，幅広く利用されるとともに，より高度な手法への改良や各種構成モデルへの適用が行われている4),5),6) ．
こうした，いわゆる計算塑性学の進展と並行して，有限変形における弾塑性理論に関する研究も盛んに行われている．現在のところ，有限変形弾塑性理論の枠組
みは大きく分けて (1) 変形速度の加算分解と亜弾性構成式を用いた定式化と，(2) 変形勾配の乗算分解と超弾性構成式を用いた定式化のふたつの系統がある．各種の商用有限要素解析コードではどちらの枠組みに基づく応力計算法も採用されており，とくに (1) は複雑な構成モデルにも対応しやすいなどの理由により，以前
から幅広く用いられている．しかし，この方法には以下に挙げるような問題点があることが知られている7) ．

1. 用いる客観応力速度によっては，明らかに不合理な挙動を示すことがある8) ．
2. 弾性変形でエネルギー散逸を生じる場合がある9) ．
3. 数値計算のために速度形構成式を増分表示に書き改める際には，客観性の保持に注意が必要である10) ．
一方，(2) の枠組みについては 1980 年代後半から急速に研究が進展し，近年主流となりつつある．この方法

近年では，地盤材料の構成モデルについても (2)の枠組みでの定式化と応力更新法が提案されている11),12),13) ．これらはいずれも現配置表示の主応力と対数主ひずみに基づく定式化14) を行っており，等方性材料に限定される．一方，異方性弾塑性モデルの定式化は現配置表示でも可能だが，中間配置を参照する諸量を用いる方が便利である．地盤材料以外では，中間配置表示で定式化を行っている研究15),16),17),18),19) もあるが，いずれも von Mises モデルなどの基本的な塑性モデルにとどまっている．地盤材料モデルを中間配置記述で定式化した研究20) もあるが，等方硬化だけを対象としており，異方性モデル（移動硬化などの誘導異方性も含む）を扱った研究はない．地盤材料の移動硬化モデルの陰的応力更新法21),22) も開発されているが，いずれも微小変形理論に基づくものである．しかしながら，地盤材料は塑性変形に伴う誘導異方性が強く発現する材料であり，そうした性質を表現できる構成モ
デルも提案されており23) ，有限変形理論に基づく定式化と高い適用性を有する応力更新法の開発が望まれる．
そこで本研究では，変形勾配の乗算分解に基づく有限変形理論の枠組みで，誘導異方性を表現しうる拡張
下負荷面 Cam-clay モデル24) の定式化を行う．その際，

1- 411 -

等方性モデルに限定した定式化ではなく，異方性弾塑性モデルにも適用しうる定式化とし，中間配置を参照
する諸量を用いて降伏関数や塑性発展方程式を記述する．また，地盤材料に特有の圧力依存性を有する体積
圧縮・せん断特性を表現しうる既存の超弾性モデルを拡張し，異方性モデルの定式化にも適用しやすい中間
配置でのテンソル表示の超弾性モデルを新たに提案する．さらに，このモデルについて，リターンマッピングを用いた陰的応力更新法を開発する．最後に，提案手法を数値解析例題に適用し，大きな変形増分を与えた場合の応力の計算精度を検証するとともに，境界値問
題の解析におけるロバスト性や解の収束性を検証する．

2.

拡張下負荷面 Cam-clay モデルの有限変形理論に基づく定式化

2.1

変形勾配の乗算分解に基づく弾塑性有限変形理論
における諸量の定義

本研究では，変形勾配の乗算分解25) に基づく弾塑性
有限変形理論に基づき定式化を進める．有限変形理論の枠組みは，等方弾塑性体について Simo26) が体系化しており，異方性弾塑性について文献27) の第 5.3 節にまとめられているので，ここでは詳細は割愛し，本論文の定式化で用いる応力や変形等の諸量の定義のみを示す．
以下ではとくに断りのない限り，原則として，基準配置を参照する諸量は大文字で，現配置を参照する諸量は小文字で，また，中間配置を参照する諸量は大文
字に上線を付して表記する．次式のように弾性部分・塑
性部分への変形勾配 F の乗算分解を仮定する．

F = Fe · Fp

1
1

ading

(1)

中間配置を参照する弾性右 Cauchy–Green 変形テン
ソル C̄ e および塑性速度勾配 L̄p を次のように定義する．

(2)

p −1

(3)

L̄ := Ḟ · (F )p

次に，中間配置を参照する応力は，現配置を参照する
応力を弾性変形勾配 F e を用いて中間配置へ引き戻すことにより，次のように与えられる．

S̄ := (F e )−1 · τ · (F e )−T ,e −TT̄ := (F ) · τ · (F )
e T

ding

図–1 ln v–ln p∗ 関係で表した等方圧密特性

ここで，τ は現配置を参照する Kirchhoﬀ 応力，S̄ は中間配置を参照する第 2 Piola–Kirchhoﬀ 応力，T̄ は中間配置を参照する Mandel 応力である．Mandel 応力 T̄ と
塑性速度勾配 L̄p は塑性散逸仕事について共役であり，これらは降伏関数や塑性流れ則の記述に用いる．なお，
Mandel 応力 T̄ は弾性モデルが等方性を有する場合を除き，一般には非対称テンソルである．

2.2 等方圧密特性の表現と等方硬化則
一般に Cam-clay モデルでは等方圧密履歴に基づいて硬化・軟化を表現する．そこで，はじめに等方圧密特性を表すモデルを示した上で，それに基づく等方硬化則を導出する．なお，本論文では一貫して，応力・変形は引張りを正と定義していることに注意されたい．

J = det(F e · F p ) = (det F e )(det F p ) = J e J p= C̄ · S̄(6)

ここで，弾性体積変化を J := det F , 塑性体積変化をJ p := det F p と定義した．式 (6) 両辺の対数を取り，対数体積ひずみ εv := ln J およびその弾性部分 εev := ln J eと塑性部分 εpv := ln J p を定義する．それらの時間微分˙ε̇ev = J˙e /J e , ε̇pv = J˙p /J p であはそれぞれ ε̇v = J/J,
る．以上より，対数体積ひずみの加算分解とその速度形が次のように導かれる．
eεv = εev + εpv ,31)e

ε̇v = ε̇ev + ε̇pv

(7)∗

ここでは Hashiguchi による体積 v と圧力 p の両対数関係を用いて等方圧密特性を表す．ここで p∗ (:=1
3 tr[τ ]) は Kirchhoﬀ 平均垂直応力（Kirchhoﬀ 圧力）である．圧密降伏応力を p∗c とすると，図–1 を参照して，
正規圧密時 (p∗ = p∗c ) はvp∗εv = ln= −λ̃∗ ln ∗c for p∗ = p∗c
(8)
v0
pc0
という関係がある．λ̃∗ は Kirchhoﬀ 圧力に対応する弾
塑性圧縮指数1 ，v0 は基準状態での圧密降伏応力 p∗c0 に

(4)
e

g

(1) 等方圧密特性の表現
式 (1) の乗算分解より，体積変化 J = det F は次のように弾性部分と塑性部分に乗算分解される．

塑性変形勾配 F は基準配置から中間配置への変形を，弾性変形勾配 F e は中間配置から現配置への変形をそれぞれ表す．一般には中間配置は剛体回転に対して任意
であり，塑性モデル諸式が満たすべき客観性について
多くの議論がある28),29) ．本論文では，材料の下部構造の方向を表す互いに直交する三つの単位方向ベクトルからなる方向三対は基準配置と中間配置で一致し，中間配置は等傾的であるとする18),19),30) ．このように中間配置を選択した場合，変形勾配 F への付加的な剛体回転に対して，中間配置を参照する諸量は不変である．

p

din

unloa

p

C̄ e := (F e )T · F e ,

loa

re-lo

1 この λ̃∗ は，微小変形理論で通常用いられる Cauchy 応力 σ(:= τ /J) の平均垂直応力 p (:= 13 tr[σ] = p∗ /J) に対応する(5)

弾塑性圧縮指数とは異なるため，上添字 (∗) を付けて示す．ま2- 412 -

対応する体積である．また，除荷・再載荷時 (p∗ > p∗c )においては，vc = v(p∗c ), vcp = v(p∗c0 ) とすると，vp∗∗p = −κ̃ ln ∗ ,
vc
pc0
p
v
p∗
εpv = ln c = −(λ̃∗ − κ̃∗ ) ln ∗cv0pc0
εev = ln(1) 超弾性ポテンシャル関数
弾性体積変化 J e = det F e = (det C̄ e )1/2 と弾性右Cauchy–Green 変形テンソル C̄ e から体積変形を除いたˆ e = 1) を用いて，
ˆ e (:= (J e )−2/3 C̄ e , det C̄
等積成分 C̄ˆ e ) = W̄ e (J e ) + Ŵ e (J e , C̄ˆ e)
W e (J e , C̄
(14)

(9)
(10)

という形式の超弾性ポテンシャル関数を考える．ここで， e −1/κ̃∗Jee∗
W̄ (J ) = −P̄0 κ̃
,
(15)
J0e
ˆ e ) = 1 μe (tr[C̄ˆ e ] − 3)(16)Ŵ e (J e , C̄
2である．J0e , P̄0 はそれぞれ基準状態での弾性体積変化と圧力である．また，次のような圧力依存型のせん断弾性係数 μe を考える．
αμe = μe0 + ∗ W̄ e (J e )(17)κ̃
ここで，μe0 はせん断弾性係数の定数部分，α はせん断弾性係数の圧力依存性を表す係数であり，μe は J e に依存して変化する．α = 0 とした場合には μe = μe0（一定）となる．

という関係がある．ここで，κ̃∗ は Kirchhoﬀ 圧力に対応する弾性圧縮指数（膨潤指数）である．式 (7) の第一式に式 (9), (10) を代入すると次式を得る．

εv = −(λ̃∗ − κ̃∗ ) ln
(2)

p∗c
p∗
− κ̃∗ ln ∗
∗
pc0
pc0

for p∗ > p∗c (11)等方硬化則

体積が v0 である基準状態での体積変化を J = 1 とすると，vcp /v0 = J p であるから，式 (10) より次式の等方硬化則を得る．

p∗c = p∗c0 (J p )−Θ(12)

ここでは Θ := 1/(λ̃∗ − κ̃∗ ) と定義した．等方弾性の場合
は Kirchhoﬀ 圧力 p∗ と Mandel 平均垂直応力（Mandel
圧力）P̄ (:= 13 tr[T̄ ]) は一致するので，等方硬化則 (12)は次のように中間配置を参照する形式に書き改めることができる．

P̄c = P̄c0 (J p )−Θ
2.3

(13)

圧力依存型の超弾性構成モデル

一般に Cam-clay モデルで用いられる弾性モデルでは，ln v–ln p 関係における膨潤線の挙動と整合するように，圧力依存型の体積弾性係数が用いられる．せん断については，せん断弾性係数あるいは Poisson 比を一定とする場合もあるが，せん断弾性係数の圧力依存性を考慮に入れる方が実験挙動を良好に表現できることが知られている．Houlsby33) は，圧力依存型のせん断弾性係数を用い，なおかつ弾性体本来のエネルギー保存を保証する弾性モデルを微小変形の枠組みで提案した．このモデルは Borja et al.34) や Tamagnini et al.35) により改良が加えられ，その力学的・数学的性質は Niemunis
and Cudny36) により検討されている．さらに，Borja
and Tamagnini11) や Callari et al.32) は対数ひずみを用いてこのモデルを有限変形に拡張した．しかし，このモデルは現配置表示で対数主ひずみと Kirchhoﬀ 主応力との関係として表されているため，本研究で行う中間
配置表示に基づく弾塑性モデルの定式化には適用しにくい．また，この形式では等方性モデルの定式化に限定され，異方性モデルに拡張できない．そこで本節では，上述の弾性モデルと同等の性質を有し，なおかつ異方性弾塑性モデルの定式化にも適用しやすい中間配
置でのテンソル表示の超弾性モデルを新たに提案する．
た，よく用いられる e–ln p 関係（e は間隙比）における弾塑性
圧縮指数とも異なる．以上の二点は文献11),12),32) で議論されており，後述の弾性圧縮指数 κ̃∗ も同様である．

(2) 超弾性構成式
式 (14) を 12 C̄ e で偏微分すると，中間配置を参照する第 2 Piola–Kirchhoﬀ 応力が次のように得られる．
ˆ e)∂W e (J e , C̄S̄ = 2∂ C̄ e e −1/κ̃∗ Jαˆ e ] − 3) (C̄ e )−1
C̄
1
+
= P̄0
(tr[
J0e
2κ̃∗

1
(18)
+ (J e )−2/3 μe I − tr[C̄ e ](C̄ e )−13
ここで I は単位テンソルである．基準状態 J e = J0e ,ˆ e = I において S̄ = (J e )−2/3 P̄ I とC̄ e = (J0e )2/3 I, C̄
00
なり，Mandel 平均垂直応力は P̄ = 13 tr[C̄ e · S̄] = P̄0となる．
(3) 超弾性接線係数
超弾性接線係数は，式 (14) を 12 C̄ e で二階偏微分することにより，中間配置表示で次のように得られる．
ˆ e)∂ 2 W e (J e , C̄eC̄ = 4∂ C̄ e ∂ C̄ e−1/κ̃∗
P̄0 J e
=− ∗
2κ̃
J0e


α
ˆ e ] − 3) (C̄ e )−1 ⊗ (C̄ e )−1
× 1 + ∗ (tr[C̄2κ̃ e −1/κ̃∗ Jαˆ e ] − 3) I
C̄
1
+
− P̄0
(tr[
J0e
2κ̃∗
 e −1/κ̃∗
P̄0 J
+ ∗
(J e )−2/3 μe
2κ̃
J0e


1
e −1
× (C̄ ) ⊗ I − tr[C̄ e ](C̄ e )−1
3


1
e
e −1
e −1
⊗ (C̄ )
+ I − tr[C̄ ](C̄ )
3
ˆe
(19)
+ C̄

3

- 413 -

ˆ e は次の通りである．
上式の右辺最終項の C̄ˆ e := 2μe (J e )−2/3 1 tr[C̄ e ]I − 1 I ⊗ (C̄ e )−1
C̄
3
3

1
1
− (C̄ e )−1 ⊗ I + (C̄ e )−1 ⊗ (C̄ e )−13
9Ā = O または R = 1 のとき，T̄mod = T̄ となる．この修正応力の平均垂直応力とせん断応力はそれぞれ

1
3 dev
P̄mod = tr[T̄mod ], Q̄mod =
(24)
T̄
3
2 mod
dev
と与えられる．ここで T̄mod:= T̄mod − P̄mod I は修正応
力の偏差成分である．上式の P̄mod と Q̄mod を用いて，下負荷面関数を次のように表す．
(Q̄mod )2f (T̄ ; Ā, P̄c , R) :=+ P̄mod (P̄mod − RP̄c ) = 0
M2
(25)

(20)

式 (19), (20) では四階テンソル I を次のように定義している．

I := −
2.4∂(C̄ e )−1, II¯J¯K̄ L̄ = −((C̄ e )−1 )I¯K̄ ((C̄ e )−1 )L̄J¯∂ C̄ e(21)

2.5 塑性流れ則・内部変数の発展則・塑性適合条件(1) 塑性流れ則
塑性流れ則は次式で与える．
∂g(T̄ ; Ā, P̄c , R)L̄p = λ̇M̄ ,
M̄ =
(26)
∂ T̄
ここで λ̇ (≥ 0) は塑性乗数である．また，g は塑性ポテ
ンシャル関数であり，本研究では式 (25) の下負荷面関数と同じ関数を用いる (f = g)．この場合，式 (26) は関連流れ則となり，塑性散逸に関する T̄ と L̄p との共役性より，最大塑性散逸の原理を満たすモデルとなる．

修正応力の導入と拡張下負荷面 Cam-clay モデ
ルの正規降伏関数および下負荷面関数

ここでは Hashiguchi24) が提案した拡張下負荷面モデ
ルにならい，修正 Cam-clay モデル37) の正規降伏関数と下負荷面関数の具体形を示す．その際，降伏関数が応力の等方関数に限定されず，異方性も含めた一般的な塑性モデルにも適用しうるよう，付加的な剛体回転に対して不変な中間配置での Mandel 応力を用いる．正規降伏面に対する下負荷面の相似中心を考え，応力空間内において塑性変形の進展とともに移動する相似中心を用いて中間配置を参照する修正応力を導入し，これを用いて正規降伏関数と下負荷面関数を表す．

(1)

正規降伏関数

圧密降伏応力 P̄c の変化による硬化・軟化に対応して，応力空間において楕円体の形状を有する正規降伏面は等方的に拡大・縮小するものとし，次式の修正 Cam-clay
モデルの正規降伏関数を考える．

F (T̄ ; P̄c ) :=Q̄2
+ P̄ (P̄ − P̄c ) ≤ 0
M2

(22)

ここで Q̄ (:= 32 T̄ dev ) は Mandel せん断応力，T̄ dev(:= T̄ − P̄ I) は Mandel 応力の偏差成分，M は限界状態における応力比である．本研究では Lode 角の影響は考慮しないこととする．

(2)(3) 圧密降伏応力 P̄c の発展則
圧密降伏応力 P̄c の発展則（等方硬化則）は既に式(13) に示した通り，速度形ではなく時間積分された形で与えられる．

修正応力の導入と下負荷面関数

応力空間において常に正規降伏面の内側に位置し，正規降伏面と相似の形状を有し，相似の位置にある下負荷面を考える．正規降伏面と下負荷面に関する相似中
心を Ā とし，塑性変形の進展に伴って相似中心は後述の発展則に従って正規降伏面の内部を移動するものとする．この Ā は中間配置を参照するテンソルであり，中間配置の Mandel 応力の変換式 (5) と同様にして
a = (F e )−T · Ā · (F e )T と現配置表示に変換できるものとする．また，正規降伏面の大きさに対する下負荷面の大きさの比（以下ではこれを正規降伏比とよぶ）を R(0 < R ≤ 1) とする．0 < R < 1 は下降伏状態，R = 1
は正規降伏状態に相当する．以上の準備のもと，次式の修正 Mandel 応力 T̄mod を導入する23) ．

T̄mod := T̄ − (1 − R)Ā

(2) 相似中心 Ā の発展則
相似中心 Ā の発展則については，回転硬化モデルに対する Hashiguchi and Chen23) の形式を書き改めた次式を用いる．
˙
˙ = λ̇c(T̄ − Ā) + P̄c Ā
Ā
(27)
P̄c
ここで c は相似中心の移動の程度を表す材料定数である．なお，現配置表示の a に関して発展則を記述する
場合には，発展則の客観性を保証するために何らかの
客観応力速度を用いる必要がある．一方，本研究では等傾的な中間配置を仮定し，それを参照する剛体回転に対して不変な相似中心 Ā を用いており，物質時間微分で発展則を記述している．

(4) 正規降伏比 R の発展則
Hashiguchi24) は微小変形の枠組みにおいて，R の発展の程度は塑性ひずみ速度のノルムに比例するとして発展則を規定している．有限変形においては，塑性ひず
み速度を変形速度テンソルの塑性部分に置き換えて同様の発展則が用いられることが多い38) ．これらと同等の性質を有する発展則として，ここでは R の発展の程度が塑性乗数 λ̇ に比例するものと考え，次式を用いる．
Ṙ = λ̇r(R)(28)

ここで，関数 r(R) は条件

lim r(R) = ∞, r(R) > 0 for 0 < R < 1,

R→0

(23)
4

- 414 -

r(1) = 0(29)

を，右側から (Fnp )−1 をそれぞれ乗じ，塑性流れ則 (26)を代入すると，次の F e に関するリターンマッピング式を得る．

を満たす単調減少関数であり，幾つかの関数形が提案されており39) ，ここでは次式を用いる．

r(R) = −u0 ln R

(30)e,(tr)e
Fn+1= Fn+1 · exp(−ΔλM̄n+1 )u0 は R の発展の程度を表す正の材料定数である．
(5) 塑性適合条件
古典弾塑性モデルでは，負荷の持続に関する塑性適
合条件は正規降伏関数 F に対して Ḟ = 0 と与えられる
が，非古典弾塑性モデルでは下負荷面関数 f に対して
次のように与えられる．
∂f ˙
∂f ˙
f˙(T̄ ; Ā, P̄c , R) =: T̄ +
: Ā
∂ T̄
∂ Ā
∂f
∂f ˙
Ṙ = 0
(31)
P̄c +
+
∂R
∂ P̄c

3.

ここで，Δλ :=(:=p −1
Fn+1 · (Fn ) ) は Δt 間で塑性変形が生じないと仮定
した試行的な弾性変形勾配である3 ．また，exp[ΔλM̄ ]の逆テンソルは exp[−ΔλM̄ ] であることを用いた．

(2)

内部変数の発展則の増分表示

後退差分を用いて，相似中心 Ā の発展則 (27) を近似
し，次のように増分表示する．

Ān+1 − Ān = Δλc(T̄n+1 − Ān+1 )リターンマッピングによる陰的応力更新+

第 2 章で述べた有限変形・拡張下負荷面 Cam-clay モ
デルを数値解析に適用するためのリターンマッピングを用いた陰的な応力と内部変数の更新法を述べる．はじめに塑性流れ則と内部変数の発展則を増分表示し，これらを次ステップの未知の状態諸量について陰的に解くリターンマッピングを行い，下負荷面にのった応力状態を求める．また，境界値問題を解く際の反復収束計算において Newton–Raphson 法本来の二次収束を得るために不可欠な整合接線係数についても述べる．

eP̄c,n+1 = P̄c,n (Jn+1 /Jn+1)−Θe,(tr)

e,(tr)

Rn+1 − Rn = Δλr(Rn+1 )
3.2

(32)

(33)

k=0

1 2
1
Z + Z3 + · · ·
2
6

負荷判定

(tr)

(tr)

(tr)

(tr)

(tr)

fn+1 := f (T̄n+1 ; Ān+1 , P̄c,n+1 , Rn+1 )
(tr)

(40)

の値を計算し，fn+1 ≤ 0 ならば弾性除荷（あるいは中立）と判定し，正規降伏比 R 以外の試行値を時刻 tn+1
での更新値として採用する．ただし，弾性除荷の場合も下負荷面は現応力点が常に面上に位置するように追従するので，更新値が fn+1 = 0 を満たすように Rn+1を計算する必要がある．その方法は第 3.4 節に後述す(tr)る．また，fn+1 > 0 ならば塑性負荷と判定し，次節で述べるリターンマッピングにより fn+1 = 0 を満たす諸量の更新値を求める．

2 一般に二階テンソル Z の指数関数は次式により与えられる41) ．

Zk = I + Z +(39)

負荷判定の手順は以下の通りである．まず，試行
e,(tr)
的な弾性右 Cauchy–Green 変形テンソル C̄n+1 (:=e,(tr) Te,(tr)(Fn+1 ) · Fn+1 ) を計算し，超弾性構成式 (18) と応(tr)
力の変換則 (5) から試行 Mandel 応力 T̄n+1 を計算する．また，試行状態では塑性内部変数は時刻 tn から変(tr)化していないものと仮定し，試行値を Ān+1 := Ān ,
(tr)
(tr)
P̄c,n+1 := P̄c,n , Rn+1 := Rn とする．これらの試行値
を下負荷面関数 (25) に代入し，exp[ΔtL̄p ] はテンソル ΔtL̄p の指数関数である2 ．本来は L̄p は [tn , tn+1 ] 間でも変化するが，数値計算ではいず
れかの時点での値を用いることとなる．ここでは L̄pn+1eを用いた陰的近似をする．式 (33) 両辺に左側から Fn+1k!

e,(tr)

最後に，後退差分を用いて正規降伏比 R の発展則 (28)を近似し，次のように増分表示する．

と与えられる．Δt 間で L̄ が変化しないものとし，式(32) を [tn , tn+1 ] 間で時間積分すると次式を得る4),40) ．

∞

1

(38)

となる．ここで，Jn+1 := det Fn+1 = Jn+1 /Jnp であe,(tr)pee
り，Jn+1/Jnp = (Jn+1 /Jn+1)/Jnp = Jn+1 /Jn+1の関係を用いた．

p

exp[Z] =(37)

次に，等方硬化則 (13) は近似することなく解析的に増分表示に書き換えることができ，(1) テンソル指数関数を用いた塑性流れ則の時間積分式 (3) より，塑性変形勾配の発展式は

p
Fn+1= exp[ΔtL̄p ] · Fnp

P̄c,n+1 − P̄c,nĀn+1
P̄c,n+1

前述のように，相似中心 Ā の発展則の記述には客観速度ではなく，剛体回転に対して不変な中間配置を参照
する Ā の物質時間微分を用いているので，上式の差分により適切な近似がなされる．

3.1 塑性流れ則と内部変数の発展則の増分表示
増分解析において，第 n ステップの時刻 tn での応力と内部変数が既知とし，時間 Δt (:= tn+1 − tn ) 後の第n + 1 ステップで変形勾配 Fn+1 が与えられる場合を考える．以下では，速度形で与えられる塑性流れ則 (26)と発展則 (27), (28)，および等方硬化則 (13) を時間区
間 [tn , tn+1 ] に関する増分表示に書き換える．

Ḟ p = L̄p · F p(36)e,(tr)
λ̇Δt としており，また，Fn+1(34)

数値計算では無限級数のうち有限項を採用するが，第何項まで採用するかについては，例えば次の条件により決定できる．
1
Z k  < TOL(35)k!
本研究では TOL = 1.0 × 10−18 とした．k = 1 の項までを採用した場合には，一次の差分近似となる．塑性流れ則の時間積分の精度については，文献27) の第 5.3 節で検討されている．

3 上添字 ‘(tr)’ は試行値 (trial value) を意味する．

5- 415 -

3.3 リターンマッピング
式 (36) に含まれるテンソル指数関数を Q̄pn+1 :=eの更新式
exp[−ΔλM̄n+1 ] とおくと，式 (2) より，C̄n+1e,(tr)e= (Q̄pn+1 )T · C̄n+1 · Q̄pn+1C̄n+1のように更新し，収束条件

Y (X (k+1) ) < TOL(41)を得る．リターンマッピングの反復収束計算の過程で，
eを式 (18) に代入して求めるので，応力の更新値は C̄n+1
式 (36) よりも式 (41) を用いる方が計算上都合がよい．

を満たすまで上記の反復を繰り返す．本研究ではリター
ンマッピングの収束判定基準を TOL = 1.0 × 10−12 とした．なお，反復計算の開始時には，次のように各未知数の初期値を設定する．

リターンマッピングでは，式 (25), (37), (39), (41) をe, Ān+1 , Rn+1 , Δλ について解く
連立して未知量 C̄n+1とともに，式 (38) を用いて P̄c,n+1 を更新する．解くべき式を整理すると以下のようになる．

3.4

(45)

未知量ベクトル X ，残差ベクトル Y (X) をそれぞれ
Te
X := C̄n+1(46)
Ān+1 Rn+1 Δλ ,TY (X) := YC̄ e YĀ YR Yf(47)と形式的に定義すると，解くべき式はY (X) = 0(48)となる．非線形連立方程式 (48) を Newton–Raphson 法
を用いて反復計算により解く．式 (48) を線形化すると，∂Y
Y (X (k) ) +(49)(X (k) ) · δX = 0∂Xとなる．ここで δX は未知数 X の反復ごとの修正量であり，上添字 (k) は反復回数を表す．また，Jacobi 行列∂Y∂X (X) はリターンマッピングの収束性を左右する重要なものであるが，その導出は煩雑ではあるが直接的な微分演算4 により可能であるため，本論文では具体形は省略する．式 (49) を解くと，
−1∂YδX = −· Y (X (k) )(52)(X (k) )∂Xがある．一般的な二階テンソル Z の指数関数の微分は文献42)に与えられており，∂ exp[Z]
∂Z



=
ijkl

∞
n

1 
n=1

n!

(Z m−1 )ik (Z n−m )lj (50)(55)

除荷時の処置式 (23), (24) より，下負荷面関数 (25) は Rn+1 に関する二次方程式であり，これを解くと次のようになる．
√−b ± b2 − 4ac2
aRn+1 + bRn+1 + c = 0  Rn+1 =2a(56)
ここで係数 a, b, c は

3Ādev : Ādev
n+1
2M 2 n+1+ P̄A,n+1 (P̄A,n+1 − P̄c,n+1 ),
(57)
3
dev
dev
− Ādev
b := 2 (T̄n+1
n+1 ) : Ān+1
M+ (P̄n+1 − P̄A,n+1 )(2P̄A,n+1 − P̄c,n+1 ), (58)
3
dev
dev
c :=
(T̄ dev − Ādev
n+1 ) : (T̄n+1 − Ān+1 )2M 2 n+1+ (P̄n+1 − P̄A,n+1 )2
(59)

a :=

なお，R → 0 になると下負荷面の大きさがゼロに近づき，数値解析上の不安定が生じることがある．これを回避するため，本研究の数値解析では R の下限値 Rminを設定し，上述の手続きで求められた Rn+1 が Rmin を
下回った場合には，Rn+1 = Rmin として下負荷面の内部に現応力点が位置することを許容し，下負荷面内部
では弾性応答を示すものとして扱うこととする．

m=1である．級数の第何項まで採用するかについては，例えば次の条件により決定できる．

1 Z m−1 ⊗† Z n−m  < TOL
n!

(tr)

式 (56) より，Rn+1 について異なるふたつの解が求められることもある．しかし，R の定義域は 0 < R ≤ 1
であることと，除荷過程では下負荷面が収縮するので(tr)
Rn+1 ≤ Rn+1 (= Rn ) であることを考慮すれば，適切な一方の解を選択できる．

(53)

4 式 (42) を微分する際には，テンソル指数関数の微分を行う必要(0)

Rn+1 = Rn+1 ,であり，P̄A := 13 tr[Ā], Ādev := Ā − P̄A I とおいた．

となり，この修正量 δX を用いてX (k+1) := X (k) + δX

(tr)

(tr)

除荷時における正規降伏比 R の更新について述べる．
下負荷面モデルでは，除荷過程では弾性挙動を示すが，再負荷が開始すると直ちに塑性変形が発生する．すな
わち，除荷過程において下負荷面は現応力点が常に面上に位置するように収縮しつつ追従する．従って，試行値(tr)を代入した下負荷面関数が fn+1 < 0 となって除荷と判定された場合には，R 以外の諸量 T̄ , Ā, P̄c については
試行値をそのまま更新値として採用するが，fn+1 = 0
を満たすように Rn+1 を計算する必要がある．

P̄c,n+1 − P̄c,nĀn+1 = O, (43)P̄c,n+1

Yf := f (T̄n+1 , Ān+1 , P̄c,n+1 , Rn+1 ) = 0

(0)

(0)

Δλ(0) = 0

YĀ := Ān+1 − Ān − Δλc(T̄n+1 − Ān+1 )(44)e,(tr)

P̄c,n+1 = P̄c,n+1 ,e,(tr)

YR := Rn+1 − Rn − Δλr(Rn+1 ) = 0,e,(0)C̄n+1 := C̄n+1 , Ān+1 = Ān+1 ,e
YC̄ e := C̄n+1− (Q̄pn+1 )T · C̄n+1 · Q̄pn+1 = O, (42)

−

(54)

n

(51)

m=1

ここで，ふたつの二階テンソルに対する演算子 ‘⊗† ’ を (Z m−1 ⊗†Z n−m )ijkl = (Z m−1 )ik (Z n−m )lj としている．また，本研
究では TOL = 1.0 × 10−18 とした．

6

- 416 -

3.5

整合接線係数表–1 材料定数と初期状態量（繰返しせん断解析）

境界値問題の解析において Newton–Raphson 法本来の二次収束を得るためには，速度形構成式の弾塑性接線係数ではなく，リターンマッピングアルゴリズムに整合した接線係数を用いる必要がある．本来，整合接線係数は前節で述べたリターンマッピング諸式の微分演算により解析的に導出できるが，その導出は非常に煩雑である．一方，Miehe43) は摂動変形に対する応力の変化を計算することにより数値的に整合接線係数を求める方法を提案しており，本研究ではこの方法により整合接線係数を求める．適切な摂動変形を用いれば，解析的に導出した整合接線係数を用いた場合と比べても遜色なく良好な解の収束が得られることが報告されており43),44),45) ，とくに地盤材料の複雑な構成モデル
では，煩雑な導出を要しない点で有用である．
数値的に整合接線係数を計算する方法の詳細は文献43)に述べられており，文献27) の第 5.3 節でも紹介されているため，以下では概略のみを述べ，本研究で提案した弾塑性モデルに対する適用性を次章で検証する．

材料定数・初期状態量
基準圧力 P̄0
基準弾性体積変化 J0e
せん断弾性係数（定数部分） μe0
μe の圧力依存性パラメータ α
基準圧密降伏応力 P̄c0
限界応力比 M
弾性圧縮指数 κ̃∗弾塑性圧縮指数 λ̃∗
過圧密比 OCR
解析開始時の正規降伏比 R0
下負荷面の発展係数 u0
相似中心の発展係数 c
相似中心の初期値 Ā0

値−100 kPa
1.0
20.0 MPa0.0−400 kPa
1.05
0.025
0.082
4.0
0.14286
5.0
25.0−50I kPaとし，これに変形勾配 Fn+1 を与えた場合と，さらに摂動 ΔF(KL) ( ) を加えた場合との応力の差として

基準配置表示での増分形弾塑性構成式は
ep,algo 1
ΔS = Cn+1
: ΔC
(60)
2
と表される．ここで ΔS は右 Cauchy–Green 変形テンソルの増分 ΔC に対する第 2 Piola–Kirchhoﬀ 応力の増分であり，Cep,algo
は基準配置表示での整合接線係数でn+15ある ．上式より，整合接線係数の各成分は近似的に
ΔSIJ
Cep,algo(61)
n+1,IJKL ≈ 2
ΔCKL

ΔS(KL) ( ) = S(Fn+1 + ΔF(KL) ( ); {Fn , Hn })− S(Fn+1 ; {Fn , Hn })(65)と求められる．式 (64), (65) を式 (61) に適用すると，整合接線係数の各成分が次のように求められる．
1
Cep,algo
I, J, K, L = 1, 2, 3
n+1,IJKL ≈ ΔSIJ(KL) ( ),

(62)

(66)
四階テンソル Cep,algoの全成分を求めるには，式(66)
n+1
で (K, L) = (1, 1), (1, 2), (1, 3), (2, 2), (2, 3), (3, 3) の
6 方向の摂動変形を考え，各摂動に対する応力の変化を求めればよい．

式 (62) から分かるように，整合接線係数を精度良く導出するためには，ΔF をどのように与えるかが重要となる．そこで，なお，現配置を参照した updated Lagrange 法による有限要素法では，上述の方法で求めた基準配置表示のep,algoを
整合接線係数 Cn+1と求められる．ここで，ΔC は次の通りである．
T
ΔC = Δ(Fn+1· Fn+1 )T= Fn+1· ΔFn+1 + (ΔFn+1 )T · Fn+1

ΔF(KL) ( ) :=2

−1 T
) · E(K) ) ⊗ E(L)
(Fn+1
−1 T
+ (Fn+1) · E(L) ) ⊗ E(K)ep,algocep,algo
n+1,ijkl = FiI FjJ FkK FlL Cn+1,IJKL



と現配置表示 cep,algoに変換すればよい．
n+1

(63)

4.

という変形勾配の摂動を考えると，式 (62) の ΔC は次のようになる．

ΔC(KL) ( ) = 2 sym[E(K) ⊗ E(L) ](67)

応力更新法の適用性検証

本章では，第 3 章で述べたリターンマッピングを用いた陰的応力更新法を幾つかの解析例に適用し，計算精度および境界値問題の解析におけるロバスト性や解の収束性を検証する．

(64)

ここで， は摂動の大きさを表すパラメータで，この値は数値的に得られる整合接線係数の精度に影響するが44),45) ，本研究では = 1.0×10−8 とした．また，E(I)(I = 1, 2, 3) はデカルト座標系の基準配置表示での基底ベクトルである．式 (64) の増分に対する応力増分 ΔS
は，時刻 tn での変形勾配 Fn と塑性内部変数（ここでは形式的にまとめて Hn と表す）に対する状態を基準

4.1

ひずみ振幅制御による繰返しせん断変形解析

両振幅のせん断変形を与えるひずみ制御の繰返し変形解析を行った．ここでは変形勾配 F を⎤⎡
Hs0
1 + Hv⎥⎢
F = F∗ · F0 , F∗ = ⎣ 0
0 ⎦1 + Hv
00
1 + Hv(68)

5 上添字 ‘algo’ はアルゴリズム上の（リターンマッピングによる応力更新の手続きに整合した）接線係数を意味する．

7- 417 -

: Nstep = 45
: Nstep = 90
: Nstep = 900

表–2 材料定数と初期状態量（1 ステップ変形解析）

材料定数・初期状態量00ࡏࢇ᩿ኚᙧH S

0.001

図–2 ひずみ振幅制御による繰返しせん断変形の解析結果

ṇつ㝆అ㠃
(Case B)M
1ึᮇᛂຊ (Case C)
ୗ㈇Ⲵ㠃 (Case C)ึᮇᛂຊ (Case A, B)

0ṇつ㝆అ㠃 (Case C)┦ఝ୰ᚰࡢึᮇ್
(Case A, B, C)
ୗ㈇Ⲵ㠃 (Case B)

200

HS = 0.12

HS = 0.16

0Cauchyࡏࢇ᩿ᛂຊ q (kPa)

Cauchy ᖹᆒᆶ┤ᛂຊ㸪p (kPa)

図–3 Case A, B, C の初期応力点および正規降伏面・下負荷
面の様子

100

HS = 0.04

HS = 0.000と与える．ここで，F0 は初期応力状態を与える変形勾配，F∗ は初期状態に対する付加的な変形である．体積変形に対応する Hv とせん断変形に対応する Hs をそれ
ぞれ指定することにより，変形を制御した解析を行う．

HS = 0.12

HS = 0.16

HS = 0.08

HS = 0.20

0.05

0.1

: Nstep = 1000
: Nstep = 1

HV = +0.020
200

100

HS = 0.04

HS = 0.00
00HS = 0.12

HS = 0.08
0.05

HS = 0.20

100

HS = 0.08

HS = 0.16
HS = 0.20
0.1ᑐᩘࡏࢇ᩿ࡦࡎࡳ ε S

HS = 0.00−200−100
Cauchyᖹᆒᆶ┤ᛂຊ p (kPa)1000HS = 0.16 HS = 0.20
HS = 0.12
HS = 0.08
HS = 0.00

HS = 0.04

−200
−100
Cauchyᖹᆒᆶ┤ᛂຊ p (kPa)0: Nstep = 1000
: Nstep = 1

HV = +0.020
200

100

0

: Nstep = 1000
: Nstep = 1

HV = 0.000200ᑐᩘࡏࢇ᩿ࡦࡎࡳ ε S

Cauchyࡏࢇ᩿ᛂຊ q (kPa)

解析に用いた材料定数と初期状態量を表–1 に示す．
解析開始時の正規降伏比 R0 は，指定した初期応力，過圧密比，相似中心の初期値 Ā0 に基づき，応力空間において初期応力点が下負荷面上にあるという条件から，
式 (23), (24), (25) を用いて定められる．初期において−100 kPa の等方応力状態として，F0 = I とした．式(68) で Hs を −0.0015 ≤ Hs ≤ +0.0015 の範囲で変化させることにより 2 サイクルの繰返しせん断変形を与え，体積変形は Hv = 0 とした．変形増分の大きさによる解の精度を比較するため，増分ステップ数 Nstep = 45,
90, 900 の 3 通りで解析を行った．解析結果を図–2 に示す．異なる増分ステップ数でもほぼ一致した結果が得られ，本論文で提案した構成モデルと応力更新法が
繰返し変形においても有用であることが確認された．

0

HS = 0.16
HS = 0.12

HS = 0.04: Nstep = 1000
: Nstep = 1

200

200

0

0.05
0.1
ᑐᩘࡏࢇ᩿ࡦࡎࡳ ε S

HV = 0.000

0: Nstep = 1000
: Nstep = 1

HV = −0.020

HS = 0.04
HS = 0.00
0

−200
−200

HS = 0.20

HS = 0.08100

ṇつ㝆అ㠃࡜
ୗ㈇Ⲵ㠃 (Case A)−400: Nstep = 1000
: Nstep = 1

HV = −0.020
Cauchyࡏࢇ᩿ᛂຊ q (kPa)

Cauchy ࡏࢇ᩿ᛂຊ㸪q (kPa)

200

Cauchyࡏࢇ᩿ᛂຊ q (kPa)−0.001

Cauchyࡏࢇ᩿ᛂຊ q (kPa)−20

値−100 kPa
1.0
20.0 MPa40.0（本文記載）
1.05
0.025
0.082
（本文記載）（本文記載）
5.0
25.0−50I kPa

基準圧力 P̄0
基準弾性体積変化 J0e
せん断弾性係数（定数部分） μe0
μe の圧力依存性パラメータ α
基準圧密降伏応力 P̄c0
限界応力比 M
弾性圧縮指数 κ̃∗弾塑性圧縮指数 λ̃∗
過圧密比 OCR
解析開始時の正規降伏比 R0
下負荷面の発展係数 u0
相似中心の発展係数 c
相似中心の初期値 Ā0

Cauchyࡏࢇ᩿ᛂຊ q (kPa)

Cauchyᛂຊࡢ (1, 2) ᡂศσ12 (kPa)20

HS = 0.04
HS = 0.08
HS = 0.12
HS = 0.16
HS = 0.20

0HS = 0.00−200−100
Cauchyᖹᆒᆶ┤ᛂຊ p (kPa)0

図–4 1 ステップ解析と 1000 ステップ解析との比較 (Case A)

4.2

1 ステップで変形を与えたときの誤差評価

1 ステップで大きな変形を与えたときの応力計算の結果を検討する．ここでも 4.1 節と同様に変形を制御した解析を行い，膨張せん断や圧縮せん断など様々な変形について評価するため，式 (68) において体積変形 Hv
を −0.020 ≤ Hv ≤ +0.020 の範囲で 0.004 刻みで，せ8- 418 -

HS = 0.00
00.05

HS = 0.04

600

HS = 0.12

HS = 0.04

HS = 0.16

200

HS = 0.08

HS = 0.20

0.05
0.1ᑐᩘࡏࢇ᩿ࡦࡎࡳ ε S

600: Nstep = 1000
: Nstep = 1

Cauchyࡏࢇ᩿ᛂຊ q (kPa)

HV = +0.020

400

200

0

HS = 0.04

HS = 0.12

HS = 0.08
0

0
−600

0.05

HS = 0.16
HS = 0.20
0.1ᑐᩘࡏࢇ᩿ࡦࡎࡳ ε S

0

0

−400
−200
Cauchyᖹᆒᆶ┤ᛂຊ p (kPa)

HV = +0.020

300

HS = 0.16

100

HS = 0.08

00.05

HS = 0.20

300

300

200

100

0
0

Cauchyᖹᆒᆶ┤ᛂຊ p (kPa)

HS = 0.04

HS = 0.08

HS = 0.00 HS = 0.12
00.05

HS = 0.16

HS = 0.20

0.1ᑐᩘࡏࢇ᩿ࡦࡎࡳ ε S

−300

−200
−100
Cauchyᖹᆒᆶ┤ᛂຊ p (kPa)0: Nstep = 1000
: Nstep = 1

HV = 0.000

HS = 0.12

100

HS = 0.08
HS = 0.04−300

HS = 0.00

−200

−100

0

Cauchyᖹᆒᆶ┤ᛂຊ p (kPa): Nstep = 1000
: Nstep = 1

HV = +0.020

HS = 0.00

HS = 0.20
HS = 0.16

00.1

HS = 0.20

HS = 0.16
HS = 0.12
HS = 0.08200ᑐᩘࡏࢇ᩿ࡦࡎࡳ εS

HS = 0.20
HS = 0.16
HS = 0.12−200

HS = 0.12

HS = 0.04

0HS = 0.04
HS = 0.08

−400

200

0

400

0
−600

0

0.1

: Nstep = 1000
: Nstep = 1

HV = 0.000

HS = 0.00: Nstep = 1000
: Nstep = 1

200

0.05

: Nstep = 1000
: Nstep = 1

HV = −0.020

100 H S = 0.04

ᑐᩘࡏࢇ᩿ࡦࡎࡳ ε S

HS = 0.12
HS = 0.08
HS = 0.04
HS = 0.00

600

Cauchyࡏࢇ᩿ᛂຊ q (kPa)

0

−200

HS = 0.20
HS = 0.16

200

200

HS = 0.00

400

HS = 0.00
0HS = 0.08
0: Nstep = 1000
: Nstep = 1

HV = 0.000
Cauchyࡏࢇ᩿ᛂຊ q (kPa)

Cauchyࡏࢇ᩿ᛂຊ q (kPa)

400

HS = 0.20

100

Cauchyᖹᆒᆶ┤ᛂຊ p (kPa): Nstep = 1000
: Nstep = 1

HV = 0.000−400

HS = 0.16

HS = 0.04

HS = 0.00ᑐᩘࡏࢇ᩿ࡦࡎࡳ ε S
600

HS = 0.12

200

Cauchyࡏࢇ᩿ᛂຊ q (kPa)

200

0
−600

0.1

HS = 0.16
HS = 0.12

HS = 0.08

300

Cauchyࡏࢇ᩿ᛂຊ q (kPa)0

HS = 0.08

HS = 0.20

400: Nstep = 1000
: Nstep = 1

HV = −0.020

300

Cauchyࡏࢇ᩿ᛂຊ q (kPa)

HS = 0.16

200

Cauchyࡏࢇ᩿ᛂຊ q (kPa)

HS = 0.04

300

Cauchyࡏࢇ᩿ᛂຊ q (kPa)

HS = 0.12: Nstep = 1000
: Nstep = 1

HV = −0.020
Cauchyࡏࢇ᩿ᛂຊ q (kPa)

Cauchyࡏࢇ᩿ᛂຊ q (kPa)

HS = 0.20

400

600

: Nstep = 1000
: Nstep = 1

HV = −0.020

Cauchyࡏࢇ᩿ᛂຊ q (kPa)600

HV = +0.020: Nstep = 1000
: Nstep = 1

200

100

0

HS = 0.12
HS = 0.16
HS = 0.20
HS = 0.00−300

HS = 0.04
HS = 0.08

−200
−100
Cauchyᖹᆒᆶ┤ᛂຊ p (kPa)0

図–6 1 ステップ解析と 1000 ステップ解析との比較 (Case C)

図–5 1 ステップ解析と 1000 ステップ解析との比較 (Case B)(p = −100 kPa, q = 0 kPa), 過圧密比 OCR = 1.0,
解析開始時の正規降伏比 R0 = 1.0.
• Case B: 等方過圧密状態．圧密降伏応力 P̄c0 =−400 kPa, 初期応力 σ11 = σ22 = σ33 = −100 kPa(p = −100 kPa, q = 0 kPa), 過圧密比 OCR = 4.0,
解析開始時の正規降伏比 R0 = 0.14286.
• Case C: 異方過圧密状態．圧密降伏応力 P̄c0 =−200 kPa, 初期応力 σ11 = σ22 = −50 kPa, σ33 =−100 kPa (p = −66.7 kPa, q = 50 kPa), 過圧
密比 OCR = 1.6274, 解析開始時の正規降伏比
R0 = 0.48596.
各ケースでの初期応力点および正規降伏面・下負荷面の様子を図–3 に示す．解析に用いた材料定数と初期状
態量を表–2 に示す．Case A, B では J0e = 1.0, F0 = Iとした．一方，Case C では J0e = 1.0 とし，初期の異方
過圧密状態を再現するため，F0 を次のように設定した．
⎡⎤
3.68505
00⎢⎥
F0 = ⎣ 0
3.68505
0 ⎦ × 10−3 (70)002.56835

ん断変形 Hs を 0 ≤ Hs ≤ 0.20 の範囲で 0.04 刻みでそれぞれ設定した．これらの値は対数ひずみに換算すると，εv = ±0.06, εs = 0.12 程度である．この変形勾配
F ∗ を 1 ステップで与えたときの応力計算値と正解との相対誤差を次式で評価し，誤差マップを作成した．
(p − p )2 + (q − q  )2
ERR :=
× 100 (%) (69)(p )2 + (q  )2
ここで，p, q はそれぞれ 1 ステップ解析で得られたCauchy 応力 σ (:= τ /J) の平均垂直応力とせん断応力である．また，上添字 ( ) が付いた量は正解であるが，本論文で扱うような複雑な構成モデルでは解析的方法により正解を求めるのは困難であるため，ここでは所定の変形勾配を F ∗ = f · · · f （1000 個の f の積）となるように等分割し，1 ステップあたりの変形勾配 f
を 1000 ステップ与えたときの解析結果を正解の代用とした．ただし，変形勾配 F ∗ の分割の仕方は一意的ではなく，ここでの検討は，上述のように分割した変形経路での結果と 1 ステップ解析とを比較した限定的な誤差評価であることには注意を要する．

図–4, 5, 6 に Case A, B, C それぞれの解析結果の一部をせん断ひずみ−せん断応力関係，および p–q 平面上の応力経路として示す．図中の実線が 1000 ステップ
解析の結果を，「◦」印が 1 ステップ解析の結果を示す．
応力経路の方向が複雑に変化する場合でも，1 ステップ

初期応力状態による比較を行うため，次に示す 3 つの初期状態について解析を行った．

• Case A: 等方正規圧密状態．圧密降伏応力 P̄c0 =−100 kPa, 初期応力 σ11 = σ22 = σ33 = −100 kPa
9- 419 -

ࡏࢇ᩿ኚᙧ㸪H S

20.0
18.8
17.5
16.3
15.0
13.8
12.5
11.3
10.0
8.75
7.50
6.25
5.00
3.75
2.50
1.25
+0.02 0.00

20.0
18.8
17.5
16.3
15.0
13.8
12.5
11.3
10.0
8.75
7.50
6.25
5.00
3.75
2.50
1.25
+0.02 0.00┦ᑐㄗᕪ (%)

0.00┦ᑐㄗᕪ (%)

0.0
−0.02

20.0
18.8
17.5
16.3
15.0
13.8
12.5
11.3
10.0
8.75
7.50
6.25
5.00
3.75
2.50
1.25
+0.02 0.00┦ᑐㄗᕪ (%)

0.2

表–3 材料定数と初期状態量（一面せん断解析）

材料定数・初期状態量య✚ኚᙧ㸪H V(a) Case Aࡏࢇ᩿ኚᙧ㸪H S

0.2

0.0
−0.02

0.00
య✚ኚᙧ㸪H V(b) Case Bࡏࢇ᩿ኚᙧ㸪H S

0.2

0.0
−0.02

0.00

値
−100 kPa
1.0
20.0 MPa100.0（本文記載）
1.05
0.025
0.082
（本文記載）（本文記載）
5.0
25.0−50I kPa

基準圧力 P̄0
基準弾性体積変化 J0e
せん断弾性係数（定数部分） μe0
μe の圧力依存性パラメータ α
基準圧密降伏応力 P̄c0
限界応力比 M
弾性圧縮指数 κ̃∗弾塑性圧縮指数 λ̃∗
過圧密比 OCR
解析開始時の正規降伏比 R0
下負荷面の発展係数 u0
相似中心の発展係数 c
相似中心の初期値 Ā0ୖ㍕ᅽ σn = 100 kPa㸦୍ᐃ㸧㖄┤ኚ఩ vỈᖹኚ఩ uỈᖹⲴ㔜 Fయ✚ኚᙧ㸪H V

1.8cm(c) Case C
図–7 誤差マップ

4.0cm 0.4cm
1.8cm

y

解析で計算された応力点と 1000 ステップ解析での最終応力点に大きな違いはなかった．

6.0cm

図–8 一面せん断解析の有限要素モデル

横軸に Hv を，縦軸に Hs を取り，式 (69) で評価した相対誤差 (%) を示した誤差マップを図–7 に示す．大きな膨張変形（Hv = +0.02 程度）を与えたときに，
Case C では 1 ステップ解析と 1000 ステップ解析で応力計算値に最大 20 %程度の差異が生じているのを除けば，1 ステップ解析と 1000 ステップ解析で概ね一致した応力が計算されている．なお，Case B の入力変形(Hv , Hs ) = (+0.020, 0.00) では 1 ステップ解析の過程
で解が発散したため，図に含めていない．

ỈᖹⲴ㔜㸪F (N)50

Case C

40

Case B

30
20

Case A

100㖄┤ኚ఩㸪v (cm)

1 ステップ解析と 1000 ステップ解析とで計算量を比較するために，収束条件 (54) を満たすのに要した反復計算の回数を調べる．解析条件が最も厳しいと思われる Case C の (Hv , Hs ) = (+0.020, 0.20) のケースでは，
1 ステップ解析で 14 回，1000 ステップ解析で合計 4996
回（1 ステップあたり最大 5 回，最小 4 回）であった．
大きな変形増分で解析を行った場合には，図–7 で示したように応力計算結果に差異が生じることがあり，また，境界値問題の解析ではつり合い状態の解を求めるために多くの反復計算を要することも考えられるが，その点に注意すれば，提案手法により大きな変形増分での解析が可能となり，計算効率の向上を図ることができる．

60

Case C

0.02

Case B0−0.02

Case A

−0.04
−0.06

0

0.1

0.2

Ỉᖹኚ఩㸪u (cm)

0.3

0.4

図–9 水平変位と水平荷重・鉛直変位との関係

4.3

境界値問題への適用：一面せん断の例

二次元平面ひずみ条件下での一面せん断解析を行い，提案手法の境界値問題への適用性と解の収束性を検証する．この問題は応力・変形の主軸の回転と変形の局所化

10- 420 -

Case C

Case B

Case A

に設定した．このように，荷重の低下や変形の局所化を伴う状況であっても Newton–Raphson 法本来の二次収束を示しており，本論文で提案したリターンマッピン
グを用いた応力更新法の適用性・有用性が確認できた．

0.000

0.550

5.

ᑐᩘࡏࢇ᩿ࡦࡎࡳ  εs

図–10 一面せん断解析での変形図（せん断変位 u = 0.38 cm）

10ṧᕪࣀ࣒ࣝ㸪||R ||10

本研究の成果をまとめると，以下の通りである．

1. 変形勾配の乗算分解に基づく有限変形理論の枠組
みで，拡張下負荷面 Cam-clay モデルの定式化を行った．その際，中間配置を参照する諸量を用いて降伏関数と塑性発展式を記述した．
2. 地盤材料に特有の圧力依存性を有する体積圧縮特性とせん断特性を表現できる超弾性モデルを提案した．その際，従来提案されている対数主ひずみと主応力とを関係付ける方法では等方性モデルの記述に限定されるため，これを拡張し，中間配置を参照するテンソルを用いて構成式を記述した．
3. 数値解析での利用を前提としたリターンマッピン
グによる効率的な応力更新法を開発し，これを有限要素解析コードに実装した．
4. 提案手法を解析例題に適用し，精度検証を行った．
その結果，比較的大きな変形増分を与えても良好な精度で応力計算が可能であることを示した．また，境界値問題の解析に提案手法を適用し，厳しい解析条件においても良好な解の収束がみられ，本手法のロバスト性と高い適用性が確認された．

−2

−5

10

−10

10

−15

཰᮰ุᐃᇶ‽

02

4඲య཯᚟ᅇᩘ k68

図–11 Case C の荷重極大点付近での解の収束状況を伴うものであり，提案手法の妥当性検証に適している．
有限要素モデルを図–8 に，解析に用いた材料定数と初期状態量を表–3 にそれぞれ示す．モデルの奥行き方向の寸法は 1 cm とした．有限要素は 8 節点二次要素（3×3 点完
全 Gauss 積分）を用い，要素数は 300 である．初期状態
では等方応力状態として σ11 = σ22 = σ33 = −100 kPa(p = −100 kPa, q = 0 kPa) と設定した．土の自重は無視した．上載圧は σn = 100 kPa で一定として，水平せん断変位 u を与えたときの水平荷重 F と鉛直変位 v を求めた．比較のため，次に示す 3 つの初期状態について解析を行った．

• Case A: 等方正規圧状態 (圧密降伏応力 P̄c0 =−100 kPa, 過圧密比 OCR = 1.0)
• Case B: 等方過圧密状態 (圧密降伏応力 P̄c0 =−200 kPa, 過圧密比 OCR = 2.0)
• Case C: 等方過圧密状態 (圧密降伏応力 P̄c0 =−400 kPa, 過圧密比 OCR = 4.0)
図–9 に水平変位と水平荷重および鉛直変位との関係を示す．また，せん断変位 u = 0.38 cm 時の変形図（コンターは対数せん断ひずみ）を図–10 に示す．Case A では荷重は単調に増加しており，せん断変形の分布は拡散的であるが，Case B と Case C では荷重極大点とその後の荷重低下がみられ，せん断変形は狭い帯状領域に局所化している．とくに Case C では荷重極大点の直後から急激に荷重が低下しており，厳しい解析条件になっている．そこで，図–9 図中の「◦」印を示した付近の 3ステップでのつり合い点の反復求解における解の収束状
況を図–11 に示す．収束判定基準は TOL = 1.0 × 10−10

まとめ

謝辞： 本研究の遂行に際し，科学研究費補助金（若手研究 (B), 研究代表者: 山川優樹，平成 20〜22 年度，
課題番号 20760311）の補助を受けた．また，（株）豊田
中央研究所の渡邊育夢氏には有限変形理論に関して貴重な助言を頂いた．ここに記して謝意を表する．
参考文献
1) Simo, J.C., Ortiz, M.: A uniﬁed approach to ﬁnite
deformation elastoplastic analysis based on the use of
hyperelastic constitutive equations. Comput. Methods
Appl. Mech. Engrg. 49, 221–245, 1985.
2) Simo, J.C., Taylor, R.L.: Consistent tangent operatorsfor rate-independent elastoplasticity. Comput. Meth.
Appl. Mech. Engrg. 48, 101–118, 1985.
3) Ortiz, M., Simo, J.C.: An analysis of a new class of integration algorithms for elastoplastic constitutive relations. Int. J. Numer. Methods Engrg. 23, 353–366, 1986.
4) Weber, G., Anand, L.: Finite deformation constitutive equations and a time integration procedurefor isotropic, hyperelastic–viscoplastic solids. Comput.
Methods Appl. Mech. Engrg. 79, 173–202, 1990.
5) Papadopoulos, P., Lu, J.: A general framework for thenumerical solution of problems in ﬁnite elasto-plasticity.
Comput. Methods Appl. Mech. Engrg. 159, 1–18, 1998.
6) Amirkhizi, A.V., Nemat-Nasser, S.: A framework fornumerical integration of crystal elasto-plastic constitutive equations compatible with explicit ﬁnite element
codes. Int. J. Plasticity 23, 1918–1937, 2007.
7) de Souza Neto, E.A., Perić, D., Owen, D.R.J.: Computational Methods for Plasticity: Theory and Applications. John Wiley & Sons, 2008.

11- 421 -

1989.
25) Lee, E.H.: Elastic–plastic decomposition at ﬁnite
strains. ASME J. Appl. Mech. 36, 1–6, 1969.
26) Simo, J.C.: Numerical analysis and simulation of plasticity. in: Ciarlet, P.G., Lions, J.L. (eds.), Handbook of
Numerical Analysis, Vol. VI, Elsevier Science, 183–499,
1998.
27) 地盤工学会編: 土の弾塑性構成モデル（地盤工学・基礎
理論シリーズ 3）．地盤工学会, 2009.
28) Casey, J., Naghdi, P.M.: A remark on the use of thedecomposition F = Fe Fp in plasticity. ASME J. Appl.
Mech. 47, 672–675, 1980.
29) Dashner, P.A.: Invariance considerations in large
strain elasto-plasticity. ASME J. Appl. Mech. 53, 55–
60, 1986.
30) Mandel, J.: Plasticité Classique et Viscoplasticité.
CISM Course 97, Springer, Udine, 1971.
31) Hashiguchi, K.: On the linear relations of V –ln p andln v–ln p for isotropic consolidation of soils. Int. J. Numer. Anal. Methods Geomech. 19, 367–376, 1995.
32) Callari, C., Auricchio, F., Sacco, E.: A ﬁnite-strain
Cam-clay model in the framework of multiplicative
elasto-plasticity. Int. J. Plasticity 14, 1155–1187, 1998.
33) Houlsby, G.T.: The use of variable shear modulus inelastoplastic model for clays. Comput. Geotech. 1, 3–13,
1985.
34) Borja, R.I., Tamagnini, C., Amorosi, A.: Coupling
plasticity and energy-conserving elasticity model forclays. ASCE J. Geotech. Geoenvironmental Engrg. 123,
948–957, 1997.
35) Tamagnini, C., Castellanza, R., Nova, R.: A Generalized Backward Euler algorithm for the numerical integration of an isotropic hardening elastoplastic model
for mechanical and chemical degradation of bonded geomaterials. Int. J. Numer. Anal. Methods Geomech. 26,
963–1004, 2002.
36) Niemunis, A., Cudny, M.: On hyperelasticity for clays.
Comput. Geotech. 23, 221–236, 1998.
37) Roscoe, K.H., Burland, J.B.: On the generalized
stress–strain behavior of ‘wet’ clay. in: Heyman, J.,
Leckie, F.A. (eds.), Engineering Plasticity. Cambridge
University Press, 535–609, 1968.
38) Hashiguchi, K.: The extended ﬂow rule in plasticity.
Int. J. Plasticity 13, 37–58, 1997.
39) 橋口公一: 最新弾塑性学. 朝倉書店, 1990.
40) Miehe, C.: Exponential map algorithm for stress updates in anisotropic multiplicative elastoplasticity forsingle crystals. Int. J. Numer. Methods Engrg. 39, 3367–
3390, 1996.
41) 梁成吉: キーポイント 行列と変換群（理工系数学のキーポイント 8）．岩波書店, 1996.
42) de Souza Neto, E.A.: The exact derivative of the exponential of an unsymmetric tensor. Comput. Methods
Appl. Mech. Engrg. 190, 2377–2383, 2001.
43) Miehe, C.: Numerical computation of algorithmic(consistent) tangent moduli in large-strain computational inelasticity. Comput. Methods Appl. Mech. Engrg. 134, 223–240, 1996.
44) Pérez-Foguet, A., Rodrı́guez-Ferran, A., Huerta, A.:
Numerical diﬀerentiation for local and global tangent
operators in computational plasticity, Comput. Methods
Appl. Mech. Engrg. 189, 277–296, 2000.
45) Fellin, W., Ostermann, A.: Consistent tangent operators for constitutive rate equations. Int. J. Numer. Anal.
Methods Geomech. 26, 1213–1233, 2002.

8) Nagtegaal, J.C., de Jong, J.E.: Some aspects of nonisotropic work-hardening in ﬁnite strain plasticity. in:
Lee, E.H., Mallet, R.L. (eds.), Plasticity of Metals at Finite Strain: Theory, Experiment and computation. Division of applied mechanics, Stanford Univ. and Dept.
of engineering & mechanics, Rensselaer polytechnic institute, 65–102, 1982.
9) Kojić, M., Bathe, K.-J.: Studies of ﬁnite element
procedures—Stress solution of a closed elastic strain
path with stretching and shearing using the update
Lagrangian Jaumann formulation. Comput. Struct. 26,
175–179, 1987.
10) Hughes, T.J.R.: Numerical implementation of constitutive models: Rate-independent deviatoric plasticity.
in: Nemat-Nasser, S. et al. (eds.), Theoretical Foundations for Large-scale Computations for Nonlinear Material Behavior. Dordrecht, The Netherlands, Martinus
Nijhoﬀ, 29–63, 1984.
11) Borja, R.I., Tamagnini, C.: Cam–Clay plasticity,
Part III: Extension of the inﬁnitesimal model to include ﬁnite strains. Comput. Methods Appl. Mech. Engrg. 155, 73–95, 1998.
12) Meschke, G., Liu, W.N.: A re-formulation of the exponential algorithm for ﬁnite strain plasticity in terms ofcauchy stresses. Comput. Methods Appl. Mech. Engrg.
173, 167–187, 1999.
13) Yamakawa, Y., Hashiguchi, K., Ikeda, K.: Implicit
stress-update algorithm for isotropic Cam-clay model
based on the subloading surface concept at ﬁnite strains.
Int. J. Plasticity 26, 634–658, 2010.
14) Simo, J.C.: Algorithms for static and dynamic multiplicative plasticity that preserve the classical return
mapping schemes of the inﬁnitesimal theory. Comput.
Methods Appl. Mech. Engrg. 99, 61–112, 1992.
15) Miehe, C.: On the representation of Prandtl–Reuss
tensors within the framework of multiplicative elastoplasticity. Int. J. Plasticity 10, 609–621, 1994.
16) Ekh, M., Runesson, K.: Modeling and numerical issues in hyperelasto-plasticity with anisotropy. Int. J.
Solids Struct. 38, 9461–9478, 2001.
17) Menzel, A., Ekh, M., Runesson, K. and Steinmann, P.:
A framework for multiplicative elastoplasticity withkinematic hardening coupled to anisotropic damage. Int.
J. Plasticity 21, 397–434, 2005.
18) Wallin, M., Ristinmaa, M. and Ottosen, N.S.: Kinematic hardening in large strain plasticity. European J.
Mech. A/Solids 22, 341–356, 2003.
19) Wallin, M., Ristinmaa, M.: Deformation gradient
based kinematic hardening model. Int. J. Plasticity 21,
2025–2050, 2005.
20) Jeremić, B., Runesson, K., Sture, S.: Finite deformation analysis of geomaterials. Int. J. Numer. Anal.
Methods Geomech. 25, 809–840, 2001.
21) Rouainia, M., Muir Wood, D.: Implicit numerical integration for a kinematic hardening soil plasticity model.
Int. J. Numer. Anal. Methods Geomech. 25, 1305–1325,
2001.
22) Borja, R.I., Lin, C.-H., Montáns, F.J.: Cam–Clayplasticity, Part IV: Implicit integration of anisotropic
bounding surface model with nonlinear hyperelasticity and ellipsoidal loading function. Comput. Methods
Appl. Mech. Engrg. 190, 3293–3323, 2001.
23) Hashiguchi, K., Chen, Z.-P.: Elastoplastic constitutive equation of soils with the subloading surface andthe rotational hardening. Int. J. Numer. Anal. Methods
Geomech. 22, 197–227, 1998.
24) Hashiguchi, K.: Subloading surface model in unconventional plasticity. Int. J. Solids Struct. 25, 917–945,（2010 年 3 月 9 日 受付）

12- 422 -