JP6664690B2

JP6664690B2 - シミュレーション装置及びプログラム

Info

Publication number: JP6664690B2
Application number: JP2016052739A
Authority: JP
Inventors: 幸人中野; 雄一松村
Original assignee: Gifu University NUC; Toyota Central R&D Labs Inc
Current assignee: Gifu University NUC; Toyota Central R&D Labs Inc
Priority date: 2016-03-16
Filing date: 2016-03-16
Publication date: 2020-03-13
Anticipated expiration: 2036-03-16
Also published as: JP2017166992A

Description

本発明は、シミュレーション装置及びプログラムに係り、特に、３次元はり構造への入力振幅に対する応答としての波動振幅を計算するシミュレーション装置及びプログラムに関する。

従来より、数値解析によって周波数応答関数を求める手法として、質量行列、剛性行列を利用した２つの手法が知られている（非特許文献１）。１つ目の手法は、直接法と呼ばれ、質量行列と剛性行列とを組み合わせた行列の逆行列を直接求めることにより、周波数応答関数が求まる。もうひとつの手法は、理論モード解析と呼ばれ、質量行列と剛性行列とを組み合わせた行列を固有値解析にかけることにより、モード特性を求め、各モードの応答の重ね合わせで周波数応答関数を求める手法である。この手法では周波数応答関数における各モードの寄与を評価することが可能である。

また、分布定数系モデルで複雑な構造物の波動伝播解析を実施している事例として、宇宙構造物を対象にした時系列の波動伝播数値解析の事例が知られている（非特許文献２）。

長松昭男、「モード解析入門」、コロナ社、１９９３年．西田明美、「３次元フレーム構造物の波動伝播特性に関する研究：ティモシェンコ梁理論の導入」、構造工学論文集Ｂ、Ｖｏｌ．５２Ｂ、pp.119-124、日本建築学会、２００６年

上記非特許文献１の周波数応答関数を求める解析技術では、振動を構成している構造物内部の各波動の伝播を物理的に考慮しておらず、波動伝播を考慮しない別の物理モデル（離散系モデル）から周波数応答関数を求めている。

また、上記非特許文献２に記載の、波動伝播に関する研究・技術開発では、衝撃荷重負荷時の時系列での波動伝播解析が主であり、周波数応答関数などの定常状態での波動伝播の扱いに関する研究事例はない。

このように、従来法のモード解析や周波数応答関数から波動伝播寄与を求めることができない。構造物の波動伝播寄与を求めるには別の手法がいる。

本発明は、上記の事情を鑑みてなされたもので、波動振幅の成分別に周波数応答関数を求めることができるシミュレーション装置及びプログラムを提供することを目的とする。

上記の目的を達成するために本発明に係るシミュレーション装置は、複数のはり要素が結合した３次元はり構造への入力振幅に対する応答としての波動振幅を計算するシミュレーション装置であって、入力振幅を表す、はり要素の軸方向であるｘ方向の縦波前進伝播波、ｘ方向の縦波後退伝播波、ねじり前進伝播波、ねじり後退伝播波、ｙ方向の曲げ前進伝播波、ｙ方向の曲げ前進近接波、ｙ方向の曲げ後退伝播波、ｙ方向の曲げ後退近接波、ｚ方向の曲げ前進伝播波、ｚ方向の曲げ前進近接波、ｚ方向の曲げ後退伝播波、及びｚ方向の曲げ後退近接波の各成分からなる波動振幅ベクトルａ₀と、前記３次元はり構造の不連続部の透過反射性を表す散乱行列Ｓと、各成分の伝播に伴った位相の変化を表す伝播行列Ｄとに基づいて、以下の式に従って、応答としての波動振幅を表す、各成分の波動振幅ベクトルａを計算する計算手段を含んで構成されている。

ただし、Ｉは、単位行列である。

本発明に係るプログラムは、複数のはり要素が結合した３次元はり構造への入力振幅に対する応答としての波動振幅を計算するためのプログラムであって、コンピュータを、入力振幅を表す、はり要素の軸方向であるｘ方向の縦波前進伝播波、ｘ方向の縦波後退伝播波、ねじり前進伝播波、ねじり後退伝播波、ｙ方向の曲げ前進伝播波、ｙ方向の曲げ前進近接波、ｙ方向の曲げ後退伝播波、ｙ方向の曲げ後退近接波、ｚ方向の曲げ前進伝播波、ｚ方向の曲げ前進近接波、ｚ方向の曲げ後退伝播波、及びｚ方向の曲げ後退近接波の各成分からなる波動振幅ベクトルａ₀と、前記３次元はり構造の不連続部の透過反射性を表す散乱行列Ｓと、各成分の伝播に伴った位相の変化を表す伝播行列Ｄとに基づいて、上記の式に従って、応答としての波動振幅を表す、各成分の波動振幅ベクトルａを計算する計算手段として機能させるためのプログラムである。

