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

Description

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

地震時の構造物の挙動や耐震安全性等を解析・評価する地震応答解析は、周波数領域で応答解析を行う周波数応答解析と、時間領域で応答解析を行う時刻歴応答解析とに大別される。地盤の動的剛性は周波数領域の複素関数であるため、地震応答解析では地盤の動的剛性をそのまま利用可能な周波数応答解析も多用されている。しかし、大きなエネルギーが構造物に入力される大地震時等には、そのエネルギーの一部が、構造物を構成する各部材の内部に亀裂を生じさせたり各部材を部分的に塑性化させる等によって消費されると共に、この亀裂発生や部分的な塑性化等に伴い各部材の強度が低下することが繰り返されるというプロセスを経るため、構造物の挙動は非線形性を有している。このため、地震時の構造物の挙動等を高精度に予測解析するためには、時刻歴応答解析により、地震時の各時点での構造物の状態（過去にどのような力が加わり、その力によってどのような状態になっているか）を考慮して各時点での構造物の挙動を解析する必要がある。

多くの材料は、内部減衰の履歴吸収エネルギー(減衰定数ｈ)が振動数ωにあまり依存しない振動数非依存特性（減衰定数ｈが振動数ωに拘わらず一定値を示す特性：図１７(Ａ)を参照）を示すことが知られており、減衰定数ｈが振動数非依存特性を示す減衰を履歴減衰、減衰定数ｈが振動数非依存特性を示す物体の減衰モデルを履歴減衰モデルと称しているが、履歴減衰モデルはヒルベルト(Hilbert)変換対を形成しないために因果律を満たさず、時刻歴応答解析へ厳密に適用することは不可能であることが知られている。なお、ωは一般に円振動数、角振動数などと称されるが、本明細書では振動数と略称する。このため、減衰定数ｈが振動数非依存特性を示す物体の時刻歴応答解析にあたっては、履歴減衰モデルに代えて、履歴減衰を近似的に表す剛性比例型やRayleigh型、歪エネルギー比例型等の減衰モデルが用いられている（例えば非特許文献１を参照）。
勝倉弘，栗本修也，「３章 減衰行列の性質と応答結果に対する減衰の影響」，日本建築学会，建築物の減衰，２０００年，ｐ．４０−６４

剛性比例型の減衰モデルは、本来は振動数ωに拘わらず一定値を示す減衰定数ｈを、例として図１７(Ｂ)に示すように振動数ωに正比例する直線で近似するモデルであり、剛性比例型の減衰モデルを用いた時刻歴応答解析を行うにあたっては、地震動に対する構造物の固有振動数のうち一次モードでの固有振動数（一次振動数）を固有値解析によって求め、一次振動数における減衰定数ｈが目標値に一致するようにパラメータ（時間領域の運動方程式における減衰マトリクス[Ｃs]の値）を設定する事前処理が行われる。剛性比例型の減衰モデルを用いた時刻歴応答解析は、構造物が減衰定数ｈの異なる複数の要素から構成されている場合にも適用可能であると共に、構造物の規模が大きい場合にも適用できるという利点を有しているが、構造物は地震動に対して二次以上の固有モードでも振動するのに対し、剛性比例型の減衰モデルは一次振動数以外の振動数域で減衰定数ｈが目標値に一致しないので減衰定数ｈの振動数非依存性を精度良く表すことができず、時刻歴応答解析の解析精度が低いという問題がある。

また、Rayleigh型の減衰モデルは、例として図１７(Ｃ)に示すように、振動数ωの変化に対する減衰定数ｈの推移を曲線で近似するモデルであり、Rayleigh型の減衰モデルを用いた時刻歴応答解析を行うにあたっては、地震動に対する構造物の固有振動数のうち、例えば一次モードにおける固有振動数(一次振動数)及び二次モードでの固有振動数(二次振動数)を固有値解析によって求め、一次振動数及び二次振動数における減衰定数ｈが目標値に各々一致するようにパラメータ（時間領域の運動方程式における減衰マトリクス[Ｃs]の値）を設定する事前処理が行われる。Rayleigh型の減衰モデルを用いた時刻歴応答解析は、剛性比例型の減衰モデルを用いた時刻歴応答解析よりは解析精度が向上するものの、Rayleigh型の減衰モデルは一次振動数及び二次振動数以外の振動数域で減衰定数ｈが目標値に一致しないので解析精度は十分ではなく、構造物が減衰定数ｈの異なる複数の要素から構成されている場合や構造物の規模が大きい場合への適用にも困難さが増大するという問題がある。

更に、歪エネルギー比例型の減衰モデルは、例として図１７(Ｄ)に示すように、振動数ωの変化に伴って変化する減衰定数ｈの値を全ての固有振動数で一致させるモデルであり、歪エネルギー比例型の減衰モデルを用いた時刻歴応答解析では、時間領域の運動方程式における減衰マトリクス[Ｃs]が、全ての固有モードに対応する多数の要素を含むフルマトリクスとなるため、歪エネルギー比例型の減衰モデルを用いた時刻歴応答解析を行うにあたっては、固有値解析によって地震動に対する全ての固有モードを求め、全ての固有モードにおける固有振動数を求めた後に、各固有モードでの固有振動数における減衰定数ｈが目標値に各々一致するように、上記の減衰マトリクス[Ｃs]の全ての要素の値を設定する事前処理を行う必要がある。このため、歪エネルギー比例型の減衰モデルを用いた時刻歴応答解析は、解析精度は高いものの固有値解析を含む事前処理に長い時間がかかるという問題があり、また、構造物の規模が大きい場合への適用も極めて困難である。これらの各減衰モデルの特徴を纏めると次の表１のようになる。

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

