JP2010112769A - 時刻歴応答シミュレーション方法、そのプログラム及び装置 - Google Patents
時刻歴応答シミュレーション方法、そのプログラム及び装置 Download PDFInfo
- Publication number
- JP2010112769A JP2010112769A JP2008283798A JP2008283798A JP2010112769A JP 2010112769 A JP2010112769 A JP 2010112769A JP 2008283798 A JP2008283798 A JP 2008283798A JP 2008283798 A JP2008283798 A JP 2008283798A JP 2010112769 A JP2010112769 A JP 2010112769A
- Authority
- JP
- Japan
- Prior art keywords
- time
- displacement
- value
- vibration
- amplitude
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Images
Landscapes
- Aerodynamic Tests, Hydrodynamic Tests, Wind Tunnels, And Water Tanks (AREA)
Abstract
【課題】空力弾性振動の時刻歴応答解析シミュレーション方法を提供すること。
【解決手段】弾性振動体が受ける非定常空気力の振動変位成分faR及び振動変位速度成分faIを変数として含む振動方程式を用いた、弾性振動体の時刻歴応答シミュレーション方法であって、弾性振動体の所定部位の時刻tの変位及び変位速度の値から時刻tの振幅y0を計算し、振幅y0の値をから空力減衰率δa及び円振動数nの値を計算し、空力減衰率δaの値から振動変位速度成分faIを計算すると共に円振動数nの値から振動変位成分faRを計算し、振動変位速度成分faI及び振動変位成分faRの値が代入された振動方程式を解いて、時刻t+Δtの所定部位の変位、変位速度及び加速度を計算する時刻歴応答計算ルーチンを繰り返し実行する方法である。
【選択図】図1
【解決手段】弾性振動体が受ける非定常空気力の振動変位成分faR及び振動変位速度成分faIを変数として含む振動方程式を用いた、弾性振動体の時刻歴応答シミュレーション方法であって、弾性振動体の所定部位の時刻tの変位及び変位速度の値から時刻tの振幅y0を計算し、振幅y0の値をから空力減衰率δa及び円振動数nの値を計算し、空力減衰率δaの値から振動変位速度成分faIを計算すると共に円振動数nの値から振動変位成分faRを計算し、振動変位速度成分faI及び振動変位成分faRの値が代入された振動方程式を解いて、時刻t+Δtの所定部位の変位、変位速度及び加速度を計算する時刻歴応答計算ルーチンを繰り返し実行する方法である。
【選択図】図1
Description
本発明は、時刻歴応答解析に関するものであり、特に空力弾性振動の時刻歴応答解析シミュレーション方法、そのプログラム及び装置に関するものである。
時刻歴応答解析は、耐震工学の分野では古い歴史がある。耐震工学の分野における時刻歴応答解析とは、地震振動(地震加速度)などの外力が与えられた構造物の所定位置の応答加速度、速度、変位などを計算することである。そして、計算された解析結果は、構造物の耐震性の判断に用いられている。
ところで、橋梁の耐風工学の分野では、橋梁、特に長大橋梁の建設の際、耐風特性の検討が行われる(特許文献1参照)。耐風特性の検討では、例えば、風に代表される外力(空気力)を受ける橋梁の所定位置の応答(変位)が検討される。そして、検討結果が橋梁の耐風設計に用いられる。このような耐風特性の検討のために、時刻歴応答解析を用いることが考えられるが、現実には時刻歴応答解析を提案する流れはなかった。
時刻歴応答解析を行うためには、構造物に作用する外力の特性を十分に把握して評価できる必要があるところ、橋梁等の構造物が受ける外力である空気力の特性を把握することは困難であるということが原因であると考えられる。つまり、空力弾性振動のような自励振動における空気力の評価は、振動振幅に対して行わざるを得ず、過渡応答的な振動における各瞬間の空気力を評価することは困難であるため、耐風特性の検討に時刻歴応答解析を用いることは困難であると考えられた。
このようなことから、現在、橋梁全体の挙動は、一般に全橋模型を用いた風洞実験によって把握されており、橋梁の耐風特性もこの全橋模型を用いた風洞実験によって検討されている。
特開平8−184525号公報
ところが、風洞実験の場合、構造物の構造特性を種々変更して実験を行うには限界がある。そして、風洞実験では、乱流中での構造物の応答挙動を推定することが非常に困難であり、耐風特性の検討に関しても統計的な検討にならざるを得ないという問題がある。また、風洞実験に使用できる空間は限られた空間であるので、風洞実験では、実物の橋梁に対して大幅に縮小された全橋模型を使用せざるを得ず、実験精度が低下しやすいという問題もある。
なお、対風応答に関するシミュレーションとしては、有限要素法(FEM)を用いたフラッタ解析がある。フラッタとは、元来、航空機の翼などの構造物に作用する外力(空力)、弾性力及び慣性力が組み合わさって生じる自励振動のことである。フラッタ解析によれば、このようなフラッタ(自励振動)の発生の有無を判断することはできるが、外力(例えば、空気力、減衰力、弾性力)を受けた構造物のフラッタ(自励振動)の振幅(応答振幅)を追跡することはできない。
本発明は、上記問題点に鑑みてなされたものであり、空力弾性振動の時刻歴応答解析シミュレーション方法およびその装置を提供することを課題とする。
本願発明の発明者は、上記課題に鑑み、減衰状態や発散状態における構造物の時々刻々の振幅に基づいて非定常空気力(瞬時作用空気力)を求める手法について、後述するとおり鋭意研究を重ねた結果、次のような本願発明に想到するに至った。
つまり、本願発明は、外力として非定常空気力を受ける弾性振動体の時刻歴応答シミュレーション方法であって、非定常空気力を受ける前記弾性振動体の所定部位の時刻tにおける変位yt及び変位速度y’tの値を用いて前記所定部位の振幅y0の時刻tにおける値を計算する第1ステップと、振幅y0の時刻tにおける値を用いて、振幅y0の関数である空力減衰率δa及び応答円振動数nの時刻tにおける値を計算する第2ステップと、空力減衰率δaの時刻tにおける値を用いて前記非定常空気力の振動変位速度成分faIの時刻tにおける値を計算すると共に、応答円振動数nの時刻tにおける値を用いて前記非定常空気力の振動変位成分faRの時刻tにおける値を計算する第3ステップと、振動変位成分faR及び振動変位速度成分faIを変数として含む弾性振動体の振動方程式に振動変位速度成分faI及び振動変位成分faRの時刻tにおける値を代入して得られた振動方程式を解いて、次の時刻t+Δtにおける前記所定部位の変位yt+Δt、変位速度y’t+Δt及び加速度y’’t+Δtの値を計算する第4ステップと、を備え、これらのステップを含む時刻歴応答計算ルーチンを各時刻について繰り返し実行して計算された値を基に弾性振動体の時刻歴応答シミュレーションを行うものであり、2回目以降の時刻歴応答計算ルーチンの第1ステップでは、当該第1ステップにおける変位yt及び変位速度y’tの値として、直前の時刻歴応答計算ルーチンの第4ステップで計算された変位y及び変位速度y’の値を用いて振幅y0の対応時刻における値を計算するものである。
ここで、非定常空気力の振動変位速度成分faI及び非定常空気力の振動変位成分faRを、次の式(5)のようにおく。そうすると、上記式(1)の振動方程式を、下記の式(6)のように表すことができる。このように、振動変位速度成分faI及び非定常空気力の振動変位成分faRを式(5)のようにおいて振動方程式を式(6)のように表す考え方は好ましいものである。
ところで、式(6)では、非定常空気力の振動変位速度成分faIが減衰項の一部を構成し、非定常空気力の振動変位成分faRが復元力項の一部を構成している。この振動方程式を解くためには、振動変位速度成分faI及び振動変位成分faRの具体的な数値を把握する必要がある。
そこで、非定常空気力の振動変位速度成分faI及び振動変位成分faRの具体的な数値を把握する方法について検討した。その結果、空力弾性振動体の所定部位の空力減衰率δaが解かれば振動変位速度成分faIの具体的数値を把握することができ、空力弾性振動体の所定部位の応答円振動数nが解かれば振動変位成分faRの具体的数値を把握することができることが解かった。ただし、空力弾性振動体に加わる外力は時々刻々と変化する非定常空気力であるので、簡単に所定部位の空力減衰率δa及び応答円振動数nを把握することはできない。
そこで、空力減衰率δa及び応答円振動数nを把握する方法について、さらに検討し、その結果、振動体の所定部位に作用する各時刻の非定常空気力は、時々刻々と変化する所定部位の振幅y0(図7参照)に対応しているという考え方を用いることによって、空力弾性振動の応答を追跡できることを見出した。より具体的に言えば、空力減衰率δa及び応答円振動数nの具体的数値を、時々刻々と変化する所定部位の振幅y0に基づいて把握することができることを見出すと共に、各時刻における振幅y0を、当該時刻における変位及び変位速度に基づいて求めることができることを見出した。これを実際のシミュレーション方法の処理の流れに沿って説明すれば、各時刻における変位及び変位速度を用いて当該時刻の振幅y0を求め、当該振幅y0から空力減衰率δa及び応答円振動数nの具体的数値を把握し、これらを用いて非定常空気力の振動変位速度成分faI及び振動変位成分faRの具体的数値を把握し、把握した非定常空気力に関する具体的数値を用いてシミュレーションに用いる振動方程式を確定することができることを見出したということになる。
なお、時刻tにおける所定部位の振幅y0とは、各時刻tにおける変位ytではなく、その時刻tにおいて取り得た最大変位のことである。振幅y0は、図7の破線で示されるように変化する。
これらのことを見出したことで、最終的に、本願発明に係る空力弾性振動体の時刻歴応答シミュレーション方法に想到するに至った。
なお、時刻tにおける所定部位の振幅y0とは、各時刻tにおける変位ytではなく、その時刻tにおいて取り得た最大変位のことである。振幅y0は、図7の破線で示されるように変化する。
これらのことを見出したことで、最終的に、本願発明に係る空力弾性振動体の時刻歴応答シミュレーション方法に想到するに至った。
本発明では、まず、ある時刻の変位及び変位速度に基づいて当該時刻の振幅y0を求めるために、非定常空気力(外力)を受けた空力弾性振動体の所定部位の振動は、正弦波的な振動で減衰しあるいは発散するものであると考える。このように考えると、空力弾性振動の変位yを、下記の式(7)のように表すことができ、変位速度y’を、下記の式(8)のように表すことができる。そして、瞬間瞬間の振動に対応する振幅(正弦波振幅)y0を各時刻tにおける「変位」の自乗と「変位速度/振動数」の自乗の和の平方根で表現することができるということを見出して、式(7),(8)から式(9)を得た。さらに、橋梁の対風応答を検討する場合、減衰定数ξは0.01未満であるので、式(9)中のルートの中は、ほぼ「1」と考えることができ、式(9)中の「e」の乗数(-ξωt)はほぼ「0」と考えることができる。そうすると、空力弾性振動の式(9)から式(10)が得られる。
ここで、上記の式(7)から式(10)に関連して図7を参照する。図7において、実線は、空力弾性振動体の所定部位の時刻tにおける変位yを示すものである。そして、破線は、上記式(10)によって計算される時々刻々の振幅y0である。つまり、式(10)を用いることによって、弾性振動体の所定部位の各時刻における変位yt及び変位速度y’tの値から当該所定部位の振幅y0の当該時刻における値を計算することができる(第1ステップ)。
図7から解かるように、空力弾性振動体の所定部位の変位が振幅に一致するとき(変位速度が0のとき)は、式(10)がなくても、振幅y0を求めることができる。ところが、振幅y0が時々刻々と過渡応答的に変化する空力弾性振動の場合、式(10)がなければ、変位速度が0のとき以外の時刻について、振幅y0を求めることはできない。
空力弾性振動のような自励振動における空気力(非定常空気力)の評価は、振動振幅に対して行わざるを得ないものであるが、これまで、過渡応答的な振動における各瞬間の振幅y0を評価することができなかったので、振幅y0の変化に対応しつつ時々刻々の非定常空気力を評価することはできなかった。このようなことから、これまでは、時刻歴応答解析を用いて耐風特性を検討することはできず、空力弾性振動の時刻歴応答シミュレーションを行うことができなかった。
空力弾性振動のような自励振動における空気力(非定常空気力)の評価は、振動振幅に対して行わざるを得ないものであるが、これまで、過渡応答的な振動における各瞬間の振幅y0を評価することができなかったので、振幅y0の変化に対応しつつ時々刻々の非定常空気力を評価することはできなかった。このようなことから、これまでは、時刻歴応答解析を用いて耐風特性を検討することはできず、空力弾性振動の時刻歴応答シミュレーションを行うことができなかった。
この点、本発明によれば、上記式(10)を用いることによって、空力弾性振動体の所定部位の時刻tにおける変位yt及び変位速度y’tに基づき、時々刻々と変化する振幅y0を、いずれの時刻においても求めることができる(第1ステップ)。
なお、変位yの自乗値と、変位速度y’を応答円振動数nで除した値の自乗値との和の平方根の値に基づいて振幅y0を求める方法であれば、種々の方法によって振幅y0を求めることができるが、式(10)を用いる方法は、簡単であり好ましい。
なお、変位yの自乗値と、変位速度y’を応答円振動数nで除した値の自乗値との和の平方根の値に基づいて振幅y0を求める方法であれば、種々の方法によって振幅y0を求めることができるが、式(10)を用いる方法は、簡単であり好ましい。
次に、先に算出した振幅y0の時刻tにおける値を用いて、振幅y0の関数である空力減衰率δa及び応答円振動数nの時刻tにおける値を計算する(第2ステップ)。なお、振幅y0の関数である空力減衰率δa(=fI(y0))及び振幅y0の関数である応答円振動数n(=fR(y0))は、例えば弾性振動体に対して振動実験を行うことによって求められる。
続いて、先に算出された空力減衰率δaの時刻tにおける値を用いて前記非定常空気力の振動変位速度成分faIの時刻tにおける値を計算すると共に、先に算出された応答円振動数nの時刻tにおける値を用いて前記非定常空気力の振動変位成分faRの時刻tにおける値を計算する(第3ステップ)。
なお、空力減衰率δaの具体的数値を用いて非定常空気力(瞬時作用空気力)の振動変位速度成分faIを計算する方法は、種々考えられるかも知れないが、本発明において振動変位速度成分faIを計算する際に用いる空力減衰率δaの関数(式(16))は、次にようにして求められる。つまり、式(6)として示される振動方程式中の振動変位速度成分faIが減衰項の一部であり、振動方程式における定常振動の条件が次の式(11)であり、振幅−空力減衰率図における定常振動の条件が式(12)であることに基づき、求められる。
上記式(11)、(12)によって示される2つの関係が共通した現象における条件である。ここで、−faI/m及びδaは空気力に起因するものであり、2ξωがδsに対応(2ξω=δs)し、また、−faI/mがδaに対応する。そして、構造減衰率δs=2πξであるので、これから導かれる「ξ=δs/2π」を用いて「2ξω」を変形すると、次式(13)のようになり、式(11)及び(12)の関係から次の式(14)が導かれる。
したがって、式(6)の減衰項は、次の式(15)のようになり、δaを振動変位速度成分faIの式になるように変形すると、式(16)のようになる。これにより、二次元剛体モデルとして表される橋梁の振動方程式の振動変位速度成分faIの具体的な形が求められたことになる。また、あらかじめ得た空力減衰率δaを表す近似関数を用いて各時刻における空力減衰率δaを計算することができるので、計算した空力減衰率δaの値を式(16)に代入することで、その時刻の振幅y0に対応した振動変位速度成分faIが得られる。
また、応答円振動数nの具体的数値を用いて非定常空気力(瞬時作用空気力)の振動変位成分faRを計算する方法は、種々考えられるかも知れないが、本発明において振動変位成分faRを計算する際に用いられる応答円振動数nの関数(式(18))は、次にようにして求められる。つまり、式(6)として示される振動方程式中の振動変位成分faRが復元力項の一部であり、復元力項の係数が橋梁モデルの応答円振動数の2乗に等しいと考えると、次の式(17)が得られる。そして、この式を変形することで、次の式(18)が得られる。また、あらかじめ得た振動数nを表す近似関数を用いて各時刻における振動数nを計算することができるので、計算した振動数nの値を式(18)に代入することで、その時刻の振動変位成分faRが得られる。
このようにして求められた振動変位成分faR及び振動変位速度成分faIの時刻tにおける具体的数値を、振動変位成分faR及び振動変位速度成分faIを変数として含む弾性振動体の振動方程式に代入して得られた振動方程式を解いて、次の時刻t+Δtにおける前記所定部位の変位yt+Δt、変位速度y’t+Δt及び加速度y’’t+Δtの値を計算する(第4ステップ)。これにより、時刻歴応答シミュレーションに用いられる数値が計算される。
本願発明に係る空力弾性振動の時刻歴応答解析シミュレーション方法では、これらの第1ステップから第4ステップを含む時刻歴応答計算ルーチンを各時刻について繰り返し実行し、計算された値を基に弾性振動体の時刻歴応答シミュレーションを行うものである。なお、2回目以降の時刻歴応答計算ルーチンの第1ステップでは、当該第1ステップにおける変位yt及び変位速度y’tの値として、直前の時刻歴応答計算ルーチンの第4ステップ(次の時刻t+Δtにおける各値を計算するステップ)で計算された変位y及び変位速度y’の値を用いて振幅y0の対応時刻における値を計算する。また、2回目以降の時刻歴応答計算ルーチンの第4ステップは、当該第4ステップにおける変位yt、変位速度y’t及び加速度y’’tの値として、直前の時刻歴応答計算ルーチンの第4ステップ(次の時刻t+Δtにおける各値を計算するステップ)で計算された変位y、変位速度y’及び加速度y’’の値を用いて、次の時刻t+Δtにおける前記所定部位の変位yt+Δt、変位速度y’t+Δt及び加速度y’’t+Δtの値を計算する。
このように、本願発明によれば、減衰状態あるいは発散状態にある所定部位の空力弾性振動について、時々刻々の振幅y0を計算することができ、計算された振幅y0を用いて非定常空気力(より具体的にはの振動変位成分faR及び振動変位速度成分faI)を計算することができる。時々刻々の非定常空気力を計算することができれば、計算された非定常空気力を振動方程式に代入することによって時々刻々の振動方程式を定め、定められた振動方程式を解くことによって、空力弾性振動の時刻歴応答解析シミュレーションをすることができる。
なお、前記振幅y0の関数である空力減衰率δa(=fI(y0))及び応答円振動数n(=fR(y0))は、前記空力弾性振動体について行う振動実験によって求めたものが好ましい。例えば、シミュレーションの対象の空力弾性振動体について自由振動実験等の振動実験を行うことによって、上記関数(δa=fI(y0))及び関数(n=fR(y0))を求めることができる。
このように、シミュレーション対象の弾性振動体について実際に実験を行うことによって、上記関数である空力減衰率δa(=fI(y0))及び応答円振動数n(=fR(y0))を正確に求めることができる。そして、このような関数があれば、時々刻々と変化する空力減衰率や応答円振動数を、迅速、簡単かつ確実に求めることができることとなり、シミュレーションをより迅速、簡単かつ確実に実行することができる。
また、前記空力減衰率δa(=fI(y0))及び応答円振動数n(=fR(y0))は、前記振幅y0及び風速の関数であることがより好ましい。振幅y0だけでなく風速が加味された関数を用いることにより、より正確なシミュレーションを行うことができる。
また、本願の別の発明は、外力として非定常空気力を受ける弾性振動体の時刻歴応答シミュレーションをするためのコンピュータを、シミュレーションの対象である前記弾性振動体の所定部位の条件を設定する手段と、風速条件を設定する手段、
前記弾性振動体の所定部位の初期時刻における変位y及び変位速度y’の値を設定する手段、前記弾性振動体の所定部位の時刻tにおける変位yt及び変位速度y’tの値を用いて当該所定部位の振幅y0の時刻tにおける値を計算する手段、
前記風速条件に基づいて特定された時刻tにおける風速と、振幅y0の時刻tにおける値に基づいて空力減衰率δa及び応答円振動数nの時刻tにおける値を計算する手段、空力減衰率δaの時刻tにおける値を用いて前記非定常空気力の振動変位速度成分faIの時刻tにおける値を計算すると共に、応答円振動数nの時刻tにおける値を用いて前記非定常空気力の振動変位成分faRの時刻tにおける値を計算する手段、振動変位成分faR及び振動変位速度成分faIを変数として含む弾性振動体の振動方程式に、振動変位速度成分faI及び振動変位成分faRの時刻tにおける値を代入して得られた振動方程式を解いて、次の時刻t+Δtにおける前記所定部位の変位yt+Δt、変位速度y’t+Δt及び加速度y’’t+Δtの値を計算する手段、前記計算する手段を順次実行する時刻歴応答計算ルーチンを各時刻について繰り返し実行させる手段であって、2回目以降の時刻歴応答計算ルーチンの振幅y0の値の計算で用いられる変位yt及び変位速度y’tの値として、直前の時刻歴応答計算ルーチンで計算された変位y及び変位速度y’の値を用いて振幅y0の値を計算させる手段、前記時刻歴応答計算ルーチンを各時刻について繰り返し実行して計算された値を用いて弾性振動体の時刻歴応答シミュレーションを出力する手段として機能させる時刻歴応答シミュレーションプログラムである。このようなプログラムによって上述したような時刻歴応答シミュレーションを実現することができ、当該プログラムがインストールされたコンピュータなどの時刻歴応答シミュレーション装置によって上述したような時刻歴応答シミュレーションを実現することができる。
前記弾性振動体の所定部位の初期時刻における変位y及び変位速度y’の値を設定する手段、前記弾性振動体の所定部位の時刻tにおける変位yt及び変位速度y’tの値を用いて当該所定部位の振幅y0の時刻tにおける値を計算する手段、
前記風速条件に基づいて特定された時刻tにおける風速と、振幅y0の時刻tにおける値に基づいて空力減衰率δa及び応答円振動数nの時刻tにおける値を計算する手段、空力減衰率δaの時刻tにおける値を用いて前記非定常空気力の振動変位速度成分faIの時刻tにおける値を計算すると共に、応答円振動数nの時刻tにおける値を用いて前記非定常空気力の振動変位成分faRの時刻tにおける値を計算する手段、振動変位成分faR及び振動変位速度成分faIを変数として含む弾性振動体の振動方程式に、振動変位速度成分faI及び振動変位成分faRの時刻tにおける値を代入して得られた振動方程式を解いて、次の時刻t+Δtにおける前記所定部位の変位yt+Δt、変位速度y’t+Δt及び加速度y’’t+Δtの値を計算する手段、前記計算する手段を順次実行する時刻歴応答計算ルーチンを各時刻について繰り返し実行させる手段であって、2回目以降の時刻歴応答計算ルーチンの振幅y0の値の計算で用いられる変位yt及び変位速度y’tの値として、直前の時刻歴応答計算ルーチンで計算された変位y及び変位速度y’の値を用いて振幅y0の値を計算させる手段、前記時刻歴応答計算ルーチンを各時刻について繰り返し実行して計算された値を用いて弾性振動体の時刻歴応答シミュレーションを出力する手段として機能させる時刻歴応答シミュレーションプログラムである。このようなプログラムによって上述したような時刻歴応答シミュレーションを実現することができ、当該プログラムがインストールされたコンピュータなどの時刻歴応答シミュレーション装置によって上述したような時刻歴応答シミュレーションを実現することができる。
そして、本願に係るさらに別の発明は、外力として非定常空気力を受ける弾性振動体の時刻歴応答シミュレーションを行う時刻歴応答シミュレーション装置であって、非定常空気力を受ける前記弾性振動体の所定部位の時刻tにおける変位yt及び変位速度y’tの値を用いて前記所定部位の振幅y0の時刻tにおける値を計算する第1計算手段と、振幅y0の時刻tにおける値を用いて、振幅y0の関数である空力減衰率δa及び応答円振動数nの時刻tにおける値を計算する第2計算手段と、空力減衰率δaの時刻tにおける値を用いて前記非定常空気力の振動変位速度成分faIの時刻tにおける値を計算すると共に、応答円振動数nの時刻tにおける値を用いて前記非定常空気力の振動変位成分faRの時刻tにおける値を計算する第3計算手段と、振動変位成分faR及び振動変位速度成分faIを変数として含む弾性振動体の振動方程式に、振動変位速度成分faI及び振動変位成分faRの時刻tにおける値を代入して得られた振動方程式を解いて、次の時刻t+Δtにおける前記所定部位の変位yt+Δt、変位速度y’t+Δt及び加速度y’’t+Δtの値を計算する第4計算手段と、前記計算する手段を順次実行する時刻歴応答計算ルーチンを各時刻について繰り返し実行させる手段であって、2回目以降の時刻歴応答計算ルーチンの振幅y0の値の計算で用いられる変位yt及び変位速度y’tの値として、直前の時刻歴応答計算ルーチンで計算された変位y及び変位速度y’値を用いて振幅y0の値を計算させる手段と、前記時刻歴応答計算ルーチンを各時刻について繰り返し実行して計算された値を用いて弾性振動体の時刻歴応答シミュレーションを出力する手段と、を備えるものである。このような時刻歴応答シミュレーション装置によって上述したような時刻歴応答シミュレーションを実現することができる。
本願発明に係る空力弾性振動体の時刻歴応答シミュレーション方法、そのプログラム及び装置によれば、橋梁などの空力弾性振動体の所定部位の空力減衰率δaを表す関数(δa=f(y0))と、応答円振動数nを表す関数(n=f(y0))とをあらかじめ求めておくことにより、非定常空気力(外力)を受ける当該構造物の耐風特性(例えば、空力弾性振動)の検討を、時刻歴応答解析シミュレーションを用いて行うことができる。
そして、風洞実験を行うことなくシミュレーションによって空力弾性振動の応答推定を行うことができるので、風洞実験を行うことが難しい長大橋梁等の構造物について、空力弾性振動の応答推定を容易、迅速に行うことができる。
また、本願発明によれば、構造物の構造様式などの条件や構造物が受ける風(風速や風速の変動)など、シミュレーションのパラメータの設定を容易に変更することができるので、多様な条件下のシミュレーションを容易、迅速に行うことができる。これにより、多様な条件下での構造物の空力弾性振動の応答推定を、容易、迅速に行うことができる。
さらに、本願発明によれば、時刻歴応答シミュレーションを用いるので、耐震設計と同等の精度で、耐風設計を行うことが可能である。つまり、構造物の空力弾性振動の応答推定を行う際に、本願発明に係る空力弾性振動体の時刻歴応答シミュレーション方法を利用することによって、より正確な耐風設計を行うことが可能である。これにより、構造物の耐風特性に関する設計の品質を飛躍的に向上させることができる。
また、時刻tにおける変位ytの自乗値と、時刻tにおける変位速度y’tを応答円振動数nで除した値の自乗値との和の平方根を求めることによって時刻tにおける所定部位の振幅y0tの値を求めることができれば、より容易、迅速かつ正確に時刻歴応答シミュレーションを行うことができる。
また、シミュレーションの各ステップでの計算に用いることができる具体的な関数が用意されているので、それらの関数を用いることによって、容易、迅速かつ正確にシミュレーションを行うことができる。
また、シミュレーション対象の弾性振動体については、実際に実験を行うことによって、空力減衰率δaを表す振幅y0の関数(δa=fI(y0))及び応答円振動数nを表す振幅y0の関数(n=fR(y0))を、シミュレーションを正確に求めることができる。そして、このような関数があれば、時々刻々と変化する空力減衰率や応答円振動数を、迅速、簡単かつ確実に求めることができることとなり、シミュレーションをより迅速、簡単かつ確実に実行することができる。
また、前記空力減衰率δaを表す振幅y0の関数(δa=fI(y0))及び前記応答円振動数nを表す振幅y0の関数(n=fR(y0))を風速ごとに求めたものを用いることにより、正確なシミュレーションを行うことができる。
本願の発明の時刻歴応答シミュレーション方法は、外力として非定常空気力を受ける弾性振動体の時刻歴応答シミュレーション方法であって、非定常空気力を受ける前記弾性振動体の所定部位の時刻tにおける変位yt及び変位速度y’tの値を用いて前記所定部位の振幅y0の時刻tにおける値を計算する第1ステップと、振幅y0の時刻tにおける値を用いて、振幅y0の関数である空力減衰率δa及び応答円振動数nの時刻tにおける値を計算する第2ステップと、空力減衰率δaの時刻tにおける値を用いて前記非定常空気力の振動変位速度成分faIの時刻tにおける値を計算すると共に、応答円振動数nの時刻tにおける値を用いて前記非定常空気力の振動変位成分faRの時刻tにおける値を計算する第3ステップと、振動変位成分faR及び振動変位速度成分faIを変数として含む弾性振動体の振動方程式に振動変位速度成分faI及び振動変位成分faRの時刻tにおける値を代入して得られた振動方程式を解いて、次の時刻t+Δtにおける前記所定部位の変位yt+Δt、変位速度y’t+Δt及び加速度y’’t+Δtの値を計算する第4ステップと、を備え、これらのステップを含む時刻歴応答計算ルーチンを各時刻について繰り返し実行して計算された値を基に弾性振動体の時刻歴応答シミュレーションを行うものであり、2回目以降の時刻歴応答計算ルーチンの第1ステップでは、当該第1ステップにおける変位yt及び変位速度y’tの値として、直前の時刻歴応答計算ルーチンの第4ステップで計算された変位y及び変位速度y’の値を用いて振幅y0の対応時刻における値を計算するものである。
そして、2回目以降の時刻歴応答計算ルーチンの第4ステップは、当該第4ステップにおける変位yt、変位速度y’t及び加速度y’’tの値として、直前の時刻歴応答計算ルーチンの第4ステップで計算された変位y、変位速度y’及び加速度y’’の値を用いて、次の時刻t+Δtにおける前記所定部位の変位yt+Δt、変位速度y’t+Δt及び加速度y’’t+Δtの値を計算するものである。
また、前記所定部位の振幅y0の時刻tにおける値を求めるための前記変位y及び変位速度y’の関数は、時刻tにおける変位ytの自乗値と、時刻tにおける変位速度y’tを応答円振動数nで除した値の自乗値との和の平方根を振幅y0の時刻tにおける値として求めることができるものである。
また、前記振動方程式は、下記の式(6)であり、前記第4ステップで振動変位速度成分faI及び振動変位成分faRの値が代入された振動方程式の解法は直接時間積分法であり、前記第1ステップで前記振幅y0の値を計算するために用いられる変位y及び変位速度y’の関数は、下記の式(10)であり、前記第3ステップで前記振動変位速度成分faIの値を計算するために用いられる前記空力減衰率δaの関数は、下記の式(16)であり、前記第3ステップで前記振動変位成分faRの値を計算するために用いられる前記応答円振動数nの関数は、下記の式(18)であることが好ましい。
また、前記空力減衰率δaを求めるための前記振幅y0の関数(δa=fI)及び前記応答円振動数nを求めるための前記振幅y0の関数(n=fR)は、前記弾性振動体について行う振動実験によって求められるものであってもよい。さらに、前記関数(δa=fI)及び前記関数(n=fR)は、前記振幅y0及び風速の関数であってもよい。
また、本願に係る別の発明は、外力として非定常空気力を受ける弾性振動体の時刻歴応答シミュレーションをするためのコンピュータを、シミュレーションの対象である前記弾性振動体の所定部位の条件を設定する手段と、風速条件を設定する手段、前記弾性振動体の所定部位の初期時刻における変位y及び変位速度y’の値を設定する手段、前記弾性振動体の所定部位の時刻tにおける変位yt及び変位速度y’tの値を用いて当該所定部位の振幅y0の時刻tにおける値を計算する手段、前記風速条件に基づいて特定された時刻tにおける風速と、振幅y0の時刻tにおける値に基づいて空力減衰率δa及び応答円振動数nの時刻tにおける値を計算する手段、空力減衰率δaの時刻tにおける値を用いて前記非定常空気力の振動変位速度成分faIの時刻tにおける値を計算すると共に、応答円振動数nの時刻tにおける値を用いて前記非定常空気力の振動変位成分faRの時刻tにおける値を計算する手段、振動変位成分faR及び振動変位速度成分faIを変数として含む弾性振動体の振動方程式に、振動変位速度成分faI及び振動変位成分faRの時刻tにおける値を代入して得られた振動方程式を解いて、次の時刻t+Δtにおける前記所定部位の変位yt+Δt、変位速度y’t+Δt及び加速度y’’t+Δtの値を計算する手段、前記計算する手段を順次実行する時刻歴応答計算ルーチンを各時刻について繰り返し実行させる手段であって、2回目以降の時刻歴応答計算ルーチンの振幅y0の値の計算で用いられる変位yt及び変位速度y’tの値として、直前の時刻歴応答計算ルーチンで計算された変位y及び変位速度y’の値を用いて振幅y0の値を計算させる手段、前記時刻歴応答計算ルーチンを各時刻について繰り返し実行して計算された値を用いて弾性振動体の時刻歴応答シミュレーションを出力する手段、として機能させる時刻歴応答シミュレーションプログラムである。
また、本願に係るさらに別の発明は、外力として非定常空気力を受ける弾性振動体の時刻歴応答シミュレーションを行う時刻歴応答シミュレーション装置であって、非定常空気力を受ける前記弾性振動体の所定部位の時刻tにおける変位yt及び変位速度y’tの値を用いて前記所定部位の振幅y0の時刻tにおける値を計算する第1計算手段と、振幅y0の時刻tにおける値を用いて、振幅y0の関数である空力減衰率δa及び応答円振動数nの時刻tにおける値を計算する第2計算手段と、空力減衰率δaの時刻tにおける値を用いて前記非定常空気力の振動変位速度成分faIの時刻tにおける値を計算すると共に、応答円振動数nの時刻tにおける値を用いて前記非定常空気力の振動変位成分faRの時刻tにおける値を計算する第3計算手段と、振動変位成分faR及び振動変位速度成分faIを変数として含む弾性振動体の振動方程式に、振動変位速度成分faI及び振動変位成分faRの時刻tにおける値を代入して得られた振動方程式を解いて、次の時刻t+Δtにおける前記所定部位の変位yt+Δt、変位速度y’t+Δt及び加速度y’’t+Δtの値を計算する第4計算手段と、前記計算する手段を順次実行する時刻歴応答計算ルーチンを各時刻について繰り返し実行させる手段であって、2回目以降の時刻歴応答計算ルーチンの振幅y0の値の計算で用いられる変位yt及び変位速度y’tの値として、直前の時刻歴応答計算ルーチンで計算された変位y及び変位速度y’の値を用いて振幅y0の値を計算させる手段と、前記時刻歴応答計算ルーチンを各時刻について繰り返し実行して計算された値を用いて弾性振動体の時刻歴応答シミュレーションを出力する手段と、を備える時刻歴応答シミュレーション装置である。
以下、本発明に係る時刻歴応答シミュレーション方法、そのプログラム及び装置の実施例について、図面を参照しつつ詳細に説明する。ここでは、図1及び図2に示される箱桁断面の橋梁模型(空力弾性振動体)を例に、バネ吊り1自由度振動系の空力弾性応答のシミュレーションについて説明する。
まず、シミュレーションの前にあらかじめ、図1に示される形状及び寸法(mm)の橋梁模型の桁部(所定部位)について自由振動実験を行った。より具体的には、橋梁模型の桁部の一部に該当する、図2に示される二次元部分剛体模型をバネで吊った状態で風洞内に設置し、風を吹かせて応答を測定するという自由振動実験を行った(図3参照)。自由振動実験に用いた二次元部分剛体模型は、実物の1/65の縮尺のものであり、長さが800mm、質量mが1.56kg、構造減衰率δsが0.0019229、固有振動数ωが7.756Hzというものであった。自由振動実験では、風速を徐々に変化させ、各風速における発散振動状況、減衰振動状況の振動波形をそのつど記録し、記録した振動波形から各風速における振幅y0ごとの減衰率δを算出した。ここでは風速0.2m/sごとに上記の記録を行った。そして、算出した減衰率δから構造減衰率δsを差し引いて、振幅y0の関数である空力減衰率δa(=fI(y0))を求めた。また、記録した振動波形を基に、振幅y0の関数である後述の応答円振動数n(=fR(y0))を求めた。所定部位とは、橋梁の桁部、主塔部など、それぞれ固有の断面形状を有する全ての部分のことである。ここでは、連続高架橋の桁部という部位の時刻歴応答シミュレーションについて説明する。なお、必要に応じて各部位について自由振動実験を行い、各部位の空力減衰率δa、応答円振動数nを算出することになるが、その手順はここで説明するものと同様であるので、説明を省略する。
ここで、実験によって求めたデータを基に作成した1自由度たわみ応答(二次元模型の対風振動応答図)を示す(図4の実験値のグラフ参照)。この図では、横軸が換算風速Vr、縦軸が無次元化した振幅(2A/D)である。
また、自由振動実験を行って求めた空力減衰率δaのグラフの一例として、渦励振時のピークを示す換算風速Vr=5.76における空力減衰率δaの応答図を示す(図5の実験値のグラフ参照)。そして、関数近似を行って、次の式(21)で示されるような空力減衰率δaの近似式(近似関数)を得た。
下記の式から解かるように、ここでは、振幅値について区間を複数設定し(実施例は3区間)、各区間ごとに近似式を求めた。このようにして、振幅値の値に応じた近似式を用いて空力減衰率δaを求めるようにすると、より正確に空力減衰率δaを求めることができ、より正確なシミュレーションを行うことができる。また、この点は、後述の応答円振動数についても同様である。なお、他の換算風速においても同様に、必要に応じて区間を区切って空力減衰率δaの近似式(近似関数)を得たが、手順等は同様であるので、ここではその詳細な説明を省略する。
なお、近似化に用いた数値を直接用いてシミュレーションを行ってもよいが、シミュレーションのためには、関数化しておく方が都合が良く、好ましい。これは、空力減衰率についてだけでなく、後述の応答円振動数nについても同様である。
また、自由振動実験を行って求めた空力減衰率δaのグラフの一例として、渦励振時のピークを示す換算風速Vr=5.76における空力減衰率δaの応答図を示す(図5の実験値のグラフ参照)。そして、関数近似を行って、次の式(21)で示されるような空力減衰率δaの近似式(近似関数)を得た。
下記の式から解かるように、ここでは、振幅値について区間を複数設定し(実施例は3区間)、各区間ごとに近似式を求めた。このようにして、振幅値の値に応じた近似式を用いて空力減衰率δaを求めるようにすると、より正確に空力減衰率δaを求めることができ、より正確なシミュレーションを行うことができる。また、この点は、後述の応答円振動数についても同様である。なお、他の換算風速においても同様に、必要に応じて区間を区切って空力減衰率δaの近似式(近似関数)を得たが、手順等は同様であるので、ここではその詳細な説明を省略する。
なお、近似化に用いた数値を直接用いてシミュレーションを行ってもよいが、シミュレーションのためには、関数化しておく方が都合が良く、好ましい。これは、空力減衰率についてだけでなく、後述の応答円振動数nについても同様である。
また、空力減衰率δaを求める際に、同時に各風速における応答円振動数nを求めた。一例として換算風速Vr=5.76に対する応答円振動数nの応答図を示す(図6の実験値のグラフ参照)。そして、関数近似を行って、次の式(22)で示されるような応答円振動数nの近似式を得た。下記の式から解かるように、空力減衰率δaの場合と同様、振幅値について区間を複数設定し(ここでは2区間)、各区間ごとに近似式を求めた。この近似式を用いて、いわゆる非定常空気力の変位比例成分を得ることができる。なお、他の換算風速においても同様に、必要に応じて区間を区切って応答円振動数nの近似式を得たが、手順等は同様であるので、ここではその詳細な説明を省略する。
このようにして、各風速および振幅に対応した橋梁の所定部位の空力減衰率δaを表す近似式(近似関数)及び応答円振動数nを表す近似式(近似関数)を求めた上で、次に説明しいている各ステップから構成される時刻歴応答計算ルーチンを時々刻々と変化する各時刻について繰り返し実行し、時刻歴応答シミュレーションを行った。
以下、時刻歴応答シミュレーション方法、そのプログラム及び装置について説明する。
図9に示されるように、本シミュレーション方法を実行するシミュレーション装置10は、CPU11等の演算プロセッサやメインメモリ12等のメモリを備え、OSが例えばWINDOWS(登録商標)XPエディションであるパソコン本体13と、パソコン本体13に入出力インターフェース14を介して接続された出力手段であるディスプレイ15と、プリンタ(不図示)と、外部記憶装置であるハードディスク16と、CD−ROMの記録を読み取るためのCD−ROMドライブ17(または各種記憶媒体のためのマルチドライブ)と、パソコン本体13に入出力インターフェース14を介して接続された入力手段であるマウス18やキーボード19とを備える周知のコンピュータである。なお、符号「20」はバスである(図9参照)。
そして、本発明にかかる時刻歴応答シミュレーションプログラムがハードディスク15又はマルチドライブにより読み取り可能なCD−ROM等の記憶媒体に記憶されており、本シミュレーション装置10にインストールされている。そして、主にCPU11やメインメモリ12が計算などのステップ処理部や出力するデータを生成する手段として機能し、メインメモリ12が入力された設定条件の記憶手段として機能する。また、計算によって算出された数値データや生成された画像データ等がハードディスク16やディスプレイ15などに出力される。なお、シミュレーション装置10として用いられているコンピュータ自体は周知の装置であるので、ここでは、その詳細な説明を省略する。
そして、本発明にかかる時刻歴応答シミュレーションプログラムがハードディスク15又はマルチドライブにより読み取り可能なCD−ROM等の記憶媒体に記憶されており、本シミュレーション装置10にインストールされている。そして、主にCPU11やメインメモリ12が計算などのステップ処理部や出力するデータを生成する手段として機能し、メインメモリ12が入力された設定条件の記憶手段として機能する。また、計算によって算出された数値データや生成された画像データ等がハードディスク16やディスプレイ15などに出力される。なお、シミュレーション装置10として用いられているコンピュータ自体は周知の装置であるので、ここでは、その詳細な説明を省略する。
次に、時刻歴応答シミュレーション方法について説明する。本シミュレーション方法は、上述したように、時刻tにおける所定部位の振幅y0tに対応する空気力が所定部位に瞬間瞬間作用しているという考え方を用いて空力弾性振動の応答を追跡し、シミュレーションするものである。
なお、時刻tにおける所定部位の振幅y0とは、各時刻tにおける所定部位の実際の変位y(図7の実線参照)ではなく、時刻tにおいて取り得た最大変位(正弦波振幅、図7の破線参照)のことである。別言すれば、時刻tにおける所定部位の時々刻々変化する振幅のことである。この振幅y0は、図7の破線で示されるように変化する。
なお、時刻tにおける所定部位の振幅y0とは、各時刻tにおける所定部位の実際の変位y(図7の実線参照)ではなく、時刻tにおいて取り得た最大変位(正弦波振幅、図7の破線参照)のことである。別言すれば、時刻tにおける所定部位の時々刻々変化する振幅のことである。この振幅y0は、図7の破線で示されるように変化する。
まず、シミュレーションの対象である弾性振動体の所定部位の条件を設定する。本実施例では、所定部位として桁部(図1参照)が設定され、桁部の条件として上述したような質量等の条件が設定される。そして、シミュレーションしようとしている状況に応じた風速条件(換算風速Vrの条件)を設定する。外力として非定常空気力を受ける状態を検討する場合は、時々刻々と変化する条件で設定がなされる。各時刻における具体的風速値(および/または換算風速Vr)をあらかじめ設定しておき必要に応じて読み取って用いることができるようにしてもよいし、各時刻における風速を算出するための関数を用意しておき必要に応じて各時刻における風速を算出して用いることができるようにしてもよい。これにより、シミュレーションの対象である橋梁の所定部位の所定時刻における風速(換算風速Vr)を定めることができる。また、最初の時刻における所定部位の変位、変位速度及び変位加速度を設定する。なお、各種条件設定については、キーボード等の入力手段を用いた入力や、あらかじめ記憶させておいた種々の条件の中からの選択など、種々の方法によって行うことができる。各種設定条件は、必要に応じて記憶手段に記憶されて用いられる。
そして、インストールされたプログラムに基づいて、以降に説明するような計算を行う。
まず、所定時刻t1(最初の時刻歴応答計算ルーチンでは最初の時刻)における風速を特定する。具体的には、先に設定した風速条件に基づいて、風速または換算風速Vrの具体的値が特定される。また、時刻t1の所定部位の変位yt1及び変位速度y’t1を式(10)中の変位y及び変位速度y’に代入して振幅y0の計算する。これにより、非定常空気力(外力)を受ける橋梁の桁部の振幅y0の時刻t1における値が求められる(第1ステップ)。時刻歴応答シミュレーションプログラムは、このような計算をする計算ルーチンを備えている(第1計算手段)。
まず、所定時刻t1(最初の時刻歴応答計算ルーチンでは最初の時刻)における風速を特定する。具体的には、先に設定した風速条件に基づいて、風速または換算風速Vrの具体的値が特定される。また、時刻t1の所定部位の変位yt1及び変位速度y’t1を式(10)中の変位y及び変位速度y’に代入して振幅y0の計算する。これにより、非定常空気力(外力)を受ける橋梁の桁部の振幅y0の時刻t1における値が求められる(第1ステップ)。時刻歴応答シミュレーションプログラムは、このような計算をする計算ルーチンを備えている(第1計算手段)。
なお、変位速度が0のときは、所定部位の実際の変位は振幅y0に一致するので、例えば時刻t1における変位速度が0である場合は、式(10)を用いずに、実際の変位yt1を時刻t1の振幅y0としてもよい。
次に、空力減衰率δa及び応答円振動数nの値を計算する(第2ステップ)。まず、第1ステップで定められた風速条件及び振幅y0の値に応じて、ここでの計算で用いる空力減衰率δaの近似式および応答円振動数nの近似式を特定する。例えば、風速が「換算風速Vr=5.76」に対する風速であり、振幅y0が0.03であれば、空力減衰率δaの計算に用いられる近似式は、上記式(21)中の近似式2に特定され、応答円振動数nの計算に用いられる近似式は、上記の式(22)中の近似式1に特定される。なお、完全に一致する条件に対応した近似式がない場合は、最も近い条件に対応した近似式を、用いる近似式として特定してもよい。このようにして特定された式に、時刻t1における振幅y0t1の値を代入して、当該時刻t1における所定部位の空力減衰率δa及び応答円振動数nの値を計算する(第2ステップ)。時刻歴応答シミュレーションプログラムは、このような計算をする計算ルーチンを備えている(第2計算手段)。
続いて、時刻t1の時点で所定部位に加わっている外力を求める。具体的には、第2ステップで求めた空力減衰率δaの具体的数値を式(16)に代入して振動変位速度成分faIの具体的数値を計算し、第2ステップで求めた応答円振動数nの具体的数値を式(18)に代入して振動変位成分faRの具体的数値を計算する(第3ステップ)。時刻歴応答シミュレーションプログラムは、このような計算をする計算ルーチンを備えている(第3計算手段)。
そして、計算された振動変位速度成分faI及び振動変位成分faRの具体的数値を式(6)に代入することによって、時刻t1における振動方程式が特定される。
この特定された振動方程式を、直接時間積分法で解いて、時間刻みに応じた振幅の時刻歴を計算する(第4ステップ)。時刻歴応答シミュレーションプログラムは、このような計算をする計算ルーチンを備えている(第4計算手段)。
本実施例では、不規則外力を受ける構造物の応答計算を数値計算によって行う際に用いられる直接時間積分法の1つであるニューマークβ法を用いて、振幅の時刻歴を計算する。
ニューマークβ法は、不規則外力を受ける構造物の応答計算を行う際に用いられる数値計算方法の1つである。この方法を用い、振動方程式の外力項を空気力とすることで、時系列での応答推定を行う。
ニューマークβ法は、時刻t+Δt(=例えば時刻t2)のときの所定部位の実際の変位を求めるために、時刻t(例えば時刻t1)のときの値とΔt時間後である時刻t+Δt(例えばt2)のときの値を使用して繰り返し計算を行う解法である。使用する方程式は、次の式(24)である。
ニューマークβ法は、不規則外力を受ける構造物の応答計算を行う際に用いられる数値計算方法の1つである。この方法を用い、振動方程式の外力項を空気力とすることで、時系列での応答推定を行う。
ニューマークβ法は、時刻t+Δt(=例えば時刻t2)のときの所定部位の実際の変位を求めるために、時刻t(例えば時刻t1)のときの値とΔt時間後である時刻t+Δt(例えばt2)のときの値を使用して繰り返し計算を行う解法である。使用する方程式は、次の式(24)である。
つまり、第4ステップでは、上記式(26)を用いて、時刻t2(=t+Δt)における前記所定部位の変位yt+Δt、変位速度y’t+Δt及び変位加速度y’’t+Δtを計算することができる(第4ステップ、第4計算手段)。
これにより、最初の時刻歴応答計算ルーチンR1において、所定部位の時刻t2(=t+Δt)における実際に変位が求められたことになる。
これにより、最初の時刻歴応答計算ルーチンR1において、所定部位の時刻t2(=t+Δt)における実際に変位が求められたことになる。
また、シミュレーションプログラムは、計算された変位yt2、変位速度y’t2及び変位加速度y’’t2の値を、次の時刻t3における所定部位の変位等を求める次の時刻歴応答計算ルーチンR2において、変位yt及び変位速度y’tの値として用いて計算を行わせるようになっている。例えば、計算された時刻t2(=t+Δt)の変位yの値及び変位速度y’の値は、時刻歴応答計算ルーチンR2の第1ステップにおいて、式(10)の変位y及び変位速度y’の値として代入され、これにより時刻t3(=t+2×Δt)における所定部位の振幅y0が計算される(第1ステップ)。また、計算された変位yt2、変位速度y’t2及び変位加速度y’’t2の値は、時刻歴応答計算ルーチンR2の第4ステップにおいて、式(26)の変位yt、変位速度y’t及び変位加速度y’’tの値として代入され、これにより、時刻t3における所定部位の変位yt+Δt、変位速度y’t+Δt及び変位加速度y’’t+Δtが計算される(第4ステップ)。
このように、2回目以降のx回目の時刻歴応答計算ルーチンRxの第1ステップでは、直前の時刻歴応答計算ルーチンRx-1の第4ステップで計算された変位ytx及び変位速度y’txの値を、変位y及び変位速度y’の値として用いて(式(10)等参照)、次の時刻tx+1における振幅y0を計算する。
また、2回目以降の時刻歴応答計算ルーチンRxの第4ステップでは、直前の時刻歴応答計算ルーチンRx-1の第4ステップで計算された変位ytx、変位速度y’tx及び変位加速度y’’txの値を、当該第4ステップにおける変位yt、変位速度y’t及び変位加速度y’’tの値として用いて次の時刻tx+1における変位y、変位速度y’及び変位加速度y’’の値を計算する。
時刻歴応答シミュレーションプログラムは、このように時刻歴応答計算ルーチンを各時刻について繰り返し実行させる計算ルーチンであって、2回目以降の時刻歴応答計算ルーチンの各計算で用いられる変位yt、変位速度y’t及び変位加速度y’’の値として、直前の時刻歴応答計算ルーチンで計算された変位y、変位速度y’及び変位加速度y’’の値を用いて振幅y0の値を計算させる計算ルーチンを、計算手段として備えている。
また、2回目以降の時刻歴応答計算ルーチンRxの第4ステップでは、直前の時刻歴応答計算ルーチンRx-1の第4ステップで計算された変位ytx、変位速度y’tx及び変位加速度y’’txの値を、当該第4ステップにおける変位yt、変位速度y’t及び変位加速度y’’tの値として用いて次の時刻tx+1における変位y、変位速度y’及び変位加速度y’’の値を計算する。
時刻歴応答シミュレーションプログラムは、このように時刻歴応答計算ルーチンを各時刻について繰り返し実行させる計算ルーチンであって、2回目以降の時刻歴応答計算ルーチンの各計算で用いられる変位yt、変位速度y’t及び変位加速度y’’の値として、直前の時刻歴応答計算ルーチンで計算された変位y、変位速度y’及び変位加速度y’’の値を用いて振幅y0の値を計算させる計算ルーチンを、計算手段として備えている。
このようなステップからなる時刻歴応答計算ルーチンを各時刻について繰り返し実行することによって、各時刻における応答(実際の変位等)を計算し、計算された変位等の値を用いて、Δt時間刻みの時刻歴応答シミュレーションを行う。シミュレーションの結果は、ディスプレイ15やプリンタなどに出力される。ハードディスク16などの記憶媒体に記録される態様で出力することもできる。時刻歴応答シミュレーションプログラムは、このような機器に時刻歴応答シミュレーションの結果を出力させる処理ルーチンを備えている。
そして、上記のようなシミュレーションを行った結果、図8(a)に示されるような時刻歴応答解析波形が得られた。また、比較のために、実物の1/130の縮尺の全橋模型について、シミュレーションの条件と同じ条件で風洞実験を行って、時刻歴応答解析波形を得た。図8(b)に、風洞実験での振幅の測定値に基づいて作成した振幅波形を示す。両波形を比較するとわかるように、両波形はほぼ一致していた。
そして、図5中の推定値のグラフは、図8(a)に示される時刻歴応答解析波形を基に算出した空力減衰率の推定値のグラフである。図5から明らかなように、空力減衰率の近似値及び推定値は、いずれも、実験値とほぼ一致していた。
また、時刻歴解析を各風速で解析し、定常振動応答をプロットして、測定値との比較を行ったものが図4である。この図に示されるように、時刻歴応答解析から求めた定常振動応答推定値は、若干のずれはあるものの実験での測定値をほぼ再現していた。
このように、本願発明によれば、減衰状態あるいは発散状態にある所定部位の空力弾性振動について、時々刻々の振幅y0を計算することができ、計算された振幅y0を用いて非定常空気力を計算することができる。時々刻々の非定常空気力を計算することができれば、計算された非定常空気力を振動方程式に代入することによって時々刻々の振動方程式を定めて、定められた振動方程式を解くことによって、空力弾性振動の時刻歴応答解析シミュレーションをすることができる。
本願の発明は、以上の実施例に限定されるものではなく、その要旨を逸脱しない範囲において、種々の変形が可能である。
上記実施例では、箱桁断面の橋梁模型(空力弾性振動体)を例にシミュレーションについて説明したが、本願の発明にかかるシミュレーション方法及び装置は、例えば、橋梁の桁部、主塔部など、弾性振動する種々の部位についての時刻歴応答シミュレーションに用いることができる。
また、上記実施例では、箱桁断面の橋梁模型を例に、バネ吊り1自由度振動系の空力弾性応答のシミュレーションについて説明したが、本願の発明にかかるシミュレーション方法及び装置は、その他にも、バネ吊り2自由度振動系の空力弾性応答や全橋の空力弾性振動応答、乱流中での橋桁の空力弾性振動応答のシミュレーションなど、種々の振動に関する空力弾性応答のシミュレーションに用いることができる。
また、上記実施例の自由振動実験では、風速を徐々に変化させ、0.2m/sごとに発散振動状況、減衰振動状況の振動波形を記録し、記録した振動波形から各風速における振幅ごとの減衰率を算出したが、波形を測定する風速の間隔は必ずしも0.2m/sごとでなければならないわけではなく、これにより大きくても良いし、小さくても良い。
Claims (8)
- 外力として非定常空気力を受ける弾性振動体の時刻歴応答シミュレーション方法であって、
非定常空気力を受ける前記弾性振動体の所定部位の時刻tにおける変位yt及び変位速度y’tの値を用いて前記所定部位の振幅y0の時刻tにおける値を計算する第1ステップと、
振幅y0の時刻tにおける値を用いて、振幅y0の関数である空力減衰率δa及び応答円振動数nの時刻tにおける値を計算する第2ステップと、
空力減衰率δaの時刻tにおける値を用いて前記非定常空気力の振動変位速度成分faIの時刻tにおける値を計算すると共に、応答円振動数nの時刻tにおける値を用いて前記非定常空気力の振動変位成分faRの時刻tにおける値を計算する第3ステップと、
振動変位成分faR及び振動変位速度成分faIを変数として含む弾性振動体の振動方程式に振動変位速度成分faI及び振動変位成分faRの時刻tにおける値を代入して得られた振動方程式を解いて、次の時刻t+Δtにおける前記所定部位の変位yt+Δt、変位速度y’t+Δt及び加速度y’’t+Δtの値を計算する第4ステップと、を備え、
これらのステップを含む時刻歴応答計算ルーチンを各時刻について繰り返し実行して計算された値を基に弾性振動体の時刻歴応答シミュレーションを行うものであり、
2回目以降の時刻歴応答計算ルーチンの第1ステップでは、当該第1ステップにおける変位yt及び変位速度y’tの値として、直前の時刻歴応答計算ルーチンの第4ステップで計算された変位y及び変位速度y’の値を用いて振幅y0の対応時刻における値を計算するものである時刻歴応答シミュレーション方法。 - 2回目以降の時刻歴応答計算ルーチンの第4ステップは、当該第4ステップにおける変位yt、変位速度y’t及び加速度y’’tの値として、直前の時刻歴応答計算ルーチンの第4ステップで計算された変位y、変位速度y’及び加速度y’’の値を用いて、次の時刻t+Δtにおける前記所定部位の変位yt+Δt、変位速度y’t+Δt及び加速度y’’t+Δtの値を計算するものである請求項1に記載の時刻歴応答シミュレーション方法。
- 前記所定部位の振幅y0の時刻tにおける値を求めるための前記変位y及び変位速度y’の関数は、時刻tにおける変位ytの自乗値と、時刻tにおける変位速度y’tを応答円振動数nで除した値の自乗値との和の平方根を振幅y0の時刻tにおける値として求めることができるものである請求項1又は請求項2に記載の時刻歴応答シミュレーション方法。
- 前記振動方程式は、下記の式(6)であり、前記第4ステップで振動変位速度成分faI及び振動変位成分faRの値が代入された振動方程式の解法は直接時間積分法であり、
前記第1ステップで前記振幅y0の値を計算するために用いられる変位y及び変位速度y’の関数は、下記の式(10)であり、
前記第3ステップで前記振動変位速度成分faIの値を計算するために用いられる前記空力減衰率δaの関数は、下記の式(16)であり、
前記第3ステップで前記振動変位成分faRの値を計算するために用いられる前記応答円振動数nの関数は、下記の式(18)である請求項1から請求項3のいずれか一項に記載の時刻歴応答シミュレーション方法。
- 前記空力減衰率δaを求めるための前記振幅y0の関数(δa=fI)及び前記応答円振動数nを求めるための前記振幅y0の関数(n=fR)は、前記弾性振動体について行う振動実験によって求められるものである請求項1から請求項4のいずれか一項に記載の時刻歴応答シミュレーション方法。
- 前記関数(δa=fI)及び前記関数(n=fR)は、前記振幅y0及び風速の関数である請求項5に記載の時刻歴応答シミュレーション方法。
- 外力として非定常空気力を受ける弾性振動体の時刻歴応答シミュレーションをするためのコンピュータを、
シミュレーションの対象である前記弾性振動体の所定部位の条件を設定する手段と、
風速条件を設定する手段、
前記弾性振動体の所定部位の初期時刻における変位y及び変位速度y’の値を設定する手段、
前記弾性振動体の所定部位の時刻tにおける変位yt及び変位速度y’tの値を用いて当該所定部位の振幅y0の時刻tにおける値を計算する手段、
前記風速条件に基づいて特定された時刻tにおける風速と、振幅y0の時刻tにおける値に基づいて空力減衰率δa及び応答円振動数nの時刻tにおける値を計算する手段、
空力減衰率δaの時刻tにおける値を用いて前記非定常空気力の振動変位速度成分faIの時刻tにおける値を計算すると共に、応答円振動数nの時刻tにおける値を用いて前記非定常空気力の振動変位成分faRの時刻tにおける値を計算する手段、
振動変位成分faR及び振動変位速度成分faIを変数として含む弾性振動体の振動方程式に、振動変位速度成分faI及び振動変位成分faRの時刻tにおける値を代入して得られた振動方程式を解いて、次の時刻t+Δtにおける前記所定部位の変位yt+Δt、変位速度y’t+Δt及び加速度y’’t+Δtの値を計算する手段、
前記計算する手段を順次実行する時刻歴応答計算ルーチンを各時刻について繰り返し実行させる手段であって、2回目以降の時刻歴応答計算ルーチンの振幅y0の値の計算で用いられる変位yt及び変位速度y’tの値として、直前の時刻歴応答計算ルーチンで計算された変位y及び変位速度y’の値を用いて振幅y0の値を計算させる手段、
前記時刻歴応答計算ルーチンを各時刻について繰り返し実行して計算された値を用いて弾性振動体の時刻歴応答シミュレーションを出力する手段
として機能させる時刻歴応答シミュレーションプログラム。 - 外力として非定常空気力を受ける弾性振動体の時刻歴応答シミュレーションを行う時刻歴応答シミュレーション装置であって、
非定常空気力を受ける前記弾性振動体の所定部位の時刻tにおける変位yt及び変位速度y’tの値を用いて前記所定部位の振幅y0の時刻tにおける値を計算する第1計算手段と、
振幅y0の時刻tにおける値を用いて、振幅y0の関数である空力減衰率δa及び応答円振動数nの時刻tにおける値を計算する第2計算手段と、
空力減衰率δaの時刻tにおける値を用いて前記非定常空気力の振動変位速度成分faIの時刻tにおける値を計算すると共に、応答円振動数nの時刻tにおける値を用いて前記非定常空気力の振動変位成分faRの時刻tにおける値を計算する第3計算手段と、
振動変位成分faR及び振動変位速度成分faIを変数として含む弾性振動体の振動方程式に、振動変位速度成分faI及び振動変位成分faRの時刻tにおける値を代入して得られた振動方程式を解いて、次の時刻t+Δtにおける前記所定部位の変位yt+Δt、変位速度y’t+Δt及び加速度y’’t+Δtの値を計算する第4計算手段と、
前記計算する手段を順次実行する時刻歴応答計算ルーチンを各時刻について繰り返し実行させる手段であって、2回目以降の時刻歴応答計算ルーチンの振幅y0の値の計算で用いられる変位yt及び変位速度y’tの値として、直前の時刻歴応答計算ルーチンで計算された変位y及び変位速度y’値を用いて振幅y0の値を計算させる手段と、
前記時刻歴応答計算ルーチンを各時刻について繰り返し実行して計算された値を用いて弾性振動体の時刻歴応答シミュレーションを出力する手段と、を備える時刻歴応答シミュレーション装置。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2008283798A JP2010112769A (ja) | 2008-11-05 | 2008-11-05 | 時刻歴応答シミュレーション方法、そのプログラム及び装置 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2008283798A JP2010112769A (ja) | 2008-11-05 | 2008-11-05 | 時刻歴応答シミュレーション方法、そのプログラム及び装置 |
Publications (1)
Publication Number | Publication Date |
---|---|
JP2010112769A true JP2010112769A (ja) | 2010-05-20 |
Family
ID=42301440
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
JP2008283798A Pending JP2010112769A (ja) | 2008-11-05 | 2008-11-05 | 時刻歴応答シミュレーション方法、そのプログラム及び装置 |
Country Status (1)
Country | Link |
---|---|
JP (1) | JP2010112769A (ja) |
Cited By (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103335858A (zh) * | 2013-06-06 | 2013-10-02 | 湖南大学 | 一种桥梁结构动态位移和振动频率的测量方法 |
CN105758602A (zh) * | 2016-03-01 | 2016-07-13 | 西南交通大学 | 一种桁架梁桥断面抖振力同步测量方法 |
CN105910779A (zh) * | 2016-05-09 | 2016-08-31 | 中国地震局工程力学研究所 | 一种桥梁梁体形变传感***及传感器 |
CN109579969A (zh) * | 2018-11-29 | 2019-04-05 | 上海交通大学 | 叶轮在加减速瞬态工况下最大振动幅值的获取方法及*** |
CN114088322A (zh) * | 2021-11-18 | 2022-02-25 | 天元建设集团有限公司 | 一种桥梁连接组件地震力响应数值模拟用实验平台 |
CN114112292A (zh) * | 2021-11-22 | 2022-03-01 | 中国建筑第五工程局有限公司 | 一种大跨度桥梁高阶涡振风洞试验模型***及其试验方法 |
CN114398712A (zh) * | 2022-03-25 | 2022-04-26 | 西南交通大学 | 斜拉索实时索力的计算方法、装置、设备及可读存储介质 |
-
2008
- 2008-11-05 JP JP2008283798A patent/JP2010112769A/ja active Pending
Cited By (10)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103335858A (zh) * | 2013-06-06 | 2013-10-02 | 湖南大学 | 一种桥梁结构动态位移和振动频率的测量方法 |
CN105758602A (zh) * | 2016-03-01 | 2016-07-13 | 西南交通大学 | 一种桁架梁桥断面抖振力同步测量方法 |
CN105910779A (zh) * | 2016-05-09 | 2016-08-31 | 中国地震局工程力学研究所 | 一种桥梁梁体形变传感***及传感器 |
CN109579969A (zh) * | 2018-11-29 | 2019-04-05 | 上海交通大学 | 叶轮在加减速瞬态工况下最大振动幅值的获取方法及*** |
CN109579969B (zh) * | 2018-11-29 | 2020-05-22 | 上海交通大学 | 叶轮在加减速瞬态工况下最大振动幅值的获取方法及*** |
CN114088322A (zh) * | 2021-11-18 | 2022-02-25 | 天元建设集团有限公司 | 一种桥梁连接组件地震力响应数值模拟用实验平台 |
CN114088322B (zh) * | 2021-11-18 | 2023-04-28 | 天元建设集团有限公司 | 一种桥梁连接组件地震力响应数值模拟用实验平台 |
CN114112292A (zh) * | 2021-11-22 | 2022-03-01 | 中国建筑第五工程局有限公司 | 一种大跨度桥梁高阶涡振风洞试验模型***及其试验方法 |
CN114112292B (zh) * | 2021-11-22 | 2023-08-08 | 中国建筑第五工程局有限公司 | 一种大跨度桥梁高阶涡振风洞试验模型***及其试验方法 |
CN114398712A (zh) * | 2022-03-25 | 2022-04-26 | 西南交通大学 | 斜拉索实时索力的计算方法、装置、设备及可读存储介质 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
JP2010112769A (ja) | 時刻歴応答シミュレーション方法、そのプログラム及び装置 | |
Siedziako et al. | An enhanced forced vibration rig for wind tunnel testing of bridge deck section models in arbitrary motion | |
Šarkić et al. | Bridge flutter derivatives based on computed, validated pressure fields | |
Zhu et al. | Vortex-induced vibration analysis of long-span bridges with twin-box decks under non-uniformly distributed turbulent winds | |
Ma et al. | Aerodynamic characteristics of a long-span cable-stayed bridge under construction | |
Shrivastava et al. | Identification of unbalance in a rotor system using a joint input-state estimation technique | |
Sampaio et al. | A simple damage detection indicator using operational deflection shapes | |
Mukhopadhyay | Structural Dynamics | |
Liu et al. | Influence of alternating loads on nonlinear vibration characteristics of cracked blade in rotor system | |
Brewick et al. | An investigation of the effects of traffic induced local dynamics on global damping estimates using operational modal analysis | |
Valverde-Marcos et al. | Influence of the rotation speed on the dynamic behaviour of a cracked rotating beam | |
Barni et al. | Time-variant self-excited force model based on 2D rational function approximation | |
Chowdhury et al. | Experimental identification of rational function coefficients for time-domain flutter analysis | |
Sinou | On the use of non-linear vibrations and the anti-resonances of Higher-Order Frequency Response Functions for crack detection in pipeline beam | |
Salvatori et al. | Effects of structural nonlinearity and along-span wind coherence on suspension bridge aerodynamics: Some numerical simulation results | |
Caracoglia et al. | Comparative and sensitivity study of flutter derivatives of selected bridge deck sections, Part 2: Implications on the aerodynamic stability of long-span bridges | |
Santos et al. | Strain-based dynamic measurements and modal testing | |
Sauder et al. | Real-time prediction of aeroelastic loads of wind turbine blades in gusty and turbulent wind using an improved load model | |
CN115046727A (zh) | 一种基于希尔伯特变换的桥梁非线性涡激力精确识别方法 | |
Kim et al. | Development of artificial neural network model for estimation of cable tension of cable-stayed bridge | |
Al-Anbaki et al. | Duality between time and frequency domains for vibration serviceability analysis of floor structures | |
Markowski et al. | Modelling and vibration analysis of some complex mechanical systems | |
Chiulan et al. | Experimental study on an aeroelastic model of a television tower in the boundary layer wind tunnel | |
Siedziako et al. | A new experimental approach to validation of aerodynamic derivatives based model for self-excited forces | |
Caracoglia et al. | Nonlinear computer model for the simulation of lock‐in vibration on long‐span bridges |