本発明によれば、計算手段が、入力振幅を表す、はり要素の軸方向であるｘ方向の縦波前進伝播波、ｘ方向の縦波後退伝播波、ねじり前進伝播波、ねじり後退伝播波、ｙ方向の曲げ前進伝播波、ｙ方向の曲げ前進近接波、ｙ方向の曲げ後退伝播波、ｙ方向の曲げ後退近接波、ｚ方向の曲げ前進伝播波、ｚ方向の曲げ前進近接波、ｚ方向の曲げ後退伝播波、及びｚ方向の曲げ後退近接波の各成分からなる波動振幅ベクトルａ₀と、前記３次元はり構造の不連続部の透過反射性を表す散乱行列Ｓと、各成分の伝播に伴った位相の変化を表す伝播行列Ｄとに基づいて、上記の式に従って、応答としての波動振幅を表す、各成分の波動振幅ベクトルａを計算する。

このように、３次元はり構造への入力振幅に対する応答として、ｘ方向の縦波前進伝播波、ｘ方向の縦波後退伝播波、ねじり前進伝播波、ねじり後退伝播波、ｙ方向の曲げ前進伝播波、ｙ方向の曲げ前進近接波、ｙ方向の曲げ後退伝播波、ｙ方向の曲げ後退近接波、ｚ方向の曲げ前進伝播波、ｚ方向の曲げ前進近接波、ｚ方向の曲げ後退伝播波、及びｚ方向の曲げ後退近接波の各成分からなる波動振幅ベクトルを計算することにより、波動振幅の成分別に周波数応答関数を求めることができる。

上記の計算手段は、周波数毎に、各成分の波動振幅ベクトルａを計算し、各成分について、周波数応答関数を求めることができる。

上記の計算手段は、前記周波数毎に、前記複数のはり要素の各々について、前記はり要素の縦弾性係数、体積質量密度、横弾性係数、断面二次極モーメント、ねじり定数、及び断面積に基づいて、ｘ方向の縦波の波数ｋ_L、ねじり波の波数ｋ_T、ｙ方向の曲げ波の波数ｋ_By、及びｚ方向の曲げ波の波数ｋ_Bzを計算する波数計算手段と、前記周波数毎に、前記複数のはり要素の各々について、前記ｘ方向の縦波の波数ｋ_L、前記ねじり波の波数ｋ_T、前記ｙ方向の曲げ波の波数ｋ_By、及び前記ｚ方向の曲げ波の波数ｋ_Bzに基づいて、前進伝播波成分及び後退伝播波成分の各々についての変位波動振幅変換行列Ψ₁、Ψ₂と、前進波及び後退波の各々についての伝播距離ｘによる位相変化を表す伝播行列Ｇ₁（ｘ）、Ｇ₂（ｘ）と、前進伝播波成分及び後退伝播波成分の各々についての内力波動振幅変換行列Φ₁、Φ₂とを計算する変換行列伝播行列計算手段と、前記はり要素の各々について、前記はり要素の全体座標系の基底ベクトルｅi及び要素座標系の基底ベクトルＥiに基づいて、不連続部において結合している前記はり要素の各々についての座標変換行列Ｃ⁽ⁱ⁾を計算する座標変換行列計算手段と、前記はり要素の各々について、前記はり要素の前進伝播波の、前記不連続部への入射方向ｄ⁽ⁱ⁾を計算する入射方向計算手段と、前記周波数毎に、前記はり要素の各々についての、前記変位波動振幅変換行列Ψ₁、Ψ₂と、前記伝播行列Ｇ₁（ｘ）、Ｇ₂（ｘ）と、前記内力波動振幅変換行列Φ₁、Φ₂と、前記座標変換行列Ｃ⁽ⁱ⁾と、前記入射方向ｄ⁽ⁱ⁾とに基づいて、前記はり要素の各々についての反射係数ｒと、不連続部において結合している前記はり要素のペアの各々について透過係数ｔとを計算する係数計算手段と、前記周波数毎に、前記はり要素の各々についての反射係数ｒと、前記はり要素のペアの各々について透過係数ｔとに基づいて、前記３次元はり構造の不連続部の透過反射性を表す散乱行列Ｓを計算する散乱行列計算手段と、前記周波数毎に、前記はり要素の各々についての前記伝播行列Ｇ₁（ｘ）、Ｇ₂（ｘ）と、前記はり要素の各々の導波路長さＬ⁽ⁱ⁾とに基づいて、前記はり要素の各々の位相情報を表す伝播行列Ｄを計算する伝播行列計算手段と、前記周波数毎に、前記波動振幅ベクトルａ₀と、前記散乱行列Ｓと、前記伝播行列Ｄとに基づいて、各成分の波動振幅ベクトルａを計算し、各成分について、周波数応答関数を求める周波数応答関数計算手段とを含むことができる。

