(57)【特許請求の範囲】
【請求項１】
多自由度系における剛性行列の対称部及び質量行列と、剛性行列の非対称部及び減衰行列の少なくとも一方と、実部が減衰率を表し且つ虚部が減衰固有角振動数を表す複素固有値及び複素固有モードと、から構成される特性方程式をコンピュータを含む計算装置で複素固有値解析し、該解析で算出された前記特性方程式を満たす前記複素固有値及び前記複素固有モードを用いて解析対象物の不安定振動を評価する振動評価方法であって、
前記計算装置が前記複素固有値解析により前記特性方程式を満たす前記複素固有値及び前記複素固有モードを算出する際、
前記特性方程式における前記剛性行列の非対称部の摩擦に関連する成分及び前記減衰行列の減衰に関連する成分を前記多自由度系に作用する外力とみなし、該外力に含まれる複素固有モードが前記多自由度系における前記剛性行列の非対称部及び前記減衰行列を省略した運動方程式で算出される実固有モードの重ね合わせで近似可能と仮定し、該仮定に従って近似された複素固有モードを含む外力及び高次の剰余剛性の前記不安定振動に対する影響を表す剰余行列を含んで前記特性方程式を満たす前記複素固有モードを表し、
前記剰余行列を算出する計算周波数を、前記多自由度系における固有角振動数のうち、前記解析の目標次数における前記複素固有値の前記減衰固有角振動数よりも小さく且つ前記減衰固有角振動数に最も近い第１固有角振動数と、前記減衰固有角振動数よりも大きな第２固有角振動数と、の中間の周波数に設定し、
前記剛性行列の非対称部、前記実固有モード及び前記剰余行列の乗算、及び、前記減衰行列、前記実固有モード及び前記剰余行列の乗算で表される行列の列成分を表す剰余ベクトルのうち、数値計算誤差の大きさよりも大きな剰余ベクトルを抽出して、
前記特性方程式を満たす前記複素固有値及び前記複素固有モードが算出され、
前記計算装置が、
前記算出された複素固有モードに対応する振動モードにおいて前記算出された複素固有値を構成する前記減衰率と前記減衰固有角振動数との比で表される減衰比を用いて前記解析対象物における不安定振動を評価する、
ように構成された、振動評価方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、複素固有値解析で得られた複素固有値及び複素固有モードを用いて不安定振動を評価する振動評価方法に関する。
【背景技術】
【０００２】
従来から、複素固有値解析で得られた複素固有値及び複素固有モードを用いて、機械、機器及び構造物等の解析対象物に生じる不安定振動を評価することが行われている。例えば、下記特許文献１に開示された部品間接触最適化設計方法では、複素固有値解析により複素固有値及び複素固有モードを算出し、得られた複素固有値及び複素固有モードを用いて不安定振動である自励振動の発生有無を評価するようになっている。従来の部品間接触最適化設計方法では、複素固有値解析にモード縮退を行う不拘束モード法を採用することで複素固有値解析における解析自由度を低減させて、複素固有値及び複素固有モードを演算するようになっている。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開２００６−３５０４７５号公報
【発明の概要】
【０００４】
ところで、上記従来の部品間接触最適化設計方法では、不拘束モード法（不拘束モード合成法）を用いて単に解析モデル規模の縮小を目的とした近似が行われるので、求められた複素固有値及び複素固有モードの計算誤差が大きくなるという問題がある。以下、この理由を具体的に説明する。
【０００５】
ディスクブレーキ装置におけるディスクロータ及びブレーキパッド間で自励振動として発生するブレーキ鳴きを解析する場合、その特性方程式は下記式５０で表される。但し、下記式５０において、Ｋはディスクロータ及びブレーキパッドからなる多自由度系（以下、多自由度系を単に「系」とも称呼する。）の剛性行列の対称部を表し、Ｍは系の質量行列を表し、Ｋμは剛性行列のうちで系における摩擦の影響を表す非対称部を表す。
【数１】
0006308117.tif 000002
【０００６】
前記式５０で表された特性方程式には剛性行列の摩擦の影響を表す非対称部Ｋμが含まれているので、前記式５０の特性方程式を満たす固有値及び固有モードは複素数になる。この場合、前記式５０におけるλは複素固有値を表し、ψは複素固有モードを表す。複素固有値λは下記数５１で表される。
【数２】
0006308117.tif 000003
但し、前記式５１中で表された複素固有値λの実部σは系の減衰率を表し、複素固有値λの虚部ωdは系の減衰固有角振動数を表す。
【０００７】
不拘束モード合成法では、複素固有値λ及び素固有モードψを算出する場合、剛性行列の非対称部Ｋμを省略した下記式５２で表される運動方程式に着目し、実固有値Ω2と実固有モードφを算出する。
【数３】
0006308117.tif 000004
前記式５０で表された複素固有モードψは、下記式５３で表されるように、前記式５２で算出された実固有モードφiを重ね合わせることにより近似される。
【数４】
0006308117.tif 000005
但し、前記式５３におけるφiはｉ次の実固有モードを表し、ｎは実固有モード数（以下、「採用モード数」と称呼する。）を表し、ξはモード寄与ベクトルを表す。
【０００８】
複素固有モードψが前記式５３で表されたように近似されると、前記式５０で表された特性方程式は下記式５４で表されるようにモード座標上の特性方程式に変換される。
【数５】
0006308117.tif 000006
【０００９】
ところで、前記式５４で表されたモード座標上の特性方程式は、前記式５０で表された特性方程式がｉ次の実固有モードφiの重ね合わせによりモード縮退され、その結果、モード規模が縮小された特性方程式である。これにより、複素固有値λ及び複素固有モードψを算出する際の計算負荷が低減されるので、例えば、計算時間が低減される等、計算効率は確かに向上する。
【００１０】
しかし、不安定振動を複素固有値解析して評価する場合、ｎ次よりも高次のモードの振動に対する影響を無視できない場合がある。この点に関し、上記従来の不拘束モード合成法では、前記式５３で表されたように複素固有モードψをｎ次までの実固有モードφiの重ね合わせで近似するので、ｎ次よりも高次のモードの振動に対する影響が考慮されない。尚、ｎ次よりも高次のモードは、一般に、剰余剛性と呼ばれる。従って、上記従来の不拘束モード合成法では、剰余剛性の不安定振動に対する影響を反映させた計算を行うことができず、その結果、計算誤差が大きくなるので得られた複素固有値λ及び複素固有モードψの計算精度が悪化する。
【００１１】
尚、上記従来の不拘束モード合成法でも、論理上は複素固有モードψをｎ次の実固有モードφi以上の実固有モードの重ね合わせで近似することは可能であるかもしれない。しかしながら、この場合には、計算対象となるモード規模が拡大するので計算負荷が著しく増加し、その結果、膨大な計算時間が必要となる。従って、上記従来の不拘束モード合成法のままでモード規模を増加させることは現実的ではない。
【００１２】
本発明は、上記課題を解決するためになされたものであり、その目的の一つは、複素固有値解析において剰余剛性を考慮し、複素固有値及び複素固有モードを高い精度で算出して不安定振動を評価することができる振動評価方法を提供することにある。
【００１３】
上記目的を達成するための本発明の振動評価方法は、多自由度系（２０）における剛性行列の対称部（Ｋ）及び質量行列（Ｍ）と、剛性行列の非対称部（Ｋμ）及び減衰行列（Ｄ）の少なくとも一方と、実部が減衰率（σ）を表し且つ虚部が減衰固有角振動数（ωd）を表す複素固有値（λ）及び複素固有モード（ψ）と、から構成される特性方程式（式１）をコンピュータ（１２）を含む計算装置（１０）で複素固有値解析し、該解析で算出された前記特性方程式を満たす前記複素固有値（λ）及び前記複素固有モード（ψ）を用いて解析対象物（２０）の不安定振動を評価するものである。
【００１４】
本発明の振動評価方法では、
前記計算装置（１０）が前記複素固有値解析により前記特性方程式（式１）を満たす前記複素固有値（λ）及び前記複素固有モード（ψ）を算出する際、
前記特性方程式（式１）における前記剛性行列の非対称部（Ｋμ）の摩擦に関連する成分及び前記減衰行列（Ｄ）の減衰に関連する成分を前記多自由度系に作用する外力（fμ，fD）とみなし（式２）、該外力（fμ，fD）に含まれる複素固有モード（ψ）（式３，４）が、前記多自由度系における前記剛性行列の非対称部（Ｋμ）及び前記減衰行列（Ｄ）を省略した運動方程式（式５２）で算出される実固有モード（Φ）の重ね合わせで近似可能と仮定し（式６）、該仮定に従って近似された複素固有モード（ψ）を含む外力（fμ，fD）及び高次の剰余剛性の前記不安定振動に対する影響を表す剰余行列（Gc）を含んで前記特性方程式（式１）を満たす前記複素固有モード（ψ）を表し（式５，７）、
前記剰余行列（Gc）を算出する計算周波数（ωc）を、前記多自由度系における固有角振動数（Ωi）のうち、解析の目標次数における前記複素固有値（λ）の前記減衰固有角振動数（ωd）よりも小さく且つ前記減衰固有角振動数（ωd）に最も近い第１固有角振動数（Ωj）と、前記減衰固有角振動数（ωd）よりも大きな第２固有角振動数（Ωj+1）と、の中間の周波数に設定し、
前記剛性行列の非対称部（Ｋμ）、前記実固有モード（Φ）及び前記剰余行列（Gc）の乗算（式８）、及び、前記減衰行列（Ｄ）、前記実固有モード（Φ）及び前記剰余行列（Gc）の乗算（式９）で表される行列（Rμ，RD）の列成分を表す剰余ベクトル（Rμ,i，RD,i）のうち、数値計算誤差の大きさよりも大きな剰余ベクトル（Rμ,i，RD,i）を抽出（式２７）して、
前記特性方程式（式１）を満たす前記複素固有値（λ）及び前記複素固有モード（ψ）が算出され、
前記計算装置（１０）が、
前記算出された複素固有モード（ψ）に対応する振動モードにおいて前記算出された複素固有値（λ）を構成する前記減衰率（σ）と前記減衰固有角振動数（ωd）との比で表される減衰比（σ／ωd）を用いて前記解析対象物（２０）における不安定振動を評価する（ステップＳ１８）、ように構成される。
【００１５】
本発明の振動評価方法によれば、摩擦力及び減衰力を多自由度系に作用する外力とみなし、これらの外力に含まれる複素固有モードを実固有モードの重ね合わせで近似することができる。これらにより、近似された複素固有モードを含む外力及び高次の剰余剛性の不安定振動に対する影響を表す剰余行列を含んで、複素固有値解析における特性方程式を満たす複素固有モードを表すことができる。従って、例えば、ブレーキ鳴き解析のような複素固有値解析において剰余剛性の影響を考慮することができる。これにより、複素固有値解析の解である複素固有値及び複素固有モードの計算精度を向上させることができる。
【００１６】
加えて、第１固有角振動数と第２固有角振動数との中間の周波数となるように剰余行列を計算する計算周波数を設定することができる。これにより、数値計算誤差を抑制しつつ精度の良い剰余行列を計算することができる。従って、数値計算誤差の影響が小さくなり、前記特性方程式を満たす複素固有値及び複素固有モードを精度良く算出することができる。
【００１７】
更に、剰余ベクトルのうち、数値計算誤差よりも大きな剰余ベクトルを抽出することができる。これにより、数値計算誤差に埋もれてしまうことなく、剰余ベクトルを構成する剰余行列を複素固有値及び複素固有モードの算出に適切に反映させることができる。その結果、数値計算誤差の影響が小さくなり、前記特性方程式を満たす複素固有値及び複素固有モードを精度よく算出することができる。
【００１８】
このように、複素固有値解析の解である複素固有値及び複素固有モードを高次の剰余剛性を考慮して精度良く算出することができるので、これらの複素固有値及び複素固有モードを用いて不安定振動を精度良く評価することができる。従って、不安定振動である自励振動により発生するブレーキ鳴きのような現象の発生有無を正確に判断することができる。
【００１９】
尚、上記説明においては、発明の理解を助けるために、実施形態に対応する発明の構成に対して、実施形態で用いた符号及び式番号を括弧書きで添えているが、発明の各構成要件は上記符号によって規定される実施形態に限定されるものではない。
【図面の簡単な説明】
【００２０】
【図１】本発明の実施形態に係る振動評価方法が実行される計算装置を示すブロック図である。
【図２】剰余モード合成法で用いられる剰余ベクトルの抽出を説明するためのグラフである。
【図３】図３（Ａ）はブレーキ鳴き解析を行う場合に想定するブレーキモデルにおいてブレーキパッド及びディスクロータに作用する力の関係を示す概略図、（Ｂ）はばねがブレーキパッド及びディスクロータから受ける力の関係を示す概略図である。
【図４】ブレーキ鳴き解析プログラムを示すフローチャートである。
【図５】直接法で算出された複素固有値を基準とした場合、従来の不拘束モード合成法で算出された複素固有値の計算誤差と、剰余モード合成法で算出された複素固有値の計算誤差と、の差を説明するための表である。
【発明を実施するための形態】
【００２１】
以下、本発明の実施形態に係る振動評価方法について図面を参照しながら説明する。この実施形態に係る振動評価方法は、図１に示した計算装置１０を利用して実行される。
【００２２】
計算装置１０は、有限要素法（ＦＥＭ）等における複素固有値解析を用いて解析対象物（機械、機器及び構造物等）に発生する不安定振動を評価する。計算装置１０は、入力部１１、演算部１２、記憶部１３及び出力部１４を含んで構成される。入力部１１は、解析対象物の形状データ、材料定数、境界条件及び最適化検討条件等の各種データを入力する。
【００２３】
演算部１２は、ＣＰＵ、ＲＯＭ、ＲＡＭ等からなるマイクロコンピュータを主要構成部品として構成される。演算部１２は、各種インターフェースを有しており、このインターフェースを介して入力部１１、記憶部１３及び出力部１４と信号を入出力可能に接続されている。
【００２４】
演算部１２は、入力部１１を介して入力された各種データを用いて解析対象物における複素固有値解析の解である複素固有値及び複素固有モードを算出する。演算部１２は、算出した複素固有値及び複素固有モードを用いて解析対象物における不安定振動の発生有無を評価する。尚、本実施形態では、後述するように、演算部１２は、不安定振動である自励振動の解析、例えば、ブレーキ鳴き解析を行い、振動モード毎に自励振動の発生有無を評価する。
【００２５】
記憶部１３は、入力部１１で入力された各種データを演算部１２を介して記憶するとともに、記憶した各種データを演算部１２に供給する。記憶部１３は、演算部１２によって算出された複素固有値及び複素固有モードと、不安定振動の発生有無の評価結果と、を関連付けて記憶する。出力部１４は、演算部１２で算出された、或いは、記憶部１３に記憶された複素固有値及び複素固有モードと、不安定振動の発生有無の評価結果と、を出力する。
【００２６】
（剰余剛性を考慮した不拘束モード合成法）
計算装置１０は、複素固有値解析を行なって複素固有値及び複素固有モードを算出する際、上記従来の不拘束モード合成法では考慮されていない高次のモードである、剰余剛性の影響をも考慮した不拘束モード合成法を利用する。尚、以下の説明においては、剰余剛性の影響を考慮した不拘束モード合成法を、単に「剰余モード合成法」と称呼する。剰余モード合成法を利用することにより、複素固有値及び複素固有モードの計算精度を従来の不拘束モード合成法に比べて大幅に改善することができる。以下、この剰余モード合成法を説明する。
【００２７】
剰余モード合成法では、下記式１により表される特性方程式を考える。
【数６】
0006308117.tif 000007
前記式１は、前述した不拘束モード合成法の前記式５０と比べて、多自由度系に生じる減衰及びジャイロ効果等に関連する減衰行列Ｄに関する項が追加されている点で異なる。従って、前記式１で表された特性方程式にも、前記式５０と同様に、剛性行列の非対称部Ｋμが含まれているので、固有値及び固有モードが複素数になる。
【００２８】
前記式１において、非対称部Ｋμの摩擦に関連する成分及び減衰行列Ｄの減衰に関連する成分を右辺に移項すると、下記式２が得られる。
【数７】
0006308117.tif 000008
前記式２の右辺におけるfμは外力として多自由度系に作用する摩擦力であり、複素固有モードψを用いた下記式３で表される。前記式２の右辺におけるfDは外力として多自由度系に作用する減衰力であり、複素固有モードψを用いた下記式４で表される。これにより、前記式２は、その左辺が対称行列で表される不減衰系（多自由度系）であり、この不減衰系（多自由度系）に非対称行列で表される摩擦力及び減衰力が外力として作用することを表している。
【数８】
0006308117.tif 000009
【数９】
0006308117.tif 000010
【００２９】
下記式５で表されるように、剰余モード合成法でも、複素固有モードψは、不減衰系の１次からｎ次までの、換言すれば、採用モード数ｎまでの実固有モードφi（実固有モードΦ）の重ね合わせを含んでいる（前述した不拘束モード合成法の前記式５３を参照）。加えて、剰余モード合成法では、前述した不拘束モード合成法において採用されなかったｎ＋１次以上の剰余剛性の影響を表す剰余行列Gcを導入することで、複素固有モードψの実固有モードΦを用いた近似精度の改善を図る。
【数１０】
0006308117.tif 000011
前記式５の右辺の摩擦力fμ及び減衰力fDには、前記式３，４で示したように複素固有モードψが含まれているので、前記式５をそのまま利用することができない。加えて、剰余行列Gcは、物理空間における自由度（物理自由度）がＮ×Ｎの密行列である。このため、剰余行列Gcを物理座標のまま保持することは現実的ではない。
【００３０】
そこで、摩擦力fμ及び減衰力fDを表す前記式３，４における複素固有モードψも、不拘束モード合成法の場合（前記式５３を参照）と同様に、下記式６に従って実固有モードΦ（実固有モードφi）の重ね合わせで表現できると仮定する。この場合、複素固有モードψを近似するときの実固有モードΦ（実固有モードφi）の寄与をηとする。
【数１１】
0006308117.tif 000012
これにより、摩擦力fμ及び減衰力fDを求める場合には、前記式６に従って複素固有モードψを実固有モードΦ（実固有モードφi）の重ね合わせで近似する。一方で、系全体の複素固有モードψを求める場合には、前記式５に従って実固有モードΦ（実固有モードφi）及び剰余行列Gcで近似する。
【００３１】
ここで、摩擦力fμを求める場合の実固有モードφiの寄与をημとし、減衰力fDを求める場合の実固有モードφiの寄与をηDとする。そして、前記式６を前記式３，４，５に代入して整理すると、下記式７，８，９が得られる。
【数１２】
0006308117.tif 000013
【数１３】
0006308117.tif 000014
【数１４】
0006308117.tif 000015
但し、前記式８，９で表されたRμ，RDは系の物理自由度Ｎ×採用モード数ｎの行列であり、それぞれ、-KμΦ，-DΦの力が作用したときの剰余行列Gcの変形（変位）を表す。尚、以下の説明のおいてはRμ，RDの各列を剰余ベクトルと称呼し、ｉ行の成分をｉ次の剰余ベクトルRμ,i，RD,i（或いは、単に、剰余ベクトルRi）と称呼する。前記式７によれば、寄与ημ，ηDは行列Rμ，RD（即ち、剰余ベクトルRμ,i，RD,i）の複素固有モードψに対する寄与を表していることが理解できる。
【００３２】
剰余モード合成法では、摩擦力fμ及び減衰力fDを実固有モードΦ（実固有モードφi）の重ね合わせで近似する場合、実固有モードΦ（実固有モードφi）の寄与を新たな未知数ημ，ηDとする。この場合、前記式３〜７で表される各定義に基づけば、物理的にはモード寄与ベクトルξ、寄与ημ及び寄与ηDの間には、ξ＝ημ＝ηDとなる従属関係が成立する。一方で、近似を行う場合には、モード寄与ベクトルξ、寄与ημ及び寄与ηDを独立に扱うことで近似精度が向上する。何故ならば、前述した従来の不拘束モード合成法のように、モード寄与ベクトルξが系全体の複素固有モードψを近似するのに対して、寄与ημは摩擦力fμが作用している点のみを近似し、寄与ηDは減衰力fDが作用している点のみを近似する（前記式７を参照）。これにより、摩擦力fμ及び減衰力fDの近似精度が改善され、その結果として複素固有モードψの近似精度も改善されるためである。
【００３３】
複素固有モードψを前記式７に従って実固有モードΦ（実固有モードφi）と行列Rμ，RD（剰余ベクトルRμ,i，RD,i）で表すことにより、前記式１で表された特性方程式は、下記式１０で表されるモード座標上の特性方程式に変換される。
【数１５】
0006308117.tif 000016
但し、前記式１０中のkはモード座標上における剛性行列の対称部であり下記式１１で表される。前記式１０中のkμはモード座標上における剛性行列の非対称部であり下記式１２で表される。前記式１０中のmはモード座標上における質量行列であり下記式１３で表される。更に、前記式１０中のdはモード座標上における減衰行列であり下記式１４で表される。
【００３４】
【数１６】
0006308117.tif 000017
【数１７】
0006308117.tif 000018
【数１８】
0006308117.tif 000019
【数１９】
0006308117.tif 000020
【００３５】
前記式１０で表された特性方程式におけるモード寄与ベクトルξ、寄与ημ及び寄与ηDが一意に決定されるためには、実固有モードΦと行列Rμ，RD（剰余ベクトルRμ,i，RD,i）の各列が線形独立でなければならない。この場合、実固有モードΦは互いに線形独立であり、且つ、剰余行列Gcは下記式１５で定義されるように１次からｎ次までの実固有モード成分φiを含まない。よって、剰余行列Gcを用いて表される行列Rμ，RDでは（前記式８，９を参照）、剰余ベクトルRμ,i，RD,iの列のうちで線形独立な列は実固有モード成分φiとも線形独立となる。尚、線形独立な剰余ベクトルRμ,i，RD,iの列を取り出す簡易的な方法は、例えば、周知のＱＲ分解等により抽出が可能となる。
【数２０】
0006308117.tif 000021
【００３６】
但し、前記式１５中のωcは、剰余行列Gcを計算する周波数（以下、この周波数を「剰余計算周波数」と称呼する。）を表す。更に、前記式１５中のΩiは、多自由度系におけるｉ次の固有角振動数を表す。前記式１５は、周波数ωの関数である剰余行列G(ω)を一つの周波数ωcで計算した定数行列である剰余行列Gcで近似することを示している。
【００３７】
更に、剰余モード合成法における特性方程式は、前記式１０〜１４からも明らかなように、３ｎ×３ｎの自由度であるので、採用モード数ｎが大きい場合には物理自由度Ｎ＜３ｎになる可能性がある。しかし、上述したように、線形独立な剰余ベクトルRμ,i，RD,iを選択することにより、採用モード数ｎが大きい場合でもモード座標上の特定方程式の自由度は物理自由度Ｎを超えない。
【００３８】
ここで、剰余計算周波数ωcと剰余行列Gcの近似精度について説明する。前記式１５によれば、定数行列の剰余行列Gcは剰余計算周波数ωcの値により決定される。従って、剰余計算周波数ωcの値により、周波数ωの関数である剰余行列G(ω)の剰余行列Gcによる近似精度が変化する。この場合、剰余計算周波数ωcの値が目標とする（求めようとする）周波数ωに近くなるほど近似精度は改善されると言われており、特に、剰余計算周波数ωcが周波数ωと等しいとき、剰余行列Gc＝剰余行列G(ω)が成立する。
【００３９】
前記式１５における剰余計算周波数ωcに代えて複素固有値λを入力すると、下記式１６が成立する。尚、複素固有値λは前記式５１からも明らかなように減衰固有角振動数ωdを含む周波数である。この場合、剰余モード合成法では摩擦力fμに加えて減衰力fDをも外力とするので、下記式１６左辺の特性行列の逆行列は下記式１６の第２行のように近似されると考えることができる。
【数２１】
0006308117.tif 000022
【００４０】
更に、剰余行列G(ω)において、ｎ＋１次以上の高次の実固有モード成分φjの影響は、下記式１７で表すことができる。
【数２２】
0006308117.tif 000023
この場合、前記式１７に周波数ωとして複素固有値λ（周波数）を代入した剰余行列G(λ)と、前記式１７に周波数ωとして剰余計算周波数ωcを代入した剰余行列Gcと、を考える。ｎ＋１次以上の高次において、前記式１６で表されたように剰余行列G(λ)を剰余行列Gcで近似したときの誤差行列E(λ)は、下記式１８で表される。尚、この場合、剰余計算周波数ωcの値は実数であるとする。
【数２３】
0006308117.tif 000024
【００４１】
ここで、下記式１９を満たすように剰余計算周波数ωc、及び、ｎ＋１次以上の高次の固有角振動数Ωjを選択すると、前記式１８中の成分は、複素固有値λを構成する任意の減衰比σ及び任意の減衰固有角振動数ωdに対して下記式２０を満たすようになる。
【数２４】
0006308117.tif 000025
【数２５】
0006308117.tif 000026
【００４２】
前記式２０が成立する場合には、前記式１８で表される誤差行列E(λ)が発散することはない。更に、剰余計算周波数ωcは、複素固有値λ（周波数）に近い程、且つ、高次の固有角振動数Ωjから離れる程、誤差行列E(λ)の値が小さくなることが理解できる。尚、本実施形態の剰余モード合成法では、計算装置１０の演算部１２の計算負荷を低減する観点から剰余計算周波数ωcを実数とし、剰余計算周波数ωcを減衰固有角振動数ωdに近づけるように計算することで近似精度の改善を図る。
【００４３】
ところで、上述した剰余行列Gc及び行列Rμ，RD（剰余ベクトルRμ,i，RD,i）を数値計算で求める場合、
（１）剰余計算周波数ωcの決定、及び、
（２）剰余ベクトルRμ,i，RD,iに含まれる数値計算誤差
の２点について留意する必要がある。以下、この２つの留意点を順に説明する。
【００４４】
（１）剰余計算周波数ωcの決定について
前述したように、剰余計算周波数ωcを解析の目標次数としてｉ次の減衰固有角振動数ωdに近づける程、誤差行列E(λ)の値が小さくなる。その結果、剰余行列Gcによる近似精度が向上するので、前記式７〜９から明らかなように複素固有モードψを表現する剰余ベクトルRμ，RDの近似精度も向上する。従って、不安定振動を生じさせる複素固有モードψが特定できている場合、剰余計算周波数ωcと減衰固有角振動数ωdとを等しくすればよいが、数値計算の観点では以下の点を考慮する必要がある。
【００４５】
剰余行列Gcの計算では、前記式１５を変形して得られる下記式２１で表されるように、不減衰系の応答計算が必要となる。
【数２６】
0006308117.tif 000027
この場合、剰余計算周波数ωcが不減衰系の固有角振動数近傍Ωiになると、前記式２１右辺第１項及び第２項が無限大に近づく場合があり、第１項及び第２項の差分で定義された剰余行列Gcが数値計算誤差に埋もれてしまう可能性がある。よって、図２に概略的に示したように、ｉ次（解析の目標次数）の減衰固有角振動数ωdに近い固有角振動数Ωj（第１固有角振動数）と、固有角振動数Ωj+1（第１固有角振動数よりも大きな第２固有角振動数）と、の中間の周波数を剰余計算周波数ωcとする。具体的には、減衰固有角振動数ωdに近い固有角振動数Ωjと、固有角振動数Ωj+1と、を用いて表される下記式２２に従い、剰余計算周波数ωcは決定される。これにより、数値計算誤差を抑制しつつ精度の良い剰余行列Gcを計算することができる。
【数２７】
0006308117.tif 000028
【００４６】
（２）剰余ベクトルRμ,i，RD,iに含まれる数値計算誤差について
前述したように、行列Rμ，RD（剰余ベクトルRμ,i，RD,i）は、それぞれ、-KμΦ，-DΦの力が作用したときの剰余行列Gcの変形（変位）を表す。このため、剰余ベクトルRμ,i，RD,iは、直接法として例えばGaussの消去法等で求めた変形（変位）と、モードの重ね合わせで求めた変形（変位）と、の差分で求められる。尚、数値計算誤差については、剰余ベクトルRμ,i，RD,i共に同様に生じるので、以下では摩擦力fμによる剰余ベクトルRμ,iを例示して説明する。
【００４７】
摩擦力fμによるi次の剰余ベクトルRμ,iは、下記式２３で表されるように、直接法で求めた変形xd,iと、モードの重ね合わせで求めた変形xm,iと、の差分で求められる。
【数２８】
0006308117.tif 000029
【００４８】
ここで、前記式２３に対して前記式２１を代入することにより、変形xd,iは下記式２４で表され、変形xm,iは下記式２５で表される。更に、前記式３及び前記式６から、摩擦力は下記式２６で表すことができる。
【数２９】
0006308117.tif 000030
【数３０】
0006308117.tif 000031
【数３１】
0006308117.tif 000032
【００４９】
図２に示したように、直接法で求めた変形xd,i（実線）とモードの重ね合わせで求めた変形xm,i（破線）との差分が剰余ベクトルRμ,iとなる。前述したように、剰余計算周波数ωcを前記式２２に従って決定した場合であっても、前記式２６に従って計算される摩擦力fμ,iの成分によっては前記式２３で求まる剰余ベクトルRμ,iの大きさ数値計算誤差に埋もれてしまう場合がある。
【００５０】
このため、計算した剰余ベクトルRμ,iのうち、数値計算誤差よりも大きく、数値計算誤差に埋もれていない列だけを抽出する。列の抽出は、下記式２７で表されるように、直接法で求めた変形xd,i（変位xd,i）によって生じる運動エネルギーと、ｉ次の剰余ベクトルRμ,i（剰余変形Rμ,i）によって生じる運動エネルギーと、の比を判定して行う。具体的には、前記運動エネルギー比が一定値ε以上（例えば、ε＝１×１０−４）になる剰余ベクトルRμ,iのみを採用する。このことは、変位によって生じる運動エネルギーの０．１％以上の寄与を有する剰余ベクトルRμ,iを採用することを意味している。
【数３２】
0006308117.tif 000033
【００５１】
＜剰余モード合成法のブレーキ鳴き解析への適用＞
次に、複素固有値解析において前述した剰余モード合成法を採用して複素固有値及び複素固有モードを算出し、この複素固有値及び複素固有モードを用いて不安定振動である自励振動の発生有無を評価する場合を具体的に説明する。以下、自励振動により生じるブレーキ鳴き解析において、前述した剰余モード合成法を利用して自励振動の発生有無を評価する場合を説明する。
【００５２】
ブレーキ鳴き解析においては、図３（Ａ）に示したディスクブレーキ装置のモデル２０（以下、単に「ブレーキモデル２０」と称呼する。）に基づき、運動方程式を導出する。ブレーキモデル２０は、ディスクロータ２１とブレーキパッド２２の接触力を表現するためにディスクロータ２１とブレーキパッド２２との間にばね定数ｋのばね２３を配置する。この場合、ばね２３のディスクロータ２１側の節点を２３１とし、ブレーキパッド２２側の節点を２３２とする。
【００５３】
加えて、制動時にディスクロータ２１とブレーキパッド２２との間に発生する押し付け力Ｆに関し、押し付け力Ｆがブレーキパッド２２からディスクロータ２１に向けて発生するように、節点２３１，２３２のｙ方向の変位関係をｙ１＞ｙ２と仮定する。更に、ディスクロータ２１はｘ方向に変位するものと仮定する。この場合、摩擦係数をμとすると、ディスクロータ２１に発生する摩擦力μＦ及びブレーキパッド２２に発生する摩擦力μＦは、図３（Ａ）に示した方向にてそれぞれ作用する。
【００５４】
ディスクロータ２１及びブレーキパッド２２に摩擦力μＦが発生している状態において、図３（Ｂ）に示したばね２３の力の釣り合いを考える。この釣り合いにおいては、ディスクロータ２１側のばね２３の節点２３１には上方向の反力が作用し、ブレーキパッド２２側のばね２３の節点２３２には下方向の反力が作用する。これらの反力は、それぞれ、下記式２８，２９で表される。
【数３３】
0006308117.tif 000034
【数３４】
0006308117.tif 000035
尚、前記式２８で表されるＦy1は節点２３１のｙ方向の力であり、前記式２９で表されるＦy2は節点２３２のｙ方向の力である。
【００５５】
更に、押し付け力Ｆにより、ディスクロータ２１側のばね２３の節点２３１に生じる摩擦力Ｆx1、及び、ブレーキパッド２２側のばね２３の節点２３２に生じる摩擦力Ｆx2は、それぞれ、下記式３０，３１で表される。
【数３５】
0006308117.tif 000036
【数３６】
0006308117.tif 000037
【００５６】
このように構成されたブレーキモデル２０について、計算装置１０の演算部１２は、図４のフローチャートで示されるブレーキ鳴き解析プログラムを実行する。具体的に、演算部１２は、ブレーキ鳴き解析プログラムの実行をステップＳ１０にて開始し、続くステップＳ１１にて、ブレーキモデル２０即ち多自由度系における剛性行列の対称部Ｋ、剛性行列の非対称部Ｋμ及び質量行列Ｍを組み立てる。この場合、演算部１２は、入力部１１を介して入力されたディスクロータ２１及びブレーキパッド２２の形状データ、材料定数及び境界条件等に基づき、例えば、有限要素法（ＦＥＭ）に従って各行列Ｋ、Ｋμ及びＭを組み立てる。これにより、ブレーキモデル２０について、慣性力を考慮した下記式３２で表される運動方程式が成立する。
【数３７】
0006308117.tif 000038
尚、前記式３２で表される運動方程式は、時間領域の変位x(t)及び外力f(t)をそれぞれ下記３３のようにおいた周波数領域の運動方程式である。
【数３８】
0006308117.tif 000039
【００５７】
前記式３２で表された運動方程式の右辺を０とすることにより、前記式５０で表された特性方程式が成立する。これにより、前記式５０で表された特性方程式を解けば、ブレーキモデル２０のように、摩擦が存在する多自由度系における固有値及び固有モードを算出することができる。
【００５８】
そこで、演算部１２は、ブレーキ鳴き解析プログラムのステップＳ１２にて、前記式５２で表される前述の運動方程式を成立させる実固有値Ω2及び実固有モードφを算出する。実固有値Ω2及び実固有モードφを算出すると、演算部１２はステップＳ１３を実行する。
【００５９】
ステップＳ１３では、演算部１２は、前記式１５に従って剰余行列G(ω)を一つの周波数ωcで計算した定数行列である剰余行列Gcで近似する。この場合、前述したように、剰余行列Gcの計算では、前記式２１で表されたように不減衰系の応答計算が必要となるので、剰余計算周波数ωc及び減衰固有角振動数ωdは前記式２２を満たすように決定される。このように決定された剰余計算周波数ωc及び減衰固有角振動数ωdを用いることにより、数値計算誤差を抑制しつつ精度のよい剰余行列Gcが計算される。
【００６０】
剰余行列Gcを計算すると、演算部１２はステップＳ１４にて剰余ベクトルを抽出する。ブレーキモデル２０においては、前記式５０で表された特性方程式が成立する。つまり、ブレーキモデル２０においては、前記式１で表された特性方程式における減衰行列Ｄに関連する項の省略された特性方程式が成立する。このため、演算部１２は、ステップＳ１４において、減衰力fDによるi次の剰余ベクトルRD,iの計算を省略し、前述したように摩擦力fμによるi次の剰余ベクトルRμ,iのみを抽出する。この場合、演算部１２は、前記式２７に従い、直接法で求めた変形xd,iによって生じる運動エネルギーとｉ次の剰余ベクトルRμ,iによって生じる運動エネルギーの比が一定値ε以上になる剰余ベクトルRμ,iのみを抽出する。
【００６１】
剰余ベクトルRμ,iを抽出すると、演算部１２はステップＳ１５にて高次のモード即ち剰余剛性を考慮した行列の縮退（モード縮退）を行う。具体的に、ブレーキモデル２０では、減衰行列Ｄに関連する項が省略される。従って、演算部１２は、前記式１１〜１３に従った下記式３４，３５，３６で表されるように、剰余ベクトルRμ,iを含む行列Rμを用いてモード座標上における剛性行列の対称部k、剛性行列の非対称部kμ及び質量行列mを決定する。
【数３９】
0006308117.tif 000040
【数４０】
0006308117.tif 000041
【数４１】
0006308117.tif 000042
【００６２】
モード座標上における剛性行列の対称部k、剛性行列の非対称部kμ及び質量行列mを決定して縮退を行うと、演算部１２は、ステップＳ１６にて、前記式１０及び前記式７に従った下記式３７及び下記式３８を解き、複素固有値λ及び複素固有モードψを算出する。
【数４２】
0006308117.tif 000043
【数４３】
0006308117.tif 000044
【００６３】
複素固有値λ及び複素固有モードψを算出すると、演算部１２は、ステップＳ１７にてブレーキ鳴きの発生有無、即ち、自励振動の発生有無を評価する。具体的に、演算部１２は、複素固有モードψに対応する振動モード（Ｌ次のモード）において、複素固有値λの実部である減衰率σと複素固有値λの虚部である減衰固有角振動数ωdとの比で表される減衰比σ／ωdが正であるか否かを判定する。
【００６４】
この判定に従い、演算部１２は、あるＬ次のモードにおいて減衰比σ／ωdが正である場合には、ディスクロータ２１及びブレーキパッド２２間で発生した振動が減衰されるので、「自励振動によるブレーキ鳴きは生じない」と評価する。一方、演算部１２は、あるＬ次のモードにおいて減衰比σ／ωdが負である場合には、ディスクロータ２１及びブレーキパッド２２間で振動が励振されるので、「自励振動によるブレーキ鳴きが生じる」と評価する。この場合には、複素固有モードψに対応して自励振動の生じる振動モード（Ｌ次のモード）が特定されるので、この振動モードが生じないように、例えば、ブレーキパッド２２の形状変更及びディスクロータ２１の剛性を大きくする等の対策を施すことができる。
【００６５】
このように、ブレーキ鳴き解析（自励振動の発生有無）を評価すると、演算部１２はステップＳ１８にてブレーキ鳴き解析プログラムの実行を終了する。
【００６６】
ところで、剰余モード合成法を用いた場合には、複素固有値λ及び複素固有モードψを算出する際に、高次のモードである剰余剛性を考慮するので、複素固有値λ及び複素固有モードψの計算精度を向上させることができる。このことに関し、節点数の少ない簡易的なモデルについて直接法で算出した複素固有値と、上述した従来の不拘束モード合成法を用いて算出した複素固有値及び剰余モード合成法を用いて算出した複素固有値と、を比較する。
【００６７】
図５に採用モード数ｎに対応して求められる採用モード周波数毎に示したように、直接法による複素固有値を基準とした場合の計算誤差（％）は、剰余モード合成法による複素固有値λに比べて従来の不拘束モード合成法による複素固有値λの方が明らかに大きい。剰余モード合成法では、剰余剛性の影響も考慮しているので、計算精度を大幅に向上させることができる。不拘束モード合成法においても採用モード数ｎ（採用モード周波数）を増加させることにより論理上は複素固有値λの計算精度を向上させることができる。しかし、この場合には、計算負荷が著しく増大するとともに膨大な計算時間が必要となり、現実的ではない。これに対して、上述した剰余モード合成法を採用した場合では、１０ＫＨｚ以下の採用モード周波数、即ち、少ない採用モード数ｎであっても十分高い計算精度で複素固有値λを算出することができる。これにより、高次のモードである剰余剛性を考慮した場合であっても、計算負荷の増加を抑えることができ、且つ、計算時間を短縮することができる。
【００６８】
以上の説明からも理解できるように、上記実施形態によれば、複素固有値解析の解である複素固有値λ及び複素固有モードψの計算精度を向上させるために、剰余モード合成法を採用することができる。剰余モード合成法を採用する場合、前記式２で表されたように摩擦力fμ及び減衰力fDが多自由度系に作用する外力とみなすことができる。更に、前記式３，４で表されたように、これらの摩擦力fμ及び減衰力fDに含まれる複素固有モードψが実固有モードΦ（実固有モードφi）の重ね合わせで近似することができる。これらにより、例えば、ブレーキ鳴き解析のような複素固有値解析において、剰余剛性を考慮する剰余モード合成法を採用することができる。
【００６９】
複素固有値解析に剰余モード合成法を採用した場合において、数値計算誤差の影響を小さくするために、前記式２２を満たすように剰余計算周波数ωcを決定することができる。これにより、数値計算誤差を抑制しつつ精度の良い剰余行列Gcを計算することができ、その結果、複素固有値λ及び複素固有モードψを精度良く算出することができる。
【００７０】
更に、複素固有値解析に剰余モード合成法を採用した場合において、数値計算誤差の影響を小さくするために、前記式２７を満たすような剰余ベクトルRiを抽出することができる。これにより、数値計算誤差に埋もれてしまうことなく、剰余ベクトルRiを構成する剰余行列Gcを複素固有値λ及び複素固有モードψの算出に適切に反映させることができる。その結果、複素固有値λ及び複素固有モードψを精度よく算出することができる。
【００７１】
本発明は、上記実施形態に限定されることなく、本発明の範囲内において種々の変更例を採用することができる。
【００７２】
例えば、上記実施形態では剰余モード合成法を採用した複素固有値解析を摩擦力fμに起因して発生する自励振動であるブレーキ鳴きの評価に用いたが、剰余モード合成法を採用した複素固有値解析は減衰力fDに起因して発生する不安定振動を評価するために用いられても良い。更に、剰余モード合成法を採用した複素固有値解析は、摩擦力fμ及び減衰力fDに起因して発生する不安定振動の評価に用いられても良い。
【符号の説明】
【００７３】
１０…計算装置、１１…入力部、１２…演算部（コンピュータ）、１３…記憶部、１４…出力部。
【図１】
0006308117.tif 000045
【図２】
0006308117.tif 000046
【図３】
0006308117.tif 000047
【図４】
0006308117.tif 000048
【図５】
0006308117.tif 000049