JP2007263623A - 地震応答解析方法 - Google Patents

Abstract

【課題】簡便で斉一性のある解が得られると共にサイクリックモビリティの影響を考慮した地震応答解析が可能な地震応答解析方法を提供する。
【解決手段】平均有効応力が変相線を越える前は、３次元応力空間における除荷点からの増分相当応力に基づいて累積損傷度を逐次計算する（ステップ２０８〜２１８）。変相線を越えた後は、変相線を越えた載荷状態か否かを判断し、変相線を越えた載荷状態の場合とそうでない場合とで異なる方法により接線剛性や平均有効応力の変動分を求める（ステップ２２０〜２３６）。
【選択図】図３

Description

本発明は、地震応答解析方法に係り、特に、地盤の液状化を考慮した地震応答解析方法に関する。

従来、地震時の地盤の液状化を考慮した時々刻々の地震応答の解析方法が提案されている。このような地盤の液状化を考慮した地震応答解析を行う際、液状化特性を示す実際の試験結果の傾向を再現するために、地震応答解析の数学モデルのパラメータを試行錯誤的に設定して計算し、この計算結果と試験結果との一致度を確認することを繰り返すことにより、数学モデルの最適なパラメータを確定する技術がある（例えば、特許文献１等参照）。

また、液状化特性を示す試験結果そのもの又はその回帰結果を数学モデルへの入力データとして用い、累積損傷度により液状化の影響を考慮した地震応答解析を行う技術がある（例えば、特許文献２等参照）。この方法によれば、試験結果との傾向の一致度が高く、また解析結果の斉一性も高い。
特許第３５６７４０７号公報 特開２００４−２４５６８４号公報

しかしながら、特許文献１に記載された技術では、多大な時間と高度な工学的経験が必要となり、実務に用いるには不便であると共に、斉一性のある解が得られず、また最適解であるという保証がない、という問題があった。

また、特許文献２に記載された技術では、完全な液状化に至るまでの挙動の解析方法であり、サイクリックモビリティ現象の影響が考慮されない、という問題があった。なお、サイクリックモビリティ現象とは、地盤の土等の密度が高い場合、一旦有効応力を喪失した後でもその後の繰返し載荷によって正のダイレイタンシーを示し、有効応力が回復するために破壊を生じない現象である。

本発明は、上記問題を解決すべく成されたものであり、簡便で斉一性のある解が得られると共にサイクリックモビリティの影響を考慮した地震応答解析が可能な地震応答解析方法を提供することを目的とする。

上記課題を解決するため、請求項１記載の発明は、解析対象モデルの質量行列情報、剛性行列情報、減衰行列情報、及び入力地震情報による運動方程式を用いて時々刻々と変化する前記解析対象モデルの変位情報を求め、前記変位情報に基づいて、時々刻々と変化する少なくとも応力に関する応答値の時間特性を求める地震応答解析方法において、前記運動方程式を時間積分することにより求めた変位情報に基づいて応力を求め、求めた応力に基づいてサイクリックモビリティ状態か否かを判断し、前記解析対象モデルの剛性を、サイクリックモビリティ状態になる前と、サイクリックモビリティ状態となった場合又は過去にサイクリックモビリティ状態となったことがある場合とで、異なる規則により求め、求めた剛性により前記剛性情報を更新する処理を、予め定めた所定時間毎に繰り返すことにより前記応答値の時間特性を求めることを特徴とする。

この発明によれば、解析対象モデルの剛性を、サイクリックモビリティ状態になる前と、サイクリックモビリティ状態となった場合又は過去にサイクリックモビリティ状態となったことがある場合とで、異なる規則により求めるため、サイクリックモビリティの影響を考慮した地震応答解析が可能となる。

なお、請求項２に記載したように、前記サイクリックモビリティ状態となった場合又は過去にサイクリックモビリティ状態となったことがある場合は、前記求めた応力が変相線を越えた載荷状態であるか否かを判断し、前記解析対象モデルの剛性を、変相線を越えた載荷状態である場合と、そうでない場合とで、異なる規則により求めるようにしてもよい。

具体的には、例えば請求項３に記載したように、前記求めた応力が変相線を越えた載荷状態である場合は、前記変相線を越えるときの接線剛性と前記変相線上の平均有効応力とに基づいて、接線剛性を求め、求めた接線剛性に基づいて前記応力の増分を求め、求めた応力の増分に基づいて、当該増分が予め定めた破壊線に沿って増加するように定められた所定式により平均有効応力の増分を求め、求めた平均有効応力の増分に基づいて応力を求めることができる。

また、請求項４に記載したように、前記求めた応力が変相線を越えた載荷状態でない場合は、ひずみと応力との関係を示す骨格曲線に基づいて接線剛性を求め、求めた接線剛性に基づいて前記応力の増分を求め、求めた応力の増分、前記応力の予め定めた半波開始点における平均有効応力、及び前記変相線上の平均有効応力に基づいて、平均有効応力の増分を求め、求めた平均有効応力の増分に基づいて応力を求めることができる。なお、請求項７に記載したように、前記半波開始点は、前記応力の除荷点としてもよい。