以上説明したように、本発明のシミュレーション装置及びプログラムによれば、３次元はり構造への入力振幅に対する応答として、ｘ方向の縦波前進伝播波、ｘ方向の縦波後退伝播波、ねじり前進伝播波、ねじり後退伝播波、ｙ方向の曲げ前進伝播波、ｙ方向の曲げ前進近接波、ｙ方向の曲げ後退伝播波、ｙ方向の曲げ後退近接波、ｚ方向の曲げ前進伝播波、ｚ方向の曲げ前進近接波、ｚ方向の曲げ後退伝播波、及びｚ方向の曲げ後退近接波の各成分からなる波動振幅ベクトルを計算することにより、波動振幅の成分別に周波数応答関数を求めることができる、という効果が得られる。

本発明の実施の形態に係るシミュレーション装置を示すブロック図である。本発明の実施の形態に係るシミュレーション装置のシミュレーション処理ルーチンの内容を示すフローチャートである。３次元はり構造の一例を示す図である。周波数応答関数の一例を示す図である。周波数応答関数の一例を示す図である。

以下、図面を参照して本発明の実施の形態を詳細に説明する。

＜本発明の実施の形態の概要＞
透過・反射係数と形状情報（導波路長さ等）と波動振幅とでレイトレースモデルと呼ばれるモデルを構成できる。このモデルは波動が構造物の不連続部から別の不連続部まで伝播する際の伝播前後の波動振幅の関係を数値的に記述したモデルである。入力が１Ｎの場合の波動振幅を計算し、波動伝播が無限回起きた時点のある任意位置での各種の波動振幅を個別にみることにより、周波数応答関数における各種波動振幅寄与を評価する。

＜本発明の実施の形態の原理＞
次に、波動伝播による波動振幅解析の原理について説明する。

＜3次元はり要素における波動振幅表現＞
まず、3次元のはり要素における波動理論を整理する。縦波、ねじり波、曲げ波についての波動振幅と変位および内力の関係を整理し、透過・反射係数の算出に必要な波動振幅から変位および内力への変換行列を導出する。

＜縦波についての波動理論＞
はりを伝播する縦波についての運動方程式は次式となる。ここでρは体積質量密度[kg/m³]、Aは断面積[m²]、Eは縦弾性係数[N/m²]を示す。

上式の運動方程式の一般解は、次式の通りである。

ここで、a_1Lはxの正の方向へ伝播する前進伝播波の波動振幅、a_2Lはxの負の方向へ伝播する後退伝播波の波動振幅を表す。また、内力F(x)は、次式より求まる。

k_Lは、縦波の波数[rad/m]であり、次式で定義される。

＜ねじり波についての波動理論＞
はり要素の軸方向をx軸とし、断面のx軸周りの変位をs(x,t)とする。また、ねじり波伝播のねじりモーメントをh(x,t)とする。ここでGは横弾性係数[N/m²]、IPは断面二次極モーメント[m⁴]、Jはねじり定数を示す。線形振動領域では、断面変形が小さく、ねじり波についての運動方程式は次式で表される。

上式の運動方程式の一般解は、次式の通りである。

また、内力H(x)は、次式より求まる。

ここで、a_1Tはxの正の方向へ伝播するねじり前進伝播波の振幅、a_2Tはxの負の方向へ伝播するねじり後退伝播波の振幅を表す。このとき、k_Tはねじり波の波数[rad/m]である。

＜曲げ波についての波動理論＞
曲げ波に関する運動方程式は次式で表すことができる。

ここで、Iは断面二次極モーメント[m⁴]である。本式の運動方程式の解は、次式で得られる。

ここで、a_1Pは曲げ前進伝播波の振幅、a_2Pは曲げ後退伝播波の振幅、a_1Nは曲げ前進近接波の振幅、a_2Nは曲げ後退近接波の振幅を表す。ここで、曲げ波の波数k_Bは次式の通りである。

このとき、たわみ角θ(x,t)とたわみv(x,t)との関係は次式である。

