JP4850132B2 - 時刻歴応答解析方法、装置及びプログラム - Google Patents

Description

本発明は時刻歴応答解析方法、装置及びプログラムに係り、特に、減衰定数ｈが振動数非依存特性を示しかつ剛性低下率α及び減衰定数ｈが歪振幅依存特性を示す物体の時刻歴応答解析を行う時刻歴応答解析方法、当該時刻歴応答解析方法を適用可能な時刻歴応答解析装置、及び、コンピュータを時刻歴応答解析装置として機能させるための時刻歴応答解析プログラムに関する。

地震時に構造物に入力される地震動は、構造物が建設されている地盤の特性（地震時の地盤の挙動）によって大きく相違する。このため、地震応答解析によって地震時の構造物の挙動や耐震安全性等を正確に解析・評価するためには、地震時の地盤の挙動も併せて解析・評価することが望ましい。地震応答解析は、周波数領域で応答解析を行う周波数応答解析と、時間領域で応答解析を行う時刻歴応答解析とに大別される。地盤の動的剛性を振動数非依存の周波数領域の複素関数として表し、これを直接用いた周波数応答解析も多用されている（例えば公知のＳＨＡＫＥを用いた周波数領域での等価線形解析）。

しかし、大きなエネルギーが入力される大地震時等には、そのエネルギーの一部が、構造物を構成する各部材の内部に亀裂を生じさせたり各部材を部分的に塑性化させる等によって消費されると共に、この亀裂発生や部分的な塑性化等に伴い各部材の強度が低下することが繰り返されるというプロセスを経るため、構造物の挙動は非線形性を示す。このため、地震時の構造物の挙動等を高精度に予測解析するためには、時刻歴応答解析により、地震時の各時点での構造物の状態（過去にどのような力が加わり、その力によってどのような状態になっているか）を考慮して各時点での構造物の挙動を解析する必要がある。地盤の挙動を時刻歴応答解析によって解析する際に適用可能なモデルとしてはRamberg-Osgoodモデル（以下Ｒ−Ｏモデル）や双曲線モデルが知られている。

上記に関連して本願出願人は、振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)−(２ω／ωm)なる値(但しｎは整数)を示す虚数部の正則成分、及び、振動数ωに拘わらず(２ω／ωm)なる値を示す虚数部の特異成分の和で表され、振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)なる値を示す虚数部と、前記虚数部の正則成分よりヒルベルト変換を用いて算出した因果律を満たす実数部と、から成る因果的単位虚数関数を定義し、当該因果的単位虚数関数を時間領域へ変換することで因果的単位虚数関数のインパルス応答値を演算し、この因果的単位虚数関数のインパルス応答値を用いることで、減衰定数ｈが振動数非依存特性を示す物体の時刻歴応答解析を行う技術を提案している（特許文献１を参照）。なお、ωは一般に円振動数、角振動数などと称されるが、本明細書では振動数と略称する。
特開２００７−０７１７５４号公報

地盤を含む多くの材料は、大地震時等のように大きなエネルギーが大振幅の外力として入力された場合、内部減衰の履歴吸収エネルギー(減衰定数ｈ)が振動数ωにあまり依存しない振動数非依存特性（減衰定数ｈが振動数ωに拘わらず一定値を示す特性）を示すことが知られており、更に、剛性低下率α（＝剛性Ｇ／初期剛性Ｇ_０）及び減衰定数ｈが物体の歪レベルγに応じて変化する歪振幅依存特性も示すことが知られている。

これに対し、前述のＳＨＡＫＥを用いた解析では、解析対象の物体の減衰定数ｈの振動数特性に関しては任意の特性（振動数非依存特性や振動数ωの変化に対して減衰定数ｈが所定の変化を示す特性）で解析を行うことができるものの、当該解析は周波数領域での解析であるため、解析対象の物体の歪レベルγに関しては大エネルギーの入力に拘わらず一定と見做して解析せざるを得ず、解析対象の物体の剛性低下率α及び減衰定数ｈが歪振幅依存特性を示していたとしても、これを考慮した解析を行うことはできない。

また、前述のＲ−Ｏモデルや双曲線モデルは、解析対象の物体の剛性低下率α及び減衰定数ｈの歪振幅依存特性を表すことは可能であるものの、減衰定数ｈの振動数特性については規定できないので、減衰定数ｈの振動数特性についてどのような特性で解析を行っているのか不明確であり、減衰定数ｈの振動数特性として任意の特性で解析を行うことは不可能である。また、Ｒ−Ｏモデルや双曲線モデルは、剛性低下率α及び減衰定数ｈの歪振幅依存特性についても明確に規定できるものではなく、チャート上に物体の歪と応力の関係を表すループ状の軌跡を描画することで解析結果を出力するにあたり、ループ状の軌跡の動きを規定するルールを設定しておくことで、剛性低下率α及び減衰定数ｈの歪振幅依存特性を間接的に表すことができるに過ぎず、物体の歪振幅依存特性を上記のルールへ置き換えることは試行錯誤的な作業であると共に、置き換えたルールの妥当性を確認する手法も確立されていないので、ルールの設定に多大な手間が掛り解析精度も低いという問題がある。

また、特許文献１に記載の技術を適用すれば、減衰定数ｈが振動数非依存特性を示す物体の時刻歴応答解析を行うことができるものの、この技術は剛性低下率α及び減衰定数ｈの歪振幅依存特性について考慮しておらず、物体の歪レベルの変化に拘わらず剛性低下率α及び減衰定数ｈを一定とみなして解析を行うので、歪振幅依特性を示す物体の時刻歴応答解析を精度良く行うことは困難である。

本発明は上記事実を考慮して成されたもので、減衰定数ｈが振動数非依存特性を示しかつ剛性低下率α及び減衰定数ｈが歪振幅依存特性を示す物体の時刻歴応答解析を容易かつ高精度に行うことができる時刻歴応答解析方法、時刻歴応答解析装置及び時刻歴応答解析プログラムを得ることが目的である。