上記目的を達成するために請求項１記載の発明に係る時刻歴応答解析方法は、減衰定数ｈが物体を振動させる外力の振動数に依存しない振動数非依存特性を示す物体の時刻歴応答解析を行う時刻歴応答解析方法であって、振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)−(２ω／ωm)なる値(但しｎは整数)を示す虚数部の正則成分、及び、振動数ωに拘わらず(２ω／ωm)なる値を示す虚数部の特異成分の和で表され、振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)なる値を示す虚数部と、前記虚数部の正則成分よりヒルベルト変換を用いて算出した因果律を満たす実数部と、から成る因果的単位虚数関数を時間領域へ変換することで前記因果的単位虚数関数のインパルス応答値を演算し、前記演算によって得られた前記因果的単位虚数関数のインパルス応答値を用いて前記物体の時刻歴応答解析を行うことを特徴としている。

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

Ｓ(ω)＝Ｋ₀(１＋２ｈ・Ｚ(ω)) …(７)
複素剛性（周波数領域の複素関数）が因果律を満たすためには、複素剛性の実数部と虚数部がヒルベルト(Hilbert)変換対を形成する必要があるが、上記(６)式によって規定される単位虚数関数Ｚ(ω)は、実数部Ｚ_R(ω)と虚数部Ｚ_I(ω)がヒルベルト(Hilbert)変換対を形成しないために因果律を満たさず、時間領域への厳密な変換は不可能である。

ここで本願発明者は、構造物に対する地震応答解析等の時刻歴応答解析では、解析対象の振動数範囲の多くが10Hz以下、原子力発電所のような剛な構造物でも20Hz以下であり、これよりも高い振動数範囲は殆どの構造物の応答に影響を与えないことに着目し、一定の振動数範囲(−ωm＜ω＜ωm)内で単位虚数関数Ｚ(ω)と同様の値を示し((６)式の条件をおおよそ満たし)、かつ因果律を満たすように単位虚数関数を修正すれば、高い解析精度が得られる可能性があることに想到した。そして本願発明者は、単位虚数関数Ｚ(ω)の虚数部Ｚ_I(ω)を、因果律を満たさない成分と因果律を満たす成分の和で表すと共に、因果律を満たさない成分とヒルベルト変換対を形成する実数部を導出すれば、上記の条件を満たす単位虚数関数が得られることに想到した。

このため、請求項１記載の発明では、目的の単位虚数関数の虚数部全体Ｚ'_I(ω)を、振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)−(２ω／ωm)なる値(但しｎは整数)を示す虚数部の正則成分Ｚ'_IR(ω)(図１(Ｂ)参照)と、振動数ωに拘わらず(２ω／ωm)なる値を示す虚数部の特異成分Ｚ'_IS(ω)(図１(Ｃ)参照)の和として規定している(Ｚ'_I(ω)＝Ｚ'_IR(ω)＋Ｚ'_IS(ω))。これにより、虚数部全体Ｚ'_I(ω)は振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)なる値を示し、より詳しくは、例として図２に示すように振動数範囲(ω＞ωm)で３以上の値を示し、図示は省略するが振動数範囲(ω＜−ωm)では−３以下の値を示すことから、振動数範囲(ω＞ωm)及び(ω＜−ωm)で単位虚数関数Ｚ(ω)の虚数部Ｚ_I(ω)の値とは相違する((６)式の条件から外れる)ものの、図１(Ａ)に示すように振動数範囲（０＜ω＜ωm）内で＋１の値を示すと共に、振動数範囲（−ωm＜ω＜０）内で−１の値を示すので、振動数範囲（−ωm＜ω＜ωm）内では単位虚数関数Ｚ(ω) の虚数部Ｚ_I(ω)と同一の値を示す(虚数部Ｚ_I(ω)に関する(６)式の条件を満たす)ことになる。

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

また、一例として振動数ωm＝20Hz（40π）の場合に(８)式によって算出される因果的単位虚数関数Ｚ'(ω)の実数部Ｚ'_R(ω)を、因果的単位虚数関数Ｚ'(ω)の虚数部Ｚ'_I(ω)と共に図３に示す。

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

そして請求項１記載の発明のように、上述した因果的単位虚数関数Ｚ'(ω)を時間領域へ変換することで因果的単位虚数関数Ｚ'(ω)のインパルス応答値を演算し、前記演算によって得られた因果的単位虚数関数Ｚ'(ω)のインパルス応答値を用いて物体の時刻歴応答解析を行った場合、因果的単位虚数関数Ｚ'(ω)自体が、２ｈ倍することで振動数範囲（−ωm＜ω＜ωm）内で履歴減衰モデルと精度良く一致する減衰モデルが得られる関数であることから、従来の各減衰モデルの何れかを用いて時刻歴応答解析を行う場合のように、各固有モードでの固有振動数を求める固有値解析を行い、求めた固有振動数において減衰定数ｈが目標値に一致するようにパラメータを設定する事前処理を行う必要がなくなり、減衰定数ｈが振動数非依存特性を示す物体の時刻歴応答解析を簡単かつ高精度に行うことができる。

また請求項１記載の発明は、時刻歴応答解析を簡単かつ高精度に行うことができることに伴い、解析対象の物体が規模の大きい構造物であっても適用できると共に、解析対象の物体が減衰定数ｈの異なる複数の要素から構成されている構造物である場合にも、個々の要素毎に部材力を演算し、演算した部材力を全要素について重ね合わせる等の演算を行うことで適用することができる。