また、請求項５に記載したように、前記サイクリックモビリティ状態になる前は、ひずみと応力との関係を示す骨格曲線に基づいて接線剛性を求め、求めた接線剛性に基づいて前記応力の増分を求め、求めた応力の増分に基づいて、前記応力の予め定めた半波開始点からの累積損傷度の増分を求め、求めた累積損傷度の増分に前記半波開始点における累積損傷度を加算することにより累積損傷度を求め、求めた累積損傷度に基づいて過剰間隙水圧比を求め、求めた過剰間隙水圧比に基づいて平均有効応力の増分を求め、求めた平均有効応力の増分に基づいて応力を求めることができる。

この場合、請求項６に記載したように、前記累積損傷度は、３次元応力空間における相当応力の増分に基づいて求めるようにしてもよい。このように、３次元応力空間における相当応力の増分に基づいて累積損傷度を計算することにより、簡便かつ精度良く地震の３次元挙動を解析することができる。

以上説明したように、本発明によれば、簡便で斉一性のある解が得られると共にサイクリックモビリティの影響を考慮した地震応答解析が可能になる、という効果を有する。

以下、図面を参照して本発明の実施形態について説明する。

図１には、地震応答解析装置１０が示されている。地震応答解析装置１０は、操作部１２、記憶部１４、演算部１６、表示部１８で構成されている。

操作部１２は、オペレータが表示部１８に表示されたメニューに従って所望の解析モデルについての地震応答解析を演算部１６に実行させるための指示や必要なパラメータを指定するためのものである。記憶部１４は、演算部１６において実行される制御ルーチンのプログラムや、様々な解析モデルのパラメータ、地震応答解析に必要な各種演算式が記憶されている。また、記憶部１４には、演算部１６による地震応答解析の解析結果が格納される。

演算部１６は、操作部１２からの指示に従って記憶部１４から必要なデータを読み出して地震応答解析を行うと共に、出力結果を記憶部１４へ記憶すると共に、表示部１８へ出力する。演算部１６では、サイクリックモビリティを考慮した地震応答解析を実行する。

地震応答解析は、次式で示される運動方程式を時間積分することにより行う。

なお、本実施の形態では、時間積分には、一般によく知られた動的非線形問題の時間積分法を用いることができる。本実施の形態では、一例として、Ｎｅｗｍａｒｋ−β法による増分解法を用いた。

次に、本実施形態の作用として、演算部１６で実行される地震応答解析の制御ルーチンについて図２に示すフローチャートを参照して説明する。

まず、ステップ１００では、オペレータが表示部１８に表示されたメニューに従って操作部１２を操作し、地震応答解析を行うべき解析モデルのパラメータを指定すると、演算部１６では、指定された解析モデルに関するデータを記憶部１４から読み込む。

この解析モデルに関するデータには、例えば解析対象の地層、材料定数、時間積分定数等がある。

次のステップ１０２では、以下の地震応答を求める演算において用いる各種データの初期値の設定を行う。このデータには、例えば加速度、速度、変位、応力、ひずみ、上記（１）式における質量行列Ｍ、剛性行列Ｋ、減衰行列Ｃなどがある。演算部１６では、これらのデータに初期値を設定する。

次に、ステップ１０４において、ｔ＋Δｔ時刻の応答変位、応答速度、応答加速度の予測値ｕ^*を下記（２）式により算出すると共に、外力Ｆ（ｔ＋Δｔ）、粘性力Ｆ_C、慣性力Ｆ_Mの算出を行う。外力は上記（１）式の右辺の算出を（ｔ＋Δｔ）時刻について行うことにより得られる。粘性力Ｆ_Cは下記（３）式により、慣性力Ｆ_Mは下記（４）式により算出することができる。

ｕ^*＝Ａｕ（ｔ）＋ＢΛ …（２）
ここで、Ａは予測オペレーター行列、Ｂは修正オペレーター行列、Λは未知数である修正値ベクトルである。

そして、次のステップ１０６で、図３に示すような処理により内力Ｒａを求める。この処理では、サイクリックモビリティを考慮した応力を計算し、この応力から内力Ｒａを求める。以下、図３に示す内力Ｒａの計算について説明する。

まず、ステップ２００では、まず増分相当ひずみを計算する。現在のステップのひずみε_(n)は、変位から次式で求めることができる。なお、添え字の（ｎ）は現在のステップを表す。

ただし、Ｂは変位ひずみ関係マトリックス、Ｕ^e _(n)は要素の節点並びに定義された変位ベクトルである。なお、添え字のｅは要素の代表点の変位であることを表す。また、１節点での変位成分は、Ｕ^T＝｛ｕ ｖ ｗ｝とする。添え字のＴは転置行列であることを表す。

