この内容は，私たちの論文¹⁾を分かり易く説明したものです。内容は，3つのコラムで構成されています。一つ目がSHAKEの減衰でここです。その他，Lysmerの提案とその問題，合理的な複素剛性に記述があります。

　SHAKE²⁾については，皆さんよくご存じと思います。SHAKEの開発は1972年で，もう，50年近く使われて来ていますが，以下に示すような基本的な所に疑問を持つ人がいなかったようです。また，文献2)にある式の誤りについても，指摘されたものを見たことがありません。なお，SHAKEには以下に示すような，基本的な理解に欠けるところがあるのですが，（意識したのかしていないのかわかりませんが）それを地盤の地震応答解析に使えるようにしたのは非常に大きな成果です。その業績にケチを付けるつもりは全くありません。ただ，ちゃんと理解しておく必要はあると思います。

　SHAKE¹⁾によると，運動方程式は次の様に書いてあります。

\(\rho \cfrac{\partial ^{2}u}{\partial t^2}=G\cfrac{\partial ^{2}u}{\partial x^2}+\eta\cfrac{\partial ^{3}u}{\partial x^2\partial t} \)(1)

ここで，\(u=u(x,t)\)は変位，\(\rho\)は密度，\(G\)はせん断剛性，\(\eta\)は減衰係数，\(x\)は鉛直下方を向く座標軸，\(t\)は時間です。この式の誘導についても少し手間が掛かるのですが，それは後の話にします。

　この式を解くのに，\(u(x,t)=U(x)e^{i\omega t}\)と変数分離し，空間に関する解を求めると

\(U(x)=Ee^{ikx}+Fe^{-ikx}\)(2)

の様になります。\(E\)と\(F\)が入射波と反射波の振幅です。ここで，

\(k^2=\cfrac{\rho \omega ^2}{G+i\omega \eta}=\cfrac{\rho \omega^2}{G^*}\)(3)

は（複素の）波数です。この式の中央部分\(G+i\omega \eta\)の複素数の部分\(\omega \eta\)を次の様に変換します。

\(\omega\eta=2G\beta\)(4)

と置くと，有名な複素剛性に対する定義

\(G^*=G(1+2i\beta)\)(5)

が得られます。問題は式(4)の\(\beta\)です。SHAKEのマニュアルでは，臨界減衰比（critical damping ratio）と書いてあります。臨界減衰比については減衰に関する用語でも説明しましたし，ほとんどの振動に関する本にも書いてあります。臨界減衰比を\(h\)と表すと，次の式になります。

\(h=\eta/(2\sqrt{mk})\)(6)

です。式(4)と式(6)の定義，全く違いますよね。式(4)には質量\(m\)は含まれていません。これは当然で，式(6)は1自由度系の振動から得られたものですし，式(4)は地盤の釣合式です。つまり，式(4)には臨界減衰（critical damping）に対応する量はありませんし，従って，比もありません。

　筆者は，これはSHAKEの開発者らの勘違いと思います。減衰に関する用語でも説明しましたように，\(\beta\)の値が小さいときには，履歴減衰により求められる減衰定数\(h\)と臨界減衰比\(\beta\)は同じになります。一方，式(4)で定義した\(\beta\)も，SHAKEのモデル化（式(5)）を使えば，減衰定数と同じになります。従って，減衰定数に対応する用語として\(\beta\)を臨界減衰比と書いたものと考えられます。

　ところで，減衰に関する用語では，1自由度系の非線形挙動の伴う履歴減衰が

\(h=\cfrac{1}{4\pi}\cfrac{S}{W}\)(7)

が成立する，すなわち，\(h=\beta\)となるのは，\(h\)が小さく，次の式が成り立つとき

\(e^{-h\pi}\approx1-h\pi\) (8)

という事を示しました。この式はマクローリン展開の第2項までをとったもので，その次の項は\(\pi h^2/2\)です。1に対してこの項が無視できるというのが，式(7)が成立する条件です。建物などでは\(h\)の最大値は5%程度ですから，\(\pi h^2/2=0.004\)つまり0.4%で無視するのに問題は無いと考えられます。しかし，地盤では\(h\)の最大値は30%程度ですので，\(\pi h^2/2=0.14\)つまり，14%も誤差があるわけで，無視できるとも思えません。この点もおかしいですね。

　土の応力－ひずみ関係の問題ですから，応力－ひずみ関係モデルから話を進めます。粘弾性のモデルとして，ダッシュポットとばねで構成されるモデルには，

の二つがあります。ここでは，Voigtモデルを使います。ここで，\(C\)は粘性係数，\(G\)はばね定数です。このモデルの応力－ひずみ関係は次の式(9)です。

\(\tau=G\gamma + C\dot{\gamma}\)(9)

