(57)【特許請求の範囲】
【請求項１】
複数のはり要素が結合した３次元はり構造への入力振幅に対する応答としての波動振幅を計算するシミュレーション装置であって、
入力振幅を表す、はり要素の軸方向であるｘ方向の縦波前進伝播波、ｘ方向の縦波後退伝播波、ねじり前進伝播波、ねじり後退伝播波、ｙ方向の曲げ前進伝播波、ｙ方向の曲げ前進近接波、ｙ方向の曲げ後退伝播波、ｙ方向の曲げ後退近接波、ｚ方向の曲げ前進伝播波、ｚ方向の曲げ前進近接波、ｚ方向の曲げ後退伝播波、及びｚ方向の曲げ後退近接波の各成分からなる波動振幅ベクトルａ0と、前記３次元はり構造の不連続部の透過反射性を表す散乱行列Ｓと、各成分の伝播に伴った位相の変化を表す伝播行列Ｄとに基づいて、以下の式に従って、応答としての波動振幅を表す、各成分の波動振幅ベクトルａを計算する計算手段
を含むシミュレーション装置。
【数１】
0006664690.tif 000046
ただし、Ｉは、単位行列である。
【請求項２】
前記計算手段は、周波数毎に、各成分の波動振幅ベクトルａを計算し、各成分について、周波数応答関数を求める請求項１記載のシミュレーション装置。
【請求項３】
前記計算手段は、
前記周波数毎に、前記複数のはり要素の各々について、前記はり要素の縦弾性係数、体積質量密度、横弾性係数、断面二次極モーメント、ねじり定数、及び断面積に基づいて、ｘ方向の縦波の波数ｋL、ねじり波の波数ｋT、ｙ方向の曲げ波の波数ｋBy、及びｚ方向の曲げ波の波数ｋBzを計算する波数計算手段と、
前記周波数毎に、前記複数のはり要素の各々について、前記ｘ方向の縦波の波数ｋL、前記ねじり波の波数ｋT、前記ｙ方向の曲げ波の波数ｋBy、及び前記ｚ方向の曲げ波の波数ｋBzに基づいて、前進伝播波成分及び後退伝播波成分の各々についての変位波動振幅変換行列Ψ1、Ψ2と、前進波及び後退波の各々についての伝播距離ｘによる位相変化を表す伝播行列Ｇ1（ｘ）、Ｇ2（ｘ）と、前進伝播波成分及び後退伝播波成分の各々についての内力波動振幅変換行列Φ1、Φ2とを計算する変換行列伝播行列計算手段と、
前記はり要素の各々について、前記はり要素の全体座標系の基底ベクトルｅi及び要素座標系の基底ベクトルＥiに基づいて、不連続部において結合している前記はり要素の各々についての座標変換行列Ｃ(i)を計算する座標変換行列計算手段と、
前記はり要素の各々について、前記はり要素の前進伝播波の、前記不連続部への入射方向ｄ(i)を計算する入射方向計算手段と、
前記周波数毎に、前記はり要素の各々についての、前記変位波動振幅変換行列Ψ1、Ψ2と、前記伝播行列Ｇ1（ｘ）、Ｇ2（ｘ）と、前記内力波動振幅変換行列Φ1、Φ2と、前記座標変換行列Ｃ(i)と、前記入射方向ｄ(i)とに基づいて、前記はり要素の各々についての反射係数ｒと、不連続部において結合している前記はり要素のペアの各々について透過係数ｔとを計算する係数計算手段と、
前記周波数毎に、前記はり要素の各々についての反射係数ｒと、前記はり要素のペアの各々について透過係数ｔとに基づいて、前記３次元はり構造の不連続部の透過反射性を表す散乱行列Ｓを計算する散乱行列計算手段と、
前記周波数毎に、前記はり要素の各々についての前記伝播行列Ｇ1（ｘ）、Ｇ2（ｘ）と、前記はり要素の各々の導波路長さＬ(i)とに基づいて、前記はり要素の各々の位相情報を表す伝播行列Ｄを計算する伝播行列計算手段と、
前記周波数毎に、前記波動振幅ベクトルａ0と、前記散乱行列Ｓと、前記伝播行列Ｄとに基づいて、各成分の波動振幅ベクトルａを計算し、各成分について、周波数応答関数を求める周波数応答関数計算手段とを含む請求項２記載のシミュレーション装置。
【請求項４】
複数のはり要素が結合した３次元はり構造への入力振幅に対する応答としての波動振幅を計算するためのプログラムであって、
コンピュータを、
入力振幅を表す、はり要素の軸方向であるｘ方向の縦波前進伝播波、ｘ方向の縦波後退伝播波、ねじり前進伝播波、ねじり後退伝播波、ｙ方向の曲げ前進伝播波、ｙ方向の曲げ前進近接波、ｙ方向の曲げ後退伝播波、ｙ方向の曲げ後退近接波、ｚ方向の曲げ前進伝播波、ｚ方向の曲げ前進近接波、ｚ方向の曲げ後退伝播波、及びｚ方向の曲げ後退近接波の各成分からなる波動振幅ベクトルａ0と、前記３次元はり構造の不連続部の透過反射性を表す散乱行列Ｓと、各成分の伝播に伴った位相の変化を表す伝播行列Ｄとに基づいて、以下の式に従って、応答としての波動振幅を表す、各成分の波動振幅ベクトルａを計算する計算手段として機能させるためのプログラム。
【数２】
0006664690.tif 000047
ただし、Ｉは、単位行列である。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、シミュレーション装置及びプログラムに係り、特に、３次元はり構造への入力振幅に対する応答としての波動振幅を計算するシミュレーション装置及びプログラムに関する。
【背景技術】
【０００２】
従来より、数値解析によって周波数応答関数を求める手法として、質量行列、剛性行列を利用した２つの手法が知られている（非特許文献１）。１つ目の手法は、直接法と呼ばれ、質量行列と剛性行列とを組み合わせた行列の逆行列を直接求めることにより、周波数応答関数が求まる。もうひとつの手法は、理論モード解析と呼ばれ、質量行列と剛性行列とを組み合わせた行列を固有値解析にかけることにより、モード特性を求め、各モードの応答の重ね合わせで周波数応答関数を求める手法である。この手法では周波数応答関数における各モードの寄与を評価することが可能である。
【０００３】
また、分布定数系モデルで複雑な構造物の波動伝播解析を実施している事例として、宇宙構造物を対象にした時系列の波動伝播数値解析の事例が知られている（非特許文献２）。
【先行技術文献】
【非特許文献】
【０００４】
【非特許文献１】長松昭男、「モード解析入門」、コロナ社、１９９３年．
【非特許文献２】西田明美、「３次元フレーム構造物の波動伝播特性に関する研究：ティモシェンコ梁理論の導入」、構造工学論文集Ｂ、Ｖｏｌ．５２Ｂ、pp.119-124、日本建築学会、２００６年
【発明の概要】
【発明が解決しようとする課題】
【０００５】
上記非特許文献１の周波数応答関数を求める解析技術では、振動を構成している構造物内部の各波動の伝播を物理的に考慮しておらず、波動伝播を考慮しない別の物理モデル（離散系モデル）から周波数応答関数を求めている。
【０００６】
また、上記非特許文献２に記載の、波動伝播に関する研究・技術開発では、衝撃荷重負荷時の時系列での波動伝播解析が主であり、周波数応答関数などの定常状態での波動伝播の扱いに関する研究事例はない。
【０００７】
このように、従来法のモード解析や周波数応答関数から波動伝播寄与を求めることができない。構造物の波動伝播寄与を求めるには別の手法がいる。
【０００８】
本発明は、上記の事情を鑑みてなされたもので、波動振幅の成分別に周波数応答関数を求めることができるシミュレーション装置及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【０００９】
上記の目的を達成するために本発明に係るシミュレーション装置は、複数のはり要素が結合した３次元はり構造への入力振幅に対する応答としての波動振幅を計算するシミュレーション装置であって、入力振幅を表す、はり要素の軸方向であるｘ方向の縦波前進伝播波、ｘ方向の縦波後退伝播波、ねじり前進伝播波、ねじり後退伝播波、ｙ方向の曲げ前進伝播波、ｙ方向の曲げ前進近接波、ｙ方向の曲げ後退伝播波、ｙ方向の曲げ後退近接波、ｚ方向の曲げ前進伝播波、ｚ方向の曲げ前進近接波、ｚ方向の曲げ後退伝播波、及びｚ方向の曲げ後退近接波の各成分からなる波動振幅ベクトルａ0と、前記３次元はり構造の不連続部の透過反射性を表す散乱行列Ｓと、各成分の伝播に伴った位相の変化を表す伝播行列Ｄとに基づいて、以下の式に従って、応答としての波動振幅を表す、各成分の波動振幅ベクトルａを計算する計算手段を含んで構成されている。
【００１０】
【数１】
0006664690.tif 000002
ただし、Ｉは、単位行列である。
【００１１】
本発明に係るプログラムは、複数のはり要素が結合した３次元はり構造への入力振幅に対する応答としての波動振幅を計算するためのプログラムであって、コンピュータを、入力振幅を表す、はり要素の軸方向であるｘ方向の縦波前進伝播波、ｘ方向の縦波後退伝播波、ねじり前進伝播波、ねじり後退伝播波、ｙ方向の曲げ前進伝播波、ｙ方向の曲げ前進近接波、ｙ方向の曲げ後退伝播波、ｙ方向の曲げ後退近接波、ｚ方向の曲げ前進伝播波、ｚ方向の曲げ前進近接波、ｚ方向の曲げ後退伝播波、及びｚ方向の曲げ後退近接波の各成分からなる波動振幅ベクトルａ0と、前記３次元はり構造の不連続部の透過反射性を表す散乱行列Ｓと、各成分の伝播に伴った位相の変化を表す伝播行列Ｄとに基づいて、上記の式に従って、応答としての波動振幅を表す、各成分の波動振幅ベクトルａを計算する計算手段として機能させるためのプログラムである。
【００１２】
本発明によれば、計算手段が、入力振幅を表す、はり要素の軸方向であるｘ方向の縦波前進伝播波、ｘ方向の縦波後退伝播波、ねじり前進伝播波、ねじり後退伝播波、ｙ方向の曲げ前進伝播波、ｙ方向の曲げ前進近接波、ｙ方向の曲げ後退伝播波、ｙ方向の曲げ後退近接波、ｚ方向の曲げ前進伝播波、ｚ方向の曲げ前進近接波、ｚ方向の曲げ後退伝播波、及びｚ方向の曲げ後退近接波の各成分からなる波動振幅ベクトルａ0と、前記３次元はり構造の不連続部の透過反射性を表す散乱行列Ｓと、各成分の伝播に伴った位相の変化を表す伝播行列Ｄとに基づいて、上記の式に従って、応答としての波動振幅を表す、各成分の波動振幅ベクトルａを計算する。
【００１３】
このように、３次元はり構造への入力振幅に対する応答として、ｘ方向の縦波前進伝播波、ｘ方向の縦波後退伝播波、ねじり前進伝播波、ねじり後退伝播波、ｙ方向の曲げ前進伝播波、ｙ方向の曲げ前進近接波、ｙ方向の曲げ後退伝播波、ｙ方向の曲げ後退近接波、ｚ方向の曲げ前進伝播波、ｚ方向の曲げ前進近接波、ｚ方向の曲げ後退伝播波、及びｚ方向の曲げ後退近接波の各成分からなる波動振幅ベクトルを計算することにより、波動振幅の成分別に周波数応答関数を求めることができる。
【００１４】
上記の計算手段は、周波数毎に、各成分の波動振幅ベクトルａを計算し、各成分について、周波数応答関数を求めることができる。
【００１５】
上記の計算手段は、前記周波数毎に、前記複数のはり要素の各々について、前記はり要素の縦弾性係数、体積質量密度、横弾性係数、断面二次極モーメント、ねじり定数、及び断面積に基づいて、ｘ方向の縦波の波数ｋL、ねじり波の波数ｋT、ｙ方向の曲げ波の波数ｋBy、及びｚ方向の曲げ波の波数ｋBzを計算する波数計算手段と、前記周波数毎に、前記複数のはり要素の各々について、前記ｘ方向の縦波の波数ｋL、前記ねじり波の波数ｋT、前記ｙ方向の曲げ波の波数ｋBy、及び前記ｚ方向の曲げ波の波数ｋBzに基づいて、前進伝播波成分及び後退伝播波成分の各々についての変位波動振幅変換行列Ψ1、Ψ2と、前進波及び後退波の各々についての伝播距離ｘによる位相変化を表す伝播行列Ｇ1（ｘ）、Ｇ2（ｘ）と、前進伝播波成分及び後退伝播波成分の各々についての内力波動振幅変換行列Φ1、Φ2とを計算する変換行列伝播行列計算手段と、前記はり要素の各々について、前記はり要素の全体座標系の基底ベクトルｅi及び要素座標系の基底ベクトルＥiに基づいて、不連続部において結合している前記はり要素の各々についての座標変換行列Ｃ(i)を計算する座標変換行列計算手段と、前記はり要素の各々について、前記はり要素の前進伝播波の、前記不連続部への入射方向ｄ(i)を計算する入射方向計算手段と、前記周波数毎に、前記はり要素の各々についての、前記変位波動振幅変換行列Ψ1、Ψ2と、前記伝播行列Ｇ1（ｘ）、Ｇ2（ｘ）と、前記内力波動振幅変換行列Φ1、Φ2と、前記座標変換行列Ｃ(i)と、前記入射方向ｄ(i)とに基づいて、前記はり要素の各々についての反射係数ｒと、不連続部において結合している前記はり要素のペアの各々について透過係数ｔとを計算する係数計算手段と、前記周波数毎に、前記はり要素の各々についての反射係数ｒと、前記はり要素のペアの各々について透過係数ｔとに基づいて、前記３次元はり構造の不連続部の透過反射性を表す散乱行列Ｓを計算する散乱行列計算手段と、前記周波数毎に、前記はり要素の各々についての前記伝播行列Ｇ1（ｘ）、Ｇ2（ｘ）と、前記はり要素の各々の導波路長さＬ(i)とに基づいて、前記はり要素の各々の位相情報を表す伝播行列Ｄを計算する伝播行列計算手段と、前記周波数毎に、前記波動振幅ベクトルａ0と、前記散乱行列Ｓと、前記伝播行列Ｄとに基づいて、各成分の波動振幅ベクトルａを計算し、各成分について、周波数応答関数を求める周波数応答関数計算手段とを含むことができる。
【発明の効果】
【００１６】
以上説明したように、本発明のシミュレーション装置及びプログラムによれば、３次元はり構造への入力振幅に対する応答として、ｘ方向の縦波前進伝播波、ｘ方向の縦波後退伝播波、ねじり前進伝播波、ねじり後退伝播波、ｙ方向の曲げ前進伝播波、ｙ方向の曲げ前進近接波、ｙ方向の曲げ後退伝播波、ｙ方向の曲げ後退近接波、ｚ方向の曲げ前進伝播波、ｚ方向の曲げ前進近接波、ｚ方向の曲げ後退伝播波、及びｚ方向の曲げ後退近接波の各成分からなる波動振幅ベクトルを計算することにより、波動振幅の成分別に周波数応答関数を求めることができる、という効果が得られる。
【図面の簡単な説明】
【００１７】
【図１】本発明の実施の形態に係るシミュレーション装置を示すブロック図である。
【図２】本発明の実施の形態に係るシミュレーション装置のシミュレーション処理ルーチンの内容を示すフローチャートである。
【図３】３次元はり構造の一例を示す図である。
【図４】周波数応答関数の一例を示す図である。
【図５】周波数応答関数の一例を示す図である。
【発明を実施するための形態】
【００１８】
以下、図面を参照して本発明の実施の形態を詳細に説明する。
【００１９】
＜本発明の実施の形態の概要＞
透過・反射係数と形状情報（導波路長さ等）と波動振幅とでレイトレースモデルと呼ばれるモデルを構成できる。このモデルは波動が構造物の不連続部から別の不連続部まで伝播する際の伝播前後の波動振幅の関係を数値的に記述したモデルである。入力が１Ｎの場合の波動振幅を計算し、波動伝播が無限回起きた時点のある任意位置での各種の波動振幅を個別にみることにより、周波数応答関数における各種波動振幅寄与を評価する。
【００２０】
＜本発明の実施の形態の原理＞
次に、波動伝播による波動振幅解析の原理について説明する。
【００２１】
＜3次元はり要素における波動振幅表現＞
まず、3次元のはり要素における波動理論を整理する。縦波、ねじり波、曲げ波についての波動振幅と変位および内力の関係を整理し、透過・反射係数の算出に必要な波動振幅から変位および内力への変換行列を導出する。
【００２２】
＜縦波についての波動理論＞
はりを伝播する縦波についての運動方程式は次式となる。ここでρは体積質量密度[kg/m3]、Aは断面積[m2]、Eは縦弾性係数[N/m2]を示す。
【００２３】
【数２】
0006664690.tif 000003
【００２４】
上式の運動方程式の一般解は、次式の通りである。
【００２５】
【数３】
0006664690.tif 000004
【００２６】
ここで、a1Lはxの正の方向へ伝播する前進伝播波の波動振幅、a2Lはxの負の方向へ伝播する後退伝播波の波動振幅を表す。また、内力F(x)は、次式より求まる。
【００２７】
【数４】
0006664690.tif 000005
【００２８】
kLは、縦波の波数[rad/m]であり、次式で定義される。
【００２９】
【数５】
0006664690.tif 000006
【００３０】
＜ねじり波についての波動理論＞
はり要素の軸方向をx軸とし、断面のx軸周りの変位をs(x,t)とする。また、ねじり波伝播のねじりモーメントをh(x,t)とする。ここでGは横弾性係数[N/m2]、IPは断面二次極モーメント[m4]、Jはねじり定数を示す。線形振動領域では、断面変形が小さく、ねじり波についての運動方程式は次式で表される。
【００３１】
【数６】
0006664690.tif 000007
【００３２】
上式の運動方程式の一般解は、次式の通りである。
【００３３】
【数７】
0006664690.tif 000008
【００３４】
また、内力H(x)は、次式より求まる。
【００３５】
【数８】
0006664690.tif 000009
【００３６】
ここで、a1Tはxの正の方向へ伝播するねじり前進伝播波の振幅、a2Tはxの負の方向へ伝播するねじり後退伝播波の振幅を表す。このとき、kTはねじり波の波数[rad/m]である。
【００３７】
【数９】
0006664690.tif 000010
【００３８】
＜曲げ波についての波動理論＞
曲げ波に関する運動方程式は次式で表すことができる。
【００３９】
【数１０】
0006664690.tif 000011
【００４０】
ここで、Iは断面二次極モーメント[m4]である。本式の運動方程式の解は、次式で得られる。
【００４１】
【数１１】
0006664690.tif 000012
【００４２】
ここで、a1Pは曲げ前進伝播波の振幅、a2Pは曲げ後退伝播波の振幅、a1Nは曲げ前進近接波の振幅、a2Nは曲げ後退近接波の振幅を表す。ここで、曲げ波の波数kBは次式の通りである。
【００４３】
【数１２】
0006664690.tif 000013
【００４４】
このとき、たわみ角θ(x,t)とたわみv(x,t)との関係は次式である。
【００４５】
【数１３】
0006664690.tif 000014
【００４６】
また、モーメントm(x,t)とせん断力w(x,t)は、それぞれ次式で求まる。
【００４７】
【数１４】
0006664690.tif 000015
【００４８】
＜波動振幅変換行列＞
はりの反射・透過係数を求める際に、はりの変位、内力を波動振幅で表すための波動振幅変換行列が必要である。次に、波動振幅変換行列の定義を述べる。
【００４９】
まず、波動振幅ベクトルの前進伝播波成分a1および後退伝播波成分a2を次式の通り定義する。
【００５０】
【数１５】
0006664690.tif 000016
【００５１】
縦波による変位の一般解である式(2)、ねじり波による変位の一般解である式(6)、曲げ波についての変位の一般解である式(10)、式(12)を、前進波a1による変位成分および後退波a2による変位成分に分ける。縦波の変位をU1、U2、 y軸方向の曲げ波のたわみをV1y、V2y、z軸方向の曲げ波のたわみをV1z、V2z、ねじり波の変位をS1、S2、y軸方向のたわみを角Θ1y、Θ2y、z軸方向のたわみ角をΘ1z、Θ2zとする。ここで、変位の一般解についてのベクトルをU、前進波および後退波による変位ベクトルをU1、U2として次式で定義する。
【００５２】
【数１６】
0006664690.tif 000017
【００５３】
これらの定義より、式(2)、式(6)、式(10)、式(12)を用いて、波動振幅の関係をまとめると、次式で表すことができる。
【００５４】
【数１７】
0006664690.tif 000018
【００５５】
このとき、Ψ1、 Ψ2はそれぞれ波動振幅の前進伝播波成分および後退伝播波成分についての変位-波動振幅変換行列であり、それぞれ次の２つの式で表すことができる。
【００５６】
【数１８】
0006664690.tif 000019
【００５７】
また、G1(x)、 G2(x)はそれぞれ、前進波および後退波の伝播距離xによる位相変化をあらわす伝播行列であり、式(20)および式(21)で表される。
【００５８】
【数１９】
0006664690.tif 000020
【００５９】
次に、縦波による内力の解である式(3)、ねじり波による内力の解である式(7)、曲げ波による内力の解である式(13)、式(14)を、前進波および後退波による成分に分ける。縦波の内力F1、F2、y軸方向の曲げ波のせん断力W1y、W2y、z軸方向の曲げ波のせん断力W1z、W2z、ねじり波の内力H1、H2、y軸方向の曲げモーメントM1y、M2y、z軸方向の曲げモーメントM1z、M2zとする。ここで、内力の一般解についてのベクトルをF、前進波および後退波による内力ベクトルをF1、F2として次式で定義する。
【００６０】
【数２０】
0006664690.tif 000021
【００６１】
これらの定義より、式(3)、式(7)、式(13)、式(14)をまとめると、内力と波動振幅の関係をまとめると、次式の通り表すことができる。
【００６２】
【数２１】
0006664690.tif 000022
【００６３】
このときの行列Φ1、Φ2はそれぞれ波動振幅の前進波成分および後退波成分についての内力−波動振幅変換行列であり、それぞれ次の２つの式で表される。ここで、Iyはy軸方向の断面二次極モーメント[m4]、Izはz軸方向の断面二次極モーメント[m4]を示す。
【００６４】
【数２２】
0006664690.tif 000023
【００６５】
以上より、縦波、ねじり波、曲げ波の一般解を整理し、波動振幅から変位および内力への変換行列を導出した。
【００６６】
＜不連続部の透過・反射係数の算出＞
次に、Nc個の要素が結合する不連続部における反射・透過係数の算出について定式化する。a1(i)、 a2(i)は、それぞれ結合要素i ( = 1〜Nc)の要素座標系のx軸の正方向に伝播する前進波と、負の向きに伝播する後退波についての波動振幅ベクトルである。
【００６７】
不連続部において結合するはり要素i ( = 1〜Nc)の要素座標系から全体座標系への座標変換行列C(i)は、全体座標系の基底ベクトルをe、要素iの要素座標系の基底ベクトルをE(i)として、次式で表される。
【００６８】
【数２３】
0006664690.tif 000024
【００６９】
式(17)において不連続部位置をx=0とすると、不連続部で結合するはり要素iの不連続部での変位U(i)は、次式で表される。
【００７０】
【数２４】
0006664690.tif 000025
【００７１】
不連続部における変位の連続性より、次式が成り立つ。
【００７２】
【数２５】
0006664690.tif 000026
【００７３】
また、式(23)において不連続部の位置をx = 0とすると、不連続部で結合するはり要素iの不連続部での内力F(i)は、次式で表される。
【００７４】
【数２６】
0006664690.tif 000027
【００７５】
不連続部における力の釣り合いより、次式が成立する。
【００７６】
【数２７】
0006664690.tif 000028
【００７７】
本式におけるd(i)は要素iの前進伝播波の不連続部への入射方向の正負を表しており、次式であらわされる。
【００７８】
【数２８】
0006664690.tif 000029
【００７９】
ここで、wI(i)およびwO(i)は、不連続部に対する入射波および反射波を定義する変数であり、式(32)、式(33)で表される。なお、xcは、不連続部の座標、x1(i)およびx2(i)は、要素iの構成節点座標であり、x1(i)からx2(i)の方向が要素座標系のx軸の正方向である。wI(i)は、不連続部に対する入射する波動振幅ベクトルが要素座標系のx軸の正方向に伝播する前進波成分であるか負の方向に伝播する後退波成分であるかを定める変数であり、1のとき要素iの前進波が不連続部に入射し、2のとき後退波が入射するように要素座標系が定義されていることになる。一方、wO(i)は、不連続部から反射する波動振幅ベクトルを定める変数で、1のとき前進波が不連続部から反射し、2のとき後退波が不連続部から反射するように、要素座標系が定義されていることになる。
【００８０】
【数２９】
0006664690.tif 000030
【００８１】
不連続部における変位の連続を表す式(28)および、力のつり合いを表す式(30)を、左辺に不連続部に波動が入射する波動振幅成分、右辺に不連続部から波動振幅が反射される波動振幅成分に整理し、行列表記すると次式で表される。
【００８２】
【数３０】
0006664690.tif 000031
0006664690.tif 000032
【００８３】
また、透過・反射係数を用いて、不連続部の波動振幅の変化を表すと次式の通りである。
【００８４】
【数３１】
0006664690.tif 000033
【００８５】
したがって、式(34)および式(35)より、(34)反射係数rおよび透過係数tは、変位−波動振幅変換行列Ψ、内力−波動振幅変換行列Φ、座標変換行列Cを用いて、次式で算出することが可能である。
【００８６】
【数３２】
0006664690.tif 000034
【００８７】
＜レイトレースモデルの一般化＞
次にN個の要素から構成される解析対象のレイトレースモデルを定式化する。はり要素i (=1〜N)を伝播する前進波の波動振幅a1(i)および後退波の波動振幅成分a2(i)からなる波動振幅ベクトルa(i)を次式の通り定義する。
【００８８】
【数３３】
0006664690.tif 000035
【００８９】
不連続部にて結合されるはり要素iとはり要素jの間の透過・反射の関係について定式化すると、次式のように表すことができる。
【００９０】
【数３４】
0006664690.tif 000036
【００９１】
ここで、R(i)は、不連続部におけるはり要素iから入射した波動に対する反射係数行列、T(i,j)ははり要素iからはり要素jに波動が透過する際の透過係数行列である。透過係数行列Tおよび反射係数行列Rは、前節の式(36)で算出された透過反射係数と、式(32)で示したwIおよびwOを用いて表すと、不連続部に対する波動ベクトルの入射および反射方向の定義を考慮し、それぞれ、式(39)および式(40)で表される。
【００９２】
【数３５】
0006664690.tif 000037
【００９３】
以上より、透過・反射の関係が定式化される。系を構成するはり要素iについての反射係数R(i)および、不連続部にて結合されているはり要素iとはり要素jの間の透過係数行列T(i,j)を使用して、次式のような系全体の散乱行列Sを構築することが可能である。
【００９４】
【数３６】
0006664690.tif 000038
【００９５】
ここで、系全体の要素についての波動振幅ベクトルを並べたaを、次式の通り定義する。
【００９６】
【数３７】
0006664690.tif 000039
【００９７】
また、波動が不連続部から発生し、はり要素内をもう一方の不連続部まで伝播することを表す伝播行列Dを次式の通り定義する。
【００９８】
【数３８】
0006664690.tif 000040
【００９９】
ある波動振幅ベクトルanがはり要素内を伝播し、不連続部に到達し、次に伝播する波動振幅ベクトルan+1に変化する波動モデルをレイトレースモデルと呼び、伝播行列Dと散乱行列Sを使用して次式で表すことができる。
【０１００】
【数３９】
0006664690.tif 000041
【０１０１】
レイトレースモデルを用いると、初期波動伝播前後の波動振幅ａ1の関係は以下のようになる
【０１０２】
【数４０】
0006664690.tif 000042
【０１０３】
定常状態の波動振幅ａを考えると、複数回伝播した波動振幅(ａ0，ａ1，ａ2…ａ∞)が混在する状況であり以下のように定式化できる。
【０１０４】
【数４１】
0006664690.tif 000043
【０１０５】
これをまとめると以下の式になる。
【０１０６】
【数４２】
0006664690.tif 000044
【０１０７】
ａ0が構造物の任意箇所に１Nの振幅の加振力を付加された時の初期波動振幅だとすると、上記のａは周波数応答関数を各種波動振幅別に分解した波動振幅ベクトルであり、周波数応答関数の波動振幅寄与を表現している。よって、上式を利用すれば周波数応答関数の波動振幅寄与を求めることが可能である。
【０１０８】
＜システム構成＞
図１に示すように、本発明の実施の形態に係るシミュレーション装置１０は、ＣＰＵ１２、ＲＯＭ１４、ＲＡＭ１６、ＨＤＤ１８、通信インタフェース２０、及びこれらを相互に接続するためのバス２２を備えている。
【０１０９】
ＣＰＵ１２は、各種プログラムを実行する。ＲＯＭ１４には、各種プログラムやパラメータ等が記憶されている。ＲＡＭ１６は、ＣＰＵ１２による各種プログラムの実行時におけるワークエリア等として用いられる。記録媒体としてのＨＤＤ１８には、後述する微分計算処理ルーチンを実行するためのプログラムを含む各種プログラムや各種データが記憶されている。
【０１１０】
本実施の形態におけるシミュレーション装置１０は、入力として、入力振幅を表す、はり要素の軸方向であるｘ方向の縦波前進伝播波、ｘ方向の縦波後退伝播波、ねじり前進伝播波、ねじり後退伝播波、ｙ方向の曲げ前進伝播波、ｙ方向の曲げ前進近接波、ｙ方向の曲げ後退伝播波、ｙ方向の曲げ後退近接波、ｚ方向の曲げ前進伝播波、ｚ方向の曲げ前進近接波、ｚ方向の曲げ後退伝播波、及びｚ方向の曲げ後退近接波の各成分からなる波動振幅ベクトルａ0を受け付ける。また、シミュレーション装置１０は、複数のはり要素が結合した３次元はり構造の形状情報として、各はり要素についての縦弾性係数Ｅ、体積質量密度ρ、横弾性係数Ｇ、断面二次極モーメントＩp、ねじり定数Ｊ、断面積Ａ、全体座標系の基底ベクトルｅ、要素座標系の基底ベクトルＥ(i)、及び導波路長さＬ(i)を受け付ける。
【０１１１】
シミュレーション装置１０は、周波数毎に、波動振幅ベクトルａ0と、３次元はり構造の不連続部の透過反射性を表す散乱行列Ｓと、はり要素の各々の位相情報を表す伝播行列Ｄとに基づいて、上記式（４７）に従って、応答としての波動振幅を表す、各成分からなる波動振幅ベクトルａを計算し、各成分について、周波数応答関数を求める。
【０１１２】
＜シミュレーション装置の動作＞
次に、本実施の形態に係るシミュレーション装置１０の作用について説明する。
【０１１３】
まず、シミュレーション装置１０に、波動振幅ベクトルａ0と、上述した３次元はり構造の形状情報が入力されると、シミュレーション装置１０によって、３次元はり構造の形状モデルが生成される。そして、シミュレーション装置１０によって、図２に示すシミュレーション処理ルーチンが実行される。
【０１１４】
まず、ステップＳ１００において、周波数ω毎に、各はり要素の縦弾性係数Ｅ、体積質量密度ρ、横弾性係数Ｇ、断面二次極モーメントＩp、ねじり定数Ｊ、及び断面積Ａに基づいて、上記式（４）、（８）、（１１）に従って、各はり要素について、ｘ方向の縦波の波数ｋL、ねじり波の波数ｋT、ｙ方向の曲げ波の波数ｋBy、及びｚ方向の曲げ波の波数ｋBzを計算する。
【０１１５】
ステップＳ１０２では、周波数ω毎に、各はり要素について、ｘ方向の縦波の波数ｋL、ねじり波の波数ｋT、ｙ方向の曲げ波の波数ｋBy、及びｚ方向の曲げ波の波数ｋBzに基づいて、上記式（１８）〜（２１）、（２４）、（２５）に従って、前進伝播波成分及び後退伝播波成分の各々についての変位波動振幅変換行列Ψ1、Ψ2と、前進波及び後退波の各々についての伝播距離ｘによる位相変化を表す伝播行列Ｇ1（ｘ）、Ｇ2（ｘ）と、前進伝播波成分及び後退伝播波成分の各々についての内力波動振幅変換行列Φ1、Φ2とを計算する。
【０１１６】
ステップＳ１０４では、はり要素の各々についての当該はり要素の全体座標系の基底ベクトルｅ及び要素座標系の基底ベクトルＥ(i)に基づいて、上記式（２６）に従って、不連続部において結合しているはり要素の各々についての座標変換行列Ｃ(i)を計算する。
【０１１７】
また、はり要素の各々について、３次元はり構造の形状情報に基づいて、上記式（３１）に従って、当該はり要素の前進伝播波の、前記不連続部への入射方向ｄ(i)を計算する。
【０１１８】
ステップＳ１０６では、周波数ω毎に、上記ステップＳ１０２で計算された、前進伝播波成分及び後退伝播波成分の各々についての変位波動振幅変換行列Ψ1、Ψ2と、前進波及び後退波の各々についての伝播距離ｘによる位相変化を表す伝播行列Ｇ1（ｘ）、Ｇ2（ｘ）と、前進伝播波成分及び後退伝播波成分の各々についての内力波動振幅変換行列Φ1、Φ2と、上記ステップＳ１０４で計算された、はり要素の各々についての前記座標変換行列Ｃ(i)と、前記はり要素の各々についての前記入射方向ｄ(i)とに基づいて、上記式（３６）に従って、はり要素の各々についての反射係数ｒと、不連続部において結合しているはり要素のペアの各々について透過係数ｔとを計算する。
【０１１９】
そして、ステップＳ１０８では、周波数ω毎に、上記ステップＳ１０６で計算された、はり要素の各々についての反射係数ｒと、はり要素のペアの各々について透過係数ｔとに基づいて、上記式（４１）に従って、３次元はり構造の不連続部の透過反射性を表す散乱行列Ｓを計算する。
【０１２０】
ステップＳ１１０では、周波数ω毎に、上記ステップＳ１０２で計算された前進波及び後退波の各々についての伝播距離ｘによる位相変化を表す伝播行列Ｇ1（ｘ）、Ｇ2（ｘ）と、はり要素の各々の導波路長さＬ(i)とに基づいて、上記式（４３）に従って、はり要素の各々の位相情報を表す伝播行列Ｄを計算する。
【０１２１】
そして、ステップＳ１１２では、周波数ω毎に、波動振幅ベクトルａ0と、上記ステップＳ１０８で計算された散乱行列Ｓと、上記ステップＳ１１０で計算された伝播行列Ｄとに基づいて、波動振幅ベクトルａを計算する。次のステップＳ１１４では、上記ステップＳ１１２で周波数ω毎に計算された波動振幅ベクトルａに基づいて、各成分について、周波数応答関数を求め、出力部（図示省略）により出力して、シミュレーション処理ルーチンを終了する。
【０１２２】
＜実施例＞
断面特性が異なる３つのはり要素を結合した３次元はり構造のモデルで波動振幅寄与の例を示す。
【０１２３】
図３に示すように、３次元はり構造モデルに対して、入力波動振幅として、拘束条件なしで左端に１Ｎの調和荷重を負荷している。評価点は加振点（１）とはり中央部（２）である。はり要素の諸元を以下の表１に記す。
【０１２４】
【表１】
0006664690.tif 000045
【０１２５】
解析結果を図４、図５に示す。３次元はり構造のすべての成分の波動の波動振幅寄与を表示すると複雑になるため、ここでは低周波の上下方向曲げ振動に着目し、前進伝播波、前進近接波、後退伝播波、後退近接波の周波数応答関数に対する寄与を表示する。図４は、評価点（１）における各成分の周波数応答関数を示しており、図５は、評価点（２）における各成分の周波数応答関数を示している。
【０１２６】
図４、図５のグラフで示すように、周波数応答関数に対して各種成分の寄与を解析することが可能である。
【０１２７】
以上説明したように、本発明の実施の形態に係るシミュレーション装置によれば、３次元はり構造への入力振幅に対する応答として、ｘ方向の縦波前進伝播波、ｘ方向の縦波後退伝播波、ねじり前進伝播波、ねじり後退伝播波、ｙ方向の曲げ前進伝播波、ｙ方向の曲げ前進近接波、ｙ方向の曲げ後退伝播波、ｙ方向の曲げ後退近接波、ｚ方向の曲げ前進伝播波、ｚ方向の曲げ前進近接波、ｚ方向の曲げ後退伝播波、及びｚ方向の曲げ後退近接波の各成分からなる波動振幅ベクトルを計算することにより、波動振幅の成分別に周波数応答関数を求めることができる。
【０１２８】
また、周波数応答関数を構造物内の波動振幅成分の寄与として評価することで、振動応答の主原因となる波動を特定でき、効果的な対策を検討することが可能となる。
【０１２９】
また、位相の異なる各種波動振幅の重ね合わせにより、任意位置の振動応答を低減する構造を提案することなども可能となる。
【０１３０】
なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【０１３１】
本発明のプログラムは、記憶媒体に格納して提供するようにしてもよい。
【符号の説明】
【０１３２】
１０ シミュレーション装置
１２ ＣＰＵ
１４ ＲＯＭ
１６ ＲＡＭ
１８ ＨＤＤ
【図１】
0006664690.tif 000048
【図２】
0006664690.tif 000049
【図３】
0006664690.tif 000050
【図４】
0006664690.tif 000051
【図５】
0006664690.tif 000052