ところで、複素剛性Ｓ(ω)のＮ個のデータ（互いに異なるＮ種類の振動数ωでの複素剛性Ｓ(ω)の値を表すデータ）からインパルス応答値を演算した場合、インパルス応答の構成成分（インパルス応答値）として、変位に依存する剛性項の同時成分ｋ(t₀)、速度に依存する減衰項の同時成分ｃ(t₀)、加速度に依存する質量項の同時成分ｍ(t₀)が得られると共に、剛性項の時間遅れ成分ｋ(t₁),ｋ(t₂),…,ｋ(t_n)及び減衰項の時間遅れ成分ｃ(t₁),ｃ(t₂),…,ｃ(t_n) がΔt刻みでｎ個（ｎ＝Ｎ−１）得られる。なお、「同時成分」は現在の状態量（変位・速度・加速度）に依存して反力を生じる量(成分)を、「時間遅れ成分」は過去の状態量に依存して反力を生じる量(成分)を意味している。図４には、インパルス応答の各成分を複素剛性Ｓ(ω)の各成分と対応させて示す。複素剛性Ｓ(ω)は上記のインパルス応答値を用いて次の(９)式で表すことができ、複素剛性Ｓ(ω)に対応する時間領域の反力Ｆ(t)は変位ｕ(t)、速度ｕ'(t)及び加速度ｕ"(t)を用いて次の(10)式で表すことができる。

但し、因果的単位虚数関数Ｚ'(ω)は、振動数範囲ω＝０〜ωmの中央に相当する振動数(ωm／２)に関して実数部Ｚ'_R(ω)が線対称、虚数部Ｚ'_I(ω)が点対称の変化を示すため、上記各データのうち、質量項の同時成分ｍ(t₀)及び減衰項の時間遅れ成分ｃ(t₁)〜ｃ(t_n-1)の値は何れも０となる。このため、因果的単位虚数関数Ｚ'(ω)及び対応する時間領域の反力Ｆz(t)は次の(11),(12)式で表される。