また、ひずみベクトルは下記並びで定義する。

増分ひずみΔε（ｄε）は、現在のｎステップと直前の（ｎ−１）ステップのひずみから次式で求める。
ｄε＝ε_(n)−ε_(n-1) ・・・（７）
次に、ステップ２０２で、平均有効応力σ’_mの経路（以下、有効応力経路という）が過去に変相線を越えているか否かを判断する。図４（Ｂ）に示すように、変相線２０は、平均有効応力σ’_mを横軸、せん断応力ηを縦軸にとった二次元空間において予め定められた境界線である。有効応力経路２２がこれを越えると、すなわち破壊線２４側の領域に入ると、サイクリックモビリティが発生する。変相線２０、破壊線２４は、試験結果の有効応力経路から予め定める。

そして、有効応力経路２２が過去に変相線２０を一度でも越えている場合には、ステップ２２０へ移行し、有効応力経路が過去に一度も変相線２０を越えていない場合には、ステップ２０６へ移行する。

ステップ２０６では、今回有効応力経路２２が変相線２０を越えたか否か、すなわち初めてサイクリックモビリティ状態に入ったか否かを判断し、有効応力経路２２が変相線２０を越えていない場合はステップ２０８へ移行し、越えている場合はステップ２２０へ移行する。

ステップ２０８では、せん断による剛性低下を考慮した応力のせん断成分を図５に示す処理により計算する。ここでは、３次元応力空間の各方向成分に共通なせん断弾性係数の接線勾配を求め、破壊曲面とその内側に降伏曲面を有するモデルで計算を行うと共に、有効応力に対応させるために無次元化した応力ひずみ関係で計算する。また、３次元に対応するために、π平面上に降伏曲面を定義する。

まず、ステップ３００において、無次元化ひずみξを計算する。無次元化ひずみξは、上記（５）式で求められるひずみεから、次式で求めることができる。

ここで、添え字のｉ，ｊは自由度を表す。また、ｅ_(n)は偏差ひずみであり、次式で表される。

δ_ijはクロネッカーのデルタ（ｉ≠ｊのときは０，ｉ＝ｊのときは１）である。

また、Ｇ_0(n-1)は平均有効応力に依存した初期せん断弾性係数であり、次式で表される。

また、σ_refは参照応力を、σ’_m(n-1)は（ｎ−１）ステップにおける平均有効応力を表す。Ｇ_refは参照せん断弾性係数を表し、σ’_m＝σ_refのときの初期せん断弾性係数を表す。ｍはパラメータ（一例として０．５）を表す。τ_max(n-1)は破壊時のせん断応力を表し、次式で表される。
τ_max(n-1)＝σ’_m(n-1)ｓｉｎφ＋ｃ ・・・(１１)
ここで、無次元化相当ひずみξを次式で定義する。

ここで、ｅ（ｎ）は相当ひずみであり、次式で表される。

無次元化ひずみ増分ｄｅは、次式で求めることができるが、これに限らず他の方法で求めることもできる。

次に、ステップ３０２で除荷点か否かを判定する。まず、除荷判定をするために、仮増分応力を求める。除荷判定は無次元化応力ηで行う。

無次元化応力ηは、次式で表される。

ここで、ｓ_(n-1)は偏差応力であり、次式で表される。

ここで、δ_ijはクロネッカーのデルタを表す。

また、無次元化相当応力η_(n-1)は次式で表される。
η_(n-1)＝σ_e(n-1)／τ_max(n-1) ・・・(１８)
ここで、σ_e(n-1)は相当応力であり、次式で表される。

なお、仮増分応力は、線形の無次元化増分応力を用い、次式で表される。
ｄη^e＝ｇ₀ｄξ ・・・(２０)
ここで、ｇ０は初期無次元化剛性を表し、次式で表される。

次に、除荷判定について具体的に説明する。

（ｎ−１）ステップの無次元化応力は、降伏曲面上にあるとする。図６に降伏曲面と応力の定義を示した。

現在の降伏曲面の大きさは、次式で求められる。

そして、現在の降伏曲面と無次元化増分応力ｄη^eを加えた後の降伏曲面の大きさを比較し、小さくなる点を除荷点とする。すなわち、次式を満足した場合に（ｎ−１）ステップ目が除荷点であると判定する。

そして、除荷点であると判定した場合には、ステップ３０４で、除荷点の無次元化応力・ひずみを更新する。無次元化応力が除荷した場合は、（ｎ−１）ステップの偏差応力、偏差ひずみ、除荷曲面の中心点、平均有効応力、破壊時のせん断応力を新たな除荷点情報として記憶する。すなわち、除荷点数ｉを１インクリメントし、次式により各パラメータを更新する。