　SHAKEでは，それぞれの振動数で調和振動（三角関数の外力を受け，定常振動している状態）しています。そこで，次の様にひずみを設定します。

\(\gamma=\gamma_0\cos{\omega t}\)　または　\(\cos{\omega t}=\gamma/\gamma_0\)(10)

これを式(9)に代入し，\(\sin{\omega t}=\pm\sqrt{1-\cos{^2\omega t}}\)の関係を用いると応力－ひずみ関係は次の様になります。

\(\tau=G\gamma_0\cos{\omega t}-C\gamma_0\omega\sin{\omega t}=G\gamma\pm C\omega\sqrt{\gamma_0^2-\gamma^2}\)(11)

この式は右辺に\(\omega\)が入っていますので，応力－ひずみ関係は振動数に依存していることになります。しかし，土では（地震応答解析に用いる程度の振動数では）振動数やひずみ速度に依存しないという性質がありますから，都合が悪いです。そこで，数学上のトリックとして，減衰係数が振動数に依存していると考えます。つまり，減衰係数の代わりに次の様なパラメータ\(\beta\)を導入します。

\(\omega C=2\beta G\)　または　\(C=2\beta G/\omega\)(12)

これを式(11)に代入すると，振動数に依存しない応力－ひずみ関係が次の様に得られます。

\(\tau=G\gamma\pm 2\beta G\sqrt{\gamma_0^2-\gamma^2}\)(13)

この式の第1項は直線，第2項は円の式で，それらを合わせたものは，

の様な楕円になります。

　ところで，式(12)と式(4)は同じです。つまり，式(4)で定義した\(\beta\)は，critical damping ratioというのではなく，Voigtモデルの応力－ひずみ関係が周波数に依存しないという設定で求まるものです。過去の文献を調べたのですが，これまでこういうアプローチはなかったようです。従って式(12)の\(\beta\)にも名前がありません。そこで，適切かどうかわからないのですが，私たちはこれを減衰パラメータと名付けました。

　これだけでは，SHAKEの減衰との関係が十分わかったと感じられない方もいらっしゃるかもしれません。そこで，もう少し議論を進めてみます。

　まず，式(12)から得られる減衰定数を求めます。最初に，履歴吸収によるエネルギー\(\Delta W\)を求めますが，これは式(12)の第2項（円の部分）の面積です。式のルートの中は半径\(\gamma_0\)の円ですので，面積は\(\pi\gamma_0^2\)ですから，面積は

\(\Delta W=2\beta G\pi\gamma_0^2\)(14)

となります。一方，ひずみエネルギー（弾性エネルギーと呼ばれることもあります）\(W\)は

\(W=G\gamma_0^2/2\)(15)

で計算できますので，減衰定数\(h\)は次の様になります。

\(h=\cfrac{1}{4\pi}\cfrac{\Delta W}{W}=\beta\)(16)

つまり，これまで使われてきたように，室内試験で得られた減衰定数\(h\)を使ってよいということがわかります。また，式(8)で示したような近似式も使っていませんので，減衰定数が小さいときに成立するというような縛りもありません。

　ついでに，もう少し詳しく，SHAKEで使われる応力－ひずみ関係を見てみましょう。SHAKEでは複素剛性を用いています。上の例ではVoigtモデルを使っていました。この様な問題に複素剛性を用いることを最初に提案したのは，Sorokin³⁾と思われます。ただ，色々調べたのですが，この論文はロシア語で書かれていることから見つけることができませんでした。この論文は多くの所で引用されています。ただし，いずれも年代が古く，著者とコンタクトすることができそうにありませんでした。最近では文献4)で引用されています。文献4)には著者のメールアドレスが書いてあったので，早速連絡を取りました。すると，彼は前職の時の図書館で見たので，現物は手元にないと言うことでした。そこで，その図書館に論文のコピーを御願いしたのですが，返事はいただけませんでした。このような複素剛性は，昔の論文を見るとSorokinの仮説⁵⁾とかSorokinの減衰⁶⁾と呼ばれているようです。

　Sorokinさんの方法に従って，応力－ひずみ関係を構築するのですが，そのためには複素数と実数とを表現上でも区別しておくのがよいと思います。そこで，複素数はオーバーバーを付けることにします。つまり複素数のひずみは\(\overline{\gamma}\)の様に表すわけです。これまで同様，定常状態の調和振動を考え，ひずみを次の様に与えます。

\(\overline{\gamma}=\gamma_0 e^{i\omega t}\)(17)

また，応力－ひずみ関係を次の様に設定します。

\(\overline{\tau}=G\overline{\gamma}+C\dot{\overline{\gamma}}\)(18)

　式(17)を式(18)に代入し，式(12)の関係を用いれば，次の応力－ひずみ関係が得られます。