また、モーメントm(x,t)とせん断力w(x,t)は、それぞれ次式で求まる。

＜波動振幅変換行列＞
はりの反射・透過係数を求める際に、はりの変位、内力を波動振幅で表すための波動振幅変換行列が必要である。次に、波動振幅変換行列の定義を述べる。

まず、波動振幅ベクトルの前進伝播波成分a₁および後退伝播波成分a₂を次式の通り定義する。

縦波による変位の一般解である式(2)、ねじり波による変位の一般解である式(6)、曲げ波についての変位の一般解である式(10)、式(12)を、前進波a₁による変位成分および後退波a₂による変位成分に分ける。縦波の変位をU₁、U₂、 y軸方向の曲げ波のたわみをV_1y、V_2y、z軸方向の曲げ波のたわみをV_1z、V_2z、ねじり波の変位をS₁、S₂、y軸方向のたわみを角Θ_1y、Θ_2y、z軸方向のたわみ角をΘ_1z、Θ_2zとする。ここで、変位の一般解についてのベクトルをU、前進波および後退波による変位ベクトルをU₁、U₂として次式で定義する。

これらの定義より、式(2)、式(6)、式(10)、式(12)を用いて、波動振幅の関係をまとめると、次式で表すことができる。

このとき、Ψ₁、 Ψ₂はそれぞれ波動振幅の前進伝播波成分および後退伝播波成分についての変位-波動振幅変換行列であり、それぞれ次の２つの式で表すことができる。

また、G₁(x)、 G₂(x)はそれぞれ、前進波および後退波の伝播距離xによる位相変化をあらわす伝播行列であり、式(20)および式(21)で表される。

次に、縦波による内力の解である式(3)、ねじり波による内力の解である式(7)、曲げ波による内力の解である式(13)、式(14)を、前進波および後退波による成分に分ける。縦波の内力F₁、F₂、y軸方向の曲げ波のせん断力W_1y、W_2y、z軸方向の曲げ波のせん断力W_1z、W_2z、ねじり波の内力H₁、H₂、y軸方向の曲げモーメントM_1y、M_2y、z軸方向の曲げモーメントM_1z、M_2zとする。ここで、内力の一般解についてのベクトルをF、前進波および後退波による内力ベクトルをF₁、F₂として次式で定義する。

これらの定義より、式(3)、式(7)、式(13)、式(14)をまとめると、内力と波動振幅の関係をまとめると、次式の通り表すことができる。

このときの行列Φ₁、Φ₂はそれぞれ波動振幅の前進波成分および後退波成分についての内力−波動振幅変換行列であり、それぞれ次の２つの式で表される。ここで、I_yはy軸方向の断面二次極モーメント[m⁴]、I_zはz軸方向の断面二次極モーメント[m⁴]を示す。

以上より、縦波、ねじり波、曲げ波の一般解を整理し、波動振幅から変位および内力への変換行列を導出した。

＜不連続部の透過・反射係数の算出＞
次に、N_c個の要素が結合する不連続部における反射・透過係数の算出について定式化する。a₁ ⁽ⁱ⁾、 a₂ ⁽ⁱ⁾は、それぞれ結合要素i ( = 1〜N_c)の要素座標系のx軸の正方向に伝播する前進波と、負の向きに伝播する後退波についての波動振幅ベクトルである。

不連続部において結合するはり要素i ( = 1〜N_c)の要素座標系から全体座標系への座標変換行列C⁽ⁱ⁾は、全体座標系の基底ベクトルをe、要素iの要素座標系の基底ベクトルをE⁽ⁱ⁾として、次式で表される。

式(17)において不連続部位置をx=0とすると、不連続部で結合するはり要素iの不連続部での変位U⁽ⁱ⁾は、次式で表される。

不連続部における変位の連続性より、次式が成り立つ。

また、式(23)において不連続部の位置をx = 0とすると、不連続部で結合するはり要素iの不連続部での内力F⁽ⁱ⁾は、次式で表される。

不連続部における力の釣り合いより、次式が成立する。

本式におけるd⁽ⁱ⁾は要素iの前進伝播波の不連続部への入射方向の正負を表しており、次式であらわされる。