従って、請求項１記載の発明における時刻歴応答解析は、具体的には、例えば請求項２に記載したように、因果的単位虚数関数のインパルス応答値として、物体の速度に依存する同時成分ｃ(t₀)、物体の変位に依存する同時成分ｋ(t₀)、物体の変位に依存するΔt刻みの時間遅れ成分ｋ(t_j)（但し、ｊは自然数でt_j＝Δt・j）を演算し、前記物体の質量マトリクスを[Ｍs]、前記物体の剛性マトリクスを[Ｋs]、時間領域での物体の変位をｕ(t)、速度をｕ'(t)、反力をＦz(t)、物体を振動させる外力の時間領域での加速度をｙ₀"(t)、時間遅れ成分ｋ(t_j)の総数をｎとしたときに、求めたインパルス応答値を、
[Ｍs][ｕ"(t)]＋[Ｋs]([ｕ(t)]＋２ｈ[Ｆz(t)])＝−ｙ₀"(t) [Ｍs][1] …(１)
但し、

上式に代入し、物体の反力Ｆz(t)、変位ｕ(t) 及び、速度ｕ'(t)をΔt刻みで順次演算することで行うことができる。

上記(１)式は、本発明に係る因果的単位虚数関数Ｚ'(ω)のインパルス応答値を用いて物体の時刻歴応答解析を行うための時間領域の運動方程式であるが、従来の各減衰モデルの何れかを用いて物体の時刻歴応答解析を行う場合、時間領域の運動方程式は次の(13)式のようになる。
[Ｍs][ｕ"(t)]＋[Ｃs][ｕ'(t)]＋[Ｋs][ｕ(t)]＝−ｙ₀"(t) [Ｍs][1] …(13)
但し、(13)式における [Ｃs]は物体の減衰マトリクスである。(13)式によって時刻歴応答解析を行う場合には、減衰マトリクス[Ｃs]の各要素に、特定の固有モード(剛性比例型は一次モードのみ、Rayleigh型は一次及び二次モード、歪エネルギー比例型は全固有モード）での固有振動数における減衰定数ｈが目標値に各々一致するように値を設定する必要があり、従来の各減衰モデルの何れかを用いて物体の時刻歴応答解析を行う場合には、減衰マトリクス[Ｃs]の各要素に設定する値を求めるために固有値解析を含む事前処理が必要となっていた。

これに対して本発明では、因果的単位虚数関数Ｚ'(ω)自体が、２ｈ倍することで振動数範囲（−ωm＜ω＜ωm）内で履歴減衰モデルと精度良く一致する減衰モデルが得られる関数であるために、(１)式を(９)式と比較しても明らかなように、(13)式における[Ｃs][ｕ'(t)]が(１)式では２ｈ[Ｋs][Ｆz(t)]に置き換わっており、減衰マトリクス[Ｃs]が省略されている。これにより、減衰マトリクス[Ｃs]の各要素に設定する値を求めるための固有値解析等の事前処理を行う必要がなくなるので、減衰定数ｈが振動数非依存特性を示す物体の時刻歴応答解析を簡単に行うことができる。

なお、請求項２記載の発明において、解析対象の物体が減衰定数ｈの異なる複数の要素から構成されている場合の時刻歴応答解析は、具体的には、例えば請求項３に記載したように、個々の要素毎の部材力[Ｆ_K(t)]_E(但し、Ｅは個々の要素を識別するための符号)を次の(３)式に従って各々演算し、
[Ｆ_K(t)]_E＝[Ｋs]_E([ｕ(t)]_E＋２ｈ_E[Ｆz(t)]_E) …(３)
個々の要素毎の部材力[Ｆ_K(t)]_Eを重ね合わせ、その結果を前記(１)式に代入することで行うことができる。これにより、解析対象の物体が減衰定数ｈの異なる複数の要素から構成されている場合であっても、当該物体の時刻歴応答解析を行うことができる。

また本願発明者は、本発明に係る因果的単位虚数関数Ｚ'(ω)を時間領域に変換する際のサンプリング点の数(データ点数)がインパルス応答値の精度に及ぼす影響を確認するために、解析対象振動数範囲(検討振動数域)の上限振動数ωm＝20Hzの因果的単位虚数関数Ｚ'(ω)に対し、次の表２に示すように、互いに異なるデータ点数で複素データのサンプリングを行って時間領域へ変換することでインパルス応答値を各々求めた。変換によって得られたインパルス応答値を剛性項ｋ(t_j)と減衰項ｃ(t_j)(但しｊ＝0,1,…,ｎ)に分けて図５に示す。

図５(Ａ)からも明らかなように、剛性項の同時成分ｋ(t₀)は、複素データの数Ｎ＝４と少なく精度が不足していると推定される「３項モデル」を除いて、各モデル共におよそ０になっている。また、剛性項の時間遅れ成分ｋ(t_j)(但しｊ＝1,…,ｎ)については、何れのモデルでも最初の項(ｔ＝0.05秒)で約-0.6を示し、その後ｔが増加するに従って徐々に値が増加して０に近づく変化を示しており、時間遅れ成分の数ｎが少ないモデル程早く０に近づく傾向がある。また、減衰項については３モデルとも同一で、同時成分ｃ(t₀)のみで構成されている。従って、因果的単位虚数関数Ｚ'(ω)のインパルス応答は、実質的に、図４に太線で囲んで示す剛性項の時間遅れ成分と減衰項の同時成分、すなわち現在の速度と過去の変位に依存する量（成分）で構成されていると見なすことができる。

また、上記で得られた各モデルのインパルス応答値を前出の(11)式に代入することで再現される複素剛性を、インパルス応答値の演算に用いた元の因果的単位虚数関数Ｚ'(ω)上でのデータ点（18項モデルのデータ点）と共に図６に示す。図６(Ａ)に示すように、何れのモデルにおいても再現された複素剛性は同様の変化を示し、元の因果的単位虚数関数Ｚ'(ω)（のデータ点）と精度良く一致している。各モデルから再現された複素剛性の差異は、解析対象振動数範囲０〜ωmの下限振動数(０Hz)付近及び上限振動数(ωm＝20Hz)付近に生じている。図６(Ａ)における０〜４Hzの振動数範囲を拡大して示す図６(Ｂ)からも明らかなように、時間遅れ成分の数ｎが多いモデル程、下限の０Hz近く迄虚数部(図では"Imag"と表記)が１に近い良好な値を示し高い変換精度が得られるものの、実数部(図では"Real"と表記)が０Hz付近で元の因果的単位虚数関数Ｚ'(ω)に対して大きく低下する傾向を示している。

なお、本発明に係る因果的単位虚数関数Ｚ'(ω)を用いるモデルは無次元振動数で定式化可能であり、解析対象振動数範囲の上限振動数ωmは任意に変更可能である。参考までに、解析対象振動数範囲をω＝０〜４Hzとし、解析対象上限振動数ωm＝４Hzの因果的単位虚数関数Ｚ'(ω)を上記と同様に時間領域へ変換することで得られたインパルス応答値を図７に、このインパルス応答値を(11)式に代入することで再現される複素剛性を図８に示す。上記のインパルス応答の演算に用いた複素データの数Ｎは先の18項モデルと同じであり、詳しくはデータ点を0.2,0.4,0.6,…,3.8Hz、時間刻みを0.25秒とした。なお、図７には解析対象振動数範囲ω＝０〜20Hzの18項モデルから得られたインパルス応答値も縦軸・横軸共に５倍に拡大して示しており、図８には解析対象振動数範囲ω＝０〜20Hzの18項モデルのインパルス応答値から再現された複素剛性を横軸のみ1/5に縮小して示している。図７，８からも明らかなように、解析対象振動数範囲ωが相違してもインパルス応答値及び再現された複素剛性の変化は同一であることが理解できる。

また本願発明者は、減衰定数ｈの値の大小が解析精度に及ぼす影響を確認するために、前出の(７)式における係数Ｋ₀をＫ₀＝１とし、単位虚数関数Ｚ(ω)に代えて因果的単位虚数関数Ｚ'(ω)を適用した因果的単位複素剛性Ｓ'₁(ω)を用い(次の(13)式を参照)、
Ｓ'₁(ω)＝１＋２ｈ・Ｚ'(ω) …(13)
上限振動数ωm＝20Hzの因果的単位虚数関数Ｚ'(ω)を、変換精度の高い18項モデル及びより簡便な８項モデルによって時間領域へ変換することで得られたインパルス応答値を、前出の(９)式へ代入することで複素剛性Ｓ(ω)を再現し、この複素剛性Ｓ(ω)を(13)式の因果的単位虚数関数Ｚ'(ω)として(13)式へ代入し、(13)式の減衰定数ｈをｈ＝１〜10％の範囲で変化させた場合の因果的単位複素剛性Ｓ'₁(ω)を比較した。結果を図９に示す。

図９(Ａ)〜(Ｄ)からも明らかなように、虚数部については減衰定数ｈが何れの値の場合にも目標値との一致度が高い結果が得られている。一方、実数部については、減衰定数ｈの値が小さい場合には目標値との一致度が高いものの、減衰定数ｈの値が高くなるに従って目標値との一致度が低下している。これは、(13)式における実数部は(１＋２ｈ・Ｚ'_R(ω))であり、減衰定数ｈの値が大きくなるに従ってＺ'_R(ω)の影響が相対的に大きくなり、目標値(＝1.0)から値が離れるためと推察される。従って本発明は、減衰定数ｈの値が大きくなるに従って解析精度が低下するという特性を有している。なお、例えばコンクリートの減衰定数ｈ＝３％程度、鉄骨の減衰定数＝１％程度であるので、これらの材料から成る物体の時刻歴応答解析を行う際には高い解析精度が得られる。

また、図９に示す結果に基づき、剛性の精度として目標値(＝1.0)に対する因果的単位複素剛性Ｓ'₁(ω)の実数部Re(Ｓ'₁(ω))の比率を、減衰定数の精度として目標とする減衰定数ｈの値に対するｈ'(ω)(次の(14)式参照)の比率を演算した結果を、18項モデルと８項モデルに分けて図１０，１１に示す。この図１０，１１においても、減衰定数ｈの値が大きくなるに従って精度が低下する傾向が見られる。

また図１２には、(Ａ)に減衰定数ｈ＝３％の場合、(Ｂ)に減衰定数ｈ＝10％の場合の振動数範囲０〜２Hz及び18〜20Hzにおける結果を一覧表として示す。図１２では精度の目安として目標値に対する誤差が±10％以上のデータにハッチングを付して示している。減衰定数ｈ＝３％の場合、剛性及び減衰定数の誤差が±10％未満となる振動数範囲は、18項モデルでω＝0.5〜19.0Hz、８項モデルでω＝1.0〜18.5Hzとなっている（解析対象振動数範囲をω＝０〜４Hzとした場合に換算すると、18項モデルでω＝0.1〜3.8Hz、８項モデルでω＝0.2〜3.7Hzに相当する）。また、減衰定数ｈ＝10％の場合、剛性及び減衰定数の誤差が±10％未満となる振動数範囲は、18項モデルでω＝1.5〜18.0Hz、８項モデルでω＝2.0〜18.0Hzとなっている（解析対象振動数範囲をω＝０〜４Hzとした場合に換算すると、18項モデルでω＝0.3〜3.6Hz、８項モデルでω＝0.4〜3.6Hzに相当する）。以上の結果に基づき、解析対象の物体の主要な固有振動数が、解析対象振動数範囲のうち精度が低い下限振動数付近及び上限振動数付近にかからないように、解析対象振動数範囲を適切に定めることで、主要な固有振動数が異なる種々の物体の時刻歴応答解析を精度良く行うことが可能となる。

地震応答解析等の時刻歴応答解析では、前述のように構造物に対する解析対象の振動数範囲の多くがおおよそ定まっていることから、解析対象の振動数範囲が振動数ω＝０〜ωmの範囲内に入るように解析対象上限振動数ωmを予め固定的に定めておくことも可能であるが、上記のように振動数ω＝０〜ωmの範囲内における剛性及び減衰定数の誤差が一定でないことを考慮すると、請求項１又は請求項２記載の発明において、例えば請求項４に記載したように、物体の主要な固有振動数を把握する実固有値解析を行い、実固有値解析によって把握した前記物体の主要な固有振動数が、振動数ω＝０〜ωmの範囲内のうち誤差が所定値未満となる振動数範囲（例えば図１２に示すように剛性及び減衰定数の誤差が10％未満となる振動数範囲）内に入るように、解析対象上限振動数ωmを設定又は選択することが好ましい。これにより、時刻歴応答解析の精度を更に向上させることができる。なお、上記の実固有値解析は、従来の各減衰モデルの何れかを用いて時刻歴応答解析を行う場合の固有値解析のように固有振動数を厳密に求める必要はなく、主要な固有振動数が何れの振動数範囲内に存在しているかを把握できれば目的を達成できるので簡単な処理で済む。

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

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

請求項７記載の発明に係る時刻歴応答解析装置は、減衰定数ｈが物体を振動させる外力の振動数に依存しない振動数非依存特性を示す物体の時刻歴応答解析を行う時刻歴応答解析装置であって、振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)−(２ω／ωm)なる値(但しｎは整数)を示す虚数部の正則成分、及び、振動数ωに拘わらず(２ω／ωm)なる値を示す虚数部の特異成分の和で表され、振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)なる値を示す虚数部と、前記虚数部の正則成分よりヒルベルト変換を用いて算出した因果律を満たす実数部と、から成る因果的単位虚数関数を時間領域へ変換することで得られた前記因果的単位虚数関数のインパルス応答値を記憶する記憶手段と、前記記憶手段に記憶されている前記因果的単位虚数関数のインパルス応答値を用いて前記物体の時刻歴応答解析を行う解析手段と、を備えたことを特徴としているので、請求項１記載の発明と同様に、減衰定数ｈが振動数非依存特性を示す物体の時刻歴応答解析を簡単かつ高精度に行うことができる。