ここで、ｓ^Ri _ijはｉ番目の除荷点における無次元化応力、ｅ^Ri _ijはｉ番目の除荷点における無次元化ひずみ、η^ci _ijはｉ番目の除荷点における除荷曲面の中心点、τ^Ri _maxはｉ番目の除荷点における破壊時のせん断応力、σ’^Ri _mはｉ番目の除荷点における平均有効応力である。

ステップ３０６では、除荷点からの無次元化相当ひずみξ^L _(n)を次式により求める。なお、添え字のＬは除荷点からの距離であることを表す。

ここで、ｅ^L _(n)は除荷点からの相当ひずみであり、次式で表される。

ここで、ｅ_(n)ijはｎステップの偏差ひずみ、ｅ^Ri _ijはｉ番目の除荷点の偏差ひずみを表す。

ステップ３０８では、骨格曲線／履歴曲線による接線剛性を計算する。骨格曲線/履歴曲線は、次式で示すように無次元化相当応力および無次元化相当ひずみの関係で定義される。

η＝Ｈ（ξ） ・・・（３３）
履歴曲線は、次式で示すように、除荷点からの距離の半分η^L／２及びξ^L／２が骨格曲線上にあるとして、Ｍａｓｉｎｇ則に対応させる。図７に骨格曲線及び履歴曲線の一例を示した。

ここで、η^L _(n)はｎステップの除荷点からの無次元化相当応力である。

そして、骨格曲線（上記３３式）及び無次元化ひずみ（上記３１式）より、無次元化接線剛性ｇを求める。

骨格曲線は、図８（ａ）に示すような予め与えられた動的変形特性としての剛性低下率Ｇ／Ｇ₀とせん断ひずみγとの関係を用いて作成する。一方、履歴曲線は、予め与えられた減衰率ｈとせん断ひずみγとの関係における減衰率ｈが等しくなるような仮想の骨格曲線を作成し、これにＭａｓｉｎｇ則を適用する。なお、これらの関係は、例えばテーブルデータ（図８（ａ）の○点のデータ）として記憶部１４に予め記憶される。

骨格曲線の計算では、上記の動的変形特性から無次元化ひずみξ_iと無次元化剛性ｇ_iとの関係に変更する。ここで、添字のｉはｉ番目のデータであることを示す。

まず、無次元化剛性ｇ_iの計算方法を示し、次に無次元化ひずみξ_iの計算方法を示す。

無次元化剛性ｇ_iは、剛性低下率Ｇ／Ｇ₀とせん断ひずみγとの関係を表すテーブルデータに基づいて求める。各区間（上記○点と○点との間）の接線剛性Ｇ_Tiは次式で求められる。

この接線剛性Ｇ_Tiから、各区間の無次元化剛性ｇ_iを次式で求める。

ただし、テーブルデータの１つ目の点における無次元化剛性ｇ_iは次式で求める。

上記（３７）式を図で表したものを図９に示す。同図に示すように、各区間では接線剛性が一定値となる。

次に、無次元化ひずみξ_iの求め方について説明する。

テーブルデータのひずみデータは、次式により無次元化ひずみに変換する。

なお、ｈ〜γ曲線のひずみも同様に無次元化ひずみに変換する。

次に、履歴曲線の計算について説明する。履歴曲線は、次式に示す双曲線モデルを用いて定義する。

ここで、Ｇ^R ₀はｈ〜γ曲線を満足するように次式により求める。そして、そのときのひずみを除荷点のひずみγ^Rとする。

ただし、χ^Rはｈより次式で表される。

無次元化接線剛性ｇは、双曲線モデルを用いて次式で求められる。

ただし、ξ^Hは除荷点からの無次元化ひずみであり、次式で求められる。

次のステップ３１０では、無次元化応力を求める。具体的には、無次元化接線剛性ｇより次式によって無次元化増分応力ｄηを求める。
ｄη＝ｇｄξ ・・・(４４)
そして、無次元化増分応力ｄηから偏差応力増分ｄｓを次式により求める。
ｄｓ＝ｄη・τ_max(n-1) ・・・(４５)
次に、無次元化応力η_(n)及び除荷点からの無次元化応力η^L _(n)を次式により更新する。
η_(n)＝η_(n-1)＋ｄη ・・・(４６)
η^L _(n)＝η_(n)−η^Ri ・・・(４７)
次のステップ３１２では、（ｎ−１）ステップの降伏曲面が直前の除荷曲面よりも大きいか否かを判断する。そして、（ｎ−１）ステップの降伏曲面が直前の除荷曲面よりも大きい場合はステップ３１４へ移行し、そうでない場合は本ルーチンを終了する。具体的には、次式が成り立つか否かを判断する。なお、この判断の概念を図１０に示す。

