時刻歴応答シミュレーション方法、そのプログラム及び装置
【課題】空力弾性振動の時刻歴応答解析シミュレーション方法を提供すること。
【解決手段】弾性振動体が受ける非定常空気力の振動変位成分ｆaR及び振動変位速度成分ｆaIを変数として含む振動方程式を用いた、弾性振動体の時刻歴応答シミュレーション方法であって、弾性振動体の所定部位の時刻ｔの変位及び変位速度の値から時刻ｔの振幅ｙ0を計算し、振幅ｙ0の値をから空力減衰率δa及び円振動数ｎの値を計算し、空力減衰率δaの値から振動変位速度成分ｆaIを計算すると共に円振動数ｎの値から振動変位成分ｆaRを計算し、振動変位速度成分ｆaI及び振動変位成分ｆaRの値が代入された振動方程式を解いて、時刻ｔ+Δｔの所定部位の変位、変位速度及び加速度を計算する時刻歴応答計算ルーチンを繰り返し実行する方法である。
【解決手段】弾性振動体が受ける非定常空気力の振動変位成分ｆaR及び振動変位速度成分ｆaIを変数として含む振動方程式を用いた、弾性振動体の時刻歴応答シミュレーション方法であって、弾性振動体の所定部位の時刻ｔの変位及び変位速度の値から時刻ｔの振幅ｙ0を計算し、振幅ｙ0の値をから空力減衰率δa及び円振動数ｎの値を計算し、空力減衰率δaの値から振動変位速度成分ｆaIを計算すると共に円振動数ｎの値から振動変位成分ｆaRを計算し、振動変位速度成分ｆaI及び振動変位成分ｆaRの値が代入された振動方程式を解いて、時刻ｔ+Δｔの所定部位の変位、変位速度及び加速度を計算する時刻歴応答計算ルーチンを繰り返し実行する方法である。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、時刻歴応答解析に関するものであり、特に空力弾性振動の時刻歴応答解析シミュレーション方法、そのプログラム及び装置に関するものである。
【背景技術】
【０００２】
時刻歴応答解析は、耐震工学の分野では古い歴史がある。耐震工学の分野における時刻歴応答解析とは、地震振動（地震加速度）などの外力が与えられた構造物の所定位置の応答加速度、速度、変位などを計算することである。そして、計算された解析結果は、構造物の耐震性の判断に用いられている。
【０００３】
ところで、橋梁の耐風工学の分野では、橋梁、特に長大橋梁の建設の際、耐風特性の検討が行われる（特許文献１参照）。耐風特性の検討では、例えば、風に代表される外力（空気力）を受ける橋梁の所定位置の応答（変位）が検討される。そして、検討結果が橋梁の耐風設計に用いられる。このような耐風特性の検討のために、時刻歴応答解析を用いることが考えられるが、現実には時刻歴応答解析を提案する流れはなかった。
【０００４】
時刻歴応答解析を行うためには、構造物に作用する外力の特性を十分に把握して評価できる必要があるところ、橋梁等の構造物が受ける外力である空気力の特性を把握することは困難であるということが原因であると考えられる。つまり、空力弾性振動のような自励振動における空気力の評価は、振動振幅に対して行わざるを得ず、過渡応答的な振動における各瞬間の空気力を評価することは困難であるため、耐風特性の検討に時刻歴応答解析を用いることは困難であると考えられた。
【０００５】
このようなことから、現在、橋梁全体の挙動は、一般に全橋模型を用いた風洞実験によって把握されており、橋梁の耐風特性もこの全橋模型を用いた風洞実験によって検討されている。
【特許文献１】特開平８−１８４５２５号公報
【発明の開示】
【発明が解決しようとする課題】
【０００６】
ところが、風洞実験の場合、構造物の構造特性を種々変更して実験を行うには限界がある。そして、風洞実験では、乱流中での構造物の応答挙動を推定することが非常に困難であり、耐風特性の検討に関しても統計的な検討にならざるを得ないという問題がある。また、風洞実験に使用できる空間は限られた空間であるので、風洞実験では、実物の橋梁に対して大幅に縮小された全橋模型を使用せざるを得ず、実験精度が低下しやすいという問題もある。
【０００７】
なお、対風応答に関するシミュレーションとしては、有限要素法（ＦＥＭ）を用いたフラッタ解析がある。フラッタとは、元来、航空機の翼などの構造物に作用する外力（空力）、弾性力及び慣性力が組み合わさって生じる自励振動のことである。フラッタ解析によれば、このようなフラッタ（自励振動）の発生の有無を判断することはできるが、外力（例えば、空気力、減衰力、弾性力）を受けた構造物のフラッタ（自励振動）の振幅（応答振幅）を追跡することはできない。
【０００８】
本発明は、上記問題点に鑑みてなされたものであり、空力弾性振動の時刻歴応答解析シミュレーション方法およびその装置を提供することを課題とする。
【課題を解決するための手段】
【０００９】
本願発明の発明者は、上記課題に鑑み、減衰状態や発散状態における構造物の時々刻々の振幅に基づいて非定常空気力（瞬時作用空気力）を求める手法について、後述するとおり鋭意研究を重ねた結果、次のような本願発明に想到するに至った。
【００１０】
つまり、本願発明は、外力として非定常空気力を受ける弾性振動体の時刻歴応答シミュレーション方法であって、非定常空気力を受ける前記弾性振動体の所定部位の時刻ｔにおける変位ｙt及び変位速度ｙ’tの値を用いて前記所定部位の振幅ｙ0の時刻ｔにおける値を計算する第１ステップと、振幅ｙ0の時刻ｔにおける値を用いて、振幅ｙ0の関数である空力減衰率δa及び応答円振動数ｎの時刻ｔにおける値を計算する第２ステップと、空力減衰率δaの時刻ｔにおける値を用いて前記非定常空気力の振動変位速度成分ｆaIの時刻ｔにおける値を計算すると共に、応答円振動数ｎの時刻ｔにおける値を用いて前記非定常空気力の振動変位成分ｆaRの時刻ｔにおける値を計算する第３ステップと、振動変位成分ｆaR及び振動変位速度成分ｆaIを変数として含む弾性振動体の振動方程式に振動変位速度成分ｆaI及び振動変位成分ｆaRの時刻ｔにおける値を代入して得られた振動方程式を解いて、次の時刻ｔ+Δｔにおける前記所定部位の変位ｙt+Δt、変位速度ｙ’t+Δt及び加速度ｙ’’t+Δtの値を計算する第４ステップと、を備え、これらのステップを含む時刻歴応答計算ルーチンを各時刻について繰り返し実行して計算された値を基に弾性振動体の時刻歴応答シミュレーションを行うものであり、２回目以降の時刻歴応答計算ルーチンの第１ステップでは、当該第１ステップにおける変位ｙt及び変位速度ｙ’tの値として、直前の時刻歴応答計算ルーチンの第４ステップで計算された変位ｙ及び変位速度ｙ’の値を用いて振幅ｙ0の対応時刻における値を計算するものである。
【００１１】
本発明のシミュレーション方法は、基本的には、次の式（１）に示される振動方程式に係るものである。そして、式（１）は、式（２）、（３）を適用することにより、式（４）のように表される。
【数２】
【００１２】
ここで、非定常空気力の振動変位速度成分ｆaI及び非定常空気力の振動変位成分ｆaRを、次の式（５）のようにおく。そうすると、上記式（１）の振動方程式を、下記の式（６）のように表すことができる。このように、振動変位速度成分ｆaI及び非定常空気力の振動変位成分ｆaRを式（５）のようにおいて振動方程式を式（６）のように表す考え方は好ましいものである。
【数３】
【００１３】
ところで、式（６）では、非定常空気力の振動変位速度成分ｆaIが減衰項の一部を構成し、非定常空気力の振動変位成分ｆaRが復元力項の一部を構成している。この振動方程式を解くためには、振動変位速度成分ｆaI及び振動変位成分ｆaRの具体的な数値を把握する必要がある。
【００１４】
そこで、非定常空気力の振動変位速度成分ｆaI及び振動変位成分ｆaRの具体的な数値を把握する方法について検討した。その結果、空力弾性振動体の所定部位の空力減衰率δaが解かれば振動変位速度成分ｆaIの具体的数値を把握することができ、空力弾性振動体の所定部位の応答円振動数ｎが解かれば振動変位成分ｆaRの具体的数値を把握することができることが解かった。ただし、空力弾性振動体に加わる外力は時々刻々と変化する非定常空気力であるので、簡単に所定部位の空力減衰率δa及び応答円振動数ｎを把握することはできない。
【００１５】
そこで、空力減衰率δa及び応答円振動数ｎを把握する方法について、さらに検討し、その結果、振動体の所定部位に作用する各時刻の非定常空気力は、時々刻々と変化する所定部位の振幅ｙ0（図７参照）に対応しているという考え方を用いることによって、空力弾性振動の応答を追跡できることを見出した。より具体的に言えば、空力減衰率δa及び応答円振動数ｎの具体的数値を、時々刻々と変化する所定部位の振幅ｙ0に基づいて把握することができることを見出すと共に、各時刻における振幅ｙ0を、当該時刻における変位及び変位速度に基づいて求めることができることを見出した。これを実際のシミュレーション方法の処理の流れに沿って説明すれば、各時刻における変位及び変位速度を用いて当該時刻の振幅ｙ0を求め、当該振幅ｙ0から空力減衰率δa及び応答円振動数ｎの具体的数値を把握し、これらを用いて非定常空気力の振動変位速度成分ｆaI及び振動変位成分ｆaRの具体的数値を把握し、把握した非定常空気力に関する具体的数値を用いてシミュレーションに用いる振動方程式を確定することができることを見出したということになる。
なお、時刻ｔにおける所定部位の振幅ｙ0とは、各時刻ｔにおける変位ｙtではなく、その時刻ｔにおいて取り得た最大変位のことである。振幅ｙ0は、図７の破線で示されるように変化する。
これらのことを見出したことで、最終的に、本願発明に係る空力弾性振動体の時刻歴応答シミュレーション方法に想到するに至った。
【００１６】
本発明では、まず、ある時刻の変位及び変位速度に基づいて当該時刻の振幅ｙ0を求めるために、非定常空気力（外力）を受けた空力弾性振動体の所定部位の振動は、正弦波的な振動で減衰しあるいは発散するものであると考える。このように考えると、空力弾性振動の変位ｙを、下記の式（７）のように表すことができ、変位速度ｙ’を、下記の式（８）のように表すことができる。そして、瞬間瞬間の振動に対応する振幅（正弦波振幅）ｙ0を各時刻ｔにおける「変位」の自乗と「変位速度／振動数」の自乗の和の平方根で表現することができるということを見出して、式（７），（８）から式（９）を得た。さらに、橋梁の対風応答を検討する場合、減衰定数ξは０．０１未満であるので、式（９）中のルートの中は、ほぼ「１」と考えることができ、式（９）中の「e」の乗数（-ξωt）はほぼ「０」と考えることができる。そうすると、空力弾性振動の式（９）から式（１０）が得られる。
【数４】
【００１７】
ここで、上記の式（７）から式（１０）に関連して図７を参照する。図７において、実線は、空力弾性振動体の所定部位の時刻ｔにおける変位ｙを示すものである。そして、破線は、上記式（１０）によって計算される時々刻々の振幅ｙ0である。つまり、式（１０）を用いることによって、弾性振動体の所定部位の各時刻における変位ｙt及び変位速度ｙ’tの値から当該所定部位の振幅ｙ0の当該時刻における値を計算することができる（第１ステップ）。
【００１８】
図７から解かるように、空力弾性振動体の所定部位の変位が振幅に一致するとき（変位速度が０のとき）は、式（１０）がなくても、振幅ｙ0を求めることができる。ところが、振幅ｙ0が時々刻々と過渡応答的に変化する空力弾性振動の場合、式（１０）がなければ、変位速度が０のとき以外の時刻について、振幅ｙ0を求めることはできない。
空力弾性振動のような自励振動における空気力（非定常空気力）の評価は、振動振幅に対して行わざるを得ないものであるが、これまで、過渡応答的な振動における各瞬間の振幅ｙ0を評価することができなかったので、振幅ｙ0の変化に対応しつつ時々刻々の非定常空気力を評価することはできなかった。このようなことから、これまでは、時刻歴応答解析を用いて耐風特性を検討することはできず、空力弾性振動の時刻歴応答シミュレーションを行うことができなかった。
【００１９】
この点、本発明によれば、上記式（１０）を用いることによって、空力弾性振動体の所定部位の時刻ｔにおける変位ｙt及び変位速度ｙ’tに基づき、時々刻々と変化する振幅ｙ0を、いずれの時刻においても求めることができる（第１ステップ）。
なお、変位ｙの自乗値と、変位速度ｙ’を応答円振動数ｎで除した値の自乗値との和の平方根の値に基づいて振幅ｙ0を求める方法であれば、種々の方法によって振幅ｙ0を求めることができるが、式（１０）を用いる方法は、簡単であり好ましい。
【００２０】
次に、先に算出した振幅ｙ0の時刻ｔにおける値を用いて、振幅ｙ0の関数である空力減衰率δa及び応答円振動数ｎの時刻ｔにおける値を計算する（第２ステップ）。なお、振幅ｙ0の関数である空力減衰率δa（＝ｆI（ｙ0））及び振幅ｙ0の関数である応答円振動数ｎ（＝ｆR（ｙ0））は、例えば弾性振動体に対して振動実験を行うことによって求められる。
【００２１】
続いて、先に算出された空力減衰率δaの時刻ｔにおける値を用いて前記非定常空気力の振動変位速度成分ｆaIの時刻ｔにおける値を計算すると共に、先に算出された応答円振動数ｎの時刻ｔにおける値を用いて前記非定常空気力の振動変位成分ｆaRの時刻ｔにおける値を計算する（第３ステップ）。
【００２２】
なお、空力減衰率δaの具体的数値を用いて非定常空気力（瞬時作用空気力）の振動変位速度成分ｆaIを計算する方法は、種々考えられるかも知れないが、本発明において振動変位速度成分ｆaIを計算する際に用いる空力減衰率δaの関数（式（１６））は、次にようにして求められる。つまり、式（６）として示される振動方程式中の振動変位速度成分ｆaIが減衰項の一部であり、振動方程式における定常振動の条件が次の式（１１）であり、振幅−空力減衰率図における定常振動の条件が式（１２）であることに基づき、求められる。
【数５】
【００２３】
上記式（１１）、（１２）によって示される２つの関係が共通した現象における条件である。ここで、−ｆaI／ｍ及びδaは空気力に起因するものであり、２ξωがδsに対応（２ξω＝δs）し、また、−ｆaI／ｍがδaに対応する。そして、構造減衰率δs＝２πξであるので、これから導かれる「ξ＝δs／２π」を用いて「２ξω」を変形すると、次式（１３）のようになり、式（１１）及び（１２）の関係から次の式（１４）が導かれる。
【数６】
【００２４】
したがって、式（６）の減衰項は、次の式（１５）のようになり、δaを振動変位速度成分ｆaIの式になるように変形すると、式（１６）のようになる。これにより、二次元剛体モデルとして表される橋梁の振動方程式の振動変位速度成分ｆaIの具体的な形が求められたことになる。また、あらかじめ得た空力減衰率δaを表す近似関数を用いて各時刻における空力減衰率δaを計算することができるので、計算した空力減衰率δaの値を式（１６）に代入することで、その時刻の振幅ｙ0に対応した振動変位速度成分ｆaIが得られる。
【数７】
【００２５】
また、応答円振動数ｎの具体的数値を用いて非定常空気力（瞬時作用空気力）の振動変位成分ｆaRを計算する方法は、種々考えられるかも知れないが、本発明において振動変位成分ｆaRを計算する際に用いられる応答円振動数ｎの関数（式（１８））は、次にようにして求められる。つまり、式（６）として示される振動方程式中の振動変位成分ｆaRが復元力項の一部であり、復元力項の係数が橋梁モデルの応答円振動数の２乗に等しいと考えると、次の式（１７）が得られる。そして、この式を変形することで、次の式（１８）が得られる。また、あらかじめ得た振動数ｎを表す近似関数を用いて各時刻における振動数ｎを計算することができるので、計算した振動数ｎの値を式（１８）に代入することで、その時刻の振動変位成分ｆaRが得られる。
【数８】
【００２６】
このようにして求められた振動変位成分ｆaR及び振動変位速度成分ｆaIの時刻ｔにおける具体的数値を、振動変位成分ｆaR及び振動変位速度成分ｆaIを変数として含む弾性振動体の振動方程式に代入して得られた振動方程式を解いて、次の時刻ｔ+Δｔにおける前記所定部位の変位ｙt+Δt、変位速度ｙ’t+Δt及び加速度ｙ’’t+Δtの値を計算する（第４ステップ）。これにより、時刻歴応答シミュレーションに用いられる数値が計算される。
【００２７】
本願発明に係る空力弾性振動の時刻歴応答解析シミュレーション方法では、これらの第１ステップから第４ステップを含む時刻歴応答計算ルーチンを各時刻について繰り返し実行し、計算された値を基に弾性振動体の時刻歴応答シミュレーションを行うものである。なお、２回目以降の時刻歴応答計算ルーチンの第１ステップでは、当該第１ステップにおける変位ｙt及び変位速度ｙ’tの値として、直前の時刻歴応答計算ルーチンの第４ステップ（次の時刻ｔ+Δｔにおける各値を計算するステップ）で計算された変位ｙ及び変位速度ｙ’の値を用いて振幅ｙ0の対応時刻における値を計算する。また、２回目以降の時刻歴応答計算ルーチンの第４ステップは、当該第４ステップにおける変位ｙt、変位速度ｙ’t及び加速度ｙ’’tの値として、直前の時刻歴応答計算ルーチンの第４ステップ（次の時刻ｔ+Δｔにおける各値を計算するステップ）で計算された変位ｙ、変位速度ｙ’及び加速度ｙ’’の値を用いて、次の時刻ｔ+Δｔにおける前記所定部位の変位ｙt+Δt、変位速度ｙ’t+Δt及び加速度ｙ’’t+Δtの値を計算する。
【００２８】
このように、本願発明によれば、減衰状態あるいは発散状態にある所定部位の空力弾性振動について、時々刻々の振幅ｙ0を計算することができ、計算された振幅ｙ0を用いて非定常空気力（より具体的にはの振動変位成分ｆaR及び振動変位速度成分ｆaI）を計算することができる。時々刻々の非定常空気力を計算することができれば、計算された非定常空気力を振動方程式に代入することによって時々刻々の振動方程式を定め、定められた振動方程式を解くことによって、空力弾性振動の時刻歴応答解析シミュレーションをすることができる。
【００２９】
なお、前記振幅ｙ0の関数である空力減衰率δa（＝ｆI（ｙ0））及び応答円振動数ｎ（＝ｆR（ｙ0））は、前記空力弾性振動体について行う振動実験によって求めたものが好ましい。例えば、シミュレーションの対象の空力弾性振動体について自由振動実験等の振動実験を行うことによって、上記関数（δa＝ｆI（ｙ0））及び関数（ｎ＝ｆR（ｙ0））を求めることができる。
【００３０】
このように、シミュレーション対象の弾性振動体について実際に実験を行うことによって、上記関数である空力減衰率δa（＝ｆI（ｙ0））及び応答円振動数ｎ（＝ｆR（ｙ0））を正確に求めることができる。そして、このような関数があれば、時々刻々と変化する空力減衰率や応答円振動数を、迅速、簡単かつ確実に求めることができることとなり、シミュレーションをより迅速、簡単かつ確実に実行することができる。
【００３１】
また、前記空力減衰率δa（＝ｆI（ｙ0））及び応答円振動数ｎ（＝ｆR（ｙ0））は、前記振幅ｙ0及び風速の関数であることがより好ましい。振幅ｙ0だけでなく風速が加味された関数を用いることにより、より正確なシミュレーションを行うことができる。
【００３２】
また、本願の別の発明は、外力として非定常空気力を受ける弾性振動体の時刻歴応答シミュレーションをするためのコンピュータを、シミュレーションの対象である前記弾性振動体の所定部位の条件を設定する手段と、風速条件を設定する手段、
前記弾性振動体の所定部位の初期時刻における変位ｙ及び変位速度ｙ’の値を設定する手段、前記弾性振動体の所定部位の時刻ｔにおける変位ｙt及び変位速度ｙ’tの値を用いて当該所定部位の振幅ｙ0の時刻ｔにおける値を計算する手段、
前記風速条件に基づいて特定された時刻ｔにおける風速と、振幅ｙ0の時刻ｔにおける値に基づいて空力減衰率δa及び応答円振動数ｎの時刻ｔにおける値を計算する手段、空力減衰率δaの時刻ｔにおける値を用いて前記非定常空気力の振動変位速度成分ｆaIの時刻ｔにおける値を計算すると共に、応答円振動数ｎの時刻ｔにおける値を用いて前記非定常空気力の振動変位成分ｆaRの時刻ｔにおける値を計算する手段、振動変位成分ｆaR及び振動変位速度成分ｆaIを変数として含む弾性振動体の振動方程式に、振動変位速度成分ｆaI及び振動変位成分ｆaRの時刻ｔにおける値を代入して得られた振動方程式を解いて、次の時刻ｔ+Δｔにおける前記所定部位の変位ｙt+Δt、変位速度ｙ’t+Δt及び加速度ｙ’’t+Δtの値を計算する手段、前記計算する手段を順次実行する時刻歴応答計算ルーチンを各時刻について繰り返し実行させる手段であって、２回目以降の時刻歴応答計算ルーチンの振幅ｙ0の値の計算で用いられる変位ｙt及び変位速度ｙ’tの値として、直前の時刻歴応答計算ルーチンで計算された変位ｙ及び変位速度ｙ’の値を用いて振幅ｙ0の値を計算させる手段、前記時刻歴応答計算ルーチンを各時刻について繰り返し実行して計算された値を用いて弾性振動体の時刻歴応答シミュレーションを出力する手段として機能させる時刻歴応答シミュレーションプログラムである。このようなプログラムによって上述したような時刻歴応答シミュレーションを実現することができ、当該プログラムがインストールされたコンピュータなどの時刻歴応答シミュレーション装置によって上述したような時刻歴応答シミュレーションを実現することができる。
【００３３】
そして、本願に係るさらに別の発明は、外力として非定常空気力を受ける弾性振動体の時刻歴応答シミュレーションを行う時刻歴応答シミュレーション装置であって、非定常空気力を受ける前記弾性振動体の所定部位の時刻ｔにおける変位ｙt及び変位速度ｙ’tの値を用いて前記所定部位の振幅ｙ0の時刻ｔにおける値を計算する第１計算手段と、振幅ｙ0の時刻ｔにおける値を用いて、振幅ｙ0の関数である空力減衰率δa及び応答円振動数ｎの時刻ｔにおける値を計算する第２計算手段と、空力減衰率δaの時刻ｔにおける値を用いて前記非定常空気力の振動変位速度成分ｆaIの時刻ｔにおける値を計算すると共に、応答円振動数ｎの時刻ｔにおける値を用いて前記非定常空気力の振動変位成分ｆaRの時刻ｔにおける値を計算する第３計算手段と、振動変位成分ｆaR及び振動変位速度成分ｆaIを変数として含む弾性振動体の振動方程式に、振動変位速度成分ｆaI及び振動変位成分ｆaRの時刻ｔにおける値を代入して得られた振動方程式を解いて、次の時刻ｔ+Δｔにおける前記所定部位の変位ｙt+Δt、変位速度ｙ’t+Δt及び加速度ｙ’’t+Δtの値を計算する第４計算手段と、前記計算する手段を順次実行する時刻歴応答計算ルーチンを各時刻について繰り返し実行させる手段であって、２回目以降の時刻歴応答計算ルーチンの振幅ｙ0の値の計算で用いられる変位ｙt及び変位速度ｙ’tの値として、直前の時刻歴応答計算ルーチンで計算された変位ｙ及び変位速度ｙ’値を用いて振幅ｙ0の値を計算させる手段と、前記時刻歴応答計算ルーチンを各時刻について繰り返し実行して計算された値を用いて弾性振動体の時刻歴応答シミュレーションを出力する手段と、を備えるものである。このような時刻歴応答シミュレーション装置によって上述したような時刻歴応答シミュレーションを実現することができる。
【発明の効果】
【００３４】
本願発明に係る空力弾性振動体の時刻歴応答シミュレーション方法、そのプログラム及び装置によれば、橋梁などの空力弾性振動体の所定部位の空力減衰率δaを表す関数（δa＝ｆ（ｙ0））と、応答円振動数ｎを表す関数（ｎ＝ｆ（ｙ0））とをあらかじめ求めておくことにより、非定常空気力（外力）を受ける当該構造物の耐風特性（例えば、空力弾性振動）の検討を、時刻歴応答解析シミュレーションを用いて行うことができる。
【００３５】
そして、風洞実験を行うことなくシミュレーションによって空力弾性振動の応答推定を行うことができるので、風洞実験を行うことが難しい長大橋梁等の構造物について、空力弾性振動の応答推定を容易、迅速に行うことができる。
【００３６】
また、本願発明によれば、構造物の構造様式などの条件や構造物が受ける風（風速や風速の変動）など、シミュレーションのパラメータの設定を容易に変更することができるので、多様な条件下のシミュレーションを容易、迅速に行うことができる。これにより、多様な条件下での構造物の空力弾性振動の応答推定を、容易、迅速に行うことができる。
【００３７】
さらに、本願発明によれば、時刻歴応答シミュレーションを用いるので、耐震設計と同等の精度で、耐風設計を行うことが可能である。つまり、構造物の空力弾性振動の応答推定を行う際に、本願発明に係る空力弾性振動体の時刻歴応答シミュレーション方法を利用することによって、より正確な耐風設計を行うことが可能である。これにより、構造物の耐風特性に関する設計の品質を飛躍的に向上させることができる。
【００３８】
また、時刻ｔにおける変位ｙtの自乗値と、時刻ｔにおける変位速度ｙ’tを応答円振動数ｎで除した値の自乗値との和の平方根を求めることによって時刻ｔにおける所定部位の振幅ｙ0tの値を求めることができれば、より容易、迅速かつ正確に時刻歴応答シミュレーションを行うことができる。
【００３９】
また、シミュレーションの各ステップでの計算に用いることができる具体的な関数が用意されているので、それらの関数を用いることによって、容易、迅速かつ正確にシミュレーションを行うことができる。
【００４０】
また、シミュレーション対象の弾性振動体については、実際に実験を行うことによって、空力減衰率δaを表す振幅ｙ0の関数（δa＝ｆI（ｙ0））及び応答円振動数ｎを表す振幅ｙ0の関数（ｎ＝ｆR（ｙ0））を、シミュレーションを正確に求めることができる。そして、このような関数があれば、時々刻々と変化する空力減衰率や応答円振動数を、迅速、簡単かつ確実に求めることができることとなり、シミュレーションをより迅速、簡単かつ確実に実行することができる。
【００４１】
また、前記空力減衰率δaを表す振幅ｙ0の関数（δa＝ｆI（ｙ0））及び前記応答円振動数ｎを表す振幅ｙ0の関数（ｎ＝ｆR（ｙ0））を風速ごとに求めたものを用いることにより、正確なシミュレーションを行うことができる。
【発明を実施するための最良の形態】
【００４２】
本願の発明の時刻歴応答シミュレーション方法は、外力として非定常空気力を受ける弾性振動体の時刻歴応答シミュレーション方法であって、非定常空気力を受ける前記弾性振動体の所定部位の時刻ｔにおける変位ｙt及び変位速度ｙ’tの値を用いて前記所定部位の振幅ｙ0の時刻ｔにおける値を計算する第１ステップと、振幅ｙ0の時刻ｔにおける値を用いて、振幅ｙ0の関数である空力減衰率δa及び応答円振動数ｎの時刻ｔにおける値を計算する第２ステップと、空力減衰率δaの時刻ｔにおける値を用いて前記非定常空気力の振動変位速度成分ｆaIの時刻ｔにおける値を計算すると共に、応答円振動数ｎの時刻ｔにおける値を用いて前記非定常空気力の振動変位成分ｆaRの時刻ｔにおける値を計算する第３ステップと、振動変位成分ｆaR及び振動変位速度成分ｆaIを変数として含む弾性振動体の振動方程式に振動変位速度成分ｆaI及び振動変位成分ｆaRの時刻ｔにおける値を代入して得られた振動方程式を解いて、次の時刻ｔ+Δｔにおける前記所定部位の変位ｙt+Δt、変位速度ｙ’t+Δt及び加速度ｙ’’t+Δtの値を計算する第４ステップと、を備え、これらのステップを含む時刻歴応答計算ルーチンを各時刻について繰り返し実行して計算された値を基に弾性振動体の時刻歴応答シミュレーションを行うものであり、２回目以降の時刻歴応答計算ルーチンの第１ステップでは、当該第１ステップにおける変位ｙt及び変位速度ｙ’tの値として、直前の時刻歴応答計算ルーチンの第４ステップで計算された変位ｙ及び変位速度ｙ’の値を用いて振幅ｙ0の対応時刻における値を計算するものである。
【００４３】
そして、２回目以降の時刻歴応答計算ルーチンの第４ステップは、当該第４ステップにおける変位ｙt、変位速度ｙ’t及び加速度ｙ’’tの値として、直前の時刻歴応答計算ルーチンの第４ステップで計算された変位ｙ、変位速度ｙ’及び加速度ｙ’’の値を用いて、次の時刻ｔ+Δｔにおける前記所定部位の変位ｙt+Δt、変位速度ｙ’t+Δt及び加速度ｙ’’t+Δtの値を計算するものである。
【００４４】
また、前記所定部位の振幅ｙ0の時刻ｔにおける値を求めるための前記変位ｙ及び変位速度ｙ’の関数は、時刻ｔにおける変位ｙtの自乗値と、時刻ｔにおける変位速度ｙ’tを応答円振動数ｎで除した値の自乗値との和の平方根を振幅ｙ0の時刻ｔにおける値として求めることができるものである。
【００４５】
また、前記振動方程式は、下記の式（６）であり、前記第４ステップで振動変位速度成分ｆaI及び振動変位成分ｆaRの値が代入された振動方程式の解法は直接時間積分法であり、前記第１ステップで前記振幅ｙ0の値を計算するために用いられる変位ｙ及び変位速度ｙ’の関数は、下記の式（１０）であり、前記第３ステップで前記振動変位速度成分ｆaIの値を計算するために用いられる前記空力減衰率δaの関数は、下記の式（１６）であり、前記第３ステップで前記振動変位成分ｆaRの値を計算するために用いられる前記応答円振動数ｎの関数は、下記の式（１８）であることが好ましい。
【数９】
【００４６】
また、前記空力減衰率δaを求めるための前記振幅ｙ0の関数（δa＝ｆI）及び前記応答円振動数ｎを求めるための前記振幅ｙ0の関数（ｎ＝ｆR）は、前記弾性振動体について行う振動実験によって求められるものであってもよい。さらに、前記関数（δa＝ｆI）及び前記関数（ｎ＝ｆR）は、前記振幅ｙ0及び風速の関数であってもよい。
【００４７】
また、本願に係る別の発明は、外力として非定常空気力を受ける弾性振動体の時刻歴応答シミュレーションをするためのコンピュータを、シミュレーションの対象である前記弾性振動体の所定部位の条件を設定する手段と、風速条件を設定する手段、前記弾性振動体の所定部位の初期時刻における変位ｙ及び変位速度ｙ’の値を設定する手段、前記弾性振動体の所定部位の時刻ｔにおける変位ｙt及び変位速度ｙ’tの値を用いて当該所定部位の振幅ｙ0の時刻ｔにおける値を計算する手段、前記風速条件に基づいて特定された時刻ｔにおける風速と、振幅ｙ0の時刻ｔにおける値に基づいて空力減衰率δa及び応答円振動数ｎの時刻ｔにおける値を計算する手段、空力減衰率δaの時刻ｔにおける値を用いて前記非定常空気力の振動変位速度成分ｆaIの時刻ｔにおける値を計算すると共に、応答円振動数ｎの時刻ｔにおける値を用いて前記非定常空気力の振動変位成分ｆaRの時刻ｔにおける値を計算する手段、振動変位成分ｆaR及び振動変位速度成分ｆaIを変数として含む弾性振動体の振動方程式に、振動変位速度成分ｆaI及び振動変位成分ｆaRの時刻ｔにおける値を代入して得られた振動方程式を解いて、次の時刻ｔ+Δｔにおける前記所定部位の変位ｙt+Δt、変位速度ｙ’t+Δt及び加速度ｙ’’t+Δtの値を計算する手段、前記計算する手段を順次実行する時刻歴応答計算ルーチンを各時刻について繰り返し実行させる手段であって、２回目以降の時刻歴応答計算ルーチンの振幅ｙ0の値の計算で用いられる変位ｙt及び変位速度ｙ’tの値として、直前の時刻歴応答計算ルーチンで計算された変位ｙ及び変位速度ｙ’の値を用いて振幅ｙ0の値を計算させる手段、前記時刻歴応答計算ルーチンを各時刻について繰り返し実行して計算された値を用いて弾性振動体の時刻歴応答シミュレーションを出力する手段、として機能させる時刻歴応答シミュレーションプログラムである。
【００４８】
また、本願に係るさらに別の発明は、外力として非定常空気力を受ける弾性振動体の時刻歴応答シミュレーションを行う時刻歴応答シミュレーション装置であって、非定常空気力を受ける前記弾性振動体の所定部位の時刻ｔにおける変位ｙt及び変位速度ｙ’tの値を用いて前記所定部位の振幅ｙ0の時刻ｔにおける値を計算する第１計算手段と、振幅ｙ0の時刻ｔにおける値を用いて、振幅ｙ0の関数である空力減衰率δa及び応答円振動数ｎの時刻ｔにおける値を計算する第２計算手段と、空力減衰率δaの時刻ｔにおける値を用いて前記非定常空気力の振動変位速度成分ｆaIの時刻ｔにおける値を計算すると共に、応答円振動数ｎの時刻ｔにおける値を用いて前記非定常空気力の振動変位成分ｆaRの時刻ｔにおける値を計算する第３計算手段と、振動変位成分ｆaR及び振動変位速度成分ｆaIを変数として含む弾性振動体の振動方程式に、振動変位速度成分ｆaI及び振動変位成分ｆaRの時刻ｔにおける値を代入して得られた振動方程式を解いて、次の時刻ｔ+Δｔにおける前記所定部位の変位ｙt+Δt、変位速度ｙ’t+Δt及び加速度ｙ’’t+Δtの値を計算する第４計算手段と、前記計算する手段を順次実行する時刻歴応答計算ルーチンを各時刻について繰り返し実行させる手段であって、２回目以降の時刻歴応答計算ルーチンの振幅ｙ0の値の計算で用いられる変位ｙt及び変位速度ｙ’tの値として、直前の時刻歴応答計算ルーチンで計算された変位ｙ及び変位速度ｙ’の値を用いて振幅ｙ0の値を計算させる手段と、前記時刻歴応答計算ルーチンを各時刻について繰り返し実行して計算された値を用いて弾性振動体の時刻歴応答シミュレーションを出力する手段と、を備える時刻歴応答シミュレーション装置である。
【実施例】
【００４９】
以下、本発明に係る時刻歴応答シミュレーション方法、そのプログラム及び装置の実施例について、図面を参照しつつ詳細に説明する。ここでは、図１及び図２に示される箱桁断面の橋梁模型（空力弾性振動体）を例に、バネ吊り１自由度振動系の空力弾性応答のシミュレーションについて説明する。
【００５０】
まず、シミュレーションの前にあらかじめ、図１に示される形状及び寸法（ｍｍ）の橋梁模型の桁部（所定部位）について自由振動実験を行った。より具体的には、橋梁模型の桁部の一部に該当する、図２に示される二次元部分剛体模型をバネで吊った状態で風洞内に設置し、風を吹かせて応答を測定するという自由振動実験を行った（図３参照）。自由振動実験に用いた二次元部分剛体模型は、実物の１／６５の縮尺のものであり、長さが８００ｍｍ、質量ｍが１．５６ｋｇ、構造減衰率δsが０．００１９２２９、固有振動数ωが７．７５６Ｈｚというものであった。自由振動実験では、風速を徐々に変化させ、各風速における発散振動状況、減衰振動状況の振動波形をそのつど記録し、記録した振動波形から各風速における振幅ｙ0ごとの減衰率δを算出した。ここでは風速０．２ｍ／ｓごとに上記の記録を行った。そして、算出した減衰率δから構造減衰率δsを差し引いて、振幅ｙ0の関数である空力減衰率δa（＝ｆI（ｙ0））を求めた。また、記録した振動波形を基に、振幅ｙ0の関数である後述の応答円振動数ｎ（＝ｆR（ｙ0））を求めた。所定部位とは、橋梁の桁部、主塔部など、それぞれ固有の断面形状を有する全ての部分のことである。ここでは、連続高架橋の桁部という部位の時刻歴応答シミュレーションについて説明する。なお、必要に応じて各部位について自由振動実験を行い、各部位の空力減衰率δa、応答円振動数ｎを算出することになるが、その手順はここで説明するものと同様であるので、説明を省略する。
【００５１】
ここで、実験によって求めたデータを基に作成した１自由度たわみ応答（二次元模型の対風振動応答図）を示す（図４の実験値のグラフ参照）。この図では、横軸が換算風速Ｖr、縦軸が無次元化した振幅（２Ａ／Ｄ）である。
また、自由振動実験を行って求めた空力減衰率δaのグラフの一例として、渦励振時のピークを示す換算風速Ｖr＝５．７６における空力減衰率δaの応答図を示す（図５の実験値のグラフ参照）。そして、関数近似を行って、次の式（２１）で示されるような空力減衰率δaの近似式（近似関数）を得た。
下記の式から解かるように、ここでは、振幅値について区間を複数設定し（実施例は３区間）、各区間ごとに近似式を求めた。このようにして、振幅値の値に応じた近似式を用いて空力減衰率δaを求めるようにすると、より正確に空力減衰率δaを求めることができ、より正確なシミュレーションを行うことができる。また、この点は、後述の応答円振動数についても同様である。なお、他の換算風速においても同様に、必要に応じて区間を区切って空力減衰率δaの近似式（近似関数）を得たが、手順等は同様であるので、ここではその詳細な説明を省略する。
なお、近似化に用いた数値を直接用いてシミュレーションを行ってもよいが、シミュレーションのためには、関数化しておく方が都合が良く、好ましい。これは、空力減衰率についてだけでなく、後述の応答円振動数ｎについても同様である。
【数１０】
【００５２】
また、空力減衰率δaを求める際に、同時に各風速における応答円振動数ｎを求めた。一例として換算風速Ｖr＝５．７６に対する応答円振動数ｎの応答図を示す（図６の実験値のグラフ参照）。そして、関数近似を行って、次の式（２２）で示されるような応答円振動数ｎの近似式を得た。下記の式から解かるように、空力減衰率δaの場合と同様、振幅値について区間を複数設定し（ここでは２区間）、各区間ごとに近似式を求めた。この近似式を用いて、いわゆる非定常空気力の変位比例成分を得ることができる。なお、他の換算風速においても同様に、必要に応じて区間を区切って応答円振動数ｎの近似式を得たが、手順等は同様であるので、ここではその詳細な説明を省略する。
【数１１】
【００５３】
このようにして、各風速および振幅に対応した橋梁の所定部位の空力減衰率δaを表す近似式（近似関数）及び応答円振動数ｎを表す近似式（近似関数）を求めた上で、次に説明しいている各ステップから構成される時刻歴応答計算ルーチンを時々刻々と変化する各時刻について繰り返し実行し、時刻歴応答シミュレーションを行った。
【００５４】
以下、時刻歴応答シミュレーション方法、そのプログラム及び装置について説明する。
【００５５】
図９に示されるように、本シミュレーション方法を実行するシミュレーション装置10は、ＣＰＵ11等の演算プロセッサやメインメモリ12等のメモリを備え、ＯＳが例えばＷＩＮＤＯＷＳ（登録商標）ＸＰエディションであるパソコン本体13と、パソコン本体13に入出力インターフェース14を介して接続された出力手段であるディスプレイ15と、プリンタ（不図示）と、外部記憶装置であるハードディスク16と、ＣＤ−ＲＯＭの記録を読み取るためのＣＤ−ＲＯＭドライブ17（または各種記憶媒体のためのマルチドライブ）と、パソコン本体13に入出力インターフェース14を介して接続された入力手段であるマウス18やキーボード19とを備える周知のコンピュータである。なお、符号「20」はバスである（図９参照）。
そして、本発明にかかる時刻歴応答シミュレーションプログラムがハードディスク15又はマルチドライブにより読み取り可能なＣＤ−ＲＯＭ等の記憶媒体に記憶されており、本シミュレーション装置10にインストールされている。そして、主にＣＰＵ11やメインメモリ12が計算などのステップ処理部や出力するデータを生成する手段として機能し、メインメモリ12が入力された設定条件の記憶手段として機能する。また、計算によって算出された数値データや生成された画像データ等がハードディスク16やディスプレイ15などに出力される。なお、シミュレーション装置10として用いられているコンピュータ自体は周知の装置であるので、ここでは、その詳細な説明を省略する。
【００５６】
次に、時刻歴応答シミュレーション方法について説明する。本シミュレーション方法は、上述したように、時刻ｔにおける所定部位の振幅ｙ0tに対応する空気力が所定部位に瞬間瞬間作用しているという考え方を用いて空力弾性振動の応答を追跡し、シミュレーションするものである。
なお、時刻ｔにおける所定部位の振幅ｙ0とは、各時刻ｔにおける所定部位の実際の変位ｙ（図７の実線参照）ではなく、時刻ｔにおいて取り得た最大変位（正弦波振幅、図７の破線参照）のことである。別言すれば、時刻ｔにおける所定部位の時々刻々変化する振幅のことである。この振幅ｙ0は、図７の破線で示されるように変化する。
【００５７】
まず、シミュレーションの対象である弾性振動体の所定部位の条件を設定する。本実施例では、所定部位として桁部（図１参照）が設定され、桁部の条件として上述したような質量等の条件が設定される。そして、シミュレーションしようとしている状況に応じた風速条件（換算風速Ｖrの条件）を設定する。外力として非定常空気力を受ける状態を検討する場合は、時々刻々と変化する条件で設定がなされる。各時刻における具体的風速値（および／または換算風速Ｖr）をあらかじめ設定しておき必要に応じて読み取って用いることができるようにしてもよいし、各時刻における風速を算出するための関数を用意しておき必要に応じて各時刻における風速を算出して用いることができるようにしてもよい。これにより、シミュレーションの対象である橋梁の所定部位の所定時刻における風速（換算風速Ｖr）を定めることができる。また、最初の時刻における所定部位の変位、変位速度及び変位加速度を設定する。なお、各種条件設定については、キーボード等の入力手段を用いた入力や、あらかじめ記憶させておいた種々の条件の中からの選択など、種々の方法によって行うことができる。各種設定条件は、必要に応じて記憶手段に記憶されて用いられる。
【００５８】
そして、インストールされたプログラムに基づいて、以降に説明するような計算を行う。
まず、所定時刻ｔ1（最初の時刻歴応答計算ルーチンでは最初の時刻）における風速を特定する。具体的には、先に設定した風速条件に基づいて、風速または換算風速Ｖrの具体的値が特定される。また、時刻ｔ1の所定部位の変位ｙt1及び変位速度ｙ’t1を式（１０）中の変位ｙ及び変位速度ｙ’に代入して振幅ｙ0の計算する。これにより、非定常空気力（外力）を受ける橋梁の桁部の振幅ｙ0の時刻ｔ1における値が求められる（第１ステップ）。時刻歴応答シミュレーションプログラムは、このような計算をする計算ルーチンを備えている（第１計算手段）。
【数１２】
【００５９】
なお、変位速度が０のときは、所定部位の実際の変位は振幅ｙ0に一致するので、例えば時刻ｔ1における変位速度が０である場合は、式（１０）を用いずに、実際の変位ｙt1を時刻ｔ1の振幅ｙ0としてもよい。
【００６０】
次に、空力減衰率δa及び応答円振動数ｎの値を計算する（第２ステップ）。まず、第１ステップで定められた風速条件及び振幅ｙ0の値に応じて、ここでの計算で用いる空力減衰率δaの近似式および応答円振動数ｎの近似式を特定する。例えば、風速が「換算風速Ｖr＝５．７６」に対する風速であり、振幅ｙ0が０．０３であれば、空力減衰率δaの計算に用いられる近似式は、上記式（２１）中の近似式２に特定され、応答円振動数ｎの計算に用いられる近似式は、上記の式（２２）中の近似式１に特定される。なお、完全に一致する条件に対応した近似式がない場合は、最も近い条件に対応した近似式を、用いる近似式として特定してもよい。このようにして特定された式に、時刻ｔ1における振幅ｙ0t1の値を代入して、当該時刻ｔ1における所定部位の空力減衰率δa及び応答円振動数ｎの値を計算する（第２ステップ）。時刻歴応答シミュレーションプログラムは、このような計算をする計算ルーチンを備えている（第２計算手段）。
【００６１】
続いて、時刻ｔ1の時点で所定部位に加わっている外力を求める。具体的には、第２ステップで求めた空力減衰率δaの具体的数値を式（１６）に代入して振動変位速度成分ｆaIの具体的数値を計算し、第２ステップで求めた応答円振動数ｎの具体的数値を式（１８）に代入して振動変位成分ｆaRの具体的数値を計算する（第３ステップ）。時刻歴応答シミュレーションプログラムは、このような計算をする計算ルーチンを備えている（第３計算手段）。
【数１３】
【００６２】
そして、計算された振動変位速度成分ｆaI及び振動変位成分ｆaRの具体的数値を式（６）に代入することによって、時刻ｔ1における振動方程式が特定される。
【００６３】
この特定された振動方程式を、直接時間積分法で解いて、時間刻みに応じた振幅の時刻歴を計算する（第４ステップ）。時刻歴応答シミュレーションプログラムは、このような計算をする計算ルーチンを備えている（第４計算手段）。
【００６４】
本実施例では、不規則外力を受ける構造物の応答計算を数値計算によって行う際に用いられる直接時間積分法の１つであるニューマークβ法を用いて、振幅の時刻歴を計算する。
ニューマークβ法は、不規則外力を受ける構造物の応答計算を行う際に用いられる数値計算方法の１つである。この方法を用い、振動方程式の外力項を空気力とすることで、時系列での応答推定を行う。
ニューマークβ法は、時刻ｔ+Δｔ（＝例えば時刻ｔ2）のときの所定部位の実際の変位を求めるために、時刻ｔ（例えば時刻ｔ1）のときの値とΔｔ時間後である時刻ｔ+Δｔ（例えばｔ2）のときの値を使用して繰り返し計算を行う解法である。使用する方程式は、次の式（２４）である。
【数１４】
【００６５】
ここで、上記式（２４）を式（６）の形に置き換えると、次の式（２５）のようになり、さらにマトリクス表示すると、次の式（２６）のようになる。
【数１５】
【００６６】
つまり、第４ステップでは、上記式（２６）を用いて、時刻ｔ2（＝ｔ+Δｔ）における前記所定部位の変位ｙt+Δt、変位速度ｙ’t+Δt及び変位加速度ｙ’’t+Δtを計算することができる（第４ステップ、第４計算手段）。
これにより、最初の時刻歴応答計算ルーチンＲ1において、所定部位の時刻ｔ2（＝ｔ+Δｔ）における実際に変位が求められたことになる。
【００６７】
また、シミュレーションプログラムは、計算された変位ｙt2、変位速度ｙ’t2及び変位加速度ｙ’’t2の値を、次の時刻t3における所定部位の変位等を求める次の時刻歴応答計算ルーチンＲ2において、変位ｙt及び変位速度ｙ’tの値として用いて計算を行わせるようになっている。例えば、計算された時刻ｔ2（＝ｔ+Δｔ）の変位ｙの値及び変位速度ｙ’の値は、時刻歴応答計算ルーチンＲ2の第１ステップにおいて、式（１０）の変位ｙ及び変位速度ｙ’の値として代入され、これにより時刻ｔ3（＝ｔ+２×Δｔ）における所定部位の振幅ｙ0が計算される（第１ステップ）。また、計算された変位ｙt2、変位速度ｙ’t2及び変位加速度ｙ’’t2の値は、時刻歴応答計算ルーチンＲ2の第４ステップにおいて、式（２６）の変位ｙt、変位速度ｙ’t及び変位加速度ｙ’’tの値として代入され、これにより、時刻t3における所定部位の変位ｙt+Δt、変位速度ｙ’t+Δt及び変位加速度ｙ’’t+Δtが計算される（第４ステップ）。
【００６８】
このように、２回目以降のｘ回目の時刻歴応答計算ルーチンＲxの第１ステップでは、直前の時刻歴応答計算ルーチンＲx-1の第４ステップで計算された変位ｙtx及び変位速度ｙ’txの値を、変位ｙ及び変位速度ｙ’の値として用いて（式（１０）等参照）、次の時刻ｔx+1における振幅ｙ0を計算する。
また、２回目以降の時刻歴応答計算ルーチンＲxの第４ステップでは、直前の時刻歴応答計算ルーチンＲx-1の第４ステップで計算された変位ｙtx、変位速度ｙ’tx及び変位加速度ｙ’’txの値を、当該第４ステップにおける変位ｙt、変位速度ｙ’t及び変位加速度ｙ’’tの値として用いて次の時刻ｔx+1における変位ｙ、変位速度ｙ’及び変位加速度ｙ’’の値を計算する。
時刻歴応答シミュレーションプログラムは、このように時刻歴応答計算ルーチンを各時刻について繰り返し実行させる計算ルーチンであって、２回目以降の時刻歴応答計算ルーチンの各計算で用いられる変位ｙt、変位速度ｙ’t及び変位加速度ｙ’’の値として、直前の時刻歴応答計算ルーチンで計算された変位ｙ、変位速度ｙ’及び変位加速度ｙ’’の値を用いて振幅ｙ0の値を計算させる計算ルーチンを、計算手段として備えている。
【００６９】
このようなステップからなる時刻歴応答計算ルーチンを各時刻について繰り返し実行することによって、各時刻における応答（実際の変位等）を計算し、計算された変位等の値を用いて、Δｔ時間刻みの時刻歴応答シミュレーションを行う。シミュレーションの結果は、ディスプレイ15やプリンタなどに出力される。ハードディスク16などの記憶媒体に記録される態様で出力することもできる。時刻歴応答シミュレーションプログラムは、このような機器に時刻歴応答シミュレーションの結果を出力させる処理ルーチンを備えている。
【００７０】
そして、上記のようなシミュレーションを行った結果、図８（ａ）に示されるような時刻歴応答解析波形が得られた。また、比較のために、実物の１／１３０の縮尺の全橋模型について、シミュレーションの条件と同じ条件で風洞実験を行って、時刻歴応答解析波形を得た。図８（ｂ）に、風洞実験での振幅の測定値に基づいて作成した振幅波形を示す。両波形を比較するとわかるように、両波形はほぼ一致していた。
【００７１】
そして、図５中の推定値のグラフは、図８（ａ）に示される時刻歴応答解析波形を基に算出した空力減衰率の推定値のグラフである。図５から明らかなように、空力減衰率の近似値及び推定値は、いずれも、実験値とほぼ一致していた。
【００７２】
また、時刻歴解析を各風速で解析し、定常振動応答をプロットして、測定値との比較を行ったものが図４である。この図に示されるように、時刻歴応答解析から求めた定常振動応答推定値は、若干のずれはあるものの実験での測定値をほぼ再現していた。
【００７３】
このように、本願発明によれば、減衰状態あるいは発散状態にある所定部位の空力弾性振動について、時々刻々の振幅ｙ0を計算することができ、計算された振幅ｙ0を用いて非定常空気力を計算することができる。時々刻々の非定常空気力を計算することができれば、計算された非定常空気力を振動方程式に代入することによって時々刻々の振動方程式を定めて、定められた振動方程式を解くことによって、空力弾性振動の時刻歴応答解析シミュレーションをすることができる。
【００７４】
本願の発明は、以上の実施例に限定されるものではなく、その要旨を逸脱しない範囲において、種々の変形が可能である。
【００７５】
上記実施例では、箱桁断面の橋梁模型（空力弾性振動体）を例にシミュレーションについて説明したが、本願の発明にかかるシミュレーション方法及び装置は、例えば、橋梁の桁部、主塔部など、弾性振動する種々の部位についての時刻歴応答シミュレーションに用いることができる。
【００７６】
また、上記実施例では、箱桁断面の橋梁模型を例に、バネ吊り１自由度振動系の空力弾性応答のシミュレーションについて説明したが、本願の発明にかかるシミュレーション方法及び装置は、その他にも、バネ吊り２自由度振動系の空力弾性応答や全橋の空力弾性振動応答、乱流中での橋桁の空力弾性振動応答のシミュレーションなど、種々の振動に関する空力弾性応答のシミュレーションに用いることができる。
【００７７】
また、上記実施例の自由振動実験では、風速を徐々に変化させ、０．２ｍ／ｓごとに発散振動状況、減衰振動状況の振動波形を記録し、記録した振動波形から各風速における振幅ごとの減衰率を算出したが、波形を測定する風速の間隔は必ずしも０．２ｍ／ｓごとでなければならないわけではなく、これにより大きくても良いし、小さくても良い。
【図面の簡単な説明】
【００７８】
【図１】本発明のシミュレーション方法のシミュレーション対象である連続高架橋を示すものであり、橋脚上に配置された桁部の概略形状を示す斜視図である。
【図２】図１に示される橋梁の桁部の断面を示す断面図である。
【図３】本実施例の自由振動実験で用いた装置の概略構成を示す構成図である。
【図４】１自由度たわみ応答図（二次元模型の対風振動応答図）である。
【図５】換算風速Ｖr＝５．７６における空力減衰率δaの応答図である。
【図６】換算風速Ｖr＝５．７６に対する応答円振動数ｎの応答図である。
【図７】非定常空気力を受けた弾性振動体の過渡応答時の振動振幅を示すグラフである。
【図８】（ａ）は、換算風速Ｖr＝５．７６における時刻歴応答解析波形を示す図であり、（ｂ）は、実際の振幅の測定値に基づいて作成した振幅波形を示す図である。
【図９】本発明のシミュレーション装置の構成を示す構成図である。
【技術分野】
【０００１】
本発明は、時刻歴応答解析に関するものであり、特に空力弾性振動の時刻歴応答解析シミュレーション方法、そのプログラム及び装置に関するものである。
【背景技術】
【０００２】
時刻歴応答解析は、耐震工学の分野では古い歴史がある。耐震工学の分野における時刻歴応答解析とは、地震振動（地震加速度）などの外力が与えられた構造物の所定位置の応答加速度、速度、変位などを計算することである。そして、計算された解析結果は、構造物の耐震性の判断に用いられている。
【０００３】
ところで、橋梁の耐風工学の分野では、橋梁、特に長大橋梁の建設の際、耐風特性の検討が行われる（特許文献１参照）。耐風特性の検討では、例えば、風に代表される外力（空気力）を受ける橋梁の所定位置の応答（変位）が検討される。そして、検討結果が橋梁の耐風設計に用いられる。このような耐風特性の検討のために、時刻歴応答解析を用いることが考えられるが、現実には時刻歴応答解析を提案する流れはなかった。
【０００４】
時刻歴応答解析を行うためには、構造物に作用する外力の特性を十分に把握して評価できる必要があるところ、橋梁等の構造物が受ける外力である空気力の特性を把握することは困難であるということが原因であると考えられる。つまり、空力弾性振動のような自励振動における空気力の評価は、振動振幅に対して行わざるを得ず、過渡応答的な振動における各瞬間の空気力を評価することは困難であるため、耐風特性の検討に時刻歴応答解析を用いることは困難であると考えられた。
【０００５】
このようなことから、現在、橋梁全体の挙動は、一般に全橋模型を用いた風洞実験によって把握されており、橋梁の耐風特性もこの全橋模型を用いた風洞実験によって検討されている。
【特許文献１】特開平８−１８４５２５号公報
【発明の開示】
【発明が解決しようとする課題】
【０００６】
ところが、風洞実験の場合、構造物の構造特性を種々変更して実験を行うには限界がある。そして、風洞実験では、乱流中での構造物の応答挙動を推定することが非常に困難であり、耐風特性の検討に関しても統計的な検討にならざるを得ないという問題がある。また、風洞実験に使用できる空間は限られた空間であるので、風洞実験では、実物の橋梁に対して大幅に縮小された全橋模型を使用せざるを得ず、実験精度が低下しやすいという問題もある。
【０００７】
なお、対風応答に関するシミュレーションとしては、有限要素法（ＦＥＭ）を用いたフラッタ解析がある。フラッタとは、元来、航空機の翼などの構造物に作用する外力（空力）、弾性力及び慣性力が組み合わさって生じる自励振動のことである。フラッタ解析によれば、このようなフラッタ（自励振動）の発生の有無を判断することはできるが、外力（例えば、空気力、減衰力、弾性力）を受けた構造物のフラッタ（自励振動）の振幅（応答振幅）を追跡することはできない。
【０００８】
本発明は、上記問題点に鑑みてなされたものであり、空力弾性振動の時刻歴応答解析シミュレーション方法およびその装置を提供することを課題とする。
【課題を解決するための手段】
【０００９】
本願発明の発明者は、上記課題に鑑み、減衰状態や発散状態における構造物の時々刻々の振幅に基づいて非定常空気力（瞬時作用空気力）を求める手法について、後述するとおり鋭意研究を重ねた結果、次のような本願発明に想到するに至った。
【００１０】
つまり、本願発明は、外力として非定常空気力を受ける弾性振動体の時刻歴応答シミュレーション方法であって、非定常空気力を受ける前記弾性振動体の所定部位の時刻ｔにおける変位ｙt及び変位速度ｙ’tの値を用いて前記所定部位の振幅ｙ0の時刻ｔにおける値を計算する第１ステップと、振幅ｙ0の時刻ｔにおける値を用いて、振幅ｙ0の関数である空力減衰率δa及び応答円振動数ｎの時刻ｔにおける値を計算する第２ステップと、空力減衰率δaの時刻ｔにおける値を用いて前記非定常空気力の振動変位速度成分ｆaIの時刻ｔにおける値を計算すると共に、応答円振動数ｎの時刻ｔにおける値を用いて前記非定常空気力の振動変位成分ｆaRの時刻ｔにおける値を計算する第３ステップと、振動変位成分ｆaR及び振動変位速度成分ｆaIを変数として含む弾性振動体の振動方程式に振動変位速度成分ｆaI及び振動変位成分ｆaRの時刻ｔにおける値を代入して得られた振動方程式を解いて、次の時刻ｔ+Δｔにおける前記所定部位の変位ｙt+Δt、変位速度ｙ’t+Δt及び加速度ｙ’’t+Δtの値を計算する第４ステップと、を備え、これらのステップを含む時刻歴応答計算ルーチンを各時刻について繰り返し実行して計算された値を基に弾性振動体の時刻歴応答シミュレーションを行うものであり、２回目以降の時刻歴応答計算ルーチンの第１ステップでは、当該第１ステップにおける変位ｙt及び変位速度ｙ’tの値として、直前の時刻歴応答計算ルーチンの第４ステップで計算された変位ｙ及び変位速度ｙ’の値を用いて振幅ｙ0の対応時刻における値を計算するものである。
【００１１】
本発明のシミュレーション方法は、基本的には、次の式（１）に示される振動方程式に係るものである。そして、式（１）は、式（２）、（３）を適用することにより、式（４）のように表される。
【数２】
【００１２】
ここで、非定常空気力の振動変位速度成分ｆaI及び非定常空気力の振動変位成分ｆaRを、次の式（５）のようにおく。そうすると、上記式（１）の振動方程式を、下記の式（６）のように表すことができる。このように、振動変位速度成分ｆaI及び非定常空気力の振動変位成分ｆaRを式（５）のようにおいて振動方程式を式（６）のように表す考え方は好ましいものである。
【数３】
【００１３】
ところで、式（６）では、非定常空気力の振動変位速度成分ｆaIが減衰項の一部を構成し、非定常空気力の振動変位成分ｆaRが復元力項の一部を構成している。この振動方程式を解くためには、振動変位速度成分ｆaI及び振動変位成分ｆaRの具体的な数値を把握する必要がある。
【００１４】
そこで、非定常空気力の振動変位速度成分ｆaI及び振動変位成分ｆaRの具体的な数値を把握する方法について検討した。その結果、空力弾性振動体の所定部位の空力減衰率δaが解かれば振動変位速度成分ｆaIの具体的数値を把握することができ、空力弾性振動体の所定部位の応答円振動数ｎが解かれば振動変位成分ｆaRの具体的数値を把握することができることが解かった。ただし、空力弾性振動体に加わる外力は時々刻々と変化する非定常空気力であるので、簡単に所定部位の空力減衰率δa及び応答円振動数ｎを把握することはできない。
【００１５】
そこで、空力減衰率δa及び応答円振動数ｎを把握する方法について、さらに検討し、その結果、振動体の所定部位に作用する各時刻の非定常空気力は、時々刻々と変化する所定部位の振幅ｙ0（図７参照）に対応しているという考え方を用いることによって、空力弾性振動の応答を追跡できることを見出した。より具体的に言えば、空力減衰率δa及び応答円振動数ｎの具体的数値を、時々刻々と変化する所定部位の振幅ｙ0に基づいて把握することができることを見出すと共に、各時刻における振幅ｙ0を、当該時刻における変位及び変位速度に基づいて求めることができることを見出した。これを実際のシミュレーション方法の処理の流れに沿って説明すれば、各時刻における変位及び変位速度を用いて当該時刻の振幅ｙ0を求め、当該振幅ｙ0から空力減衰率δa及び応答円振動数ｎの具体的数値を把握し、これらを用いて非定常空気力の振動変位速度成分ｆaI及び振動変位成分ｆaRの具体的数値を把握し、把握した非定常空気力に関する具体的数値を用いてシミュレーションに用いる振動方程式を確定することができることを見出したということになる。
なお、時刻ｔにおける所定部位の振幅ｙ0とは、各時刻ｔにおける変位ｙtではなく、その時刻ｔにおいて取り得た最大変位のことである。振幅ｙ0は、図７の破線で示されるように変化する。
これらのことを見出したことで、最終的に、本願発明に係る空力弾性振動体の時刻歴応答シミュレーション方法に想到するに至った。
【００１６】
本発明では、まず、ある時刻の変位及び変位速度に基づいて当該時刻の振幅ｙ0を求めるために、非定常空気力（外力）を受けた空力弾性振動体の所定部位の振動は、正弦波的な振動で減衰しあるいは発散するものであると考える。このように考えると、空力弾性振動の変位ｙを、下記の式（７）のように表すことができ、変位速度ｙ’を、下記の式（８）のように表すことができる。そして、瞬間瞬間の振動に対応する振幅（正弦波振幅）ｙ0を各時刻ｔにおける「変位」の自乗と「変位速度／振動数」の自乗の和の平方根で表現することができるということを見出して、式（７），（８）から式（９）を得た。さらに、橋梁の対風応答を検討する場合、減衰定数ξは０．０１未満であるので、式（９）中のルートの中は、ほぼ「１」と考えることができ、式（９）中の「e」の乗数（-ξωt）はほぼ「０」と考えることができる。そうすると、空力弾性振動の式（９）から式（１０）が得られる。
【数４】
【００１７】
ここで、上記の式（７）から式（１０）に関連して図７を参照する。図７において、実線は、空力弾性振動体の所定部位の時刻ｔにおける変位ｙを示すものである。そして、破線は、上記式（１０）によって計算される時々刻々の振幅ｙ0である。つまり、式（１０）を用いることによって、弾性振動体の所定部位の各時刻における変位ｙt及び変位速度ｙ’tの値から当該所定部位の振幅ｙ0の当該時刻における値を計算することができる（第１ステップ）。
【００１８】
図７から解かるように、空力弾性振動体の所定部位の変位が振幅に一致するとき（変位速度が０のとき）は、式（１０）がなくても、振幅ｙ0を求めることができる。ところが、振幅ｙ0が時々刻々と過渡応答的に変化する空力弾性振動の場合、式（１０）がなければ、変位速度が０のとき以外の時刻について、振幅ｙ0を求めることはできない。
空力弾性振動のような自励振動における空気力（非定常空気力）の評価は、振動振幅に対して行わざるを得ないものであるが、これまで、過渡応答的な振動における各瞬間の振幅ｙ0を評価することができなかったので、振幅ｙ0の変化に対応しつつ時々刻々の非定常空気力を評価することはできなかった。このようなことから、これまでは、時刻歴応答解析を用いて耐風特性を検討することはできず、空力弾性振動の時刻歴応答シミュレーションを行うことができなかった。
【００１９】
この点、本発明によれば、上記式（１０）を用いることによって、空力弾性振動体の所定部位の時刻ｔにおける変位ｙt及び変位速度ｙ’tに基づき、時々刻々と変化する振幅ｙ0を、いずれの時刻においても求めることができる（第１ステップ）。
なお、変位ｙの自乗値と、変位速度ｙ’を応答円振動数ｎで除した値の自乗値との和の平方根の値に基づいて振幅ｙ0を求める方法であれば、種々の方法によって振幅ｙ0を求めることができるが、式（１０）を用いる方法は、簡単であり好ましい。
【００２０】
次に、先に算出した振幅ｙ0の時刻ｔにおける値を用いて、振幅ｙ0の関数である空力減衰率δa及び応答円振動数ｎの時刻ｔにおける値を計算する（第２ステップ）。なお、振幅ｙ0の関数である空力減衰率δa（＝ｆI（ｙ0））及び振幅ｙ0の関数である応答円振動数ｎ（＝ｆR（ｙ0））は、例えば弾性振動体に対して振動実験を行うことによって求められる。
【００２１】
続いて、先に算出された空力減衰率δaの時刻ｔにおける値を用いて前記非定常空気力の振動変位速度成分ｆaIの時刻ｔにおける値を計算すると共に、先に算出された応答円振動数ｎの時刻ｔにおける値を用いて前記非定常空気力の振動変位成分ｆaRの時刻ｔにおける値を計算する（第３ステップ）。
【００２２】
なお、空力減衰率δaの具体的数値を用いて非定常空気力（瞬時作用空気力）の振動変位速度成分ｆaIを計算する方法は、種々考えられるかも知れないが、本発明において振動変位速度成分ｆaIを計算する際に用いる空力減衰率δaの関数（式（１６））は、次にようにして求められる。つまり、式（６）として示される振動方程式中の振動変位速度成分ｆaIが減衰項の一部であり、振動方程式における定常振動の条件が次の式（１１）であり、振幅−空力減衰率図における定常振動の条件が式（１２）であることに基づき、求められる。
【数５】
【００２３】
上記式（１１）、（１２）によって示される２つの関係が共通した現象における条件である。ここで、−ｆaI／ｍ及びδaは空気力に起因するものであり、２ξωがδsに対応（２ξω＝δs）し、また、−ｆaI／ｍがδaに対応する。そして、構造減衰率δs＝２πξであるので、これから導かれる「ξ＝δs／２π」を用いて「２ξω」を変形すると、次式（１３）のようになり、式（１１）及び（１２）の関係から次の式（１４）が導かれる。
【数６】
【００２４】
したがって、式（６）の減衰項は、次の式（１５）のようになり、δaを振動変位速度成分ｆaIの式になるように変形すると、式（１６）のようになる。これにより、二次元剛体モデルとして表される橋梁の振動方程式の振動変位速度成分ｆaIの具体的な形が求められたことになる。また、あらかじめ得た空力減衰率δaを表す近似関数を用いて各時刻における空力減衰率δaを計算することができるので、計算した空力減衰率δaの値を式（１６）に代入することで、その時刻の振幅ｙ0に対応した振動変位速度成分ｆaIが得られる。
【数７】
【００２５】
また、応答円振動数ｎの具体的数値を用いて非定常空気力（瞬時作用空気力）の振動変位成分ｆaRを計算する方法は、種々考えられるかも知れないが、本発明において振動変位成分ｆaRを計算する際に用いられる応答円振動数ｎの関数（式（１８））は、次にようにして求められる。つまり、式（６）として示される振動方程式中の振動変位成分ｆaRが復元力項の一部であり、復元力項の係数が橋梁モデルの応答円振動数の２乗に等しいと考えると、次の式（１７）が得られる。そして、この式を変形することで、次の式（１８）が得られる。また、あらかじめ得た振動数ｎを表す近似関数を用いて各時刻における振動数ｎを計算することができるので、計算した振動数ｎの値を式（１８）に代入することで、その時刻の振動変位成分ｆaRが得られる。
【数８】
【００２６】
このようにして求められた振動変位成分ｆaR及び振動変位速度成分ｆaIの時刻ｔにおける具体的数値を、振動変位成分ｆaR及び振動変位速度成分ｆaIを変数として含む弾性振動体の振動方程式に代入して得られた振動方程式を解いて、次の時刻ｔ+Δｔにおける前記所定部位の変位ｙt+Δt、変位速度ｙ’t+Δt及び加速度ｙ’’t+Δtの値を計算する（第４ステップ）。これにより、時刻歴応答シミュレーションに用いられる数値が計算される。
【００２７】
本願発明に係る空力弾性振動の時刻歴応答解析シミュレーション方法では、これらの第１ステップから第４ステップを含む時刻歴応答計算ルーチンを各時刻について繰り返し実行し、計算された値を基に弾性振動体の時刻歴応答シミュレーションを行うものである。なお、２回目以降の時刻歴応答計算ルーチンの第１ステップでは、当該第１ステップにおける変位ｙt及び変位速度ｙ’tの値として、直前の時刻歴応答計算ルーチンの第４ステップ（次の時刻ｔ+Δｔにおける各値を計算するステップ）で計算された変位ｙ及び変位速度ｙ’の値を用いて振幅ｙ0の対応時刻における値を計算する。また、２回目以降の時刻歴応答計算ルーチンの第４ステップは、当該第４ステップにおける変位ｙt、変位速度ｙ’t及び加速度ｙ’’tの値として、直前の時刻歴応答計算ルーチンの第４ステップ（次の時刻ｔ+Δｔにおける各値を計算するステップ）で計算された変位ｙ、変位速度ｙ’及び加速度ｙ’’の値を用いて、次の時刻ｔ+Δｔにおける前記所定部位の変位ｙt+Δt、変位速度ｙ’t+Δt及び加速度ｙ’’t+Δtの値を計算する。
【００２８】
このように、本願発明によれば、減衰状態あるいは発散状態にある所定部位の空力弾性振動について、時々刻々の振幅ｙ0を計算することができ、計算された振幅ｙ0を用いて非定常空気力（より具体的にはの振動変位成分ｆaR及び振動変位速度成分ｆaI）を計算することができる。時々刻々の非定常空気力を計算することができれば、計算された非定常空気力を振動方程式に代入することによって時々刻々の振動方程式を定め、定められた振動方程式を解くことによって、空力弾性振動の時刻歴応答解析シミュレーションをすることができる。
【００２９】
なお、前記振幅ｙ0の関数である空力減衰率δa（＝ｆI（ｙ0））及び応答円振動数ｎ（＝ｆR（ｙ0））は、前記空力弾性振動体について行う振動実験によって求めたものが好ましい。例えば、シミュレーションの対象の空力弾性振動体について自由振動実験等の振動実験を行うことによって、上記関数（δa＝ｆI（ｙ0））及び関数（ｎ＝ｆR（ｙ0））を求めることができる。
【００３０】
このように、シミュレーション対象の弾性振動体について実際に実験を行うことによって、上記関数である空力減衰率δa（＝ｆI（ｙ0））及び応答円振動数ｎ（＝ｆR（ｙ0））を正確に求めることができる。そして、このような関数があれば、時々刻々と変化する空力減衰率や応答円振動数を、迅速、簡単かつ確実に求めることができることとなり、シミュレーションをより迅速、簡単かつ確実に実行することができる。
【００３１】
また、前記空力減衰率δa（＝ｆI（ｙ0））及び応答円振動数ｎ（＝ｆR（ｙ0））は、前記振幅ｙ0及び風速の関数であることがより好ましい。振幅ｙ0だけでなく風速が加味された関数を用いることにより、より正確なシミュレーションを行うことができる。
【００３２】
また、本願の別の発明は、外力として非定常空気力を受ける弾性振動体の時刻歴応答シミュレーションをするためのコンピュータを、シミュレーションの対象である前記弾性振動体の所定部位の条件を設定する手段と、風速条件を設定する手段、
前記弾性振動体の所定部位の初期時刻における変位ｙ及び変位速度ｙ’の値を設定する手段、前記弾性振動体の所定部位の時刻ｔにおける変位ｙt及び変位速度ｙ’tの値を用いて当該所定部位の振幅ｙ0の時刻ｔにおける値を計算する手段、
前記風速条件に基づいて特定された時刻ｔにおける風速と、振幅ｙ0の時刻ｔにおける値に基づいて空力減衰率δa及び応答円振動数ｎの時刻ｔにおける値を計算する手段、空力減衰率δaの時刻ｔにおける値を用いて前記非定常空気力の振動変位速度成分ｆaIの時刻ｔにおける値を計算すると共に、応答円振動数ｎの時刻ｔにおける値を用いて前記非定常空気力の振動変位成分ｆaRの時刻ｔにおける値を計算する手段、振動変位成分ｆaR及び振動変位速度成分ｆaIを変数として含む弾性振動体の振動方程式に、振動変位速度成分ｆaI及び振動変位成分ｆaRの時刻ｔにおける値を代入して得られた振動方程式を解いて、次の時刻ｔ+Δｔにおける前記所定部位の変位ｙt+Δt、変位速度ｙ’t+Δt及び加速度ｙ’’t+Δtの値を計算する手段、前記計算する手段を順次実行する時刻歴応答計算ルーチンを各時刻について繰り返し実行させる手段であって、２回目以降の時刻歴応答計算ルーチンの振幅ｙ0の値の計算で用いられる変位ｙt及び変位速度ｙ’tの値として、直前の時刻歴応答計算ルーチンで計算された変位ｙ及び変位速度ｙ’の値を用いて振幅ｙ0の値を計算させる手段、前記時刻歴応答計算ルーチンを各時刻について繰り返し実行して計算された値を用いて弾性振動体の時刻歴応答シミュレーションを出力する手段として機能させる時刻歴応答シミュレーションプログラムである。このようなプログラムによって上述したような時刻歴応答シミュレーションを実現することができ、当該プログラムがインストールされたコンピュータなどの時刻歴応答シミュレーション装置によって上述したような時刻歴応答シミュレーションを実現することができる。
【００３３】
そして、本願に係るさらに別の発明は、外力として非定常空気力を受ける弾性振動体の時刻歴応答シミュレーションを行う時刻歴応答シミュレーション装置であって、非定常空気力を受ける前記弾性振動体の所定部位の時刻ｔにおける変位ｙt及び変位速度ｙ’tの値を用いて前記所定部位の振幅ｙ0の時刻ｔにおける値を計算する第１計算手段と、振幅ｙ0の時刻ｔにおける値を用いて、振幅ｙ0の関数である空力減衰率δa及び応答円振動数ｎの時刻ｔにおける値を計算する第２計算手段と、空力減衰率δaの時刻ｔにおける値を用いて前記非定常空気力の振動変位速度成分ｆaIの時刻ｔにおける値を計算すると共に、応答円振動数ｎの時刻ｔにおける値を用いて前記非定常空気力の振動変位成分ｆaRの時刻ｔにおける値を計算する第３計算手段と、振動変位成分ｆaR及び振動変位速度成分ｆaIを変数として含む弾性振動体の振動方程式に、振動変位速度成分ｆaI及び振動変位成分ｆaRの時刻ｔにおける値を代入して得られた振動方程式を解いて、次の時刻ｔ+Δｔにおける前記所定部位の変位ｙt+Δt、変位速度ｙ’t+Δt及び加速度ｙ’’t+Δtの値を計算する第４計算手段と、前記計算する手段を順次実行する時刻歴応答計算ルーチンを各時刻について繰り返し実行させる手段であって、２回目以降の時刻歴応答計算ルーチンの振幅ｙ0の値の計算で用いられる変位ｙt及び変位速度ｙ’tの値として、直前の時刻歴応答計算ルーチンで計算された変位ｙ及び変位速度ｙ’値を用いて振幅ｙ0の値を計算させる手段と、前記時刻歴応答計算ルーチンを各時刻について繰り返し実行して計算された値を用いて弾性振動体の時刻歴応答シミュレーションを出力する手段と、を備えるものである。このような時刻歴応答シミュレーション装置によって上述したような時刻歴応答シミュレーションを実現することができる。
【発明の効果】
【００３４】
本願発明に係る空力弾性振動体の時刻歴応答シミュレーション方法、そのプログラム及び装置によれば、橋梁などの空力弾性振動体の所定部位の空力減衰率δaを表す関数（δa＝ｆ（ｙ0））と、応答円振動数ｎを表す関数（ｎ＝ｆ（ｙ0））とをあらかじめ求めておくことにより、非定常空気力（外力）を受ける当該構造物の耐風特性（例えば、空力弾性振動）の検討を、時刻歴応答解析シミュレーションを用いて行うことができる。
【００３５】
そして、風洞実験を行うことなくシミュレーションによって空力弾性振動の応答推定を行うことができるので、風洞実験を行うことが難しい長大橋梁等の構造物について、空力弾性振動の応答推定を容易、迅速に行うことができる。
【００３６】
また、本願発明によれば、構造物の構造様式などの条件や構造物が受ける風（風速や風速の変動）など、シミュレーションのパラメータの設定を容易に変更することができるので、多様な条件下のシミュレーションを容易、迅速に行うことができる。これにより、多様な条件下での構造物の空力弾性振動の応答推定を、容易、迅速に行うことができる。
【００３７】
さらに、本願発明によれば、時刻歴応答シミュレーションを用いるので、耐震設計と同等の精度で、耐風設計を行うことが可能である。つまり、構造物の空力弾性振動の応答推定を行う際に、本願発明に係る空力弾性振動体の時刻歴応答シミュレーション方法を利用することによって、より正確な耐風設計を行うことが可能である。これにより、構造物の耐風特性に関する設計の品質を飛躍的に向上させることができる。
【００３８】
また、時刻ｔにおける変位ｙtの自乗値と、時刻ｔにおける変位速度ｙ’tを応答円振動数ｎで除した値の自乗値との和の平方根を求めることによって時刻ｔにおける所定部位の振幅ｙ0tの値を求めることができれば、より容易、迅速かつ正確に時刻歴応答シミュレーションを行うことができる。
【００３９】
また、シミュレーションの各ステップでの計算に用いることができる具体的な関数が用意されているので、それらの関数を用いることによって、容易、迅速かつ正確にシミュレーションを行うことができる。
【００４０】
また、シミュレーション対象の弾性振動体については、実際に実験を行うことによって、空力減衰率δaを表す振幅ｙ0の関数（δa＝ｆI（ｙ0））及び応答円振動数ｎを表す振幅ｙ0の関数（ｎ＝ｆR（ｙ0））を、シミュレーションを正確に求めることができる。そして、このような関数があれば、時々刻々と変化する空力減衰率や応答円振動数を、迅速、簡単かつ確実に求めることができることとなり、シミュレーションをより迅速、簡単かつ確実に実行することができる。
【００４１】
また、前記空力減衰率δaを表す振幅ｙ0の関数（δa＝ｆI（ｙ0））及び前記応答円振動数ｎを表す振幅ｙ0の関数（ｎ＝ｆR（ｙ0））を風速ごとに求めたものを用いることにより、正確なシミュレーションを行うことができる。
【発明を実施するための最良の形態】
【００４２】
本願の発明の時刻歴応答シミュレーション方法は、外力として非定常空気力を受ける弾性振動体の時刻歴応答シミュレーション方法であって、非定常空気力を受ける前記弾性振動体の所定部位の時刻ｔにおける変位ｙt及び変位速度ｙ’tの値を用いて前記所定部位の振幅ｙ0の時刻ｔにおける値を計算する第１ステップと、振幅ｙ0の時刻ｔにおける値を用いて、振幅ｙ0の関数である空力減衰率δa及び応答円振動数ｎの時刻ｔにおける値を計算する第２ステップと、空力減衰率δaの時刻ｔにおける値を用いて前記非定常空気力の振動変位速度成分ｆaIの時刻ｔにおける値を計算すると共に、応答円振動数ｎの時刻ｔにおける値を用いて前記非定常空気力の振動変位成分ｆaRの時刻ｔにおける値を計算する第３ステップと、振動変位成分ｆaR及び振動変位速度成分ｆaIを変数として含む弾性振動体の振動方程式に振動変位速度成分ｆaI及び振動変位成分ｆaRの時刻ｔにおける値を代入して得られた振動方程式を解いて、次の時刻ｔ+Δｔにおける前記所定部位の変位ｙt+Δt、変位速度ｙ’t+Δt及び加速度ｙ’’t+Δtの値を計算する第４ステップと、を備え、これらのステップを含む時刻歴応答計算ルーチンを各時刻について繰り返し実行して計算された値を基に弾性振動体の時刻歴応答シミュレーションを行うものであり、２回目以降の時刻歴応答計算ルーチンの第１ステップでは、当該第１ステップにおける変位ｙt及び変位速度ｙ’tの値として、直前の時刻歴応答計算ルーチンの第４ステップで計算された変位ｙ及び変位速度ｙ’の値を用いて振幅ｙ0の対応時刻における値を計算するものである。
【００４３】
そして、２回目以降の時刻歴応答計算ルーチンの第４ステップは、当該第４ステップにおける変位ｙt、変位速度ｙ’t及び加速度ｙ’’tの値として、直前の時刻歴応答計算ルーチンの第４ステップで計算された変位ｙ、変位速度ｙ’及び加速度ｙ’’の値を用いて、次の時刻ｔ+Δｔにおける前記所定部位の変位ｙt+Δt、変位速度ｙ’t+Δt及び加速度ｙ’’t+Δtの値を計算するものである。
【００４４】
また、前記所定部位の振幅ｙ0の時刻ｔにおける値を求めるための前記変位ｙ及び変位速度ｙ’の関数は、時刻ｔにおける変位ｙtの自乗値と、時刻ｔにおける変位速度ｙ’tを応答円振動数ｎで除した値の自乗値との和の平方根を振幅ｙ0の時刻ｔにおける値として求めることができるものである。
【００４５】
また、前記振動方程式は、下記の式（６）であり、前記第４ステップで振動変位速度成分ｆaI及び振動変位成分ｆaRの値が代入された振動方程式の解法は直接時間積分法であり、前記第１ステップで前記振幅ｙ0の値を計算するために用いられる変位ｙ及び変位速度ｙ’の関数は、下記の式（１０）であり、前記第３ステップで前記振動変位速度成分ｆaIの値を計算するために用いられる前記空力減衰率δaの関数は、下記の式（１６）であり、前記第３ステップで前記振動変位成分ｆaRの値を計算するために用いられる前記応答円振動数ｎの関数は、下記の式（１８）であることが好ましい。
【数９】
【００４６】
また、前記空力減衰率δaを求めるための前記振幅ｙ0の関数（δa＝ｆI）及び前記応答円振動数ｎを求めるための前記振幅ｙ0の関数（ｎ＝ｆR）は、前記弾性振動体について行う振動実験によって求められるものであってもよい。さらに、前記関数（δa＝ｆI）及び前記関数（ｎ＝ｆR）は、前記振幅ｙ0及び風速の関数であってもよい。
【００４７】
また、本願に係る別の発明は、外力として非定常空気力を受ける弾性振動体の時刻歴応答シミュレーションをするためのコンピュータを、シミュレーションの対象である前記弾性振動体の所定部位の条件を設定する手段と、風速条件を設定する手段、前記弾性振動体の所定部位の初期時刻における変位ｙ及び変位速度ｙ’の値を設定する手段、前記弾性振動体の所定部位の時刻ｔにおける変位ｙt及び変位速度ｙ’tの値を用いて当該所定部位の振幅ｙ0の時刻ｔにおける値を計算する手段、前記風速条件に基づいて特定された時刻ｔにおける風速と、振幅ｙ0の時刻ｔにおける値に基づいて空力減衰率δa及び応答円振動数ｎの時刻ｔにおける値を計算する手段、空力減衰率δaの時刻ｔにおける値を用いて前記非定常空気力の振動変位速度成分ｆaIの時刻ｔにおける値を計算すると共に、応答円振動数ｎの時刻ｔにおける値を用いて前記非定常空気力の振動変位成分ｆaRの時刻ｔにおける値を計算する手段、振動変位成分ｆaR及び振動変位速度成分ｆaIを変数として含む弾性振動体の振動方程式に、振動変位速度成分ｆaI及び振動変位成分ｆaRの時刻ｔにおける値を代入して得られた振動方程式を解いて、次の時刻ｔ+Δｔにおける前記所定部位の変位ｙt+Δt、変位速度ｙ’t+Δt及び加速度ｙ’’t+Δtの値を計算する手段、前記計算する手段を順次実行する時刻歴応答計算ルーチンを各時刻について繰り返し実行させる手段であって、２回目以降の時刻歴応答計算ルーチンの振幅ｙ0の値の計算で用いられる変位ｙt及び変位速度ｙ’tの値として、直前の時刻歴応答計算ルーチンで計算された変位ｙ及び変位速度ｙ’の値を用いて振幅ｙ0の値を計算させる手段、前記時刻歴応答計算ルーチンを各時刻について繰り返し実行して計算された値を用いて弾性振動体の時刻歴応答シミュレーションを出力する手段、として機能させる時刻歴応答シミュレーションプログラムである。
【００４８】
また、本願に係るさらに別の発明は、外力として非定常空気力を受ける弾性振動体の時刻歴応答シミュレーションを行う時刻歴応答シミュレーション装置であって、非定常空気力を受ける前記弾性振動体の所定部位の時刻ｔにおける変位ｙt及び変位速度ｙ’tの値を用いて前記所定部位の振幅ｙ0の時刻ｔにおける値を計算する第１計算手段と、振幅ｙ0の時刻ｔにおける値を用いて、振幅ｙ0の関数である空力減衰率δa及び応答円振動数ｎの時刻ｔにおける値を計算する第２計算手段と、空力減衰率δaの時刻ｔにおける値を用いて前記非定常空気力の振動変位速度成分ｆaIの時刻ｔにおける値を計算すると共に、応答円振動数ｎの時刻ｔにおける値を用いて前記非定常空気力の振動変位成分ｆaRの時刻ｔにおける値を計算する第３計算手段と、振動変位成分ｆaR及び振動変位速度成分ｆaIを変数として含む弾性振動体の振動方程式に、振動変位速度成分ｆaI及び振動変位成分ｆaRの時刻ｔにおける値を代入して得られた振動方程式を解いて、次の時刻ｔ+Δｔにおける前記所定部位の変位ｙt+Δt、変位速度ｙ’t+Δt及び加速度ｙ’’t+Δtの値を計算する第４計算手段と、前記計算する手段を順次実行する時刻歴応答計算ルーチンを各時刻について繰り返し実行させる手段であって、２回目以降の時刻歴応答計算ルーチンの振幅ｙ0の値の計算で用いられる変位ｙt及び変位速度ｙ’tの値として、直前の時刻歴応答計算ルーチンで計算された変位ｙ及び変位速度ｙ’の値を用いて振幅ｙ0の値を計算させる手段と、前記時刻歴応答計算ルーチンを各時刻について繰り返し実行して計算された値を用いて弾性振動体の時刻歴応答シミュレーションを出力する手段と、を備える時刻歴応答シミュレーション装置である。
【実施例】
【００４９】
以下、本発明に係る時刻歴応答シミュレーション方法、そのプログラム及び装置の実施例について、図面を参照しつつ詳細に説明する。ここでは、図１及び図２に示される箱桁断面の橋梁模型（空力弾性振動体）を例に、バネ吊り１自由度振動系の空力弾性応答のシミュレーションについて説明する。
【００５０】
まず、シミュレーションの前にあらかじめ、図１に示される形状及び寸法（ｍｍ）の橋梁模型の桁部（所定部位）について自由振動実験を行った。より具体的には、橋梁模型の桁部の一部に該当する、図２に示される二次元部分剛体模型をバネで吊った状態で風洞内に設置し、風を吹かせて応答を測定するという自由振動実験を行った（図３参照）。自由振動実験に用いた二次元部分剛体模型は、実物の１／６５の縮尺のものであり、長さが８００ｍｍ、質量ｍが１．５６ｋｇ、構造減衰率δsが０．００１９２２９、固有振動数ωが７．７５６Ｈｚというものであった。自由振動実験では、風速を徐々に変化させ、各風速における発散振動状況、減衰振動状況の振動波形をそのつど記録し、記録した振動波形から各風速における振幅ｙ0ごとの減衰率δを算出した。ここでは風速０．２ｍ／ｓごとに上記の記録を行った。そして、算出した減衰率δから構造減衰率δsを差し引いて、振幅ｙ0の関数である空力減衰率δa（＝ｆI（ｙ0））を求めた。また、記録した振動波形を基に、振幅ｙ0の関数である後述の応答円振動数ｎ（＝ｆR（ｙ0））を求めた。所定部位とは、橋梁の桁部、主塔部など、それぞれ固有の断面形状を有する全ての部分のことである。ここでは、連続高架橋の桁部という部位の時刻歴応答シミュレーションについて説明する。なお、必要に応じて各部位について自由振動実験を行い、各部位の空力減衰率δa、応答円振動数ｎを算出することになるが、その手順はここで説明するものと同様であるので、説明を省略する。
【００５１】
ここで、実験によって求めたデータを基に作成した１自由度たわみ応答（二次元模型の対風振動応答図）を示す（図４の実験値のグラフ参照）。この図では、横軸が換算風速Ｖr、縦軸が無次元化した振幅（２Ａ／Ｄ）である。
また、自由振動実験を行って求めた空力減衰率δaのグラフの一例として、渦励振時のピークを示す換算風速Ｖr＝５．７６における空力減衰率δaの応答図を示す（図５の実験値のグラフ参照）。そして、関数近似を行って、次の式（２１）で示されるような空力減衰率δaの近似式（近似関数）を得た。
下記の式から解かるように、ここでは、振幅値について区間を複数設定し（実施例は３区間）、各区間ごとに近似式を求めた。このようにして、振幅値の値に応じた近似式を用いて空力減衰率δaを求めるようにすると、より正確に空力減衰率δaを求めることができ、より正確なシミュレーションを行うことができる。また、この点は、後述の応答円振動数についても同様である。なお、他の換算風速においても同様に、必要に応じて区間を区切って空力減衰率δaの近似式（近似関数）を得たが、手順等は同様であるので、ここではその詳細な説明を省略する。
なお、近似化に用いた数値を直接用いてシミュレーションを行ってもよいが、シミュレーションのためには、関数化しておく方が都合が良く、好ましい。これは、空力減衰率についてだけでなく、後述の応答円振動数ｎについても同様である。
【数１０】
【００５２】
また、空力減衰率δaを求める際に、同時に各風速における応答円振動数ｎを求めた。一例として換算風速Ｖr＝５．７６に対する応答円振動数ｎの応答図を示す（図６の実験値のグラフ参照）。そして、関数近似を行って、次の式（２２）で示されるような応答円振動数ｎの近似式を得た。下記の式から解かるように、空力減衰率δaの場合と同様、振幅値について区間を複数設定し（ここでは２区間）、各区間ごとに近似式を求めた。この近似式を用いて、いわゆる非定常空気力の変位比例成分を得ることができる。なお、他の換算風速においても同様に、必要に応じて区間を区切って応答円振動数ｎの近似式を得たが、手順等は同様であるので、ここではその詳細な説明を省略する。
【数１１】
【００５３】
このようにして、各風速および振幅に対応した橋梁の所定部位の空力減衰率δaを表す近似式（近似関数）及び応答円振動数ｎを表す近似式（近似関数）を求めた上で、次に説明しいている各ステップから構成される時刻歴応答計算ルーチンを時々刻々と変化する各時刻について繰り返し実行し、時刻歴応答シミュレーションを行った。
【００５４】
以下、時刻歴応答シミュレーション方法、そのプログラム及び装置について説明する。
【００５５】
図９に示されるように、本シミュレーション方法を実行するシミュレーション装置10は、ＣＰＵ11等の演算プロセッサやメインメモリ12等のメモリを備え、ＯＳが例えばＷＩＮＤＯＷＳ（登録商標）ＸＰエディションであるパソコン本体13と、パソコン本体13に入出力インターフェース14を介して接続された出力手段であるディスプレイ15と、プリンタ（不図示）と、外部記憶装置であるハードディスク16と、ＣＤ−ＲＯＭの記録を読み取るためのＣＤ−ＲＯＭドライブ17（または各種記憶媒体のためのマルチドライブ）と、パソコン本体13に入出力インターフェース14を介して接続された入力手段であるマウス18やキーボード19とを備える周知のコンピュータである。なお、符号「20」はバスである（図９参照）。
そして、本発明にかかる時刻歴応答シミュレーションプログラムがハードディスク15又はマルチドライブにより読み取り可能なＣＤ−ＲＯＭ等の記憶媒体に記憶されており、本シミュレーション装置10にインストールされている。そして、主にＣＰＵ11やメインメモリ12が計算などのステップ処理部や出力するデータを生成する手段として機能し、メインメモリ12が入力された設定条件の記憶手段として機能する。また、計算によって算出された数値データや生成された画像データ等がハードディスク16やディスプレイ15などに出力される。なお、シミュレーション装置10として用いられているコンピュータ自体は周知の装置であるので、ここでは、その詳細な説明を省略する。
【００５６】
次に、時刻歴応答シミュレーション方法について説明する。本シミュレーション方法は、上述したように、時刻ｔにおける所定部位の振幅ｙ0tに対応する空気力が所定部位に瞬間瞬間作用しているという考え方を用いて空力弾性振動の応答を追跡し、シミュレーションするものである。
なお、時刻ｔにおける所定部位の振幅ｙ0とは、各時刻ｔにおける所定部位の実際の変位ｙ（図７の実線参照）ではなく、時刻ｔにおいて取り得た最大変位（正弦波振幅、図７の破線参照）のことである。別言すれば、時刻ｔにおける所定部位の時々刻々変化する振幅のことである。この振幅ｙ0は、図７の破線で示されるように変化する。
【００５７】
まず、シミュレーションの対象である弾性振動体の所定部位の条件を設定する。本実施例では、所定部位として桁部（図１参照）が設定され、桁部の条件として上述したような質量等の条件が設定される。そして、シミュレーションしようとしている状況に応じた風速条件（換算風速Ｖrの条件）を設定する。外力として非定常空気力を受ける状態を検討する場合は、時々刻々と変化する条件で設定がなされる。各時刻における具体的風速値（および／または換算風速Ｖr）をあらかじめ設定しておき必要に応じて読み取って用いることができるようにしてもよいし、各時刻における風速を算出するための関数を用意しておき必要に応じて各時刻における風速を算出して用いることができるようにしてもよい。これにより、シミュレーションの対象である橋梁の所定部位の所定時刻における風速（換算風速Ｖr）を定めることができる。また、最初の時刻における所定部位の変位、変位速度及び変位加速度を設定する。なお、各種条件設定については、キーボード等の入力手段を用いた入力や、あらかじめ記憶させておいた種々の条件の中からの選択など、種々の方法によって行うことができる。各種設定条件は、必要に応じて記憶手段に記憶されて用いられる。
【００５８】
そして、インストールされたプログラムに基づいて、以降に説明するような計算を行う。
まず、所定時刻ｔ1（最初の時刻歴応答計算ルーチンでは最初の時刻）における風速を特定する。具体的には、先に設定した風速条件に基づいて、風速または換算風速Ｖrの具体的値が特定される。また、時刻ｔ1の所定部位の変位ｙt1及び変位速度ｙ’t1を式（１０）中の変位ｙ及び変位速度ｙ’に代入して振幅ｙ0の計算する。これにより、非定常空気力（外力）を受ける橋梁の桁部の振幅ｙ0の時刻ｔ1における値が求められる（第１ステップ）。時刻歴応答シミュレーションプログラムは、このような計算をする計算ルーチンを備えている（第１計算手段）。
【数１２】
【００５９】
なお、変位速度が０のときは、所定部位の実際の変位は振幅ｙ0に一致するので、例えば時刻ｔ1における変位速度が０である場合は、式（１０）を用いずに、実際の変位ｙt1を時刻ｔ1の振幅ｙ0としてもよい。
【００６０】
次に、空力減衰率δa及び応答円振動数ｎの値を計算する（第２ステップ）。まず、第１ステップで定められた風速条件及び振幅ｙ0の値に応じて、ここでの計算で用いる空力減衰率δaの近似式および応答円振動数ｎの近似式を特定する。例えば、風速が「換算風速Ｖr＝５．７６」に対する風速であり、振幅ｙ0が０．０３であれば、空力減衰率δaの計算に用いられる近似式は、上記式（２１）中の近似式２に特定され、応答円振動数ｎの計算に用いられる近似式は、上記の式（２２）中の近似式１に特定される。なお、完全に一致する条件に対応した近似式がない場合は、最も近い条件に対応した近似式を、用いる近似式として特定してもよい。このようにして特定された式に、時刻ｔ1における振幅ｙ0t1の値を代入して、当該時刻ｔ1における所定部位の空力減衰率δa及び応答円振動数ｎの値を計算する（第２ステップ）。時刻歴応答シミュレーションプログラムは、このような計算をする計算ルーチンを備えている（第２計算手段）。
【００６１】
続いて、時刻ｔ1の時点で所定部位に加わっている外力を求める。具体的には、第２ステップで求めた空力減衰率δaの具体的数値を式（１６）に代入して振動変位速度成分ｆaIの具体的数値を計算し、第２ステップで求めた応答円振動数ｎの具体的数値を式（１８）に代入して振動変位成分ｆaRの具体的数値を計算する（第３ステップ）。時刻歴応答シミュレーションプログラムは、このような計算をする計算ルーチンを備えている（第３計算手段）。
【数１３】
【００６２】
そして、計算された振動変位速度成分ｆaI及び振動変位成分ｆaRの具体的数値を式（６）に代入することによって、時刻ｔ1における振動方程式が特定される。
【００６３】
この特定された振動方程式を、直接時間積分法で解いて、時間刻みに応じた振幅の時刻歴を計算する（第４ステップ）。時刻歴応答シミュレーションプログラムは、このような計算をする計算ルーチンを備えている（第４計算手段）。
【００６４】
本実施例では、不規則外力を受ける構造物の応答計算を数値計算によって行う際に用いられる直接時間積分法の１つであるニューマークβ法を用いて、振幅の時刻歴を計算する。
ニューマークβ法は、不規則外力を受ける構造物の応答計算を行う際に用いられる数値計算方法の１つである。この方法を用い、振動方程式の外力項を空気力とすることで、時系列での応答推定を行う。
ニューマークβ法は、時刻ｔ+Δｔ（＝例えば時刻ｔ2）のときの所定部位の実際の変位を求めるために、時刻ｔ（例えば時刻ｔ1）のときの値とΔｔ時間後である時刻ｔ+Δｔ（例えばｔ2）のときの値を使用して繰り返し計算を行う解法である。使用する方程式は、次の式（２４）である。
【数１４】
【００６５】
ここで、上記式（２４）を式（６）の形に置き換えると、次の式（２５）のようになり、さらにマトリクス表示すると、次の式（２６）のようになる。
【数１５】
【００６６】
つまり、第４ステップでは、上記式（２６）を用いて、時刻ｔ2（＝ｔ+Δｔ）における前記所定部位の変位ｙt+Δt、変位速度ｙ’t+Δt及び変位加速度ｙ’’t+Δtを計算することができる（第４ステップ、第４計算手段）。
これにより、最初の時刻歴応答計算ルーチンＲ1において、所定部位の時刻ｔ2（＝ｔ+Δｔ）における実際に変位が求められたことになる。
【００６７】
また、シミュレーションプログラムは、計算された変位ｙt2、変位速度ｙ’t2及び変位加速度ｙ’’t2の値を、次の時刻t3における所定部位の変位等を求める次の時刻歴応答計算ルーチンＲ2において、変位ｙt及び変位速度ｙ’tの値として用いて計算を行わせるようになっている。例えば、計算された時刻ｔ2（＝ｔ+Δｔ）の変位ｙの値及び変位速度ｙ’の値は、時刻歴応答計算ルーチンＲ2の第１ステップにおいて、式（１０）の変位ｙ及び変位速度ｙ’の値として代入され、これにより時刻ｔ3（＝ｔ+２×Δｔ）における所定部位の振幅ｙ0が計算される（第１ステップ）。また、計算された変位ｙt2、変位速度ｙ’t2及び変位加速度ｙ’’t2の値は、時刻歴応答計算ルーチンＲ2の第４ステップにおいて、式（２６）の変位ｙt、変位速度ｙ’t及び変位加速度ｙ’’tの値として代入され、これにより、時刻t3における所定部位の変位ｙt+Δt、変位速度ｙ’t+Δt及び変位加速度ｙ’’t+Δtが計算される（第４ステップ）。
【００６８】
このように、２回目以降のｘ回目の時刻歴応答計算ルーチンＲxの第１ステップでは、直前の時刻歴応答計算ルーチンＲx-1の第４ステップで計算された変位ｙtx及び変位速度ｙ’txの値を、変位ｙ及び変位速度ｙ’の値として用いて（式（１０）等参照）、次の時刻ｔx+1における振幅ｙ0を計算する。
また、２回目以降の時刻歴応答計算ルーチンＲxの第４ステップでは、直前の時刻歴応答計算ルーチンＲx-1の第４ステップで計算された変位ｙtx、変位速度ｙ’tx及び変位加速度ｙ’’txの値を、当該第４ステップにおける変位ｙt、変位速度ｙ’t及び変位加速度ｙ’’tの値として用いて次の時刻ｔx+1における変位ｙ、変位速度ｙ’及び変位加速度ｙ’’の値を計算する。
時刻歴応答シミュレーションプログラムは、このように時刻歴応答計算ルーチンを各時刻について繰り返し実行させる計算ルーチンであって、２回目以降の時刻歴応答計算ルーチンの各計算で用いられる変位ｙt、変位速度ｙ’t及び変位加速度ｙ’’の値として、直前の時刻歴応答計算ルーチンで計算された変位ｙ、変位速度ｙ’及び変位加速度ｙ’’の値を用いて振幅ｙ0の値を計算させる計算ルーチンを、計算手段として備えている。
【００６９】
このようなステップからなる時刻歴応答計算ルーチンを各時刻について繰り返し実行することによって、各時刻における応答（実際の変位等）を計算し、計算された変位等の値を用いて、Δｔ時間刻みの時刻歴応答シミュレーションを行う。シミュレーションの結果は、ディスプレイ15やプリンタなどに出力される。ハードディスク16などの記憶媒体に記録される態様で出力することもできる。時刻歴応答シミュレーションプログラムは、このような機器に時刻歴応答シミュレーションの結果を出力させる処理ルーチンを備えている。
【００７０】
そして、上記のようなシミュレーションを行った結果、図８（ａ）に示されるような時刻歴応答解析波形が得られた。また、比較のために、実物の１／１３０の縮尺の全橋模型について、シミュレーションの条件と同じ条件で風洞実験を行って、時刻歴応答解析波形を得た。図８（ｂ）に、風洞実験での振幅の測定値に基づいて作成した振幅波形を示す。両波形を比較するとわかるように、両波形はほぼ一致していた。
【００７１】
そして、図５中の推定値のグラフは、図８（ａ）に示される時刻歴応答解析波形を基に算出した空力減衰率の推定値のグラフである。図５から明らかなように、空力減衰率の近似値及び推定値は、いずれも、実験値とほぼ一致していた。
【００７２】
また、時刻歴解析を各風速で解析し、定常振動応答をプロットして、測定値との比較を行ったものが図４である。この図に示されるように、時刻歴応答解析から求めた定常振動応答推定値は、若干のずれはあるものの実験での測定値をほぼ再現していた。
【００７３】
このように、本願発明によれば、減衰状態あるいは発散状態にある所定部位の空力弾性振動について、時々刻々の振幅ｙ0を計算することができ、計算された振幅ｙ0を用いて非定常空気力を計算することができる。時々刻々の非定常空気力を計算することができれば、計算された非定常空気力を振動方程式に代入することによって時々刻々の振動方程式を定めて、定められた振動方程式を解くことによって、空力弾性振動の時刻歴応答解析シミュレーションをすることができる。
【００７４】
本願の発明は、以上の実施例に限定されるものではなく、その要旨を逸脱しない範囲において、種々の変形が可能である。
【００７５】
上記実施例では、箱桁断面の橋梁模型（空力弾性振動体）を例にシミュレーションについて説明したが、本願の発明にかかるシミュレーション方法及び装置は、例えば、橋梁の桁部、主塔部など、弾性振動する種々の部位についての時刻歴応答シミュレーションに用いることができる。
【００７６】
また、上記実施例では、箱桁断面の橋梁模型を例に、バネ吊り１自由度振動系の空力弾性応答のシミュレーションについて説明したが、本願の発明にかかるシミュレーション方法及び装置は、その他にも、バネ吊り２自由度振動系の空力弾性応答や全橋の空力弾性振動応答、乱流中での橋桁の空力弾性振動応答のシミュレーションなど、種々の振動に関する空力弾性応答のシミュレーションに用いることができる。
【００７７】
また、上記実施例の自由振動実験では、風速を徐々に変化させ、０．２ｍ／ｓごとに発散振動状況、減衰振動状況の振動波形を記録し、記録した振動波形から各風速における振幅ごとの減衰率を算出したが、波形を測定する風速の間隔は必ずしも０．２ｍ／ｓごとでなければならないわけではなく、これにより大きくても良いし、小さくても良い。
【図面の簡単な説明】
【００７８】
【図１】本発明のシミュレーション方法のシミュレーション対象である連続高架橋を示すものであり、橋脚上に配置された桁部の概略形状を示す斜視図である。
【図２】図１に示される橋梁の桁部の断面を示す断面図である。
【図３】本実施例の自由振動実験で用いた装置の概略構成を示す構成図である。
【図４】１自由度たわみ応答図（二次元模型の対風振動応答図）である。
【図５】換算風速Ｖr＝５．７６における空力減衰率δaの応答図である。
【図６】換算風速Ｖr＝５．７６に対する応答円振動数ｎの応答図である。
【図７】非定常空気力を受けた弾性振動体の過渡応答時の振動振幅を示すグラフである。
【図８】（ａ）は、換算風速Ｖr＝５．７６における時刻歴応答解析波形を示す図であり、（ｂ）は、実際の振幅の測定値に基づいて作成した振幅波形を示す図である。
【図９】本発明のシミュレーション装置の構成を示す構成図である。
【特許請求の範囲】
【請求項１】
外力として非定常空気力を受ける弾性振動体の時刻歴応答シミュレーション方法であって、
非定常空気力を受ける前記弾性振動体の所定部位の時刻ｔにおける変位ｙt及び変位速度ｙ’tの値を用いて前記所定部位の振幅ｙ0の時刻ｔにおける値を計算する第１ステップと、
振幅ｙ0の時刻ｔにおける値を用いて、振幅ｙ0の関数である空力減衰率δa及び応答円振動数ｎの時刻ｔにおける値を計算する第２ステップと、
空力減衰率δaの時刻ｔにおける値を用いて前記非定常空気力の振動変位速度成分ｆaIの時刻ｔにおける値を計算すると共に、応答円振動数ｎの時刻ｔにおける値を用いて前記非定常空気力の振動変位成分ｆaRの時刻ｔにおける値を計算する第３ステップと、
振動変位成分ｆaR及び振動変位速度成分ｆaIを変数として含む弾性振動体の振動方程式に振動変位速度成分ｆaI及び振動変位成分ｆaRの時刻ｔにおける値を代入して得られた振動方程式を解いて、次の時刻ｔ+Δｔにおける前記所定部位の変位ｙt+Δt、変位速度ｙ’t+Δt及び加速度ｙ’’t+Δtの値を計算する第４ステップと、を備え、
これらのステップを含む時刻歴応答計算ルーチンを各時刻について繰り返し実行して計算された値を基に弾性振動体の時刻歴応答シミュレーションを行うものであり、
２回目以降の時刻歴応答計算ルーチンの第１ステップでは、当該第１ステップにおける変位ｙt及び変位速度ｙ’tの値として、直前の時刻歴応答計算ルーチンの第４ステップで計算された変位ｙ及び変位速度ｙ’の値を用いて振幅ｙ0の対応時刻における値を計算するものである時刻歴応答シミュレーション方法。
【請求項２】
２回目以降の時刻歴応答計算ルーチンの第４ステップは、当該第４ステップにおける変位ｙt、変位速度ｙ’t及び加速度ｙ’’tの値として、直前の時刻歴応答計算ルーチンの第４ステップで計算された変位ｙ、変位速度ｙ’及び加速度ｙ’’の値を用いて、次の時刻ｔ+Δｔにおける前記所定部位の変位ｙt+Δt、変位速度ｙ’t+Δt及び加速度ｙ’’t+Δtの値を計算するものである請求項１に記載の時刻歴応答シミュレーション方法。
【請求項３】
前記所定部位の振幅ｙ0の時刻ｔにおける値を求めるための前記変位ｙ及び変位速度ｙ’の関数は、時刻ｔにおける変位ｙtの自乗値と、時刻ｔにおける変位速度ｙ’tを応答円振動数ｎで除した値の自乗値との和の平方根を振幅ｙ0の時刻ｔにおける値として求めることができるものである請求項１又は請求項２に記載の時刻歴応答シミュレーション方法。
【請求項４】
前記振動方程式は、下記の式（６）であり、前記第４ステップで振動変位速度成分ｆaI及び振動変位成分ｆaRの値が代入された振動方程式の解法は直接時間積分法であり、
前記第１ステップで前記振幅ｙ0の値を計算するために用いられる変位ｙ及び変位速度ｙ’の関数は、下記の式（１０）であり、
前記第３ステップで前記振動変位速度成分ｆaIの値を計算するために用いられる前記空力減衰率δaの関数は、下記の式（１６）であり、
前記第３ステップで前記振動変位成分ｆaRの値を計算するために用いられる前記応答円振動数ｎの関数は、下記の式（１８）である請求項１から請求項３のいずれか一項に記載の時刻歴応答シミュレーション方法。
【数１】
【請求項５】
前記空力減衰率δaを求めるための前記振幅ｙ0の関数（δa＝ｆI）及び前記応答円振動数ｎを求めるための前記振幅ｙ0の関数（ｎ＝ｆR）は、前記弾性振動体について行う振動実験によって求められるものである請求項１から請求項４のいずれか一項に記載の時刻歴応答シミュレーション方法。
【請求項６】
前記関数（δa＝ｆI）及び前記関数（ｎ＝ｆR）は、前記振幅ｙ0及び風速の関数である請求項５に記載の時刻歴応答シミュレーション方法。
【請求項７】
外力として非定常空気力を受ける弾性振動体の時刻歴応答シミュレーションをするためのコンピュータを、
シミュレーションの対象である前記弾性振動体の所定部位の条件を設定する手段と、
風速条件を設定する手段、
前記弾性振動体の所定部位の初期時刻における変位ｙ及び変位速度ｙ’の値を設定する手段、
前記弾性振動体の所定部位の時刻ｔにおける変位ｙt及び変位速度ｙ’tの値を用いて当該所定部位の振幅ｙ0の時刻ｔにおける値を計算する手段、
前記風速条件に基づいて特定された時刻ｔにおける風速と、振幅ｙ0の時刻ｔにおける値に基づいて空力減衰率δa及び応答円振動数ｎの時刻ｔにおける値を計算する手段、
空力減衰率δaの時刻ｔにおける値を用いて前記非定常空気力の振動変位速度成分ｆaIの時刻ｔにおける値を計算すると共に、応答円振動数ｎの時刻ｔにおける値を用いて前記非定常空気力の振動変位成分ｆaRの時刻ｔにおける値を計算する手段、
振動変位成分ｆaR及び振動変位速度成分ｆaIを変数として含む弾性振動体の振動方程式に、振動変位速度成分ｆaI及び振動変位成分ｆaRの時刻ｔにおける値を代入して得られた振動方程式を解いて、次の時刻ｔ+Δｔにおける前記所定部位の変位ｙt+Δt、変位速度ｙ’t+Δt及び加速度ｙ’’t+Δtの値を計算する手段、
前記計算する手段を順次実行する時刻歴応答計算ルーチンを各時刻について繰り返し実行させる手段であって、２回目以降の時刻歴応答計算ルーチンの振幅ｙ0の値の計算で用いられる変位ｙt及び変位速度ｙ’tの値として、直前の時刻歴応答計算ルーチンで計算された変位ｙ及び変位速度ｙ’の値を用いて振幅ｙ0の値を計算させる手段、
前記時刻歴応答計算ルーチンを各時刻について繰り返し実行して計算された値を用いて弾性振動体の時刻歴応答シミュレーションを出力する手段
として機能させる時刻歴応答シミュレーションプログラム。
【請求項８】
外力として非定常空気力を受ける弾性振動体の時刻歴応答シミュレーションを行う時刻歴応答シミュレーション装置であって、
非定常空気力を受ける前記弾性振動体の所定部位の時刻ｔにおける変位ｙt及び変位速度ｙ’tの値を用いて前記所定部位の振幅ｙ0の時刻ｔにおける値を計算する第１計算手段と、
振幅ｙ0の時刻ｔにおける値を用いて、振幅ｙ0の関数である空力減衰率δa及び応答円振動数ｎの時刻ｔにおける値を計算する第２計算手段と、
空力減衰率δaの時刻ｔにおける値を用いて前記非定常空気力の振動変位速度成分ｆaIの時刻ｔにおける値を計算すると共に、応答円振動数ｎの時刻ｔにおける値を用いて前記非定常空気力の振動変位成分ｆaRの時刻ｔにおける値を計算する第３計算手段と、
振動変位成分ｆaR及び振動変位速度成分ｆaIを変数として含む弾性振動体の振動方程式に、振動変位速度成分ｆaI及び振動変位成分ｆaRの時刻ｔにおける値を代入して得られた振動方程式を解いて、次の時刻ｔ+Δｔにおける前記所定部位の変位ｙt+Δt、変位速度ｙ’t+Δt及び加速度ｙ’’t+Δtの値を計算する第４計算手段と、
前記計算する手段を順次実行する時刻歴応答計算ルーチンを各時刻について繰り返し実行させる手段であって、２回目以降の時刻歴応答計算ルーチンの振幅ｙ0の値の計算で用いられる変位ｙt及び変位速度ｙ’tの値として、直前の時刻歴応答計算ルーチンで計算された変位ｙ及び変位速度ｙ’値を用いて振幅ｙ0の値を計算させる手段と、
前記時刻歴応答計算ルーチンを各時刻について繰り返し実行して計算された値を用いて弾性振動体の時刻歴応答シミュレーションを出力する手段と、を備える時刻歴応答シミュレーション装置。
【請求項１】
外力として非定常空気力を受ける弾性振動体の時刻歴応答シミュレーション方法であって、
非定常空気力を受ける前記弾性振動体の所定部位の時刻ｔにおける変位ｙt及び変位速度ｙ’tの値を用いて前記所定部位の振幅ｙ0の時刻ｔにおける値を計算する第１ステップと、
振幅ｙ0の時刻ｔにおける値を用いて、振幅ｙ0の関数である空力減衰率δa及び応答円振動数ｎの時刻ｔにおける値を計算する第２ステップと、
空力減衰率δaの時刻ｔにおける値を用いて前記非定常空気力の振動変位速度成分ｆaIの時刻ｔにおける値を計算すると共に、応答円振動数ｎの時刻ｔにおける値を用いて前記非定常空気力の振動変位成分ｆaRの時刻ｔにおける値を計算する第３ステップと、
振動変位成分ｆaR及び振動変位速度成分ｆaIを変数として含む弾性振動体の振動方程式に振動変位速度成分ｆaI及び振動変位成分ｆaRの時刻ｔにおける値を代入して得られた振動方程式を解いて、次の時刻ｔ+Δｔにおける前記所定部位の変位ｙt+Δt、変位速度ｙ’t+Δt及び加速度ｙ’’t+Δtの値を計算する第４ステップと、を備え、
これらのステップを含む時刻歴応答計算ルーチンを各時刻について繰り返し実行して計算された値を基に弾性振動体の時刻歴応答シミュレーションを行うものであり、
２回目以降の時刻歴応答計算ルーチンの第１ステップでは、当該第１ステップにおける変位ｙt及び変位速度ｙ’tの値として、直前の時刻歴応答計算ルーチンの第４ステップで計算された変位ｙ及び変位速度ｙ’の値を用いて振幅ｙ0の対応時刻における値を計算するものである時刻歴応答シミュレーション方法。
【請求項２】
２回目以降の時刻歴応答計算ルーチンの第４ステップは、当該第４ステップにおける変位ｙt、変位速度ｙ’t及び加速度ｙ’’tの値として、直前の時刻歴応答計算ルーチンの第４ステップで計算された変位ｙ、変位速度ｙ’及び加速度ｙ’’の値を用いて、次の時刻ｔ+Δｔにおける前記所定部位の変位ｙt+Δt、変位速度ｙ’t+Δt及び加速度ｙ’’t+Δtの値を計算するものである請求項１に記載の時刻歴応答シミュレーション方法。
【請求項３】
前記所定部位の振幅ｙ0の時刻ｔにおける値を求めるための前記変位ｙ及び変位速度ｙ’の関数は、時刻ｔにおける変位ｙtの自乗値と、時刻ｔにおける変位速度ｙ’tを応答円振動数ｎで除した値の自乗値との和の平方根を振幅ｙ0の時刻ｔにおける値として求めることができるものである請求項１又は請求項２に記載の時刻歴応答シミュレーション方法。
【請求項４】
前記振動方程式は、下記の式（６）であり、前記第４ステップで振動変位速度成分ｆaI及び振動変位成分ｆaRの値が代入された振動方程式の解法は直接時間積分法であり、
前記第１ステップで前記振幅ｙ0の値を計算するために用いられる変位ｙ及び変位速度ｙ’の関数は、下記の式（１０）であり、
前記第３ステップで前記振動変位速度成分ｆaIの値を計算するために用いられる前記空力減衰率δaの関数は、下記の式（１６）であり、
前記第３ステップで前記振動変位成分ｆaRの値を計算するために用いられる前記応答円振動数ｎの関数は、下記の式（１８）である請求項１から請求項３のいずれか一項に記載の時刻歴応答シミュレーション方法。
【数１】
【請求項５】
前記空力減衰率δaを求めるための前記振幅ｙ0の関数（δa＝ｆI）及び前記応答円振動数ｎを求めるための前記振幅ｙ0の関数（ｎ＝ｆR）は、前記弾性振動体について行う振動実験によって求められるものである請求項１から請求項４のいずれか一項に記載の時刻歴応答シミュレーション方法。
【請求項６】
前記関数（δa＝ｆI）及び前記関数（ｎ＝ｆR）は、前記振幅ｙ0及び風速の関数である請求項５に記載の時刻歴応答シミュレーション方法。
【請求項７】
外力として非定常空気力を受ける弾性振動体の時刻歴応答シミュレーションをするためのコンピュータを、
シミュレーションの対象である前記弾性振動体の所定部位の条件を設定する手段と、
風速条件を設定する手段、
前記弾性振動体の所定部位の初期時刻における変位ｙ及び変位速度ｙ’の値を設定する手段、
前記弾性振動体の所定部位の時刻ｔにおける変位ｙt及び変位速度ｙ’tの値を用いて当該所定部位の振幅ｙ0の時刻ｔにおける値を計算する手段、
前記風速条件に基づいて特定された時刻ｔにおける風速と、振幅ｙ0の時刻ｔにおける値に基づいて空力減衰率δa及び応答円振動数ｎの時刻ｔにおける値を計算する手段、
空力減衰率δaの時刻ｔにおける値を用いて前記非定常空気力の振動変位速度成分ｆaIの時刻ｔにおける値を計算すると共に、応答円振動数ｎの時刻ｔにおける値を用いて前記非定常空気力の振動変位成分ｆaRの時刻ｔにおける値を計算する手段、
振動変位成分ｆaR及び振動変位速度成分ｆaIを変数として含む弾性振動体の振動方程式に、振動変位速度成分ｆaI及び振動変位成分ｆaRの時刻ｔにおける値を代入して得られた振動方程式を解いて、次の時刻ｔ+Δｔにおける前記所定部位の変位ｙt+Δt、変位速度ｙ’t+Δt及び加速度ｙ’’t+Δtの値を計算する手段、
前記計算する手段を順次実行する時刻歴応答計算ルーチンを各時刻について繰り返し実行させる手段であって、２回目以降の時刻歴応答計算ルーチンの振幅ｙ0の値の計算で用いられる変位ｙt及び変位速度ｙ’tの値として、直前の時刻歴応答計算ルーチンで計算された変位ｙ及び変位速度ｙ’の値を用いて振幅ｙ0の値を計算させる手段、
前記時刻歴応答計算ルーチンを各時刻について繰り返し実行して計算された値を用いて弾性振動体の時刻歴応答シミュレーションを出力する手段
として機能させる時刻歴応答シミュレーションプログラム。
【請求項８】
外力として非定常空気力を受ける弾性振動体の時刻歴応答シミュレーションを行う時刻歴応答シミュレーション装置であって、
非定常空気力を受ける前記弾性振動体の所定部位の時刻ｔにおける変位ｙt及び変位速度ｙ’tの値を用いて前記所定部位の振幅ｙ0の時刻ｔにおける値を計算する第１計算手段と、
振幅ｙ0の時刻ｔにおける値を用いて、振幅ｙ0の関数である空力減衰率δa及び応答円振動数ｎの時刻ｔにおける値を計算する第２計算手段と、
空力減衰率δaの時刻ｔにおける値を用いて前記非定常空気力の振動変位速度成分ｆaIの時刻ｔにおける値を計算すると共に、応答円振動数ｎの時刻ｔにおける値を用いて前記非定常空気力の振動変位成分ｆaRの時刻ｔにおける値を計算する第３計算手段と、
振動変位成分ｆaR及び振動変位速度成分ｆaIを変数として含む弾性振動体の振動方程式に、振動変位速度成分ｆaI及び振動変位成分ｆaRの時刻ｔにおける値を代入して得られた振動方程式を解いて、次の時刻ｔ+Δｔにおける前記所定部位の変位ｙt+Δt、変位速度ｙ’t+Δt及び加速度ｙ’’t+Δtの値を計算する第４計算手段と、
前記計算する手段を順次実行する時刻歴応答計算ルーチンを各時刻について繰り返し実行させる手段であって、２回目以降の時刻歴応答計算ルーチンの振幅ｙ0の値の計算で用いられる変位ｙt及び変位速度ｙ’tの値として、直前の時刻歴応答計算ルーチンで計算された変位ｙ及び変位速度ｙ’値を用いて振幅ｙ0の値を計算させる手段と、
前記時刻歴応答計算ルーチンを各時刻について繰り返し実行して計算された値を用いて弾性振動体の時刻歴応答シミュレーションを出力する手段と、を備える時刻歴応答シミュレーション装置。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【公開番号】特開２０１０−１１２７６９（Ｐ２０１０−１１２７６９Ａ）
【公開日】平成２２年５月２０日（２０１０．５．２０）
【国際特許分類】
【出願番号】特願２００８−２８３７９８（Ｐ２００８−２８３７９８）
【出願日】平成２０年１１月５日（２００８．１１．５）
【出願人】（５０４１７４１３５）国立大学法人九州工業大学 (489)
【Ｆターム（参考）】
【公開日】平成２２年５月２０日（２０１０．５．２０）
【国際特許分類】
【出願日】平成２０年１１月５日（２００８．１１．５）
【出願人】（５０４１７４１３５）国立大学法人九州工業大学 (489)
【Ｆターム（参考）】
[ Back to top ]