請求項８記載の発明に係る時刻歴応答解析プログラムは、記憶手段を備えたコンピュータを、減衰定数ｈが物体を振動させる外力の振動数に依存しない振動数非依存特性を示す物体の時刻歴応答解析を行う時刻歴応答解析装置として機能させるための時刻歴応答解析プログラムであって、前記記憶手段には、振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)−(２ω／ωm)なる値(但しｎは整数)を示す虚数部の正則成分、及び、振動数ωに拘わらず(２ω／ωm)なる値を示す虚数部の特異成分の和で表され、振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)なる値を示す虚数部と、前記虚数部の正則成分よりヒルベルト変換を用いて算出した因果律を満たす実数部と、から成る因果的単位虚数関数を時間領域へ変換することで得られた前記因果的単位虚数関数のインパルス応答値が記憶されており、前記コンピュータを、前記記憶手段に記憶されている前記因果的単位虚数関数のインパルス応答値を用いて前記物体の時刻歴応答解析を行う解析手段として機能させることを特徴としている。

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

以上説明したように本発明は、振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)−(２ω／ωm)なる値を示す虚数部の正則成分、及び、振動数ωに拘わらず(２ω／ωm)なる値を示す虚数部の特異成分の和で表され、振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)なる値を示す虚数部と、前記虚数部の正則成分よりヒルベルト変換を用いて算出した因果律を満たす実数部と、から成る因果的単位虚数関数を時間領域へ変換することで因果的単位虚数関数のインパルス応答値を演算し、当該因果的単位虚数関数のインパルス応答値を用いて物体の時刻歴応答解析を行うようにしたので、減衰定数ｈが振動数非依存特性を示す物体の時刻歴応答解析を簡単かつ高精度に行うことができる、という優れた効果を有する。

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