ステップ３１４では、（ｎ−１）ステップの降伏曲面を削除し、直前の除荷曲面の情報に置換する。具体的には、次のように直前の除荷曲面の情報を取り出す。まずｉ番目の除荷曲面データをクリアし、除荷点数ｉを１デクリメントし、ｉ＝０のときは骨格曲線に戻り、原点からの計算とし、以降は更新されたｉ番目の除荷点の情報を用いる。この処理の後は、ステップ３０６へ戻って再度接線剛性を計算し、無次元化応力を求め直す。

図３へ戻って、ステップ２１０では、応力の体積成分を求める。まず、過圧密を仮定し、膨潤線の膨張係数κを用いて体積弾性係数Ｋを求める。κと体積圧縮係数ｍ_vの関係は次式で表される。

ここで、σ’_mは拘束圧（ただし、ここでは圧縮が正）、ｅ₀は初期間隙比である。

なお、ｄσ’_m／σ’_mが十分に小さい場合は、次式を用いることが出来る。

また、体積圧縮指数ｍ_vと体積弾性係数Ｋとの関係は次式で表される。

上記（５０）式、（５１）式より、体積弾性係数Ｋは次式で求めることができる。

なお、Ｇは無次元化接線剛性ｇより次式により計算する。

Ｇ＝ｇＧ₀ ・・・(５３)
次に、応力の体積成分を次式で求める。

ｄσ_ij＝ｄｓ_ij＋δＫｄε_v ・・・(５４)
ここで、ｄε_vは体積ひずみ増分であり、ｄε_v＝ｄε_kkである。

次のステップ２１２では、増分応力を求めると共に次式により応力を更新する。

σ_n＝σ_n-1＋ｄσ ・・・(５５)
次に、ステップ２１４において、累積損傷度Ｄを求める。

累積損傷度Ｄは、本実施形態では一例として除荷点から次の除荷点までを半波の区切りとし、除荷点からの増分相当応力に対して逐次計算する。除荷点の偏差応力Ｓ^R _ijからの偏差応力増分について相当応力σ^L _e(n)を次式により求める。

なお、除荷点からの相当応力σ^L _e(n)の１／２を半波のせん断応力とし、応力比Ｒは次式で計算することができる。

このＲを用いて、累積損傷度増分ΔＤを求める。すなわち、図１１に示すように液状化強度曲線から応力比Ｒに対応する繰り返し回数Ｎ_ifを求め、この繰り返し回数Ｎ_ifから、除荷点からの累積損傷度増分ΔＤを計算する。なお、累積損傷度Ｄは、前述した特許文献２と同様に求めることができる。すなわち、累積損傷度Ｄは、同図に示すように、半波開始点ｔ_S、すなわち、除荷点における累積損傷度をＤ_Rとして次式で表すことができる。

Ｄ＝Ｄ_R＋ΔＤ …（５８）
また、累積損傷度Ｄ_Rは、相当応力σ^L _e(n)が除荷点となる毎に更新される。すなわち、図１１の場合には、半波開始点ｔ_S、半波終了点ｔ_eの時点において、そのときの累積損傷度Ｄが累積損傷度Ｄ_Rに設定される。累積損傷度増分ΔＤは、半波開始点ｔ_s、半波終了点ｔ_eの時点においてゼロに初期化される。

ステップ２１６では、過剰間隙水圧比増分を求める。累積損傷度Ｄから過剰間隙水圧比ｒ_uを求める方法は、例えば従来のＳｅｅｄの次式による方法や、地盤調査より得られる実測値のテーブルデータを用いる方法を用いることができる。ここで、過剰間隙水圧比ｒ_uを求めた後、過剰間隙水圧比の増分Δｒ_uを計算する。

次のステップ２１８では、平均有効応力の変動ｄσ’_mを次式により求める。
ｄσ’_m＝−ｄｒ_u・｜σ’_m0｜ ・・・(６０)
一方、ステップ２０２で過去にサイクリックモビリティ状態に入ったことがあると判断された場合、２０６でサイクリックモビリティ状態に入っていると判断された場合には、ステップ２２０において、有効応力経路が変相線２０を越えた載荷状態か否かを判断する。すなわち、せん断応力τが前回のステップで計算したそれより増加しているか否かを判断する。

そして、載荷状態であると判断した場合には、ステップ２２２へ移行し、載荷状態でない、すなわち除荷状態であると判断した場合は、ステップ２３０へ移行する。

ステップ２２２では、図１２（ａ）、（ｂ）に示すように変相線２０を越えるときのせん断剛性（せん断弾性係数の接線勾配）Ｇ_CMと平均有効応力σ’_mとから接線せん断剛性Ｇ_Tを求める。

せん断剛性Ｇ_CMは、初期剛性Ｇ_m0と応力経路が変相線２０と交差した時の平均有効応力σ’_mpを用いた、例えば次式により求めることができるが、これに限られるものではない。

ここで、Ｒ_smは有効応力低下率であり、次式で表される。

なお、上記（６１）式は、有効応力低下率Ｒ_smと接線剛性低下率（Ｇ_CM／Ｇ_m0）との関係を示す試験結果から導き出される。

