振動解析方法及び振動解析装置
【課題】実際に測定を行う測定点を減少させても、同定精度をほとんど低下させることなく剛体特性を同定することができる振動解析方法又は振動解析装置を提供する。
【解決手段】加振試験で測定点において実際に測定される測定自由度よりも同定する特性行列の自由度を大きくし、これにより測定データに基づいて算出されるモード特性の値の数を特性行列に基づいて算出されるモード特性の値の数よりも少なくする。特性行列を同定する際に、特性行列に基づいて算出されるモード特性のうち測定データに基づいて算出されたモード特性に対応するモード特性についてのみ、特性行列に基づいて算出されるモード特性の値が測定データに基づいて算出されたモード特性の値に一致するように特性行列の各成分の値を変更していくことで特性行列を同定する。
【解決手段】加振試験で測定点において実際に測定される測定自由度よりも同定する特性行列の自由度を大きくし、これにより測定データに基づいて算出されるモード特性の値の数を特性行列に基づいて算出されるモード特性の値の数よりも少なくする。特性行列を同定する際に、特性行列に基づいて算出されるモード特性のうち測定データに基づいて算出されたモード特性に対応するモード特性についてのみ、特性行列に基づいて算出されるモード特性の値が測定データに基づいて算出されたモード特性の値に一致するように特性行列の各成分の値を変更していくことで特性行列を同定する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、振動解析方法及び振動解析装置に関する。
【背景技術】
【０００２】
自動車やエンジン等の複雑な構造物では、その剛体特性（質量、質量中心位置、三つの主慣性モーメント及びこれらに対応する三本の主軸の向き等）はその構造物の運動性能や振動騒音性能及び振動騒音制御の性能を左右する最大因子の一つである。事実、剛体特性の値は、設計支援のための運動解析、振動解析、防振機構解析、制御性能解析等、多岐にわたる動的挙動の解析、設計、最適化を開始するときに最も必要なパラメータである。従って、複雑な構造物の剛体特性を容易に実用的な精度で同定（推定）することは非常に重要である。
【０００３】
自動車等の構造物の剛体特性を同定する場合には、構造物をワイヤーで吊り下げることによる古典的な剛体特性の計測法を適用することは困難である。このため、現在では振り子ベンチ形式の計測装置が開発され、利用されている。しかし、この計測装置では、自動車を搭載する振り子ベンチ部分及び計測対象の自動車を剛体として扱う必要があるため、セッティングには手間がかかる。また、慣性モーメント計測時には、予め自動車の質量中心位置を同定しておいて、その質量中心が計測装置の回転中心軸線上に正確に位置するように自動車を設置する必要がある。このように、上記振り子ベンチ形式の計測装置には剛体特性を計測するにあたり、いくつかの問題点がある。
【０００４】
本願の発明者は、既に上記課題に対する解決方法として、「実験的特性行列同定法」と名付けた新しい方法を提案している（特許文献１）。この方法では、同定対象物を剛体として扱う必要は無く、むしろ弾性体としての第１次から２次又は３次までの適当な数の共振振動特性から剛体特性を求めることができる。また、加振試験として単点加振多点応答試験を１度実行するだけで、すべての剛体特性を同定することができる。
【０００５】
【特許文献１】特開２００１−３５０７４１号公報
【発明の開示】
【発明が解決しようとする課題】
【０００６】
ところで、上記本願発明者の提案した方法では、加振試験として単点加振多点応答試験が行われる。ここで、複雑な構造物に対して剛体特性を高い精度で同定するためには、その複雑さや大きさに応じて或る程度多くの測定点を設定する必要がある。
【０００７】
ところが、実際に測定を行う測定点の数が増大すると、測定用のセンサのキャリブレーションや計測システムのセッティング等を含む加振試験の測定作業にかかる時間が増大する。また、各測定点毎に周波数応答関数の同定が行われるため、実際に測定を行う測定点の数が増大するとコンピュータによる解析負荷が大きくなる。一方、測定点数を減らして加振試験を行うと、剛体特性の同定精度が低くなってしまう。
【０００８】
従って、上記問題点に鑑み、本発明の目的は、実際に測定を行う測定点を減少させても、同定精度をほとんど低下させることなく剛体特性を同定することができる振動解析方法又は振動解析装置を提供することにある。
【課題を解決するための手段】
【０００９】
上記課題を解決するために、第１の発明では、解析対象物上に設定された複数の測定点の座標データを要求する工程と、上記解析対象物上の少なくとも一点において加振する加振試験が行われた時に得られた各測定点における測定データに基づいてモード特性を算出する工程と、上記座標データに基づいて、上記解析対象物の振動を表す特性行列を構成する成分間の制約条件式を算出する工程と、上記測定データに基づいて算出されたモード特性の値に基づいて上記制約条件式を用いて上記特性行列を同定する工程とを具備する、振動解析方法において、上記加振試験で上記測定点において実際に測定される測定自由度よりも上記特性行列の自由度の方が大きく、よって上記測定データに基づいて算出されるモード特性の値の数の方が上記特性行列に基づいて算出されるモード特性の値の数よりも少なく、上記特性行列を同定する工程では、上記特性行列に基づいて算出されるモード特性のうち上記測定データに基づいて算出されたモード特性に対応するモード特性についてのみ、特性行列に基づいて算出されるモード特性の値が上記測定データに基づいて算出されたモード特性の値に一致するように上記特性行列の各成分の値を変更していくことで、上記特性行列を同定する。
【００１０】
上記課題を解決するために、第２の発明では、解析対象物上に設定された複数の実測定点及び仮想測定点の座標データを要求する工程と、上記解析対象物上の少なくとも一点において加振する加振試験が行われた時に実測定点においてのみ測定を行い、加振により得られた各実測定点における測定データに基づいてモード特性を算出する工程と、上記座標データに基づいて上記解析対象物の振動を表す特性行列を構成する成分間の制約条件式を算出する工程と、上記測定データに基づいて、算出されたモード特性の値に基づいて上記制約条件式を用いて上記特性行列を同定する工程とを具備する、振動解析方法において、上記特性行列は実測定点及び仮想測定点の両測定点における自由度に等しい自由度で作成され、上記特性行列を同定する工程では、上記特性行列に基づいて算出されるモード特性のうち上記測定データに基づいて算出されたモード特性に対応するモード特性についてのみ、特性行列に基づいて算出されるモード特性の値が上記測定データに基づいて算出されたモード特性の値に一致するように上記特性行列の各成分の値を変更していくことで、上記特性行列を同定する。
【００１１】
第３の発明では、第２の発明において、上記仮想測定点はその仮想測定点に最も近接する測定点が実測定点となるように設けられる。
【００１２】
第４の発明では、第２又は第３の発明において、測定点間の構造的結合状態の入力を要求する工程をさらに具備し、上記制約条件式を算出する工程では、上記座標データに加えて測定点間の構造的結合状態に基づいて制約条件式を算出する。
【００１３】
第５の発明では、第２〜第４のいずれか一つの発明において、上記特性行列を同定する工程で同定される特性行列は質量行列であり、上記特性行列を同定する工程において同定された質量行列に基づいて剛体特性を算出する工程をさらに具備する。
【００１４】
上記課題を解決するために、第６の発明では、解析対象物上の測定点における位置、速度又は加速度の変位を測定する測定器と、上記解析対象物上の少なくとも一点において加振を行う加振器と、上記測定器及び加振器からそれぞれ測定データ及び加振データが入力されると共にこれらデータに基づいて演算処理を行う解析部とを具備する振動解析装置において、上記解析部は、上記加振器での加振試験により得られた各測定点における測定データに基づいてモード特性を算出し、各測定点の座標データに基づいて上記解析対象物の振動を表す特性行列を構成する成分間の制約条件式を算出し、上記測定データに基づいて算出されたモード特性の値に基づいて上記制約条件式を用いて上記特性行列を同定し、上記加振試験で上記測定点において実際に測定される測定自由度よりも上記特性行列の自由度の方が大きく、よって上記測定データに基づいて算出されるモード特性の値の数の方が上記特性行列に基づいて算出されるモード特性の値の数よりも少なく、上記特性行列を同定する際には、上記特性行列に基づいて算出されるモード特性のうち上記測定データに基づいて算出されたモード特性に対応するモード特性についてのみ、特性行列に基づいて算出されるモード特性の値が上記測定データに基づいて算出されたモード特性の値に一致するように上記特性行列の各成分の値を変更していくことで、上記特性行列を同定する。
【００１５】
上記課題を解決するために、第７の発明では、解析対象物上に設定された実測定点及び仮想測定点のうち実測定点における位置、速度又は加速度の変位を測定する測定器と、上記解析対象物上の少なくとも一点において加振を行う加振器と、上記測定器及び加振器からそれぞれ測定データ及び加振データが入力されると共にこれらデータに基づいて演算処理を行う解析部とを具備する振動解析装置において、上記解析部は、上記加振器による加振試験により得られた実測定点における測定データに基づいてモード特性を算出し、上記実測定点及び仮想測定点の座標データに基づいて特性行列を構成する成分間の制約条件式を算出し、上記測定データに基づいて算出されたモード特性の値に基づいて上記制約条件式を用いて上記特性行列を同定し、上記特性行列は実測定点及び仮想測定点の両測定点における自由度に等しい自由度で作成され、上記特性行列を同定する際には、上記特性行列に基づいて算出されるモード特性のうち上記測定データに基づいて算出されたモード特性に対応するモード特性についてのみ、特性行列に基づいて算出されるモード特性の値が上記測定データに基づいて算出されたモード特性の値に一致するように上記特性行列の各成分の値を変更していくことで、上記特性行列を同定する。
【００１６】
第８の発明では、第７の発明において、上記仮想測定点はその仮想測定点に最も近接する測定点が実測定点となるように設けられる。
【００１７】
第９の発明では、第７又は第８の発明において、上記制約条件式を算出する際には、上記座標データに加えて測定点間の構造的結合状態に基づいて制約条件式が算出される。
【００１８】
第１０の発明では、第７〜第９のいずれか一つの発明において、上記同定される特性行列は質量行列であり、上記解析部は、上記特性行列を同定する工程において同定された質量行列に基づいて剛体特性を算出する。
【発明の効果】
【００１９】
本発明によれば、実際に測定を行う測定点を減少させても、推定精度をほとんど低下させることなく剛体特性を推定することができる。
【発明を実施するための最良の形態】
【００２０】
以下、図面を参照して本発明の実施形態について詳細に説明する。
【００２１】
図１は、本実施形態の解析装置のブロック図である。図１において、解析装置１０は、解析対象物１の剛体特性等を解析的に算出可能な装置であり、測定部２０と解析部３０とを具備する。測定部２０は、先端に力変換器２１を備えた打撃ハンマ２２と、解析対象物１に取り付けられる複数の計測器（例えば、加速度計）２３ａ〜２３ｃと、力変換器２１及び各計測器２３ａ〜２３ｃに接続された増幅器２４、２５ａ〜２５ｃと、これら増幅器２４、２５ａ〜２５ｃに接続された外部インタフェース部２６とを具備する。
【００２２】
一方、解析部３０は、リムーバブルメディアへの書き込み等を行うリムーバブルメディア装置（ＲＭ装置）３１と、リムーバブルメディア装置３１を制御するリムーバブルメディア制御部（ＲＭ制御部）３２と、ユーザからの入力を受け付ける入力装置３３と、入力装置３３を制御する入力制御部３４と、ディスプレイ等への出力を行う出力装置３５と、出力装置３５を制御する出力制御部３６と、各種演算処理を行う中央処理装置３７と、記憶装置３８と、これら外部インタフェース部２６、リムーバブルメディア制御部３２、入力制御部３４、出力制御部３６、中央処理装置３７、記憶装置３８等を相互に接続するバス３９とを具備する。
【００２３】
打撃ハンマ２２は、解析対象物１を加振するのに用いられる。加振時に加えられた力の大きさは、打撃ハンマ２２に設けられた力変換器２１によって電気信号に変換される。この電気信号は、増幅器２４で増幅された後に外部インタフェース２６に出力される。なお、本実施形態では加振器として打撃ハンマ２２を用いて解析対象物１に非定常波を与えているが、加振器や加振方法はこれに限定されるものではない。解析対象物１を加振することができると共に解析対象物１に加えられた力の大きさが分かれば、打撃ハンマ２２の代わりに如何なる加振器を用いても良い。従って、加振器として、油圧式や圧電式等の加振器、及び音圧や磁界を用いた非接触加振器等を使用することができる。また、非定常波の代わりに、定常波や周期波等の各種の加振波形の加振を行ってもよい。さらに、本実施形態では、単点加振試験を行っているが、多点加振試験を行ってもよい。
【００２４】
加速度計２３ａ〜２３ｃはそれぞれ解析対象物１上の測定点に配置され、測定点において解析対象物１に発生した加速度の大きさを測定し、これを電気信号に変換して出力する。この電気信号は、増幅器２５ａ〜２５ｃで増幅された後に外部インタフェース２６に出力される。測定点は、注目したい次数までの固有モードを幾何学的に充分に表現することができる程度に適切な分布及び数となるようにユーザにより決定される。なお、本実施形態では、計測器として３軸加速度計２３ａ〜２３ｃを用いているが、渦電流やレーザ光線を利用した非接触変位計、レーザドップラ速度計、ひずみゲージ等、測定点における各種パラメータ（位置、速度、加速度等）の変位を計測することができれば加速度計２３ａ〜２３ｃの代わりに如何なる計測器を用いても良い。
【００２５】
外部インタフェース部２６は、力変換器２１及び加速度計２３から増幅器２４、２５を介して入力された電気信号をデジタル信号に変換し、このデジタル信号をバス３９に出力する。リムーバブルメディア装置３１は、リムーバブルメディア制御部３２によって制御される。リムーバブルメディア装置３１では、リムーバブルメディア２に記録されたプログラム（例えば、後述する解析を行うプログラム）等を解析部３０にインストールしたり、解析結果をリムーバブルメディア２に記録したりすることができる。なお、リムーバブルメディア２としては、ＦＤ、ＣＤ、ＤＶＤ、ＭＯ等、如何なる記憶媒体を用いても良い。
【００２６】
入力装置３３は、キーボード、マウス及びデジタイザ等の入力機器であり、入力制御部３４によって制御される。この入力装置３３により、ユーザは解析対象物の形状や各測定点の座標データ等、各種データを入力することができる。出力装置３５は、ディスプレイやプリンタ等の出力機器であり、出力制御部３６によって制御される。出力装置３５には、入力装置３３に入力された各種データ、プログラム実行のためのメニュー、解析結果等が出力される。
【００２７】
中央処理装置３７は、バス３９を介して、外部インタフェース部２６、リムーバブルメディア制御部３２、入力制御部３４、出力制御部３６及び記憶装置３８と接続され、これら制御部等に指示を与え、また、記憶装置１２９に格納され各プログラムに従いフーリエ変換や実験的特性行列同定法等の各種処理を行う。また、記憶装置３８は、ＲＡＭ（ランダムアクセスメモリ）やＲＯＭ（リードオンリメモリ）等から構成され、プログラムや、各プログラム実行中の一時的なデータ、解析結果等が格納される。
【００２８】
次に、このように構成された解析装置１０によって行われる解析処理について説明する。
【００２９】
まず、本実施形態の解析装置１０によって行われる解析処理についての基本的な概念について説明する。例えば、図２（ａ）に示したような解析対象物を多自由度系でモデル化する。特に、図２（ａ）に示した例では、計測点がｎ個あり、各計測点における自由度がｘ、ｙ、ｚ方向の３自由度であるため、解析対象物を自由度３ｎの多自由度系でモデル化している。この場合、質量行列、粘性減衰行列及び剛性行列をそれぞれ［Ｍ］、［Ｃ］、及び［Ｋ］とおくと、モデル化した解析対象物の運動方程式は下記式（１）で表される。なお、下記式（１）では粘性減衰を仮定しているが、必要に応じて構造減衰を仮定してもよいし、両者の減衰を同時に組み込んだ運動方程式を用いてもよい。
【数１】
【００３０】
ここで、式（１）中の｛ψ｝は測定点と定めた各点におけるｘ方向変位ｘi、ｙ方向変位ｙi、ｚ方向変位ｚiを縦に並べてできるベクトル、｛ｆ｝は測定点と定めた各点におけるｘ方向の加振力ｆｘi、ｙ方向の加振力ｆｙi、ｚ方向の加振力ｆｚiを縦に並べてできるベクトルを表している（下記式（２）参照）。
【数２】
【００３１】
なお、質量行列［Ｍ］、粘性減衰行列［Ｃ］及び剛性行列［Ｋ］（以下、これらをまとめて「特性行列」という）は３ｎ行３ｎ列の対称行列（すなわち、行数及び列数が解析対象物の測定自由度３ｎに等しい対称行列）であり、ｉは測定点の番号を示している。また、式（１）中の｛ψ’｝は変位を表すベクトルの一階微分、すなわち測定点の速度を表すベクトルであり、｛ψ’’｝は変位を表すベクトルの二階微分、すなわち測定点の加速度を表すベクトルである。本明細書では、行列を［］で、ベクトルを｛｝で表す。
【００３２】
上記運動方程式（１）において、質量行列［Ｍ］等の特性行列を求めると、解析対象物１の振動の様子を把握することができると共に、解析対象物１の剛体特性（例えば、解析対象物１の質量、質量中心位置、三つの主慣性モーメント及びこれらに対応する三つの主軸の向き等）を求めることができる。そこで、本実施形態では、打撃ハンマ２２によって解析対処物１が加振されたときの各測定点での測定データに基づいて、解析対象物１の特性行列を求めることとしている。
【００３３】
加振時の測定データに基づく解析対象物１の特性行列の同定及び剛体特性の算出は、図３に示したように大別して四つの処理工程によって行われる。
【００３４】
まず、ステップＳ１０に示したように、単点加振試験時の測定データに基づいて、モード特性の同定が行われる。すなわち、打撃ハンマ２２によって解析対象物１の一点が打撃（加振）されると共に、この加振力が力変換器２１によって測定され且つこの加振によって発生した振動が各加速度計２３によって測定される。その後、このようにして力変換器２１及び加速度計２３によって測定された測定データから各測定点について周波数応答関数及びコヒーレンス関数が算出される。次いで、算出された周波数応答関数を用いて、一般的な方法で、固有振動数、固有モードベクトル及びモード減衰比（以下、これらをまとめて「モード特性」という）を同定する。
【００３５】
次いで、ステップＳ２０では、特性行列成分間の物理的制約条件式の算出が行われる。すなわち、特性行列の特定の成分の値や特定の成分間の制約条件は、各種力学原理、各測定点の座標、測定点間の構造的結合状態等に基づいて、加振試験の測定結果とは無関係に求められる。例えば、測定点間の構造的結合状態に基づいて特性行列中の零成分の配置を或る程度特定することができる。また、各種力学原理や測定点の座標等に基づいて、特性行列中の非零成分間の関係（制約条件）を或る程度特定することができる。
【００３６】
このように、特性行列成分間の制約条件式を算出することにより、上記ステップＳ１０において求められたモード特性に基づいて同定しなければならない特性行列中の成分の数が減少する。例えば質量行列［Ｍ］が３ｎ行３ｎ列の行列である場合、この質量行列［Ｍ］中には３ｎ（３ｎ＋１）／２個の成分が含まれている（実際の成分の数は３ｎ×３ｎ個であるが、行列の対称性を考慮すると３ｎ（３ｎ＋１）／２個であるということができる）。ここで、この質量行列［Ｍ］の成分間の制約条件等が求められていない場合には、上記ステップＳ１０において求められたモード特性に基づいて同定しなければならない質量行列［Ｍ］中の成分の数は３ｎ（３ｎ＋１）／２ということになり、現実的には質量行列［Ｍ］の全ての成分を適切に同定することは困難である。
【００３７】
しかしながら、本ステップＳ２０では質量行列［Ｍ］の成分間の制約条件式が求められ、よって質量行列［Ｍ］中の零成分の配置や、質量行列［Ｍ］中の特定の成分と別の特定の成分との関係式が求められる。従って、質量行列［Ｍ］成分間の制約条件式を算出することにより、上記ステップＳ１０において求められたモード特性に基づいて同定しなければならない質量行列［Ｍ］中の成分の数がかなり減少する。具体的には、質量行列［Ｍ］中の全ての成分の数３ｎ（３ｎ＋１）／２個から、質量行列［Ｍ］中の零成分の数及び作成された制約条件式の数を減算した数の成分（以下、「独立未知成分」という）を上記ステップＳ１０において求められたモード特性に基づいて同定することになる。
【００３８】
次いで、ステップＳ３０では、数学的最適化法により質量行列［Ｍ］及び剛性行列［Ｋ］の独立未知成分の同定が行われる。すなわち、質量行列［Ｍ］及び剛性行列［Ｋ］と固有角振動数Ωt及び固有モードベクトル｛ψ｝tとの関係は下記式（３）のように表せる（なお、ｔは次数を表す）。
【数３】
【００３９】
式（３）から分かるように、質量行列［Ｍ］及び剛性行列［Ｋ］の各成分を変数とみなしてその値を変化させると、これに伴って固有角振動数Ωt及び固有モードベクトル｛ψ｝tの値も変化することになる。従って、本ステップＳ３０では、上記式（３）によって求められる固有角振動数Ωt及び固有モードベクトル｛ψ｝tが上記ステップＳ１０で求められた固有角振動数及び固有モードベクトルに一致するように、ステップＳ２０で求められた質量行列［Ｍ］及び剛性行列［Ｋ］の成分間の制約条件式を満足させながらこれら質量行列［Ｍ］及び剛性行列［Ｋ］の独立未知成分の値が徐々に変化せしめられる。固有角振動数同士及び固有モードベクトル同士がほぼ一致したときの質量行列［Ｍ］及び剛性行列［Ｋ］が、図２（ａ）に示したように多自由度系でモデル化された解析対象物の質量行列［Ｍ］及び剛性行列［Ｋ］として同定される。
【００４０】
以下では、このように固有角振動数同士及び固有モードベクトル同士を一致させることを固有角振動数及び固有モードベクトルの「一致化」として説明し、固有角振動数同士及び固有モードベクトル同士を一致させることで質量行列［Ｍ］及び剛性行列［Ｋ］の各成分の最適な値を算出・同定することを質量行列［Ｍ］及び剛性行列［Ｋ］の「最適化」として説明する。
【００４１】
なお、このように質量行列［Ｍ］及び剛性行列［Ｋ］の最適化を行うための数学的手法としては、「ニュートン法類に属する最急降下法」等、多くの手法を適用することが可能であり、最適化が効率よく且つ安定して実行できれば如何なる手法を用いることも可能である。
【００４２】
次いで、ステップＳ４０では、ステップＳ３０において同定された質量行列［Ｍ］に基づいて、解析対象物１の剛体特性、すなわち解析対象物１の質量、質量中心位置、三つの主慣性モーメント及びこれらに対応する三つの主軸の向きが求められる。
【００４３】
ところで、上述した特性行列及び剛体特性の同定法では、適当な数ｎの測定点を解析対象物１に設定して、加振試験時においてこれら各測定点でｘ、ｙ、ｚ方向の振動について周波数応答関数を計測し、それに対応した自由度の特性行列を同定することで、剛体特性を推定している。解析対象物１が複雑な構造物であるような場合には、その形状の複雑さや大きさに応じて測定点の数を或る程度多く設定する必要がある。
【００４４】
ところが、上述した特性行列及び剛体特性の同定法においては、解析対象物１上への測定点の設定から周波数応答関数の算出までの間に、測定点の設定、加速度計のキャリブレーション、計測システムのセッティング等の作業が行われ、よって測定点の設定から周波数応答関数の算出までが最も手間と時間のかかる作業である。従って、測定点数が増大すると、それに伴って周波数応答関数の算出までにかかる手間と時間が大きく増大し、その結果、特性行列及び剛体特性の同定にかかる時間が増大してしまう。また、周波数応答関数は中央処理装置３７において行われる同定処理によって求められるため、測定点数が増大すると中央処理装置３７における作業負荷が増大してしまい、周波数応答関数の収束速度が大きく低下してしまう。したがって、特性行列及び剛体特性の同定を迅速に行うという観点からは、実際に測定を行う測定点の数が少ないことが好ましい。
【００４５】
そこで、本発明の実施形態では、測定点数を減少させずに、すなわち特性行列の自由度を小さくすることなく、実際に測定を行う測定点の数を減少させることとしている。すなわち、解析対象物１上に従来と同様なほどの十分な数の測定点を設定すると共に、そのうちの一部の測定点については実際に測定を行って周波数応答関数を算出し（以下、このような測定点を「実測定点」という）、残りの測定点については実際には測定を行わないようにしている（以下、このような測定点を「仮想測定点」という）。換言すると、本実施形態では、実際に測定が行われる測定点の自由度よりも大きな自由度の特性行列を同定することとしている。
【００４６】
本実施形態では、このように実測定点と仮想測定点とを設定した上で、上記ステップＳ３０において、数学的最適化法により質量行列［Ｍ］及び剛性行列［Ｋ］の独立未知成分の同定を行うにあたり、固有モードベクトルのうち実測定点に対応する固有モードベクトルの成分のみを一致化させるようにしている。
【００４７】
ここで、本願発明者は、このように仮想測定点という概念を導入した場合であっても、測定点間の構造的結合状態等に基づいて適切に制約条件式が求められている場合には、固有モードベクトルのうち実測定点に対応する成分のみを一致化させることができれば、仮想測定点に対応する成分も従属的に（または自動的に）一致化させることができ、その結果、適切な質量行列［Ｍ］及び剛性行列［Ｋ］を同定することができることを発見した。換言すると、本願発明者は、固有モードベクトルの一致化をするにあたって、必ずしも全ての固有モードベクトルの成分について一致化させなくても、十分な精度で特性行列の最適化を行うことができることを発見したのである。
【００４８】
以下では、本発明の実施形態に係る解析装置１０によって行われる解析処理の手順について詳細に説明する。なお、本実施形態における解析処理も、基本的に図３に示した処理手順に従って行われる。従って、以下では、図４〜図７を参照して、上記解析処理の処理手順に対応する形で本実施形態の解析処理の処理手順について説明する。
【００４９】
［モード特性の同定］
まず、ステップＳ１０に示したように、単点加振試験による測定結果に基づいてモード特性の同定が行われる。図４は、本発明のモード特性の同定処理の処理手順を示すフローチャートである。モード特性の同定は、加振試験を行って周波数応答関数等を算出する振動計測処理と、算出された周波数応答関数等に基づいてモード特性を同定する特性同定処理とに大別される。
【００５０】
図４に示した処理手順の実行の前に、ユーザによる計測システムのセッティングが行われる。計測システムのセッティングにあたっては、解析対象物１がユーザによりゴム紐やタイヤチューブ等の弾性体により支持されるように配置される。これにより、解析対象物１の周辺自由状態が近似的に実現される。なお、解析対象物１の支持は、擬似的に周辺自由状態が実現されれば、如何なる態様で行われてもよい。
【００５１】
図４に示したようなモード特性の同定処理においては、まず、ステップＳ１１において、中央処理装置３７は、ユーザに解析対象物１の形状・寸法等の初期データを入力するように出力画面３５を介して要求する。このような初期データには、実測定点及び仮想測定点を含む全ての測定点の座標データが含まれる。
【００５２】
ユーザは実測定点及び仮想測定点を解析対象物１上の任意の位置に設定する。このとき、仮想測定点は、できるだけ仮想測定点同士が隣り合うことの無いように、すなわちその仮想測定点に最も近接する測定点が実測定点となるように設定される。換言すると、仮想測定点は、少なくとも一つの或る方向において実測定点の間に配置されるように設定される。また、仮想測定点は、仮想測定点数が全測定点数の１／２以上にならないように（好ましくは、１／３以上にならないように）設定される。
【００５３】
図２（ｂ）は実測定点及び仮想測定点の設定の例を示している。図２（ｂ）中に丸で囲まれた測定点は実際に計測が行われる実測定点を示しており、丸で囲まれていない測定点は実際には計測が行われない仮想測定点を示している。すなわち、図示した例では、３、６、７、１１番の測定点においては実際に測定が行われない。
【００５４】
ユーザは、このような条件のもと各測定点の設定位置を決定すると共に、中央処理装置３７の上記要求に応じて全ての測定点の座標データを入力する。また、実測定点の設定位置には加速度計２３が取り付けられる。
【００５５】
次に、ユーザは、打撃ハンマ２２によって加速度計２３が取り付けられた位置の中のいずれかの１つの位置（単点）で解析対象物１に打撃（加振）を加える。
ユーザが打撃ハンマ２２によって解析対象物１に加えた打撃は力変換器２１によって測定され、その測定データは増幅器２４を介して外部インタフェース部２６に入力される。また、この打撃によって解析対象物１に発生した振動も加速度計２３によって測定され、その測定データは増幅器２５を介して外部インタフェース部２６に入力される。入力されたこれら測定データは中央処理装置３７に取り込まれる。このようにして単点加振試験が行われる。
【００５６】
次に、ステップＳ１３において、測定データに基づいて周波数応答関数等が算出される。すなわち、中央処理装置３７は、これら力変換器２１及び加速度計２３によって測定された測定データをフーリエ変換（例えば、高速フーリエ変換）して周波数スペクトルになおす。その後、中央処理装置３７は、解析対象物１への入力に関する周波数スペクトル（力変換器２１の出力に基づく周波数スペクトル）及び解析対象物１からの出力に関する周波数スペクトル（加速度計２３の出力に基づく周波数スペクトル）からパワースペクトルとクロススペクトルとを算出し、これらパワースペクトル及びクロススペクトルから周波数応答関数及びコヒーレンス関数を算出する。これら周波数応答関数及びコヒーレンス関数は各実測定点毎に算出される。
【００５７】
ここで、コヒーレンス関数（関連度関数）は、入力と出力と内積であり、入出力の相関を示すものである。その値が０である場合には、入力と出力とは無相関であり、１に近づくにつれて相関性が増す。このコヒーレンス関数は、後述する周波数応答関数に基づくモード特性の同定にも用いられる（ステップＳ１５で使用）。
【００５８】
なお、より信頼性の高いデータを得る観点から、この単点加振試験を複数回繰り返し、個々のパワースペクトルおよびクロススペクトルを平均し、その平均値を用いるようにすることが好ましい。
【００５９】
上記ステップＳ１１〜Ｓ１３の処理が振動計測処理であり、単点加振試験により得られた測定データに基づいて各実計測点における周波数応答関数及びコヒーレント関数が算出される。
【００６０】
このように、周波数応答関数が算出された後で、算出された周波数応答関数等に基づいてモード特性を同定する特性同定処理が行われる。
【００６１】
特性同定処理では、まず、ステップＳ１４に示したように、中央処理装置３７が、ユーザにモード特性の同定を行う周波数帯域（以下、「同定周波数帯域」という）を入力するように出力画面３５を介して要求する。ユーザは入力装置３３から同定周波数帯域を入力する。同定周波数帯域は、その帯域幅が広いほどその帯域内に現れる共振峰の数が多くなる。このため、同定周波数帯域を広げると、多くの次数のモード特性を同定することができる。
【００６２】
次いで、ステップＳ１５において、ステップＳ１３で算出された周波数応答関数及びコヒーレント関数に基づいて、所定の次数までのモード特性が同定される。この所定の次数は、上記ステップＳ１４で入力された同定周波数帯域応じて変化する。なお、モード特性の同定法としては一般に知られている各種のモード同定法を利用することができ、このようなモード同定法としては、例えば、多点参照法（ポリリファレンス法、poly-referencemethod）、偏分反復法、プロリー法（Prony's method）、サークルフィット法（circle fit method）等が挙げられる。
【００６３】
以下では、周波数応答関数に基づいてモード特性を同定する際の基本的な原理について簡単に説明する。減衰を比例粘性減衰とした振動系の限定された周波数帯域（すなわち、上記同定周波数帯域）に関する周波数応答関数は近似的に下記式（４）で表せる。
【数４】
ここで、式（４）中のＹは慣性項定数パラメータ、Ｚは剰余項定数パラメータであり、また同定周波数帯域内に存在する共振峰の数をａとする。さらに、そのａ個の共振峰に対応する固有角振動数をΩt、モード減衰比をζt、加振点ｑ、応答点（測定点）ｐの固有モード成分をそれぞれψqt、ψpt（ｔは次数）とする。
【００６４】
このように式（４）で表した周波数応答関数Ｈpqは、固有角振動数Ωt、モード減衰比ζt、固有モード成分ψpt、ψqtを適切に選択すれば、ステップＳ１３において振動試験に基づく周波数応答関数にほぼ一致する。そこで、本モード特性同定では、式（４）によって算出される周波数応答関数ＨpqがステップＳ１３において算出された周波数応答関数にほぼ一致するような固有角振動数Ωt、モード減衰比ζt、固有モード成分ψpt、ψqtを最適化法等によって算出し、これらモード特性を同定するようにしている。また、このモード特性の同定の過程で上記コヒーレンス関数が用いられる。
【００６５】
上記ステップＳ１４及びＳ１５の処理が特性同定処理であり、周波数応答関数に基づいてモード特性が同定される。
【００６６】
［特性行列成分間の物理的制約条件式の同定］
次に、図３にステップＳ２０に示したように、特性行列成分間の物理的制約条件式の算出が行われる。図５は、本発明の特性行列成分間の物理的制約条件式の算出処理の処理手順を示すフローチャートである。
【００６７】
特性行列成分間の物理的制約条件式の算出にあたって、中央処理装置３７は、ステップＳ２１において、測定点間の構造的結合状態（コネクティビティ）を入力するように出力画面３５を介して要求する。ここで、測定点間の構造的結合状態とは、或る測定点と別の測定点とが構造的に結びついているか否か、すなわち或る測定点と他の測定点とが構造的に隣り合う関係にあるか否かの状態を意味する。
【００６８】
例えば図２（ａ）に示した構造について考慮すると、測定点間の構造的結合状態は図２（ｃ）のように表される。図中、互いに構造的に結合しているとされる測定点同士は直線によって互いに結ばれている。図２（ａ）に示した例では、例えば測定点１は測定点２、３、４と互いに構造的に結合されており、測定点３は測定点１、２、４及び５と互いに構造的に結合されている。なお、本実施形態では、測定点間の構造的結合状態はユーザの主観により測定点の配置や測定点間の距離等に基づいて定められて、ユーザにより入力装置３３を用いて入力される。しかしながら、例えば、ユーザに解析対象物１の形状及び測定点の配置を入力させて、この入力されたデータに基づいて中央処理装置３７により測定点間の構造的結合状態が算出されるようにしてもよい。
【００６９】
次いで、ステップＳ２２では、ステップ２１において入力された測定点間の構造的結合状態に基づいて質量行列［Ｍ］及び剛性行列［Ｋ］の零成分の配置が特定される。以下、零成分の配置について説明する。
【００７０】
ここで、図６（ａ）に示したような単純な均一剛体要素と回転バネで直列に結合されたモデルを考える。各剛体要素は、均一形状で長さがそれぞれｌi、質量中心が中央にありその質量がｍi、回転バネ定数がｋiであるとして、曲げ振動表現の自由度（上記実施形態における測定点数に相当）を点１〜５番の並進変位とする。この場合の質量行列［Ｍ］及び剛性行列［Ｋ］はそれぞれ下記式（５）、（６）のように表される。
【数５】
【数６】
【００７１】
ここで、式（５）、（６）中の＊は零成分でないこと、すなわち非零成分であることを意味するものである。このように、図６に示したようなモデルでは、質量行列［Ｍ］については３重対角行列、剛性行列［Ｋ］については５重対角行列として表される。図６（ｂ）、（ｃ）を参照してこの理由について説明する。なお、図６では、剛体要素間の寸法を大きくして回転弾性対偶を描いているが、剛体要素間の寸法はないものとして考える。すなわち、剛体要素Ｉの右端と剛体要素ＩＩの左端とは同一位置であるとして考える。
【００７２】
図６（ｂ）を参照して、剛体要素Ｉの質量中心（剛体要素Ｉの図心）について、自由体図によって運動方程式を作成すると、回転運動については式（７）のように、並進運動については式（８）のように表される。これら式（７）、（８）を行列式でまとめて表示すると、式（９）のようになる。
【数７】
【数８】
【数９】
【００７３】
剛体要素Ｉの両端の並進変位自由度と質量中心での変位自由度との変換則は下記式（１０）のように表せる。そこで、式（１０）を式（９）に代入して、両辺左側から［Ｔ1］Tを乗じると、下記式（１１）のようになる。なお、［Ｔ1］TはＴ1の転置行列を表す。
【数１０】
【数１１】
ここで、剛体要素Ｉが均一形状であって均一材料密度分布であるとすると、下記式（１２）のように表せる。なお、式（１２）においてｍ1は剛体要素１の質量である。
【数１２】
【００７４】
同様に、図６（ｃ）を参照して、剛体要素ＩＩの質量中心（要素の図心）について自由体図によって上記剛体要素Ｉの場合と同様に運動方程式を作成して整理すると、下記式（１３）が導かれる。
【数１３】
【００７５】
この様に順々に剛体要素III、剛体要素IV、…と運動方程式を求めて、全ての運動方程式を行列表現でまとめて、全系の運動方程式を作成すると、下記式（１４）に示したようにその質量行列［Ｍ］は３重対角行列、剛性行列［Ｋ］は５重対角行列になる。このように、質量行列［Ｍ］及び剛性行列｛Ｋ］の非零成分の配置は、測定点間の構造的結合状態によって決定することができる。
【数１４】
【００７６】
上述の力学原理に基づいて、図２（ｃ）に示したような測定点間の構造的結合状態となっている解析対象物のモデルを考えると、例えば剛性行列［Ｋ］の零成分の配置は下記式（１５）のようになる。
【数１５】
【００７７】
上記式（１５）の剛性行列［ｋ］中の［Ｎ］は３行３列の非零成分となる可能性のある小行列部分を、［０］は３行３列の零成分となる小行列部分を表している。ここで、［Ｎ］、［０］が３行３列の小行列となっているのは各測定点における測定自由度がｘ、ｙ、ｚ方向の３自由度であるためであり、例えば各測定点における測定自由度がｘ、ｙ、ｚ方向に加えてｘ軸、ｙ軸、ｚ軸回りの回転方向を加えた６自由度である場合には［Ｎ］、［０］は６行６列の小行列となる。
【００７８】
図２（ｃ）に示した例では、測定点１は測定点２、３に物理的に結合しており、これら測定点２又は３は測定点４、５に物理的に結合している。従って、式（１７）に示したように、剛性行列［Ｋ］では、測定点１に対応する１列目には１〜５行目に非零小行列成分が配置されることになる。このように、ユーザがデザインする構造的結合モデルに従って特性行列の零成分配置が力学原理に基づいてコンピュータにより自動的に決定される。
【００７９】
次に、ステップＳ２３において、３次元構造物の質量行列［Ｍ］については、如何なるｎ自由度のモデルについての質量行列［Ｍ］も下記式（１６）左辺の演算によって、構造物中に任意に設定された剛体運動表現点としての１点の６自由度変位に関する６行６列の剛体質量行列Ｍrigid］に変換できるという力学原理に基づいて、質量行列［Ｍ］の非零成分間の連立１次方程式が物理的制約条件式として作成される。ここで、式（１６）中の行列［Φ］は下記式（１７）のように表される３ｎ行６列の行列（以下、「剛体変位モード行列」という）であり、剛体質量行列[Ｍrigid］は下記式（１８）のように表される６行６列の行列である。
【数１６】
【数１７】
【数１８】
【００８０】
なお、剛体質量行列中のＩxx、Ｉyy、Ｉzzは剛体運動表現点を通り各座標軸と並行な軸回りの慣性モーメント、Ｉyx、Ｉzx、Ｉyzはそれら３軸相互間についての慣性乗積である。Ａ、Ｂ、Ｃは、剛体運動表現点と重心とのズレに起因する値であり、解析対象物１の質量をｍ、重心座標を（ｘG、ｙG、ｚG）、剛体運動表現点座標を（ｘR、ｙR、ｚR）とすると、Ａ＝ｍ（ｘG−ｘR）、Ｂ＝ｍ（ｙG−ｙR）、Ｃ＝ｍ（ｚG−ｚR）である。［Ｍrigid］は、上記式（１８）に示すように対称行列である。
【００８１】
以下、上記剛体変位モード行列［Φ］、剛体質量行列［Ｍrigid］及び上記式（１８）の導出方法について簡単に説明する。
上記剛体変位モード行列［Φ]は、６本の互いに独立な剛体運動ベクトルで構成される。具体的に、入力データの一部である測定点座標値を利用して、剛体運動表現点がｘ軸方向、ｙ軸方向、ｚ軸方向に単位並進変位するときの構造物上の測定点の変位を表すベクトルをそれぞれ第１列、第２列、第３列に当てはめ、剛体運動表現点がその点を通ってｘ軸と平行な軸回り、ｙ軸と平行な軸回り及びｚ軸と平行な軸回りに微少角変位する剛体角変位運動に関する構造物上の測定点の並進変位を線形近似表現するベクトルを第４列、第５列、第６列に当てはめることで剛体変位モード行列［Φ］は容易に作成することができる。
【００８２】
例えば、図７を参考に、剛体中に設定した剛体運動表現点座標をＲ（ｘR，ｙR，ｚR）、測定点ｉ番のｘ、ｙ、ｚ座標を（ｘi，ｙi，ｚi）として、剛体が単にｘ軸方向へδRxだけ並進変位するとそれに伴って測定点ｉ番の変位は、下記式（１９）のように、δRxに定数で構成される変換行列（この場合は列ベクトル）を乗じることで表現することができる。
【数１９】
【００８３】
剛体が剛体運動表現点で単にｘ軸に平行な軸回りに角度θxだけ右手に回転することによる測定点ｉ番の変位は幾何学から下記式（２１）のように代数的に表現することができる。
【数２０】
【００８４】
ここで、今対象としている機械振動の並進振幅や角変位振動は小さい値であることを考慮して、ｓｉｎθｘ≒θｘ、ｃｏｓθｘ≒１の近似を十分適用することができるため、上記式（２０）は下記式（２１）と変形することができる。式（１９）と同様に座標データに基づく定数で構成される変換行列を運動表現点の角変位に乗じることで測定点ｉの変位が表現できる。
【数２１】
【００８５】
このようにして６自由度の全ての剛体運動を考えると、剛体変位モード行列［Φ］の成分構成は測定点ｉのｘ、ｙ、ｚ方向変位成分についてのみ明記した形式として、下記式（２２）のように作成することができる。
【数２２】
【００８６】
従って、各測定点の座標データを入力することにより、上記剛体変位モード行列［Φ］を求めることができる。なお、フーリエ変換、モード特性同定法等の実験的同定法では、すべてデータ処理法（信号処理法）に過ぎず、幾何学的データ（または構造的データ）としての測定点座標を本質的に必要する「実験的同定法」はこれまで皆無であったと考えることができ、本手法はこの点についても大きな特徴を有する。
【００８７】
次に、剛体質量行列［Ｍrigid］について説明する。上記式（１６）の右辺の剛体質量行列［Ｍrigid］は、如何なる対象構造物に関しても上記式（１８）のような独特な成分配置関係となることが知られている。例えば、剛体質量行列［Ｍrigid］には、１行１列成分と２行２列成分及び３行３列成分は常に同一の値でありその値は解析対象物の質量の値であること、２行１列成分と１行２列成分は常に零となること、５行１列成分と４行２列成分とは絶対値は同一で符号が互いに逆の値となること、等の特徴的な成分配置関係がある。
【００８８】
このように表される剛体変位モード行列［Φ］、剛体質量行列［Ｍrigid］及び質量行列［Ｍ］の間には、上記式（１６）のような関係がある。ここで、式（１８）の非零成分（ｍ、Ａ、Ｉxx等）の値そのものは同定前には未知であるが、これら非零成分間の関係を利用して、式（１６）からその左辺の質量行列［Ｍ］の成分間に関する１１本の等式制約条件式を自動生成することができる。
【００８９】
なお、解析対象物１の質量（重さ）等の剛体特性値が既知として与えられていれば、これを式（１８）に代入して質量行列成分間の制約条件式を生成するようにしてもよい。このように既知の剛体特性値が得られれば、この既知の剛体特性値については、他の成分で表す必要が無くなるという利点がある。
【００９０】
このようにして生成された連立方程式は、互いに独立な方程式よりも未知数である非零成分の数の方が必ず多くなる。従って、このように生成された連立方程式のみからは質量行列［Ｍ］の非零成分を全て特定することはできない。しかしながら、方程式の数だけの未知数（非零成分）を他の未知数（非零成分）の従属変数とすることができる。以下の説明では、本ステップＳ２３等において求められた制約条件式により特性行列の他の成分に従属する形で表現される成分を従属成分と、制約条件式によっても他の成分に従属する形では表現されない成分を独立未成分として説明する。
【００９１】
このような等式制約条件式を上記手法により自動生成すると、例えば下記式（２３）のような等式が作成される。
【数２３】
ここで、ｍi,jは質量行列のｉ行ｊ列の成分を表し、左辺の成分が従属成分、右辺の各項の成分が独立未知成分である。なお、どの成分が独立未知成分に割り当てられ、どの成分が従属成分に割り当てられるかはプログラミングに依存する。例えば、式（２３）の最上段の式の場合、ｍ1,1を独立未知成分とすることもでき、この場合ｍ2,1が従属成分となる。
【００９２】
次に、ステップＳ２４において、３次元構造物の剛性行列［Ｋ］と減衰行列［Ｃ］について、剛体運動状態ではその構造物のどの位置にもひずみが生じず、そのひずみに基づく応力も全ての部分において零であるという基本的な力学原理に基づいて、剛性行列［Ｋ］及び減衰行列［Ｃ］の非零成分間に連立１次方程式が物理的制約条件式として作成される。上記力学原理からは、下記式（２４）が導き出される。なお、式（２４）の右辺の［０］はｎ行６列の零行列である。
【数２４】
【００９３】
この式（２４）により、各特性行列中の非零成分に関する連立１次方程式が作成される。これにより、ステップＳ２３と同様に、剛性行列［Ｋ］及び粘性行列［Ｃ］についても、成分間に関する複数の等式制約条件式が生成され、生成された方程式の数だけ独立未知成分を減少させることができる。
【００９４】
このように、ステップＳ２１〜Ｓ２４では、物理原則等に基づいて各特性行列毎にその各成分間の関係を示す制約条件式が求められる。これにより、各特性行列中に含まれる独立未知成分の数を減少させると共に力学（物理学）的成分間関連特性を付与することができ、上記ステップＳ３０における最適化処理における最適化精度を高めることができる。
【００９５】
［特性行列の同定］
次に、図３にステップＳ３０で示したように、質量行列［Ｍ］、剛性行列［Ｋ］及び粘性行列［Ｃ］中の独立未知成分の同定、従ってこれら特性行列の同定が行われる。上述したように、質量行列［Ｍ］及び剛性行列［Ｋ］の最適化を行うための数学的手法としては多くの手法を採用可能であり、以下では図８及び図９を参照してその一つの例について説明する。図８及び図９は、本発明の特性行列の同定処理の処理手順を示すフローチャートである。
【００９６】
本実施形態の特性行列の同定処理は、三つの処理に大別される。一つ目の処理は、独立未知成分として適当な値が代入された質量行列［Ｍ］が正定値行列になるように、また、独立未知成分として適当な値が代入された剛性行列［Ｋ］が非負定値行列となるように、独立未知成分の値を修正する正定値・非負定値化処理である。二つ目の処理は、上記式（３）によって求められる固有角振動数Ωt及び固有モードベクトル｛ψ｝tが上記ステップＳ１０で求められた固有角振動数及び固有モードベクトル（以下、これらをそれぞれ「目標固有角振動数Ωtt」、「目標固有モードベクトル｛ψ｝ttという）に一致するように、質量行列［Ｍ］及び剛性行列［Ｋ］の独立未知成分の値が徐々に変化させる最適化処理である。そして、三つ目の処理は、正規化されている質量行列［Ｍ］及び剛性行列［Ｋ］にゲイン平行移動係数αを乗算して、最終的に質量行列［Ｍ］及び剛性行列［Ｋ］を同定するゲイン調整処理である。図８及び図９では、これら三つの大きな処理毎に破線で囲まれている。
【００９７】
特性行列の同定処理では、まずステップＳ３１において、質量行列［Ｍ］及び剛性行列［Ｋ］の各成分のうちステップＳ２０（具体的にはステップＳ２４）で物理的制約条件式を作成した結果独立未知成分となった成分に初期値として乱数を代入する。また、代入された乱数に基づいてステップＳ２０において作成された制約条件式を用いて従属成分を算出し、初期値としての質量行列［Ｍ］と剛性行列［Ｋ］を数値行列として作成する。
【００９８】
しかしながら、この初期値としての質量行列［Ｍ］と剛性行列［Ｋ］は乱数から作成されているため、物理学的に最低限満足しなければならない条件をも満足しない場合がある。すなわち、質量や回転慣性は必ず正の値であることから質量行列は正定値行列となる。また、同様に剛性行列も基本的に正の値となるが、剛体運動ではひずみが零であるので固有値に零のものが最大６個まで存在する（剛体変位自由度だけ固有値零が存在する）ことから剛性行列は非負定値行列となる。しかしながら、上記初期値としての質量行列［Ｍ］と剛性行列［Ｋ］はこのような条件を満足しない場合がある。
【００９９】
そこで、本実施形態では、ステップＳ３２〜３４において、質量行列については正定値行列化を行い、剛性行列については非負定値行列化を行う。換言すると、質量行列［Ｍ］に関しては全ての３ｎ個の固有値が正値となるように質量行列［Ｍ］の独立未知成分を修正し、剛性行列［Ｋ］の全ての３ｎ個の固有値が負値とならないように剛性行列［Ｋ］の独立未知成分を修正する。ここでは、代表として質量行列の正定値行列化手法について説明する。
【０１００】
まず、質量行列［Ｍ］をその成分の絶対値最大の値に対して正規化をした上で、すなわち質量行列［Ｍ］の成分の最大絶対値が１になるようにした上で、標準的固有値問題である下記式（２５）を解析して、３ｎ個の全ての固有値λと固有ベクトル｛ｘ｝を求めて、負値の固有値を選び出す。
【数２５】
このようにして算出された固有値λの中に負値のものがなければ、質量行列［Ｍ］はこの時点の各成分の値で正定値行列である。一方、固有値λに負値のものがある場合には、これら負値の固有値λを正値へ変化させる必要がある。このために、負値の固有値λを正値へ変化させるための１次感度解析を行う。この感度解析としては、例えば、最急降下法、共役勾配法といったニュートン法に代表されるニュートン法類の手法等、一般に知られている各種の感度解析の手法を利用することができる。以下では、その一例について示す。質量行列［Ｍ］内の独立未知成分ｐiに関する固有値λについての１次感度ｓt,iは上記式（２５）を独立未知成分ｐiで微分した下記式（２６）で算出される。
【数２６】
【０１０１】
ここで、式（２６）中のλtと｛ｘ｝tはそれぞれｔ次の固有値と固有ベクトルを示している。また、ｐiは質量行列［Ｍ］中にｑ個ある独立未知成分のうちの一つであるとする（ｉ＝１〜ｑ）。上記式（２６）で算出される１次感度ｓt,iを用いて、ｑ個の独立未知成分を微少量Δｐi（ｉ＝１〜ｑ）だけ増加させた場合の固有値λtの増分Δλtは近似的に下記式（２７）のように見積もることができる。
【数２７】
【０１０２】
そこで、負値の固有値λ全てについてこの式を生成する。このとき固有値の増分Δλtには、例えば、その負値の固有値の絶対値よりも大きい正値が代入される。そして、このようにして生成された式を満たすような独立未知成分の微少増分Δｐiを算出し、このようにして算出された独立未知成分の微少増分Δｐiを用いて独立未知成分の更新を行う。このような操作を繰り返すことにより、質量行列［Ｍ］の正定値化を行うことができる。なお、剛性行列［Ｋ］の非負定値化についても、上記質量行列［Ｍ］の正定値化と同様に行われる。
【０１０３】
このような正定値化の概念に基づいて、ステップＳ３２では、上記式（２５）に基づいて算出された全ての固有値λtのうち負値のものが選定される。次いで、ステップＳ３３では、上記式（２６）を用いてステップＳ３２で選定された負値の固有値λ毎に各独立未知成分についての固有値の１次感度ｓt,iが算出される。次いで、ステップＳ３４では、ステップＳ３３において算出された１次感度を用いて上記式（２７）により独立未知成分の微少増分Δｐiが算出され、算出された微少増分Δｐiを用いて下記式（２８）により独立未知成分の更新が行われる。
【数２８】
【０１０４】
次いで、ステップＳ３５では、ステップＳ３４において独立未知成分の更新が行われた結果、質量行列［Ｍ］の正定値化及び剛性行列［Ｋ］の非負定値化が実現されたか否かが判定され、質量行列［Ｍ］の正定値化及び剛性行列［Ｋ］の非負定値化が実現されていないと判定された場合にはステップＳ３２〜ステップ３４が繰り返される。一方、ステップＳ３５において、質量行列［Ｍ］の正定値化及び剛性行列［Ｋ］の非負定値化が実現されたと判定された場合には、ステップＳ３６へ、すなわち最適化処理へと進む。ステップＳ３２〜ステップＳ３５までの処理が上述した正定値・非負定値化処理に該当する。
【０１０５】
特性行列の同定処理では、次に最適化処理が行われる。最適化処理では、下記式（２９）の一般的固有値問題を解いて算出される固有角振動数（不減衰固有角振動数）Ωtと固有モードベクトル｛ψ｝tを、実験で得られた固有角振動数（目標固有角振動数Ωtt）及び固有モードベクトル（目標固有モードベクトル｛ψ｝tt）に一致化させる処理が行われる。このため、本実施形態では、これら固有角振動数同士及び固有モードベクトル同士を一致化させるべく、１次感度解析を行う。この感度解析においても、一般に知られている各種の感度解析の手法を利用することができる。
【数２９】
【０１０６】
固有角振動数Ωについての１次感度は、第ｔ次の固有角振動数Ωtと固有モードベクトル｛ψ｝tについて上記式（２９）を独立未知成分ｐiについて微分して変形することで得られる下記式（３０）に基づいて算出される。
【数３０】
【０１０７】
ここで、式（３０）中のｐiは、質量行列［Ｍ］及び剛性行列［Ｋ］中にｑ個ある独立未知成分のうちの一つであるとする（ｉ＝１〜ｑ）。固有ベクトル感度は、上記式（２９）と固有ベクトルの質量行列に関する正規化の方程式（下記式（３１））と固有モードの直交性条件に関する式（下記式（３２））を利用しながら微分操作することによって下記式（３３）のように得られる。
【数３１】
【数３２】
【数３３】
【０１０８】
このようにして得られた１次感度を用いて、ｑ個の独立未知成分を微少量Δｐi（ｉ＝１〜ｑ）だけ増加させた場合の固有角振動数を二乗した値の増分が、現在の質量行列［Ｍ］及び剛性行列［Ｋ］に基づいて上記式（２９）の一般的固有値問題を解いて算出される固有角振動数Ωtを二乗した値と加振試験に基づいてステップＳ１５において算出された目標固有角振動数Ωttを二乗した値との差分（Ωtt2−Ωt2）になるように、且つｑ個の独立未知成分を微少量Δｐiだけ増加させた場合の固有モードベクトルの増分が、現在の質量行列［Ｍ］及び剛性行列［Ｋ］に基づいて上記式（２９）の一般的固有値問題を解いて算出される固有モードベクトル｛ψ｝tと加振試験に基づいてステップＳ１５において算出された目標固有モードベクトル｛ψ｝ttとの差分（｛ψ｝tt−｛ψ｝t）になるように、上記微少量の増分Δｐiが算出される（下記式（３４））。
【数３４】
【０１０９】
次いで、このようにして算出された独立未知成分の微少増分Δｐiを用いて上記式（２８）と同様な式により独立未知成分の更新を行う。そして、このような操作を繰り返すことで、質量行列［Ｍ］及び剛性行列［Ｋ］の最適化を行い、質量行列［Ｍ］及び剛性行列［Ｋ］の成分を同定する。
【０１１０】
なお、上記式（３４）の下式において、本実施形態では、左辺の列ベクトル（目標固有モードベクトル）｛ψtt｝には、実測定点自由度に対応する固有モード成分のみが配置され、仮想測定点自由度に対応する固有モード成分は配置されない。このため、本実施形態では、実測定点自由度に対応する固有モード成分についてのみ目標値への一致化を上記式（３４）によって行い、仮想測定点自由度に対応する固有モード成分については目標値への一致化を行わない。すなわち、本実施形態では、同定すべき質量行列［Ｍ］及び剛性行列［Ｋ］の固有値解析で得られる３ｎ個の成分から成る固有モードベクトルの中から実測定点の自由度に対応する自由度の成分のみを抜き出して列ベクトル｛ψt｝を作成して上記式（３４）の最適化を図る。
【０１１１】
このように、実測定点自由度に対応する固有モード成分についてのみ目標値への一致化を行っても、質量行列［Ｍ］及び剛性行列［Ｋ］を適切に最適化することができる。すなわち、本実施形態では、上記ステップＳ２３、Ｓ２４等で求められる測定点間の構造的結合状態等に基づく制約条件式が利用されている。このため、実測定点が固有モードの形状を適切に表現することができるように解析対象物１上に十分適切に配置されれば、実測定点間に配置される仮想測定点の変位や固有モード成分は実測定点の変位や固有モード成分に基づいて内挿補間的に表現できる。よって、実測定点自由度に対応する固有モード成分についてのみ目標値への一致化を行えば、実測定点間に配置される仮想測定点自由度に対応する固有モード成分は従属的に適切な値となるように調整されることになる。
【０１１２】
なお、このように実測定点自由度に対応する固有モード成分についてのみ一致化を行うことにより、解析対象物１上に設定された全測定点の自由度に対応する固有モード成分について一致化を行う場合に比べて、質量行列［Ｍ］及び剛性行列［Ｋ］の最適化解析の計算負荷が大幅に軽減され、よって同定解析の計算速度を向上させることができる。特に、上述したように測定点の一部を仮想測定点とすることから、加速度計２３等のセッティングに必要な時間も大幅に削減できることから、本実施形態では解析時間を大幅に縮めることができる。
【０１１３】
また、仮想測定点を配置することにより十分な自由度の大きさを確保できることから固有モードの形状を滑らかに表現することができるため、解析対象物１上に設定された全測定点の自由度に対応する固有モード成分について一致化を行う場合に対して、質量行列［Ｍ］及び剛性行列［Ｋ］の同定精度はほとんど低下しない。したがって、本実施形態によれば、質量行列［Ｍ］及び剛性行列［Ｋ］の同定精度をほとんど低下させることなく解析時間を極めて短縮することができる。
【０１１４】
このような最適化の概念に基づいて、ステップＳ３６では、上記式（３３）により固有各振動数Ω毎に各独立未知成分についての固有各振動数の１次感度が算出される。次いで、ステップＳ３７では、ステップＳ３６において算出された１次感度を用いて上記式（３４）により独立未知成分の微少増分Δｐiが算出され、算出された微少増分Δｐiを用いて上記式（２８）と同様な式により独立未知成分の更新が行われる。
【０１１５】
ステップＳ３８では、ステップＳ３７で更新された独立未知成分を含む質量行列［Ｍ］及び剛性行列［Ｋ］に基づいて算出された固有角振動数Ωt及び固有モードベクトル｛ψ｝tが目標固有角振動数Ωtt及び固有モードベクトル｛ψ｝ttと許容範囲内で一致しているか否かが判定される。このステップＳ３８での判断は、実験結果から得られたモード特性と計算によって得られたモード特性との相関性を示す値を計算し、この値が許容範囲を満足するか否かで判断する。
【０１１６】
ステップＳ３８において固有角振動数同士及び固有モードベクトル同士が許容範囲内で一致していないと判定された場合には、ステップＳ３６及びＳ３７が繰り返される。一方、固有角振動数同士及び固有モード同士が許容範囲内で一致していると判定された場合には、ステップＳ３９へ、すなわちゲイン調整処理へと進む。ステップＳ３６〜ステップＳ３８までの処理が上述した最適化処理に該当する。
【０１１７】
特性行列の同定処理では、次にゲイン調整処理が行われる。上述した反復計算による質量行列［Ｍ］と剛性行列［Ｋ］の同定法では、最適化の数学処理上の都合から、質量行列［Ｍ］の成分の最大絶対値を１とする正規化条件を保持しながら最適化処理が行われている。このため、上記最適化処理で得られた特性行列を用いて計算される加振試験の測定結果に対応する周波数応答関数は、加振試験で得られている周波数応答関数と共振峰の周波数（すなわち固有角振動数）と固有モードは一致するものの、ゲイン（振動振幅レベル）は一致しない。ただし、図１０に示したように、最適化処理で得られた特性行列を用いて計算された周波数応答関数（図１０中の実線Ｘ）と加振試験で得られている周波数応答関数（図１０中の破線Ｙ）とは、上記最適化処理で得られた質量行列［Ｍ］を適切に定数倍すると両者が一致するようなずれの関係となっている。そこで、ゲインも最適に一致するように、最適化操作を行う。
【０１１８】
ここで、特性行列に基づく周波数応答関数Ｘを加振試験に基づく周波数応答関数Ｙに一致させるためには、特定の周波数における両周波数応答関数の値を算出し、算出された値の比率に基づいてゲインを算出することが考えられる。ただし、共振周波数付近（図１０中のＺ）の周波数帯域では周波数応答関数の値の変化が大きいため、共振周波数付近以外の周波数帯域における両周波数応答関数の値の比率に基づいてゲインを算出するのが好ましい。
【０１１９】
そこで、本実施形態では、共振周波数付近を除く周波数帯域における下記式（３５）によって得られる周波数応答関数の値が、加振試験で実際に得られている周波数応答関数の値に最適に一致するようにゲイン調整パラメータαの値を決定する。なお、式（３５）において、｛ｆ（ω）｝は周波数ωにおける周波数応答関数を計算するための外力ベクトルであり、｛ｈ（ω）｝は周波数ωにおける周波数応答関数を表すベクトルである。
【数３５】
【０１２０】
ここで、加振試験では全ての測定点について実際に測定が行われるわけではないため、特性行列に基づく周波数応答関数が全ての測定点について求められるのにたいして、加振試験に基づく周波数応答関数は実測定点についてのみしか求められない。そこで、本実施形態では、特性行列に基づく周波数応答関数のうち実測定点に対応するものについてのみ加振試験に基づく周波数応答関数と一致するようにゲイン調整パラメータαの値を決定する。
【０１２１】
最後に、このようにして算出されたゲイン調整パラメータαを上記最適化処理で算出された質量行列［Ｍ］及び剛性行列［Ｋ］に乗算し、これによって算出された行列（α［Ｍ］及びα［Ｋ］）をそれぞれ最終同定結果としての質量行列と剛性行列とする。
【０１２２】
このようなゲイン調整の概念に基づいて、ステップＳ３９では、ステップ３６〜３８で同定された質量行列［Ｍ］及び剛性行列［Ｋ］に基づいて式（３５）を用いて周波数応答関数が算出される。次いで、ステップＳ４０では、ステップ３９で算出された周波数応答関数と、加振試験で得られている周波数応答関数とを比較して、ゲイン調整パラメータαを算出する。次いで、ステップＳ４１では、ステップ３６〜３８で同定された質量行列［Ｍ］及び剛性行列［Ｋ］にαを乗算して、最終同定結果としての質量行列と剛性行列を算出している（［Ｍ］＝α［Ｍ］、［Ｋ］＝α［Ｋ］）。
【０１２３】
なお、後述するように剛体特性の同定にあたっては質量行列［Ｍ］のみが必要であるため、上記実施形態では減衰行列［Ｃ］の同定を行っていない。しかしながら、減衰行列［Ｃ］の同定も可能であるため、以下では減衰行列［Ｃ］の同定法について簡単に説明する。
【０１２４】
減衰行列［Ｃ］の同定にあたっては、まず、減衰行列［Ｃ］の初期値を剛性行列［Ｋ］をコピーして［Ｃ］＝β［Ｋ］により作成する。そして、最適に一致化させたい第１次から所定の次数までのモード減衰比ζtについての最小二乗法によって最適なβを決定する。すなわち、下記式（３６）で計算されるモード減衰比ζtが、適合させるべき全ての次数について、振動試験に基づいてステップＳ１５で同定されたモード減衰比ζttに適合するように、補正係数βを算出する。
【数３６】
【０１２５】
このように算出された補正係数βを用いて、［Ｃ］＝β［Ｋ］により暫定的な減衰行列を作成する。次に、更なる最適化を図るために、このようにして作成された減衰行列を初期値として、減衰行列に基づいて算出されるモード減衰比ζtが振動試験に基づいて同定されたモード減衰比ζttに一致するように、１次感度解析を行う。この感度解析においても、一般に知られている各種の感度解析の手法を利用することができる。
【０１２６】
モード減衰比についての１次感度は、第ｔ次のモード減衰比ζtについて、減衰行列を構成する独立未知成分ｐi（ｉ＝１〜ｑ）に関して下記式（３７）で算出される。このように算出された１次感度を用いて、ｑ個の独立未知成分を微少量Δｐi（ｉ＝１〜ｑ）だけ増加させた場合のモード減衰比の増分が、現在の減衰行列［Ｃ］に基づいて算出されるモード減衰比ζtと加振試験に基づいて算出されたモード減衰比ζttとの差分（ζtt−ζt）になるように、上記微少量の増分Δｐiが算出される（下記式（３８））。
【数３７】
【数３８】
【０１２７】
次いで、このようにして算出された独立未知成分の微少増分Δｐiを用いて上記式（２８）と同様な式により独立未知成分の更新を行う。そして、このような操作を繰り返すことで、減衰行列［Ｃ］の最適化を行い、減衰行列［Ｃ］の成分を同定する。
【０１２８】
［剛体特性の同定］
次に、図３にステップ５０で示したように、上記同定処理に基づいて同定された質量行列［Ｍ］に基づいて剛体特性の算出が行われる。図１１は、本発明の質量行列［Ｍ］に基づく剛体特性の算出処理の処理手順を示すフローチャートである。
【０１２９】
ところで、上記同定処理に基づいて同定された質量行列［Ｍ］を上記式（１６）の左辺に代入して計算すると、その右辺に示される成分構成の６行６列の行列（剛体質量行列）が数値行列として得られる。この剛体質量行列の各成分は上記式（１８）に示した成分構成となっている。従って、このように算出された剛体質量行列に基づいて質量ｍを直接的に求めることができる。また、剛体運動表現点を座標原点にして剛体変位モード行列［Φ］を作成して上記式（１６）の左辺を計算すると、質量中心の座標（Ａ／ｍ、Ｂ／ｍ、Ｃ／ｍ）も直接算出される。
【０１３０】
一方、主慣性モーメントと慣性主軸の向きの算出は以下のように行われる。剛体変位モード行列［Φ］の第４〜第６列目の剛体回転運動ベクトルを、剛体運動表現点ではなく上述したように算出された質量中心の座標を座標原点として、上記式（１６）の左辺の演算を改めて行う。すると、演算結果の剛体質量行列［Ｍrigid］は下記式（３９）のようになる。
【数３９】
【０１３１】
式（３９）の４行４列から６行６列目までの、いわゆる慣性テンソル（３行３列）について、下記式（４０）の標準的固有値問題を解く。これにより算出された三つの固有値が解析対象物１の３つの主慣性モーメントの値、それらに対応して得られる長さ１に正規化された固有ベクトルがそれぞれの主慣性モーメントに対応する慣性主軸の方向余弦ベクトル（慣性主軸の向きを表す）である。すなわち、式（４０）の固有値問題の解として、下記式（４１）が得られる。
【数４０】
【数４１】
【０１３２】
このようなモード特性同定の概念に基づいて、ステップＳ５１では、上記同定処理に基づいて同定された質量行列［Ｍ］及び剛体運動表現点を座標原点にして作成された剛体変位モード行列［Φ］に基づいて、上記式（１６）を用いて質量ｍ及び質量中心の座標（Ａ／ｍ、Ｂ／ｍ、Ｃ／ｍ）が算出される。次いで、ステップＳ５２では、上記同定処理に基づいて同定された質量行列［Ｍ］及びステップＳ５１で算出された質量中心の座標を原点にして作成された剛体変位モード行列［Φ］に基づいて、上記式（１６）を用いて剛体質量行列［Ｍrigid］が算出される。ステップＳ５３では、ステップＳ５２で算出された剛体質量行列［Ｍrigid］に基づいて解析対象物１の３つの主慣性モーメントの値、及び各主慣性モーメントに対応する慣性主軸の方向余弦ベクトルが算出される。
【０１３３】
本発明の実施形態によれば、上述したように解析対象物の単点加振試験の測定データに基づいて多点の周波数応答関数から、解析対象物１の振動に関する特性行列の同定及び剛体特性の同定を行うことができる。したがって、本発明の実施形態では、実験的に特性行列の同定を行うことができるので、直接振動解析を行いたい解析対象物の物理モデルを表現することができる。したがって、部分構造合成法、最適設計、有限要素法との一体化及び各種シミュレーション等に直接的に適用することができる。
【０１３４】
特に、本発明の実施形態では、測定点のうち一部の測定点を仮想測定点とすることにより、実際に測定を行う実測定点の測定自由度よりも大きい自由度の特性行列が同定される。これにより、上述したように特性行列や剛体特性の同定精度をほとんど低下させることなく、解析時間を極めて短縮することができる。逆に考えると、実際に測定を行う実測定点の測定自由度と同じ自由度の特性行列を同定する場合に比べて、特性行列や剛体特性の同定精度を高めることができる。
【０１３５】
なお、上記実施形態では、測定点のうち一部の測定点を仮想測定点とすることで、実際の測定自由度よりも大きい自由度の特性行列を同定している。しかしながら、加振試験で実際に測定される測定自由度よりも同定される特性行列の自由度の方が大きければ本発明を適用可能である。
【０１３６】
従って、例えば、各測定点の測定自由度が２自由度（例えばｘ軸方向及びｙ軸方向のみ）であるｎ個の測定点に基づいて３ｎ行３ｎ列の特性行列を同定したり、各測定点の測定自由度が３自由度（例えば、ｘ軸方向、ｙ軸方向及びｚ軸方向）であるｎ個の測定点に基づいて６ｎ行６ｎ列の特性行列を同定したりすることも可能である。なお、６ｎ行６ｎ列の特性行列は、ｘ、ｙ、ｚ軸方向の並進運動に加えてｘ、ｙ、ｚ軸回りの回転運動を仮定して作成される特性行列である。
【０１３７】
換言すると、本発明は、実際に測定が行われる実測定自由度と実際には測定が行われない仮想測定自由度とを合わせた自由度の特性行列を、実際に測定が行われる実測定自由度の測定データに基づいて同定する場合にも適用可能である。
【図面の簡単な説明】
【０１３８】
【図１】本実施形態の解析装置のブロック図を示す図である。
【図２】解析対象物の多自由度系でのモデル化の例を示す図である。
【図３】解析処理のフローチャートである。
【図４】モード特性の同定処理のフローチャートである。
【図５】物理的制約条件式の算出処理の処理手順を示すフローチャートである。
【図６】構造的結合状態と零成分の配置との関係を説明するための図である。
【図７】剛体変位モード行列及び剛体質量行列を説明するための図である。
【図８】特性行列の同定処理の処理手順を示すフローチャートの一部である。
【図９】特性行列の同定処理の処理手順を示すフローチャートの一部である。
【図１０】ゲイン調整処理を説明するための図である。
【図１１】剛体特性の算出処理の処理手順を示すフローチャートである。
【符号の説明】
【０１３９】
１ 解析対象物
１０ 解析装置
２０ 測定部
２１ 力変換器
２２ 打撃ハンマ
２３ 加速度計
２４、２５ 増幅器
２６ 外部インタフェース部
３０ 解析部
３１ リムーバブルメディア装置
３３ 入力装置
３５ 出力装置
３７ 中央処理装置
３８ 記憶装置
【技術分野】
【０００１】
本発明は、振動解析方法及び振動解析装置に関する。
【背景技術】
【０００２】
自動車やエンジン等の複雑な構造物では、その剛体特性（質量、質量中心位置、三つの主慣性モーメント及びこれらに対応する三本の主軸の向き等）はその構造物の運動性能や振動騒音性能及び振動騒音制御の性能を左右する最大因子の一つである。事実、剛体特性の値は、設計支援のための運動解析、振動解析、防振機構解析、制御性能解析等、多岐にわたる動的挙動の解析、設計、最適化を開始するときに最も必要なパラメータである。従って、複雑な構造物の剛体特性を容易に実用的な精度で同定（推定）することは非常に重要である。
【０００３】
自動車等の構造物の剛体特性を同定する場合には、構造物をワイヤーで吊り下げることによる古典的な剛体特性の計測法を適用することは困難である。このため、現在では振り子ベンチ形式の計測装置が開発され、利用されている。しかし、この計測装置では、自動車を搭載する振り子ベンチ部分及び計測対象の自動車を剛体として扱う必要があるため、セッティングには手間がかかる。また、慣性モーメント計測時には、予め自動車の質量中心位置を同定しておいて、その質量中心が計測装置の回転中心軸線上に正確に位置するように自動車を設置する必要がある。このように、上記振り子ベンチ形式の計測装置には剛体特性を計測するにあたり、いくつかの問題点がある。
【０００４】
本願の発明者は、既に上記課題に対する解決方法として、「実験的特性行列同定法」と名付けた新しい方法を提案している（特許文献１）。この方法では、同定対象物を剛体として扱う必要は無く、むしろ弾性体としての第１次から２次又は３次までの適当な数の共振振動特性から剛体特性を求めることができる。また、加振試験として単点加振多点応答試験を１度実行するだけで、すべての剛体特性を同定することができる。
【０００５】
【特許文献１】特開２００１−３５０７４１号公報
【発明の開示】
【発明が解決しようとする課題】
【０００６】
ところで、上記本願発明者の提案した方法では、加振試験として単点加振多点応答試験が行われる。ここで、複雑な構造物に対して剛体特性を高い精度で同定するためには、その複雑さや大きさに応じて或る程度多くの測定点を設定する必要がある。
【０００７】
ところが、実際に測定を行う測定点の数が増大すると、測定用のセンサのキャリブレーションや計測システムのセッティング等を含む加振試験の測定作業にかかる時間が増大する。また、各測定点毎に周波数応答関数の同定が行われるため、実際に測定を行う測定点の数が増大するとコンピュータによる解析負荷が大きくなる。一方、測定点数を減らして加振試験を行うと、剛体特性の同定精度が低くなってしまう。
【０００８】
従って、上記問題点に鑑み、本発明の目的は、実際に測定を行う測定点を減少させても、同定精度をほとんど低下させることなく剛体特性を同定することができる振動解析方法又は振動解析装置を提供することにある。
【課題を解決するための手段】
【０００９】
上記課題を解決するために、第１の発明では、解析対象物上に設定された複数の測定点の座標データを要求する工程と、上記解析対象物上の少なくとも一点において加振する加振試験が行われた時に得られた各測定点における測定データに基づいてモード特性を算出する工程と、上記座標データに基づいて、上記解析対象物の振動を表す特性行列を構成する成分間の制約条件式を算出する工程と、上記測定データに基づいて算出されたモード特性の値に基づいて上記制約条件式を用いて上記特性行列を同定する工程とを具備する、振動解析方法において、上記加振試験で上記測定点において実際に測定される測定自由度よりも上記特性行列の自由度の方が大きく、よって上記測定データに基づいて算出されるモード特性の値の数の方が上記特性行列に基づいて算出されるモード特性の値の数よりも少なく、上記特性行列を同定する工程では、上記特性行列に基づいて算出されるモード特性のうち上記測定データに基づいて算出されたモード特性に対応するモード特性についてのみ、特性行列に基づいて算出されるモード特性の値が上記測定データに基づいて算出されたモード特性の値に一致するように上記特性行列の各成分の値を変更していくことで、上記特性行列を同定する。
【００１０】
上記課題を解決するために、第２の発明では、解析対象物上に設定された複数の実測定点及び仮想測定点の座標データを要求する工程と、上記解析対象物上の少なくとも一点において加振する加振試験が行われた時に実測定点においてのみ測定を行い、加振により得られた各実測定点における測定データに基づいてモード特性を算出する工程と、上記座標データに基づいて上記解析対象物の振動を表す特性行列を構成する成分間の制約条件式を算出する工程と、上記測定データに基づいて、算出されたモード特性の値に基づいて上記制約条件式を用いて上記特性行列を同定する工程とを具備する、振動解析方法において、上記特性行列は実測定点及び仮想測定点の両測定点における自由度に等しい自由度で作成され、上記特性行列を同定する工程では、上記特性行列に基づいて算出されるモード特性のうち上記測定データに基づいて算出されたモード特性に対応するモード特性についてのみ、特性行列に基づいて算出されるモード特性の値が上記測定データに基づいて算出されたモード特性の値に一致するように上記特性行列の各成分の値を変更していくことで、上記特性行列を同定する。
【００１１】
第３の発明では、第２の発明において、上記仮想測定点はその仮想測定点に最も近接する測定点が実測定点となるように設けられる。
【００１２】
第４の発明では、第２又は第３の発明において、測定点間の構造的結合状態の入力を要求する工程をさらに具備し、上記制約条件式を算出する工程では、上記座標データに加えて測定点間の構造的結合状態に基づいて制約条件式を算出する。
【００１３】
第５の発明では、第２〜第４のいずれか一つの発明において、上記特性行列を同定する工程で同定される特性行列は質量行列であり、上記特性行列を同定する工程において同定された質量行列に基づいて剛体特性を算出する工程をさらに具備する。
【００１４】
上記課題を解決するために、第６の発明では、解析対象物上の測定点における位置、速度又は加速度の変位を測定する測定器と、上記解析対象物上の少なくとも一点において加振を行う加振器と、上記測定器及び加振器からそれぞれ測定データ及び加振データが入力されると共にこれらデータに基づいて演算処理を行う解析部とを具備する振動解析装置において、上記解析部は、上記加振器での加振試験により得られた各測定点における測定データに基づいてモード特性を算出し、各測定点の座標データに基づいて上記解析対象物の振動を表す特性行列を構成する成分間の制約条件式を算出し、上記測定データに基づいて算出されたモード特性の値に基づいて上記制約条件式を用いて上記特性行列を同定し、上記加振試験で上記測定点において実際に測定される測定自由度よりも上記特性行列の自由度の方が大きく、よって上記測定データに基づいて算出されるモード特性の値の数の方が上記特性行列に基づいて算出されるモード特性の値の数よりも少なく、上記特性行列を同定する際には、上記特性行列に基づいて算出されるモード特性のうち上記測定データに基づいて算出されたモード特性に対応するモード特性についてのみ、特性行列に基づいて算出されるモード特性の値が上記測定データに基づいて算出されたモード特性の値に一致するように上記特性行列の各成分の値を変更していくことで、上記特性行列を同定する。
【００１５】
上記課題を解決するために、第７の発明では、解析対象物上に設定された実測定点及び仮想測定点のうち実測定点における位置、速度又は加速度の変位を測定する測定器と、上記解析対象物上の少なくとも一点において加振を行う加振器と、上記測定器及び加振器からそれぞれ測定データ及び加振データが入力されると共にこれらデータに基づいて演算処理を行う解析部とを具備する振動解析装置において、上記解析部は、上記加振器による加振試験により得られた実測定点における測定データに基づいてモード特性を算出し、上記実測定点及び仮想測定点の座標データに基づいて特性行列を構成する成分間の制約条件式を算出し、上記測定データに基づいて算出されたモード特性の値に基づいて上記制約条件式を用いて上記特性行列を同定し、上記特性行列は実測定点及び仮想測定点の両測定点における自由度に等しい自由度で作成され、上記特性行列を同定する際には、上記特性行列に基づいて算出されるモード特性のうち上記測定データに基づいて算出されたモード特性に対応するモード特性についてのみ、特性行列に基づいて算出されるモード特性の値が上記測定データに基づいて算出されたモード特性の値に一致するように上記特性行列の各成分の値を変更していくことで、上記特性行列を同定する。
【００１６】
第８の発明では、第７の発明において、上記仮想測定点はその仮想測定点に最も近接する測定点が実測定点となるように設けられる。
【００１７】
第９の発明では、第７又は第８の発明において、上記制約条件式を算出する際には、上記座標データに加えて測定点間の構造的結合状態に基づいて制約条件式が算出される。
【００１８】
第１０の発明では、第７〜第９のいずれか一つの発明において、上記同定される特性行列は質量行列であり、上記解析部は、上記特性行列を同定する工程において同定された質量行列に基づいて剛体特性を算出する。
【発明の効果】
【００１９】
本発明によれば、実際に測定を行う測定点を減少させても、推定精度をほとんど低下させることなく剛体特性を推定することができる。
【発明を実施するための最良の形態】
【００２０】
以下、図面を参照して本発明の実施形態について詳細に説明する。
【００２１】
図１は、本実施形態の解析装置のブロック図である。図１において、解析装置１０は、解析対象物１の剛体特性等を解析的に算出可能な装置であり、測定部２０と解析部３０とを具備する。測定部２０は、先端に力変換器２１を備えた打撃ハンマ２２と、解析対象物１に取り付けられる複数の計測器（例えば、加速度計）２３ａ〜２３ｃと、力変換器２１及び各計測器２３ａ〜２３ｃに接続された増幅器２４、２５ａ〜２５ｃと、これら増幅器２４、２５ａ〜２５ｃに接続された外部インタフェース部２６とを具備する。
【００２２】
一方、解析部３０は、リムーバブルメディアへの書き込み等を行うリムーバブルメディア装置（ＲＭ装置）３１と、リムーバブルメディア装置３１を制御するリムーバブルメディア制御部（ＲＭ制御部）３２と、ユーザからの入力を受け付ける入力装置３３と、入力装置３３を制御する入力制御部３４と、ディスプレイ等への出力を行う出力装置３５と、出力装置３５を制御する出力制御部３６と、各種演算処理を行う中央処理装置３７と、記憶装置３８と、これら外部インタフェース部２６、リムーバブルメディア制御部３２、入力制御部３４、出力制御部３６、中央処理装置３７、記憶装置３８等を相互に接続するバス３９とを具備する。
【００２３】
打撃ハンマ２２は、解析対象物１を加振するのに用いられる。加振時に加えられた力の大きさは、打撃ハンマ２２に設けられた力変換器２１によって電気信号に変換される。この電気信号は、増幅器２４で増幅された後に外部インタフェース２６に出力される。なお、本実施形態では加振器として打撃ハンマ２２を用いて解析対象物１に非定常波を与えているが、加振器や加振方法はこれに限定されるものではない。解析対象物１を加振することができると共に解析対象物１に加えられた力の大きさが分かれば、打撃ハンマ２２の代わりに如何なる加振器を用いても良い。従って、加振器として、油圧式や圧電式等の加振器、及び音圧や磁界を用いた非接触加振器等を使用することができる。また、非定常波の代わりに、定常波や周期波等の各種の加振波形の加振を行ってもよい。さらに、本実施形態では、単点加振試験を行っているが、多点加振試験を行ってもよい。
【００２４】
加速度計２３ａ〜２３ｃはそれぞれ解析対象物１上の測定点に配置され、測定点において解析対象物１に発生した加速度の大きさを測定し、これを電気信号に変換して出力する。この電気信号は、増幅器２５ａ〜２５ｃで増幅された後に外部インタフェース２６に出力される。測定点は、注目したい次数までの固有モードを幾何学的に充分に表現することができる程度に適切な分布及び数となるようにユーザにより決定される。なお、本実施形態では、計測器として３軸加速度計２３ａ〜２３ｃを用いているが、渦電流やレーザ光線を利用した非接触変位計、レーザドップラ速度計、ひずみゲージ等、測定点における各種パラメータ（位置、速度、加速度等）の変位を計測することができれば加速度計２３ａ〜２３ｃの代わりに如何なる計測器を用いても良い。
【００２５】
外部インタフェース部２６は、力変換器２１及び加速度計２３から増幅器２４、２５を介して入力された電気信号をデジタル信号に変換し、このデジタル信号をバス３９に出力する。リムーバブルメディア装置３１は、リムーバブルメディア制御部３２によって制御される。リムーバブルメディア装置３１では、リムーバブルメディア２に記録されたプログラム（例えば、後述する解析を行うプログラム）等を解析部３０にインストールしたり、解析結果をリムーバブルメディア２に記録したりすることができる。なお、リムーバブルメディア２としては、ＦＤ、ＣＤ、ＤＶＤ、ＭＯ等、如何なる記憶媒体を用いても良い。
【００２６】
入力装置３３は、キーボード、マウス及びデジタイザ等の入力機器であり、入力制御部３４によって制御される。この入力装置３３により、ユーザは解析対象物の形状や各測定点の座標データ等、各種データを入力することができる。出力装置３５は、ディスプレイやプリンタ等の出力機器であり、出力制御部３６によって制御される。出力装置３５には、入力装置３３に入力された各種データ、プログラム実行のためのメニュー、解析結果等が出力される。
【００２７】
中央処理装置３７は、バス３９を介して、外部インタフェース部２６、リムーバブルメディア制御部３２、入力制御部３４、出力制御部３６及び記憶装置３８と接続され、これら制御部等に指示を与え、また、記憶装置１２９に格納され各プログラムに従いフーリエ変換や実験的特性行列同定法等の各種処理を行う。また、記憶装置３８は、ＲＡＭ（ランダムアクセスメモリ）やＲＯＭ（リードオンリメモリ）等から構成され、プログラムや、各プログラム実行中の一時的なデータ、解析結果等が格納される。
【００２８】
次に、このように構成された解析装置１０によって行われる解析処理について説明する。
【００２９】
まず、本実施形態の解析装置１０によって行われる解析処理についての基本的な概念について説明する。例えば、図２（ａ）に示したような解析対象物を多自由度系でモデル化する。特に、図２（ａ）に示した例では、計測点がｎ個あり、各計測点における自由度がｘ、ｙ、ｚ方向の３自由度であるため、解析対象物を自由度３ｎの多自由度系でモデル化している。この場合、質量行列、粘性減衰行列及び剛性行列をそれぞれ［Ｍ］、［Ｃ］、及び［Ｋ］とおくと、モデル化した解析対象物の運動方程式は下記式（１）で表される。なお、下記式（１）では粘性減衰を仮定しているが、必要に応じて構造減衰を仮定してもよいし、両者の減衰を同時に組み込んだ運動方程式を用いてもよい。
【数１】
【００３０】
ここで、式（１）中の｛ψ｝は測定点と定めた各点におけるｘ方向変位ｘi、ｙ方向変位ｙi、ｚ方向変位ｚiを縦に並べてできるベクトル、｛ｆ｝は測定点と定めた各点におけるｘ方向の加振力ｆｘi、ｙ方向の加振力ｆｙi、ｚ方向の加振力ｆｚiを縦に並べてできるベクトルを表している（下記式（２）参照）。
【数２】
【００３１】
なお、質量行列［Ｍ］、粘性減衰行列［Ｃ］及び剛性行列［Ｋ］（以下、これらをまとめて「特性行列」という）は３ｎ行３ｎ列の対称行列（すなわち、行数及び列数が解析対象物の測定自由度３ｎに等しい対称行列）であり、ｉは測定点の番号を示している。また、式（１）中の｛ψ’｝は変位を表すベクトルの一階微分、すなわち測定点の速度を表すベクトルであり、｛ψ’’｝は変位を表すベクトルの二階微分、すなわち測定点の加速度を表すベクトルである。本明細書では、行列を［］で、ベクトルを｛｝で表す。
【００３２】
上記運動方程式（１）において、質量行列［Ｍ］等の特性行列を求めると、解析対象物１の振動の様子を把握することができると共に、解析対象物１の剛体特性（例えば、解析対象物１の質量、質量中心位置、三つの主慣性モーメント及びこれらに対応する三つの主軸の向き等）を求めることができる。そこで、本実施形態では、打撃ハンマ２２によって解析対処物１が加振されたときの各測定点での測定データに基づいて、解析対象物１の特性行列を求めることとしている。
【００３３】
加振時の測定データに基づく解析対象物１の特性行列の同定及び剛体特性の算出は、図３に示したように大別して四つの処理工程によって行われる。
【００３４】
まず、ステップＳ１０に示したように、単点加振試験時の測定データに基づいて、モード特性の同定が行われる。すなわち、打撃ハンマ２２によって解析対象物１の一点が打撃（加振）されると共に、この加振力が力変換器２１によって測定され且つこの加振によって発生した振動が各加速度計２３によって測定される。その後、このようにして力変換器２１及び加速度計２３によって測定された測定データから各測定点について周波数応答関数及びコヒーレンス関数が算出される。次いで、算出された周波数応答関数を用いて、一般的な方法で、固有振動数、固有モードベクトル及びモード減衰比（以下、これらをまとめて「モード特性」という）を同定する。
【００３５】
次いで、ステップＳ２０では、特性行列成分間の物理的制約条件式の算出が行われる。すなわち、特性行列の特定の成分の値や特定の成分間の制約条件は、各種力学原理、各測定点の座標、測定点間の構造的結合状態等に基づいて、加振試験の測定結果とは無関係に求められる。例えば、測定点間の構造的結合状態に基づいて特性行列中の零成分の配置を或る程度特定することができる。また、各種力学原理や測定点の座標等に基づいて、特性行列中の非零成分間の関係（制約条件）を或る程度特定することができる。
【００３６】
このように、特性行列成分間の制約条件式を算出することにより、上記ステップＳ１０において求められたモード特性に基づいて同定しなければならない特性行列中の成分の数が減少する。例えば質量行列［Ｍ］が３ｎ行３ｎ列の行列である場合、この質量行列［Ｍ］中には３ｎ（３ｎ＋１）／２個の成分が含まれている（実際の成分の数は３ｎ×３ｎ個であるが、行列の対称性を考慮すると３ｎ（３ｎ＋１）／２個であるということができる）。ここで、この質量行列［Ｍ］の成分間の制約条件等が求められていない場合には、上記ステップＳ１０において求められたモード特性に基づいて同定しなければならない質量行列［Ｍ］中の成分の数は３ｎ（３ｎ＋１）／２ということになり、現実的には質量行列［Ｍ］の全ての成分を適切に同定することは困難である。
【００３７】
しかしながら、本ステップＳ２０では質量行列［Ｍ］の成分間の制約条件式が求められ、よって質量行列［Ｍ］中の零成分の配置や、質量行列［Ｍ］中の特定の成分と別の特定の成分との関係式が求められる。従って、質量行列［Ｍ］成分間の制約条件式を算出することにより、上記ステップＳ１０において求められたモード特性に基づいて同定しなければならない質量行列［Ｍ］中の成分の数がかなり減少する。具体的には、質量行列［Ｍ］中の全ての成分の数３ｎ（３ｎ＋１）／２個から、質量行列［Ｍ］中の零成分の数及び作成された制約条件式の数を減算した数の成分（以下、「独立未知成分」という）を上記ステップＳ１０において求められたモード特性に基づいて同定することになる。
【００３８】
次いで、ステップＳ３０では、数学的最適化法により質量行列［Ｍ］及び剛性行列［Ｋ］の独立未知成分の同定が行われる。すなわち、質量行列［Ｍ］及び剛性行列［Ｋ］と固有角振動数Ωt及び固有モードベクトル｛ψ｝tとの関係は下記式（３）のように表せる（なお、ｔは次数を表す）。
【数３】
【００３９】
式（３）から分かるように、質量行列［Ｍ］及び剛性行列［Ｋ］の各成分を変数とみなしてその値を変化させると、これに伴って固有角振動数Ωt及び固有モードベクトル｛ψ｝tの値も変化することになる。従って、本ステップＳ３０では、上記式（３）によって求められる固有角振動数Ωt及び固有モードベクトル｛ψ｝tが上記ステップＳ１０で求められた固有角振動数及び固有モードベクトルに一致するように、ステップＳ２０で求められた質量行列［Ｍ］及び剛性行列［Ｋ］の成分間の制約条件式を満足させながらこれら質量行列［Ｍ］及び剛性行列［Ｋ］の独立未知成分の値が徐々に変化せしめられる。固有角振動数同士及び固有モードベクトル同士がほぼ一致したときの質量行列［Ｍ］及び剛性行列［Ｋ］が、図２（ａ）に示したように多自由度系でモデル化された解析対象物の質量行列［Ｍ］及び剛性行列［Ｋ］として同定される。
【００４０】
以下では、このように固有角振動数同士及び固有モードベクトル同士を一致させることを固有角振動数及び固有モードベクトルの「一致化」として説明し、固有角振動数同士及び固有モードベクトル同士を一致させることで質量行列［Ｍ］及び剛性行列［Ｋ］の各成分の最適な値を算出・同定することを質量行列［Ｍ］及び剛性行列［Ｋ］の「最適化」として説明する。
【００４１】
なお、このように質量行列［Ｍ］及び剛性行列［Ｋ］の最適化を行うための数学的手法としては、「ニュートン法類に属する最急降下法」等、多くの手法を適用することが可能であり、最適化が効率よく且つ安定して実行できれば如何なる手法を用いることも可能である。
【００４２】
次いで、ステップＳ４０では、ステップＳ３０において同定された質量行列［Ｍ］に基づいて、解析対象物１の剛体特性、すなわち解析対象物１の質量、質量中心位置、三つの主慣性モーメント及びこれらに対応する三つの主軸の向きが求められる。
【００４３】
ところで、上述した特性行列及び剛体特性の同定法では、適当な数ｎの測定点を解析対象物１に設定して、加振試験時においてこれら各測定点でｘ、ｙ、ｚ方向の振動について周波数応答関数を計測し、それに対応した自由度の特性行列を同定することで、剛体特性を推定している。解析対象物１が複雑な構造物であるような場合には、その形状の複雑さや大きさに応じて測定点の数を或る程度多く設定する必要がある。
【００４４】
ところが、上述した特性行列及び剛体特性の同定法においては、解析対象物１上への測定点の設定から周波数応答関数の算出までの間に、測定点の設定、加速度計のキャリブレーション、計測システムのセッティング等の作業が行われ、よって測定点の設定から周波数応答関数の算出までが最も手間と時間のかかる作業である。従って、測定点数が増大すると、それに伴って周波数応答関数の算出までにかかる手間と時間が大きく増大し、その結果、特性行列及び剛体特性の同定にかかる時間が増大してしまう。また、周波数応答関数は中央処理装置３７において行われる同定処理によって求められるため、測定点数が増大すると中央処理装置３７における作業負荷が増大してしまい、周波数応答関数の収束速度が大きく低下してしまう。したがって、特性行列及び剛体特性の同定を迅速に行うという観点からは、実際に測定を行う測定点の数が少ないことが好ましい。
【００４５】
そこで、本発明の実施形態では、測定点数を減少させずに、すなわち特性行列の自由度を小さくすることなく、実際に測定を行う測定点の数を減少させることとしている。すなわち、解析対象物１上に従来と同様なほどの十分な数の測定点を設定すると共に、そのうちの一部の測定点については実際に測定を行って周波数応答関数を算出し（以下、このような測定点を「実測定点」という）、残りの測定点については実際には測定を行わないようにしている（以下、このような測定点を「仮想測定点」という）。換言すると、本実施形態では、実際に測定が行われる測定点の自由度よりも大きな自由度の特性行列を同定することとしている。
【００４６】
本実施形態では、このように実測定点と仮想測定点とを設定した上で、上記ステップＳ３０において、数学的最適化法により質量行列［Ｍ］及び剛性行列［Ｋ］の独立未知成分の同定を行うにあたり、固有モードベクトルのうち実測定点に対応する固有モードベクトルの成分のみを一致化させるようにしている。
【００４７】
ここで、本願発明者は、このように仮想測定点という概念を導入した場合であっても、測定点間の構造的結合状態等に基づいて適切に制約条件式が求められている場合には、固有モードベクトルのうち実測定点に対応する成分のみを一致化させることができれば、仮想測定点に対応する成分も従属的に（または自動的に）一致化させることができ、その結果、適切な質量行列［Ｍ］及び剛性行列［Ｋ］を同定することができることを発見した。換言すると、本願発明者は、固有モードベクトルの一致化をするにあたって、必ずしも全ての固有モードベクトルの成分について一致化させなくても、十分な精度で特性行列の最適化を行うことができることを発見したのである。
【００４８】
以下では、本発明の実施形態に係る解析装置１０によって行われる解析処理の手順について詳細に説明する。なお、本実施形態における解析処理も、基本的に図３に示した処理手順に従って行われる。従って、以下では、図４〜図７を参照して、上記解析処理の処理手順に対応する形で本実施形態の解析処理の処理手順について説明する。
【００４９】
［モード特性の同定］
まず、ステップＳ１０に示したように、単点加振試験による測定結果に基づいてモード特性の同定が行われる。図４は、本発明のモード特性の同定処理の処理手順を示すフローチャートである。モード特性の同定は、加振試験を行って周波数応答関数等を算出する振動計測処理と、算出された周波数応答関数等に基づいてモード特性を同定する特性同定処理とに大別される。
【００５０】
図４に示した処理手順の実行の前に、ユーザによる計測システムのセッティングが行われる。計測システムのセッティングにあたっては、解析対象物１がユーザによりゴム紐やタイヤチューブ等の弾性体により支持されるように配置される。これにより、解析対象物１の周辺自由状態が近似的に実現される。なお、解析対象物１の支持は、擬似的に周辺自由状態が実現されれば、如何なる態様で行われてもよい。
【００５１】
図４に示したようなモード特性の同定処理においては、まず、ステップＳ１１において、中央処理装置３７は、ユーザに解析対象物１の形状・寸法等の初期データを入力するように出力画面３５を介して要求する。このような初期データには、実測定点及び仮想測定点を含む全ての測定点の座標データが含まれる。
【００５２】
ユーザは実測定点及び仮想測定点を解析対象物１上の任意の位置に設定する。このとき、仮想測定点は、できるだけ仮想測定点同士が隣り合うことの無いように、すなわちその仮想測定点に最も近接する測定点が実測定点となるように設定される。換言すると、仮想測定点は、少なくとも一つの或る方向において実測定点の間に配置されるように設定される。また、仮想測定点は、仮想測定点数が全測定点数の１／２以上にならないように（好ましくは、１／３以上にならないように）設定される。
【００５３】
図２（ｂ）は実測定点及び仮想測定点の設定の例を示している。図２（ｂ）中に丸で囲まれた測定点は実際に計測が行われる実測定点を示しており、丸で囲まれていない測定点は実際には計測が行われない仮想測定点を示している。すなわち、図示した例では、３、６、７、１１番の測定点においては実際に測定が行われない。
【００５４】
ユーザは、このような条件のもと各測定点の設定位置を決定すると共に、中央処理装置３７の上記要求に応じて全ての測定点の座標データを入力する。また、実測定点の設定位置には加速度計２３が取り付けられる。
【００５５】
次に、ユーザは、打撃ハンマ２２によって加速度計２３が取り付けられた位置の中のいずれかの１つの位置（単点）で解析対象物１に打撃（加振）を加える。
ユーザが打撃ハンマ２２によって解析対象物１に加えた打撃は力変換器２１によって測定され、その測定データは増幅器２４を介して外部インタフェース部２６に入力される。また、この打撃によって解析対象物１に発生した振動も加速度計２３によって測定され、その測定データは増幅器２５を介して外部インタフェース部２６に入力される。入力されたこれら測定データは中央処理装置３７に取り込まれる。このようにして単点加振試験が行われる。
【００５６】
次に、ステップＳ１３において、測定データに基づいて周波数応答関数等が算出される。すなわち、中央処理装置３７は、これら力変換器２１及び加速度計２３によって測定された測定データをフーリエ変換（例えば、高速フーリエ変換）して周波数スペクトルになおす。その後、中央処理装置３７は、解析対象物１への入力に関する周波数スペクトル（力変換器２１の出力に基づく周波数スペクトル）及び解析対象物１からの出力に関する周波数スペクトル（加速度計２３の出力に基づく周波数スペクトル）からパワースペクトルとクロススペクトルとを算出し、これらパワースペクトル及びクロススペクトルから周波数応答関数及びコヒーレンス関数を算出する。これら周波数応答関数及びコヒーレンス関数は各実測定点毎に算出される。
【００５７】
ここで、コヒーレンス関数（関連度関数）は、入力と出力と内積であり、入出力の相関を示すものである。その値が０である場合には、入力と出力とは無相関であり、１に近づくにつれて相関性が増す。このコヒーレンス関数は、後述する周波数応答関数に基づくモード特性の同定にも用いられる（ステップＳ１５で使用）。
【００５８】
なお、より信頼性の高いデータを得る観点から、この単点加振試験を複数回繰り返し、個々のパワースペクトルおよびクロススペクトルを平均し、その平均値を用いるようにすることが好ましい。
【００５９】
上記ステップＳ１１〜Ｓ１３の処理が振動計測処理であり、単点加振試験により得られた測定データに基づいて各実計測点における周波数応答関数及びコヒーレント関数が算出される。
【００６０】
このように、周波数応答関数が算出された後で、算出された周波数応答関数等に基づいてモード特性を同定する特性同定処理が行われる。
【００６１】
特性同定処理では、まず、ステップＳ１４に示したように、中央処理装置３７が、ユーザにモード特性の同定を行う周波数帯域（以下、「同定周波数帯域」という）を入力するように出力画面３５を介して要求する。ユーザは入力装置３３から同定周波数帯域を入力する。同定周波数帯域は、その帯域幅が広いほどその帯域内に現れる共振峰の数が多くなる。このため、同定周波数帯域を広げると、多くの次数のモード特性を同定することができる。
【００６２】
次いで、ステップＳ１５において、ステップＳ１３で算出された周波数応答関数及びコヒーレント関数に基づいて、所定の次数までのモード特性が同定される。この所定の次数は、上記ステップＳ１４で入力された同定周波数帯域応じて変化する。なお、モード特性の同定法としては一般に知られている各種のモード同定法を利用することができ、このようなモード同定法としては、例えば、多点参照法（ポリリファレンス法、poly-referencemethod）、偏分反復法、プロリー法（Prony's method）、サークルフィット法（circle fit method）等が挙げられる。
【００６３】
以下では、周波数応答関数に基づいてモード特性を同定する際の基本的な原理について簡単に説明する。減衰を比例粘性減衰とした振動系の限定された周波数帯域（すなわち、上記同定周波数帯域）に関する周波数応答関数は近似的に下記式（４）で表せる。
【数４】
ここで、式（４）中のＹは慣性項定数パラメータ、Ｚは剰余項定数パラメータであり、また同定周波数帯域内に存在する共振峰の数をａとする。さらに、そのａ個の共振峰に対応する固有角振動数をΩt、モード減衰比をζt、加振点ｑ、応答点（測定点）ｐの固有モード成分をそれぞれψqt、ψpt（ｔは次数）とする。
【００６４】
このように式（４）で表した周波数応答関数Ｈpqは、固有角振動数Ωt、モード減衰比ζt、固有モード成分ψpt、ψqtを適切に選択すれば、ステップＳ１３において振動試験に基づく周波数応答関数にほぼ一致する。そこで、本モード特性同定では、式（４）によって算出される周波数応答関数ＨpqがステップＳ１３において算出された周波数応答関数にほぼ一致するような固有角振動数Ωt、モード減衰比ζt、固有モード成分ψpt、ψqtを最適化法等によって算出し、これらモード特性を同定するようにしている。また、このモード特性の同定の過程で上記コヒーレンス関数が用いられる。
【００６５】
上記ステップＳ１４及びＳ１５の処理が特性同定処理であり、周波数応答関数に基づいてモード特性が同定される。
【００６６】
［特性行列成分間の物理的制約条件式の同定］
次に、図３にステップＳ２０に示したように、特性行列成分間の物理的制約条件式の算出が行われる。図５は、本発明の特性行列成分間の物理的制約条件式の算出処理の処理手順を示すフローチャートである。
【００６７】
特性行列成分間の物理的制約条件式の算出にあたって、中央処理装置３７は、ステップＳ２１において、測定点間の構造的結合状態（コネクティビティ）を入力するように出力画面３５を介して要求する。ここで、測定点間の構造的結合状態とは、或る測定点と別の測定点とが構造的に結びついているか否か、すなわち或る測定点と他の測定点とが構造的に隣り合う関係にあるか否かの状態を意味する。
【００６８】
例えば図２（ａ）に示した構造について考慮すると、測定点間の構造的結合状態は図２（ｃ）のように表される。図中、互いに構造的に結合しているとされる測定点同士は直線によって互いに結ばれている。図２（ａ）に示した例では、例えば測定点１は測定点２、３、４と互いに構造的に結合されており、測定点３は測定点１、２、４及び５と互いに構造的に結合されている。なお、本実施形態では、測定点間の構造的結合状態はユーザの主観により測定点の配置や測定点間の距離等に基づいて定められて、ユーザにより入力装置３３を用いて入力される。しかしながら、例えば、ユーザに解析対象物１の形状及び測定点の配置を入力させて、この入力されたデータに基づいて中央処理装置３７により測定点間の構造的結合状態が算出されるようにしてもよい。
【００６９】
次いで、ステップＳ２２では、ステップ２１において入力された測定点間の構造的結合状態に基づいて質量行列［Ｍ］及び剛性行列［Ｋ］の零成分の配置が特定される。以下、零成分の配置について説明する。
【００７０】
ここで、図６（ａ）に示したような単純な均一剛体要素と回転バネで直列に結合されたモデルを考える。各剛体要素は、均一形状で長さがそれぞれｌi、質量中心が中央にありその質量がｍi、回転バネ定数がｋiであるとして、曲げ振動表現の自由度（上記実施形態における測定点数に相当）を点１〜５番の並進変位とする。この場合の質量行列［Ｍ］及び剛性行列［Ｋ］はそれぞれ下記式（５）、（６）のように表される。
【数５】
【数６】
【００７１】
ここで、式（５）、（６）中の＊は零成分でないこと、すなわち非零成分であることを意味するものである。このように、図６に示したようなモデルでは、質量行列［Ｍ］については３重対角行列、剛性行列［Ｋ］については５重対角行列として表される。図６（ｂ）、（ｃ）を参照してこの理由について説明する。なお、図６では、剛体要素間の寸法を大きくして回転弾性対偶を描いているが、剛体要素間の寸法はないものとして考える。すなわち、剛体要素Ｉの右端と剛体要素ＩＩの左端とは同一位置であるとして考える。
【００７２】
図６（ｂ）を参照して、剛体要素Ｉの質量中心（剛体要素Ｉの図心）について、自由体図によって運動方程式を作成すると、回転運動については式（７）のように、並進運動については式（８）のように表される。これら式（７）、（８）を行列式でまとめて表示すると、式（９）のようになる。
【数７】
【数８】
【数９】
【００７３】
剛体要素Ｉの両端の並進変位自由度と質量中心での変位自由度との変換則は下記式（１０）のように表せる。そこで、式（１０）を式（９）に代入して、両辺左側から［Ｔ1］Tを乗じると、下記式（１１）のようになる。なお、［Ｔ1］TはＴ1の転置行列を表す。
【数１０】
【数１１】
ここで、剛体要素Ｉが均一形状であって均一材料密度分布であるとすると、下記式（１２）のように表せる。なお、式（１２）においてｍ1は剛体要素１の質量である。
【数１２】
【００７４】
同様に、図６（ｃ）を参照して、剛体要素ＩＩの質量中心（要素の図心）について自由体図によって上記剛体要素Ｉの場合と同様に運動方程式を作成して整理すると、下記式（１３）が導かれる。
【数１３】
【００７５】
この様に順々に剛体要素III、剛体要素IV、…と運動方程式を求めて、全ての運動方程式を行列表現でまとめて、全系の運動方程式を作成すると、下記式（１４）に示したようにその質量行列［Ｍ］は３重対角行列、剛性行列［Ｋ］は５重対角行列になる。このように、質量行列［Ｍ］及び剛性行列｛Ｋ］の非零成分の配置は、測定点間の構造的結合状態によって決定することができる。
【数１４】
【００７６】
上述の力学原理に基づいて、図２（ｃ）に示したような測定点間の構造的結合状態となっている解析対象物のモデルを考えると、例えば剛性行列［Ｋ］の零成分の配置は下記式（１５）のようになる。
【数１５】
【００７７】
上記式（１５）の剛性行列［ｋ］中の［Ｎ］は３行３列の非零成分となる可能性のある小行列部分を、［０］は３行３列の零成分となる小行列部分を表している。ここで、［Ｎ］、［０］が３行３列の小行列となっているのは各測定点における測定自由度がｘ、ｙ、ｚ方向の３自由度であるためであり、例えば各測定点における測定自由度がｘ、ｙ、ｚ方向に加えてｘ軸、ｙ軸、ｚ軸回りの回転方向を加えた６自由度である場合には［Ｎ］、［０］は６行６列の小行列となる。
【００７８】
図２（ｃ）に示した例では、測定点１は測定点２、３に物理的に結合しており、これら測定点２又は３は測定点４、５に物理的に結合している。従って、式（１７）に示したように、剛性行列［Ｋ］では、測定点１に対応する１列目には１〜５行目に非零小行列成分が配置されることになる。このように、ユーザがデザインする構造的結合モデルに従って特性行列の零成分配置が力学原理に基づいてコンピュータにより自動的に決定される。
【００７９】
次に、ステップＳ２３において、３次元構造物の質量行列［Ｍ］については、如何なるｎ自由度のモデルについての質量行列［Ｍ］も下記式（１６）左辺の演算によって、構造物中に任意に設定された剛体運動表現点としての１点の６自由度変位に関する６行６列の剛体質量行列Ｍrigid］に変換できるという力学原理に基づいて、質量行列［Ｍ］の非零成分間の連立１次方程式が物理的制約条件式として作成される。ここで、式（１６）中の行列［Φ］は下記式（１７）のように表される３ｎ行６列の行列（以下、「剛体変位モード行列」という）であり、剛体質量行列[Ｍrigid］は下記式（１８）のように表される６行６列の行列である。
【数１６】
【数１７】
【数１８】
【００８０】
なお、剛体質量行列中のＩxx、Ｉyy、Ｉzzは剛体運動表現点を通り各座標軸と並行な軸回りの慣性モーメント、Ｉyx、Ｉzx、Ｉyzはそれら３軸相互間についての慣性乗積である。Ａ、Ｂ、Ｃは、剛体運動表現点と重心とのズレに起因する値であり、解析対象物１の質量をｍ、重心座標を（ｘG、ｙG、ｚG）、剛体運動表現点座標を（ｘR、ｙR、ｚR）とすると、Ａ＝ｍ（ｘG−ｘR）、Ｂ＝ｍ（ｙG−ｙR）、Ｃ＝ｍ（ｚG−ｚR）である。［Ｍrigid］は、上記式（１８）に示すように対称行列である。
【００８１】
以下、上記剛体変位モード行列［Φ］、剛体質量行列［Ｍrigid］及び上記式（１８）の導出方法について簡単に説明する。
上記剛体変位モード行列［Φ]は、６本の互いに独立な剛体運動ベクトルで構成される。具体的に、入力データの一部である測定点座標値を利用して、剛体運動表現点がｘ軸方向、ｙ軸方向、ｚ軸方向に単位並進変位するときの構造物上の測定点の変位を表すベクトルをそれぞれ第１列、第２列、第３列に当てはめ、剛体運動表現点がその点を通ってｘ軸と平行な軸回り、ｙ軸と平行な軸回り及びｚ軸と平行な軸回りに微少角変位する剛体角変位運動に関する構造物上の測定点の並進変位を線形近似表現するベクトルを第４列、第５列、第６列に当てはめることで剛体変位モード行列［Φ］は容易に作成することができる。
【００８２】
例えば、図７を参考に、剛体中に設定した剛体運動表現点座標をＲ（ｘR，ｙR，ｚR）、測定点ｉ番のｘ、ｙ、ｚ座標を（ｘi，ｙi，ｚi）として、剛体が単にｘ軸方向へδRxだけ並進変位するとそれに伴って測定点ｉ番の変位は、下記式（１９）のように、δRxに定数で構成される変換行列（この場合は列ベクトル）を乗じることで表現することができる。
【数１９】
【００８３】
剛体が剛体運動表現点で単にｘ軸に平行な軸回りに角度θxだけ右手に回転することによる測定点ｉ番の変位は幾何学から下記式（２１）のように代数的に表現することができる。
【数２０】
【００８４】
ここで、今対象としている機械振動の並進振幅や角変位振動は小さい値であることを考慮して、ｓｉｎθｘ≒θｘ、ｃｏｓθｘ≒１の近似を十分適用することができるため、上記式（２０）は下記式（２１）と変形することができる。式（１９）と同様に座標データに基づく定数で構成される変換行列を運動表現点の角変位に乗じることで測定点ｉの変位が表現できる。
【数２１】
【００８５】
このようにして６自由度の全ての剛体運動を考えると、剛体変位モード行列［Φ］の成分構成は測定点ｉのｘ、ｙ、ｚ方向変位成分についてのみ明記した形式として、下記式（２２）のように作成することができる。
【数２２】
【００８６】
従って、各測定点の座標データを入力することにより、上記剛体変位モード行列［Φ］を求めることができる。なお、フーリエ変換、モード特性同定法等の実験的同定法では、すべてデータ処理法（信号処理法）に過ぎず、幾何学的データ（または構造的データ）としての測定点座標を本質的に必要する「実験的同定法」はこれまで皆無であったと考えることができ、本手法はこの点についても大きな特徴を有する。
【００８７】
次に、剛体質量行列［Ｍrigid］について説明する。上記式（１６）の右辺の剛体質量行列［Ｍrigid］は、如何なる対象構造物に関しても上記式（１８）のような独特な成分配置関係となることが知られている。例えば、剛体質量行列［Ｍrigid］には、１行１列成分と２行２列成分及び３行３列成分は常に同一の値でありその値は解析対象物の質量の値であること、２行１列成分と１行２列成分は常に零となること、５行１列成分と４行２列成分とは絶対値は同一で符号が互いに逆の値となること、等の特徴的な成分配置関係がある。
【００８８】
このように表される剛体変位モード行列［Φ］、剛体質量行列［Ｍrigid］及び質量行列［Ｍ］の間には、上記式（１６）のような関係がある。ここで、式（１８）の非零成分（ｍ、Ａ、Ｉxx等）の値そのものは同定前には未知であるが、これら非零成分間の関係を利用して、式（１６）からその左辺の質量行列［Ｍ］の成分間に関する１１本の等式制約条件式を自動生成することができる。
【００８９】
なお、解析対象物１の質量（重さ）等の剛体特性値が既知として与えられていれば、これを式（１８）に代入して質量行列成分間の制約条件式を生成するようにしてもよい。このように既知の剛体特性値が得られれば、この既知の剛体特性値については、他の成分で表す必要が無くなるという利点がある。
【００９０】
このようにして生成された連立方程式は、互いに独立な方程式よりも未知数である非零成分の数の方が必ず多くなる。従って、このように生成された連立方程式のみからは質量行列［Ｍ］の非零成分を全て特定することはできない。しかしながら、方程式の数だけの未知数（非零成分）を他の未知数（非零成分）の従属変数とすることができる。以下の説明では、本ステップＳ２３等において求められた制約条件式により特性行列の他の成分に従属する形で表現される成分を従属成分と、制約条件式によっても他の成分に従属する形では表現されない成分を独立未成分として説明する。
【００９１】
このような等式制約条件式を上記手法により自動生成すると、例えば下記式（２３）のような等式が作成される。
【数２３】
ここで、ｍi,jは質量行列のｉ行ｊ列の成分を表し、左辺の成分が従属成分、右辺の各項の成分が独立未知成分である。なお、どの成分が独立未知成分に割り当てられ、どの成分が従属成分に割り当てられるかはプログラミングに依存する。例えば、式（２３）の最上段の式の場合、ｍ1,1を独立未知成分とすることもでき、この場合ｍ2,1が従属成分となる。
【００９２】
次に、ステップＳ２４において、３次元構造物の剛性行列［Ｋ］と減衰行列［Ｃ］について、剛体運動状態ではその構造物のどの位置にもひずみが生じず、そのひずみに基づく応力も全ての部分において零であるという基本的な力学原理に基づいて、剛性行列［Ｋ］及び減衰行列［Ｃ］の非零成分間に連立１次方程式が物理的制約条件式として作成される。上記力学原理からは、下記式（２４）が導き出される。なお、式（２４）の右辺の［０］はｎ行６列の零行列である。
【数２４】
【００９３】
この式（２４）により、各特性行列中の非零成分に関する連立１次方程式が作成される。これにより、ステップＳ２３と同様に、剛性行列［Ｋ］及び粘性行列［Ｃ］についても、成分間に関する複数の等式制約条件式が生成され、生成された方程式の数だけ独立未知成分を減少させることができる。
【００９４】
このように、ステップＳ２１〜Ｓ２４では、物理原則等に基づいて各特性行列毎にその各成分間の関係を示す制約条件式が求められる。これにより、各特性行列中に含まれる独立未知成分の数を減少させると共に力学（物理学）的成分間関連特性を付与することができ、上記ステップＳ３０における最適化処理における最適化精度を高めることができる。
【００９５】
［特性行列の同定］
次に、図３にステップＳ３０で示したように、質量行列［Ｍ］、剛性行列［Ｋ］及び粘性行列［Ｃ］中の独立未知成分の同定、従ってこれら特性行列の同定が行われる。上述したように、質量行列［Ｍ］及び剛性行列［Ｋ］の最適化を行うための数学的手法としては多くの手法を採用可能であり、以下では図８及び図９を参照してその一つの例について説明する。図８及び図９は、本発明の特性行列の同定処理の処理手順を示すフローチャートである。
【００９６】
本実施形態の特性行列の同定処理は、三つの処理に大別される。一つ目の処理は、独立未知成分として適当な値が代入された質量行列［Ｍ］が正定値行列になるように、また、独立未知成分として適当な値が代入された剛性行列［Ｋ］が非負定値行列となるように、独立未知成分の値を修正する正定値・非負定値化処理である。二つ目の処理は、上記式（３）によって求められる固有角振動数Ωt及び固有モードベクトル｛ψ｝tが上記ステップＳ１０で求められた固有角振動数及び固有モードベクトル（以下、これらをそれぞれ「目標固有角振動数Ωtt」、「目標固有モードベクトル｛ψ｝ttという）に一致するように、質量行列［Ｍ］及び剛性行列［Ｋ］の独立未知成分の値が徐々に変化させる最適化処理である。そして、三つ目の処理は、正規化されている質量行列［Ｍ］及び剛性行列［Ｋ］にゲイン平行移動係数αを乗算して、最終的に質量行列［Ｍ］及び剛性行列［Ｋ］を同定するゲイン調整処理である。図８及び図９では、これら三つの大きな処理毎に破線で囲まれている。
【００９７】
特性行列の同定処理では、まずステップＳ３１において、質量行列［Ｍ］及び剛性行列［Ｋ］の各成分のうちステップＳ２０（具体的にはステップＳ２４）で物理的制約条件式を作成した結果独立未知成分となった成分に初期値として乱数を代入する。また、代入された乱数に基づいてステップＳ２０において作成された制約条件式を用いて従属成分を算出し、初期値としての質量行列［Ｍ］と剛性行列［Ｋ］を数値行列として作成する。
【００９８】
しかしながら、この初期値としての質量行列［Ｍ］と剛性行列［Ｋ］は乱数から作成されているため、物理学的に最低限満足しなければならない条件をも満足しない場合がある。すなわち、質量や回転慣性は必ず正の値であることから質量行列は正定値行列となる。また、同様に剛性行列も基本的に正の値となるが、剛体運動ではひずみが零であるので固有値に零のものが最大６個まで存在する（剛体変位自由度だけ固有値零が存在する）ことから剛性行列は非負定値行列となる。しかしながら、上記初期値としての質量行列［Ｍ］と剛性行列［Ｋ］はこのような条件を満足しない場合がある。
【００９９】
そこで、本実施形態では、ステップＳ３２〜３４において、質量行列については正定値行列化を行い、剛性行列については非負定値行列化を行う。換言すると、質量行列［Ｍ］に関しては全ての３ｎ個の固有値が正値となるように質量行列［Ｍ］の独立未知成分を修正し、剛性行列［Ｋ］の全ての３ｎ個の固有値が負値とならないように剛性行列［Ｋ］の独立未知成分を修正する。ここでは、代表として質量行列の正定値行列化手法について説明する。
【０１００】
まず、質量行列［Ｍ］をその成分の絶対値最大の値に対して正規化をした上で、すなわち質量行列［Ｍ］の成分の最大絶対値が１になるようにした上で、標準的固有値問題である下記式（２５）を解析して、３ｎ個の全ての固有値λと固有ベクトル｛ｘ｝を求めて、負値の固有値を選び出す。
【数２５】
このようにして算出された固有値λの中に負値のものがなければ、質量行列［Ｍ］はこの時点の各成分の値で正定値行列である。一方、固有値λに負値のものがある場合には、これら負値の固有値λを正値へ変化させる必要がある。このために、負値の固有値λを正値へ変化させるための１次感度解析を行う。この感度解析としては、例えば、最急降下法、共役勾配法といったニュートン法に代表されるニュートン法類の手法等、一般に知られている各種の感度解析の手法を利用することができる。以下では、その一例について示す。質量行列［Ｍ］内の独立未知成分ｐiに関する固有値λについての１次感度ｓt,iは上記式（２５）を独立未知成分ｐiで微分した下記式（２６）で算出される。
【数２６】
【０１０１】
ここで、式（２６）中のλtと｛ｘ｝tはそれぞれｔ次の固有値と固有ベクトルを示している。また、ｐiは質量行列［Ｍ］中にｑ個ある独立未知成分のうちの一つであるとする（ｉ＝１〜ｑ）。上記式（２６）で算出される１次感度ｓt,iを用いて、ｑ個の独立未知成分を微少量Δｐi（ｉ＝１〜ｑ）だけ増加させた場合の固有値λtの増分Δλtは近似的に下記式（２７）のように見積もることができる。
【数２７】
【０１０２】
そこで、負値の固有値λ全てについてこの式を生成する。このとき固有値の増分Δλtには、例えば、その負値の固有値の絶対値よりも大きい正値が代入される。そして、このようにして生成された式を満たすような独立未知成分の微少増分Δｐiを算出し、このようにして算出された独立未知成分の微少増分Δｐiを用いて独立未知成分の更新を行う。このような操作を繰り返すことにより、質量行列［Ｍ］の正定値化を行うことができる。なお、剛性行列［Ｋ］の非負定値化についても、上記質量行列［Ｍ］の正定値化と同様に行われる。
【０１０３】
このような正定値化の概念に基づいて、ステップＳ３２では、上記式（２５）に基づいて算出された全ての固有値λtのうち負値のものが選定される。次いで、ステップＳ３３では、上記式（２６）を用いてステップＳ３２で選定された負値の固有値λ毎に各独立未知成分についての固有値の１次感度ｓt,iが算出される。次いで、ステップＳ３４では、ステップＳ３３において算出された１次感度を用いて上記式（２７）により独立未知成分の微少増分Δｐiが算出され、算出された微少増分Δｐiを用いて下記式（２８）により独立未知成分の更新が行われる。
【数２８】
【０１０４】
次いで、ステップＳ３５では、ステップＳ３４において独立未知成分の更新が行われた結果、質量行列［Ｍ］の正定値化及び剛性行列［Ｋ］の非負定値化が実現されたか否かが判定され、質量行列［Ｍ］の正定値化及び剛性行列［Ｋ］の非負定値化が実現されていないと判定された場合にはステップＳ３２〜ステップ３４が繰り返される。一方、ステップＳ３５において、質量行列［Ｍ］の正定値化及び剛性行列［Ｋ］の非負定値化が実現されたと判定された場合には、ステップＳ３６へ、すなわち最適化処理へと進む。ステップＳ３２〜ステップＳ３５までの処理が上述した正定値・非負定値化処理に該当する。
【０１０５】
特性行列の同定処理では、次に最適化処理が行われる。最適化処理では、下記式（２９）の一般的固有値問題を解いて算出される固有角振動数（不減衰固有角振動数）Ωtと固有モードベクトル｛ψ｝tを、実験で得られた固有角振動数（目標固有角振動数Ωtt）及び固有モードベクトル（目標固有モードベクトル｛ψ｝tt）に一致化させる処理が行われる。このため、本実施形態では、これら固有角振動数同士及び固有モードベクトル同士を一致化させるべく、１次感度解析を行う。この感度解析においても、一般に知られている各種の感度解析の手法を利用することができる。
【数２９】
【０１０６】
固有角振動数Ωについての１次感度は、第ｔ次の固有角振動数Ωtと固有モードベクトル｛ψ｝tについて上記式（２９）を独立未知成分ｐiについて微分して変形することで得られる下記式（３０）に基づいて算出される。
【数３０】
【０１０７】
ここで、式（３０）中のｐiは、質量行列［Ｍ］及び剛性行列［Ｋ］中にｑ個ある独立未知成分のうちの一つであるとする（ｉ＝１〜ｑ）。固有ベクトル感度は、上記式（２９）と固有ベクトルの質量行列に関する正規化の方程式（下記式（３１））と固有モードの直交性条件に関する式（下記式（３２））を利用しながら微分操作することによって下記式（３３）のように得られる。
【数３１】
【数３２】
【数３３】
【０１０８】
このようにして得られた１次感度を用いて、ｑ個の独立未知成分を微少量Δｐi（ｉ＝１〜ｑ）だけ増加させた場合の固有角振動数を二乗した値の増分が、現在の質量行列［Ｍ］及び剛性行列［Ｋ］に基づいて上記式（２９）の一般的固有値問題を解いて算出される固有角振動数Ωtを二乗した値と加振試験に基づいてステップＳ１５において算出された目標固有角振動数Ωttを二乗した値との差分（Ωtt2−Ωt2）になるように、且つｑ個の独立未知成分を微少量Δｐiだけ増加させた場合の固有モードベクトルの増分が、現在の質量行列［Ｍ］及び剛性行列［Ｋ］に基づいて上記式（２９）の一般的固有値問題を解いて算出される固有モードベクトル｛ψ｝tと加振試験に基づいてステップＳ１５において算出された目標固有モードベクトル｛ψ｝ttとの差分（｛ψ｝tt−｛ψ｝t）になるように、上記微少量の増分Δｐiが算出される（下記式（３４））。
【数３４】
【０１０９】
次いで、このようにして算出された独立未知成分の微少増分Δｐiを用いて上記式（２８）と同様な式により独立未知成分の更新を行う。そして、このような操作を繰り返すことで、質量行列［Ｍ］及び剛性行列［Ｋ］の最適化を行い、質量行列［Ｍ］及び剛性行列［Ｋ］の成分を同定する。
【０１１０】
なお、上記式（３４）の下式において、本実施形態では、左辺の列ベクトル（目標固有モードベクトル）｛ψtt｝には、実測定点自由度に対応する固有モード成分のみが配置され、仮想測定点自由度に対応する固有モード成分は配置されない。このため、本実施形態では、実測定点自由度に対応する固有モード成分についてのみ目標値への一致化を上記式（３４）によって行い、仮想測定点自由度に対応する固有モード成分については目標値への一致化を行わない。すなわち、本実施形態では、同定すべき質量行列［Ｍ］及び剛性行列［Ｋ］の固有値解析で得られる３ｎ個の成分から成る固有モードベクトルの中から実測定点の自由度に対応する自由度の成分のみを抜き出して列ベクトル｛ψt｝を作成して上記式（３４）の最適化を図る。
【０１１１】
このように、実測定点自由度に対応する固有モード成分についてのみ目標値への一致化を行っても、質量行列［Ｍ］及び剛性行列［Ｋ］を適切に最適化することができる。すなわち、本実施形態では、上記ステップＳ２３、Ｓ２４等で求められる測定点間の構造的結合状態等に基づく制約条件式が利用されている。このため、実測定点が固有モードの形状を適切に表現することができるように解析対象物１上に十分適切に配置されれば、実測定点間に配置される仮想測定点の変位や固有モード成分は実測定点の変位や固有モード成分に基づいて内挿補間的に表現できる。よって、実測定点自由度に対応する固有モード成分についてのみ目標値への一致化を行えば、実測定点間に配置される仮想測定点自由度に対応する固有モード成分は従属的に適切な値となるように調整されることになる。
【０１１２】
なお、このように実測定点自由度に対応する固有モード成分についてのみ一致化を行うことにより、解析対象物１上に設定された全測定点の自由度に対応する固有モード成分について一致化を行う場合に比べて、質量行列［Ｍ］及び剛性行列［Ｋ］の最適化解析の計算負荷が大幅に軽減され、よって同定解析の計算速度を向上させることができる。特に、上述したように測定点の一部を仮想測定点とすることから、加速度計２３等のセッティングに必要な時間も大幅に削減できることから、本実施形態では解析時間を大幅に縮めることができる。
【０１１３】
また、仮想測定点を配置することにより十分な自由度の大きさを確保できることから固有モードの形状を滑らかに表現することができるため、解析対象物１上に設定された全測定点の自由度に対応する固有モード成分について一致化を行う場合に対して、質量行列［Ｍ］及び剛性行列［Ｋ］の同定精度はほとんど低下しない。したがって、本実施形態によれば、質量行列［Ｍ］及び剛性行列［Ｋ］の同定精度をほとんど低下させることなく解析時間を極めて短縮することができる。
【０１１４】
このような最適化の概念に基づいて、ステップＳ３６では、上記式（３３）により固有各振動数Ω毎に各独立未知成分についての固有各振動数の１次感度が算出される。次いで、ステップＳ３７では、ステップＳ３６において算出された１次感度を用いて上記式（３４）により独立未知成分の微少増分Δｐiが算出され、算出された微少増分Δｐiを用いて上記式（２８）と同様な式により独立未知成分の更新が行われる。
【０１１５】
ステップＳ３８では、ステップＳ３７で更新された独立未知成分を含む質量行列［Ｍ］及び剛性行列［Ｋ］に基づいて算出された固有角振動数Ωt及び固有モードベクトル｛ψ｝tが目標固有角振動数Ωtt及び固有モードベクトル｛ψ｝ttと許容範囲内で一致しているか否かが判定される。このステップＳ３８での判断は、実験結果から得られたモード特性と計算によって得られたモード特性との相関性を示す値を計算し、この値が許容範囲を満足するか否かで判断する。
【０１１６】
ステップＳ３８において固有角振動数同士及び固有モードベクトル同士が許容範囲内で一致していないと判定された場合には、ステップＳ３６及びＳ３７が繰り返される。一方、固有角振動数同士及び固有モード同士が許容範囲内で一致していると判定された場合には、ステップＳ３９へ、すなわちゲイン調整処理へと進む。ステップＳ３６〜ステップＳ３８までの処理が上述した最適化処理に該当する。
【０１１７】
特性行列の同定処理では、次にゲイン調整処理が行われる。上述した反復計算による質量行列［Ｍ］と剛性行列［Ｋ］の同定法では、最適化の数学処理上の都合から、質量行列［Ｍ］の成分の最大絶対値を１とする正規化条件を保持しながら最適化処理が行われている。このため、上記最適化処理で得られた特性行列を用いて計算される加振試験の測定結果に対応する周波数応答関数は、加振試験で得られている周波数応答関数と共振峰の周波数（すなわち固有角振動数）と固有モードは一致するものの、ゲイン（振動振幅レベル）は一致しない。ただし、図１０に示したように、最適化処理で得られた特性行列を用いて計算された周波数応答関数（図１０中の実線Ｘ）と加振試験で得られている周波数応答関数（図１０中の破線Ｙ）とは、上記最適化処理で得られた質量行列［Ｍ］を適切に定数倍すると両者が一致するようなずれの関係となっている。そこで、ゲインも最適に一致するように、最適化操作を行う。
【０１１８】
ここで、特性行列に基づく周波数応答関数Ｘを加振試験に基づく周波数応答関数Ｙに一致させるためには、特定の周波数における両周波数応答関数の値を算出し、算出された値の比率に基づいてゲインを算出することが考えられる。ただし、共振周波数付近（図１０中のＺ）の周波数帯域では周波数応答関数の値の変化が大きいため、共振周波数付近以外の周波数帯域における両周波数応答関数の値の比率に基づいてゲインを算出するのが好ましい。
【０１１９】
そこで、本実施形態では、共振周波数付近を除く周波数帯域における下記式（３５）によって得られる周波数応答関数の値が、加振試験で実際に得られている周波数応答関数の値に最適に一致するようにゲイン調整パラメータαの値を決定する。なお、式（３５）において、｛ｆ（ω）｝は周波数ωにおける周波数応答関数を計算するための外力ベクトルであり、｛ｈ（ω）｝は周波数ωにおける周波数応答関数を表すベクトルである。
【数３５】
【０１２０】
ここで、加振試験では全ての測定点について実際に測定が行われるわけではないため、特性行列に基づく周波数応答関数が全ての測定点について求められるのにたいして、加振試験に基づく周波数応答関数は実測定点についてのみしか求められない。そこで、本実施形態では、特性行列に基づく周波数応答関数のうち実測定点に対応するものについてのみ加振試験に基づく周波数応答関数と一致するようにゲイン調整パラメータαの値を決定する。
【０１２１】
最後に、このようにして算出されたゲイン調整パラメータαを上記最適化処理で算出された質量行列［Ｍ］及び剛性行列［Ｋ］に乗算し、これによって算出された行列（α［Ｍ］及びα［Ｋ］）をそれぞれ最終同定結果としての質量行列と剛性行列とする。
【０１２２】
このようなゲイン調整の概念に基づいて、ステップＳ３９では、ステップ３６〜３８で同定された質量行列［Ｍ］及び剛性行列［Ｋ］に基づいて式（３５）を用いて周波数応答関数が算出される。次いで、ステップＳ４０では、ステップ３９で算出された周波数応答関数と、加振試験で得られている周波数応答関数とを比較して、ゲイン調整パラメータαを算出する。次いで、ステップＳ４１では、ステップ３６〜３８で同定された質量行列［Ｍ］及び剛性行列［Ｋ］にαを乗算して、最終同定結果としての質量行列と剛性行列を算出している（［Ｍ］＝α［Ｍ］、［Ｋ］＝α［Ｋ］）。
【０１２３】
なお、後述するように剛体特性の同定にあたっては質量行列［Ｍ］のみが必要であるため、上記実施形態では減衰行列［Ｃ］の同定を行っていない。しかしながら、減衰行列［Ｃ］の同定も可能であるため、以下では減衰行列［Ｃ］の同定法について簡単に説明する。
【０１２４】
減衰行列［Ｃ］の同定にあたっては、まず、減衰行列［Ｃ］の初期値を剛性行列［Ｋ］をコピーして［Ｃ］＝β［Ｋ］により作成する。そして、最適に一致化させたい第１次から所定の次数までのモード減衰比ζtについての最小二乗法によって最適なβを決定する。すなわち、下記式（３６）で計算されるモード減衰比ζtが、適合させるべき全ての次数について、振動試験に基づいてステップＳ１５で同定されたモード減衰比ζttに適合するように、補正係数βを算出する。
【数３６】
【０１２５】
このように算出された補正係数βを用いて、［Ｃ］＝β［Ｋ］により暫定的な減衰行列を作成する。次に、更なる最適化を図るために、このようにして作成された減衰行列を初期値として、減衰行列に基づいて算出されるモード減衰比ζtが振動試験に基づいて同定されたモード減衰比ζttに一致するように、１次感度解析を行う。この感度解析においても、一般に知られている各種の感度解析の手法を利用することができる。
【０１２６】
モード減衰比についての１次感度は、第ｔ次のモード減衰比ζtについて、減衰行列を構成する独立未知成分ｐi（ｉ＝１〜ｑ）に関して下記式（３７）で算出される。このように算出された１次感度を用いて、ｑ個の独立未知成分を微少量Δｐi（ｉ＝１〜ｑ）だけ増加させた場合のモード減衰比の増分が、現在の減衰行列［Ｃ］に基づいて算出されるモード減衰比ζtと加振試験に基づいて算出されたモード減衰比ζttとの差分（ζtt−ζt）になるように、上記微少量の増分Δｐiが算出される（下記式（３８））。
【数３７】
【数３８】
【０１２７】
次いで、このようにして算出された独立未知成分の微少増分Δｐiを用いて上記式（２８）と同様な式により独立未知成分の更新を行う。そして、このような操作を繰り返すことで、減衰行列［Ｃ］の最適化を行い、減衰行列［Ｃ］の成分を同定する。
【０１２８】
［剛体特性の同定］
次に、図３にステップ５０で示したように、上記同定処理に基づいて同定された質量行列［Ｍ］に基づいて剛体特性の算出が行われる。図１１は、本発明の質量行列［Ｍ］に基づく剛体特性の算出処理の処理手順を示すフローチャートである。
【０１２９】
ところで、上記同定処理に基づいて同定された質量行列［Ｍ］を上記式（１６）の左辺に代入して計算すると、その右辺に示される成分構成の６行６列の行列（剛体質量行列）が数値行列として得られる。この剛体質量行列の各成分は上記式（１８）に示した成分構成となっている。従って、このように算出された剛体質量行列に基づいて質量ｍを直接的に求めることができる。また、剛体運動表現点を座標原点にして剛体変位モード行列［Φ］を作成して上記式（１６）の左辺を計算すると、質量中心の座標（Ａ／ｍ、Ｂ／ｍ、Ｃ／ｍ）も直接算出される。
【０１３０】
一方、主慣性モーメントと慣性主軸の向きの算出は以下のように行われる。剛体変位モード行列［Φ］の第４〜第６列目の剛体回転運動ベクトルを、剛体運動表現点ではなく上述したように算出された質量中心の座標を座標原点として、上記式（１６）の左辺の演算を改めて行う。すると、演算結果の剛体質量行列［Ｍrigid］は下記式（３９）のようになる。
【数３９】
【０１３１】
式（３９）の４行４列から６行６列目までの、いわゆる慣性テンソル（３行３列）について、下記式（４０）の標準的固有値問題を解く。これにより算出された三つの固有値が解析対象物１の３つの主慣性モーメントの値、それらに対応して得られる長さ１に正規化された固有ベクトルがそれぞれの主慣性モーメントに対応する慣性主軸の方向余弦ベクトル（慣性主軸の向きを表す）である。すなわち、式（４０）の固有値問題の解として、下記式（４１）が得られる。
【数４０】
【数４１】
【０１３２】
このようなモード特性同定の概念に基づいて、ステップＳ５１では、上記同定処理に基づいて同定された質量行列［Ｍ］及び剛体運動表現点を座標原点にして作成された剛体変位モード行列［Φ］に基づいて、上記式（１６）を用いて質量ｍ及び質量中心の座標（Ａ／ｍ、Ｂ／ｍ、Ｃ／ｍ）が算出される。次いで、ステップＳ５２では、上記同定処理に基づいて同定された質量行列［Ｍ］及びステップＳ５１で算出された質量中心の座標を原点にして作成された剛体変位モード行列［Φ］に基づいて、上記式（１６）を用いて剛体質量行列［Ｍrigid］が算出される。ステップＳ５３では、ステップＳ５２で算出された剛体質量行列［Ｍrigid］に基づいて解析対象物１の３つの主慣性モーメントの値、及び各主慣性モーメントに対応する慣性主軸の方向余弦ベクトルが算出される。
【０１３３】
本発明の実施形態によれば、上述したように解析対象物の単点加振試験の測定データに基づいて多点の周波数応答関数から、解析対象物１の振動に関する特性行列の同定及び剛体特性の同定を行うことができる。したがって、本発明の実施形態では、実験的に特性行列の同定を行うことができるので、直接振動解析を行いたい解析対象物の物理モデルを表現することができる。したがって、部分構造合成法、最適設計、有限要素法との一体化及び各種シミュレーション等に直接的に適用することができる。
【０１３４】
特に、本発明の実施形態では、測定点のうち一部の測定点を仮想測定点とすることにより、実際に測定を行う実測定点の測定自由度よりも大きい自由度の特性行列が同定される。これにより、上述したように特性行列や剛体特性の同定精度をほとんど低下させることなく、解析時間を極めて短縮することができる。逆に考えると、実際に測定を行う実測定点の測定自由度と同じ自由度の特性行列を同定する場合に比べて、特性行列や剛体特性の同定精度を高めることができる。
【０１３５】
なお、上記実施形態では、測定点のうち一部の測定点を仮想測定点とすることで、実際の測定自由度よりも大きい自由度の特性行列を同定している。しかしながら、加振試験で実際に測定される測定自由度よりも同定される特性行列の自由度の方が大きければ本発明を適用可能である。
【０１３６】
従って、例えば、各測定点の測定自由度が２自由度（例えばｘ軸方向及びｙ軸方向のみ）であるｎ個の測定点に基づいて３ｎ行３ｎ列の特性行列を同定したり、各測定点の測定自由度が３自由度（例えば、ｘ軸方向、ｙ軸方向及びｚ軸方向）であるｎ個の測定点に基づいて６ｎ行６ｎ列の特性行列を同定したりすることも可能である。なお、６ｎ行６ｎ列の特性行列は、ｘ、ｙ、ｚ軸方向の並進運動に加えてｘ、ｙ、ｚ軸回りの回転運動を仮定して作成される特性行列である。
【０１３７】
換言すると、本発明は、実際に測定が行われる実測定自由度と実際には測定が行われない仮想測定自由度とを合わせた自由度の特性行列を、実際に測定が行われる実測定自由度の測定データに基づいて同定する場合にも適用可能である。
【図面の簡単な説明】
【０１３８】
【図１】本実施形態の解析装置のブロック図を示す図である。
【図２】解析対象物の多自由度系でのモデル化の例を示す図である。
【図３】解析処理のフローチャートである。
【図４】モード特性の同定処理のフローチャートである。
【図５】物理的制約条件式の算出処理の処理手順を示すフローチャートである。
【図６】構造的結合状態と零成分の配置との関係を説明するための図である。
【図７】剛体変位モード行列及び剛体質量行列を説明するための図である。
【図８】特性行列の同定処理の処理手順を示すフローチャートの一部である。
【図９】特性行列の同定処理の処理手順を示すフローチャートの一部である。
【図１０】ゲイン調整処理を説明するための図である。
【図１１】剛体特性の算出処理の処理手順を示すフローチャートである。
【符号の説明】
【０１３９】
１ 解析対象物
１０ 解析装置
２０ 測定部
２１ 力変換器
２２ 打撃ハンマ
２３ 加速度計
２４、２５ 増幅器
２６ 外部インタフェース部
３０ 解析部
３１ リムーバブルメディア装置
３３ 入力装置
３５ 出力装置
３７ 中央処理装置
３８ 記憶装置
【特許請求の範囲】
【請求項１】
解析対象物上に設定された複数の測定点の座標データを要求する工程と、
上記解析対象物上の少なくとも一点において加振する加振試験が行われた時に得られた各測定点における測定データに基づいてモード特性を算出する工程と、
上記座標データに基づいて、上記解析対象物の振動を表す特性行列を構成する成分間の制約条件式を算出する工程と、
上記測定データに基づいて算出されたモード特性の値に基づいて上記制約条件式を用いて上記特性行列を同定する工程とを具備する、振動解析方法において、
上記加振試験で上記測定点において実際に測定される測定自由度よりも上記特性行列の自由度の方が大きく、よって上記測定データに基づいて算出されるモード特性の値の数の方が上記特性行列に基づいて算出されるモード特性の値の数よりも少なく、
上記特性行列を同定する工程では、上記特性行列に基づいて算出されるモード特性のうち上記測定データに基づいて算出されたモード特性に対応するモード特性についてのみ、特性行列に基づいて算出されるモード特性の値が上記測定データに基づいて算出されたモード特性の値に一致するように上記特性行列の各成分の値を変更していくことで、上記特性行列を同定する、振動解析方法。
【請求項２】
解析対象物上に設定された複数の実測定点及び仮想測定点の座標データを要求する工程と、
上記解析対象物上の少なくとも一点において加振する加振試験が行われた時に実測定点においてのみ測定を行い、加振により得られた各実測定点における測定データに基づいてモード特性を算出する工程と、
上記座標データに基づいて上記解析対象物の振動を表す特性行列を構成する成分間の制約条件式を算出する工程と、
上記測定データに基づいて、算出されたモード特性の値に基づいて上記制約条件式を用いて上記特性行列を同定する工程とを具備する、振動解析方法において、
上記特性行列は実測定点及び仮想測定点の両測定点における自由度に等しい自由度で作成され、
上記特性行列を同定する工程では、上記特性行列に基づいて算出されるモード特性のうち上記測定データに基づいて算出されたモード特性に対応するモード特性についてのみ、特性行列に基づいて算出されるモード特性の値が上記測定データに基づいて算出されたモード特性の値に一致するように上記特性行列の各成分の値を変更していくことで、上記特性行列を同定する、振動解析方法。
【請求項３】
上記仮想測定点はその仮想測定点に最も近接する測定点が実測定点となるように設けられる、請求項２に記載の振動解析方法。
【請求項４】
測定点間の構造的結合状態の入力を要求する工程をさらに具備し、
上記制約条件式を算出する工程では、上記座標データに加えて測定点間の構造的結合状態に基づいて制約条件式を算出する、請求項２又は３に記載の振動解析方法。
【請求項５】
上記特性行列を同定する工程で同定される特性行列は質量行列であり、
上記特性行列を同定する工程において同定された質量行列に基づいて剛体特性を算出する工程をさらに具備する、請求項２〜４のいずれか１項に記載の振動解析方法。
【請求項６】
解析対象物上の測定点における位置、速度又は加速度の変位を測定する測定器と、上記解析対象物上の少なくとも一点において加振を行う加振器と、上記測定器及び加振器からそれぞれ測定データ及び加振データが入力されると共にこれらデータに基づいて演算処理を行う解析部とを具備する振動解析装置において、
上記解析部は、上記加振器での加振試験により得られた各測定点における測定データに基づいてモード特性を算出し、各測定点の座標データに基づいて上記解析対象物の振動を表す特性行列を構成する成分間の制約条件式を算出し、上記測定データに基づいて算出されたモード特性の値に基づいて上記制約条件式を用いて上記特性行列を同定し、
上記加振試験で上記測定点において実際に測定される測定自由度よりも上記特性行列の自由度の方が大きく、よって上記測定データに基づいて算出されるモード特性の値の数の方が上記特性行列に基づいて算出されるモード特性の値の数よりも少なく、
上記特性行列を同定する際には、上記特性行列に基づいて算出されるモード特性のうち上記測定データに基づいて算出されたモード特性に対応するモード特性についてのみ、特性行列に基づいて算出されるモード特性の値が上記測定データに基づいて算出されたモード特性の値に一致するように上記特性行列の各成分の値を変更していくことで、上記特性行列を同定する、振動解析装置。
【請求項７】
解析対象物上に設定された実測定点及び仮想測定点のうち実測定点における位置、速度又は加速度の変位を測定する測定器と、上記解析対象物上の少なくとも一点において加振を行う加振器と、上記測定器及び加振器からそれぞれ測定データ及び加振データが入力されると共にこれらデータに基づいて演算処理を行う解析部とを具備する振動解析装置において、
上記解析部は、上記加振器による加振試験により得られた実測定点における測定データに基づいてモード特性を算出し、上記実測定点及び仮想測定点の座標データに基づいて特性行列を構成する成分間の制約条件式を算出し、上記測定データに基づいて算出されたモード特性の値に基づいて上記制約条件式を用いて上記特性行列を同定し、
上記特性行列は実測定点及び仮想測定点の両測定点における自由度に等しい自由度で作成され、
上記特性行列を同定する際には、上記特性行列に基づいて算出されるモード特性のうち上記測定データに基づいて算出されたモード特性に対応するモード特性についてのみ、特性行列に基づいて算出されるモード特性の値が上記測定データに基づいて算出されたモード特性の値に一致するように上記特性行列の各成分の値を変更していくことで、上記特性行列を同定する、振動解析装置。
【請求項８】
上記仮想測定点はその仮想測定点に最も近接する測定点が実測定点となるように設けられる、請求項７に記載の振動解析装置。
【請求項９】
上記制約条件式を算出する際には、上記座標データに加えて測定点間の構造的結合状態に基づいて制約条件式が算出される、請求項７又は８に記載の振動解析装置。
【請求項１０】
上記同定される特性行列は質量行列であり、
上記解析部は、上記特性行列を同定する工程において同定された質量行列に基づいて剛体特性を算出する、請求項７〜９のいずれか１項に記載の振動解析装置。
【請求項１】
解析対象物上に設定された複数の測定点の座標データを要求する工程と、
上記解析対象物上の少なくとも一点において加振する加振試験が行われた時に得られた各測定点における測定データに基づいてモード特性を算出する工程と、
上記座標データに基づいて、上記解析対象物の振動を表す特性行列を構成する成分間の制約条件式を算出する工程と、
上記測定データに基づいて算出されたモード特性の値に基づいて上記制約条件式を用いて上記特性行列を同定する工程とを具備する、振動解析方法において、
上記加振試験で上記測定点において実際に測定される測定自由度よりも上記特性行列の自由度の方が大きく、よって上記測定データに基づいて算出されるモード特性の値の数の方が上記特性行列に基づいて算出されるモード特性の値の数よりも少なく、
上記特性行列を同定する工程では、上記特性行列に基づいて算出されるモード特性のうち上記測定データに基づいて算出されたモード特性に対応するモード特性についてのみ、特性行列に基づいて算出されるモード特性の値が上記測定データに基づいて算出されたモード特性の値に一致するように上記特性行列の各成分の値を変更していくことで、上記特性行列を同定する、振動解析方法。
【請求項２】
解析対象物上に設定された複数の実測定点及び仮想測定点の座標データを要求する工程と、
上記解析対象物上の少なくとも一点において加振する加振試験が行われた時に実測定点においてのみ測定を行い、加振により得られた各実測定点における測定データに基づいてモード特性を算出する工程と、
上記座標データに基づいて上記解析対象物の振動を表す特性行列を構成する成分間の制約条件式を算出する工程と、
上記測定データに基づいて、算出されたモード特性の値に基づいて上記制約条件式を用いて上記特性行列を同定する工程とを具備する、振動解析方法において、
上記特性行列は実測定点及び仮想測定点の両測定点における自由度に等しい自由度で作成され、
上記特性行列を同定する工程では、上記特性行列に基づいて算出されるモード特性のうち上記測定データに基づいて算出されたモード特性に対応するモード特性についてのみ、特性行列に基づいて算出されるモード特性の値が上記測定データに基づいて算出されたモード特性の値に一致するように上記特性行列の各成分の値を変更していくことで、上記特性行列を同定する、振動解析方法。
【請求項３】
上記仮想測定点はその仮想測定点に最も近接する測定点が実測定点となるように設けられる、請求項２に記載の振動解析方法。
【請求項４】
測定点間の構造的結合状態の入力を要求する工程をさらに具備し、
上記制約条件式を算出する工程では、上記座標データに加えて測定点間の構造的結合状態に基づいて制約条件式を算出する、請求項２又は３に記載の振動解析方法。
【請求項５】
上記特性行列を同定する工程で同定される特性行列は質量行列であり、
上記特性行列を同定する工程において同定された質量行列に基づいて剛体特性を算出する工程をさらに具備する、請求項２〜４のいずれか１項に記載の振動解析方法。
【請求項６】
解析対象物上の測定点における位置、速度又は加速度の変位を測定する測定器と、上記解析対象物上の少なくとも一点において加振を行う加振器と、上記測定器及び加振器からそれぞれ測定データ及び加振データが入力されると共にこれらデータに基づいて演算処理を行う解析部とを具備する振動解析装置において、
上記解析部は、上記加振器での加振試験により得られた各測定点における測定データに基づいてモード特性を算出し、各測定点の座標データに基づいて上記解析対象物の振動を表す特性行列を構成する成分間の制約条件式を算出し、上記測定データに基づいて算出されたモード特性の値に基づいて上記制約条件式を用いて上記特性行列を同定し、
上記加振試験で上記測定点において実際に測定される測定自由度よりも上記特性行列の自由度の方が大きく、よって上記測定データに基づいて算出されるモード特性の値の数の方が上記特性行列に基づいて算出されるモード特性の値の数よりも少なく、
上記特性行列を同定する際には、上記特性行列に基づいて算出されるモード特性のうち上記測定データに基づいて算出されたモード特性に対応するモード特性についてのみ、特性行列に基づいて算出されるモード特性の値が上記測定データに基づいて算出されたモード特性の値に一致するように上記特性行列の各成分の値を変更していくことで、上記特性行列を同定する、振動解析装置。
【請求項７】
解析対象物上に設定された実測定点及び仮想測定点のうち実測定点における位置、速度又は加速度の変位を測定する測定器と、上記解析対象物上の少なくとも一点において加振を行う加振器と、上記測定器及び加振器からそれぞれ測定データ及び加振データが入力されると共にこれらデータに基づいて演算処理を行う解析部とを具備する振動解析装置において、
上記解析部は、上記加振器による加振試験により得られた実測定点における測定データに基づいてモード特性を算出し、上記実測定点及び仮想測定点の座標データに基づいて特性行列を構成する成分間の制約条件式を算出し、上記測定データに基づいて算出されたモード特性の値に基づいて上記制約条件式を用いて上記特性行列を同定し、
上記特性行列は実測定点及び仮想測定点の両測定点における自由度に等しい自由度で作成され、
上記特性行列を同定する際には、上記特性行列に基づいて算出されるモード特性のうち上記測定データに基づいて算出されたモード特性に対応するモード特性についてのみ、特性行列に基づいて算出されるモード特性の値が上記測定データに基づいて算出されたモード特性の値に一致するように上記特性行列の各成分の値を変更していくことで、上記特性行列を同定する、振動解析装置。
【請求項８】
上記仮想測定点はその仮想測定点に最も近接する測定点が実測定点となるように設けられる、請求項７に記載の振動解析装置。
【請求項９】
上記制約条件式を算出する際には、上記座標データに加えて測定点間の構造的結合状態に基づいて制約条件式が算出される、請求項７又は８に記載の振動解析装置。
【請求項１０】
上記同定される特性行列は質量行列であり、
上記解析部は、上記特性行列を同定する工程において同定された質量行列に基づいて剛体特性を算出する、請求項７〜９のいずれか１項に記載の振動解析装置。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【公開番号】特開２００９−２０４５１７（Ｐ２００９−２０４５１７Ａ）
【公開日】平成２１年９月１０日（２００９．９．１０）
【国際特許分類】
【出願番号】特願２００８−４８２７７（Ｐ２００８−４８２７７）
【出願日】平成２０年２月２８日（２００８．２．２８）
【出願人】（３０４０２１４１７）国立大学法人東京工業大学 (1,821)
【Ｆターム（参考）】
【公開日】平成２１年９月１０日（２００９．９．１０）
【国際特許分類】
【出願日】平成２０年２月２８日（２００８．２．２８）
【出願人】（３０４０２１４１７）国立大学法人東京工業大学 (1,821)
【Ｆターム（参考）】
[ Back to top ]