Vol.26, No.1, 1-17, (2019)
時間領域差分法による表層地盤の非線形応答の計算
Computation of non-linear site response
by the time-domain finite-difference method
鳥
越 友 輔 (Yusuke T
ORIGOE)
*小
松 正 直 (Masanao K
OMATSU)
**竹
中 博 士 (Hiroshi T
AKENAKA)
**Abstract
We have developed a time-domain staggered-grid finite-difference code for modeling non-linear response of a one-dimensionally inhomogeneous subsurface structure to a SH plane-wave incidence. It employs the velocity-stress formulation of elastodynamic equation for the linear part, and adopts a elastoplastic rheology model for the non-linear relation between the stress and strain. In this paper, we apply this code to four constitutive models from linear-elastic to nonlinear: (1) linear elastic model, (2) linear viscoelastic model, (3) elastoplastic model, and (4) viscoelastoplastic model, which simulate shallow sand and clay structures and are vibrated by a vertically incident SH plane-wave of Ricker wavelet, to compare the linear and the non-linear soil behaviors including low strains damping (viscoelastic effect) and/or hysteretic attenuation (non-linear effect). We also apply it to a local strong-motion record of the 2000 Western-Tottori earthquake (MW6.8). We then simulate
characteristics of non-linear site response such as reduction of the spectral amplitude in the high frequency band and shift of the peak frequencies to lower frequencies. Keywords: finite-difference method, non-linear site response, strong motion
1. はじめに 強震動の評価をするためには，地震の震源特性， 伝播経路特性，地盤特性の 3 つの特性を把握する 必要がある．その内，地盤特性とは地震基盤上に存 在する地盤(堆積層)が地震動に与える影響のこと である．特に，やわらかく厚い地盤では，地震波の 増幅特性により大きな被害を引き起こすこともあ る．このような地盤特性は強震動に大きな影響を 与えるため，地盤特性を評価することは重要であ る． 軟弱な表層地盤では，弱震動によって生じる小 さなひずみでは応力－ひずみ関係は線形であるが， 強震動によって生じる大きなひずみでは応力－ひ * 岡山大学理学部地球科学科 (現所属：岡山大学大学院環境生命科学研究科), 〒700-8530 岡山市北区津島中三丁目１－１
Department of Earth Sciences, Faculty of Science, Okayama University, Okayama 700-8530, Japan (Present affiliation: Graduate School of Environmental and Life Science, Okayama University)
** 岡山大学大学院自然科学研究科，〒700-8530 岡山市北区津島中三丁目１－１
Graduate School of Natural Science and Technology, Okayama University, Okayama 700-8530, Japan.
ずみ関係が線形から外れ，剛性率が低下し減衰が 増大する(例えば，野口, 2009)．これを表層地盤の 非線形応答と呼ぶ．地震の非線形応答の解析法に は大きく分けて等価線形化法，非線形逐次積分法 がある(例えば，吉田, 2010)．このうちの等価線形 化法は運動方程式を周波数領域で解析する方法で， 繰返しせん断特性から得られる剛性率 𝐺 および 減衰定数 ℎ を用いて，非線形応答と等価な線形応 答を計算するものである(例えば，Schnabel et al., 1972)．それに対して，非線形逐次積分法は運動方 程式を逐次積分に基づいて時間領域で解析する方 法で，時々刻々と変化する材料の接線剛性を用い て計算するものである．等価線形化法は簡便であ
るため実務でよく用いられている．大ひずみ領域 では，最大加速度を過大評価すること(吉田, 1994) や ， 高 周 波 数 成 分 の 増 幅 比 を 過 小 評 価 す る (Masuda et al., 2001)という問題があったが，改良 が進められている(例えば，杉戸・他, 1994; Kausel and Assimaki, 2002; Yoshida et al., 2002)．一方， 非線形逐次積分法は，原理的に小ひずみ領域から 大ひずみ領域まで広い範囲のひずみレベルでも計 算可能である．本研究では非線形逐次積分法を用 いて計算を行う． 表層地盤の非線形挙動による影響としては，ピ ーク周波数の低周波数側へのシフト，および高周 波数帯の増幅率の低下が挙げられる(例えば，時松・ 翠川，1988；野口・笹谷, 2011)．地震動の振幅が増 加するとピーク周波数が低周波数側へシフトする のは，剛性率 𝐺 が低下することによって，S 波速度 が遅くなることに起因すると考えられる．また，高 周波数帯の増幅率の低下は減衰定数 ℎ の増大に 関連付けられる(Fig. 1)．地盤の非線形挙動は一般 に表層地盤が軟弱なほど生じやすく，解放基盤で は生じないと考えられている． 本研究の目的は，表層地盤の非線形応答による 地盤挙動の変化を確認することである．そのため に，まず非線形応答の計算を行うための非線形逐 次積分法の計算プログラムを作成した．次に，単純 な入射波形を用いて，入射波の強さや中心周波数 の変化，また地盤材料の違いが非線形性の発達に 及ぼす影響を数値的に調べた．さらに，実際に観測 された地震波形を入射波として用いて，複雑な入 射波が表層地盤の非線形現象に及ぼす影響を調べ た． 2. 計算手法 今回，線形応答の計算には時間領域差分法(田中・ 竹中, 2005)を用いた．非線形応答の計算について はせん断応力の評価に Bardet and Tobita (2001) により提案された非線形 1 次元地盤応答解析法 (NERA, Nonlinear Earthquake site Response Analysis)のアルゴリズムを用いた．これらの計算 手法について以下で説明する．座標系は，𝑧 軸を鉛 直下向きを正とする右手系の 3 次元デカルト座標 系を採用する． 2.1. 線形の支配方程式と差分解法 田中・竹中 (2005)より，深さ方向にのみ物性値 が変化する媒質の平面波入射問題における SH 波 の速度－応力型の支配方程式は 𝜕𝑣 𝜕𝑡 1 ⊿ 𝜕𝜏 𝜕𝑧 , 1 𝜕𝜏 𝜕𝑡 𝐺 𝜕𝑣 𝜕𝑧 . 2 ただし， ⊿ 𝜌 𝐺𝑝 . 3 ここで，𝑡は時間，𝑣 は y 方向の地動速度，𝜏 はせ ん断応力，𝐺 は剛性率，𝜌 は密度，𝑝 は水平スロー ネスである．ただし，本研究では対象が表層の軟弱 地盤であることから鉛直入射を仮定し，𝑝 0 とす る．この方程式を数値的に解くために，スタガード Fig. 1. 土質材料の非線形せん断特性の概要． Fig. 2. スタガード格子の配置. 空間配置(左)， 時間配置(右). 田中・竹中 (2005)を一部改変．
格子差分法(Virieux, 1984; Levander, 1988)を用い て離散化する．スタガード格子の配列を Fig. 2 に 示す．スタガード格子は，地動速度とせん断応力の 成分を半格子分ずらして配置する離散化の手法で ある．従来の差分法と比べて数値的安定性が優れ ているので，スタガード格子差分法は地震動の計 算で多用されている(竹中, 1998)． 式(1)，(2)を時間 𝑡 に関して 2 次精度，空間 𝑧 に 関して 4 次精度の差分で離散化するとそれぞれ以 下のようになる(田中・竹中, 2005)． 𝑣 𝑣 ⊿𝑡 1 ⊿ 𝑐 𝜏 ⁄⁄ 𝜏 ⁄⁄ ⊿𝑧 𝑐 𝜏 ⁄ ⁄ 𝜏 ⁄⁄ ⊿𝑧 , 4 𝜏 ⁄⁄ 𝜏 ⁄ ⁄ ⊿𝑡 𝐺 𝑐 𝑣 𝑣 ⊿𝑧 𝑐 𝑣 𝑣 ⊿𝑧 . 5 ここで，式中の添字について 𝑣 を例にとって説明 すると，深さ𝑧 𝑖 1 ⊿𝑧 ，時刻𝑡 𝑛 1 ⊿𝑡 に おける 𝑣 を表している．また，定数𝑐 ，𝑐 は等間隔 格子の場合，4 次精度では𝑐 9 8⁄ ，𝑐 1 24⁄ である．一方 𝑐 1 ，𝑐 0 とすると 2 次精度と なる． 2.2. 非線形地盤応答の計算アルゴリズム 非線形応答の計算を行うために，式(5)のせん断 応力の計算を Bardet and Tobita (2001)によって 提案されたアルゴリズムに置き換える．このアル O A B C D E F G H O R1 R2 R3 2R1 2R2 2R3 2R1 Strain St re ss R1 R2 R3 A B C D H1 H2 H3 Strain St re ss
𝑘
𝑘
𝑘
𝑘
𝑅
𝑅
𝑅
𝛾
𝜏
Fig. 3. Iwan(1967)と Mróz(1967)が提案したレオロジーモデル．ここでは，地震学の慣習 に従い，応力は引っ張りを正としている． Fig. 4. (a)載荷中の骨格曲線と(b)繰り返しの載荷・除荷中の履歴的な応力-ひずみ関係． O (a) (b)ゴリズムでは土の非線形的な応力－ひずみ関係を モデル化するために，Iwan (1967) と Mróz (1967) によって提案されたレオロジーモデルを導入して いる．以降，このモデルを IM モデルと呼ぶ．IM モデルはFig. 3 に示すように，バネとスライダー から成る複数の要素を直列に連結することで構成 さ れ て い る ．𝑖 番 目 の バ ネ は バ ネ 定 数 𝑘 𝑖 1, 2, ⋯ , 𝑛 を有し，𝑖 番目のスライダーは降伏応力 𝑅 𝑅 𝑅 ⋯ 𝑅 )を有する．全てのスライダ ーの始めの残留応力は0 である．単調載荷中は，𝑖 番目のスライダー 𝑖 はせん断応力 𝜏 が 𝑅 に達し たときに降伏し，降伏後は 𝑅 に等しい正の残留応 力を持つ．3 個のバネ，3 個のスライダーで構成さ れた IM モデル( 𝑛 3 )によって描かれる応力－ひ ずみ関係を Fig. 4 に示す．応力―ひずみ関係は折 れ線で表現され，その傾き𝐻 は以下のようにバネ 定数 𝑘 の調和平均に等しい．
Fig. 5. ひずみ増分から応力を計算する Bardet and Tobita (2001)のアルゴリズムの フローチャート．
𝐻 ⎩ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎧ 𝐻 𝑘 if 0 𝜏 𝑅 , 𝐻 𝑘 𝑘 if 𝑅 𝜏 𝑅 , ⋮ 𝐻 𝑘 𝑘 ⋯ 𝑘 if 𝑅 𝜏 𝑅 , 𝐻 𝑘 𝑘 ⋯ 𝑘 𝑘 if 𝑅 𝜏 𝑅 , 0 if 𝜏 𝑅 . 6
Fig. 4(a)の折れ線 OABCD は，単調載荷を受ける 際の応力－ひずみ関係(骨格曲線)である．載荷過程 で変化するバックストレス𝛼 を導入すると，載荷 中は𝜏 が 𝛼 𝑅 ，除荷中は 𝜏 が 𝛼 𝑅 に達した ときにスライダー𝑖 は降伏する．繰返し載荷を生じ るときの応力－ひずみ関係は Masing 則(Masing, 1926)に従い，Fig. 4(b)に示すような履歴曲線を描 く．Fig. 4(b)の折れ線 DEFGH は骨格曲線(折れ線 OABCD)を相似形に 2 倍大きくし，始点を除荷点 に移動したものである．
Bardet and Tobita (2001)のアルゴリズムのフロ ーチャートを Fig. 5 に示す．彼らのアルゴリズム では，ひずみ増分 ⊿𝛾 から 𝜏 を求める．⊿𝛾 は 以下の差分式によって得ることができる． 𝑣 𝑣 ∆𝑡𝑣 , 7 𝛾 𝑣 𝑣 ⊿𝑧 , 8 ⊿𝛾 𝛾 𝛾 . 9 ここで，𝑣 は y 方向の地動変位である．フローチャ ート中では，載荷か除荷かを変数 𝑥 を用いて表し， 載荷のとき 𝑥 1 ，除荷のとき 𝑥 1 としている． Fig. 4 の骨格曲線は，剛性率𝐺 とせん断ひずみ 𝛾 の 関 係 を 𝑛 個 の デ ー タ 点 ， 𝐺 , 𝛾 𝑖 1, . . . , 𝑛 を用いて表したものである．この場合，傾 き𝐻 は以下のように求まる． 𝐻 𝐺 , 𝐻 𝐺 𝛾 𝐺 𝛾 𝛾 𝛾 , 𝑖 2, … , 𝑛 1 , 𝐻 0 . (10) ここで，𝐺 は線形のときの剛性率である．バックス トレス𝛼 が最初は 0 であると仮定すると，𝑅 は 𝑅 𝐺 𝛾 , 𝑖 1, … , 𝑛 . 11 式(10)，式(11)は，𝐺 𝐺⁄ ― 𝛾 関係を用いると 𝐻 𝐺 𝐺 𝛾 𝐺 𝛾 𝛾 𝛾 , 𝑖 2, … , 𝑛 1, 12 𝑅 𝐺 𝐺 𝛾 , 𝑖 1, … , 𝑛. 13 ここで，𝐺′ 𝐺 𝐺 ⁄ である．以上のアルゴリズムを 実装したFortran 77 のサブルーチンの例を付録に 示す． 2.3. 非弾性減衰の導入 実際の地盤を考えるとき，応力－ひずみ関係の 非線形性に起因する履歴減衰だけでなく，非弾性 減衰(粘性減衰)の効果を考慮する必要がある． Graves (1996)により提案された方法を用いると 非弾性減衰の効果を簡便に導入することができる (例えば, 林田・他, 1999)．林田・他 (1999)では， 各時間ステップの速度と応力を計算する際に，以 下の係数 𝐴 を 1 つ前の時間ステップの値にかける ことによって非弾性減衰の効果を導入する． 𝐴 exp 𝜋⊿𝑡𝑄 𝑇 . 14 ここで，𝑄 は Q 値の逆数，𝑇 は減衰性を考慮す る中心周期，⊿𝑡 は時間ステップである．本研究で は，𝑄 については S 波に対する Q 値(𝑄 )を用い る．(14)式の Q の値は周波数1/𝑇 の Q 値と考える ことができ，(14)式の適用は周波数の 1 乗に比例す る周波数依存性を持った Q と等価な効果を与える． 非線形性を考慮しない場合の計算では，林田・他 (1999)の方法に従う．非線形性を考慮する場合の計 算では，2.2 節に示したアルゴリズムによる計算に おいて，応力 𝜏 とバックストレス α の両方に 𝐴 をかけるスキームと， 𝜏 のみに𝐴 かけるスキーム の両方を試したが，両者の結果にほとんど違いが 見られなかったため，以下，𝜏 のみに 𝐴 をかける
スキームを採用する． 2.4. 境界条件 今回，境界条件として，田中・竹中（2005）と同 じ手法を採用し，地表で自由表面条件，最下層で無 反射境界条件を適用した．無反射境界条件とは，計 算領域の端において反射しないで波を透過させる 条件のことである．無反射境界条件には，計算領域 の下端からの人工的な反射波を最小限に抑えるた め，one-way 波動方程式に基づく吸収境界条件 (Clayton and Engquist, 1977) と Cerjan et al. (1985)の吸収境界条件を併用した．
3. プログラムの検証
本研究で作成した計算プログラムと，Bardet and Tobita (2001) により提案された非線形 1 次元 地盤応答解析法(NERA, Nonlinear Earthquake site Response Analysis)において，同一条件での解 析結果を比較し，作成した計算プログラムの検証 を行う．NERA のプログラムは飛田氏のサイト [https://sites.google.com/site/tt60898/home/softw are (最終閲覧日：2017 年 12 月 13 日)]で公開され ており，表計算プログラムの Excel を用いて作ら れている．比較にはNERA for Windows 10-64bit 版を用い，計算にはExcel 2016 を使用した．本研 究で作成したプログラムとNERA では差分法のス キームが異なる．時間に関してはどちらも 2 次精 度のスキームで離散化しているが，空間に関して は，今回作成した計算プログラムでは 4 次精度で あるのに対し，NERA では 2 次精度のスキームで 離散化している． 本節の計算及び次節の計算例のケース1～3 では， リッカー・ウェーブレット(Ricker wavelet)を用い て非線形性の発達に及ぼす影響を調べる．リッカ ー・ウェーブレットは載荷・除荷のサイクルが少な いため現実的ではないが，波の大きさや中心周波 数を変化させることが容易である．また計算にか かる時間も短いため，非線形挙動の素過程を調べ るためには適している．本研究では波の大きさの 指標として解放基盤における PGA を使用する． PGA は一般に地震動の指標として非線形の研究で 使用されている(例えば, Beresnev and Wen, 1996)． リッカー・ウェーブレットの時系列の加速度波形 と加速度フーリエスペクトルを Fig. 6 に示す．こ の波形は時間領域では， 𝑤 𝑡 PGA 2 2𝜋 𝑓 𝑡 1 exp 𝜋 𝑓 𝑡 , 15 周波数領域では， 𝑊 𝑓 PGA √𝜋 𝑓 𝑓 exp 𝑓 𝑓 , 16 と表せる．ここで， 𝑓は周波数，𝑓 は中心周波数で ある． プログラムの検証に用いる地盤構造モデルを Fig. 7 に示す．厚さ 30 m の砂層からなる表層地盤 と半無限の工学的基盤の2 層モデルで，物性値は， S 波速度𝑉 は表層地盤 300 m/s，工学的基盤 700 m/s ， 密 度 𝜌 は 表 層 地 盤 18.0 kN/m3 (=1834.9 kg/m3)，工学的基盤 21.0 kN/m3 (=2140.7 kg/m3)と した．このモデルの工学的基盤からリッカー・ウェ ーブレットの加速度波形を積分した速度波形を鉛 直入射して差分法計算(空間 4 次精度，時間 2 次精 度)を実施し，その結果を数値微分(差分)して出力 Fig. 6. (a)入射加速度波形とその(b)フーリエス ペクトル． (a) (b)
波形(地動加速度)を得た．解析に使用した差分法計 算のパラメーターの値をTable 1 に示す．⊿𝑡 は時 間間隔，𝑁 は時間ステップ数，⊿𝑧 は深さ方向の格 子間隔，𝑁 は深さ方向の格子数である． 解析には土質材料の剛性比𝐺 𝐺⁄ とせん断ひずみ 𝛾の関係を表すデータが必要である．ここでは， Bardet and Tobita (2001) に 砂 の サ ン プ ル (“material type No. 2”)として掲載されている値 (Table 2)を用いた．このデータは Seed and Idriss (1970)に基づいている． 計算した地表の加速度波形をFig. 8，深さ 5.1 m おける応力－ひずみ関係を Fig. 9 に示す．NERA と本計算プログラムの両者の結果がよく一致して いることから，本研究で作成した計算プログラム による結果は正しいと考えられる． ひずみ G/G0 0.000001 1.000 0.000003 1.000 0.00001 0.990 0.00003 0.960 0.0001 0.850 0.0003 0.640 0.001 0.370 0.003 0.180 0.01 0.080 0.03 0.050 0.1 0.035
数値はBardet and Tobita (2001)による Table 2. 使用した砂のせん断特性の値 30 m 𝑉 300 m/s 𝜌 18.0 kN/m3 𝑉 700 m/s 𝜌 21.0 kN/m3 表層地盤 工学的 基盤 Fig. 7. 検証に用いた地盤構造モデル． -5 -4 -3 -2 -1 0 1 2 3 0 0.5 1 1.5 2 2.5 A ccelera tion (m/ s/ s) Time (s) NERA 作成したプログラム -20 -15 -10 -5 0 5 10 15 -0.02 -0.01 0 0.01 0.02 Stre ss (k P a) Strain (%) NERA 作成したプログラム Fig. 8. NERA (青)と作成した計算プログラム(赤) による加速度波形の解析結果. Fig. 9. NERA (青)と作成した計算プログラム(赤) で求まった深さ 5.1 m における応力－ひずみ関 係. (s) 12500 0.2 (m) 5000 Table 1. 検証に用いた差分パラメーター
∆𝑡
∆𝑧
𝑁
𝑁
2.0 104. 解析ケースと結果 作成したプログラムを用いて，非線形応答の計 算を行う．本研究では4 つのケースの計算を行う． ケース 1 として，まず入射波の強さが非線形性の 発達に及ぼす影響を調べる．基盤への入射加速度 波形にリッカー・ウェーブレットを採用し，解放基 盤における地動最大加速度(PGA)を変化させて表 層地盤の応答を計算する．次に，ケース2 として， 入射波の周波数の違いによる影響を調べるために， リッカー・ウェーブレットの中心周波数𝑓 を変化 させた結果の比較を行う．さらに，ケース 3 とし て，砂で構成された表層地盤と粘土で構成された 表層地盤の非線形挙動の比較を行う．最後にケー ス 4 として，実波形を入射加速度波形として用い た計算を行う． 4.1. ケース 1 ケース1 では，リッカー・ウェーブレットを用い て入射波の強さの違いが非線形性の発達に及ぼす 影響を調べる．ケース 1 で用いる地盤構造モデル をFig. 10 に示す．厚さ 10 m の 3 つの砂層からな る表層地盤と半無限の工学的基盤の4 層モデルで， 物性値は，S 波速度𝑉 が表層 1 層目で 150 m/s，2 層目で200 m/s，3 層目で 300 m/s，工学的基盤で 700 m/s，密度𝜌 が表層地盤で 1800 kg/m3，工学的 基盤で2300 kg/m3である．粘性を考慮する際の𝑄 は表層地盤で 15，基盤で∞，減衰の中心周期𝑇 は 1/𝑓 とした．砂の剛性比𝐺 𝐺⁄ とせん断ひずみ𝛾の関 係は前節と同じデータを用いた．解析に使用した 差分法計算のパラメーター値をTable. 3 に示す． 本節では(ケース 4 を除いて)表層地盤に以下の 4 つの材料構成モデルを用いる． ① 弾性モデル：応力－ひずみ関係がフックの法則 に従って線形関係である．非弾性減衰の効果は 考慮しない． ② 粘弾性モデル：2.3 節で述べた非弾性減衰の効果 を考慮する． ③ 弾塑性モデル：表層地盤の非線形性を考慮する． ④ 粘弾塑性モデル：表層地盤の非弾性減衰の効果 と非線形性の両方を考慮する． リッカー・ウェーブレットの中心周波数𝑓 を 10 Hz とし，解放基盤における PGA を(a) 0.05g，(b) 0.10g，(c) 0.25g，(d) 0.50g とした場合の，砂の表 層地盤における，地表での計算加速度波形を Fig. 11，そのフーリエ振幅スペクトルを Fig. 12 に示す． ここで，g は重力加速度( 9.80 m/s2)である．PGA が0.05g のとき(Fig. 11(a))，ピークでの振幅は弾性 モデルと弾塑性モデルに大きな違いはないが，3 秒 付近からピークの生じる時間にズレが生じ始め， 弾塑性モデルで波の位相が遅れている．また，粘弾 性モデル及び粘弾塑性モデルでピークの振幅が， 弾性モデルに比べて小さくなっていることが確認 できる．Fig. 12(a)のフーリエスペクトルでは粘弾 性モデルで全周波数帯域にわたって振幅が弾性モ デルよりも小さくなっているが，これは非弾性減 衰の効果によるものである．また，弾塑性モデルで は，5 Hz 付近から高周波数側で振幅が減衰してい る．さらに，5 Hz 付近から高周波数の帯域で，ピ ークの低周波数側へのシフトが見られる．同様の ことが粘弾塑性モデルにも言えるが，粘弾塑性モ デルの場合は，非弾性減衰の効果によって最初に エネルギーの一部を散逸するため，非線形性の影 響は弾塑性モデルよりも小さくなり，高周波数部 分における振幅のピークの低周波数側へシフトは 弾塑性モデルよりも小さくなる．このように小さ なPGA の場合であっても励起された特定の周波数 帯域で非線形性が発達し得ることがわかる． 𝑉 150 m/s 𝜌 1800 kg/m3 10 m 10 m 10 m 𝑉 200 m/s 𝜌 1800 kg/m3 𝑉 300 m/s 𝜌 1800 kg/m3 𝑉 700 m/s 𝜌 2300 kg/m3
表層地盤
工学的
基盤
Fig. 10. リッカー・ウェーブレットによる解析 に用いる地盤構造モデル． (s) 12500 0.2 (m) 5000 Table 3. ケース1で用いた差分パラメーター∆𝑡
∆𝑧
𝑁
𝑁
2.0 10次に，入射波のPGA を大きくした場合の結果に ついて見る．Fig. 11(b)～(d)の時刻歴波形では，入 射波のPGA が増加するにつれて，非線形のモデル のピーク加速度が線形のモデルよりも小さくなっ ている．同様に，非線形のモデルでは入射波のPGA が増加するにつれて，振幅のピークが生じる時間 がより遅くなっている．この位相の遅れは，時間の 経過に伴って顕著である．また，Fig. 12(b)～(d)の Fig. 11. リッカー・ウェーブレットの𝑓 を 10 Hz とし，PGA を(a) 0.05g，(b) 0.10g，(c) 0.25g， (d) 0.50g とした場合の地表での加速度波形. (a) PGA=0.05g (b) PGA=0.10g (c) PGA=0.25g (d) PGA=0.50g 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0 5 10 15 20 25 Fo u ri er s pe ctr u m (m /s ) Frequency (Hz) 弾性 粘弾性 弾塑性 粘弾塑性 0.00 0.20 0.40 0.60 0.80 1.00 1.20 0 5 10 15 20 25 Fo u ri er s pe ctr u m (m /s ) Frequency (Hz) 弾性 粘弾性 弾塑性 粘弾塑性 Fig. 12. リッカー・ウェーブレットの𝑓 を 10 Hz，PGA を(a) 0.05g，(b) 0.10g，(c) 0.25g， (d) 0.50g とした場合の地表での加速度フーリ エ振幅スペクトル. (a) PGA=0.05g (b) PGA=0.10g (c) PGA=0.25g (d) PGA=0.50g
フーリエ振幅スペクトルでは，線形のモデルの場 合，PGA が 0.05g の Fig. 12(a)のスペクトルと比 較してそれぞれ，2 倍，5 倍，10 倍と大きくなって いる．これは，伝播する波形が入力波形のPGA の 大きさと線形に比例するという事実を示している． 非線形のモデルでは，入力波形のPGA が増加する につれてピーク加速度も振幅スペクトルも減衰し ている．Fig. 12(c)，(d)では，振幅スペクトルが 2 Hz 付近から高周波数側で減衰し，ピークが低周波 数側へシフトしていることが確認できる．また，周 波数が高くなるほど減衰が顕著となっている． 次に，Fig. 13 に PGA=0.25g のときの，深さに 対する最大せん断ひずみと最大せん断応力のプロ ファイル，さらに各構成モデルにおける深さ5.1 m での応力－ひずみ関係を示す．どのモデルにおい ても最大せん断ひずみと最大せん断応力は深さ 4 m 付近で最も大きくなっている．そして，弾塑性 及び粘弾塑性の非線形モデルでは，弾性及び粘弾 性の線形モデルよりも深い位置にピークが生じて おり，プロファイル曲線が尖っている．また，地盤 が粘弾性モデルである場合(Fig. 13②)の深さに対 する最大せん断ひずみと最大せん断応力のプロフ ァイルは弾性モデルの場合と比べて曲線の形は似 ているが，値は小さい．これは，振幅が非弾性減衰 の効果によって減衰しているためである．応力－ ひずみ関係は履歴曲線を描いており，この履歴曲 線に囲まれた面積はエネルギーが非弾性減衰を介 して散逸されたことを示している．地盤が弾塑性 モデルである場合(Fig. 13③)は，最大せん断ひずみ は弾性モデルの場合よりも小さいが，粘弾性モデ ルの場合よりは大きくなっている．しかし，弾塑性 モデルの最大せん断応力は弾性モデルおよび粘弾 性モデルよりも小さい．これは非線形性による減 衰の効果である．応力－ひずみ関係は粘弾性モデ ルと同様に履歴曲線を描いているが，履歴曲線に 囲まれた面積は粘弾性モデルよりも大きくなって おり，より減衰が大きいことがわかる．地盤が粘弾 塑性モデルの場合(Fig. 13④)では，最大せん断ひず みと最大せん断応力の両方で最も小さくなってい る．これは非弾性減衰と非線形性による履歴減衰 の両方が生じているためである． 4.2. ケース 2 ケース 2 では，解析条件についてはケース 1 と 同じものを用い，入射波リッカー・ウェーブレット の中心周波数 𝑓 が非線形性の発達に及ぼす影響 を調べる． リッカー・ウェーブレットの解放基盤における PGA が 0.10g のとき，中心周波数𝑓 を(a) 5 Hz， (b) 10 Hz，(c) 20 Hz と変化させた場合の砂の表層 地盤における，地表加速度のフーリエ振幅スペク トルをFig. 14 に示す．地盤が粘弾性モデル及び粘 弾塑性モデルの場合，𝑓 が大きくなるほど振幅が 小さくなっている．これは粘弾性による非弾性減 衰の効果が周波数に依存しているためである．周 波数が大きくなるほど非弾性減衰の効果は大きく なる．また，地盤が弾塑性の場合，𝑓 が大きくな -30 -25 -20 -15 -10 -5 0 0 0.02 0.04 0.06 De pt h ( m )
Maximum Shear strain (%) 弾性 粘弾性 弾塑性 粘弾塑性 -30 -25 -20 -15 -10 -5 0 0 10 20 30 De pt h ( m )
Maximum Shear stress (kPa) 弾性 粘弾性 弾塑性 粘弾塑性 -18 -12 -6 0 6 12 18 -0.06 -0.03 0 0.03 0.06 St re ss (k P a) Strain (%) -18 -12 -6 0 6 12 18 -0.06 -0.03 0 0.03 0.06 St re ss (k P a) Strain (%) -18 -12 -6 0 6 12 18 -0.06 -0.03 0 0.03 0.06 St re ss (k P a) Strain (%) -18 -12 -6 0 6 12 18 -0.06 -0.03 0 0.03 0.06 St re ss (k P a) Strain (%) Fig. 13. 深さに対する最大せん断ひずみと最大 せん断応力と，地盤が①弾性，②粘弾性，③弾塑 性，④粘弾塑性モデルの場合の深さ 5.1 m におけ る応力－ひずみ関係(𝑓 =10 Hz, PGA=0.25g の場 合). ①弾性モデル ②粘弾性モデル ③弾塑性モデル ④粘弾塑性モデル 最大せん断ひずみ 最大せん断応力
るほど減衰及びピークの低周波数側へのシフトが 生じる周波数が高くなっている．(a)では 2 Hz 付近 から高周波数側で減衰および低周波数側へシフト が現れているが，(b)では 5 Hz 付近，(c)では 8 Hz 付近から見られる．これは，リッカー・ウェーブレ ットの 𝑓 を大きくすると，振動継続時間が短くな ることに起因していると考えられる． 4.3. ケース 3 ケース 3 では，地盤の土質材料の違いによる影 響を調べるために，表層地盤が砂で構成されてい る場合と，粘土で構成されている場合の解析結果 の比較を行う．解析条件はケース1 と同じである． 土質材料の剛性比𝐺 𝐺⁄ とせん断ひずみ𝛾 の関係に ついては，砂はTable 2，粘土は Table 4 の値を用 いた．Table 4 の値は Bardet and Tobita (2001)に 粘土のサンプル(“material type No. 1”)として掲載 されており，Seed and Sun (1989)に基づいている． また，物性値のS 波速度𝑉 と密度𝜌 は砂層と粘土層 で同じ値とした． リッカー・ウェーブレットの中心周波数𝑓 を 10 Hz，解放基盤における PGA を(a) 0.05g，(b) 0.10g， (c) 0.25g，(d) 0.50g とした場合について，粘土表層 地盤上の地表加速度波形をFig. 15，そのフーリエ 振幅スペクトルをFig. 16 に示す．物性値が砂地盤 と同じであるので，線形モデルでは加速度波形も 振幅スペクトルも表層地盤が砂の場合(Fig. 11, 12) と同じ結果となっている．以下，PGA が 0.05g の ときから順に表層地盤が砂の場合の結果と比較す る．Fig. 15(a)では，弾性モデルと弾塑性モデルで 大きな違いはなく，表層地盤が砂の場合(Fig. 11(a)) と異なり，ピークの生じる時間のずれも生じてい ひずみ G/G0 0.000001 1.000 0.000003 1.000 0.00001 1.000 0.00003 0.981 0.0001 0.941 0.0003 0.847 0.001 0.656 0.003 0.438 0.01 0.238 0.03 0.144 0.1 0.110
数値はBardet and Tobita (2001)による Table 4. 使用した粘土のせん断特性の値 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0 5 10 15 F our ie r s pe ct rum ( m /s ) Frequency (Hz) 弾性 粘弾性 弾塑性 粘弾塑性 0.00 0.05 0.10 0.15 0.20 0.25 0 5 10 15 20 25 F our ie r s pe ct rum (m /s ) Frequency (Hz) 弾性 粘弾性 弾塑性 粘弾塑性 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0 5 10 15 20 25 30 35 40 45 F our ie r s pe ct rum (m /s ) Frequency (Hz) 弾性 粘弾性 弾塑性 粘弾塑性 (a) 𝑓 = 5 (Hz) (b) 𝑓 = 10 (Hz) (c) 𝑓 = 20 (Hz) Fig. 14. リッカー・ウェーブレットの PGA を 0.10g，𝑓 を(a) 5 Hz，(b) 10 Hz，(c) 20 Hz と した場合の地表での加速度フーリエ振幅スペク トル.
ない．Fig. 16(a)についても，弾塑性モデルの場合， 減衰が生じるのは表層地盤が砂の場合も粘土の場 合も周波数が5 Hz 付近からであるが，粘土は砂の 場合に比べて減衰が小さく，ピークの低周波数側 へのシフトも見られない．Fig. 15(b)～(d)では，入 射波のPGA が増加するにつれて砂の場合と同じよ -1.5 -1.0 -0.5 0.0 0.5 1.0 0 0.5 1 1.5 2 2.5 A cce le ra ti on ( m /s /s ) Time (s) 弾性 粘弾性 弾塑性 粘弾塑性 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 0 0.5 1 1.5 2 2.5 A cce le ra ti on ( m /s /s ) Time (s) 弾性 粘弾性 弾塑性 粘弾塑性 -6.0 -4.0 -2.0 0.0 2.0 4.0 0 0.5 1 1.5 2 2.5 A cce le ra ti on ( m /s /s ) Time (s) 弾性 粘弾性 弾塑性 粘弾塑性 -12.0 -7.0 -2.0 3.0 8.0 0 0.5 1 1.5 2 2.5 A cce le ra ti on ( m /s /s ) Time (s) 弾性 粘弾性 弾塑性 粘弾塑性 Fig. 15. リッカー・ウェーブレットの𝑓 を 10 Hz，PGA を(a) 0.05g，(b) 0.10g，(c) 0.25g， (d) 0.50g とした場合の，粘土の表層地盤におけ る地表での加速度波形. (a) PGA=0.05g (b) PGA=0.10g (c) PGA=0.25g (d) PGA=0.50g Fig. 16. リッカー・ウェーブレットの𝑓 を 10 Hz とし，PGA を(a) 0.05g，(b) 0.10g，(c) 0.25g， (d) 0.50g とした場合の，粘土の表層地盤にお け る 地 表 で の 加 速 度 フ ー リ エ 振 幅 ス ペ ク ト ル. 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0 5 10 15 20 25 Fo ur ie r s pe ct rum (m /s ) Frequency (Hz) 弾性 粘弾性 弾塑性 粘弾塑性 0.00 0.05 0.10 0.15 0.20 0.25 0 5 10 15 20 25 Fo uri er sp ec tr um (m /s ) Frequency (Hz) 弾性 粘弾性 弾塑性 粘弾塑性 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0 5 10 15 20 25 Fo uri er sp ec tr um (m /s ) Frequency (Hz) 弾性 粘弾性 弾塑性 粘弾塑性 0.00 0.20 0.40 0.60 0.80 1.00 1.20 0 5 10 15 20 25 Fo uri er sp ec tr um (m /s ) Frequency (Hz) 弾性 粘弾性 弾塑性 粘弾塑性 (a) PGA=0.05g (b) PGA=0.10g (c) PGA=0.25g (d) PGA=0.50g
うに振幅が小さくなり，反射波に遅れが生じてい るが，砂の場合(Fig. 12(b)～(d))と比較して明らか に振幅の減少が小さい．同様に，Fig. 16(b)～(d)で も減衰は大きくなっているが，どのケースでも砂 の場合(Fig. 12(b)～(d))ほどの減衰は生じていない． 次に，Fig. 17 に PGA が 0.25g のときの，最大 せん断ひずみ及び最大せん断応力の深さプロファ イルと，各構成モデルにおける深さ 5.1 m での応 力－ひずみ関係を示す．砂の場合(Fig. 13)と同様に， 深さ4 m 付近で最大せん断ひずみと最大せん断応 力が最大となっている．また，4 m 付近の最大せん 断ひずみは，砂の場合と異なり，弾性モデルよりも 弾塑性モデルの方がピークの最大値が大きくなっ ている．一方，最大せん断応力は，弾塑性モデルが 弾性モデルよりも小さくなっているが，砂の場合 ほどの減衰は見られず，粘弾性モデルより値が大 きくなっている．応力－ひずみ曲線を見ると，弾塑 性で砂の場合に比べて履歴曲線は細く，曲線に囲 -30 -25 -20 -15 -10 -5 0 0 0.02 0.04 0.06 De pt h ( m )
Maximum Shear strain (%) 弾性 粘弾性 弾塑性 粘弾塑性 -30 -25 -20 -15 -10 -5 0 0 10 20 30 De pt h ( m )
Maximum Shear stress (kPa) 弾性 粘弾性 弾塑性 粘弾塑性 -18 -12 -6 0 6 12 18 -0.06 -0.03 0 0.03 0.06 St re ss (k P a) Strain (%) -18 -12 -6 0 6 12 18 -0.06 -0.03 0 0.03 0.06 St re ss (k P a) Strain (%) -18 -12 -6 0 6 12 18 -0.06 -0.03 0 0.03 0.06 St re ss (k P a) Strain (%) -18 -12 -6 0 6 12 18 -0.06 -0.03 0 0.03 0.06 St re ss (k P a) Strain (%) Fig. 17. 表層地盤が粘土のときの，深さに対する 最大せん断ひずみと最大せん断応力と，地盤が① 弾性，②粘弾性，③弾塑性，④粘弾塑性モデルの 場合の深さ 5.1 m における応力－ひずみ関係 (𝑓 =10 Hz, PGA=0.25g の場合). ①弾性モデル ③弾塑性モデル ②粘弾性モデル ④粘弾塑性モデル 最大せん断応力 最大せん断ひずみ Fig. 18. 震央と観測点 OKYH07(神郷)の位置. Fig. 19. 入射加速度波形(上)と加速度フーリエ スペクトル(下).波形の時刻ゼロは P 波の到達時 間．
まれた面積も小さくなっている．これは，非線形性 による履歴減衰によって生じたエネルギー散逸は， 粘土の場合の方が砂の場合よりも小さいことを示 している． 4.4. ケース 4 ケース4 では，実波形による解析を行う．実波形 は ， 防 災 科 学 技 術 研 究 所 の KiK-net 観 測 点 OKYH07(神郷)で得られた 2000 年鳥取県西部地震 の加速度波形を解放基盤上の記録として用いる． 震央と観測点の位置をFig. 18 に示す．この地震は 2000 年 10 月 6 日 13 時 30 分，鳥取県西部で発生 し，マグニチュード𝑀 7.3 (𝑀 6.8)であった．震 央は北緯35 度 16 分 27 秒，東経 133 度 20 分 56 秒，深さは9 km と推定されている(気象庁一元化 震源による)．OKYH07 で観測された地表の水平加 速度記録を回転して得たtransverse 成分について， 振幅を半分にした波形と，そのフーリエ振幅スペ クトルをFig. 19 に示す．これを表層地盤モデルの 基盤への入射波とする． 実波形による解析に用いた地盤構造モデルを Fig. 20 に示す．深さ 30 m の表層地盤と半無限の 基盤の 2 層構造を仮定し，物性値は，S 波速度𝑉 (s) 150000 0.2 (m) 100673 Table 5. ケース4で用いた差分パラメーター
∆𝑡
∆𝑧
𝑁
𝑁
2.0 10 Fig. 20. 実波形による解析に用いる地盤構造モ デル. Fig. 21. 地表での加速度波形(上)とそのフーリ エ振幅スペクトル(下). ①弾性モデル ②弾塑性モデル Fig. 22. 深さに対する最大せん断ひずみ(左上) と最大せん断応力(右上)と，深さ 5.1 m における 応力－ひずみ関係(下). 最大せん断応力 最大せん断ひずみは表層で150 m/s，基盤で 510 m/s，密度𝜌 は表層 で1800 kg/m3，基盤で1900 kg/m3とした．基盤の 物性値は，観測点のPS 検層の最上層の値を参考に した．ここでは，構成モデルは弾性モデルと弾塑性 モデルのみ，地盤の土質材料については砂のみの 解析を行う．砂の𝐺 𝐺⁄ 𝛾 関係については，これ までのケースと同じデータを用いた．差分法計算 に使用した離散化パラメーターの値を Table 5 に 示す． Fig. 21 に実波形による解析で得られた，弾性モ デルと弾塑性モデルの地表加速度波形と，フーリ エ振幅スペクトルを示す．加速度波形は，明らかに 弾塑性モデルの方が弾性モデルよりも振幅が小さ い．このことは，特に4 秒付近から 16 秒付近の間 で顕著であり，非線形性による履歴減衰が生じて いることを示している．振幅スペクトルは，3 Hz 付近から高周波数側で大きな減衰が見られ，6 Hz 以上ではほぼ連続的に減衰している．さらに，13 Hz よりも高周波数側では，低周波数側へピークの シフトも見られる． 次に，Fig. 22 に最大せん断ひずみと最大せん断 応力の深さプロファイルと，両構成モデルにおけ る深さ 5.1 m での応力－ひずみ曲線を示す．最大 せん断ひずみは，弾性モデル，弾塑性モデルの両方 で，表層と基盤層のインピーダンス比により基盤 近くの深さ30 m 付近で大きく増幅しているが，弾 塑性モデルの方が増幅が大きくなっている．また， 地表から深さ5 m 付近までは弾性モデルの方が大 きくなっているが，深さ5 m 以深では弾塑性モデ ルの方が大きい．これらは非線形性により剛性率𝐺 が低下したことに起因している．最大せん断応力 は，弾性モデルが地表から基盤面まで弾塑性モデ ルよりも大きい．これは，弾塑性モデルで非線形性 による減衰が生じているためであり，弾塑性モデ ルの応力－ひずみ関係が履歴曲線を描いているこ とから明らかである． 5. 結論 以上の結果より，入射地震動が強いほど，そして 表層地盤が粘土よりも砂で構成されている方が， 非線形性が顕著に表れ，さらに，入射地震動が強く なるとともに，スペクトルのピークの減衰と低周 波数側へのシフトが現れ始める周波数も低くなる ことがわかった． 本研究では，非線形応答の計算を行うための非 線形逐次積分法の計算プログラムを作成し，それ を用いて解析を行った．入射波リッカー・ウェーブ レットの解放基盤におけるPGA を変化させた場合， PGA の増加に伴って表層地盤の非線形性による影 響が大きくなり，振幅スペクトルのピークの低周 波数側へのシフトと高周波数帯のスペクトルの低 下が顕著となることが確認できた．また，リッカ ー・ウェーブレットの中心周波数𝑓 を変化させた 場合，𝑓 の増加に伴って，粘弾性モデルでは，周波 数に依存する非弾性減衰により，全周波数帯で大 きな減衰が生じた．一方，弾塑性モデルでは𝑓 が 大きくなるほど非線形性の影響が小さくなった． また，表層地盤を粘土とした場合では，表層地盤を 砂とした場合に比べて非線形性は小さくなった． つまり，粘土地盤に比べ砂地盤の方が非線形性の 影響が現れやすいことが確認できた．また，2000 年鳥取県西部地震の実波形を用いて解析を行った 結果，リッカー・ウェーブレットによる解析と同様 に，非線形性の影響を示し，スペクトルのピークの 低周波数側へシフトと高周波数帯のスペクトルの 低下を確認できた．この結果は，強震動の評価にお いて非線形性を考慮することの必要性を示してい る．この実波形による解析では，2 層の単純な地盤 モデルで解析を行ったが，実際の複雑な地盤モデ ルを用いることで，より現実的な非線形性の影響 を予測することが可能であろう． 謝辞 本論文は第一著者の岡山大学理学部地球科学科 平成29 年度課題研究に基づいています．本研究で は防災科学技術研究所の KiK-net の強震波形記録 を使用しました．一部の作図は Generic Mapping Tools (Wessel and Smith, 1998)を使用しました． 記して感謝申し上げます．
引用文献
Bardet, J. P. and T. Tobita, 2001, NERA: A com-puter program for nonlinear earthquake site response analyses of layered soil deposits, De-partment of Civil Engineering, University of Southern California, Los Angeles,
https://sites.google.com/site/tt60898/home/sof tware.
Beresnev, I. A. and K. L. Wen, 1996, Nonlinear soil response: a reality?, Bulletin of the Seis-mological Society of America, 97, 2180-2196. Cerjan, C., D. Kosloff, R. Kosloff, and M. Reshef,
1985, A nonreflecting boundary condition for discrete acoustic and elastic wave equation, Geophysics, 50, 705-708.
Clayton, R. and B. Engquist, 1977, Absorbing boundary conditions for acoustic and elastic wave equations, Bulletin of the Seismological Society of America, 67, 1529-1540.
Graves, R. W., 1996, Simulating seismic wave propagation in 3D elastic media using stag-gered-grid finite differences, Bulletin of the Seismological Society of America, 86, 1091-1106.
林田智宏・竹中博士・岡元太郎, 1999, 速度－応力 型のスタガード格子差分法を用いた2 次元およ び3 次元地震波動計算コードの作成, 九州大学 理学部研究報告地球惑星科学, 20, 99-110. Iwan, W. D., 1967, On a class of models for the
yielding behavior of continuous and composite systems, Journal of Applied Mechanics, 34(3), 612-617.
Kausel, E. and D. Assimaki, 2002, Seismic simu-lation of inelastic soils via frequency-depend-ent moduli and damping, Journal of Engineer-ing Mechanics, 128(1), 34-47.
Levander A. R., 1988, Fourth-order finite-differ-ence P-SV seismograms, Geophysics, 53, 1425-1436.
Masing, G., 1926, Eigenspannungen und Verfes-tigung beim Messing, Proceedings of the Sec-ond International Congress of Applied Me-chanics, 332-335.
Masuda, T., S. Yasuda, N. Yoshida, and M. Sato, 2001, Field investigations and laboratory soil tests on heterogeneous nature of alluvial de-posits, Soil and Foundations, 41 (4), 1-16. Mróz, Z., 1967, On the description of anisotropic
workhardening, Journal of the Mechanics and Physics of Solids, 15, 163-175. 野口科子, 2009, 地震記録に基づく表層地盤の非 線形応答に関する研究, 北海道大学学位論文, pp.109. 野口科子・笹谷 努, 2011, 2003 年宮城県沖スラブ 内地震における表層地盤の非線形応答とその強 震動への影響, 地震第 2 輯, 63(3), 165-187.
Schnabel, P. B., J. Lysmer, and H. B. Seed, 1972, SHAKE A Computer program for earthquake
response analysis of horizontally layered sites, Report No. EERC72-12, University of Califor-nia Berkeley.
Seed, H. B. and I. M. Idriss, 1970, Soil moduli and damping Factors for dynamic response analyses, Report No. EERC70-10, Earth-quake Engineering Research Center, Univer-sity of California, Berkeley, pp.40.
Seed, H. B. and J. H. Sun, 1989, Implication of site effects in the Mexico City earthquake of September 19, 1985 for earthquake-re-sistance-design criteria in the San Francisco Bay Area of California, No. UCB/EERC-89/03, University of California, Berkeley, California, pp.140. 杉戸真太・会田尚義・増田民夫, 1994, 周波数特性 を考慮した等価ひずみによる地盤の地震応答解 析法に関する一考察, 土木学会論文集, No. 493/ Ⅲ-27, 49-58. 竹中博士, 1998, 第 16 章「シミュレーション」第 3.1 節「弾性波探査」, 物理探査ハンドブック, 物理探査学会, 851-862 (2016 年出版の物理探査 ハンドブック増補改訂版では938-949). 田中宏樹・竹中博士, 2005, 鉛直方向任意不均質弾 性媒質における平面波入射問題の時間領域差分 解法, 地震第 2 輯, 57, 343-354. 時松孝次・翠川三郎, 1988, 地表で観測された強震 記録から推定した表層地盤の非線形性状, 日本 建築学会論文報告集, 388, 131-137. 吉田 望, 1994, 実用プログラム SHAKE の適用性, 軟弱地盤における地震動増幅シンポジウム発表 論文集, 土質工学会, 14-31.
Yoshida, N., S. Kobayashi, I. Suetomi, and K. Miura, 2002, Equivalent linear method con-sidering frequency dependent characteristics of stiffness and damping, Soil Dynamics and Earthquake Engineering, 22, 205-222.
吉田 望, 2010, 地盤の地震応答解析, 鹿島出版会, pp. 272.
Virieux, J., 1984, SH-wave propagation in heter-ogeneous media: Velocity-stress finite-differ-ence method, Geophysics, 49, 1933-1957. Wessel, P. and W. H. F. Smith, 1998, New
im-proved version of the Generic Mapping Tools released, EOS Transaction American Geo-physical Union, 79, 579.
付録 IM モデルによるせん断応力計算コードの例 IM モ デ ル に よ り せ ん 断 応 力 𝜏 を 求 め る Fortran のサブルーチンを示す．コードは Fortran 77 形式で書かれており，引数は以下の通りである． nz ：深さ方向の全格子数． nbz ：地表から工学的基盤上面の深さまでの格子 数 nim ：バネの数(=スライダーの数) syz ：せん断応力𝜏 dgyz ：式(9)から計算されるひずみ増分∆𝛾 alp ：バックストレス𝛼 Ri ：式(13)から計算される降伏応力𝑅 Hi ：式(12)から計算される骨格曲線の傾き𝐻 ad ：式(14)から計算される非弾性減衰の項𝐴 ． 要素数 nz の 1 次元配列で，格子ごとに設 定．非弾性減衰がない場合の値は1.0．
