(57)【特許請求の範囲】
【請求項１】
解析対象モデルの質量行列情報、剛性行列情報、減衰行列情報、及び入力地震情報による運動方程式を用いて時々刻々と変化する前記解析対象モデルの変位情報を求め、前記変位情報に基づいて、時々刻々と変化する少なくとも応力に関する応答値の時間特性を求める地震応答解析方法において、
前記運動方程式を時間積分することにより求めた変位情報に基づいて応力を求め、
求めた応力に基づいてサイクリックモビリティ状態か否かを判断し、
前記解析対象モデルの剛性を、サイクリックモビリティ状態になる前と、サイクリックモビリティ状態となった場合又は過去にサイクリックモビリティ状態となったことがある場合とで、異なる規則により求め、
求めた剛性により前記剛性情報を更新する処理を、
予め定めた所定時間毎に繰り返すことにより前記応答値の時間特性を求める
ことを特徴とする地震応答解析方法。
【請求項２】
前記サイクリックモビリティ状態となった場合又は過去にサイクリックモビリティ状態となったことがある場合は、
前記求めた応力が変相線を越えた載荷状態であるか否かを判断し、
前記解析対象モデルの剛性を、変相線を越えた載荷状態である場合と、そうでない場合とで、異なる規則により求める
ことを特徴とする請求項１記載の地震応答解析方法。
【請求項３】
前記求めた応力が変相線を越えた載荷状態である場合は、
前記変相線を越えるときの接線剛性と前記変相線上の平均有効応力とに基づいて、接線剛性を求め、
求めた接線剛性に基づいて前記応力の増分を求め、
求めた応力の増分に基づいて、当該増分が予め定めた破壊線に沿って増加するように定められた所定式により平均有効応力の増分を求め、
求めた平均有効応力の増分に基づいて応力を求める
ことを特徴とする請求項２記載の地震応答解析方法。
【請求項４】
前記求めた応力が変相線を越えた載荷状態でない場合は、
ひずみと応力との関係を示す骨格曲線に基づいて接線剛性を求め、
求めた接線剛性に基づいて前記応力の増分を求め、
求めた応力の増分、前記応力の予め定めた半波開始点における平均有効応力、及び前記変相線上の平均有効応力に基づいて、平均有効応力の増分を求め、
求めた平均有効応力の増分に基づいて応力を求める
ことを特徴とする請求項２又は請求項３記載の地震応答解析方法。
【請求項５】
前記サイクリックモビリティ状態になる前は、
ひずみと応力との関係を示す骨格曲線に基づいて接線剛性を求め、
求めた接線剛性に基づいて前記応力の増分を求め、
求めた応力の増分に基づいて、前記応力の予め定めた半波開始点からの累積損傷度の増分を求め、
求めた累積損傷度の増分に前記半波開始点における累積損傷度を加算することにより累積損傷度を求め、
求めた累積損傷度に基づいて過剰間隙水圧比を求め、
求めた過剰間隙水圧比に基づいて平均有効応力の増分を求め、
求めた平均有効応力の増分に基づいて応力を求める
ことを特徴とする請求項１乃至請求項４の何れか１項に記載の地震応答解析方法。
【請求項６】
前記累積損傷度は、３次元応力空間における相当応力の増分に基づいて求めることを特徴とする請求項５記載の地震応答解析方法。
【請求項７】
前記半波開始点は、前記応力の除荷点であることを特徴とする請求項４乃至請求項６の何れか１項に記載の地震応答解析方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、地震応答解析方法に係り、特に、地盤の液状化を考慮した地震応答解析方法に関する。
【背景技術】
【０００２】
従来、地震時の地盤の液状化を考慮した時々刻々の地震応答の解析方法が提案されている。このような地盤の液状化を考慮した地震応答解析を行う際、液状化特性を示す実際の試験結果の傾向を再現するために、地震応答解析の数学モデルのパラメータを試行錯誤的に設定して計算し、この計算結果と試験結果との一致度を確認することを繰り返すことにより、数学モデルの最適なパラメータを確定する技術がある（例えば、特許文献１等参照）。
【０００３】
また、液状化特性を示す試験結果そのもの又はその回帰結果を数学モデルへの入力データとして用い、累積損傷度により液状化の影響を考慮した地震応答解析を行う技術がある（例えば、特許文献２等参照）。この方法によれば、試験結果との傾向の一致度が高く、また解析結果の斉一性も高い。
【特許文献１】特許第３５６７４０７号公報
【特許文献２】特開２００４−２４５６８４号公報
【発明の開示】
【発明が解決しようとする課題】
【０００４】
しかしながら、特許文献１に記載された技術では、多大な時間と高度な工学的経験が必要となり、実務に用いるには不便であると共に、斉一性のある解が得られず、また最適解であるという保証がない、という問題があった。
【０００５】
また、特許文献２に記載された技術では、完全な液状化に至るまでの挙動の解析方法であり、サイクリックモビリティ現象の影響が考慮されない、という問題があった。なお、サイクリックモビリティ現象とは、地盤の土等の密度が高い場合、一旦有効応力を喪失した後でもその後の繰返し載荷によって正のダイレイタンシーを示し、有効応力が回復するために破壊を生じない現象である。
【０００６】
本発明は、上記問題を解決すべく成されたものであり、簡便で斉一性のある解が得られると共にサイクリックモビリティの影響を考慮した地震応答解析が可能な地震応答解析方法を提供することを目的とする。
【課題を解決するための手段】
【０００７】
上記課題を解決するため、請求項１記載の発明は、解析対象モデルの質量行列情報、剛性行列情報、減衰行列情報、及び入力地震情報による運動方程式を用いて時々刻々と変化する前記解析対象モデルの変位情報を求め、前記変位情報に基づいて、時々刻々と変化する少なくとも応力に関する応答値の時間特性を求める地震応答解析方法において、前記運動方程式を時間積分することにより求めた変位情報に基づいて応力を求め、求めた応力に基づいてサイクリックモビリティ状態か否かを判断し、前記解析対象モデルの剛性を、サイクリックモビリティ状態になる前と、サイクリックモビリティ状態となった場合又は過去にサイクリックモビリティ状態となったことがある場合とで、異なる規則により求め、求めた剛性により前記剛性情報を更新する処理を、予め定めた所定時間毎に繰り返すことにより前記応答値の時間特性を求めることを特徴とする。
【０００８】
この発明によれば、解析対象モデルの剛性を、サイクリックモビリティ状態になる前と、サイクリックモビリティ状態となった場合又は過去にサイクリックモビリティ状態となったことがある場合とで、異なる規則により求めるため、サイクリックモビリティの影響を考慮した地震応答解析が可能となる。
【０００９】
なお、請求項２に記載したように、前記サイクリックモビリティ状態となった場合又は過去にサイクリックモビリティ状態となったことがある場合は、前記求めた応力が変相線を越えた載荷状態であるか否かを判断し、前記解析対象モデルの剛性を、変相線を越えた載荷状態である場合と、そうでない場合とで、異なる規則により求めるようにしてもよい。
【００１０】
具体的には、例えば請求項３に記載したように、前記求めた応力が変相線を越えた載荷状態である場合は、前記変相線を越えるときの接線剛性と前記変相線上の平均有効応力とに基づいて、接線剛性を求め、求めた接線剛性に基づいて前記応力の増分を求め、求めた応力の増分に基づいて、当該増分が予め定めた破壊線に沿って増加するように定められた所定式により平均有効応力の増分を求め、求めた平均有効応力の増分に基づいて応力を求めることができる。
【００１１】
また、請求項４に記載したように、前記求めた応力が変相線を越えた載荷状態でない場合は、ひずみと応力との関係を示す骨格曲線に基づいて接線剛性を求め、求めた接線剛性に基づいて前記応力の増分を求め、求めた応力の増分、前記応力の予め定めた半波開始点における平均有効応力、及び前記変相線上の平均有効応力に基づいて、平均有効応力の増分を求め、求めた平均有効応力の増分に基づいて応力を求めることができる。なお、請求項７に記載したように、前記半波開始点は、前記応力の除荷点としてもよい。
【００１２】
また、請求項５に記載したように、前記サイクリックモビリティ状態になる前は、ひずみと応力との関係を示す骨格曲線に基づいて接線剛性を求め、求めた接線剛性に基づいて前記応力の増分を求め、求めた応力の増分に基づいて、前記応力の予め定めた半波開始点からの累積損傷度の増分を求め、求めた累積損傷度の増分に前記半波開始点における累積損傷度を加算することにより累積損傷度を求め、求めた累積損傷度に基づいて過剰間隙水圧比を求め、求めた過剰間隙水圧比に基づいて平均有効応力の増分を求め、求めた平均有効応力の増分に基づいて応力を求めることができる。
【００１３】
この場合、請求項６に記載したように、前記累積損傷度は、３次元応力空間における相当応力の増分に基づいて求めるようにしてもよい。このように、３次元応力空間における相当応力の増分に基づいて累積損傷度を計算することにより、簡便かつ精度良く地震の３次元挙動を解析することができる。
【発明の効果】
【００１４】
以上説明したように、本発明によれば、簡便で斉一性のある解が得られると共にサイクリックモビリティの影響を考慮した地震応答解析が可能になる、という効果を有する。
【発明を実施するための最良の形態】
【００１５】
以下、図面を参照して本発明の実施形態について説明する。
【００１６】
図１には、地震応答解析装置１０が示されている。地震応答解析装置１０は、操作部１２、記憶部１４、演算部１６、表示部１８で構成されている。
【００１７】
操作部１２は、オペレータが表示部１８に表示されたメニューに従って所望の解析モデルについての地震応答解析を演算部１６に実行させるための指示や必要なパラメータを指定するためのものである。記憶部１４は、演算部１６において実行される制御ルーチンのプログラムや、様々な解析モデルのパラメータ、地震応答解析に必要な各種演算式が記憶されている。また、記憶部１４には、演算部１６による地震応答解析の解析結果が格納される。
【００１８】
演算部１６は、操作部１２からの指示に従って記憶部１４から必要なデータを読み出して地震応答解析を行うと共に、出力結果を記憶部１４へ記憶すると共に、表示部１８へ出力する。演算部１６では、サイクリックモビリティを考慮した地震応答解析を実行する。
【００１９】
地震応答解析は、次式で示される運動方程式を時間積分することにより行う。
【００２０】
【数１】
0004513776.tif 000002
【００２１】
なお、本実施の形態では、時間積分には、一般によく知られた動的非線形問題の時間積分法を用いることができる。本実施の形態では、一例として、Ｎｅｗｍａｒｋ−β法による増分解法を用いた。
【００２２】
次に、本実施形態の作用として、演算部１６で実行される地震応答解析の制御ルーチンについて図２に示すフローチャートを参照して説明する。
【００２３】
まず、ステップ１００では、オペレータが表示部１８に表示されたメニューに従って操作部１２を操作し、地震応答解析を行うべき解析モデルのパラメータを指定すると、演算部１６では、指定された解析モデルに関するデータを記憶部１４から読み込む。
【００２４】
この解析モデルに関するデータには、例えば解析対象の地層、材料定数、時間積分定数等がある。
【００２５】
次のステップ１０２では、以下の地震応答を求める演算において用いる各種データの初期値の設定を行う。このデータには、例えば加速度、速度、変位、応力、ひずみ、上記（１）式における質量行列Ｍ、剛性行列Ｋ、減衰行列Ｃなどがある。演算部１６では、これらのデータに初期値を設定する。
【００２６】
次に、ステップ１０４において、ｔ＋Δｔ時刻の応答変位、応答速度、応答加速度の予測値ｕ*を下記（２）式により算出すると共に、外力Ｆ（ｔ＋Δｔ）、粘性力ＦC、慣性力ＦMの算出を行う。外力は上記（１）式の右辺の算出を（ｔ＋Δｔ）時刻について行うことにより得られる。粘性力ＦCは下記（３）式により、慣性力ＦMは下記（４）式により算出することができる。
【００２７】
ｕ*＝Ａｕ（ｔ）＋ＢΛ    …（２）
ここで、Ａは予測オペレーター行列、Ｂは修正オペレーター行列、Λは未知数である修正値ベクトルである。
【００２８】
【数２】
0004513776.tif 000003
【００２９】
【数３】
0004513776.tif 000004
【００３０】
そして、次のステップ１０６で、図３に示すような処理により内力Ｒａを求める。この処理では、サイクリックモビリティを考慮した応力を計算し、この応力から内力Ｒａを求める。以下、図３に示す内力Ｒａの計算について説明する。
【００３１】
まず、ステップ２００では、まず増分相当ひずみを計算する。現在のステップのひずみε(n)は、変位から次式で求めることができる。なお、添え字の（ｎ）は現在のステップを表す。
【００３２】
【数４】
0004513776.tif 000005
【００３３】
ただし、Ｂは変位ひずみ関係マトリックス、Ｕe(n)は要素の節点並びに定義された変位ベクトルである。なお、添え字のｅは要素の代表点の変位であることを表す。また、１節点での変位成分は、ＵT＝｛ｕ  ｖ  ｗ｝とする。添え字のＴは転置行列であることを表す。
【００３４】
また、ひずみベクトルは下記並びで定義する。
【００３５】
【数５】
0004513776.tif 000006
【００３６】
増分ひずみΔε（ｄε）は、現在のｎステップと直前の（ｎ−１）ステップのひずみから次式で求める。
ｄε＝ε(n)−ε(n-1)    ・・・（７）
次に、ステップ２０２で、平均有効応力σ’mの経路（以下、有効応力経路という）が過去に変相線を越えているか否かを判断する。図４（Ｂ）に示すように、変相線２０は、平均有効応力σ’mを横軸、せん断応力ηを縦軸にとった二次元空間において予め定められた境界線である。有効応力経路２２がこれを越えると、すなわち破壊線２４側の領域に入ると、サイクリックモビリティが発生する。変相線２０、破壊線２４は、試験結果の有効応力経路から予め定める。
【００３７】
そして、有効応力経路２２が過去に変相線２０を一度でも越えている場合には、ステップ２２０へ移行し、有効応力経路が過去に一度も変相線２０を越えていない場合には、ステップ２０６へ移行する。
【００３８】
ステップ２０６では、今回有効応力経路２２が変相線２０を越えたか否か、すなわち初めてサイクリックモビリティ状態に入ったか否かを判断し、有効応力経路２２が変相線２０を越えていない場合はステップ２０８へ移行し、越えている場合はステップ２２０へ移行する。
【００３９】
ステップ２０８では、せん断による剛性低下を考慮した応力のせん断成分を図５に示す処理により計算する。ここでは、３次元応力空間の各方向成分に共通なせん断弾性係数の接線勾配を求め、破壊曲面とその内側に降伏曲面を有するモデルで計算を行うと共に、有効応力に対応させるために無次元化した応力ひずみ関係で計算する。また、３次元に対応するために、π平面上に降伏曲面を定義する。
【００４０】
まず、ステップ３００において、無次元化ひずみξを計算する。無次元化ひずみξは、上記（５）式で求められるひずみεから、次式で求めることができる。
【００４１】
【数６】
0004513776.tif 000007
【００４２】
ここで、添え字のｉ，ｊは自由度を表す。また、ｅ(n)は偏差ひずみであり、次式で表される。
【００４３】
【数７】
0004513776.tif 000008
【００４４】
δijはクロネッカーのデルタ（ｉ≠ｊのときは０，ｉ＝ｊのときは１）である。
【００４５】
また、Ｇ0(n-1)は平均有効応力に依存した初期せん断弾性係数であり、次式で表される。
【００４６】
【数８】
0004513776.tif 000009
【００４７】
また、σrefは参照応力を、σ’m(n-1)は（ｎ−１）ステップにおける平均有効応力を表す。Ｇrefは参照せん断弾性係数を表し、σ’m＝σrefのときの初期せん断弾性係数を表す。ｍはパラメータ（一例として０．５）を表す。τmax(n-1)は破壊時のせん断応力を表し、次式で表される。
τmax(n-1)＝σ’m(n-1)ｓｉｎφ＋ｃ    ・・・(１１)
ここで、無次元化相当ひずみξを次式で定義する。
【００４８】
【数９】
0004513776.tif 000010
【００４９】
ここで、ｅ（ｎ）は相当ひずみであり、次式で表される。
【００５０】
【数１０】
0004513776.tif 000011
【００５１】
無次元化ひずみ増分ｄｅは、次式で求めることができるが、これに限らず他の方法で求めることもできる。
【００５２】
【数１１】
0004513776.tif 000012
【００５３】
次に、ステップ３０２で除荷点か否かを判定する。まず、除荷判定をするために、仮増分応力を求める。除荷判定は無次元化応力ηで行う。
【００５４】
無次元化応力ηは、次式で表される。
【００５５】
【数１２】
0004513776.tif 000013
【００５６】
ここで、ｓ(n-1)は偏差応力であり、次式で表される。
【００５７】
【数１３】
0004513776.tif 000014
【００５８】
ここで、δijはクロネッカーのデルタを表す。
【００５９】
また、無次元化相当応力η(n-1)は次式で表される。
η(n-1)＝σe(n-1)／τmax(n-1)      ・・・(１８)
ここで、σe(n-1)は相当応力であり、次式で表される。
【００６０】
【数１４】
0004513776.tif 000015
【００６１】
なお、仮増分応力は、線形の無次元化増分応力を用い、次式で表される。
ｄηe＝ｇ0ｄξ      ・・・(２０)
ここで、ｇ０は初期無次元化剛性を表し、次式で表される。
【００６２】
【数１５】
0004513776.tif 000016
【００６３】
次に、除荷判定について具体的に説明する。
【００６４】
（ｎ−１）ステップの無次元化応力は、降伏曲面上にあるとする。図６に降伏曲面と応力の定義を示した。
【００６５】
現在の降伏曲面の大きさは、次式で求められる。
【００６６】
【数１６】
0004513776.tif 000017
【００６７】
そして、現在の降伏曲面と無次元化増分応力ｄηeを加えた後の降伏曲面の大きさを比較し、小さくなる点を除荷点とする。すなわち、次式を満足した場合に（ｎ−１）ステップ目が除荷点であると判定する。
【００６８】
【数１７】
0004513776.tif 000018
【００６９】
【数１８】
0004513776.tif 000019
【００７０】
そして、除荷点であると判定した場合には、ステップ３０４で、除荷点の無次元化応力・ひずみを更新する。無次元化応力が除荷した場合は、（ｎ−１）ステップの偏差応力、偏差ひずみ、除荷曲面の中心点、平均有効応力、破壊時のせん断応力を新たな除荷点情報として記憶する。すなわち、除荷点数ｉを１インクリメントし、次式により各パラメータを更新する。
【００７１】
【数１９】
0004513776.tif 000020
【００７２】
ここで、ｓRiijはｉ番目の除荷点における無次元化応力、ｅRiijはｉ番目の除荷点における無次元化ひずみ、ηciijはｉ番目の除荷点における除荷曲面の中心点、τRimaxはｉ番目の除荷点における破壊時のせん断応力、σ’Rimはｉ番目の除荷点における平均有効応力である。
【００７３】
ステップ３０６では、除荷点からの無次元化相当ひずみξL(n)を次式により求める。なお、添え字のＬは除荷点からの距離であることを表す。
【００７４】
【数２０】
0004513776.tif 000021
【００７５】
ここで、ｅL(n)は除荷点からの相当ひずみであり、次式で表される。
【００７６】
【数２１】
0004513776.tif 000022
【００７７】
ここで、ｅ(n)ijはｎステップの偏差ひずみ、ｅRiijはｉ番目の除荷点の偏差ひずみを表す。
【００７８】
ステップ３０８では、骨格曲線／履歴曲線による接線剛性を計算する。骨格曲線/履歴曲線は、次式で示すように無次元化相当応力および無次元化相当ひずみの関係で定義される。
【００７９】
η＝Ｈ（ξ）      ・・・（３３）
履歴曲線は、次式で示すように、除荷点からの距離の半分ηL／２及びξL／２が骨格曲線上にあるとして、Ｍａｓｉｎｇ則に対応させる。図７に骨格曲線及び履歴曲線の一例を示した。
【００８０】
【数２２】
0004513776.tif 000023
【００８１】
ここで、ηL(n)はｎステップの除荷点からの無次元化相当応力である。
【００８２】
そして、骨格曲線（上記３３式）及び無次元化ひずみ（上記３１式）より、無次元化接線剛性ｇを求める。
【００８３】
骨格曲線は、図８（ａ）に示すような予め与えられた動的変形特性としての剛性低下率Ｇ／Ｇ0とせん断ひずみγとの関係を用いて作成する。一方、履歴曲線は、予め与えられた減衰率ｈとせん断ひずみγとの関係における減衰率ｈが等しくなるような仮想の骨格曲線を作成し、これにＭａｓｉｎｇ則を適用する。なお、これらの関係は、例えばテーブルデータ（図８（ａ）の○点のデータ）として記憶部１４に予め記憶される。
【００８４】
骨格曲線の計算では、上記の動的変形特性から無次元化ひずみξiと無次元化剛性ｇiとの関係に変更する。ここで、添字のｉはｉ番目のデータであることを示す。
【００８５】
まず、無次元化剛性ｇiの計算方法を示し、次に無次元化ひずみξiの計算方法を示す。
【００８６】
無次元化剛性ｇiは、剛性低下率Ｇ／Ｇ0とせん断ひずみγとの関係を表すテーブルデータに基づいて求める。各区間（上記○点と○点との間）の接線剛性ＧTiは次式で求められる。
【００８７】
【数２３】
0004513776.tif 000024
【００８８】
この接線剛性ＧTiから、各区間の無次元化剛性ｇiを次式で求める。
【００８９】
【数２４】
0004513776.tif 000025
【００９０】
ただし、テーブルデータの１つ目の点における無次元化剛性ｇiは次式で求める。
【００９１】
【数２５】
0004513776.tif 000026
【００９２】
上記（３７）式を図で表したものを図９に示す。同図に示すように、各区間では接線剛性が一定値となる。
【００９３】
次に、無次元化ひずみξiの求め方について説明する。
【００９４】
テーブルデータのひずみデータは、次式により無次元化ひずみに変換する。
【００９５】
【数２６】
0004513776.tif 000027
【００９６】
なお、ｈ〜γ曲線のひずみも同様に無次元化ひずみに変換する。
【００９７】
次に、履歴曲線の計算について説明する。履歴曲線は、次式に示す双曲線モデルを用いて定義する。
【００９８】
【数２７】
0004513776.tif 000028
【００９９】
ここで、ＧR0はｈ〜γ曲線を満足するように次式により求める。そして、そのときのひずみを除荷点のひずみγRとする。
【０１００】
【数２８】
0004513776.tif 000029
【０１０１】
ただし、χRはｈより次式で表される。
【０１０２】
【数２９】
0004513776.tif 000030
【０１０３】
無次元化接線剛性ｇは、双曲線モデルを用いて次式で求められる。
【０１０４】
【数３０】
0004513776.tif 000031
【０１０５】
ただし、ξHは除荷点からの無次元化ひずみであり、次式で求められる。
【０１０６】
【数３１】
0004513776.tif 000032
【０１０７】
次のステップ３１０では、無次元化応力を求める。具体的には、無次元化接線剛性ｇより次式によって無次元化増分応力ｄηを求める。
ｄη＝ｇｄξ    ・・・(４４)
そして、無次元化増分応力ｄηから偏差応力増分ｄｓを次式により求める。
ｄｓ＝ｄη・τmax(n-1)      ・・・(４５)
次に、無次元化応力η(n)及び除荷点からの無次元化応力ηL(n)を次式により更新する。
η(n)＝η(n-1)＋ｄη    ・・・(４６)
ηL(n)＝η(n)−ηRi    ・・・(４７)
次のステップ３１２では、（ｎ−１）ステップの降伏曲面が直前の除荷曲面よりも大きいか否かを判断する。そして、（ｎ−１）ステップの降伏曲面が直前の除荷曲面よりも大きい場合はステップ３１４へ移行し、そうでない場合は本ルーチンを終了する。具体的には、次式が成り立つか否かを判断する。なお、この判断の概念を図１０に示す。
【０１０８】
【数３２】
0004513776.tif 000033
【０１０９】
ステップ３１４では、（ｎ−１）ステップの降伏曲面を削除し、直前の除荷曲面の情報に置換する。具体的には、次のように直前の除荷曲面の情報を取り出す。まずｉ番目の除荷曲面データをクリアし、除荷点数ｉを１デクリメントし、ｉ＝０のときは骨格曲線に戻り、原点からの計算とし、以降は更新されたｉ番目の除荷点の情報を用いる。この処理の後は、ステップ３０６へ戻って再度接線剛性を計算し、無次元化応力を求め直す。
【０１１０】
図３へ戻って、ステップ２１０では、応力の体積成分を求める。まず、過圧密を仮定し、膨潤線の膨張係数κを用いて体積弾性係数Ｋを求める。κと体積圧縮係数ｍvの関係は次式で表される。
【０１１１】
【数３３】
0004513776.tif 000034
【０１１２】
ここで、σ’mは拘束圧（ただし、ここでは圧縮が正）、ｅ0は初期間隙比である。
【０１１３】
なお、ｄσ’m／σ’mが十分に小さい場合は、次式を用いることが出来る。
【０１１４】
【数３４】
0004513776.tif 000035
【０１１５】
また、体積圧縮指数ｍvと体積弾性係数Ｋとの関係は次式で表される。
【０１１６】
【数３５】
0004513776.tif 000036
【０１１７】
上記（５０）式、（５１）式より、体積弾性係数Ｋは次式で求めることができる。
【０１１８】
【数３６】
0004513776.tif 000037
【０１１９】
なお、Ｇは無次元化接線剛性ｇより次式により計算する。
【０１２０】
Ｇ＝ｇＧ0        ・・・(５３)
次に、応力の体積成分を次式で求める。
【０１２１】
ｄσij＝ｄｓij＋δＫｄεv    ・・・(５４)
ここで、ｄεvは体積ひずみ増分であり、ｄεv＝ｄεkkである。
【０１２２】
次のステップ２１２では、増分応力を求めると共に次式により応力を更新する。
【０１２３】
σn＝σn-1＋ｄσ    ・・・(５５)
次に、ステップ２１４において、累積損傷度Ｄを求める。
【０１２４】
累積損傷度Ｄは、本実施形態では一例として除荷点から次の除荷点までを半波の区切りとし、除荷点からの増分相当応力に対して逐次計算する。除荷点の偏差応力ＳRijからの偏差応力増分について相当応力σLe(n)を次式により求める。
【０１２５】
【数３７】
0004513776.tif 000038
【０１２６】
なお、除荷点からの相当応力σLe(n)の１／２を半波のせん断応力とし、応力比Ｒは次式で計算することができる。
【０１２７】
【数３８】
0004513776.tif 000039
【０１２８】
このＲを用いて、累積損傷度増分ΔＤを求める。すなわち、図１１に示すように液状化強度曲線から応力比Ｒに対応する繰り返し回数Ｎifを求め、この繰り返し回数Ｎifから、除荷点からの累積損傷度増分ΔＤを計算する。なお、累積損傷度Ｄは、前述した特許文献２と同様に求めることができる。すなわち、累積損傷度Ｄは、同図に示すように、半波開始点ｔS、すなわち、除荷点における累積損傷度をＤRとして次式で表すことができる。
【０１２９】
Ｄ＝ＤR＋ΔＤ      …（５８）
また、累積損傷度ＤRは、相当応力σLe(n)が除荷点となる毎に更新される。すなわち、図１１の場合には、半波開始点ｔS、半波終了点ｔeの時点において、そのときの累積損傷度Ｄが累積損傷度ＤRに設定される。累積損傷度増分ΔＤは、半波開始点ｔs、半波終了点ｔeの時点においてゼロに初期化される。
【０１３０】
ステップ２１６では、過剰間隙水圧比増分を求める。累積損傷度Ｄから過剰間隙水圧比ｒuを求める方法は、例えば従来のＳｅｅｄの次式による方法や、地盤調査より得られる実測値のテーブルデータを用いる方法を用いることができる。ここで、過剰間隙水圧比ｒuを求めた後、過剰間隙水圧比の増分Δｒuを計算する。
【０１３１】
【数３９】
0004513776.tif 000040
【０１３２】
次のステップ２１８では、平均有効応力の変動ｄσ’mを次式により求める。
ｄσ’m＝−ｄｒu・｜σ’m0｜  ・・・(６０)
一方、ステップ２０２で過去にサイクリックモビリティ状態に入ったことがあると判断された場合、２０６でサイクリックモビリティ状態に入っていると判断された場合には、ステップ２２０において、有効応力経路が変相線２０を越えた載荷状態か否かを判断する。すなわち、せん断応力τが前回のステップで計算したそれより増加しているか否かを判断する。
【０１３３】
そして、載荷状態であると判断した場合には、ステップ２２２へ移行し、載荷状態でない、すなわち除荷状態であると判断した場合は、ステップ２３０へ移行する。
【０１３４】
ステップ２２２では、図１２（ａ）、（ｂ）に示すように変相線２０を越えるときのせん断剛性（せん断弾性係数の接線勾配）ＧCMと平均有効応力σ’mとから接線せん断剛性ＧTを求める。
【０１３５】
せん断剛性ＧCMは、初期剛性Ｇm0と応力経路が変相線２０と交差した時の平均有効応力σ’mpを用いた、例えば次式により求めることができるが、これに限られるものではない。
【０１３６】
【数４０】
0004513776.tif 000041
【０１３７】
ここで、Ｒsmは有効応力低下率であり、次式で表される。
【０１３８】
【数４１】
0004513776.tif 000042
【０１３９】
なお、上記（６１）式は、有効応力低下率Ｒsmと接線剛性低下率（ＧCM／Ｇm0）との関係を示す試験結果から導き出される。
【０１４０】
そして、接線せん断剛性ＧTは、接線剛性ＧCMを初期値として平均有効応力σ’mpからの平均有効応力比（σ’m(n-1)−σ’mp）／σ’m0により次式で表される。
【０１４１】
【数４２】
0004513776.tif 000043
【０１４２】
ここで、ｆb（γmax）は接線勾配増加率を表す関数であり、最大せん断ひずみγmaxを用いて次式で表される。
【０１４３】
【数４３】
0004513776.tif 000044
【０１４４】
ここで、ａ，ｂ，ｃは所定のパラメータであり、ｂはｆb（γmax）の最小値、ａ，ｂは試験結果のデータを用いた回帰分析により求めることができる。
【０１４５】
なお、上記（６３）式は、平均有効応力σ’mpからの平均有効応力比と接線剛性低下率との関係を表す試験結果から導き出され、上記（６４）式は、最大せん断ひずみγmaxと接線勾配増加率ｆｂとの関係を示す試験結果から導き出される。
【０１４６】
次に、ステップ２２４において、体積弾性係数Ｋを算出する。これは、サイクリックモビリティ状態に入る前と同様に、すなわちステップ２１０と同様に（５２）式により求めることができる。
【０１４７】
次に、ステップ２２６において、増分応力の計算を行うと共に応力を更新する。具体的には、まず、３次元の応力ひずみ関係マトリックスＤを組み立て、増分応力計算に用いる。Ｄマトリックスを下式に示す。
【０１４８】
【数４４】
0004513776.tif 000045
【０１４９】
増分応力ｄσ、応力σnは次式により求めることができる。
ｄσ＝Ｄｄε    ・・・(６６)
σn＝σn-1＋ｄσ    ・・・(６７)
次に、ステップ２２８において、サイクリックモビリティによる平均有効応力の増分Δσ’m（ｄσ’m）を次式により求める。これにより、破壊線に沿って平均有効応力が増加する。
【０１５０】
【数４５】
0004513776.tif 000046
【０１５１】
ここで、σeは相当応力、τmaxは破壊時のせん断応力（τmax＝σ’m（ｓｉｎφ＋ｃｃｏｓφ）、Ｍ0は変相線の破壊線との比（＝ｔａｎφp／ｔａｎφf）、φfは破壊角、φpは変相角、ｄσeは増分相当応力を示す。
【０１５２】
一方、ステップ２２０で載荷状態ではなく除荷状態であると判断した場合は、ステップ２３０において、せん断ひずみγによる剛性低下Ｇ／Ｇ0と平均有効応力σ’mとから接線せん断剛性ＧTを求める。
【０１５３】
具体的には、除荷点からの無次元化ひずみξの１／２を用いて骨格曲線／履歴曲線Ｈ（ξ）で無次元化剛性ｇを求め、せん断弾性係数の接線勾配ＧTを計算する。ただし、ここで用いる初期剛性Ｇ0Uは、平均有効応力σ’mに比例して低下するものとする。
【０１５４】
ＧT＝ｇＧ0U      ・・・（６９）
ここで、初期剛性Ｇ0Uは次式で表される。
【０１５５】
【数４６】
0004513776.tif 000047
【０１５６】
なお、Ｇm0は初期平均有効応力σ’m0のときの初期剛性である。
【０１５７】
ＧTが図４（ａ）に示す最小剛性（変相線２０上の接線剛性ＧCM）より小さい時は、ＧCMを用いる。ＧCMは、上記（６１）式により求めることができる。
【０１５８】
また、図１３（ａ），（ｂ）に示すように、除荷点が変相線を越えていない場合も同様に接線勾配を求める。
【０１５９】
次に、ステップ２３２において、体積弾性係数Ｋを算出する。これは、サイクリックモビリティ状態に入る前と同様に、すなわちステップ２１０と同様に（５２）式により求めることができる。
【０１６０】
次に、ステップ２３４において、ステップ２２６と同様に増分応力の計算を行うと共に応力を更新する。
【０１６１】
ステップ２３６では、サイクリックモビリティによる平均有効応力増分Δσ’mを求める。平均有効応力増分Δσ’m（ｄσ’m）は、次式で表される。この式は、経験式であり、下記（７４）式を微分して増分形式としたものである。下記（７４）式は、平均有効応力とせん断応力との関係を示す中空ねじり試験の応力経路を変相線上で１になるようにノーマライズした除荷点からの応力比と平均有効応力低下率の関係を表す。
【０１６２】
【数４７】
0004513776.tif 000048
【０１６３】
ここで、σ’mRは除荷点の平均有効応力を、σ’mpは変相線上の平均有効応力を示す。また、ηLは除荷点からの応力比であり、次式で表される。
【０１６４】
【数４８】
0004513776.tif 000049
【０１６５】
ここで、ηRは除荷点上の応力比を、ηp（＝ｔａｎφp：φpは変相角）は変相線上の応力比を示す。ｄηは応力比増分であり、次式で表される。
【０１６６】
【数４９】
0004513776.tif 000050
【０１６７】
なお、ここでの応力比の定義はη＝τ／｜σ’m｜とする。
【０１６８】
除荷点からの平均有効応力低下率（σ’m−σ’mR）／（σ’mp−σ’mR）は次式で表される。
【０１６９】
【数５０】
0004513776.tif 000051
【０１７０】
なお、図１４（ａ），（ｂ）に示すように、変相線２０を越える前に応力が反転する場合も、その点を除荷点として上記（７１）式を適用し、平均有効応力の変動分を求める。
【０１７１】
変相線２０上の平均有効応力σ’mpはＮ／Ｎf−ｒu関係より求める。Ｎ／Ｎf−ｒu関係は図１５に示すように、サイクリックモビリティで通るｒuの最大値の包絡線を定義したものである。
【０１７２】
次のステップ２３８では、平均有効応力の変動分ｄσ’mを考慮した応力を次式により求める。
σij＝σ’ij＋δijｄσ’m    ・・・(７５)
ステップ２４０では、求めた応力σから内力Ｒａを次式により求める。
【０１７３】
【数５１】
0004513776.tif 000052
【０１７４】
ここで、Ｖは対象とする体積であり、Ｎは節点（計算位置）への変換関数である。
【０１７５】
図２に戻って、ステップ１０８では、外部から作用する外力Ｆ、粘性力ＦC、慣性力ＦMなどからなる外力項とステップ１０６で得られた内力Ｒａとの残差力（不釣り合い力）ΔＲを求める。
【０１７６】
次に、ステップ１１０において、上記剛性を用いて要素行列の算出とこれを用いた全体剛性の算出を行う。
【０１７７】
次に、ステップ１１２において、動的解析用全体行列Ｋ*を算出する。動的解析用全体行列Ｋ*は次式で示される。
Ｋ*＝Ｍ＋γΔｔＣ＋βΔｔ2Ｋ    …（７７）
ここで、β、γは、Ｎｅｗｍａｒｋ−β法における係数（一定）である。また、Ｋは接線剛性である。
【０１７８】
次に、ステップ１１４において、残差力より（ｔ＋Δｔ）時刻の修正値Λを下記（７８）式により算出し、この修正値Λにより予測値ｕ*を修正し、（ｔ＋Δｔ）時刻の応答変位、応答速度、応答加速度を求め、これにより外力Ｆ、慣性力ＦM及び粘性力ＦCを算出する。
【０１７９】
Δｔは計算時間間隔、すなわち時刻歴ループの実行間隔である。なお、Δｔは一定の値でもよいし、計算毎に変化させてもよい。
【０１８０】
【数５２】
0004513776.tif 000053
【０１８１】
次に、ステップ１１６において、ステップ１０６と同様に内力Ｒａを算出する。次に、ステップ１１８において、ステップ１０８と同様に外力と内力の差（残差力）ΔＲを算出する。
【０１８２】
次に、ステップ１２０において、この残差力ΔＲが許容値以下であるか否かを判断し、許容値を越えていれば、ステップ１１０から反復して計算する。この残差力ΔＲが許容値以下であれば、次のステップ１２２へ進む。
【０１８３】
次に、ステップ１２２において、変位、速度、加速度等を更新すると共に、算出した各応答値を表示部１８や記憶部１４へ出力する。出力値としては、変位、速度、加速度、ひずみ、応力等がある。
【０１８４】
このようにして、地震時間全体についてステップ１０４からステップ１２２までの処理（時刻歴ループ）を行い、各時刻毎の応答値を算出する。
【０１８５】
以上説明したように、本実施形態では、除荷点から次の除荷点までを半波の区切りとし、３次元応力空間における除荷点からの増分相当応力に基づいて累積損傷度を逐次計算するので、簡便かつ精度良く地震の３次元挙動を解析することができる。
【０１８６】
また、有効応力経路が変相線を越える前と後とで異なる手法により平均有効応力の変動を求めるので、サイクリックモビリティの影響が反映された地震応答解析を行うことができる。
【０１８７】
なお、本実施形態では、３次元応力空間における除荷点からの増分相当応力に基づいて累積損傷度を逐次計算する場合について説明したが、従来（例えば前記特許文献２記載の技術）のように１次元の応力、すなわちせん断応力に基づいて累積損傷度を逐次計算するようにしてもよい。この場合、増分ひずみの計算において、上記の１節点での変位成分はｕT＝｛ｕ｝として、上記（６）式は下記（７９）式に、上記（７）式は下記（８０）式に置き換える。
εT＝｛γ｝    ・・・（７９）
ｄγ＝γ(n)−γ(n-1)    ・・・（８０）
なお、γはせん断ひずみである。
【０１８８】
また、累積損傷度の計算において、上記（５６）式はτ（せん断応力）に、上記（５７）式は下記（８１）式に置き換える。
【０１８９】
【数５３】
0004513776.tif 000054
【０１９０】
また、増分応力の計算、応力の更新処理において、上記（６５）式は下記（８２）式に、上記（６６）式は下記（８３）式に、上記（６７）式は下記（８４）式に置き換える。
Ｄ＝［ＧT］      ・・・（８２）
ｄτ＝ＧTｄγ    ・・・（８３）
τn＝τ(n-1)＋ｄτ  ・・・（８４）
また、本実施形態では、除荷点から次の除荷点までを半波の区切りとした場合について説明したが、これに限らず、半波の開始点及び終了点は任意の位置に設定することができる。
【図面の簡単な説明】
【０１９１】
【図１】地震応答解析装置の概略ブロック図である。
【図２】演算部において実行される制御ルーチンのフローチャートである。
【図３】内力の計算処理のフローチャートである。
【図４】（ａ）は接線剛性について説明するためのせん断ひずみとせん断応力との関係を示す図、（ｂ）は有効応力経路について説明するための図である。
【図５】剛性低下の計算処理のフローチャートである。
【図６】降伏曲面と除荷曲面について説明するための図である。
【図７】骨格曲線及び履歴曲線について説明するための図である。
【図８】（ａ）は剛性低下率とせん断ひずみ、減衰との関係を示す図、（ｂ）はせん断ひずみとせん断応力との関係を示す図である。
【図９】各区間の接線剛性について説明するためのせん断ひずみとせん断応力との関係を示す図である。
【図１０】降伏曲面が直前の除荷曲面よりも大きいか否かの判断について説明するための図である。
【図１１】相当応力から液状化強度曲線を用いて累積損傷度を求める場合について説明するための図である。
【図１２】（ａ）は接線剛性について説明するためのせん断ひずみとせん断応力との関係を示す図、（ｂ）は有効応力経路について説明するための図である。
【図１３】（ａ）は除荷点が変相線を越えていない場合において接線剛性を求める場合について説明するためのせん断ひずみとせん断応力との関係を示す図、（ｂ）は有効応力経路について説明するための図である。
【図１４】（ａ）、（ｂ）は変相線を越える前に応力が反転する場合において平均有効応力の変動分を求める場合について説明するための有効応力経路を示す図である。
【図１５】Ｎ／Ｎf−ｒu関係を示す図である。
【符号の説明】
【０１９２】
１０      地震応答解析装置
１２      操作部
１４      記憶部
１６      演算部
２０      変相線
２２      有効応力経路
２４      破壊線
【図１】
0004513776.tif 000055
【図２】
0004513776.tif 000056
【図３】
0004513776.tif 000057
【図４】
0004513776.tif 000058
【図５】
0004513776.tif 000059
【図６】
0004513776.tif 000060
【図７】
0004513776.tif 000061
【図８】
0004513776.tif 000062
【図９】
0004513776.tif 000063
【図１０】
0004513776.tif 000064
【図１１】
0004513776.tif 000065
【図１２】
0004513776.tif 000066
【図１３】
0004513776.tif 000067
【図１４】
0004513776.tif 000068
【図１５】
0004513776.tif 000069