ここで、w_I ⁽ⁱ⁾およびw_O ⁽ⁱ⁾は、不連続部に対する入射波および反射波を定義する変数であり、式(32)、式(33)で表される。なお、x_cは、不連続部の座標、x₁ ⁽ⁱ⁾およびx₂ ⁽ⁱ⁾は、要素iの構成節点座標であり、x₁ ⁽ⁱ⁾からx₂ ⁽ⁱ⁾の方向が要素座標系のx軸の正方向である。w_I ⁽ⁱ⁾は、不連続部に対する入射する波動振幅ベクトルが要素座標系のx軸の正方向に伝播する前進波成分であるか負の方向に伝播する後退波成分であるかを定める変数であり、1のとき要素iの前進波が不連続部に入射し、2のとき後退波が入射するように要素座標系が定義されていることになる。一方、w_O ⁽ⁱ⁾は、不連続部から反射する波動振幅ベクトルを定める変数で、1のとき前進波が不連続部から反射し、2のとき後退波が不連続部から反射するように、要素座標系が定義されていることになる。

不連続部における変位の連続を表す式(28)および、力のつり合いを表す式(30)を、左辺に不連続部に波動が入射する波動振幅成分、右辺に不連続部から波動振幅が反射される波動振幅成分に整理し、行列表記すると次式で表される。

また、透過・反射係数を用いて、不連続部の波動振幅の変化を表すと次式の通りである。

したがって、式(34)および式(35)より、(34)反射係数rおよび透過係数tは、変位−波動振幅変換行列Ψ、内力−波動振幅変換行列Φ、座標変換行列Cを用いて、次式で算出することが可能である。

＜レイトレースモデルの一般化＞
次にN個の要素から構成される解析対象のレイトレースモデルを定式化する。はり要素i (=1〜N)を伝播する前進波の波動振幅a₁ ⁽ⁱ⁾および後退波の波動振幅成分a₂ ⁽ⁱ⁾からなる波動振幅ベクトルa⁽ⁱ⁾を次式の通り定義する。

不連続部にて結合されるはり要素iとはり要素jの間の透過・反射の関係について定式化すると、次式のように表すことができる。

ここで、R⁽ⁱ⁾は、不連続部におけるはり要素iから入射した波動に対する反射係数行列、T(i,j)ははり要素iからはり要素jに波動が透過する際の透過係数行列である。透過係数行列Tおよび反射係数行列Rは、前節の式(36)で算出された透過反射係数と、式(32)で示したw_Iおよびw_Oを用いて表すと、不連続部に対する波動ベクトルの入射および反射方向の定義を考慮し、それぞれ、式(39)および式(40)で表される。

以上より、透過・反射の関係が定式化される。系を構成するはり要素iについての反射係数R⁽ⁱ⁾および、不連続部にて結合されているはり要素iとはり要素jの間の透過係数行列T(i,j)を使用して、次式のような系全体の散乱行列Sを構築することが可能である。

ここで、系全体の要素についての波動振幅ベクトルを並べたaを、次式の通り定義する。

また、波動が不連続部から発生し、はり要素内をもう一方の不連続部まで伝播することを表す伝播行列Dを次式の通り定義する。

ある波動振幅ベクトルa_nがはり要素内を伝播し、不連続部に到達し、次に伝播する波動振幅ベクトルa_n+1に変化する波動モデルをレイトレースモデルと呼び、伝播行列Dと散乱行列Sを使用して次式で表すことができる。

レイトレースモデルを用いると、初期波動伝播前後の波動振幅ａ₁の関係は以下のようになる

定常状態の波動振幅ａを考えると、複数回伝播した波動振幅(ａ₀，ａ₁，ａ₂…ａ_∞)が混在する状況であり以下のように定式化できる。

これをまとめると以下の式になる。

ａ₀が構造物の任意箇所に１Nの振幅の加振力を付加された時の初期波動振幅だとすると、上記のａは周波数応答関数を各種波動振幅別に分解した波動振幅ベクトルであり、周波数応答関数の波動振幅寄与を表現している。よって、上式を利用すれば周波数応答関数の波動振幅寄与を求めることが可能である。

＜システム構成＞
図１に示すように、本発明の実施の形態に係るシミュレーション装置１０は、ＣＰＵ１２、ＲＯＭ１４、ＲＡＭ１６、ＨＤＤ１８、通信インタフェース２０、及びこれらを相互に接続するためのバス２２を備えている。

ＣＰＵ１２は、各種プログラムを実行する。ＲＯＭ１４には、各種プログラムやパラメータ等が記憶されている。ＲＡＭ１６は、ＣＰＵ１２による各種プログラムの実行時におけるワークエリア等として用いられる。記録媒体としてのＨＤＤ１８には、後述する微分計算処理ルーチンを実行するためのプログラムを含む各種プログラムや各種データが記憶されている。

