(57)【要約】【課題】 本発明は、単点加振実験的特性行列同定法の方法、装置およびコンピュータ読み取り可能な記録媒体に関する。【解決手段】 解析対象物の測定点の座標データを求める処理と、解析対象物を単点加振して周波数応答関数を求める処理と、任意の弾性構造物モデルの質量行列が剛性質量行列に変換できること、剛性質量行列が或る特定の要素構成となること、および、剛体運動状態では解析対象物がそのすべての領域でひずみと応力とを生じないこと、という力学原理から特性行列成分間の制約条件式を求める処理と、同定周波数帯域内の共振周波数を設定しモード特性を同定する処理と、モード特性の同定から得られる不減衰固有振動数および固有モードを目標値として質量行列および剛性行列を求める処理を含むプログラムを記憶媒体であるＦＤ１０１から記憶装置１２９に格納して、中央処理装置１２８で実行するようにする。【選択図】        図１
2007147634.tif 000002
【特許請求の範囲】
【請求項１】
解析対象物の単点加振によって前記解析対象物の複数測定点に生じる変位の測定データ（Δｘｉ，Δｙｉ，Δｚｉ）の入力を受け付け（ただしｉ＝１・・ｎ）、前記測定データ（Δｘｉ，Δｙｉ，Δｚｉ）と重力加速度ｇとを行列要素として、（Ａ１）式の行列［Ψ］を作成するステップと、
【数１】
2007147634.tif 000020
既知の前記行列［Ψ］と、既知の行列要素を予め含む（Ａ２）式の剛体質量行列［Ｍrigid］とを、弾性構造物モデルの未知の質量行列［Ｍ］を剛体質量行列［Ｍrigid］に変換する（Ａ３）式に代入することにより、前記質量行列［Ｍ］の未知の行列要素が満足すべき連立一次の第１制約条件式を得るステップと、
【数２】
2007147634.tif 000021
【数３】
2007147634.tif 000022
剛体運動状態では前記解析対象物にひずみと応力が生じないことを示す（Ａ４）式に、既知の前記行列［Ψ］を代入することにより、剛性行列［Ｋ］の未知の行列要素が満足すべき連立一次の第２制約条件式を得るステップと、
［Ｋ］［Ψ］＝［０］  ・・・（Ａ４）式
前記解析対象物の単点加振によって前記複数測定点に生じる振動の周波数応答から求められるモード特性の解析データ（不減衰固有振動数Ω，固有モードのベクトル｛φ｝）の入力を受け付けるステップと、
前記第１制約条件式を満足する前記質量行列［Ｍ］の候補と、前記第２制約条件式を満足する前記剛性行列［Ｋ］の候補との中から、（Ａ５）式を満足するものを探索することにより、前記解析対象物の前記質量行列［Ｍ］と前記剛性行列［Ｋ］を特定するステップと
（［Ｋ］−Ω2［Ｍ］）｛φ｝＝｛０｝    ・・・（Ａ５）式
をコンピュータで実行することを特徴とする振動解析方法。
【請求項２】
請求項１に記載の振動解析方法において、
前記（Ａ２）式中の前記補正係数αは、
１．０≦補正係数α≦１．３
であること
を特徴とする振動解析方法。
【請求項３】
請求項１または請求項２に記載の振動解析方法において、
前記質量行列［Ｍ］と前記剛性行列［Ｋ］が予め定められる許容範囲内にないと判断すると、固有モード行列［Φ］、不減衰固有振動数行列［Ω］、および単位行列［Ｅ］を下式に代入し、
［Ｍ］＝［Φ］-T［Ｅ］［Φ］-1
［Ｋ］＝［Φ］-T［Ω2］［Φ］-1
前記質量行列［Ｍ］および前記剛性行列［Ｋ］を算出するステップ
を更にコンピュータで実行することを特徴とする振動解析方法。
【請求項４】
請求項１または請求項２に記載の振動解析方法において、
前記質量行列［Ｍ］および前記剛性行列［Ｋ］が予め定められる許容範囲内にないと判断すると、モード次数ｓの固有モード｛φｓ｝と不減衰固有振動数Ωｓについて（Ａ６）式の第３制約条件式を満足するか否かを判断し、
【数４】
2007147634.tif 000023
前記第３制約条件式を満足しない剰余固有振動数を同定対象の周波数帯域から高い周波数へ移動するステップ
を更にコンピュータで実行することを特徴とする振動解析方法。
【請求項５】
請求項１ないし請求項４のいずれか１項に記載の振動解析方法において、
前記解析対象物の単点加振によって前記複数測定点に生じる振動の周波数応答関数Ｈ（ω）の解析結果の入力を受け付けるステップと、
前記質量行列［Ｍ］および前記剛性行列［Ｋ］から算出される周波数応答関数ｈ（ω）が、前記周波数応答関数Ｈ（ω）と最も一致するように、前記質量行列［Ｍ］および前記剛性行列［Ｋ］を定数倍して最適化するステップと
を更にコンピュータで実行することを特徴とする振動解析方法。
【請求項６】
解析対象物の単点加振によって前記解析対象物の複数測定点に生じる変位の測定データ（Δｘｉ，Δｙｉ，Δｚｉ）の入力を受け付け（ただしｉ＝１・・ｎ）、前記測定データ（Δｘｉ，Δｙｉ，Δｚｉ）と重力加速度ｇとを行列要素として、（Ａ１）式の行列［Ψ］を作成する第１処理手段と、
【数５】
2007147634.tif 000024
既知の前記行列［Ψ］と、既知の行列要素を予め含む（Ａ２）式の剛体質量行列［Ｍrigid］とを、弾性構造物モデルの未知の質量行列［Ｍ］を剛体質量行列［Ｍrigid］に変換する（Ａ３）式に代入することにより、前記質量行列［Ｍ］の未知の行列要素が満足すべき連立一次の第１制約条件式を得る第２処理手段と、
【数６】
2007147634.tif 000025
【数７】
2007147634.tif 000026
剛体運動状態では前記解析対象物にひずみと応力が生じないことを示す（Ａ４）式に、既知の前記行列［Ψ］を代入することにより、剛性行列［Ｋ］の未知の行列要素が満足すべき連立一次の第２制約条件式を得る第３処理手段と、
［Ｋ］［Ψ］＝［０］  ・・・（Ａ４）式
前記解析対象物の単点加振によって前記複数測定点に生じる振動の周波数応答から求められるモード特性の解析データ（不減衰固有振動数Ω，固有モードのベクトル｛φ｝）の入力を受け付ける第４処理手段と、
前記第１制約条件式を満足する前記質量行列［Ｍ］の候補と、前記第２制約条件式を満足する前記剛性行列［Ｋ］の候補との中から、（Ａ５）式を満足するものを探索することにより、前記解析対象物の前記質量行列［Ｍ］と前記剛性行列［Ｋ］を特定する第５処理手段と
（［Ｋ］−Ω2［Ｍ］）｛φ｝＝｛０｝    ・・・（Ａ５）式
を備えたことを特徴とする振動解析装置。
【請求項７】
請求項６に記載の振動解析装置において、
前記（Ａ２）式中の前記補正係数αは、
１．０≦補正係数α≦１．３
であること
を特徴とする振動解析装置。
【請求項８】
請求項６または請求項７に記載の振動解析装置において、
前記質量行列［Ｍ］と前記剛性行列［Ｋ］が予め定められる許容範囲内にないと判断すると、固有モード行列［Φ］、不減衰固有振動数行列［Ω］、および単位行列［Ｅ］を下式に代入し、
［Ｍ］＝［Φ］-T［Ｅ］［Φ］-1
［Ｋ］＝［Φ］-T［Ω2］［Φ］-1
前記質量行列［Ｍ］および前記剛性行列［Ｋ］を算出する第６処理手段
を更に備えたことを特徴とする振動解析装置。
【請求項９】
請求項６または請求項７に記載の振動解析装置において、
前記質量行列［Ｍ］および前記剛性行列［Ｋ］が予め定められる許容範囲内にないと判断すると、モード次数ｓの固有モード｛φｓ｝と不減衰固有振動数Ωｓについて（Ａ６）式の第３制約条件式を満足するか否かを判断し、
【数８】
2007147634.tif 000027
前記第３制約条件式を満足しない剰余固有振動数を同定対象の周波数帯域から高い周波数へ移動する第７処理手段
を更に備えたことを特徴とする振動解析装置。
【請求項１０】
請求項６ないし請求項９のいずれか１項に記載の振動解析装置において、
前記解析対象物の単点加振によって前記複数測定点に生じる振動の周波数応答関数Ｈ（ω）の解析結果の入力を受け付ける第８処理手段と、
前記質量行列［Ｍ］および前記剛性行列［Ｋ］から算出される周波数応答関数ｈ（ω）が、前記周波数応答関数Ｈ（ω）と最も一致するように、前記質量行列［Ｍ］および前記剛性行列［Ｋ］を定数倍して最適化する第９処理手段と
を更にコンピュータで実行することを特徴とする振動解析装置。
【請求項１１】
コンピュータに、請求項１ないし請求項５のいずれか１項に記載の振動解析方法を実行させるための振動解析プログラムを記録したコンピュータ読取り可能な記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、振動実験結果から直接特性行列を求める実験的特性行列同定法に関し、構造物に単一の加振を与えて得られる多点周波数応答関数を用いて構造物の振動を解析する方法および装置に関する。さらに、この方法を実行させるためのプログラムを記録したコンピュータ読み取り可能な記録媒体に関する。
【背景技術】
【０００２】
各種の機械や構造物などの振動に対して、何らかの対策を施す必要があるとき、まず、振動特性を何らかの方法で把握し、動特性を表現する数式モデルを決定する。一般に、この数式モデルを決定することを「系を同定する」と呼ぶが、この同定法には、理論的同定法と実験的同定法とがある。
【０００３】
理論的同定法では、有限要素法（FEM ）を用いて構造物を離散化し、特性行列を決定する方法が主流である。一方、実験的同定法では、解析対象の構造物に加振実験を行い周波数応答関数（伝達関数）を得てモード特性を同定するモード特性同定法が主流である。さらに、このモード特性同定法には、１自由度法、多点自由度法および多点参照法がある。１自由度法は、１つの周波数応答関数に含まれる１つの固有モードのみに着目して、これを同定する方法である。多自由度法は、異なる固有モードどうしの影響を考慮しながら複数の固有モードのモード特性を同時に同定する方法である。多点参照法は、複数の周波数応答関数を同時に参照することにより全体項（固有振動数、モード減衰比）の精度を高めるとともに一貫したモード特性を得る方法である。
【非特許文献１】大熊政明「実験的特性行列同定法の開発（開発理論と基礎的検証）」日本機械学会論文集Ｃ VOL.63 NO.616 P.4171〜P.4178（1997/12）
【発明の開示】
【発明が解決しようとする課題】
【０００４】
ところで、理論的同定法では、離散化を行っているため、複雑な構造物に対しては精度のよいモデル化が難しいという問題がある。
一方、実験的同定法では、低次の固有振動数から必要なモードのみを取り出して自由度を縮小する近似手法であるため、全自由度の特性行列を構築することが困難であるという問題がある。さらに、実験的同定法では、６自由度のうち３自由度は物理法則により他の３自由度に従属させることができるとしても、なお３自由度は独立であるため、少なくとも場所の異なる３点について加振実験を行う必要があるという問題がある。
【０００５】
そこで、本発明では、構造物に単一の加振を与えて得られる多点周波数応答関数から直接特性行列を求めることができる、振動解析の方法および装置を提供することを目的とする。さらに、この方法を実行させるためのプログラムを記録したコンピュータ読み取り可能な記録媒体を提供することを目的とする。
【課題を解決するための手段】
【０００６】
《１》  本発明の振動解析方法は、コンピュータで次のステップを実行する。
(1)  解析対象物の単点加振によって解析対象物の複数測定点に生じる変位の測定データ（Δｘｉ，Δｙｉ，Δｚｉ）の入力を受け付け（添え字ｉは測定点をあらわす。ｉ＝１・・ｎ）、測定データ（Δｘｉ，Δｙｉ，Δｚｉ）と重力加速度ｇとを行列要素として、（Ａ１）式の行列［Ψ］を作成するステップ
(2)  既知の行列［Ψ］と、既知の行列要素を予め含む（Ａ２）式の剛体質量行列［Ｍrigid］とを、弾性構造物モデルの未知の質量行列［Ｍ］を剛体質量行列［Ｍrigid］に変換する（Ａ３）式に代入することにより、質量行列［Ｍ］の未知の行列要素が満足すべき連立一次の第１制約条件式を得るステップ
(3)  剛体運動状態では解析対象物にひずみと応力が生じないことを示す（Ａ４）式に、既知の行列［Ψ］を代入することにより、剛性行列［Ｋ］の未知の行列要素が満足すべき連立一次の第２制約条件式を得るステップ
［Ｋ］［Ψ］＝［０］  ・・・（Ａ４）式
(4)  解析対象物の単点加振によって複数測定点に生じる振動の周波数応答から求められるモード特性の解析データ（不減衰固有振動数Ω，固有モードのベクトル｛φ｝）の入力を受け付けるステップ
(5)  第１制約条件式を満足する質量行列［Ｍ］の候補と、第２制約条件式を満足する剛性行列［Ｋ］の候補との中から、（Ａ５）式を満足するものを探索することにより、解析対象物の質量行列［Ｍ］と剛性行列［Ｋ］を特定するステップ
（［Ｋ］−Ω2［Ｍ］）｛φ｝＝｛０｝    ・・・（Ａ５）式
《２》  なお好ましくは、（Ａ２）式中の補正係数αは、１．０≦補正係数α≦１．３である。
《３》  また好ましくは、コンピュータは下記ステップを更に実行する。
質量行列［Ｍ］と剛性行列［Ｋ］が予め定められる許容範囲内にないと判断すると、固有モード行列［Φ］、不減衰固有振動数行列［Ω］、および単位行列［Ｅ］を下式に代入して、質量行列［Ｍ］および剛性行列［Ｋ］を改めて算出するステップ。
［Ｍ］＝［Φ］-T［Ｅ］［Φ］-1
［Ｋ］＝［Φ］-T［Ω2］［Φ］-1
《４》  なお好ましくは、コンピュータは下記ステップを更に実行する。
質量行列［Ｍ］および剛性行列［Ｋ］が予め定められる許容範囲内にないと判断すると、モード次数ｓの固有モード｛φｓ｝と不減衰固有振動数Ωｓについて（Ａ６）式の第３制約条件式を満足するか否かを判断し、この第３制約条件式を満足しない剰余固有振動数を同定対象の周波数帯域から高い周波数へ移動するステップ。
《５》  また好ましくは、コンピュータは下記ステップを更に実行する。
(1)  解析対象物の単点加振によって複数測定点に生じる振動の周波数応答関数Ｈ（ω）の解析結果の入力を受け付けるステップ
(2)  質量行列［Ｍ］および剛性行列［Ｋ］から算出される周波数応答関数ｈ（ω）が、周波数応答関数Ｈ（ω）と最も一致するように、質量行列［Ｍ］および剛性行列［Ｋ］を定数倍して最適化するステップ
《６》  本発明の振動解析装置は、上記《１》に対応する装置発明である。
《７》  本発明の振動解析装置は、上記《２》に対応する装置発明である。
《８》  本発明の振動解析装置は、上記《３》に対応する装置発明である。
《９》  本発明の振動解析装置は、上記《４》に対応する装置発明である。
《１０》  本発明の振動解析装置は、上記《５》に対応する装置発明である。
《１１》  本発明のコンピュータ読取可能な記録媒体には、上記《１》〜《５》のいずれか１項に記載の振動解析方法を実行させるための振動解析プログラムを記録する。
【発明の効果】
【０００７】
本発明によれば、構造物に単一の加振を与えることで、多点の周波数応答関数から直接的に特性行列を求めることができる。そして、本発明では、実験的に特性行列の同定を行うことができるので、直接振動解析を行いたい構造物の物理座標モデルを表現することができる。このため、本発明は、構造物のいかなる動挙動も容易に解析することができる。そして、本発明は、構造物を複数の部分構造に分け各々の動特性を求めそれらを合成して構造物全系の解析を行う部分構造合成法、動特性を設計者の望む値に変更する最適設計、有限要素法との一体化、および、各種シミュレーションなどに直接的に適用することができる。
【０００８】
もちろん、本発明は、多点加振して得られたデータを利用して目的の解析処理をすることもできる。
【発明を実施するための最良の形態】
【０００９】
以下、本発明の実施形態について図面に基づいて説明する。
【００１０】
（実施形態の構成）
図１は、本実施形態における解析装置のブロック図である。
図１において、解析装置は、測定部と解析部とに大別される。測定部は、打撃ハンマ１１１、力変換器１１２、加速度計１１３、増幅器１１４、１１５および外部インタフェイス部１２１を備えて構成される。解析部は、フレキシブルディスク装置（以下、「ＦＤ装置」と略記する。）１２２、フレキシブルディスク制御部（以下、「ＦＤ制御部」と略記する。）１２３、入力装置１２４、入力制御部１２５、出力装置１２６、出力制御部１２７、バス１３０、中央処理装置１２８および記憶装置１２９を備えて構成される。
【００１１】
打撃ハンマ１１１は、解析対象物１４０に加振を与え、その加えられた力の大きさは、打撃ハンマ１１１に設けられた力変換器１１２によって電気信号に変換される。この電気信号は、増幅器１１４で増幅された後に外部インタフェース１２１に出力される。
【００１２】
加速度計１１３は、変換器として用いられ、ピエゾ素子によって、解析対象物１４０に発生した加速度の大きさを測定してこれを電気信号に変換する。この電気信号は、増幅器１１５で増幅された後に外部インタフェース１２１に出力される。解析者は、注目したい次数までの固有モードを幾何学的に充分よく表現することができる程度に適切な分布と数で測定点を決定し、加速度計１１３は、その決定箇所に取り付けられる。
【００１３】
外部インタフェース１２１は、入力された電気信号を標本化してディジタル信号に変換し、解析部に合った信号レベルに調整して、バス１３０に信号を出力する。
ＦＤ装置１２２は、ＦＤ制御部１２３によって制御され、記録媒体の１つであるフレキシブルディスク（以下、「ＦＤ」と略記する。）１０１に記録された本発明にかかるプログラムを解析装置１００にインストールしたり、解析結果をこのＦＤ１０１または別のＦＤに記録したりする。
【００１４】
入力装置１２４は、キーボード、マウスおよびデジタイザなどの入力機器であり、入力制御部１２５によって制御される。解析者は、入力装置１２５を操作してコマンドや解析対象の形状などの各種データを入力する。
出力装置１２６は、ＣＲＴなどのディスプレイやプリンタなどの出力機器であり、出力制御部１２７によって制御される。出力装置１２６は、入力装置１２４に入力された各種データやコマンド、プログラム実行のためのメニューおよび解析結果などが出力される。
【００１５】
記憶装置１２９は、半導体メモリやハードディスなどであり、本発明にかかるプログラムおよび各プログラム実行中の一時的なデータや解析結果などが格納される。
【００１６】
中央処理装置１２８は、バス１３０を介して、外部インタフェース１２１、ＦＤ制御部１２３、入力制御部１２５、出力制御部１２７および記憶装置１２９と接続され、これらに指示を与えたり、記憶装置１２９に格納され各プログラムに従いフーリエ変換や実験的特性行列同定法などの処理を行う。
【００１７】
（本実施形態の作用効果）
次に、本実施形態の作用効果を説明する。
図２ないし図５は、本実施形態における解析装置のフローチャートを示す図である。
【００１８】
解析者は、解析したい解析対象物１４０を用意する。解析対象物１４０は、ゴムひもやタイヤチューブなどの弾性体でつるされ、周辺自由状態が近似的に実現される。なお、弾性支持状態でもよい。
図２ないし図５において、解析者は、解析装置を入力装置１２４から起動する。中央処理装置１２８は、記憶装置１２９に格納されている実験的特性行列同定法のプログラムを読み込み、実行する（Ｓ１）。
【００１９】
次に、中央処理装置１２８は、解析対象物１４０の形状・寸法などの初期データを入力するように解析者に出力画面１２６を介して要求する。解析者は、入力装置１２４から、解析対象物の測定点座標データを入力する（Ｓ２）。
【００２０】
次に、解析者は、打撃ハンマ１１１によって加速度計１１３が取り付けられた位置の中のいずれかの１つの位置（単点）で解析対象物１４０に打撃を加える。
解析者が打撃ハンマ１１１によって解析対象物１４０に加えた打撃は、力変換器１１２によって測定され、増幅器１１５および外部インタフェイス部１２１を介して、中央処理装置１２８に取り込まれる。そして、この打撃によって解析対象物１４０に発生した振動も、加速度計１１３によって測定され、増幅器１１５および外部インタフェイス部１２１を介して、中央処理装置１２８に取り込まれる。このようにして単点加振試験が行われる（Ｓ３）。
【００２１】
次に、中央処理装置１２８は、これらデータをフーリエ変換（例えば、高速フーリエ変換）して周波数スペクトルになおす。中央処理装置１２８は、入力（力変換器１１２の出力）と出力（加速度計１１３の出力）からパワースペクトルとクロススペクトルとを計算し、周波数応答関数およびコヒーレンス関数を計算する。
【００２２】
ここで、コヒーレンス関数（関連度関数）は、入力と出力との内積であり、入出力の相関を示すものである。その値が０である場合には、入力と出力とは無相関であり、１に近づくにつれて相関性が増す。コヒーレンス関数の値が測定データの信頼性を示すので、後述される重み付き最小自乗法の重み関数は、このコヒーレンス関数の値を使用する。
【００２３】
なお、より信頼性の高いデータを得る観点から、この単点加振試験を複数回繰り返し、個々のパワースペクトルおよびクロススペクトルを平均し、その平均値を周波数応答関数とコヒーレンス関数とを求めるために用いられるパワースペクトルおよびクロススペクトルとすることが好ましい。
また、予めこれらのデータを記憶装置１２９またはＦＤ１０１に格納しておき、解析装置は、起動に従って出力装置１２６に格納内容を表示して、解析者は、これから選択して座標データ、単点加振試験のデータを入力するようにしてもよい。
【００２４】
このＳ２の処理ないしＳ４の処理が振動測定の処理である。
次に、測定点間の構造的結合の定義・設定を行う。すなわち、測定点間においてどの測定点がどの測定点と構造的に結びついているかを解析者が入力装置１２４から入力することにより、測定点相互の物理的結合状態をモデル化する（Ｓ５）。
【００２５】
次に、中央処理装置１２８は、式１、式３および式４からなる特性行列の成分間の制約条件式を計算する（Ｓ６）。これによって解析対象物１４０の物理的制約条件が決定され、第１次物理モデル化の処理が行われる。
【数９】
2007147634.tif 000003
【数１０】
2007147634.tif 000004
【数１１】
2007147634.tif 000005
【数１２】
2007147634.tif 000006
これら式１ないし式４は、本発明にかかる単点加振実験的特性行列同定法の基礎方程式である。これら基礎方程式について説明する。なお、説明の都合上、座標系は、ｘｙｚ直角座標系を使用するが、どのような座標系においても同様に解析することができる。
【００２６】
同定される特性行列の質量行列および剛性行列の自由度は、（測定点×３）であるが、質量行列は、『任意の弾性構造物モデルの質量行列が剛体質量行列に変換できる』という力学原理から式１となる。
ここで、同定される質量行列を［Ｍ］と表す。［Ψ］は、式２に示すようにｎ行６列の行列であり、測定点の座標値からｘ方向単位並進変位、ｙ方向単位並進変位、ｚ方向単位方向並進変位、ｘ軸回り単位回転変位、ｙ軸回り単位回転変位およびｚ軸回り単位回転変位の剛体変位ベクトルを順に並べた行列である。これは、最も一般的な場合の［Ψ］の構成を考えて剛体自由度が６であることから、互いに直交する方向への並進変位モード３つとそれらに対する回転変位モード３つとから構成された行列である。
【００２７】
なお、△ｘ、△ｙおよび△ｚは、正確にはこれに方向余弦を乗じたものであるが、方向余弦を「１」と近似している。また、ｇは、重力加速度であり、△ｘ、△ｙおよび△ｚに添えられている数字は、測定点の番号に対応する。そして、［Ψ］T は［Ψ］の転置行列である。転置行列は、以下、行列に上付き１／４角の「Ｔ」を添えて表示する。
【００２８】
この［Ψ］を剛体変位モード行列と呼ぶことにする。
一方、式１の［Ｍrigid］ は、『剛体質量行列が或る特定の要素構成となる』という力学原理から式３となる。
ここで、ｗは、解析対象物１４０の重さ、Ｉxx、Ｉyy、Ｉzzは、それぞれ座標原点を通る各軸回りの慣性モーメント、Ｉyx、Ｉzx、Ｉyzは、慣性乗積である。Ａ、Ｂ、Ｃは、剛体表現の点と重心とのズレに起因する値であり、解析対象物１４０の質量をｍ、重心座標を（ｘg、ｙg、ｚg ）とすると、Ａ＝ｍｘg 、Ｂ＝ｍｙg 、Ｃ＝ｍｚg である。そして、αは同定精度を上げるための補正係数である。［Ｍrigid］ は、式３に示すように対称行列である。
【００２９】
なお、［Ψ］の重力加速度で表した成分を「１」とすれば、ｗは、解析対象物１４０の質量ｍとなる。
［Ｍrigid］ の成分は、質量と重心位置が未定であっても式３に示す関係が成立している。そこで、式１の関係から、未知である質量行列［Ｍ］の行列成分に関して連立１次方程式が構成され、この方程式の数に等しい未知成分を他の成分に従属させることができ、これを質量行列［Ｍ］の物理的な制約条件式として利用することができる。
【００３０】
ここで、解析対象物１４０の質量（重さ）など剛体特性値が既知として与えられていれば、これを式３に代入することにより行列成分間の制約条件式を生成するようにしてもよい。このように既知の剛体特性値が得られれば、この既知の剛体特性値については、他の成分で表す必要が無くなるという利点がある。
【００３１】
さらに、同定すべき剛性行列［Ｋ］については、『剛体運動状態では前記解析対象物がそのすべての領域でひずみと応力とを生じない』という力学原理から式４となる。
よって、この式４から剛性行列［Ｋ］の成分間の物理的制約条件式が導出される。つまり、剛性行列［Ｋ］の成分の中のいくつかを独立変数とし、他の成分は、その独立変数とした成分に従属させることができる。
【００３２】
この結果、単点加振実験的特性行列同定法は、質量行列［Ｍ］および剛性行列［Ｋ］における独立変数を、実験的に得られた固有モードおよび不減衰固有振動数に、適切に一致させる一種の最適化問題となる。この一致させるべき、固有モードを目標固有モードと呼び、不減衰固有振動数を目標不減衰固有振動数と呼ぶことにする。さらに、減衰行列を最適化する場合に用いられる、一致させるべきモード減衰比を目標モード減衰比と呼ぶことにする。
【００３３】
図２ないし図５に戻って、この最適化について説明する。
まず、中央処理装置１２８は、同定する周波数帯域の入力を指示し、解析者は、入力装置１２４から同定する周波数帯域を入力する（Ｓ７）。同定周波数帯域は、一次共振周波数より低い値から要求する或る周波数までとする。
次に、中央処理装置１２８は、所定のモード同定法に従う処理を行い、第１次共振周波数を含む同定周波数帯域の周波数応答関数を利用して、慣性項パラメータ、同定周波数範囲内の共振に関するモード特性および剰余項パラメータを同定する。この周波数応答関数には、測定で得られた結果を用いる。
【００３４】
これにより、中央処理装置１２８は、目標固有モード、目標不減衰固有振動数および目標モード減衰比を設定する（Ｓ８）。
モード特性同定法としては、一般に知られている各種のモード同定法を利用することができるが、例えば、多点参照法（ポリリファレンス法、poly-reference method）、偏分反復法、プロリー法（Prony's method）、サークルフィット法（circle fit method）などがある。
【００３５】
このＳ７およびＳ８の処理が目標モード特性の設定の処理である。
次に、中央処理装置１２８は、質量行列［Ｍ］および剛性行列［Ｋ］のそれぞれについて、独立未知数に乱数を当てはめ、初期値の行列を作成する（Ｓ９）。
次に、中央処理装置１２８は、式５により質量行列［Ｍ］の正定値化を行う（Ｓ１０）。すなわち、質量行列［Ｍ］に関し、すべてのｎ個の固有値が正値となるように特性行列の独立未知数を修正する。さらに、中央処理装置１２８は、式６により剛性行列［Ｋ］の正定値化を行う（Ｓ１１）。すなわち、剛性行列［Ｋ］に関し、既に物理的制約条件式から剛体変位自由度に相当する数の固有値は、必然的に零の値となるように設定されているが、他の固有値がすべて正値となるように特性行列の独立未知数を修正する。
【数１３】
2007147634.tif 000007
【数１４】
2007147634.tif 000008
Ｓ９において、質量行列［Ｍ］および剛性行列［Ｋ］が単に乱数で初期値を与えられるだけなので、このＳ１０およびＳ１１によって、物理的に質量行列［Ｍ］が正定値行列でなければならない条件や剛性行列［Ｋ］が剛体変位自由度だけ固有値零が存在する非負定値行列でなければならない条件を満足させる。
【００３６】
なお、剛性行列［Ｋ］は、周辺自由境界条件では非負定値行列であるが、弾性拘束境界条件では正定値行列である。
次に、中央処理装置１２８は、式７で表される不減衰系について、一般固有値問題を解き、不減衰固有振動数および固有モードを求める（Ｓ１２、Ｓ１３、Ｓ１４）。
【数１５】
2007147634.tif 000009
ここで、Ωは不減衰固有振動数であり、｛φ｝は固有モードのベクトルである。なお、本明細書ではベクトルを｛｝で示す。
より具体的に説明すると、中央処理装置１２８は、まず、Ｓ８で求めた目標不減衰固有振動数と対応するように、式７を用いて、低次側から同定周波数帯以内に存在する数の固有振動数Ωが満足な一致を示すまで反復の感度解析により質量行列［Ｍ］および剛性行列［Ｋ］の修正を行う（Ｓ１２）。
【００３７】
次に、中央処理装置１２８は、Ｓ８で求めた目標固有モードと対応するように、式７を用いて、固有モード｛φ｝に対しても同様に反復の感度解析により質量行列［Ｍ］および剛性行列［Ｋ］の修正を行う（Ｓ１３）。
ここで、Ｓ１０ないしＳ１３の処理において、処理すべき特性行列について、正規化処理を行う。これは、一般固有値問題を解くことから何らかの基準を必要とするからである。例えば、処理すべき特性行列の成分における最大絶対値を「１」とする正規化処理を行えばよい。
【００３８】
また、この反復の感度解析としては、一般に知られている各種の感度解析の手法を利用することができるが、例えば、最急降下法、共役勾配法などニュートン法に代表されるニュートン法類の手法を利用することができる。
次に、中央処理装置１２８は、Ｓ１２およびＳ１３で質量行列［Ｍ］および剛性行列［Ｋ］の修正を行ったので、これらが正定値行列であるか否かを判定し、正定値行列である場合には、Ｓ１５の処理を行い、正定値行列でない場合には、Ｓ１０に戻り、上述のＳ１０ないしＳ１３の処理を行う。
【００３９】
このＳ９ないしＳ１４の処理が第１段階の同定の処理である。
この第１段階の同定の処理によって、質量行列［Ｍ］および剛性行列［Ｋ］を求められるが、これらをより高精度に同定する観点からさらに第２段階の同定の処理を行うことが好ましい。次に、この第２段階の同定の処理、Ｓ１５ないしＳ１９について説明する。
【００４０】
中央処理装置１２８は、第１段階の同定処理で求められた質量行列［Ｍ］および剛性行列［Ｋ］が許容範囲内で一致しているか否かを判断し、一致している場合には、第２段階の同定を行うことなくＳ２０の処理を行い、一致していない場合には、Ｓ１６の処理を行う。このＳ１５の判断は、実験結果から得られた特性行列と計算によって得られた特性行列との相関性を示す値を計算し、この値が許容範囲を満足するか否かで判断する。許容範囲は、一致の精度を左右する。例えば、内積を計算しその値が０．９以上となるか否かで判断する。より精度良く一致させるためには、０．９より１に近い値、例えば、０．９５に設定すればよい。
【００４１】
次に、中央処理装置１２８は、第１次物理モデル化条件の遵守の有無を入力するように解析者に出力画面１２６を介して要求する。中央処理装置１２８は、入力された指示に基づき第１次物理モデル化条件を遵守する必要があるか否かを判断し、必要がない場合にはＳ１７の処理を行い、必要がある場合にはＳ１８の処理を行う（Ｓ１６）。
【００４２】
なお、Ｓ１５における許容値を０．９より小さい値とした場合でも強制的にＳ１６以下の処理を行うようにすれば、精度よく一致させることが可能である。
Ｓ１７の処理では、中央処理装置１２８は、固有モード｛φ｝における個々の固有モードについて、実験結果から得られた固有モードと計算によって得られた固有モードとの相関性を計算する。そして、所定の許容範囲外である固有モードを目標固有モードに置き換える。この置き換えられた固有モードベクトルを｛Φ｝とおき、式８および式９により、質量行列［Ｍ］および剛性行列［Ｋ］を修正する（Ｓ１７）。
【数１６】
2007147634.tif 000010
【数１７】
2007147634.tif 000011
ここで、［Ｅ］は単位行列であり、［Φ］-1 は［Φ］の逆行列である。逆行列は、以下、行列に上付き１／４角の「−１」を添えて表示する。
一方、Ｓ１８の処理では、中央処理装置１２８は、式１０に基づき、計算された剛性モードと目標固有モードとによる特性行列成分間における制約条件式を求める（Ｓ１８）。
【数１８】
2007147634.tif 000012
ここで、φおよびΩに添えられた下付１／４角の「ｓ」はモードの次数を示す。
次に、中央処理装置１２８は、すべての剰余固有振動数を同定周波数帯域より高い周波数帯域に移動させる（Ｓ１９）。より具体的に以下に説明する。
【００４３】
剛性行列中の或る独立変数ｋij（i行j列成分）を△ｋijだけ微少増加させた場合の第ｓ次の固有振動数における変化量の近似予測値は、式１１である。
【数１９】
2007147634.tif 000013
ここで、式１１の△ｋijの係数は、式１２で示される運動方程式をｋijについて微分して得られる感度（１階微分）、式１３である。
【数２０】
2007147634.tif 000014
【数２１】
2007147634.tif 000015
同様に、質量行列［Ｍ］中の或る独立変数ｍijを微小変化させた場合の第ｓ次の固有振動数における変化量の近似予測値は、式１４である。
【数２２】
2007147634.tif 000016
ここで、同定周波数帯域内に実在する固有振動数（実験的特性行列同定法で観察される固有振動数）については、既に制約条件が作成されているので、感度が零となる。そこで、もし感度が零ではない固有振動数が、第１段階の同定で求められた質量行列［Ｍ］および剛性行列［Ｋ］を用いた計算から同定周波数帯域内に得られたら、その感度を用いてニュートン法類の手法で同定周波数帯域より高い周波数にそれら固有振動数を移動すればよい。この計算操作は、第１段階の同定処理における不減衰固有振動数の目標対応化の計算操作と同様に行うことができる。
【００４４】
次に、中央処理装置１２８は、減衰行列［Ｃ］の初期値を設定する（Ｓ２０）。初期値は、乱数などで与えてもよいが、収束時間を短縮する観点から、剛性行列の定数倍を減衰行列の初期値としてもよい。
次に、中央処理装置１２８は、Ｓ１２やＳ１３と同様な処理により、減衰行列［Ｃ］について、一般的固有値問題を解いて、Ｓ８で求めた目標モード減衰比と対応するように、反復の感度解析によりモード減衰比の目標対応化を行う（Ｓ２１）。
【００４５】
次に、中央処理装置１２８は、Ｓ２１処理で求められた減衰行列［Ｃ］が許容範囲内で一致しているか否かを判断し、一致している場合には、Ｓ２４の処理を行い、一致していない場合には、Ｓ２３の処理を行う（Ｓ２２）。このＳ２２の判断は、Ｓ１５の処理と同様に、実験結果から得られた減衰行列［Ｃ］と計算によって得られたそれとの相関性を示す値を計算して行う。
【００４６】
Ｓ２３の処理において、中央処理装置１２８は、式１５により、減衰行列［Ｃ］を修正する（Ｓ２３）。
【数２３】
2007147634.tif 000017
ここで、ζはモード減衰比である。
このＳ２０ないしＳ２３の処理が減衰行列の同定処理である。
【００４７】
以上により、同定周波数帯域内において目標値として設定された物理的に実在する共振モード（実験結果の共振モード）に対応する固有振動数と固有モードが表現され、質量行列［Ｍ］と剛性行列［Ｋ］と減衰行列［Ｃ］との組が導き出される。
次に、中央処理装置１２８は、式１６により、周波数応答関数ｈ（ω）の振幅が実験結果Ｈ（ω）と最も一致するように特性行列を定数倍εする（Ｓ２４）。
【数２４】
2007147634.tif 000018
ここで、ｊは虚数単位である。
次に、中央処理装置１２８は、式１７により、周波数応答関数ｈ（ω）の一致度をさらに高める目的で、周波数応答関数ｈ（ω）の一致度合いを目的関数として最適化を行う（Ｓ２５）。この最適化は、統計学の最ゆう推定法の理論に基づいて同定周波数帯域において、単点加振試験結果である実験周波数応答関数Ｈ（ω）と求められた特性行列の組から得られる周波数応答関数ｈ（ω）とにより計算されるゆう度を最大にすることを目標として、重み付き最小自乗法で行う。
【数２５】
2007147634.tif 000019
ここで、［Ｂ］は重み付き最小自乗法の重み関数である。また、＊は共役転置を表す。
このＳ２４およびＳ２５の処理が第３段階の同定処理である。
次に、中央処理装置１２８は、質量行列［Ｍ］、剛性行列［Ｋ］、減衰行列［Ｃ］および周波数応答関数ｈ（ω）などを出力装置１２６に表示し、さらにこれら結果を記憶装置１２９や外部の記録媒体であるＦＤ１０１などに保存し（２６）、実験的特性行列同定法に基づく振動解析プログラムを終了する（Ｓ２７）。
【００４８】
なお、本実施形態では、加振は、打撃ハンマ１１１によって非定常波を与えたが、これに限定されるものではない。油圧式や圧電式などの加振器、および、音圧や磁界を用いた非接触加振器などを使用することができる。そして、非定常波だけでなく、定常波や周期波など各種の加振波形を使用することができる。これらは、解析対象物１４０における、重さ、大きさ、減衰の大きさ、固有モードの密集性などによって決定される。
【００４９】
そして、本実施形態では、変換器は、加速度計の場合について説明したが、これに限定されるものではない。渦電流やレーザ光線を利用した非接触変位計、レーザドップラ速度計およびひずみゲージなどを使用することができる。
さらに、本実施形態では、記録媒体としてＦＤの場合について説明したが、これに限定されるものではない。メモリチップ（例えば、ＲＯＭチップ）、コンパクトディスクのメモリ（例えば、ＣＤ−ＲＯＭ、ＣＤ−Ｒ、ＣＤ−ＲＷ）、光磁気ディスク（ＭＯ）、光ディスク（ＰＤ）、ディジタル・ビデオ・ディスクのメモリ（例えば、ＤＶＤ−ＲＯＭ、ＤＶＤ−ＲＡＭ）など各種記録媒体を用いることができる。
【００５０】
次に、具体的な実験・計算結果について説明する。
図６は、実験に使用した解析対象物の構造を示す図である。図６（ａ）は、解析対象物の斜視図であり、図６（ｂ）は、部材の断面図であり、図６（ｃ）は、測定点１ないし１０の座標値である。
実験に使用した解析対象物は、図６（ｂ）に示すように断面形状が厚さ２ｍｍで一辺３０ｍｍの正方形である鋼製筒状態の部材１〜５から構成される、図６（ａ）に示すように、漢字の部首で言えば「けいがまえ」状の構造物である。長さ０．８０ｍの部材２、３は、その一方の端を長さ１．００ｍの部材１の両端近くにそれぞれ接合され、他端を長さ０．４０ｍの部材４、５にその中央部でそれぞれ接合される。
【００５１】
加速度計は、図６（ｃ）に示す座標位置でこの解析対象物に密着させて設けられる。
単点加振試験は、この解析対象物を近似的に周辺自由状態でつるし、測定点１をｘ方向に打撃ハンマで加振した。
図７は、剛体特性を示す図である。
【００５２】
図８は、固有振動数およびモード減衰比を示す図である。
図９は、周波数応答関数を示す図である。
図９（ａ）は、周波数応答関数の位相を示し、図９（ｂ）は、アクセレランス（accelerance ）を示す。そして、破線が実験結果を示し、実線が同定結果を示す。
【００５３】
図１０は、１次モードないし５次モードのモード形を示す図である。
図１１は、６次モードないし１０次モードのモード形を示す図である。
この単点加振実験的特性行列同定では、式３に示す［Ｍrigid ］においてα＝１として同定を行った。また、参照値は、鋼製部材の材料定数と解析対象物の形状データとから手計算によって得た近似値である。
【００５４】
図９から分かるように、本発明による単点加振実験的特性行列同定法による計算結果と実験結果とは、よく一致していることが分かる。
このような事例実験を複数回行った結果、慣性モーメントについては、参照値より大きな値として同定される場合があることが分かった。そこで、式３に示す［Ｍrigid ］における補正係数αの最適化について検討した。
【００５５】
図１２は、主慣性モーメントの同定結果を示す図である。
図１２から分かるように、この同定精度を上げるために慣性モーメントを補正する補正係数αは、１．０≦α≦１．３にすることが好ましい。さらに、事例結果に統計処理を行うと、補正係数αは、約１．１７であることが好ましい。
一方、目標モード特性の設定の際に使用する共振周波数について検討を行った。
【００５６】
図１３は、同定周波数帯域６Ｈｚないし５０Ｈｚの第１次ないし第３次の共振周波数を同定対象とした場合の質量行列および剛体特性を示す図である。
図１４は、同定周波数帯域６Ｈｚないし２６Ｈｚの第１次ないし第３次の共振周波数を同定対象とした場合の質量行列および剛体特性を示す図である。
図１５は、同定周波数帯域６Ｈｚないし２０Ｈｚの第１次の共振周波数を同定対象とした場合の質量行列および剛体特性を示す図である。
【００５７】
この単点加振実験的特性行列同定では、図６に示す解析対象物に対し、式３に示す［Ｍrigid ］においてα＝１として同定を行った。
図１３ないし図１５から分かるように、同定対象とする共振周波数を第１次のみ、または第１次から第３次までを対象にして実用的に満足な精度の結果を得ることができる。より多くの共振周波数を用いた方が同定に際して有利であるから、結局、同定対象とする共振周波数を第１次のみ、または第１次から任意の次数までを対象にして実用的に満足な精度の結果を得ることができる。
【産業上の利用可能性】
【００５８】
以上説明したように、本発明は振動解析プログラムなどに利用可能な技術である。
【図面の簡単な説明】
【００５９】
【図１】本実施形態の解析装置のブロック図を示す図である。
【図２】本実施形態における解析装置のフローチャート・その１を示す図である。
【図３】本実施形態における解析装置のフローチャート・その２を示す図である。
【図４】本実施形態における解析装置のフローチャート・その３を示す図である。
【図５】本実施形態における解析装置のフローチャート・その４を示す図である。
【図６】実験に使用した解析対象物の構造を示す図である。
【図７】剛体特性を示す図である。
【図８】固有振動数およびモード減衰比を示す図である。
【図９】周波数応答関数を示す図である。
【図１０】１次モードないし５次モードのモード形を示す図である。
【図１１】１次モードないし５次モードのモード形を示す図である。
【図１２】主慣性モーメントの同定結果を示す図である。
【図１３】同定周波数帯域６Ｈｚないし５０Ｈｚの第１次ないし第３次の共振周波数を同定対象とした場合の質量行列および剛体特性を示す図である。
【図１４】同定周波数帯域６Ｈｚないし２６Ｈｚの第１次ないし第３次の共振周波数を同定対象とした場合の質量行列および剛体特性を示す図である。
【図１５】同定周波数帯域６Ｈｚないし２０Ｈｚの第１次の共振周波数を同定対象とした場合の質量行列および剛体特性を示す図である。
【符号の説明】
【００６０】
１１１  打撃ハンマ
１１２  力変換器
１１３  加速度計
１２１  外部インタフェイス
１０１  フレキシブルディスク
１２２  ＦＤ装置
１２８  中央処理装置
１２９  記憶装置
【図１】
2007147634.tif 000028
【図２】
2007147634.tif 000029
【図３】
2007147634.tif 000030
【図４】
2007147634.tif 000031
【図５】
2007147634.tif 000032
【図６】
2007147634.tif 000033
【図７】
2007147634.tif 000034
【図８】
2007147634.tif 000035
【図９】
2007147634.tif 000036
【図１０】
2007147634.tif 000037
【図１１】
2007147634.tif 000038
【図１２】
2007147634.tif 000039
【図１３】
2007147634.tif 000040
【図１４】
2007147634.tif 000041
【図１５】
2007147634.tif 000042