ＰＣ１０のＨＤＤ２０には、インパルス応答値データベース（ＤＢ）が記憶されており（詳細は後述）、後述する時刻歴応答解析処理を行うための時刻歴応答解析プログラムもインストールされている。この時刻歴応答解析プログラムは、請求項８記載の発明に係る時刻歴応答解析プログラムに対応している。時刻歴応答解析プログラムをＰＣ１０にインストール（移入）するには幾つかの方法があるが、例えば時刻歴応答解析プログラムをセットアッププログラムと共にＣＤ−ＲＯＭ２２に記録しておき、該ＣＤ−ＲＯＭ２２をＣＤ−ＲＯＭドライブ２４にセットし、ＣＰＵ１０Ａに対して前記セットアッププログラムの実行を指示すれば、ＣＤ−ＲＯＭ２２から時刻歴応答解析プログラムが順に読み出され、読み出された時刻歴応答解析プログラムがＨＤＤ２０に順に書き込まれることで、時刻歴応答解析プログラムのインストールが行われる。

ＰＣ１０は、ＣＰＵ１０Ａが時刻歴応答解析プログラムを実行することで、請求項７記載の発明に係る時刻歴応答解析装置として機能する。なお、ＰＣ１０は請求項８に記載のコンピュータにも対応しているが、請求項８に記載のコンピュータはＰＣ１０に限られるものではなく、例えばワークステーションであってもよいし、汎用の大型コンピュータであってもよい。

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

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

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

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

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

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

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

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

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

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

続いて、上述したインパルス応答値ＤＢに記憶されているインパルス応答値を用いてＰＣ１０によって行われる時刻歴応答解析処理について、図１５を参照して説明する。この時刻歴応答解析処理は請求項７，８に記載の解析手段に相当する処理であり、まずステップ１２０では、解析対象の構造物の剛性・質量のデータとして、解析対象の構造物の質量マトリクス[Ｍs]や剛性マトリクス[Ｋs]等のデータを取得する。次のステップ１２２では、ステップ１２０で取得した解析対象の構造物の質量マトリクス[Ｍs]や剛性マトリクス[Ｋs]に基づき、解析対象の構造物に対して実固有値解析を行い、解析対象の構造物の主要な固有振動数（例えば一次振動数及び二次振動数）と固有モードを算出する。

ところで、構造物に入力する地震動が比較的小さい場合には、構造物に損傷劣化が生じないので構造物の応答は線形範囲内に収まり、構造物の剛性や固有振動数の変化は生じないが、入力する地震動が比較的大きくなると、構造物に損傷劣化が生ずることで構造物の応答が非線形化し、構造物の剛性及び固有振動数が低下する。そして、構造物の応答の非線形化に伴って変化した固有振動数が、後述する時刻歴応答解析における解析対象振動数範囲０〜ωmから外れたり、解析対象振動数範囲０〜ωmのうち精度の低い振動数範囲に入ってしまうと、時刻歴応答解析の途中で処理を継続できない状態に陥ったり、解析精度が低下する恐れがある。

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

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

データ点数に関しては、図１０と図１１を比較しても明らかなように、或る程度以上のデータ点数があればデータ点数を更に増加させても精度はそれ程変化しないので、時刻歴応答解析処理の簡素化を考慮し、上記の条件を満たす範囲内でなるべく少ないデータ点数を選択すればよい。また図１２に示すように、データ点数が多くなるに従い剛性及び減衰定数の誤差が所定値未満となる振動数範囲が若干広くなるので、一部の固有振動数が、剛性及び減衰定数の誤差が所定値未満となる振動数範囲の端部付近に存在している場合には、時刻歴応答解析処理の精度向上のためにデータ点数を増加させるようにしてもよい。

解析対象上限振動数ωm及びデータ点数を選択すると、次のステップ１２８では、選択した解析対象上限振動数ωm及びデータ点数に対応するインパルス応答値(剛性項の同時成分ｋ(t₀)、減衰項の同時成分ｃ(t₀)、及び、剛性項の時間遅れ成分ｋ(t₁)〜ｋ(t_n))をインパルス応答値ＤＢから読み込む。そしてステップ１３０では、インパルス応答値ＤＢから読み込んだインパルス応答値を前出の(12)式((２)式)に代入することで、時刻ｔにおける反力Ｆz(t)を演算し、演算した反力Ｆz(t) 及び構造物を振動させる外力（地震動）の時間領域での加速度ｙ₀"(t)を(１)式に代入して時間領域での物体の変位ｕ(t)、速度ｕ'(t)を演算することを、時刻ｔ＝０から時間Δｔ刻みで繰り返すことで、解析対象の構造物の時刻歴応答解析を行う。なお、時刻歴応答解析を行うにあたり、剛性項の時間遅れ成分ｋ(t₁)〜ｋ(t_n)のうち一部を省略して用いるようにしてもよい。

これにより、解析対象の構造物の時刻歴応答解析を高精度に行うことができる。また、本実施形態に係る時刻歴応答解析処理では、(13)式の減衰マトリクス[Ｃs]の各要素に、特定の固有モードでの固有振動数における減衰定数ｈが目標値に各々一致するように値を設定するための固有値解析を含む事前処理を行う必要もなく、解析処理時間の短縮化や事前処理を行うコンピュータに加わる負荷の低減も実現することができる。

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

また、上記では本発明に係る時刻歴応答解析を、構造物に地震動が入力された際の構造物の応答を解析する地震応答解析に適用した例を説明したが、これに限定されるものではなく、減衰定数ｈが物体を振動させる外力の振動数に依存しない振動数非依存特性を示す任意の物体の時刻歴応答解析に適用可能である。

次に、本発明方式（因果的単位虚数関数を時間領域へ変換することで得られるインパルス応答値を用いた時刻歴応答解析）による解析精度を確認するために本願発明者が行った解析検討の結果について説明する。