上記目的を達成するために請求項１記載の発明に係る時刻歴応答解析方法は、減衰定数ｈが物体を振動させる外力の振動数に依存しない振動数非依存特性を示しかつ剛性低下率α及び減衰定数ｈが前記物体の歪レベルγに応じて変化する歪振幅依存特性を示す前記物体の時刻歴応答解析を行う時刻歴応答解析方法であって、振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)−(２ω／ωm)なる値(但しｎは整数)を示す虚数部の正則成分、及び、振動数ωに拘わらず(２ω／ωm)なる値を示す虚数部の特異成分の和で表され、振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)なる値を示す虚数部と、前記虚数部の正則成分のヒルベルト変換値に対応する実数部と、から成る因果的単位虚数関数を時間領域へ変換するか、又は前記虚数部のみを時間領域へ変換することで、前記因果的単位虚数関数のインパルス応答値として、前記物体の速度に依存する同時成分ｃ(t₀)、前記物体の変位に依存する同時成分ｋ(t₀)、前記物体の変位に依存するΔt刻みの時間遅れ成分ｋ(t_j)（但しｊは自然数でt_j＝Δt・j）を演算し、前記物体の質量マトリクスを[Ｍs]、前記物体の初期剛性マトリクスを[Ｋ₀]、時間領域での物体の変位ベクトルを{ｕ(t)}、速度ベクトルを{ｕ'(t)}、反力ベクトルを{Ｆ(γ,t)}、前記物体を振動させる外力の時間領域での加速度をｙ_０"(t)、時間遅れ成分ｋ(t_j)の総数をｎとしたときに、前記演算したインパルス応答値を、
[Ｍs]{ｕ"(t)}＋{Ｆ(γ,t)}＝−ｙ_０"(t)[Ｍs]{1} …(１)
但し、

上式に代入し、剛性低下率α(γ)及び減衰定数ｈ(γ)を各時刻における前記物体の歪レベルγに応じた値に切り替えながら物体の反力ベクトル{Ｆ(γ,t)}、変位ベクトル{ｕ(t)}及び速度ベクトル{ｕ'(t)}をΔt刻みで順次演算することで、前記物体の時刻歴応答解析を行うことを特徴としている。

減衰定数ｈが振動数非依存特性を示す物体の減衰（履歴減衰）を表す履歴減衰モデル、すなわち履歴減衰を有する複素剛性Ｓ(ω)は、次の(３)式で表される。なお、(３)式においてＫ_０は剛性、ｈは減衰定数、ｉは虚数単位である。
Ｓ(ω)＝Ｋ_０(１＋２ｈ・ｉ) …(３)
ここで、(３)式における虚数単位ｉに代えて次の(４)式で表される単位虚数関数Ｚ(ω)を用いれば、先の(３)式は(６)式で表される。
Ｚ(ω)＝Ｚ_Ｒ(ω)＋Ｚ_Ｉ(ω)・ｉ …(４)
なお、Ｚ_Ｒ(ω)は単位虚数関数Ｚ(ω)の実数部、Ｚ_Ｉ(ω)は単位虚数関数Ｚ(ω)の虚数部であり、それぞれ以下の(５)式で表される。

Ｓ(ω)＝Ｋ_０(１＋２ｈ・Ｚ(ω)) …(６)
但し、上記(５)式によって規定される単位虚数関数Ｚ(ω)は、実数部Ｚ_Ｒ(ω)と虚数部Ｚ_Ｉ(ω)がヒルベルト(Hilbert)変換対を形成しないために因果律を満たさず、時間領域への厳密な変換は不可能である。

これに対して請求項１記載の発明では、振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)−(２ω／ωm)なる値(但しｎは整数)を示す虚数部の正則成分Ｚ'_ＩＲ(ω)(図１(Ｂ)参照)、及び、振動数ωに拘わらず(２ω／ωm)なる値を示す虚数部の特異成分Ｚ'_ＩＳ(ω)(図１(Ｃ)参照)の和で表され、振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)なる値を示す虚数部Ｚ'_Ｉ(ω)と、虚数部の正則成分Ｚ'_ＩＲ(ω)のヒルベルト変換値に対応する（すなわち因果律を満たす）実数部Ｚ'_Ｒ(ω)と、から成る因果的単位虚数関数Ｚ'(ω)を用いている。因果的単位虚数関数Ｚ'(ω)の虚数部Ｚ'_Ｉ(ω)は、図２に示すように振動数範囲(ω＞ωm)及び(ω＜−ωm)で単位虚数関数Ｚ(ω)の虚数部Ｚ_Ｉ(ω)の値と相違する((５)式の条件から外れる)ものの、図１(Ａ)に示すように、振動数範囲（０＜ω＜ωm）内で＋１の値を示すと共に、振動数範囲（−ωm＜ω＜０）内で−１の値を示すので、振動数範囲（−ωm＜ω＜ωm）内では単位虚数関数Ｚ(ω) の虚数部Ｚ_Ｉ(ω)と同一の値を示す(虚数部Ｚ_Ｉ(ω)に関する(５)式の条件を満たす)。

また、因果的単位虚数関数Ｚ'(ω)の虚数部の正則成分Ｚ'_ＩＲ(ω)と特異成分Ｚ'_ＩＳ(ω)のうち、正則成分Ｚ'_ＩＲ(ω)は因果律を満たさない成分、特異成分Ｚ'_ＩＳ(ω)は因果律を満たす成分であるが、請求項１記載の発明では、因果的単位虚数関数Ｚ'(ω)の実数部Ｚ'_Ｒ(ω)を、虚数部の正則成分Ｚ'_ＩＲ(ω)のヒルベルト変換値に対応する値（すなわち因果律を満たす値）とすることで、一定の振動数範囲(−ωm＜ω＜ωm)内で単位虚数関数Ｚ(ω)と同様の値を示し、かつ因果律を満たす因果的単位虚数関数Ｚ'(ω)を得ている。なお、虚数部の正則成分Ｚ'_ＩＲ(ω)のヒルベルト変換値に対応する（すなわち因果律を満たす）実数部Ｚ'_Ｒ(ω)は、次の(７)式（ヒルベルト変換）を用いることで算出することができる。

一例として振動数ωm＝20Hz（40π）の場合に(７)式によって算出される因果的単位虚数関数Ｚ'(ω)の実数部Ｚ'_Ｒ(ω)を、因果的単位虚数関数Ｚ'(ω)の虚数部Ｚ'_Ｉ(ω)と共に図３に示す。図３からも明らかなように、因果的単位虚数関数Ｚ'(ω)の実数部Ｚ'_Ｒ(ω)は、単位虚数関数Ｚ(ω)の実数部Ｚ_Ｒ(ω)のように０ではなく(実数部Ｚ_Ｒ(ω)に関する(５)式の条件を満たしておらず)、振動数ωに依存した値となっている点で単位虚数関数Ｚ(ω)と相違しているが、因果的単位虚数関数Ｚ'(ω)は、実数部Ｚ'_Ｒ(ω)を上記のように定めることで実数部Ｚ'_Ｒ(ω)と虚数部Ｚ'_Ｉ(ω)がヒルベルト変換対を形成するために因果律を満たし、時間領域へ厳密に変換することが可能となる。また、因果的単位虚数関数Ｚ'(ω)の虚数部Ｚ'_Ｉ(ω)は、前述のように振動数範囲（−ωm＜ω＜ωm）内で単位虚数関数Ｚ(ω)の虚数部Ｚ_Ｉ(ω)と値が一致している。従って、(６)式における単位虚数関数Ｚ(ω)に代えて因果的単位虚数関数Ｚ'(ω)を用いる（因果的単位虚数関数Ｚ'(ω)を２ｈ倍して用いる）ことにより、振動数範囲（−ωm＜ω＜ωm）内で履歴減衰モデルと精度良く一致し（振動数範囲（−ωm＜ω＜ωm）内で物体の減衰定数ｈの振動数非依存性を精度良く表し）かつ因果律を満たす（時間領域への厳密な変換が可能な）減衰モデルを得ることができる。

また、請求項１記載の発明では、上記の因果的単位虚数関数Ｚ'(ω)を時間領域へ変換することで、因果的単位虚数関数Ｚ'(ω)のインパルス応答値として、物体の速度に依存する同時成分ｃ(t₀)、物体の変位に依存する同時成分ｋ(t₀)、物体の変位に依存するΔt刻みの時間遅れ成分ｋ(t_j)（但しｊは自然数でt_j＝Δt・j）を演算する。なお、「同時成分」は現在の状態量（変位・速度・加速度）に依存して反力を生じる量(成分)を、「時間遅れ成分」は過去の状態量に依存して反力を生じる量(成分)を意味している。

なお、因果的単位虚数関数Ｚ'(ω)のインパルス応答値は、因果的単位虚数関数Ｚ'(ω)の実数部Ｚ'_Ｒ(ω)をヒルベルト変換を用いて算出する場合は、因果的単位虚数関数Ｚ'(ω)を時間領域へ変換することで演算することができるが、請求項１記載の発明はこれに限られるものではなく、因果的単位虚数関数Ｚ'(ω)の実数部Ｚ'_Ｒ(ω)を算出することなく、因果的単位虚数関数Ｚ'(ω)の虚数部Ｚ'_Ｉ(ω)のみを時間領域へ変換することで因果的単位虚数関数Ｚ'(ω)のインパルス応答値を演算することも可能である。

すなわち、時間領域で因果律を満たす複素関数の実数部及び虚数部は周波数領域でヒルベルト変換対を形成し、複素関数の実数部及び虚数部のうちの虚数部のみが既知の場合にも、ヒルベルト変換により、当該虚数部から当該虚数部に対応しかつ時間領域で因果律を満たす実数部を一意に算出できるので、虚数部には実数部の成分も含まれていると見なすことができる。また、時間領域で因果律を満たさない（ヒルベルト変換対を形成しない）複素関数についても、虚数部が既知であれば、未知の実数部も近似によっておおよそ求まるので、時間領域で因果律を満たす複素関数と同様に、虚数部には実数部の情報が含まれていると見なすことができる。

本願発明者は上記事実に事実に着目し、特許第３８７８６２６号で提案したインパルス応答演算方法を基礎として、物体の動的剛性のうちの虚数成分（又は実数成分）のみを用いて時間領域変換（インパルス応答の演算）を行うことに想到し、物体の動的剛性Ｓ(ω)の実数成分Ｓ_Ｒ(ω)及び虚数成分Ｓ_Ｉ(ω)のうち虚数成分Ｓ_Ｉ(ω)のみを用いて時間領域変換（インパルス応答の演算）を行う数式として以下の(８)式を定め、

(６)式で表される関係を、動的剛性のＮ個の虚数部のデータ(Ｄ_Ｉ(ω_１)〜Ｄ_Ｉ(ω_Ｎ))のデータに対してマトリクス表示した、以下の(９),(10)式を導出し、

上記の(８)〜(10)式の動的剛性の虚数成分のみを用いたインパルス応答の演算（虚部変換法）を提案している。そして本願発明者は、この虚部変換法の有効性を確認する解析検討を行い、実用上十分な精度でインパルス応答を算出できることを確認している（中村尚弘，「実部もしくは虚部のデータのみを用いた複素剛性の時間領域変換法」，日本建築学会構造系論文集，２００７年２月，第６１２号，ｐ．７９−８６、及び、特願２００６−３３９９１３号）。従って、請求項１記載の発明において、上記の虚部変換法を適用し、因果的単位虚数関数Ｚ'(ω)の実数部Ｚ'_Ｒ(ω)を算出することなく、因果的単位虚数関数Ｚ'(ω)の虚数部Ｚ'_Ｉ(ω)のみを時間領域へ変換することで因果的単位虚数関数Ｚ'(ω)のインパルス応答値を演算するようにしてもよい。

一方、減衰定数ｈが振動数非依存特性を示す物体の剛性マトリクス[Ｋ(ω)]は、因果的単位虚数関数Ｚ'(ω)によって履歴減衰を表しかつ因果律を満たす因果的履歴減衰モデルを用いて次の(11)式で表される。なお、[Ｋ_０]は初期剛性マトリクスである。
[Ｋ(ω)]＝[Ｋ_０](１＋２ｈ・Ｚ'(ω)) …(11)
ここで、物体の剛性低下率α及び減衰定数ｈが物体の歪レベルγに応じて変化する歪振幅依存特性を示す場合、この歪振幅依存特性により、物体の剛性マトリクスも物体の歪レベルγに応じて変化する（次の(12)式を参照）。
[Ｋ(γ,ω)]＝α(γ)[Ｋ_０](１＋２ｈ(γ)・Ｚ'(ω)) …(12)
なお、(12)式において、[Ｋ(γ,ω)]は物体の歪レベルγ(及び振動数ω)に依存して変化する物体の剛性マトリクス、α(γ)は歪レベルγに応じて変化する剛性低下率、ｈ(γ)は歪レベルγに応じて変化する減衰定数を各々表す。

(12)式の両辺に変位{ｕ(ω)}を各々乗ずると次の(13)式が得られる。
{Ｆ(γ,ω)}＝[Ｋ(γ,ω)]・{ｕ(ω)}＝{Ｆ_１(γ,ω)}＋{Ｆ_２(γ,ω)} …(13)
但し、{Ｆ_１(γ,ω)}＝α(γ)[Ｋ_０]・{ｕ(ω)}
{Ｆ_２(γ,ω)}＝α(γ)[Ｋ_０]・２ｈ(γ)・Ｚ'(ω)・{ｕ(ω)}
上記の(13)式を時間領域へ変換すると次の(14)式が得られる。

なお、{ｕ(t)}は時間領域での物体の変位ベクトル、{ｕ'(t)}は時間領域での物体の速度ベクトル、{Ｆ(γ,ｔ)}は物体の反力ベクトル、ｎは時間遅れ成分ｋ(t_j)の総数である。そして上記の(14)式を整理することで請求項１記載の発明に係る前出の(２)式を得ることができる。なお、前出の(１)式は因果的単位虚数関数Ｚ'(ω)のインパルス応答値を用いて物体の時刻歴応答解析を行うための時間領域の運動方程式である。

請求項１記載の発明に係る前出の(１),(２)式には、物体の剛性低下率αの歪振幅依存特性がα(γ)として明示的に表現されていると共に、物体の減衰定数ｈの歪振幅依存特性がｈ(γ)として明示的に表現されている。そして請求項１記載の発明では、演算したインパルス応答値を前出の(１),(２)式に代入し、剛性低下率α(γ)及び減衰定数ｈ(γ)を各時刻における物体の歪レベルγに応じた値に切り替えながら物体の反力ベクトル{Ｆ(γ,t)}、変位ベクトル{ｕ(t)}及び速度ベクトル{ｕ'(t)}をΔt刻みで順次演算することで、物体の時刻歴応答解析を行うので、物体の時刻歴応答解析を行うにあたり、前記物体の減衰定数ｈが振動数非依存特性を示すこと、及び、前記物体の剛性低下率α及び減衰定数ｈが歪振幅依存特性を示すことを考慮した高精度の解析を行うことができる。また、請求項１記載の発明に係る (１),(２)式には物体の剛性低下率αの歪振幅依存特性がα(γ)として、減衰定数ｈの歪振幅依存特性がｈ(γ)として明示的に組み込まれているので、剛性低下率α及び減衰定数ｈの歪振幅依存特性を、妥当性の確認が困難な別のルールに置き換える等の煩雑な作業が不要となり、剛性低下率α及び減衰定数ｈの歪振幅依存特性として任意の特性を適用することも可能となる。

従って、請求項１記載の発明によれば、減衰定数ｈが振動数非依存特性を示しかつ剛性低下率α及び減衰定数ｈが歪振幅依存特性を示す物体の時刻歴応答解析を容易かつ高精度に行うことができる。

地震応答解析等の時刻歴応答解析では、構造物に対する解析対象の振動数範囲の多くがおおよそ定まっていることから、解析対象の振動数範囲が振動数ω＝０〜ωmの範囲内に入るように解析対象上限振動数ωmを予め固定的に定めておくことも可能であるが、解析対象の振動数範囲のうち下限振動数付近及び上限振動数付近では精度が低下することが本願発明者によって確認されており、これを考慮すると、請求項１記載の発明において、例えば請求項２に記載したように、物体の主要な固有振動数を把握する実固有値解析を行い、実固有値解析によって把握した前記物体の主要な固有振動数が、振動数ω＝０〜ωmの範囲内のうち誤差（例えば演算したインパルス応答値を用いて複素剛性Ｓ(ω)等を再現したときの目標値との比率）が所定値未満となる振動数範囲内に入るように、解析対象上限振動数ωmを設定又は選択することが好ましい。これにより、時刻歴応答解析の精度を更に向上させることができる。なお、上記の実固有値解析は、主要な固有振動数が何れの振動数範囲内に存在しているかを把握できれば目的を達成できるので簡単な処理で済む。

また、請求項２記載の発明において、例えば請求項３に記載したように、物体を振動させる外力と物体の挙動との関係の非線形化によって物体の固有振動数が変化するか否かを推定し、物体の固有振動数が変化すると判断した場合には、概略の減衰に基づく予備解析により非線形化後の固有振動数をおおよそ把握し、把握した非線形化後の固有振動数も振動数ω＝０〜ωmの範囲内のうち誤差が所定値未満となる振動数範囲内に入るように解析対象上限振動数ωmを設定又は選択することが好ましい。これにより、時刻歴応答解析の途中で、物体を振動させる外力と物体の挙動との関係の非線形化によって物体の固有振動数が変化する解析結果が得られたとしても、変化後の固有振動数が、誤差が所定値未満となる振動数範囲から外れてしまうことを防止することができ、時刻歴応答解析の精度を更に向上させることができると共に、解析対象上限振動数ωmを試行錯誤的に変更しながら解析を繰り返すことを回避できることで時刻歴応答解析の処理時間を短縮することができる。

また、解析対象の物体（の主要な固有振動数や非線形化後の固有振動数）に応じて解析対象上限振動数ωmを設定又は選択する態様であっても、解析対象上限振動数ωmの値の種類数は限られていることを考慮すると、請求項２又は請求項３記載の発明において、例えば請求項４に記載したように、解析対象上限振動数ωmの値が互いに異なる複数種の因果的単位虚数関数について、時間領域への変換を各々行い複数種の因果的単位虚数関数のインパルス応答値を各々演算して記憶手段に記憶しておき、物体の時刻歴応答解析に際し、記憶手段に記憶した複数種の因果的単位虚数関数のインパルス応答値のうち、把握した固有振動数が振動数ω＝０〜ωmの範囲内のうち誤差が所定値未満となる振動数範囲内に入る特定の因果的単位虚数関数のインパルス応答値を読み出して用いることが好ましい。これにより、時刻歴応答解析の度に因果的単位虚数関数のインパルス応答値を演算する必要が無くなるので、時刻歴応答解析の処理時間を短縮することができる。

請求項５記載の発明に係る時刻歴応答解析装置は、減衰定数ｈが物体を振動させる外力の振動数に依存しない振動数非依存特性を示しかつ剛性低下率α及び減衰定数ｈが前記物体の歪レベルγに応じて変化する歪振幅依存特性を示す前記物体の時刻歴応答解析を行う時刻歴応答解析装置であって、振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)−(２ω／ωm)なる値(但しｎは整数)を示す虚数部の正則成分、及び、振動数ωに拘わらず(２ω／ωm)なる値を示す虚数部の特異成分の和で表され、振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)なる値を示す虚数部と、前記虚数部の正則成分のヒルベルト変換値に対応する実数部と、から成る因果的単位虚数関数を時間領域へ変換するか、又は前記虚数部のみを時間領域へ変換することで、前記因果的単位虚数関数のインパルス応答値として演算された、前記物体の速度に依存する同時成分ｃ(t₀)、前記物体の変位に依存する同時成分ｋ(t₀)、前記物体の変位に依存するΔt刻みの時間遅れ成分ｋ(t_j)（但しｊは自然数でt_j＝Δt・j）を記憶する記憶手段と、前記物体の質量マトリクスを[Ｍs]、前記物体の初期剛性マトリクスを[Ｋ₀]、時間領域での物体の変位ベクトルを{ｕ(t)}、速度ベクトルを{ｕ'(t)}、反力ベクトルを{Ｆ(γ,t)}、前記物体を振動させる外力の時間領域での加速度をｙ_０"(t)、時間遅れ成分ｋ(t_j)の総数をｎとしたときに、前記記憶手段に記憶されているインパルス応答値を、
[Ｍs]{ｕ"(t)}＋{Ｆ(γ,t)}＝−ｙ_０"(t)[Ｍs]{1} …(１)
但し、

上式に代入し、剛性低下率α(γ)及び減衰定数ｈ(γ)を各時刻における前記物体の歪レベルγに応じた値に切り替えながら物体の反力ベクトル{Ｆ(γ,t)}、変位ベクトル{ｕ(t)}及び速度ベクトル{ｕ'(t)}をΔt刻みで順次演算することで、前記物体の時刻歴応答解析を行う解析手段と、を備えたことを特徴としているので、請求項１記載の発明と同様に、減衰定数ｈが振動数非依存特性を示しかつ剛性低下率α及び減衰定数ｈが歪振幅依存特性を示す物体の時刻歴応答解析を容易かつ高精度に行うことができる。

請求項６記載の発明に係る時刻歴応答解析プログラムは、記憶手段と接続されたコンピュータを、減衰定数ｈが物体を振動させる外力の振動数に依存しない振動数非依存特性を示しかつ剛性低下率α及び減衰定数ｈが前記物体の歪レベルγに応じて変化する歪振幅依存特性を示す前記物体の時刻歴応答解析を行う時刻歴応答解析装置として機能させるための時刻歴応答解析プログラムであって、前記記憶手段には、振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)−(２ω／ωm)なる値(但しｎは整数)を示す虚数部の正則成分、及び、振動数ωに拘わらず(２ω／ωm)なる値を示す虚数部の特異成分の和で表され、振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)なる値を示す虚数部と、前記虚数部の正則成分のヒルベルト変換値に対応する実数部と、から成る因果的単位虚数関数を時間領域へ変換するか、又は前記虚数部のみを時間領域へ変換することで、前記因果的単位虚数関数のインパルス応答値として演算された、前記物体の速度に依存する同時成分ｃ(t₀)、前記物体の変位に依存する同時成分ｋ(t₀)、前記物体の変位に依存するΔt刻みの時間遅れ成分ｋ(t_j)（但しｊは自然数でt_j＝Δt・j）が記憶されており、前記コンピュータを、前記物体の質量マトリクスを[Ｍs]、前記物体の初期剛性マトリクスを[Ｋ₀]、時間領域での物体の変位ベクトルを{ｕ(t)}、速度ベクトルを{ｕ'(t)}、反力ベクトルを{Ｆ(γ,t)}、前記物体を振動させる外力の時間領域での加速度をｙ_０"(t)、時間遅れ成分ｋ(t_j)の総数をｎとしたときに、前記演算手段によって演算されたインパルス応答値を、
[Ｍs]{ｕ"(t)}＋{Ｆ(γ,t)}＝−ｙ_０"(t)[Ｍs]{1} …(１)
但し、

上式に代入し、剛性低下率α(γ)及び減衰定数ｈ(γ)を各時刻における前記物体の歪レベルγに応じた値に切り替えながら物体の反力ベクトル{Ｆ(γ,t)}、変位ベクトル{ｕ(t)}及び速度ベクトル{ｕ'(t)}をΔt刻みで順次演算することで、前記物体の時刻歴応答解析を行う解析手段として機能させることを特徴としている。

請求項６記載の発明に係る時刻歴応答解析プログラムは、上記の記憶手段と接続されたコンピュータを上記の解析手段として機能させるためのプログラムであるので、上記のコンピュータが請求項６記載の発明に係る時刻歴応答解析プログラムを実行することで、上記のコンピュータが請求項５に記載の時刻歴応答解析装置として機能することになり、請求項１及び請求項５記載の発明と同様に、減衰定数ｈが振動数非依存特性を示しかつ剛性低下率α及び減衰定数ｈが歪振幅依存特性を示す物体の時刻歴応答解析を容易かつ高精度に行うことができる。

以上説明したように本発明は、減衰定数ｈが振動数非依存特性を示しかつ剛性低下率α及び減衰定数ｈが歪振幅依存特性を示す物体に対し、振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)−(２ω／ωm)なる値(但しｎは整数)を示す虚数部の正則成分、及び、振動数ωに拘わらず(２ω／ωm)なる値を示す虚数部の特異成分の和で表され、振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)なる値を示す虚数部と、虚数部の正則成分のヒルベルト変換値に対応する実数部と、から成る因果的単位虚数関数を時間領域へ変換するか、又は前記虚数部のみを時間領域へ変換することで、因果的単位虚数関数のインパルス応答値として、物体の速度に依存する同時成分ｃ(t₀)、物体の変位に依存する同時成分ｋ(t₀)、物体の変位に依存するΔt刻みの時間遅れ成分ｋ(t_j)（但しｊは自然数でt_j＝Δt・j）を演算し、物体の質量マトリクスを[Ｍs]、物体の初期剛性マトリクスを[Ｋ₀]、時間領域での物体の変位ベクトルを{ｕ(t)}、速度ベクトルを{ｕ'(t)}、反力ベクトルを{Ｆ(γ,t)}、物体を振動させる外力の時間領域での加速度をｙ_０"(t)、時間遅れ成分ｋ(t_j)の総数をｎとしたときに、演算したインパルス応答値を、
[Ｍs]{ｕ"(t)}＋{Ｆ(γ,t)}＝−ｙ_０"(t)[Ｍs]{1} …(１)
但し、

上式に代入し、剛性低下率α(γ)及び減衰定数ｈ(γ)を各時刻における前記物体の歪レベルγに応じた値に切り替えながら物体の反力ベクトル{Ｆ(γ,t)}、変位ベクトル{ｕ(t)}及び速度ベクトル{ｕ'(t)}をΔt刻みで順次演算することで、物体の時刻歴応答解析を行うので、減衰定数ｈが振動数非依存特性を示しかつ剛性低下率α及び減衰定数ｈが歪振幅依存特性を示す物体の時刻歴応答解析を容易かつ高精度に行うことができる、という優れた効果を有する。

以下、図面を参照して本発明の実施形態の一例を詳細に説明する。図４には本発明を適用可能なパーソナル・コンピュータ（ＰＣ）１０が示されている。ＰＣ１０は、ＣＰＵ１０Ａ、ＲＯＭ１０Ｂ、ＲＡＭ１０Ｃ及び入出力ポート１０Ｄが、データバス、制御バス、アドレスバス等から成るバス１０Ｅを介して互いに接続されて構成されている。また入出力ポート１０Ｄには、各種の周辺機器として、ＣＲＴ又はＬＣＤから成るディスプレイ１２、キーボード１４、マウス１６、プリンタ１８、ハードディスクドライブ（ＨＤＤ）２０、ＣＤ−ＲＯＭ２２からの情報の読み出しを行うＣＤ−ＲＯＭドライブ２４が各々接続されている。

ＰＣ１０のＨＤＤ２０には、インパルス応答値データベース（ＤＢ）が記憶されており（詳細は後述）、後述する時刻歴応答解析処理を行うための時刻歴応答解析プログラムもインストールされている。時刻歴応答解析プログラムは請求項６記載の発明に係る時刻歴応答解析プログラムに対応している。ＰＣ１０は、ＣＰＵ１０Ａが時刻歴応答解析プログラムを実行することで、請求項５記載の発明に係る時刻歴応答解析装置として機能する。なお、ＰＣ１０は請求項６に記載のコンピュータにも対応しているが、請求項６に記載のコンピュータはＰＣ１０に限られるものではなく、例えばワークステーションや汎用の大型コンピュータ等であってもよい。

次に本実施形態の作用として、まずＨＤＤ２０に記憶されているインパルス応答値ＤＢについて説明する。このインパルス応答値ＤＢには、因果的単位虚数関数Ｚ'(ω)の解析対象上限振動数ωmと、因果的単位虚数関数Ｚ'(ω)の時間領域への変換に用いるデータ点の数の少なくとも一方が互いに異なる複数種のインパルス応答値が各々記憶されている。インパルス応答値ＤＢは、例として図５に示すインパルス応答値演算処理を任意のコンピュータ（ＰＣ１０でもよいし、ＰＣ１０とは別のコンピュータであってもよい）で実行させることによって生成される。以下、このインパルス応答値演算処理について説明する。なお、以下で説明するインパルス応答値演算処理は、請求項１のうち因果的単位虚数関数のインパルス応答値を演算するステップに対応している。

図５に示すインパルス応答値演算処理では、予め定められた解析対象上限振動数ωmの複数種の値について因果的単位虚数関数Ｚ'(ω)のインパルス応答値を各々演算する。このため、まずステップ１００では、解析対象上限振動数ωmの複数種の値の中からインパルス応答値未演算の任意の解析対象上限振動数ωmを選択する。またステップ１０２では、ステップ１００で選択した解析対象上限振動数ωmに対応する因果的単位虚数関数Ｚ'(ω)の虚数部Ｚ'_Ｉ(ω)を規定する正則成分Ｚ'_ＩＲ(ω) 及び特異成分Ｚ'_ＩＳ(ω)を設定する。なお、図１(Ｂ)に示すように正則成分Ｚ'_ＩＲ(ω)は振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)−(２ω／ωm)なる値(但しｎは整数)を示し、図１(Ｃ)に示すように特異成分Ｚ'_ＩＳ(ω)は振動数ωに拘わらず(２ω／ωm)なる値を示すので、正則成分Ｚ'_ＩＲ(ω)と特異成分Ｚ'_ＩＳ(ω)の和である因果的単位虚数関数Ｚ'(ω)の虚数部Ｚ'_Ｉ(ω) (Ｚ'_Ｉ(ω)＝Ｚ'_ＩＲ(ω)＋Ｚ'_ＩＳ(ω))は、図１(Ａ)に示すように振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)なる値を示す。

次のステップ１０４では、ステップ１０２で設定した正則成分Ｚ'_ＩＲ(ω)に対してヒルベルト変換(前出の(７)式を参照)を行うことで、因果律を満たす因果的単位虚数関数Ｚ'(ω)の実数部Ｚ'_Ｒ(ω)を算出する。これにより、正則成分Ｚ'_ＩＲ(ω)と特異成分Ｚ'_ＩＳ(ω)の和で表される虚数部Ｚ'_Ｉ(ω)と、振動数ωの変化に対して例えば図３に示すように値が変化する実数部Ｚ'_Ｒ(ω)と、から成る因果的単位虚数関数Ｚ'(ω)が得られ、ステップ１０４では得られた因果的単位虚数関数Ｚ'(ω)のデータをメモリ（ＲＡＭ１０Ｃ）に一時記憶させる。なお、図３に示す実数部Ｚ'_Ｒ(ω)は解析対象上限振動数ωm＝20Hzの場合であり、ステップ１０４で算出される実数部Ｚ'_Ｒ(ω)は、実際には振動数ω＝０Hzからステップ１００で選択した解析対象上限振動数ωmの範囲内で図３に示すように値が変化する特性を示す。

また本実施形態に係るインパルス応答値演算処理では、単一の因果的単位虚数関数Ｚ'(ω)を時間領域へ変換してインパルス応答値を演算する際にサンプリングを行うデータ点の数として複数種の値が予め設定されており、ステップ１０６では、予め設定された複数種のデータ点数の中からインパルス応答値未演算のデータ点数を選択する。そしてステップ１０８では、ステップ１０６で選択したデータ点数に対応するデータ点で因果的単位虚数関数Ｚ'(ω)をサンプリングし、時間領域に変換することでインパルス応答値を算出する。このインパルス応答値の算出は、以下のようにして行うことができる。

すなわち、まずステップ１０２，１０４の処理によって得られた因果的単位虚数関数Ｚ'(ω)のデータから、振動数ω＝０〜ωmの範囲内で、ステップ１０６で選択したデータ点数Ｎに対応するＮ種の振動数ω_１,…,ω_Ｎにおける因果的単位虚数関数Ｚ'(ω)の値を表すＮ個の複素データＳ(ω_１),…,Ｓ(ω_Ｎ)を各々抽出（サンプリング）し、抽出した複素データをメモリ又はＨＤＤ２０に記憶させる。なお、因果的単位虚数関数Ｚ'(ω)のデータから抽出した複素データは、次に述べるようにインパルス応答の演算に用いられ、この演算により時刻ｔ＝０及び時刻ｔ＝Δt・j（j=1,2,…)の各時刻におけるインパルス応答値が得られるが、得られるインパルス応答値の個数は演算に用いる複素データの個数に応じて定まり（すなわちjの最大値jmax＝複素データの個数−１(＝ｎ))、得られるインパルス応答値によって表される因果的単位虚数関数Ｚ'(ω)のインパルス応答の時刻範囲も演算に用いる複素データの個数に応じて定まる（例えば複素データの個数が21個、Δt＝0.05秒とすると、ｔmax＝Δt・jmax＝0.05×20＝１秒となり、時刻ｔ＝０〜１秒の時刻範囲の因果的単位虚数関数Ｚ'(ω)のインパルス応答を表す21個のインパルス応答値が得られる）ことになるので、因果的単位虚数関数Ｚ'(ω)のデータから抽出する複素データの個数（複数種のデータ点数の各々の値）は、因果的単位虚数関数Ｚ'(ω)のインパルス応答を算出すべき時刻範囲の長さも勘案して予め定められている。

次に、周波数領域の動的剛性を時間領域のインパルス応答へ変換するために本願発明者が特許第3878626号で提案した次の(15)式及び(16)式の連立方程式をＨＤＤ２０から読み出し、読み出した連立方程式に、先に因果的単位虚数関数Ｚ'(ω)のデータから抽出したＮ個の複素データＳ(ω_１),…,Ｓ(ω_Ｎ)を代入し、この連立方程式の解を求めることで、因果的単位虚数関数Ｚ'(ω)のインパルス応答を表すインパルス応答値を、予め設定されたΔｔ刻みで演算する。

上記の演算により、因果的単位虚数関数Ｚ'(ω)のインパルス応答値として、変位に依存するインパルス応答の剛性項の同時成分ｋ(t_０)、速度に依存するインパルス応答の減衰項の同時成分ｃ(t_０)、加速度に依存するインパルス応答の質量項の同時成分ｍ(t_０)のデータが得られると共に、変位に依存するインパルス応答の剛性項の時間遅れ成分ｋ(t_ｊ)のデータがΔt刻みでｎ個（ｎ＝Ｎ−１）得られ、速度に依存するインパルス応答の減衰項の時間遅れ成分ｃ(t_ｊ)のデータがΔt刻みでｎ−１個得られることになる。但し、因果的単位虚数関数Ｚ'(ω)は、振動数範囲ω＝０〜ωmの中央に相当する振動数(ωm／２)に関して実数部Ｚ'_Ｒ(ω)が線対称、虚数部Ｚ'_Ｉ(ω)が点対称の変化を示すため、上記各データのうち、質量項の同時成分ｍ(t_０)及び減衰項の時間遅れ成分ｃ(t_１)〜ｃ(t_ｎ−１)の値は何れも０となる。このため、因果的単位虚数関数Ｚ'(ω)のインパルス応答値として、
剛性項の同時成分ｋ(t_０)、減衰項の同時成分ｃ(t_０)、及び、剛性項の時間遅れ成分ｋ(t_１)〜ｋ(t_ｎ)がメモリ又はＨＤＤ２０に記憶される。

上記のようにして、選択した解析対象上限振動数ωm及びデータ点数に対応する因果的単位虚数関数Ｚ'(ω)のインパルス応答値が得られると、次のステップ１１０では、得られたインパルス応答値を、選択した解析対象上限振動数ωm及びデータ点数を表す情報と対応付けてインパルス応答値ＤＢへ登録する。次のステップ１１２では、予め定められた複数種のデータ点数の中にインパルス応答値を未演算のデータ点数が有るか否か判定する。判定が肯定された場合はステップ１０６に戻り、ステップ１１２の判定が否定される迄ステップ１０６〜ステップ１１２を繰り返す。これにより、単一の因果的単位虚数関数Ｚ'(ω)から、サンプリングを行うデータ点の数が互いに異なる複数種のインパルス応答値が各々演算されてインパルス応答値ＤＢに各々登録されることになる。

また、単一の因果的単位虚数関数Ｚ'(ω)について、データ点数が互いに異なる複数種のインパルス応答値が各々演算されてインパルス応答値ＤＢに各々登録されることでステップ１１２の判定が否定されると、次のステップ１１４において、予め定められた複数種の解析対象上限振動数ωmの中にインパルス応答値が未演算の解析対象上限振動数ωmが有るか否か判定する。判定が肯定された場合はステップ１００に戻り、ステップ１１４の判定が否定される迄ステップ１００〜ステップ１１４を繰り返す。これにより、インパルス応答値ＤＢには解析対象上限振動数ωm又はデータ点数の少なくとも一方が互いに異なる因果的単位虚数関数Ｚ'(ω)の多数のインパルス応答値が各々登録されることになる。そして、ステップ１１４の判定が否定されるとインパルス応答値演算処理を終了する。

なお、上述したインパルス応答値演算処理をＰＣ１０とは別のコンピュータで実行させた場合、生成されたインパルス応答値ＤＢは、例えばＰＣ１０への時刻歴応答解析プログラムのインストール時や他のタイミング（ＰＣ１０で時刻歴応答解析処理が実行される前のタイミング）でＰＣ１０のＨＤＤ２０に複写される。このように、ＨＤＤ２０は請求項５，６に記載の記憶手段に対応している。

続いて、時刻歴応答解析の実行を所望しているオペレータによりキーボード１４又はマウス１６を介して時刻歴応答解析プログラムの実行が指示されることで、ＰＣ１０のＣＰＵ１０Ａで実行される時刻歴応答解析処理について、図６のフローチャートを参照して説明する。この時刻歴応答解析処理は、減衰定数ｈが振動数非依存特性を示しかつ剛性低下率α及び減衰定数ｈが歪振幅依存特性を示す物体を解析対象とし、上述したインパルス応答値ＤＢに記憶されているインパルス応答値を用いて時刻歴応答解析を行う処理であり、請求項１のうち物体の時刻歴応答解析を行うステップに対応していると共に、請求項５，６に記載の解析手段に対応している。なお、以下では解析対象の物体として地盤を適用した例を説明するが、本発明はこれに限定されるものではなく、地盤上に建設され地盤の挙動の影響を受ける構造物も解析対象とし、その挙動も同時に解析するようにしてもよいし、解析対象の物体としてその他の物体（例えば粘弾性体等）を適用してもよい。

本実施形態では、解析対象の物体に関する各種データ（例えば解析対象の物体の質量マトリクス[Ｍs]や剛性マトリクス[Ｋs]、初期剛性マトリクス[Ｋ₀]、剛性低下率特性γ−α、減衰定数特性γ−Ｈ等）がＰＣ１０に予め入力されてＨＤＤ２０に記憶されており、時刻歴応答解析処理では、まずステップ１２０において、解析対象の物体に関する各種データをＨＤＤ２０から取得する。なお、解析対象の物体が地盤である場合、当該地盤の質量マトリクス[Ｍs]や剛性マトリクス[Ｋs]、初期剛性マトリクス[Ｋ₀]は、例えば図７(Ａ)に示すように、解析対象の地盤の地層構成、各地層の層厚Ｈ、地震波の伝播速度Ｖs、ポアソン比ν、密度ρ、減衰率ｈ等のデータを用いて所定の演算を行うことで得ることができる。また、地盤の剛性低下率特性γ−α及び減衰定数特性γ−Ｈ（一例を図７(Ｂ)に示す）については、実験等を行って求めてもよいし、文献等に記載されている特性を用いてもよい。

次のステップ１２２では、ステップ１２０で取得した解析対象の物体に関する各種データのうち質量マトリクス[Ｍs]及び剛性マトリクス[Ｋs]に基づき、解析対象の物体に対して実固有値解析を行い、解析対象の物体の主要な固有振動数（例えば一次振動数及び二次振動数）と固有モードを算出する。物体に入力する地震動が比較的小さい場合には、物体に損傷劣化が生じないので物体の応答は線形範囲内に収まり、物体の剛性や固有振動数の変化は生じないが、入力する地震動が比較的大きくなると、物体に損傷劣化が生ずることで物体の応答が非線形化し、物体の剛性及び固有振動数が低下する。そして、物体の応答の非線形化に伴って変化した固有振動数が、後述する時刻歴応答解析における解析対象振動数範囲０〜ωmから外れたり、解析対象振動数範囲０〜ωmのうち精度の低い振動数範囲に入ってしまうと、時刻歴応答解析の途中で処理を継続できない状態に陥ったり、解析精度が低下する恐れがある。

このため、次のステップ１２４では、解析対象の物体に対して実行を予定している時刻歴応答解析で入力する地震動と同一の地震動を入力したときの解析対象の物体の概略の挙動を解析する予備解析を行い、前記地震動の入力によって解析対象の物体の応答が非線形化するか否かを判断すると共に、解析対象の物体の応答の非線形化が生ずる場合には非線形化後の固有振動数を算出する。なお、上記の予備解析としては、例えば剛性比例型の減衰モデルやRayleigh型の減衰モデルを用いた簡易的な時刻歴応答解析を行うことで実現できる。また、上記の予備解析を省略し、物体の損傷（応答の非線形化）が生ずるか否か、及び、物体の応答の非線形化が生ずる場合の非線形化後の固有振動数を、オペレータが経験的に判断して判断結果を入力するようにすることも可能である。

次のステップ１２６では、ステップ１２２で算出した解析対象の物体の主要な固有振動数に基づいて、時刻歴応答解析に用いるインパルス応答値の解析対象上限振動数ωm及びデータ点数を選択する。また、ステップ１２４の予備解析によって解析対象の物体の損傷（応答の非線形化）が生ずると判断された場合、ステップ１２６では予備解析で算出した非線形化後の固有振動数も考慮して解析対象上限振動数ωm及びデータ点数を選択する。より具体的には、例えば「解析対象の物体の主要な固有振動数（及び非線形化後の固有振動数）が、振動数ω＝０〜ωmの範囲内のうち剛性及び減衰定数の誤差が所定値未満（例えば±10％未満）となる振動数範囲内に全て入る」という条件を満たすように解析対象上限振動数ωmを選択する。

データ点数に関しては、或る程度以上のデータ点数があればデータ点数を更に増加させても精度はそれ程変化しないことが本願発明者によって確認されているので（特開２００７−０７１７５４号公報参照）、時刻歴応答解析処理の簡素化を考慮し、上記の条件を満たす範囲内でなるべく少ないデータ点数を選択すればよい。また、データ点数が多くなるに従い剛性及び減衰定数の誤差が所定値未満となる振動数範囲が若干広くなることも本願発明者によって確認されているので（特開２００７−０７１７５４号公報参照）、一部の固有振動数が、剛性及び減衰定数の誤差が所定値未満となる振動数範囲の端部付近に存在している場合には、時刻歴応答解析処理の精度向上のためにデータ点数を増加させるようにしてもよい。解析対象上限振動数ωm及びデータ点数を選択すると、次のステップ１２８では、選択した解析対象上限振動数ωm及びデータ点数に対応するインパルス応答値(剛性項の同時成分ｋ(t_０)、減衰項の同時成分ｃ(t_０)、及び、剛性項の時間遅れ成分ｋ(t_１)〜ｋ(t_ｎ))をインパルス応答値ＤＢから読み込む。

次のステップ１３０以降では、外力が加わった場合の解析対象の物体の挙動を解析する時刻歴応答解析を行う。すなわち、ステップ１３０では解析対象時刻ｔを０に初期化する。次のステップ１３２では、解析対象時刻ｔに解析対象の物体に加わる外力を表すデータを取り込み、取り込んだデータが表す外力に基づいて、解析対象時刻ｔにおける解析対象の物体の変位ベクトル{ｕ(t)}、速度ベクトル{ｕ'(t)}、加速度ベクトル{ｕ"(t)}を各々推定する。なお、解析対象の物体が地盤であり、外力が地震動である場合、前記外力を表すデータの取り込みは、予め想定した地震動のデータから解析対象時刻ｔに解析対象の地盤に加わる地震動のデータを抽出することで行うことができる。またステップ１３４では、ステップ１３２で推定した解析対象時刻ｔにおける解析対象の物体の変位ベクトル{ｕ(t)}に基づいて、解析対象時刻ｔにおける解析対象の物体の歪レベルγを演算し、演算結果をＲＡＭ１０Ｃに記憶する。

次のステップ１３６では、先のステップ１２０で取得した解析対象の物体に関する各種のデータのうち、解析対象の物体の剛性低下率特性γ−α及び減衰定数特性γ−ｈを表すデータに基づき、解析対象時刻ｔにおける応答解析に用いる解析対象の物体の剛性低下率α及び減衰定数ｈとして、ステップ１３４で演算した歪レベルγに対応する剛性低下率α及び減衰定数ｈを各々演算する。なお、剛性低下率α及び減衰定数ｈの演算は、剛性低下率特性γ−αを表すデータや減衰定数特性γ−ｈを表すデータが、歪レベルγと剛性低下率α（又は減衰定数ｈ）との関係を関数等の形態で表すデータであれば、当該関数等にステップ１３４で演算した歪レベルγを代入して演算することで実現できるが、剛性低下率特性γ−αを表すデータや減衰定数特性γ−ｈを表すデータが、歪レベルγと剛性低下率α（又は減衰定数ｈ）との関係を離散的に表すデータ（歪レベルγが互いに異なる複数種の値のときの剛性低下率α（又は減衰定数ｈ）の値を各々表すデータ等）である場合は、ステップ１３４で演算した歪レベルγの前後の歪レベルγに対応する剛性低下率α（又は減衰定数ｈ）の複数の値を抽出し、抽出した複数の値に基づいて、ステップ１３４で演算した歪レベルγに対応する剛性低下率α（又は減衰定数ｈ）を補間演算によって求めることで実現できる。

ここで、解析対象時刻ｔにおける応答解析に用いる解析対象の物体の剛性低下率α及び減衰定数ｈは、前出の(２)式に代入する剛性低下率α及び減衰定数ｈを意味しているが、(２)式の右辺第１項は解析対象時刻ｔに対応する同時成分の項であるので、右辺第１項における剛性低下率α(γ)及び減衰定数ｈ(γ)としては、解析対象時刻ｔにおける歪レベルγに対応する剛性低下率α及び減衰定数ｈをそのまま用いることができる。一方、(２)式の右辺第２項は解析対象時刻ｔよりも前の時刻ｔ−ｔj（但し、ｊは自然数でｔj＝Δｔ・ｊ)における解析対象の物体の挙動が及ぼす影響を表す時間遅れ成分の項であるのに対し、外力が入力されると解析対象の物体の歪レベルγも時々刻々変化するので、(２)式の右辺第２項における剛性低下率α(γ)及び減衰定数ｈ(γ)として、何れの時点での歪レベルγに対応する値を用いるのかという点に関しては選択の余地がある。

(２)式の右辺第２項の剛性低下率α(γ)及び減衰定数ｈ(γ)として用いる剛性低下率α及び減衰定数ｈを決定するための歪レベルγの選択肢としては、例えば解析対象時刻ｔよりも前の時刻ｔ−ｔjにおける解析対象の物体の歪レベルγ(ｔ−ｔj)を用いるケース（ケース１と称する）、解析対象時刻ｔにおける解析対象の物体の歪レベルγ(ｔ)を用いるケース（ケース３と称する）、ケース１とケース３の平均値に相当する歪レベル(γ(ｔ−ｔj)+γ(ｔ)／２）を用いるケース（ケース２と称する）が挙げられる。

本願発明者が実施した解析検討の結果、(２)式の右辺第２項の剛性低下率α(γ)及び減衰定数ｈ(γ)として、ケース１〜ケース３の何れの歪レベルγに対応する剛性低下率α及び減衰定数ｈを用いた場合にも時刻歴応答解析の精度には大きな差が無いことが確認されており、(２)式の右辺第２項の剛性低下率α(γ)及び減衰定数ｈ(γ)としては、ケース１〜ケース３の何れの歪レベルγに対応する剛性低下率α及び減衰定数ｈを用いてもよい。また、解析対象の物体の種類や条件等によっては、(２)式の右辺第２項の剛性低下率α(γ)及び減衰定数ｈ(γ)として、ケース１〜ケース３の何れの歪レベルγに対応する剛性低下率α及び減衰定数ｈを用いるかに応じて時刻歴応答解析の精度が変動する可能性もあるが、ケース１〜ケース３の何れの歪レベルγに対応する剛性低下率α及び減衰定数ｈを用いたとしても演算量自体は略同じであるので、その場合はケース１〜ケース３のうち最も高い精度が得られる何れかのケースの歪レベルγに対応する剛性低下率α及び減衰定数ｈを選択的に用いればよい。

次のステップ１３８では、ステップ１３２で仮定した解析対象の物体の変位ベクトル{ｕ(t)}、速度ベクトル{ｕ'(t)}、加速度ベクトル{ｕ"(t)}、ステップ１２０で取得した解析対象の物体の質量マトリクス[Ｍs]及び初期剛性マトリクス[Ｋ₀]、ステップ１２８で取り込んだインパルス応答値、ステップ１３６で演算した剛性低下率α及び減衰定数ｈを前出の(１),(２)式に代入し、解析対象時刻ｔにおける解析対象の物体の挙動を演算する。またステップ１４０では、ステップ１３８の演算の結果、解析対象時刻ｔに解析対象の物体に加わる外力が、解析対象時刻ｔにおける解析対象の物体の反力と釣り合っているか否か判定する。(１)式の運動方程式における右辺は解析対象時刻ｔに解析対象の物体に加わる外力を表しており、(１)式の左辺は解析対象時刻ｔにおける解析対象の物体の反力を表している。ステップ１４０の判定は、(１)式の右辺の値と(６)式の左辺の値との偏差（外力と反力との釣合の誤差）が許容範囲内か否かを判断することで行われる。

ステップ１４０の判定が否定された場合はステップ１３２へ戻り、ステップ１３２において、先に仮定した解析対象時刻ｔにおける解析対象の物体の変位ベクトル{ｕ(t)}、速度ベクトル{ｕ'(t)}、加速度ベクトル{ｕ"(t)}を修正した後にステップ１３４以降の処理を行う。上述したステップ１３２〜１４０は、ステップ１４０の判定が肯定される迄繰り返されるので、解析対象時刻ｔにおける解析対象の物体の変位ベクトル{ｕ(t)}、速度ベクトル{ｕ'(t)}、加速度ベクトル{ｕ"(t)}は、解析対象時刻ｔにおける外力に対する解析対象の物体の反力の偏差が許容範囲内となる値に収束することになる。

ステップ１４０の判定が肯定され、解析対象時刻ｔにおける外力に対する反力の偏差が許容範囲内となる解析対象の物体の変位ベクトル{ｕ(t)}、速度ベクトル{ｕ'(t)}、加速度ベクトル{ｕ"(t)}が求まるとステップ１４２へ移行し、上記の時刻歴応答解析（各時刻における解析対象の物体の変位ベクトル{ｕ(t)}、速度ベクトル{ｕ'(t)}、加速度ベクトル{ｕ"(t)}の演算）を、時刻歴応答解析の解析終了時刻ｔmax迄行ったか（解析対象時刻ｔが解析終了時刻ｔmaxに達したか）否か判定する。判定が否定された場合はステップ１４４へ移行し、解析対象時刻ｔにΔｔを加えることで解析対象時刻ｔを更新してステップ１３２に戻る。

これにより、ステップ１４２の判定が肯定される迄ステップ１３２〜１４４が繰り返され、解析対象時刻ｔ＝０から時間Δｔ刻みの各時刻について、演算に用いる剛性低下率α及び減衰定数ｈを各時刻における解析対象の物体の歪レベルγに応じた値に切り替えながら時刻歴応答解析（解析対象の各時刻における解析対象の物体の変位ベクトル{ｕ(t)}、速度ベクトル{ｕ'(t)}、加速度ベクトル{ｕ"(t)}の演算）が順次行われることになる。そして、ステップ１４２の判定が肯定されると時刻歴応答解析処理を終了する。

上述したように、本実施形態に係る時刻歴応答解析処理は、減衰定数ｈが振動数非依存特性を示す物体をモデル化して表す因果的単位虚数関数のインパルス応答値を用いると共に、物体の歪レベルγの変化に伴って変化する物体の剛性低下率α及び減衰定数ｈを明示的に表して物体の挙動をモデル化した(１),(２)式（非線形の因果的履歴減衰モデル）を用い、解析対象時刻ｔが各時刻のときの解析対象の物体の歪レベルγを演算し、演算に用いる剛性低下率α及び減衰定数ｈを各時刻における解析対象の物体の歪レベルγに応じた値に切り替えながら、各時刻における解析対象の物体の変位、速度及び加速度を演算する時刻歴応答解析を行うので、解析対象の物体の歪レベルγの変化に伴って解析対象の物体の剛性低下率α及び減衰定数ｈが変化し、この剛性低下率α及び減衰定数ｈの変化に伴って解析対象の物体の挙動が変化することを、上記の剛性低下率α及び減衰定数ｈの切替えによって各時刻の解析に反映させることができる。従って、解析対象の物体の時刻歴応答解析を、解析対象の物体の減衰定数ｈが振動数非依存特性を示しかつ解析対象の物体の剛性低下率α及び減衰定数ｈが歪振幅依存特性を示すことによる影響を考慮して高精度に行うことができる。

なお、上記では因果的単位虚数関数Ｚ'(ω)の解析対象上限振動数ωmと、因果的単位虚数関数Ｚ'(ω)の時間領域への変換に用いるデータ点の数の少なくとも一方が互いに異なる複数種のインパルス応答値をインパルス応答値ＤＢに予め記憶しておく態様を説明したが、本発明はこれに限定されるものではなく、時刻歴応答解析処理を行う都度、選択した解析対象上限振動数ωmに対応する因果的単位虚数関数Ｚ'(ω)のデータから、選択したデータ点数の複素データを抽出し、時間領域への変換を行ってインパルス応答値を演算するようにしてもよい。

また、上記では因果的単位虚数関数Ｚ'(ω)の虚数部Ｚ'_Ｉ(ω)の正則成分Ｚ'_ＩＲ(ω)に対してヒルベルト変換を行うことで、因果律を満たす因果的単位虚数関数Ｚ'(ω)の実数部Ｚ'_Ｒ(ω)を算出した後に、因果的単位虚数関数Ｚ'(ω)を時間領域へ変換することでインパルス応答値を演算する態様を説明したが、本発明はこれに限定されるものではなく、虚部変換法を適用し、因果的単位虚数関数Ｚ'(ω)の虚数部Ｚ'_Ｉ(ω)のみを時間領域へ変換することで、因果的単位虚数関数Ｚ'(ω)の実数部Ｚ'_Ｒ(ω)を算出することなく、インパルス応答値を演算するようにしてもよい。

また、上記では本発明に係る時刻歴応答解析を、地震動が入力された際の物体の応答を解析する地震応答解析に適用した例を説明したが、本発明はこれに限定されるものではなく、減衰定数ｈが振動数非依存特性を示しかつ剛性低下率α及び減衰定数ｈが歪振幅依存特性を示す任意の物体について、地震動やそれ以外の任意の外力が入力された際の時刻歴応答解析に適用可能である。

更に、上記では本発明に係る時刻歴応答解析プログラムがＰＣ１０のＨＤＤ２０に予め記憶（インストール）されている態様を説明したが、本発明に係る時刻歴応答解析プログラムは、ＣＤ−ＲＯＭやＤＶＤ−ＲＯＭ等の記録媒体に記録されている形態で提供することも可能である。

次に、本発明（に係る非線形の因果的履歴減衰モデル）の有用性を確認するために本願発明者が実施した解析検討の結果について説明する。

（解析検討の条件）
この解析検討では、解析対象の地盤として、図８(Ａ)に示すように、層厚Ｈ＝20ｍ、地震波の伝播速度Ｖs＝200ｍ／ｓ(但し弾性時)、密度ρ＝2.0ｔ／ｍ^３の地層が剛基盤上に存在している地層構成の地盤モデル（図８(Ｂ)に示す１自由度モデル）を用い、この地盤モデルに対し、Ｒ−Ｏモデル、双曲線モデル及び本発明に係る非線形の因果的履歴減衰モデルを各々適用して地震応答解析（時刻歴応答解析）を各々行った。

解析検討に用いたＲ−Ｏモデルは、剛性低下率α＝0.5のときの歪レベルγ＝0.1%、減衰定数ｈの最大値ｈmax＝26％に設定した。また、解析検討に用いた双曲線モデルは、剛性低下率α＝0.5となるときの歪レベルγ＝0.1%とした。これは、最大剪断応力τmaxを0.001Ｇ_０に設定したことに相当する。上記パラメータを設定したＲ−Ｏモデル及び双曲線モデルによって表される歪レベルγ−剛性低下率α特性、歪レベルγ−減衰定数ｈ特性を図９に示す。また、入力地震動としては、神戸位相のレベル２告示波（時間刻み0.005秒、継続時間30秒）を係数倍（0.2倍, 0.5倍, 1.0倍, 2.0倍, 4.0倍）して用いた。

また、本発明に係る非線形の因果的履歴減衰モデルについては、解析対象の振動数範囲を原則として０〜10Ｈzとし、(２)式の右辺第２項に代入する剛性低下率α及び減衰定数ｈを決定するための歪レベルγとしてはケース２を適用した（後述する解析検討の結果はケース２に対応しているが、本願発明者はケース１，３についても解析検討を行い、ケース２と解析結果の差異が小さいことを確認している）。また歪レベルγに応じて変化する各解析対象時刻ｔにおける剛性低下率α及び減衰定数ｈの値は、各解析対象時刻ｔから過去１秒間の歪レベルγの最大値に応じて変化するものとした。本願発明者は、この時間（過去の歪レベルγの最大値を記憶している時間）を0.25秒〜2.0秒の間で変化させた解析検討も行ったが、この時間を変化させても解析結果に大きな差異が生じないことを確認している。また、時間積分にはNewmark-β法（β=1/4）、収束計算には修正Newton-Raphson法を適用した。

（Ｒ−Ｏモデルとの比較検討）
Ｒ−Ｏモデルと比較検討を行うため、本発明に係る非線形の因果的履歴減衰モデルに対し、図９に示すＲ−Ｏモデルと同一の歪レベルγ−剛性低下率α特性、歪レベルγ−減衰定数ｈ特性を与え、前述の地盤モデルに対して地震応答解析（時刻歴応答解析）を行い、全ケースの全ての解析対象時刻で演算結果の収束を得た。次の表１に、非線形の因果的履歴減衰モデルを用いた地震応答解析によって得られた最大応答値と、Ｒ−Ｏモデルを用いた地震応答解析によって得られた最大応答値との比較を示す。

表１からも明らかなように、Ｒ−Ｏモデルによって得られた最大応答値に比して、非線形の因果的履歴減衰モデルによって得られた最大応答値は全体的に若干小さめの値となっているが、その差異はほぼ±10％以内である。

また図１０には、入力地震動のレベルが0.2倍、１倍、４倍のときに非線形の因果的履歴減衰モデル及びＲ−Ｏモデルによって得られた加速度波形（5〜20秒）を各々示す。図１０から明らかなように、両モデルによって得られた加速度波形は、若干の相違はあるものの全体としてはおおよそ一致している。

更に図１１には入力地震動のレベルが0.2倍、１倍、４倍のときに非線形の因果的履歴減衰モデル及びＲ−Ｏモデルによって得られた歪と応力の関係を各々示す。図１１から明らかなように、非線形の因果的履歴減衰モデルによって得られる歪と応力の関係においても、Ｒ−Ｏモデルと同様、歪の増大に伴い、ループの膨らみの増大（減衰定数ｈの増大を表す）や、正負の最大値を結ぶ中心軸の傾きの減少（剛性の低下を表す）が生じている。両モデルによる結果は、ループの形状がやや相違しているものの、ループの大きさや傾きは全体としておおよそ一致している。

以上の結果より、比較検討を行った非線形の因果的履歴減衰モデルは、入力地震動のレベルの変化に拘わらず、Ｒ−Ｏモデルとおよそ同等の精度の良好な解析結果が得られることが理解できる。

（双曲線モデルとの比較検討）
続いて、双曲線モデルと比較検討を行うため、本発明に係る非線形の因果的履歴減衰モデルに対し、図９に示す双曲線モデルと同一の歪レベルγ−剛性低下率α特性、歪レベルγ−減衰定数ｈ特性を与え、前述の地盤モデルに対して地震応答解析（時刻歴応答解析）を行った。次の表２に、非線形の因果的履歴減衰モデルを用いた地震応答解析によって得られた最大応答値と、双曲線モデルを用いた地震応答解析によって得られた最大応答値との比較を示す。

非線形の因果的履歴減衰モデルを用いた地震応答解析では、入力地震動のレベルが0.2倍〜１倍の範囲では、全ての解析対象時刻で演算結果が収束し、表２からも明らかなように、非線形の因果的履歴減衰モデルによって得られた最大応答値は、Ｒ−Ｏモデルによって得られた最大応答値に対して±15％以内に収まった。

一方、非線形の因果的履歴減衰モデルを用いた地震応答解析において、入力地震動のレベルが2.0倍以上のケースでは演算結果（解）が発散した。これは以下の原因が考えられる。すなわち、本解析検討で用いた双曲線モデルの動的変形特性（歪レベルγ−剛性低下率α特性、歪レベルγ−減衰定数ｈ特性）は、歪レベルγが0.1％以上となると減衰定数ｈが急速に上昇し、剛性も急速に低下する。入力地震動のレベルが2.0倍のケースでは、歪レベルγが最大約1.5％で、その時点の減衰定数ｈが約50％、剛性低下率αが約10％となっている。そして、歪レベルγが最大となっている時の共振振動数は初期の2.5Ｈzから0.8Ｈzに低下するものと考えられる。本発明に係る非線形の因果的履歴減衰モデルは、複素減衰Ｋ₀(１＋２ｈ・ｉ)の虚数単位ｉを時間領域で近似するものである。このため減衰定数ｈの値が小さくなるほど精度が高く、逆に減衰定数ｈの値が大きくなるほど精度が低下する傾向がある。この精度の低下は解析対象の振動数範囲の下限や上限及びその付近で大きい。本願発明者が詳細に検討したところ、減衰定数ｈが40％程度以上の場合に解の発散が生じていることが判明した。解析検討の条件にもよるが、減衰定数ｈがこのように大きな値となる場合には、モデルの設定に関して検討の余地があると考えられるが、一方で通常の地盤では減衰定数ｈがこのように大きな値となる場合はかなり限定されるものと考えられる。

また本願発明者は、入力地震動のレベルが2.0倍のケースについて、解析対象の振動数範囲の上限を10Ｈzから５Ｈzに低下させて再度解析検討を行った。この場合は解の収束を得たが、表２からも明らかなように最大応答値の精度は大きく低下した。本発明に係る非線形の因果的履歴減衰モデルでは、解析対象の振動数範囲外に対し、解析対象の振動数範囲内の３倍以上の大きな減衰定数が与えられるため、実質的に5〜10Ｈzの成分がカットされ、これが精度の低下に繋がったものと考えられる。

一方、図１２には、入力地震動のレベルが0.2倍、１倍のときに非線形の因果的履歴減衰モデル及び双曲線モデルによって得られた加速度波形（5〜20秒）を各々示す。図１２から明らかなように、両モデルによって得られた加速度波形は、若干の相違はあるものの全体としてはおおよそ一致している。

更に図１３には入力地震動のレベルが0.2倍、１倍のときに非線形の因果的履歴減衰モデル及び双曲線モデルによって得られた歪と応力の関係を各々示す。図１３から明らかなように、入力地震動のレベルが１倍のケースでループの形状がやや相違しているものの、ループの大きさや傾きは全体としておおよそ一致している。

以上の結果より、比較検討を行った非線形の因果的履歴減衰モデルは、入力地震動のレベルの変化に拘わらず、双曲線モデルとおよそ同等の精度の良好な解析結果が得られることが理解できる。

（まとめ）
上記の解析検討から、本発明に係る非線形の因果的履歴減衰モデルの特性に関し、以下の事項が明らかとなった。すなわち、Ｒ−Ｏモデルや双曲線モデルと同一の動的変形特性（歪レベルγ−剛性低下率α特性、歪レベルγ−減衰定数ｈ特性）を与えた場合、本発明に係る非線形の因果的履歴減衰モデルは、Ｒ−Ｏモデルや双曲線モデルとおよそ同等の精度の良好な解析結果が得られる。Ｒ−Ｏモデルや双曲線モデルは動的変形特性の設定に制限があり、多くの場合、実験等より得られた動的変形特性をそのまま用いることは困難である。これに対し、本発明に係る因果的履歴減衰モデルはそのような制限がなく、任意の動的変形特性で解析を行うことができるので、この点で優れている。

因果的単位虚数関数の、(Ａ)は虚数部全体、(Ｂ)は正則成分、(Ｃ)は特異成分を各々示す線図である。 虚数部の正則成分と特異成分の和で表される因果的単位虚数関数の虚数部全体の特性を示す線図である。 因果的単位虚数関数Ｚ'(ω)の一例を示す線図である。 本実施形態に係るＰＣの概略構成を示すブロック図である。 インパルス応答値演算処理の内容を示すフローチャートである。 時刻歴応答解析処理の内容を示すフローチャートである。 (Ａ)は解析対象の物体としての地盤の構成の一例を示す概念図、(Ｂ)は歪レベルγ−剛性低下率α特性、歪レベルγ−減衰定数ｈ特性の一例を示す線図である。 本願発明者が実施した解析検討に用いた地盤モデルを示す概念図である。 解析検討に適用したＲ−Ｏモデル及び双曲線モデルにおけるγ−α特性、γ−ｈ特性を示す線図である。 本願発明者が行った解析検討で得られた非線形の因果的履歴減衰モデル及びＲ−Ｏモデルの加速度波形を各々示す線図である。 上記の解析検討で得られた非線形の因果的履歴減衰モデル及びＲ−Ｏモデルの歪−応力の関係を各々示す線図である。 上記の解析検討で得られた非線形の因果的履歴減衰モデル及び双曲線モデルの加速度波形を各々示す線図である。 上記の解析検討で得られた非線形の因果的履歴減衰モデル及び双曲線モデルの歪−応力の関係を各々示す線図である。

符号の説明

１０ ＰＣ
１２ ディスプレイ
１４ キーボード
１６ マウス
２０ ＨＤＤ

Claims (6)

1. 減衰定数ｈが物体を振動させる外力の振動数に依存しない振動数非依存特性を示しかつ剛性低下率α及び減衰定数ｈが前記物体の歪レベルγに応じて変化する歪振幅依存特性を示す前記物体の時刻歴応答解析を行う時刻歴応答解析方法であって、
   振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)−(２ω／ωm)なる値(但しｎは整数)を示す虚数部の正則成分、及び、振動数ωに拘わらず(２ω／ωm)なる値を示す虚数部の特異成分の和で表され、振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)なる値を示す虚数部と、前記虚数部の正則成分のヒルベルト変換値に対応する実数部と、から成る因果的単位虚数関数を時間領域へ変換するか、又は前記虚数部のみを時間領域へ変換することで、前記因果的単位虚数関数のインパルス応答値として、前記物体の速度に依存する同時成分ｃ(t₀)、前記物体の変位に依存する同時成分ｋ(t₀)、前記物体の変位に依存するΔt刻みの時間遅れ成分ｋ(t_j)（但しｊは自然数でt_j＝Δt・j）を演算し、
   前記物体の質量マトリクスを[Ｍs]、前記物体の初期剛性マトリクスを[Ｋ₀]、時間領域での物体の変位ベクトルを{ｕ(t)}、速度ベクトルを{ｕ'(t)}、反力ベクトルを{Ｆ(γ,t)}、前記物体を振動させる外力の時間領域での加速度をｙ_０"(t)、時間遅れ成分ｋ(t_j)の総数をｎとしたときに、前記演算したインパルス応答値を、
   [Ｍs]{ｕ"(t)}＋{Ｆ(γ,t)}＝−ｙ_０"(t)[Ｍs]{1} …(１)
   但し、
   
   上式に代入し、剛性低下率α(γ)及び減衰定数ｈ(γ)を各時刻における前記物体の歪レベルγに応じた値に切り替えながら物体の反力ベクトル{Ｆ(γ,t)}、変位ベクトル{ｕ(t)}及び速度ベクトル{ｕ'(t)}をΔt刻みで順次演算することで、前記物体の時刻歴応答解析を行うことを特徴とする時刻歴応答解析方法。
2. 物体の主要な固有振動数を把握する実固有値解析を行い、前記実固有値解析によって把握した前記物体の主要な固有振動数が、振動数ω＝０〜ωmの範囲内のうち誤差が所定値未満となる振動数範囲内に入るように解析対象上限振動数ωmを設定又は選択することを特徴とする請求項１記載の時刻歴応答解析方法。
3. 物体を振動させる外力と前記物体の挙動との関係の非線形化によって前記物体の固有振動数が変化するか否かを推定し、前記物体の固有振動数が変化すると判断した場合には、概略の減衰に基づく予備解析により前記非線形化後の固有振動数をおおよそ把握し、把握した前記非線形化後の固有振動数も振動数ω＝０〜ωmの範囲内のうち誤差が所定値未満となる振動数範囲内に入るように解析対象上限振動数ωmを設定又は選択することを特徴とする請求項２記載の時刻歴応答解析方法。
4. 振動数ωmの値が互いに異なる複数種の因果的単位虚数関数について、時間領域への変換を各々行い前記複数種の因果的単位虚数関数のインパルス応答値を各々演算して記憶手段に記憶しておき、
   前記物体の時刻歴応答解析に際し、前記記憶手段に記憶した前記複数種の因果的単位虚数関数のインパルス応答値のうち、前記把握した固有振動数が振動数ω＝０〜ωmの範囲内のうち誤差が所定値未満となる振動数範囲内に入る特定の因果的単位虚数関数のインパルス応答値を読み出して用いることを特徴とする請求項２又は請求項３記載の時刻歴応答解析方法。
5. 減衰定数ｈが物体を振動させる外力の振動数に依存しない振動数非依存特性を示しかつ剛性低下率α及び減衰定数ｈが前記物体の歪レベルγに応じて変化する歪振幅依存特性を示す前記物体の時刻歴応答解析を行う時刻歴応答解析装置であって、
   振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)−(２ω／ωm)なる値(但しｎは整数)を示す虚数部の正則成分、及び、振動数ωに拘わらず(２ω／ωm)なる値を示す虚数部の特異成分の和で表され、振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)なる値を示す虚数部と、前記虚数部の正則成分のヒルベルト変換値に対応する実数部と、から成る因果的単位虚数関数を時間領域へ変換するか、又は前記虚数部のみを時間領域へ変換することで、前記因果的単位虚数関数のインパルス応答値として演算された、前記物体の速度に依存する同時成分ｃ(t₀)、前記物体の変位に依存する同時成分ｋ(t₀)、前記物体の変位に依存するΔt刻みの時間遅れ成分ｋ(t_j)（但しｊは自然数でt_j＝Δt・j）を記憶する記憶手段と、
   前記物体の質量マトリクスを[Ｍs]、前記物体の初期剛性マトリクスを[Ｋ₀]、時間領域での物体の変位ベクトルを{ｕ(t)}、速度ベクトルを{ｕ'(t)}、反力ベクトルを{Ｆ(γ,t)}、前記物体を振動させる外力の時間領域での加速度をｙ_０"(t)、時間遅れ成分ｋ(t_j)の総数をｎとしたときに、前記記憶手段に記憶されているインパルス応答値を、
   [Ｍs]{ｕ"(t)}＋{Ｆ(γ,t)}＝−ｙ_０"(t)[Ｍs]{1} …(１)
   但し、
   
   上式に代入し、剛性低下率α(γ)及び減衰定数ｈ(γ)を各時刻における前記物体の歪レベルγに応じた値に切り替えながら物体の反力ベクトル{Ｆ(γ,t)}、変位ベクトル{ｕ(t)}及び速度ベクトル{ｕ'(t)}をΔt刻みで順次演算することで、前記物体の時刻歴応答解析を行う解析手段と、
   を備えたことを特徴とする時刻歴応答解析装置。
6. 記憶手段と接続されたコンピュータを、減衰定数ｈが物体を振動させる外力の振動数に依存しない振動数非依存特性を示しかつ剛性低下率α及び減衰定数ｈが前記物体の歪レベルγに応じて変化する歪振幅依存特性を示す前記物体の時刻歴応答解析を行う時刻歴応答解析装置として機能させるための時刻歴応答解析プログラムであって、
   前記記憶手段には、振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)−(２ω／ωm)なる値(但しｎは整数)を示す虚数部の正則成分、及び、振動数ωに拘わらず(２ω／ωm)なる値を示す虚数部の特異成分の和で表され、振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)なる値を示す虚数部と、前記虚数部の正則成分のヒルベルト変換値に対応する実数部と、から成る因果的単位虚数関数を時間領域へ変換するか、又は前記虚数部のみを時間領域へ変換することで、前記因果的単位虚数関数のインパルス応答値として演算された、前記物体の速度に依存する同時成分ｃ(t₀)、前記物体の変位に依存する同時成分ｋ(t₀)、前記物体の変位に依存するΔt刻みの時間遅れ成分ｋ(t_j)（但しｊは自然数でt_j＝Δt・j）が記憶されており、
   前記コンピュータを、前記物体の質量マトリクスを[Ｍs]、前記物体の初期剛性マトリクスを[Ｋ₀]、時間領域での物体の変位ベクトルを{ｕ(t)}、速度ベクトルを{ｕ'(t)}、反力ベクトルを{Ｆ(γ,t)}、前記物体を振動させる外力の時間領域での加速度をｙ_０"(t)、時間遅れ成分ｋ(t_j)の総数をｎとしたときに、前記演算手段によって演算されたインパルス応答値を、
   [Ｍs]{ｕ"(t)}＋{Ｆ(γ,t)}＝−ｙ_０"(t)[Ｍs]{1} …(１)
   但し、
   
   上式に代入し、剛性低下率α(γ)及び減衰定数ｈ(γ)を各時刻における前記物体の歪レベルγに応じた値に切り替えながら物体の反力ベクトル{Ｆ(γ,t)}、変位ベクトル{ｕ(t)}及び速度ベクトル{ｕ'(t)}をΔt刻みで順次演算することで、前記物体の時刻歴応答解析を行う解析手段として機能させることを特徴とする時刻歴応答解析プログラム。