\(\overline{\tau}=G(1+2i\beta)\gamma_0e^{i\omega t}=\overline{G_S^*}\overline{\gamma}\)(19)

ここで，

\(\overline{G_S^*}=G(1+2i\beta)\)(20)

がSHAKEで用いられている複素剛性です。ここで，下添え字に\(S\)を付けているのは，Sorokinモデルの複素剛性であることを意味しています。Sorokinモデルの履歴吸収エネルギー\(\Delta W\)は前と同様に求められ，

\(\Delta W_S=2G\beta\pi\gamma_0^2\)(21)

で，室内試験の履歴吸収エネルギーと同じです。つまり，SHAKEでは室内試験で得られた\(G\)と\(h\)の両方を満たすモデルになっています。ただし，式(19)や上の図からわかるように，最大せん断応力は\(\sqrt{1+4\beta^2}G\gamma_0\)で，室内試験の結果より大きくなっています。

3. SHAKEの基礎式

　SHAKE¹⁾は一次元の地盤の解析を行っています。変位を\(\overline{u}=\overline{u}(z,t)\)とし，質量密度を\(\rho\)，減衰係数を\(C\)とすれば，運動方程式は，次のように書けるとしています。

\(\rho\cfrac{\partial^2\overline{u}}{\partial t^2}=G\cfrac{\partial^2\overline{u}}{\partial z^2}+C\cfrac{\partial^2\overline{u}}{\partial t\partial z^2}\)(22)

ここで，\(z\)は鉛直下方を向く座標です。変位を変数分離し

\(\overline{u}(z,t)=\overline{U}(z)e^{i\omega t}\)(23)

と置くと，式(22)は次のようになります。

\((G+i\omega C)\cfrac{\partial^2\overline{U}}{\partial z^2}=-\rho\omega^2\overline{U}(z)\)(24)

なお，SHAKEのマニュアル²⁾には右辺の負号はありません。これは虚数単位（\(i\)）の二乗が-1になるというところを間違えたものです。SHAKEではもう一度同じ間違いをしているので，結果として得られた支配方程式は正しいのですが。

　ここで式(12)の関係を用いて，複素剛性\(\overline{G}_S^*\)を式(20)のように置けば，周波数に依存しない式が得られます。

　なお，式(22)は次のように導きます。まず，応力－ひずみ関係は式(9)示されており，これに式(12)の関係とひずみ－変位関係式\(\gamma=\partial u/\partial z\)に代入すれば，次の式となります。

\(\tau=G\cfrac{\partial t}{\partial z}+\cfrac{2\beta G}{\omega}\cfrac{\partial^2u}{\partial t\partial z}\)(25)

一方，一次元地盤の微小要素に水平方向に作用する力の釣合式は次式で表されます。

\(\cfrac{\partial u}{\partial z}=\rho\cfrac{\partial^2u}{\partial t^2}\)(26)

式(26)に式(25)を\(z\)で偏微分して代入すれば，式(22)が得られます。式(25)ではSorokinモデルによる応力－ひずみ関係を用いているので，SHAKEが用いているのはSorokinモデルに対応した波動方程式であることがわかります。

　なお，現在のSHAKEをはじめ多くの等価線形地震応答解析プログラムでは，Sorokinモデルは使っていません。それは著者の一人であるLysmerさんがこちらの方が良いという新しい複素剛性を示し⁷⁾，California大学からSHAKEを購入した人に，SHAKEを修正する依頼文を出したので，修正されているわけです。しかし，この修正は力学的にはおかしいところもありますし，入力した減衰定数が正しく減衰定数として使われていないという問題があります。これについては項を改め，Lysmerの提案とその問題で説明します。

1) 吉田望，安達健司：地盤の地震応答解析のための複素剛性法，日本地震学会論文集（投稿中）
2) Schnabel, P. B., Lysmer, J. and Seed, H. B. (1972): SHAKE A Computer program for earthquake response analysis of horizontally layered sites, Report No. EERC72-12, University of California, Berkeley
3) Sorokin, F., Internal and external friction by vibrations of rigid bodies, CNIISK 162, Moscow, 1957
4) Tesar, A.: Tuned vibration control in aeroelasticity of slender wood bridges, Coupled Systems Mechanics, Vol. 1, No. 3 pp. 219-234，2012
5) 小西一郎，高岡宣喜：構造動力学，丸善，295pp.，1973
6) Hricko, J.: Modelling compliant mechanisms - comparison of models in MATLAN/SimMechanics vs. FEM, Proc., 21st International Workshop on Robotics in Alpe-Adria-Danube Region, pp. 57-62, 2012
7) Lysmer, J.: Modal damping and complex stiffness, University of California Lecture note, University of California, Berkeley, 1973