本実施の形態におけるシミュレーション装置１０は、入力として、入力振幅を表す、はり要素の軸方向であるｘ方向の縦波前進伝播波、ｘ方向の縦波後退伝播波、ねじり前進伝播波、ねじり後退伝播波、ｙ方向の曲げ前進伝播波、ｙ方向の曲げ前進近接波、ｙ方向の曲げ後退伝播波、ｙ方向の曲げ後退近接波、ｚ方向の曲げ前進伝播波、ｚ方向の曲げ前進近接波、ｚ方向の曲げ後退伝播波、及びｚ方向の曲げ後退近接波の各成分からなる波動振幅ベクトルａ₀を受け付ける。また、シミュレーション装置１０は、複数のはり要素が結合した３次元はり構造の形状情報として、各はり要素についての縦弾性係数Ｅ、体積質量密度ρ、横弾性係数Ｇ、断面二次極モーメントＩ_p、ねじり定数Ｊ、断面積Ａ、全体座標系の基底ベクトルｅ、要素座標系の基底ベクトルＥ⁽ⁱ⁾、及び導波路長さＬ⁽ⁱ⁾を受け付ける。

シミュレーション装置１０は、周波数毎に、波動振幅ベクトルａ₀と、３次元はり構造の不連続部の透過反射性を表す散乱行列Ｓと、はり要素の各々の位相情報を表す伝播行列Ｄとに基づいて、上記式（４７）に従って、応答としての波動振幅を表す、各成分からなる波動振幅ベクトルａを計算し、各成分について、周波数応答関数を求める。

＜シミュレーション装置の動作＞
次に、本実施の形態に係るシミュレーション装置１０の作用について説明する。

まず、シミュレーション装置１０に、波動振幅ベクトルａ₀と、上述した３次元はり構造の形状情報が入力されると、シミュレーション装置１０によって、３次元はり構造の形状モデルが生成される。そして、シミュレーション装置１０によって、図２に示すシミュレーション処理ルーチンが実行される。

まず、ステップＳ１００において、周波数ω毎に、各はり要素の縦弾性係数Ｅ、体積質量密度ρ、横弾性係数Ｇ、断面二次極モーメントＩp、ねじり定数Ｊ、及び断面積Ａに基づいて、上記式（４）、（８）、（１１）に従って、各はり要素について、ｘ方向の縦波の波数ｋ_L、ねじり波の波数ｋ_T、ｙ方向の曲げ波の波数ｋ_By、及びｚ方向の曲げ波の波数ｋ_Bzを計算する。

ステップＳ１０２では、周波数ω毎に、各はり要素について、ｘ方向の縦波の波数ｋ_L、ねじり波の波数ｋ_T、ｙ方向の曲げ波の波数ｋ_By、及びｚ方向の曲げ波の波数ｋ_Bzに基づいて、上記式（１８）〜（２１）、（２４）、（２５）に従って、前進伝播波成分及び後退伝播波成分の各々についての変位波動振幅変換行列Ψ₁、Ψ₂と、前進波及び後退波の各々についての伝播距離ｘによる位相変化を表す伝播行列Ｇ₁（ｘ）、Ｇ₂（ｘ）と、前進伝播波成分及び後退伝播波成分の各々についての内力波動振幅変換行列Φ₁、Φ₂とを計算する。

ステップＳ１０４では、はり要素の各々についての当該はり要素の全体座標系の基底ベクトルｅ及び要素座標系の基底ベクトルＥ⁽ⁱ⁾に基づいて、上記式（２６）に従って、不連続部において結合しているはり要素の各々についての座標変換行列Ｃ⁽ⁱ⁾を計算する。

また、はり要素の各々について、３次元はり構造の形状情報に基づいて、上記式（３１）に従って、当該はり要素の前進伝播波の、前記不連続部への入射方向ｄ⁽ⁱ⁾を計算する。

ステップＳ１０６では、周波数ω毎に、上記ステップＳ１０２で計算された、前進伝播波成分及び後退伝播波成分の各々についての変位波動振幅変換行列Ψ₁、Ψ₂と、前進波及び後退波の各々についての伝播距離ｘによる位相変化を表す伝播行列Ｇ₁（ｘ）、Ｇ₂（ｘ）と、前進伝播波成分及び後退伝播波成分の各々についての内力波動振幅変換行列Φ₁、Φ₂と、上記ステップＳ１０４で計算された、はり要素の各々についての前記座標変換行列Ｃ⁽ⁱ⁾と、前記はり要素の各々についての前記入射方向ｄ⁽ⁱ⁾とに基づいて、上記式（３６）に従って、はり要素の各々についての反射係数ｒと、不連続部において結合しているはり要素のペアの各々について透過係数ｔとを計算する。