この解析検討では図１６に示す４自由度の解析モデルを用い、本発明を含む４種類の方式、すなわち、(1)本発明に係る因果的単位虚数関数から８項モデルによって求めたインパルス応答値を用いる方式、(２)本発明に係る因果的単位虚数関数から18項モデルによって求めたインパルス応答値を用いる方式、(3)剛性比例型の減衰モデルを用いる方式、及び、(4)歪エネルギー比例型の減衰モデルを用いる方式を適用して、時刻歴応答解析（地震応答解析）を各々行った。なお、入力地震動はEl Centro1940NS(時間刻み0.0025秒、継続時間10.00秒)の最大加速度300Galとして入力し、時刻歴応答解析の時間積分法としては平均加速度法を用いた。また、解析精度の比較の基準として用いるために周波数応答解析も行い、この周波数応答解析ではΔｆ＝0.0488Hzとして振動数ω＝20Hzまで解析を行った。

まず剛性比例型及び歪エネルギー比例型の減衰モデルを適用するために、地震応答解析に先立って行った固有値解析の結果を次の表３に示す。

次に、各方式の解析によって得られた最大応答加速度を、周波数応答解析によって得られた最大応答加速度と共に次の表４に示す。

表４に示すように、周波数応答解析によって得られた最大応答加速度に対し、本発明方式(18項モデル及び８項モデル)と歪エネルギー比例型の減衰モデルを用いる方式は、全点で最大応答加速度の差異が±１％以内と高い精度が得られているが、剛性比例型の減衰モデルを用いる方式では最大加速度の差異が大きくなっている。このように、本願発明者が行った解析検討により、本発明は歪エネルギー比例型の減衰モデルを用いた場合と同等の精度で時刻歴応答解析を行えることが確認された。

一方、歪エネルギー比例型の減衰モデルを用いる時刻歴応答解析は、時間領域の運動方程式における減衰マトリクス[Ｃs]が全ての固有モードに対応する多数の要素を含むフルマトリクスとなるため、固有値解析によって地震動に対する全ての固有モードを求め、全ての固有モードにおける固有振動数を求めた後に、各固有モードでの固有振動数における減衰定数ｈが目標値に各々一致するように減衰マトリクス[Ｃs]の全ての要素の値を設定する事前処理を行う必要がある。このため、歪エネルギー比例型の減衰モデルを用いた時刻歴応答解析は事前処理に長い時間がかかり、構造物の規模が大きい場合への適用も極めて困難である。これに対して本発明方式は、固有値解析等の複雑な事前処理を行う必要がなく、構造物の規模が大きい等の場合にも計算負荷が小さい点で歪エネルギー比例型の減衰モデルを用いる時刻歴応答解析よりも優れている。

また、歪エネルギー比例型の減衰モデルを用いる時刻歴応答解析は、解析モデル全体で統一的に固有値毎の減衰定数を付与するものであり、解析モデルの一部のみにこの減衰を用いるような場合には部分構造合成法等の技術が必要になる。これに対して本発明方式では、解析モデルを要素単位で定義することも可能であり(前出の(３)式を参照)、必要に応じて個々の要素毎に解析対象振動数範囲（検討振動数域）や時間領域への変換に用いるデータ点数を変更することも可能であるので、解析処理の自由度が大きいという利点も有している。また、有限要素法のシェル要素やソリッド要素等の任意の要素に適用可能であるという特長も有している。

本発明に係る因果的単位虚数関数の、(Ａ)は虚数部全体、(Ｂ)は正則成分、(Ｃ)は特異成分を各々示す線図である。 虚数部の正則成分と特異成分の和で表される因果的単位虚数関数の虚数部全体の特性を示す線図である。 因果的単位虚数関数Ｚ'(ω)の一例を示す線図である。 インパルス応答の構成成分と複素剛性の成分との対応を示す図表である。 上限振動数ωm＝20Hzの因果的単位虚数関数を時間領域へ変換することで得られるインパルス応答値の一例を示す線図である。 図５のインパルス応答値から再現された複素剛性の一例を示す線図である。 上限振動数ωm＝４Hzの因果的単位虚数関数を時間領域へ変換することで得られるインパルス応答値の一例を示す線図である。 図７のインパルス応答値から再現された複素剛性の一例を示す線図である。 減衰定数ｈの変化に対する因果的単位複素剛性の変化の一例を示す線図である。 18項モデルでインパルス応答値を演算した場合の目標値に対する精度を示す線図である。 ８項モデルでインパルス応答値を演算した場合の目標値に対する精度を示す線図である。 18項モデル及び８項モデルで減衰定数ｈ＝3％,10％の条件でインパルス応答値を演算した場合の剛性及び減衰定数の精度を示す図表である。 本実施形態に係るＰＣの概略構成を示すブロック図である。 インパルス応答値演算処理の内容を示すフローチャートである。 時刻歴応答解析処理の内容を示すフローチャートである。 本願発明者が実施した解析検討で用いた４自由度モデルの、(Ａ)は概念図、(Ｂ)は諸元を示す図表である。 (Ａ)は多くの材料の振動数ω−減衰定数ｈ特性を示す線図、(Ｂ)は剛性比例型、(Ｃ)はRayleigh型、(Ｄ)は歪エネルギー比例型の各減衰モデルにおける振動数ω−減衰定数ｈ特性を各々示す線図である。

符号の説明

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

Claims (8)

1. 減衰定数ｈが物体を振動させる外力の振動数に依存しない振動数非依存特性を示す物体の時刻歴応答解析を行う時刻歴応答解析方法であって、
   振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)−(２ω／ωm)なる値(但しｎは整数)を示す虚数部の正則成分、及び、振動数ωに拘わらず(２ω／ωm)なる値を示す虚数部の特異成分の和で表され、振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)なる値を示す虚数部と、前記虚数部の正則成分よりヒルベルト変換を用いて算出した因果律を満たす実数部と、から成る因果的単位虚数関数を時間領域へ変換することで前記因果的単位虚数関数のインパルス応答値を演算し、
   前記演算によって得られた前記因果的単位虚数関数のインパルス応答値を用いて前記物体の時刻歴応答解析を行うことを特徴とする時刻歴応答解析方法。