そして、接線せん断剛性Ｇ_Tは、接線剛性Ｇ_CMを初期値として平均有効応力σ’_mpからの平均有効応力比（σ’_m(n-1)−σ’_mp）／σ’_m0により次式で表される。

ここで、ｆ_b（γ_max）は接線勾配増加率を表す関数であり、最大せん断ひずみγ_maxを用いて次式で表される。

ここで、ａ，ｂ，ｃは所定のパラメータであり、ｂはｆ_b（γ_max）の最小値、ａ，ｂは試験結果のデータを用いた回帰分析により求めることができる。

なお、上記（６３）式は、平均有効応力σ’_mpからの平均有効応力比と接線剛性低下率との関係を表す試験結果から導き出され、上記（６４）式は、最大せん断ひずみγ_maxと接線勾配増加率ｆｂとの関係を示す試験結果から導き出される。

次に、ステップ２２４において、体積弾性係数Ｋを算出する。これは、サイクリックモビリティ状態に入る前と同様に、すなわちステップ２１０と同様に（５２）式により求めることができる。

次に、ステップ２２６において、増分応力の計算を行うと共に応力を更新する。具体的には、まず、３次元の応力ひずみ関係マトリックスＤを組み立て、増分応力計算に用いる。Ｄマトリックスを下式に示す。

増分応力ｄσ、応力σ_nは次式により求めることができる。
ｄσ＝Ｄｄε ・・・(６６)
σ_n＝σ_n-1＋ｄσ ・・・(６７)
次に、ステップ２２８において、サイクリックモビリティによる平均有効応力の増分Δσ’_m（ｄσ’_m）を次式により求める。これにより、破壊線に沿って平均有効応力が増加する。

ここで、σ_eは相当応力、τ_maxは破壊時のせん断応力（τ_max＝σ’_m（ｓｉｎφ＋ｃｃｏｓφ）、Ｍ₀は変相線の破壊線との比（＝ｔａｎφ_p／ｔａｎφ_f）、φ_fは破壊角、φ_pは変相角、ｄσ_eは増分相当応力を示す。

一方、ステップ２２０で載荷状態ではなく除荷状態であると判断した場合は、ステップ２３０において、せん断ひずみγによる剛性低下Ｇ／Ｇ₀と平均有効応力σ’_mとから接線せん断剛性Ｇ_Tを求める。

具体的には、除荷点からの無次元化ひずみξの１／２を用いて骨格曲線／履歴曲線Ｈ（ξ）で無次元化剛性ｇを求め、せん断弾性係数の接線勾配Ｇ_Tを計算する。ただし、ここで用いる初期剛性Ｇ_0Uは、平均有効応力σ’_mに比例して低下するものとする。

Ｇ_T＝ｇＧ_0U ・・・（６９）
ここで、初期剛性Ｇ_0Uは次式で表される。

なお、Ｇ_m0は初期平均有効応力σ’_m0のときの初期剛性である。

Ｇ_Tが図４（ａ）に示す最小剛性（変相線２０上の接線剛性Ｇ_CM）より小さい時は、Ｇ_CMを用いる。Ｇ_CMは、上記（６１）式により求めることができる。

また、図１３（ａ），（ｂ）に示すように、除荷点が変相線を越えていない場合も同様に接線勾配を求める。

次に、ステップ２３２において、体積弾性係数Ｋを算出する。これは、サイクリックモビリティ状態に入る前と同様に、すなわちステップ２１０と同様に（５２）式により求めることができる。

次に、ステップ２３４において、ステップ２２６と同様に増分応力の計算を行うと共に応力を更新する。

ステップ２３６では、サイクリックモビリティによる平均有効応力増分Δσ’_mを求める。平均有効応力増分Δσ’_m（ｄσ’_m）は、次式で表される。この式は、経験式であり、下記（７４）式を微分して増分形式としたものである。下記（７４）式は、平均有効応力とせん断応力との関係を示す中空ねじり試験の応力経路を変相線上で１になるようにノーマライズした除荷点からの応力比と平均有効応力低下率の関係を表す。

ここで、σ’_mRは除荷点の平均有効応力を、σ’_mpは変相線上の平均有効応力を示す。また、η^Lは除荷点からの応力比であり、次式で表される。

ここで、η^Rは除荷点上の応力比を、η^p（＝ｔａｎφ_p：φ_pは変相角）は変相線上の応力比を示す。ｄηは応力比増分であり、次式で表される。

なお、ここでの応力比の定義はη＝τ／｜σ’_m｜とする。

除荷点からの平均有効応力低下率（σ’_m−σ’_mR）／（σ’_mp−σ’_mR）は次式で表される。

なお、図１４（ａ），（ｂ）に示すように、変相線２０を越える前に応力が反転する場合も、その点を除荷点として上記（７１）式を適用し、平均有効応力の変動分を求める。

