3.3 応力ひずみ関係の数値積分における非線形反復計算法
3.3.1 目的
前節で示した二相系の運動方程式の解法に対し、第 2 章で示した、砂の繰り返し載荷時 の挙動を表現する力学モデル（2.2 節、2.3 節参照）を適用して地震応答解析を行うと、し ばしば応答加速度時刻歴にスパイク状のピークが発生する。スパイク状の加速度は、後述 するように１次元地盤モデルを用いた地震応答解析においても発生する。この場合は、過 剰間隙水圧上昇時に応答加速度時刻歴にスパイク状のピークが多発し、この傾向は、時間 積分間隔Δtを小さくとるといっそう激しくなる（図 3.4-5①参照）。また、盛土や矢板式岸 壁などの地震応答解析を行う場合にも、特に、液状化する可能性のある基礎地盤が常に軸 差応力の作用下にあるような場合に、しばしば応答加速度時刻歴に異常振動が見られる。
変相線を超えた領域における塑性せん断仕事の負のダイレタンシーへの寄与程度を見直し た砂の力学モデル（修正モデル：2.3節参照）を使用した場合にはこの傾向は著しいものの
（図 3.4-12参照）、従来モデルを使用した場合にも同様の現象が見られる。これらの状況で 共通するのは、砂のダイレタンシーの影響でせん断応力の作用により有効応力が変動し、
これによって応力－ひずみ関係が短時間に大きく変動すると言うことである。
従来は各時間ステップ（または荷重ステップ）内では砂の応力－ひずみ関係は変動しな いと言う前提の基に（2.2.3 項参照）非線形反復計算を行っていた。この方法によれば、静 止している応力－ひずみ関係を対象に非線形反復計算を行うので、各時間ステップでの収 束性に優れていて、その応力－ひずみ関係を前提とする限りにおいては、当該ステップの 非線形反復計算収束時の要素力は、対応する荷重と平衡状態にある。次のステップに進ん だ場合、特に、過剰間隙水圧モデル（2.2.2 項参照）を適用する場合は、前のステップの収 束状態に応じて更新された液状化フロントS0や状態変数Sにより定まる応力－ひずみ関係 を前提に収束計算が開始されるが、これだと、たとえ荷重変動が無くても、応力－ひずみ 関係の変動による不平衡力が生じることになる（図 3.3-1参照）。Newmark 法などの時間積 分法は、不平衡力を解消するために慣性力も動員して釣り合いを保とうとするので、それ に対応して不自然な加速度応答が生じることになる。
上述の応答加速度時刻歴にスパイク状のピークが多発する現象は、各時間ステップ内で 応力－ひずみ関係を不変と仮定したことに起因して発生するものと考え、本節では、前章 で示した多重せん断機構モデルに即して、各時間ステップ内で応力－ひずみ関係の変動を 追跡しながら非線形反復計算を行う方法を提案し、その効果を検証する。
第 3 章 地盤・構造物系の数値解析法
図 3.3-1 不平衡力発生の概念図（各時間ステップ内では骨格曲線は不変とした場合）
骨格曲線の変動は、主に、液状化フロントパラメータS0か状態変数Sの変動による。
図 3.3-2 各時間ステップ内で状態変数Sの変動による骨格曲線の変動に 追随する非線形反復計算法の概念図
実線矢印のように、各時間ステップ内で骨格曲線を乗り移って行く。
破線矢印は従来法で、同じ骨格曲線上を移動する。
τ
τ
ｆG
0γ
ｒｅｆγ
←ｔ時の骨格曲線
←ｔ＋Δｔ時の骨格曲線
収束点 不平衡力
τ
ｆG
0γ
ｒｅｆγ
τ
←ｔ＋Δｔ時の骨格曲線
←ｔ時の骨格曲線
3.3 応力ひずみ関係の数値積分における非線形反復計算法
3.3.2 前提条件
第 2 章で述べた砂の力学モデルによると、せん断応力－せん断ひずみ関係の変動は、以 下の要因により生じる。
① 過剰間隙水圧モデル（2.2.2項）を適用しない場合は、平均有効応力σm'の変動。
※(2.2-68)式、(2.2-69)式を参照
② 過剰間隙水圧モデルを適用する場合は、液状化フロントパラメータS0の更新。
※(2.2-39)式、(2.2-40)式を参照
③ 過剰間隙水圧モデルを適用する場合は、状態変数Sの変動。
※(2.2-39)式、(2.2-40)式を参照
本節では、上記③の要因による各時間ステップ内の骨格曲線の変動に対して、その変動 を前提とした非線形反復計算手法を示す。すなわち、収束時には、Sの変動による骨格曲線 の変動は折り込み済みで、これによる次ステップへの不平衡力の持ち越しは生じないよう なアルゴリズムを示す。図 3.3-2に、この方法の概念図を示す。
土の挙動を表現するのに、過剰間隙水圧モデルを適用せず、多重せん断ばねモデルのみ を適用した場合には、応答加速度時刻歴に異常波形は見られないので、①の理由による骨 格曲線などの変動への対策は取らない。また、変相線を超えた応力空間における塑性せん 断仕事は負のダイレタンシーへ寄与しないと言う修正モデル（2.3節参照）を用いた場合で、
応力ポイントが破壊線近傍にあり、従って、液状化フロントパラメータ S0が変動しないよ うな場合でも（②の要因が考えられない場合でも）、加速度時刻歴には相当な異常波形が見 られることがあるので、③の理由による骨格曲線などの変動に対処することは必要と考え る。問題は、②の要因による骨格曲線などの変動による不平衡力発生に対して、対策を取 らなくてもよいかと言うことになるが、取扱いが困難なので、ここでは、③による変動対 策のみを行う。
従って、ここで示す方法では、各時間ステップ内の非線形反復計算においては、液状化 フロントパラメータ S0は、従来通り固定されている。また、次ステップにおいて、液状化 フロントパラメータS0の更新による不平衡力が、従来通り生じることになる。
一方、状態変数 S は、各時間ステップ内の非線形反復計算においては、随時、更新する こととし、次ステップにおいては、Sの更新による不平衡力は生じない。
第 3 章 地盤・構造物系の数値解析法
3.3.3 骨格曲線の変動を考慮に入れた応力ひずみ関係の接線勾配
(1) 定式化Ⅰ（∂σ／∂εの導出）
前章で示した多重せん断機構モデルに含まれる過剰間隙水圧モデル（2.2.2 項）を用いる 場合の非線形反復計算で必要な∂{σ'}/∂{ε}の算定式を導く。ここに、{σ'}T = (σx',σy', τxy)、{ε}T=(εx,εy,γxy)である。また、前述の通り、非線形反復計算の途中で、状態変数S の値は随時変更するが、液状化フロントS0の値は一定として取り扱う。
多重せん断機構モデルでは、偏差応力－偏差ひずみ関係は、以下のように表される。
θ Δ θ σ γ
τ σ n i
i θ
x y
d 2 F( i)cos
2
2 1 '
'
∑
=
− =
≡ (3.3-1)
θ Δ θ γ
τ n i
i θ
xy 22 F( i)sin
∑
1=
= (3.3-2)
但し、n：多重せん断ばねモデルの1/4円当たりのばねの数
Δθ＝π／2n (3.3-3)
i xy i d
θi γ θ γ θ
γ = cos + sin (3.3-4)
θ Δ
θi =(i−1) (3.3-5)
また、
x y
d ε ε
γ ≡ − (3.3-6)
である。また、偏差ひずみγd,γxy には、先行する各解析フェーズにおけるひずみの累積も 含まれている。
(3.3-1, 2, 4)式より、偏差応力τd,τxyは、ともに偏差ひずみγd,γxyの関数であり、また、
後で示すように、F(γ)の形は、初期せん断剛性G0とせん断強度τfにより決まる。従って、
τdとτxyを以下のように表す。
) , , ,
( d xy 0 f
d
d τ γ γ G τ
τ = (3.3-7)
) , , ,
( d xy 0 f
xy
xy τ γ γ G τ
τ = (3.3-8)
また、γdとγxyは、Bマトリックスを介して、節点変位ベクトルUと関連付けられる。す なわち、γdとγxyはUの関数である。
) (U
d
d γ
γ = (3.3-9)
)
xy(U
xy γ
γ = (3.3-10)
G0とτfは、過剰間隙水圧モデル（2.2.2 節）によると、液状化フロントS0と状態変数S の 関数であるが、各時間ステップ内ではS0は一定と見なすので、状態変数Sのみの関数とお
3.3 応力ひずみ関係の数値積分における非線形反復計算法 く。さらに、Sは、S0およびτdとτxyの関数であるが、やはり、各時間ステップ内では S0は一定と見なすので、τdとτxyの関数とする。
)
0(
0 G S
G = (3.3-11)
)
f(S
f τ
τ = (3.3-12)
) , ( d xy S
S = τ τ (3.3-13)
τdとτxyは、ともに偏差ひずみγd,γxyの関数であり、偏差ひずみは節点変位ベクトルU の関数である。結局、唯一の独立変数は節点変位ベクトルである。
ここで、dUに対するτdとτxyの増分を求める。
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂ +
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
+
∂
∂
∂
∂
∂
∂
∂
∂
=
xy d
xy f f xy d
f f xy
xy f f d d
f f
d
xy d
xy xy
d xy
xy d
d d
xy d
xy xy d
xy
xy d d
d
xy d
d d S dS S d
dS d
S dS S d
dS d τ
d d S dS dG G S dS dG G
S dS dG G S
dS dG G d
d d
d
τ τ τ τ τ τ τ τ τ τ
τ τ τ τ τ
τ τ
τ τ τ τ
τ τ
τ τ
τ τ
γ γ γ
τ γ
τ
γ τ γ
τ τ
τ
0 0 0
0
0 0 0
0
(3.3-14)
上式右辺の係数行列を、第一項から順に、K1、K2、K3と書くと、上式は、次のようになる。
+
+
=
xy d xy
d xy
d xy
d
d K d d
K d d
K d d
d
τ τ τ
τ γ
γ τ
τ
3 2
1 (3.3-15)
これより、次式を得る。
−
−
=
−
xy d xy
d
d K d K K d I
d
γ γ τ
τ
1 1 3
2 )
( (3.3-16)
あるいは、
=
xy d xy
d
d K d d
d
γ γ τ
τ (3.3-17)
ここに、
1 1 3
2 )
(I K K K
K = − − − (3.3-18)
なお、(3.3-14)式の各係数行列の成分の算定法については後で示す。特に、各時間ステップ 内で、状態変数Sも一定とした従来の定式化では、K2=K3=0であり、従って、K=K1となる。
一方、過剰間隙水圧モデルを適用した場合の、平均有効応力－平均ひずみ関係は、次式 で与えられる。
mK
vp v '
y ' ' x
m B
2
+ −
−
− + =
≡ 1
1
)
(ε ε
σ
σ σ (3.3-19)
第 3 章 地盤・構造物系の数値解析法
ここに、
K
K
m m ma K ma
m K B
−
−
=
1 1
) '
1 (
σ
(3.3-20)
また、εvは、全体積ひずみ（弾性体積ひずみ＋塑性体積ひずみ）であり、これには、先行 する各解析フェーズにおける体積ひずみ（但し弾性体積ひずみ）の累積も含み、次式で表 される。
x y v ≡ε +ε
ε (3.3-21)
εvpは、現解析フェーズで生じた塑性体積ひずみで、状態変数Sの関数である。
)
vp(S
vp ε
ε = (3.3-22)
なお、最終解析フェーズでのみ液状化を考慮することが出来るものとし、先行する各解析 フェーズにおける塑性体積ひずみの発生は考えていない。
ここで、偏差応力と同様、σm'の増分を求める。
) (
) (
) 1 )(
( 1
' '
1 '
vp v m
ma ma m
vp v m m vp v K m
d d K
d m d
B d
K
K K
ε σ ε
σ
ε ε ε
ε σ
−
=
− +
− −
= −
(3.3-23)
特に、(3.3-17)式を利用して、dεvpは、次のように変形出来る。
=
∂
∂
∂
= ∂
∂
∂
∂
= ∂
xy d
xy d xy
d vp
xy d xy d
vp vp
d K d
d K d S S
dS d
d S d S
dS d d
γ γ
γ γ τ
τ ε
τ τ τ τ ε ε
4
(3.3-24)
なお、
S K S
dS K d
xy d
vp
∂
∂
∂
= ∂
τ τ ε
4 (3.3-25)
従って、(3.3-23)式は、次のように書ける。
)
( 4
'
−
=
xy d v
m d
K d d a
d γ
ε γ
σ (3.3-26)
但し、