JP3335765B2 - 振動特性解析装置 - Google Patents
振動特性解析装置Info
- Publication number
- JP3335765B2 JP3335765B2 JP14008194A JP14008194A JP3335765B2 JP 3335765 B2 JP3335765 B2 JP 3335765B2 JP 14008194 A JP14008194 A JP 14008194A JP 14008194 A JP14008194 A JP 14008194A JP 3335765 B2 JP3335765 B2 JP 3335765B2
- Authority
- JP
- Japan
- Prior art keywords
- term
- equation
- parameter
- nonlinear
- linear
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Landscapes
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
Description
【０００１】
【産業上の利用分野】本発明は振動特性解析装置に関
し、特にエンジンブロック等の被試験物を加振すること
により応答を得てその応答関数から被試験物のモード特
性を同定する振動特性解析装置に関するものである。
し、特にエンジンブロック等の被試験物を加振すること
により応答を得てその応答関数から被試験物のモード特
性を同定する振動特性解析装置に関するものである。
【０００２】実際の機械構造物の動特性を解析するため
には、該構造物を加振試験することにより得られた振動
応答特性から該構造物の動特性を同定する必要がある
が、この場合、振動現象を直接表現するパラメータであ
り、現象を直接理解し易く、収束性が良好なモード特性
（固有振動数、モード減衰比、固有モード形状などのモ
ーダルパラメータ）を同定する手法が現在では主流とな
っている。
には、該構造物を加振試験することにより得られた振動
応答特性から該構造物の動特性を同定する必要がある
が、この場合、振動現象を直接表現するパラメータであ
り、現象を直接理解し易く、収束性が良好なモード特性
（固有振動数、モード減衰比、固有モード形状などのモ
ーダルパラメータ）を同定する手法が現在では主流とな
っている。
【０００３】
【従来の技術】図１は従来より知られている振動特性解
析装置の構成を概略的に示したもので、まず被試験物１
とこの被試験物１を加振するためのインパルスハンマー
２とを用意する。
析装置の構成を概略的に示したもので、まず被試験物１
とこの被試験物１を加振するためのインパルスハンマー
２とを用意する。
【０００４】そして、ハンマー２には力検出用ロードセ
ル３を取り付け、また該被試験物１の任意の場所に該ハ
ンマー２による加振力の３次元（ｘ，ｙ，ｚ）方向加速
度を検出するためのピックアップ４を取り付ける。
ル３を取り付け、また該被試験物１の任意の場所に該ハ
ンマー２による加振力の３次元（ｘ，ｙ，ｚ）方向加速
度を検出するためのピックアップ４を取り付ける。
【０００５】そして、ピックアップ４の出力信号とロー
ドセル３の出力信号とを演算装置５に与えてＦＦＴ処理
を行い周波数応答関数（伝達関数：コンプライアンス）
を求めモード特性を計算して出力する構成を有してい
る。
ドセル３の出力信号とを演算装置５に与えてＦＦＴ処理
を行い周波数応答関数（伝達関数：コンプライアンス）
を求めモード特性を計算して出力する構成を有してい
る。
【０００６】この場合、加振実験は通常、ハンマー２に
よる被試験物１の加振場所を固定して行い、ピックアッ
プ４は逐次移動させて複数の応答関数を演算装置５に与
えるようにする。なお、被試験物１は柔らかいバネ（図
示せず）等により固定されている。
よる被試験物１の加振場所を固定して行い、ピックアッ
プ４は逐次移動させて複数の応答関数を演算装置５に与
えるようにする。なお、被試験物１は柔らかいバネ（図
示せず）等により固定されている。
【０００７】このような振動特性解析装置の演算装置
（実験モード解析装置）５のアルゴリズムとしては今ま
で多くの手法が提案されているが、この中で周波数領域
における『偏分反復法』は精度の良さで現在最も広く使
用されている。
（実験モード解析装置）５のアルゴリズムとしては今ま
で多くの手法が提案されているが、この中で周波数領域
における『偏分反復法』は精度の良さで現在最も広く使
用されている。
【０００８】この偏分反復法は簡単に言えば、実験デー
タと理論値の誤差を検定しながら反復計算を行い、最小
二乗法により最も確からしいモード特性を同定するもの
である。更に周波数領域で自由に重み付けすることがで
き、注目するモード特性以外のモード特性の影響を考慮
することができる。
タと理論値の誤差を検定しながら反復計算を行い、最小
二乗法により最も確からしいモード特性を同定するもの
である。更に周波数領域で自由に重み付けすることがで
き、注目するモード特性以外のモード特性の影響を考慮
することができる。
【０００９】＜偏分反復法＞以下にこの偏分反復法につ
いて説明すると、まず、ハンマー２で被試験物１を加振
したときのピックアップ４から演算装置５へ与えられる
モード特性の応答関数（理論値）は振幅と位相を持つの
で、次式のように実部ＧR（ω）と虚部ＧI（ω）とから
成る複素数で表されることがよく知られている。
いて説明すると、まず、ハンマー２で被試験物１を加振
したときのピックアップ４から演算装置５へ与えられる
モード特性の応答関数（理論値）は振幅と位相を持つの
で、次式のように実部ＧR（ω）と虚部ＧI（ω）とから
成る複素数で表されることがよく知られている。
【００１０】
【数１】
【００１１】この場合、角振動数をωで表し、注目する
角振動数範囲内の固有モ−ド数（自由度）をｎ（後述す
る図７の共振点を示す山の数に相当する）としている。
なお、Ｕr，Ｖrは、測定点に対するｒ番目の留数の実部
及び虚部を示している。
角振動数範囲内の固有モ−ド数（自由度）をｎ（後述す
る図７の共振点を示す山の数に相当する）としている。
なお、Ｕr，Ｖrは、測定点に対するｒ番目の留数の実部
及び虚部を示している。
【００１２】そして、この複素数とピックアップ４から
得られる実験データ（ＦＦＴ演算処理した後の周波数を
横軸とし応答関数（コンプライアンス＝変位／力）を縦
軸としたデータ）との誤差の二乗和（残差二乗和）が最
小になるように式(1) 中のパラメータ（定数）を決定
し、そのパラメータからモード特性を決定する。
得られる実験データ（ＦＦＴ演算処理した後の周波数を
横軸とし応答関数（コンプライアンス＝変位／力）を縦
軸としたデータ）との誤差の二乗和（残差二乗和）が最
小になるように式(1) 中のパラメータ（定数）を決定
し、そのパラメータからモード特性を決定する。
【００１３】即ち、式(1) に対応する実験デ−タから最
も確からしいパラメータ｛γ｝T＝｛‥,σr,ωdr,Ｕr,
Ｖr,‥,Ｃ,Ｄ,Ｅ,Ｆ,‥｝（｛‥｝はベクトルを表示す
る）を求めるには、理論値に対する実験デ−タの誤差が
正規分布するとして残差二乗和Ｓが最小となるようにパ
ラメータ｛γ｝を初期値（これは図７の実験データＥの
共振点の周波数を読み取った値）から順次修正して行け
ば、図７に示すように実験データＥ（実線）の周波数応
答特性に最も近い（確からしい）理論値Ｔ（一点鎖線）
が得られることとなる。この残差二乗和Ｓは次式で表さ
れる。
も確からしいパラメータ｛γ｝T＝｛‥,σr,ωdr,Ｕr,
Ｖr,‥,Ｃ,Ｄ,Ｅ,Ｆ,‥｝（｛‥｝はベクトルを表示す
る）を求めるには、理論値に対する実験デ−タの誤差が
正規分布するとして残差二乗和Ｓが最小となるようにパ
ラメータ｛γ｝を初期値（これは図７の実験データＥの
共振点の周波数を読み取った値）から順次修正して行け
ば、図７に示すように実験データＥ（実線）の周波数応
答特性に最も近い（確からしい）理論値Ｔ（一点鎖線）
が得られることとなる。この残差二乗和Ｓは次式で表さ
れる。
【００１４】
【数２】
【００１５】なお、複数の応答関数を同時に処理する場
合にはそれらの実部、虚部を各々直列に並べて考えれば
よく、応答関数の数（測定点の数）をｍ，１つの応答関
数（１回の加振）中のデータ数をＮとするとｐ＝ｍ・Ｎ
となる。また、Ｗi は重みであり各々の実測データの分
散の逆数に比例する。
合にはそれらの実部、虚部を各々直列に並べて考えれば
よく、応答関数の数（測定点の数）をｍ，１つの応答関
数（１回の加振）中のデータ数をＮとするとｐ＝ｍ・Ｎ
となる。また、Ｗi は重みであり各々の実測データの分
散の逆数に比例する。
【００１６】残差二乗和Ｓが最小となるパラメータ
｛γ｝を求めるには、パラメータ｛γ｝の中に式(2) に
おける応答関数（理論式）ｆの分母にσr,ωdrなどの非
線形項を含むので、線形近似による簡略化された改良反
復法（ガウス−ニュートン法）を次式のとおり適用する
ことによりパラメータ｛γ｝を修正するための修正ベク
トル｛△γ｝を求めて残差二乗和Ｓが最小値に収束する
まで反復して計算を行えばよい。
｛γ｝を求めるには、パラメータ｛γ｝の中に式(2) に
おける応答関数（理論式）ｆの分母にσr,ωdrなどの非
線形項を含むので、線形近似による簡略化された改良反
復法（ガウス−ニュートン法）を次式のとおり適用する
ことによりパラメータ｛γ｝を修正するための修正ベク
トル｛△γ｝を求めて残差二乗和Ｓが最小値に収束する
まで反復して計算を行えばよい。
【００１７】
【数３】 ｛△γ｝＝([Ａ]T[Ｗ][Ａ])-1・[Ａ]T[Ｗ]｛ｙi−ｆi｝ ・・・式(3) ただし、Ａij ＝∂ｆi ／∂γj であり、これは理論式
(1) 中のγの一つを変化させたときに該理論式(1) がど
れだけ変化したかを示しているが、これにより非線形項
を含む分母を強制的に分子に変換する「線形近似」を実
現している。また、[Ｗ]はＷiを対角成分とする対角行
列を示している。
(1) 中のγの一つを変化させたときに該理論式(1) がど
れだけ変化したかを示しているが、これにより非線形項
を含む分母を強制的に分子に変換する「線形近似」を実
現している。また、[Ｗ]はＷiを対角成分とする対角行
列を示している。
【００１８】
【発明が解決しょうとする課題】しかし、上記の偏分反
復法には次の２つの欠点がある。 特に応答関数が複数の場合、求めるパラメータ｛γ｝
も応答関数の増加に応じて増加するため、式(3) の逆行
列の計算時間が急増してしまう（図３の特性Ａ参照）。 変数が増加するため反復計算が収束し難く発散し易く
なり、途中で計算不能になることがある（図８参照）。
復法には次の２つの欠点がある。 特に応答関数が複数の場合、求めるパラメータ｛γ｝
も応答関数の増加に応じて増加するため、式(3) の逆行
列の計算時間が急増してしまう（図３の特性Ａ参照）。 変数が増加するため反復計算が収束し難く発散し易く
なり、途中で計算不能になることがある（図８参照）。
【００１９】従って本発明は、被試験物を加振するため
のインパルスハンマーに取り付けられた力検出用ロード
セルと、該被試験物の任意の場所に取り付けられて該加
振力の３次元方向加速度を検出するためのピックアップ
と、該ピックアップの出力信号と該ロードセルの出力信
号とを受けて周波数応答関数を求めモード特性パラメー
タを計算することにより該被試験物のモード特性を同定
する演算装置とを備えた振動特性解析装置において、こ
の発散し易く計算時間がかかるというアルゴリズムの本
質に係わる改良を実現することを目的とする。
のインパルスハンマーに取り付けられた力検出用ロード
セルと、該被試験物の任意の場所に取り付けられて該加
振力の３次元方向加速度を検出するためのピックアップ
と、該ピックアップの出力信号と該ロードセルの出力信
号とを受けて周波数応答関数を求めモード特性パラメー
タを計算することにより該被試験物のモード特性を同定
する演算装置とを備えた振動特性解析装置において、こ
の発散し易く計算時間がかかるというアルゴリズムの本
質に係わる改良を実現することを目的とする。
【００２０】
【課題を解決するための手段】上記の目的を達成するた
め、本発明に係る振動特性解析装置は、演算装置が、実
測データとモード特性パラメータによる理論値との残差
が最小になるまで該パラメータを繰り返し変えて行く偏
分反復法を線形項と非線形項とに分けて計算すると共に
該線形項及び非線形項が変化しても該残差が発散しない
条件を設けることにより該モード特性を同定することを
特徴としたものである。
め、本発明に係る振動特性解析装置は、演算装置が、実
測データとモード特性パラメータによる理論値との残差
が最小になるまで該パラメータを繰り返し変えて行く偏
分反復法を線形項と非線形項とに分けて計算すると共に
該線形項及び非線形項が変化しても該残差が発散しない
条件を設けることにより該モード特性を同定することを
特徴としたものである。
【００２１】また上記の演算装置は、該発散しない条件
が該残差に対する該線形項及び非線形項による偏微分が
ゼロとなることであり、このときの２次偏微分項を省略
することができる。
が該残差に対する該線形項及び非線形項による偏微分が
ゼロとなることであり、このときの２次偏微分項を省略
することができる。
【００２２】更に上記の演算装置は、上記の発散しない
条件として該偏分反復法を実行するときの該非線形項の
修正ベクトルに対して閾値以内の縮小条件を付加するこ
とが好ましい。
条件として該偏分反復法を実行するときの該非線形項の
修正ベクトルに対して閾値以内の縮小条件を付加するこ
とが好ましい。
【００２３】
【作用】本発明に係る振動特性解析装置においては、モ
ード特性パラメータ線形項と非線形項とに分割する。
ード特性パラメータ線形項と非線形項とに分割する。
【００２４】即ち、非線形項パラメータの数が、応答関
数の数が増加しても図７に示す山の位置（周波数）が変
わらないため不変であり、応答関数の数の影響を受けな
いことが特徴であり、この非線形項が求まれば線形項を
容易に求めることができるからである。
数の数が増加しても図７に示す山の位置（周波数）が変
わらないため不変であり、応答関数の数の影響を受けな
いことが特徴であり、この非線形項が求まれば線形項を
容易に求めることができるからである。
【００２５】そこで非線形項だけを先ず非線形パラメー
タの初期値（暫定値）を設定し、この初期値により線形
項パラメータを計算で求める。
タの初期値（暫定値）を設定し、この初期値により線形
項パラメータを計算で求める。
【００２６】このようにして求めた線形項と初期値の非
線形項とを用いて式(1) により応答関数ＧR 及びＧI を
計算する。
線形項とを用いて式(1) により応答関数ＧR 及びＧI を
計算する。
【００２７】このようにして計算した応答関数ＧR ，Ｇ
I と実験値との残差二乗和Ｓを式(2) に従って計算し、
この残差二乗和Ｓが最小値（極小値）であると判ったと
きには、上記のようにして求めた線形項と非線形項を同
定したいモード特性の最終値であるとする。
I と実験値との残差二乗和Ｓを式(2) に従って計算し、
この残差二乗和Ｓが最小値（極小値）であると判ったと
きには、上記のようにして求めた線形項と非線形項を同
定したいモード特性の最終値であるとする。
【００２８】残差二乗和Ｓが最小値ではないと判ったと
きには残差二乗和Ｓは、この残差二乗和Ｓを最小値にす
るために新しい非線形項を求める。
きには残差二乗和Ｓは、この残差二乗和Ｓを最小値にす
るために新しい非線形項を求める。
【００２９】この非線形項だけを求めるには修正ベクト
ル｛△β｝を上記の式(3) と同様にして求め、収束する
まで反復計算を行うが、線形項のことを考慮するため
に、残差二乗和Ｓの極小点近傍において発散防止を図る
ための制約条件を設け、この制約条件により非線形項の
修正ベクトル｛△β｝が求められる。
ル｛△β｝を上記の式(3) と同様にして求め、収束する
まで反復計算を行うが、線形項のことを考慮するため
に、残差二乗和Ｓの極小点近傍において発散防止を図る
ための制約条件を設け、この制約条件により非線形項の
修正ベクトル｛△β｝が求められる。
【００３０】但し、修正ベクトル｛△β｝は求まるが、
このまま使用すると発散する場合もあるので更に側面制
約を設定する。
このまま使用すると発散する場合もあるので更に側面制
約を設定する。
【００３１】このようにして求められた修正ベクトル
｛△β｝と前回の非線形項パラメータ｛β｝（最初は初
期値）とを加算して今回の新しい非線形項パラメーを計
算し、この新しい非線形項パラメータを用いて繰り返し
実行することにより、最適なモード特性を示す｛α｝と
｛β｝を決定することができる。
｛△β｝と前回の非線形項パラメータ｛β｝（最初は初
期値）とを加算して今回の新しい非線形項パラメーを計
算し、この新しい非線形項パラメータを用いて繰り返し
実行することにより、最適なモード特性を示す｛α｝と
｛β｝を決定することができる。
【００３２】
【実施例】本発明に係る振動特性解析装置の構成は上述
した図１に示した従来例と同様のものを用いることがで
きるが、従来例との差異は演算装置５における演算処理
である。以下にこの演算処理について図２に示したフロ
ーチャートに沿って説明する。
した図１に示した従来例と同様のものを用いることがで
きるが、従来例との差異は演算装置５における演算処理
である。以下にこの演算処理について図２に示したフロ
ーチャートに沿って説明する。
【００３３】(1) パラメータの分割：複数の応答関数を
同時に考慮する時、その応答関数の数をｍ，同定する固
有モ−ドの数を上記と同様にｎとするとパラメータ
｛γ｝の数は2n+(2n+4)mである。そのパラメータ｛γ｝
のうち式(1) において分母中にある非線形項は２ｎ、分
子中にある線形項は(2n+4)m である。このパラメータ
｛γ｝を次式のように線形項｛α｝と非線形項｛β｝と
に分割する。
同時に考慮する時、その応答関数の数をｍ，同定する固
有モ−ドの数を上記と同様にｎとするとパラメータ
｛γ｝の数は2n+(2n+4)mである。そのパラメータ｛γ｝
のうち式(1) において分母中にある非線形項は２ｎ、分
子中にある線形項は(2n+4)m である。このパラメータ
｛γ｝を次式のように線形項｛α｝と非線形項｛β｝と
に分割する。
【００３４】
【数４】 ｛γ｝T＝｛｛α｝T,｛β｝T｝ ｛α｝T＝｛‥,Ｕr,Ｖr,‥,Ｃ,Ｄ,Ｅ,Ｆ,‥｝ ｛β｝T＝｛‥,σr,ωdr,‥｝ ・・・式(4)
【００３５】ここで、非線形項｛β｝中のσr は各ピー
ク（ｎ個存在し得る図７に示す山々）、つまり各モード
特性における減衰率を示し、ωdrは各ピーク、つまり各
モード特性における減衰固有角振動数を示す。この減衰
率と固有振動数は応答点（ピックアップ４の取付位置）
が変わっても変化しない。つまり非線形項のパラメータ
｛β｝の数２ｎは応答関数の数ｍが増加しても図７に示
す山の位置（周波数）が変わらないため不変であること
が特徴となっており、これにより応答関数の数の影響を
受けないことが分かる。そして、非線形項が求まれば線
形項だけは容易に求めることができる。
ク（ｎ個存在し得る図７に示す山々）、つまり各モード
特性における減衰率を示し、ωdrは各ピーク、つまり各
モード特性における減衰固有角振動数を示す。この減衰
率と固有振動数は応答点（ピックアップ４の取付位置）
が変わっても変化しない。つまり非線形項のパラメータ
｛β｝の数２ｎは応答関数の数ｍが増加しても図７に示
す山の位置（周波数）が変わらないため不変であること
が特徴となっており、これにより応答関数の数の影響を
受けないことが分かる。そして、非線形項が求まれば線
形項だけは容易に求めることができる。
【００３６】そこで非線形項だけを先ず求めることを工
夫する。このため、まず非線形パラメータ｛β｝の初期
値を設定する（図２のステップＳ１）。これは、図７に
示すような周波数応答のピークの周波数を読んでその値
をωdrとし、σrは適当な値を設定することにより行わ
れる。
夫する。このため、まず非線形パラメータ｛β｝の初期
値を設定する（図２のステップＳ１）。これは、図７に
示すような周波数応答のピークの周波数を読んでその値
をωdrとし、σrは適当な値を設定することにより行わ
れる。
【００３７】非線形項パラメータ｛β｝が初期設定され
れば、次式により線形項パラメータ｛α｝を１回の計算
で求めることができる（ステップＳ２）。
れば、次式により線形項パラメータ｛α｝を１回の計算
で求めることができる（ステップＳ２）。
【００３８】
【数５】 ｛α｝＝([Ｄ]T[Ｗ][Ｄ])-1[Ｄ]T[Ｗ]｛ｙi｝ ・・・式(5) ただし、Ｄij＝∂ｆi／∂αj （ｉ＝1〜2p,ｊ＝1〜(2
n+4)m)
n+4)m)
【００３９】ここで、式(4) に示したように、｛α｝T
＝｛‥,Ｕr,Ｖr,‥,Ｃ,Ｄ,Ｅ,Ｆ,‥｝であるので、Ｄij
は次のようになる。
＝｛‥,Ｕr,Ｖr,‥,Ｃ,Ｄ,Ｅ,Ｆ,‥｝であるので、Ｄij
は次のようになる。
【００４０】まず、ｉ＝１〜ｐのときは、ｆi＝Ｇ
R(ωi) であるので、Ｄij＝∂ＧR(ωi)／∂αj とな
り、各線形項について次式が得られる。
R(ωi) であるので、Ｄij＝∂ＧR(ωi)／∂αj とな
り、各線形項について次式が得られる。
【００４１】
【数６】
【００４２】また、上記式(6) におけるａr,ｂr も式
(1) で示した通り非線形項σr,ωdrで表すことができ、
非線形項σr,ωdrが上記のとおり初期設定されているの
で、それらの初期値を入れ且つωi を入れることにより
Ｄijが計算できる。
(1) で示した通り非線形項σr,ωdrで表すことができ、
非線形項σr,ωdrが上記のとおり初期設定されているの
で、それらの初期値を入れ且つωi を入れることにより
Ｄijが計算できる。
【００４３】同様に、ｉ＝ｐ＋１〜２ｐのときは、ｆ
i＝ＧI(ωi) であるので、Ｄij＝∂ＧI(ωi)／∂αj と
なり、各線形項について次式が得られる。
i＝ＧI(ωi) であるので、Ｄij＝∂ＧI(ωi)／∂αj と
なり、各線形項について次式が得られる。
【００４４】
【数７】
【００４５】この式においても、非線形項の初期値を入
れ且つωi を入れることによりＤijが計算できる。
れ且つωi を入れることによりＤijが計算できる。
【００４６】このようにして求めた線形項｛α｝と非線
形項｛β｝とを用いて式(1) により応答関数ＧR 及びＧ
I を計算する（ステップＳ３）。
形項｛β｝とを用いて式(1) により応答関数ＧR 及びＧ
I を計算する（ステップＳ３）。
【００４７】このようにして計算した応答関数ＧR ，Ｇ
I と実験値との残差二乗和Ｓを式(2) に従って計算する
（ステップＳ４）。
I と実験値との残差二乗和Ｓを式(2) に従って計算する
（ステップＳ４）。
【００４８】そして、この残差二乗和Ｓが前回求めた値
より小さくなっているか否かを判定する（ステップＳ
５）。なお、この場合の最初の残差二乗和Ｓは適当な値
を設定しておけばよい。
より小さくなっているか否かを判定する（ステップＳ
５）。なお、この場合の最初の残差二乗和Ｓは適当な値
を設定しておけばよい。
【００４９】この結果、残差二乗和Ｓが前回の値より小
さくなっていないときには残差二乗和Ｓが最小値に達し
たことになるので、上記のようにして求めた線形項
｛α｝と非線形項｛β｝を同定したいモード特性の最終
値であるとして出力し（ステップＳ６）、このルーチン
を終了する。
さくなっていないときには残差二乗和Ｓが最小値に達し
たことになるので、上記のようにして求めた線形項
｛α｝と非線形項｛β｝を同定したいモード特性の最終
値であるとして出力し（ステップＳ６）、このルーチン
を終了する。
【００５０】残差二乗和Ｓが前回の値より小さくなって
いるときには残差二乗和Ｓはまだ最小値に達していない
ことになるので、この残差二乗和Ｓを最小値にするため
の次の非線形項パラメータ｛β｝を以下のとおり求め
る。
いるときには残差二乗和Ｓはまだ最小値に達していない
ことになるので、この残差二乗和Ｓを最小値にするため
の次の非線形項パラメータ｛β｝を以下のとおり求め
る。
【００５１】即ち、非線形項だけを求めるには次式によ
り修正ベクトル｛△β｝を上記の式(3) における｛γ｝
と同様にして求め、収束するまで反復計算を行う。
り修正ベクトル｛△β｝を上記の式(3) における｛γ｝
と同様にして求め、収束するまで反復計算を行う。
【００５２】
【数８】 ｛△β｝＝([Ｃ]T[Ｗ][Ｃ])-1[Ｃ]T[Ｗ]・｛ｙi−ｆi｝ ・・・式(8) ただし、Ｃij＝∂ｆi／∂βj （ｉ＝1〜2p,ｊ＝1〜2
n）
n）
【００５３】しかしながら、このままでは線形項｛α｝
のことを考慮していないので、線形項｛α｝を考慮する
ために次のように工夫する。
のことを考慮していないので、線形項｛α｝を考慮する
ために次のように工夫する。
【００５４】残差二乗和Ｓの極小点において次式が成立
する。
する。
【００５５】
【数９】 ∂Ｓ／∂αj＝０ （ｊ＝1〜(2n+4)m ） ・・・式(9)
【００５６】この近傍において次式が成立すれば発散防
止を図ることができる。
止を図ることができる。
【００５７】
【数１０】 ∂2Ｓ／∂αj∂βk＝０ （ｊ＝1〜(2n+4)m,ｋ＝1〜2n） ・・・式(10)
【００５８】即ち、この式(10)の意味は、∂Ｓ／∂αj
が線形項を変えると残差二乗和Ｓがどれだけ変わるかを
示し、これを更に∂βkで偏微分することにより非線形
項を変えても∂Ｓ／∂αjが変化しないということを示
し（右辺が０のため）、以て線形項｛α｝を考慮した発
散防止を図るための制約条件となっている。
が線形項を変えると残差二乗和Ｓがどれだけ変わるかを
示し、これを更に∂βkで偏微分することにより非線形
項を変えても∂Ｓ／∂αjが変化しないということを示
し（右辺が０のため）、以て線形項｛α｝を考慮した発
散防止を図るための制約条件となっている。
【００５９】上記の式(10)の数は、2n×(2n+4)m であり
線形項｛α｝を消去するには多過ぎるが、ここでは式(1
0)を全て成立させるのではなく変数を縮小するのが狙い
であるので、式(10)の下での応答関数ｆのβk による微
分係数をＨk とすると、この微分係数Ｈk は次式のよう
に表される計算される（ステップＳ７）。
線形項｛α｝を消去するには多過ぎるが、ここでは式(1
0)を全て成立させるのではなく変数を縮小するのが狙い
であるので、式(10)の下での応答関数ｆのβk による微
分係数をＨk とすると、この微分係数Ｈk は次式のよう
に表される計算される（ステップＳ７）。
【００６０】
【数１１】
【００６１】この微分係数Ｈk を求めるためには、∂α
j／∂βk（ｊ＝1〜(2n+4)m,ｋ＝1〜2n）を計算する必要
がある。
j／∂βk（ｊ＝1〜(2n+4)m,ｋ＝1〜2n）を計算する必要
がある。
【００６２】このため、まず式(10)に式(2) の残差二乗
和Ｓを代入すると、次式のようになる。
和Ｓを代入すると、次式のようになる。
【００６３】
【数１２】
【００６４】この式(12)の左辺第１項の∂ｆi／∂βkの
代わりに式(11)のＨk を代入し、更に整理して行列表示
すると次式のようになる。
代わりに式(11)のＨk を代入し、更に整理して行列表示
すると次式のようになる。
【００６５】
【数１３】
【００６６】従って、式(13)を変形して、｛ａ｝＝
｛ｂ｝[Ｂ]-1とすることにより、∂αj／∂βkを求める
ことができる。
｛ｂ｝[Ｂ]-1とすることにより、∂αj／∂βkを求める
ことができる。
【００６７】また、［Ｂ］は(2n+4)m ×(2n+4)m の正方
行列であるが、Ｂstにおいて、ｓ，ｔが異なる応答関数
に属する時には［Ｂ］は次式に示すようにＢst＝０とな
るので(2n+4)×(2n+4)の行列ｍ個に分割できる。
行列であるが、Ｂstにおいて、ｓ，ｔが異なる応答関数
に属する時には［Ｂ］は次式に示すようにＢst＝０とな
るので(2n+4)×(2n+4)の行列ｍ個に分割できる。
【００６８】
【数１４】
【００６９】従って、｛ａ｝を求めるための上記の逆行
列[Ｂ]-1の大きさは応答関数の数ｍに無関係となる。
列[Ｂ]-1の大きさは応答関数の数ｍに無関係となる。
【００７０】そして、式(11)の微分係数Ｈk が求まれば
式(3) と同様に次式から非線形項の修正ベクトル｛△
β｝が求められる（ステップＳ８）。
式(3) と同様に次式から非線形項の修正ベクトル｛△
β｝が求められる（ステップＳ８）。
【００７１】
【数１５】 ｛△β｝＝([Ｃ]T[Ｗ][Ｃ])-1[Ｃ]T[Ｗ]・｛ｙi−ｆi｝ ・・・式(15) ただし、Ｃij＝Ｈj（ｆ＝ｆi） （ｉ＝1〜2p,ｊ＝1〜2
n）
n）
【００７２】＜側面制約の設定＞式(15)により非線形項
の修正ベクトル｛△β｝が求まるが、上記の式(10)の制
約条件にも関わらずこのまま使用すると発散する場合が
あるので、更に側面制約を設定する。
の修正ベクトル｛△β｝が求まるが、上記の式(10)の制
約条件にも関わらずこのまま使用すると発散する場合が
あるので、更に側面制約を設定する。
【００７３】発散防止の側面制約としては修正ベクトル
に縮小因子をかけてその方向は同じで大きさを縮小する
方法があるが、この場合、モ−ドが異なれば非線形項相
互の影響は少なくなるので、非線形項各々に次式の制約
を設定する（ステップＳ９）。
に縮小因子をかけてその方向は同じで大きさを縮小する
方法があるが、この場合、モ−ドが異なれば非線形項相
互の影響は少なくなるので、非線形項各々に次式の制約
を設定する（ステップＳ９）。
【００７４】
【数１６】 ｜△ωdr／ωdr｜≦ Ｋ1 ｜△σr ／σr ｜≦ Ｋ2 ， σr ＞０ ・・・ 式(16)
【００７５】今回、この係数Ｋ1 は0.05、Ｋ2 は0.8 を
使用した。Ｋ2 ＜１で且つσr の初期値が正であればσ
r は繰り返し計算後も正となる。
使用した。Ｋ2 ＜１で且つσr の初期値が正であればσ
r は繰り返し計算後も正となる。
【００７６】このようにして求められた修正ベクトル
｛△β｝と前回の非線形項パラメータ｛β｝（最初は初
期値）とを加算して今回の新しい非線形項パラメータ
｛β*｝を計算する（ステップＳ１０）。
｛△β｝と前回の非線形項パラメータ｛β｝（最初は初
期値）とを加算して今回の新しい非線形項パラメータ
｛β*｝を計算する（ステップＳ１０）。
【００７７】
【数１７】 ｛β*｝＝｛β｝＋｛Δβ｝ ・・・式(17)
【００７８】このようにして求めた新しい非線形項パラ
メータ｛β*｝を用いて上記のステップＳ３以降を繰り
返し実行することとなり、最適なモード特性を示す
｛α｝と｛β｝を出力する（ステップＳ６）こととな
る。
メータ｛β*｝を用いて上記のステップＳ３以降を繰り
返し実行することとなり、最適なモード特性を示す
｛α｝と｛β｝を出力する（ステップＳ６）こととな
る。
【００７９】なお、本発明では、式(12)の第２項目であ
る応答関数ｆの２階偏微分項を無視して同様に展開する
事もできる。
る応答関数ｆの２階偏微分項を無視して同様に展開する
事もできる。
【００８０】
【発明の効果】以上説明したように、本発明に係る振動
特性解析装置によれば、実測データとモード特性パラメ
ータによる理論値との残差が最小になるまで該パラメー
タを繰り返し変えて行く偏分反復法を線形項と非線形項
とに分けて計算すると共に該線形項及び非線形項が変化
しても該残差が発散しない条件を設けることにより該モ
ード特性を同定するように構成したので、次に示す如
く、従来技術に対する効果が得られる。 （１）計算時間について：計算時間の比較を図３と下記
の表１に示す。同定するモ−ド数ｎは５で１つのコンプ
ライアンス中のデータ数は５７３である。単点参照にお
いては従来例によるほうが計算機使用時間が少ないが、
従来例では、周波数応答関数のｍが増加すると計算機使
用時間は急増する。これは解くべき逆行列の大きさが
（２ｎ＋４）ｍとなる為である。改良型は逆行列の大き
さは２ｎで一定であり、これを使用する計算回数が増加
するだけなのでほぼ直線的に増加している。
特性解析装置によれば、実測データとモード特性パラメ
ータによる理論値との残差が最小になるまで該パラメー
タを繰り返し変えて行く偏分反復法を線形項と非線形項
とに分けて計算すると共に該線形項及び非線形項が変化
しても該残差が発散しない条件を設けることにより該モ
ード特性を同定するように構成したので、次に示す如
く、従来技術に対する効果が得られる。 （１）計算時間について：計算時間の比較を図３と下記
の表１に示す。同定するモ−ド数ｎは５で１つのコンプ
ライアンス中のデータ数は５７３である。単点参照にお
いては従来例によるほうが計算機使用時間が少ないが、
従来例では、周波数応答関数のｍが増加すると計算機使
用時間は急増する。これは解くべき逆行列の大きさが
（２ｎ＋４）ｍとなる為である。改良型は逆行列の大き
さは２ｎで一定であり、これを使用する計算回数が増加
するだけなのでほぼ直線的に増加している。
【００８１】
【表１】
【００８２】（２）収束性について：収束性を比較する
為に計算で作成したコンプライアンスを使用する。モデ
ルは５質点、５バネで１番バネの端を固定している。諸
元と不減衰固有振動数を下記の表２に示す。同定するモ
ード特性数は５で、データは50〜950Hz までの１Hz置き
で 901個である。モード特性減衰比は全モード特性で１
％とした。
為に計算で作成したコンプライアンスを使用する。モデ
ルは５質点、５バネで１番バネの端を固定している。諸
元と不減衰固有振動数を下記の表２に示す。同定するモ
ード特性数は５で、データは50〜950Hz までの１Hz置き
で 901個である。モード特性減衰比は全モード特性で１
％とした。
【００８３】
【表２】
【００８４】本発明と従来例における収束性を比較する
為に非線形項パラメータωdrの初期値を変えて繰り返し
計算を行った。モード特性減衰比の初期値は正解値の1
/ 10である0.1 ％である。
為に非線形項パラメータωdrの初期値を変えて繰り返し
計算を行った。モード特性減衰比の初期値は正解値の1
/ 10である0.1 ％である。
【００８５】これらの結果を結果を図４〜図６に示す。
縦軸は基準化された残差Ｚであり、次式で表される。
縦軸は基準化された残差Ｚであり、次式で表される。
【００８６】
【数１８】
【００８７】これらの図の内、図４は式(12)の左辺第
２項の２次微分項を省略した場合に対応し、図５はこ
の２次微分項を省略しない場合に対応し、そして、図
６は２次微分項を省略せず且つ式(16)による側面制約を
考慮した場合の特性をそれぞれ示しており、この順に収
束性が改善されて行き、図６の場合が従来例から大幅に
改善されていることが判る。
２項の２次微分項を省略した場合に対応し、図５はこ
の２次微分項を省略しない場合に対応し、そして、図
６は２次微分項を省略せず且つ式(16)による側面制約を
考慮した場合の特性をそれぞれ示しており、この順に収
束性が改善されて行き、図６の場合が従来例から大幅に
改善されていることが判る。
【図１】本発明及び従来例に共通な振動特性解析装置を
示したブロック図である。
示したブロック図である。
【図２】本発明に係る振動特性解析装置における演算装
置に格納され且つ実行される演算アルゴリズムを示した
フローチャート図である。
置に格納され且つ実行される演算アルゴリズムを示した
フローチャート図である。
【図３】本発明と従来例の演算時間を比較するためのグ
ラフ図である。
ラフ図である。
【図４】本発明に係る振動特性解析装置の収束性（その
１）を示したグラフ図である。
１）を示したグラフ図である。
【図５】本発明に係る振動特性解析装置の収束性（その
２）を示したグラフ図である。
２）を示したグラフ図である。
【図６】本発明に係る振動特性解析装置の収束性（その
３）を示したグラフ図である。
３）を示したグラフ図である。
【図７】被試験物を加振したときのコンプライアンス
（応答関数）の測定結果と理論値を周波数に対して示し
たグラフ図である。
（応答関数）の測定結果と理論値を周波数に対して示し
たグラフ図である。
【図８】従来例の収束性を示したグラフ図である。
１ 被試験物 ２ インパルスハンマー ３ ロードセル ４ 加速度ピックアップ ５ 演算装置 図中、同一符号は同一又は相当部分を示す。
───────────────────────────────────────────────────── フロントページの続き (58)調査した分野(Int.Cl.7，ＤＢ名) G01H 13/00 G01M 7/08 ＪＩＣＳＴファイル（ＪＯＩＳ)
Claims (3)
- 【請求項１】 被試験物を加振するためのインパルスハ
ンマーに取り付けられた力検出用ロードセルと、該被試
験物の任意の場所に取り付けられて該加振力の３次元方
向加速度を検出するためのピックアップと、該ピックア
ップからの実測データと該ロードセルの出力信号とを受
けて周波数応答関数を求めモード特性パラメータを計算
することにより該被試験物のモード特性を同定する演算
装置と、を備えた振動特性解析装置において、 該演算装置が、該実測データと該パラメータによる理論
値との残差が最小になるまで該パラメータを繰り返し変
えて行く偏分反復法を線形項と非線形項とに分けて計算
すると共に該線形項及び非線形項が変化しても該残差が
発散しない条件を設けることにより該モード特性を同定
することを特徴とした振動特性解析装置。 - 【請求項２】 該演算装置が、該発散しない条件が該残
差に対する該線形項及び非線形項による偏微分がゼロと
なることであり、このときの２次偏微分項を省略したこ
とを特徴とする請求項１に記載の振動特性解析装置。 - 【請求項３】 該演算装置が、該発散しない条件として
該偏分反復法を実行するときの該非線形項の修正ベクト
ルに対して閾値以内の縮小条件を付加することを特徴と
した請求項１又は２に記載の振動特性解析装置。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP14008194A JP3335765B2 (ja) | 1994-06-22 | 1994-06-22 | 振動特性解析装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP14008194A JP3335765B2 (ja) | 1994-06-22 | 1994-06-22 | 振動特性解析装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JPH085450A JPH085450A (ja) | 1996-01-12 |
| JP3335765B2 true JP3335765B2 (ja) | 2002-10-21 |
Family
ID=15260525
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP14008194A Expired - Fee Related JP3335765B2 (ja) | 1994-06-22 | 1994-06-22 | 振動特性解析装置 |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP3335765B2 (ja) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP3394989B2 (ja) * | 1999-11-01 | 2003-04-07 | 独立行政法人航空宇宙技術研究所 | 数値計算方法および数値計算装置、並びに数値計算用プログラムを記録した記録媒体 |
| CN106768784A (zh) * | 2017-01-22 | 2017-05-31 | 郑州云海信息技术有限公司 | 一种服务器检测方法及装置 |
-
1994
- 1994-06-22 JP JP14008194A patent/JP3335765B2/ja not_active Expired - Fee Related
Also Published As
| Publication number | Publication date |
|---|---|
| JPH085450A (ja) | 1996-01-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Richardson et al. | Global curve fitting of frequency response measurements using the rational fraction polynomial method | |
| US6763311B2 (en) | Shaking test apparatus and method for structures | |
| Yang et al. | Joint structural parameter identification using a subset of frequency response function measurements | |
| JPH0510846A (ja) | 構造物の振動試験装置、振動試験方法及び振動応答解析装置 | |
| JP3837099B2 (ja) | 構造物の損傷推定システムおよびプログラム | |
| Ashory | High quality modal testing methods | |
| Li et al. | Hybrid approach for damage detection in flexible structures | |
| JPH11118661A (ja) | 振動特性解析装置 | |
| JP3335765B2 (ja) | 振動特性解析装置 | |
| JP2003322585A (ja) | 常時微動計測に基づく建物の健全性診断法 | |
| JP3599882B2 (ja) | 振動特性解析装置 | |
| CN108573084B (zh) | 环境振动试验方法及*** | |
| JP3640583B2 (ja) | 地震応答解析方法 | |
| CN115859742A (zh) | 基于粒子群算法的有限元模型修正方法及其*** | |
| JPH05209805A (ja) | ばね・質点系のパラメータ同定装置およびその方法 | |
| JPH1130566A (ja) | 振動特性解析装置 | |
| JP3601907B2 (ja) | 振動特性解析装置 | |
| JP2800911B2 (ja) | 制御用の地震動の強度測定方法 | |
| JPH08297069A (ja) | 台上振動応力評価装置 | |
| Ruotolo et al. | A global smoothing technique for FRF data fitting | |
| JPH11281522A (ja) | 振動特性解析方法及び装置 | |
| Nagesh et al. | Finite element model updating of linear dynamic systems using a hybrid static and dynamic testing technique | |
| Formenti et al. | Modal Modeling and Structural Dynamics Modification | |
| Richardson et al. | Structural dynamics modification and modal modeling | |
| Echaniz Granado | Model calibration of a vehicle tailgate using frequency response functions |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20020723 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20080802 Year of fee payment: 6 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20090802 Year of fee payment: 7 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20090802 Year of fee payment: 7 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100802 Year of fee payment: 8 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110802 Year of fee payment: 9 |
|
| LAPS | Cancellation because of no payment of annual fees |