変相線２０上の平均有効応力σ’_mpはＮ／Ｎ_f−ｒ_u関係より求める。Ｎ／Ｎ_f−ｒ_u関係は図１５に示すように、サイクリックモビリティで通るｒ_uの最大値の包絡線を定義したものである。

次のステップ２３８では、平均有効応力の変動分ｄσ’_mを考慮した応力を次式により求める。
σ_ij＝σ’_ij＋δ_ijｄσ’_m ・・・(７５)
ステップ２４０では、求めた応力σから内力Ｒａを次式により求める。

ここで、Ｖは対象とする体積であり、Ｎは節点（計算位置）への変換関数である。

図２に戻って、ステップ１０８では、外部から作用する外力Ｆ、粘性力Ｆ_C、慣性力Ｆ_Mなどからなる外力項とステップ１０６で得られた内力Ｒａとの残差力（不釣り合い力）ΔＲを求める。

次に、ステップ１１０において、上記剛性を用いて要素行列の算出とこれを用いた全体剛性の算出を行う。

次に、ステップ１１２において、動的解析用全体行列Ｋ^*を算出する。動的解析用全体行列Ｋ^*は次式で示される。
Ｋ^*＝Ｍ＋γΔｔＣ＋βΔｔ²Ｋ …（７７）
ここで、β、γは、Ｎｅｗｍａｒｋ−β法における係数（一定）である。また、Ｋは接線剛性である。

次に、ステップ１１４において、残差力より（ｔ＋Δｔ）時刻の修正値Λを下記（７８）式により算出し、この修正値Λにより予測値ｕ^*を修正し、（ｔ＋Δｔ）時刻の応答変位、応答速度、応答加速度を求め、これにより外力Ｆ、慣性力Ｆ_M及び粘性力Ｆ_Cを算出する。

Δｔは計算時間間隔、すなわち時刻歴ループの実行間隔である。なお、Δｔは一定の値でもよいし、計算毎に変化させてもよい。

次に、ステップ１１６において、ステップ１０６と同様に内力Ｒａを算出する。次に、ステップ１１８において、ステップ１０８と同様に外力と内力の差（残差力）ΔＲを算出する。

次に、ステップ１２０において、この残差力ΔＲが許容値以下であるか否かを判断し、許容値を越えていれば、ステップ１１０から反復して計算する。この残差力ΔＲが許容値以下であれば、次のステップ１２２へ進む。

次に、ステップ１２２において、変位、速度、加速度等を更新すると共に、算出した各応答値を表示部１８や記憶部１４へ出力する。出力値としては、変位、速度、加速度、ひずみ、応力等がある。

このようにして、地震時間全体についてステップ１０４からステップ１２２までの処理（時刻歴ループ）を行い、各時刻毎の応答値を算出する。

以上説明したように、本実施形態では、除荷点から次の除荷点までを半波の区切りとし、３次元応力空間における除荷点からの増分相当応力に基づいて累積損傷度を逐次計算するので、簡便かつ精度良く地震の３次元挙動を解析することができる。

また、有効応力経路が変相線を越える前と後とで異なる手法により平均有効応力の変動を求めるので、サイクリックモビリティの影響が反映された地震応答解析を行うことができる。

なお、本実施形態では、３次元応力空間における除荷点からの増分相当応力に基づいて累積損傷度を逐次計算する場合について説明したが、従来（例えば前記特許文献２記載の技術）のように１次元の応力、すなわちせん断応力に基づいて累積損傷度を逐次計算するようにしてもよい。この場合、増分ひずみの計算において、上記の１節点での変位成分はｕ^T＝｛ｕ｝として、上記（６）式は下記（７９）式に、上記（７）式は下記（８０）式に置き換える。
ε^T＝｛γ｝ ・・・（７９）
ｄγ＝γ_(n)−γ_(n-1) ・・・（８０）
なお、γはせん断ひずみである。

また、累積損傷度の計算において、上記（５６）式はτ（せん断応力）に、上記（５７）式は下記（８１）式に置き換える。

また、増分応力の計算、応力の更新処理において、上記（６５）式は下記（８２）式に、上記（６６）式は下記（８３）式に、上記（６７）式は下記（８４）式に置き換える。
Ｄ＝［Ｇ_T］ ・・・（８２）
ｄτ＝Ｇ_Tｄγ ・・・（８３）
τ_n＝τ_(n-1)＋ｄτ ・・・（８４）
また、本実施形態では、除荷点から次の除荷点までを半波の区切りとした場合について説明したが、これに限らず、半波の開始点及び終了点は任意の位置に設定することができる。