2. 前記因果的単位虚数関数のインパルス応答値として、物体の速度に依存する同時成分ｃ(t₀)、物体の変位に依存する同時成分ｋ(t₀)、物体の変位に依存するΔt刻みの時間遅れ成分ｋ(t_j)（但し、ｊは自然数でt_j＝Δt・j）を演算し、
   前記物体の質量マトリクスを[Ｍs]、前記物体の剛性マトリクスを[Ｋs]、時間領域での物体の変位をｕ(t)、速度をｕ'(t)、反力をＦz(t)、物体を振動させる外力の時間領域での加速度をｙ₀"(t)、時間遅れ成分ｋ(t_j)の総数をｎとしたときに、求めたインパルス応答値を、
   [Ｍs][ｕ"(t)]＋[Ｋs]([ｕ(t)]＋２ｈ[Ｆz(t)])＝−ｙ₀"(t) [Ｍs][1] …(１)
   但し、
   

   

   上式に代入し、物体の反力Ｆz(t)、変位ｕ(t) 及び、速度ｕ'(t)をΔt刻みで順次演算することで、前記物体の時刻歴応答解析を行うことを特徴とする請求項１記載の時刻歴応答解析方法。
3. 解析対象の物体が、減衰定数ｈの異なる複数の要素から構成されている場合に、個々の要素毎の部材力[Ｆ_K(t)]_E(但し、Ｅは個々の要素を識別するための符号)を次の(３)式に従って各々演算し、
   [Ｆ_K(t)]_E＝[Ｋs]_E([ｕ(t)]_E＋２ｈ_E[Ｆz(t)]_E) …(３)
   個々の要素毎の部材力[Ｆ_K(t)]_Eを重ね合わせ、その結果を前記(１)式に代入することで、前記解析対象の物体の時刻歴応答解析を行うことを特徴とする請求項２記載の時刻歴応答解析方法。
4. 物体の主要な固有振動数を把握する実固有値解析を行い、前記実固有値解析によって把握した前記物体の主要な固有振動数が、振動数ω＝０〜ωmの範囲内のうち誤差が所定値未満となる振動数範囲内に入るように解析対象上限振動数ωmを設定又は選択することを特徴とする請求項１又は請求項２記載の時刻歴応答解析方法。
5. 物体を振動させる外力と前記物体の挙動との関係の非線形化によって前記物体の固有振動数が変化するか否かを推定し、前記物体の固有振動数が変化すると判断した場合には、概略の減衰に基づく予備解析により前記非線形化後の固有振動数をおおよそ把握し、把握した前記非線形化後の固有振動数も振動数ω＝０〜ωmの範囲内のうち誤差が所定値未満となる振動数範囲内に入るように解析対象上限振動数ωmを設定又は選択することを特徴とする請求項４記載の時刻歴応答解析方法。
6. 振動数ωmの値が互いに異なる複数種の因果的単位虚数関数について、時間領域への変換を各々行い前記複数種の因果的単位虚数関数のインパルス応答値を各々演算して記憶手段に記憶しておき、
   前記物体の時刻歴応答解析に際し、前記記憶手段に記憶した前記複数種の因果的単位虚数関数のインパルス応答値のうち、前記把握した固有振動数が振動数ω＝０〜ωmの範囲内のうち誤差が所定値未満となる振動数範囲内に入る特定の因果的単位虚数関数のインパルス応答値を読み出して用いることを特徴とする請求項４又は請求項５記載の時刻歴応答解析方法。
7. 減衰定数ｈが物体を振動させる外力の振動数に依存しない振動数非依存特性を示す物体の時刻歴応答解析を行う時刻歴応答解析装置であって、
   振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)−(２ω／ωm)なる値(但しｎは整数)を示す虚数部の正則成分、及び、振動数ωに拘わらず(２ω／ωm)なる値を示す虚数部の特異成分の和で表され、振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)なる値を示す虚数部と、前記虚数部の正則成分よりヒルベルト変換を用いて算出した因果律を満たす実数部と、から成る因果的単位虚数関数を時間領域へ変換することで得られた前記因果的単位虚数関数のインパルス応答値を記憶する記憶手段と、
   前記記憶手段に記憶されている前記因果的単位虚数関数のインパルス応答値を用いて前記物体の時刻歴応答解析を行う解析手段と、
   を備えたことを特徴とする時刻歴応答解析装置。
8. 記憶手段と接続されたコンピュータを、減衰定数ｈが物体を振動させる外力の振動数に依存しない振動数非依存特性を示す物体の時刻歴応答解析を行う時刻歴応答解析装置として機能させるための時刻歴応答解析プログラムであって、
   前記記憶手段には、振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)−(２ω／ωm)なる値(但しｎは整数)を示す虚数部の正則成分、及び、振動数ωに拘わらず(２ω／ωm)なる値を示す虚数部の特異成分の和で表され、振動数ωが(ｎ−１)・ωmからｎ・ωmの範囲で(２ｎ−１)なる値を示す虚数部と、前記虚数部の正則成分よりヒルベルト変換を用いて算出した因果律を満たす実数部と、から成る因果的単位虚数関数を時間領域へ変換することで得られた前記因果的単位虚数関数のインパルス応答値が記憶されており、
   前記コンピュータを、
   前記記憶手段に記憶されている前記因果的単位虚数関数のインパルス応答値を用いて前記物体の時刻歴応答解析を行う解析手段
   として機能させることを特徴とする時刻歴応答解析プログラム。

Images