そして、ステップＳ１０８では、周波数ω毎に、上記ステップＳ１０６で計算された、はり要素の各々についての反射係数ｒと、はり要素のペアの各々について透過係数ｔとに基づいて、上記式（４１）に従って、３次元はり構造の不連続部の透過反射性を表す散乱行列Ｓを計算する。

ステップＳ１１０では、周波数ω毎に、上記ステップＳ１０２で計算された前進波及び後退波の各々についての伝播距離ｘによる位相変化を表す伝播行列Ｇ₁（ｘ）、Ｇ₂（ｘ）と、はり要素の各々の導波路長さＬ⁽ⁱ⁾とに基づいて、上記式（４３）に従って、はり要素の各々の位相情報を表す伝播行列Ｄを計算する。

そして、ステップＳ１１２では、周波数ω毎に、波動振幅ベクトルａ₀と、上記ステップＳ１０８で計算された散乱行列Ｓと、上記ステップＳ１１０で計算された伝播行列Ｄとに基づいて、波動振幅ベクトルａを計算する。次のステップＳ１１４では、上記ステップＳ１１２で周波数ω毎に計算された波動振幅ベクトルａに基づいて、各成分について、周波数応答関数を求め、出力部（図示省略）により出力して、シミュレーション処理ルーチンを終了する。

＜実施例＞
断面特性が異なる３つのはり要素を結合した３次元はり構造のモデルで波動振幅寄与の例を示す。

図３に示すように、３次元はり構造モデルに対して、入力波動振幅として、拘束条件なしで左端に１Ｎの調和荷重を負荷している。評価点は加振点（１）とはり中央部（２）である。はり要素の諸元を以下の表１に記す。

解析結果を図４、図５に示す。３次元はり構造のすべての成分の波動の波動振幅寄与を表示すると複雑になるため、ここでは低周波の上下方向曲げ振動に着目し、前進伝播波、前進近接波、後退伝播波、後退近接波の周波数応答関数に対する寄与を表示する。図４は、評価点（１）における各成分の周波数応答関数を示しており、図５は、評価点（２）における各成分の周波数応答関数を示している。

図４、図５のグラフで示すように、周波数応答関数に対して各種成分の寄与を解析することが可能である。

以上説明したように、本発明の実施の形態に係るシミュレーション装置によれば、３次元はり構造への入力振幅に対する応答として、ｘ方向の縦波前進伝播波、ｘ方向の縦波後退伝播波、ねじり前進伝播波、ねじり後退伝播波、ｙ方向の曲げ前進伝播波、ｙ方向の曲げ前進近接波、ｙ方向の曲げ後退伝播波、ｙ方向の曲げ後退近接波、ｚ方向の曲げ前進伝播波、ｚ方向の曲げ前進近接波、ｚ方向の曲げ後退伝播波、及びｚ方向の曲げ後退近接波の各成分からなる波動振幅ベクトルを計算することにより、波動振幅の成分別に周波数応答関数を求めることができる。

また、周波数応答関数を構造物内の波動振幅成分の寄与として評価することで、振動応答の主原因となる波動を特定でき、効果的な対策を検討することが可能となる。

また、位相の異なる各種波動振幅の重ね合わせにより、任意位置の振動応答を低減する構造を提案することなども可能となる。

なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。

本発明のプログラムは、記憶媒体に格納して提供するようにしてもよい。

１０シミュレーション装置
１２ＣＰＵ
１４ＲＯＭ
１６ＲＡＭ
１８ＨＤＤ

Claims