地震応答解析装置の概略ブロック図である。 演算部において実行される制御ルーチンのフローチャートである。 内力の計算処理のフローチャートである。 （ａ）は接線剛性について説明するためのせん断ひずみとせん断応力との関係を示す図、（ｂ）は有効応力経路について説明するための図である。 剛性低下の計算処理のフローチャートである。 降伏曲面と除荷曲面について説明するための図である。 骨格曲線及び履歴曲線について説明するための図である。 （ａ）は剛性低下率とせん断ひずみ、減衰との関係を示す図、（ｂ）はせん断ひずみとせん断応力との関係を示す図である。 各区間の接線剛性について説明するためのせん断ひずみとせん断応力との関係を示す図である。 降伏曲面が直前の除荷曲面よりも大きいか否かの判断について説明するための図である。 相当応力から液状化強度曲線を用いて累積損傷度を求める場合について説明するための図である。 （ａ）は接線剛性について説明するためのせん断ひずみとせん断応力との関係を示す図、（ｂ）は有効応力経路について説明するための図である。 （ａ）は除荷点が変相線を越えていない場合において接線剛性を求める場合について説明するためのせん断ひずみとせん断応力との関係を示す図、（ｂ）は有効応力経路について説明するための図である。 （ａ）、（ｂ）は変相線を越える前に応力が反転する場合において平均有効応力の変動分を求める場合について説明するための有効応力経路を示す図である。 Ｎ／Ｎ_f−ｒ_u関係を示す図である。

符号の説明

１０ 地震応答解析装置
１２ 操作部
１４ 記憶部
１６ 演算部
２０ 変相線
２２ 有効応力経路
２４ 破壊線

Claims (7)

1. 解析対象モデルの質量行列情報、剛性行列情報、減衰行列情報、及び入力地震情報による運動方程式を用いて時々刻々と変化する前記解析対象モデルの変位情報を求め、前記変位情報に基づいて、時々刻々と変化する少なくとも応力に関する応答値の時間特性を求める地震応答解析方法において、
   前記運動方程式を時間積分することにより求めた変位情報に基づいて応力を求め、
   求めた応力に基づいてサイクリックモビリティ状態か否かを判断し、
   前記解析対象モデルの剛性を、サイクリックモビリティ状態になる前と、サイクリックモビリティ状態となった場合又は過去にサイクリックモビリティ状態となったことがある場合とで、異なる規則により求め、
   求めた剛性により前記剛性情報を更新する処理を、
   予め定めた所定時間毎に繰り返すことにより前記応答値の時間特性を求める
   ことを特徴とする地震応答解析方法。
2. 前記サイクリックモビリティ状態となった場合又は過去にサイクリックモビリティ状態となったことがある場合は、
   前記求めた応力が変相線を越えた載荷状態であるか否かを判断し、
   前記解析対象モデルの剛性を、変相線を越えた載荷状態である場合と、そうでない場合とで、異なる規則により求める
   ことを特徴とする請求項１記載の地震応答解析方法。
3. 前記求めた応力が変相線を越えた載荷状態である場合は、
   前記変相線を越えるときの接線剛性と前記変相線上の平均有効応力とに基づいて、接線剛性を求め、
   求めた接線剛性に基づいて前記応力の増分を求め、
   求めた応力の増分に基づいて、当該増分が予め定めた破壊線に沿って増加するように定められた所定式により平均有効応力の増分を求め、
   求めた平均有効応力の増分に基づいて応力を求める
   ことを特徴とする請求項２記載の地震応答解析方法。
4. 前記求めた応力が変相線を越えた載荷状態でない場合は、
   ひずみと応力との関係を示す骨格曲線に基づいて接線剛性を求め、
   求めた接線剛性に基づいて前記応力の増分を求め、
   求めた応力の増分、前記応力の予め定めた半波開始点における平均有効応力、及び前記変相線上の平均有効応力に基づいて、平均有効応力の増分を求め、
   求めた平均有効応力の増分に基づいて応力を求める
   ことを特徴とする請求項２又は請求項３記載の地震応答解析方法。
5. 前記サイクリックモビリティ状態になる前は、
   ひずみと応力との関係を示す骨格曲線に基づいて接線剛性を求め、
   求めた接線剛性に基づいて前記応力の増分を求め、
   求めた応力の増分に基づいて、前記応力の予め定めた半波開始点からの累積損傷度の増分を求め、
   求めた累積損傷度の増分に前記半波開始点における累積損傷度を加算することにより累積損傷度を求め、
   求めた累積損傷度に基づいて過剰間隙水圧比を求め、
   求めた過剰間隙水圧比に基づいて平均有効応力の増分を求め、
   求めた平均有効応力の増分に基づいて応力を求める
   ことを特徴とする請求項１乃至請求項４の何れか１項に記載の地震応答解析方法。
6. 前記累積損傷度は、３次元応力空間における相当応力の増分に基づいて求めることを特徴とする請求項５記載の地震応答解析方法。
7. 前記半波開始点は、前記応力の除荷点であることを特徴とする請求項４乃至請求項６の何れか１項に記載の地震応答解析方法。

Images