複数のはり要素が結合した３次元はり構造への入力振幅に対する応答としての波動振幅を計算するシミュレーション装置であって、
入力振幅を表す、はり要素の軸方向であるｘ方向の縦波前進伝播波、ｘ方向の縦波後退伝播波、ねじり前進伝播波、ねじり後退伝播波、ｙ方向の曲げ前進伝播波、ｙ方向の曲げ前進近接波、ｙ方向の曲げ後退伝播波、ｙ方向の曲げ後退近接波、ｚ方向の曲げ前進伝播波、ｚ方向の曲げ前進近接波、ｚ方向の曲げ後退伝播波、及びｚ方向の曲げ後退近接波の各成分からなる波動振幅ベクトルａ₀と、前記３次元はり構造の不連続部の透過反射性を表す散乱行列Ｓと、各成分の伝播に伴った位相の変化を表す伝播行列Ｄとに基づいて、以下の式に従って、応答としての波動振幅を表す、各成分の波動振幅ベクトルａを計算する計算手段
を含むシミュレーション装置。
ただし、Ｉは、単位行列である。
前記計算手段は、周波数毎に、各成分の波動振幅ベクトルａを計算し、各成分について、周波数応答関数を求める請求項１記載のシミュレーション装置。
前記計算手段は、
前記周波数毎に、前記複数のはり要素の各々について、前記はり要素の縦弾性係数、体積質量密度、横弾性係数、断面二次極モーメント、ねじり定数、及び断面積に基づいて、ｘ方向の縦波の波数ｋ_L、ねじり波の波数ｋ_T、ｙ方向の曲げ波の波数ｋ_By、及びｚ方向の曲げ波の波数ｋ_Bzを計算する波数計算手段と、
前記周波数毎に、前記複数のはり要素の各々について、前記ｘ方向の縦波の波数ｋ_L、前記ねじり波の波数ｋ_T、前記ｙ方向の曲げ波の波数ｋ_By、及び前記ｚ方向の曲げ波の波数ｋ_Bzに基づいて、前進伝播波成分及び後退伝播波成分の各々についての変位波動振幅変換行列Ψ₁、Ψ₂と、前進波及び後退波の各々についての伝播距離ｘによる位相変化を表す伝播行列Ｇ₁（ｘ）、Ｇ₂（ｘ）と、前進伝播波成分及び後退伝播波成分の各々についての内力波動振幅変換行列Φ₁、Φ₂とを計算する変換行列伝播行列計算手段と、
前記はり要素の各々について、前記はり要素の全体座標系の基底ベクトルｅi及び要素座標系の基底ベクトルＥiに基づいて、不連続部において結合している前記はり要素の各々についての座標変換行列Ｃ⁽ⁱ⁾を計算する座標変換行列計算手段と、
前記はり要素の各々について、前記はり要素の前進伝播波の、前記不連続部への入射方向ｄ⁽ⁱ⁾を計算する入射方向計算手段と、
前記周波数毎に、前記はり要素の各々についての、前記変位波動振幅変換行列Ψ₁、Ψ₂と、前記伝播行列Ｇ₁（ｘ）、Ｇ₂（ｘ）と、前記内力波動振幅変換行列Φ₁、Φ₂と、前記座標変換行列Ｃ⁽ⁱ⁾と、前記入射方向ｄ⁽ⁱ⁾とに基づいて、前記はり要素の各々についての反射係数ｒと、不連続部において結合している前記はり要素のペアの各々について透過係数ｔとを計算する係数計算手段と、
前記周波数毎に、前記はり要素の各々についての反射係数ｒと、前記はり要素のペアの各々について透過係数ｔとに基づいて、前記３次元はり構造の不連続部の透過反射性を表す散乱行列Ｓを計算する散乱行列計算手段と、
前記周波数毎に、前記はり要素の各々についての前記伝播行列Ｇ₁（ｘ）、Ｇ₂（ｘ）と、前記はり要素の各々の導波路長さＬ⁽ⁱ⁾とに基づいて、前記はり要素の各々の位相情報を表す伝播行列Ｄを計算する伝播行列計算手段と、
前記周波数毎に、前記波動振幅ベクトルａ₀と、前記散乱行列Ｓと、前記伝播行列Ｄとに基づいて、各成分の波動振幅ベクトルａを計算し、各成分について、周波数応答関数を求める周波数応答関数計算手段とを含む請求項２記載のシミュレーション装置。
複数のはり要素が結合した３次元はり構造への入力振幅に対する応答としての波動振幅を計算するためのプログラムであって、
コンピュータを、
入力振幅を表す、はり要素の軸方向であるｘ方向の縦波前進伝播波、ｘ方向の縦波後退伝播波、ねじり前進伝播波、ねじり後退伝播波、ｙ方向の曲げ前進伝播波、ｙ方向の曲げ前進近接波、ｙ方向の曲げ後退伝播波、ｙ方向の曲げ後退近接波、ｚ方向の曲げ前進伝播波、ｚ方向の曲げ前進近接波、ｚ方向の曲げ後退伝播波、及びｚ方向の曲げ後退近接波の各成分からなる波動振幅ベクトルａ₀と、前記３次元はり構造の不連続部の透過反射性を表す散乱行列Ｓと、各成分の伝播に伴った位相の変化を表す伝播行列Ｄとに基づいて、以下の式に従って、応答としての波動振幅を表す、各成分の波動振幅ベクトルａを計算する計算手段として機能させるためのプログラム。
ただし、Ｉは、単位行列である。