(57)【要約】  【課題】全体座標系上における物体の剛体変位と弾性振動とを解析する方法において、弾性自由度が大きい弾性振動において剛体変位の効果を考慮することが必要である場合であっても物体の剛体変位と弾性振動とを容易に解析する。  【解決手段】全体座標系上を運動しつつ物体を局所的に観測する局所観測座標系ＬＯＣであって、その物体の少なくとも近傍に位置して運動を行うものを全体座標系と共に用いることにより、その物体の剛体変位と弾性振動とを、剛体変位は剛体モードにより、弾性振動は弾性モードによりそれぞれ表すとともに、それら剛体モードと弾性モードとを区別しないで一律に取り扱うモード解析により、解析する。【選択図】図３
2005091370.tif 000002
【特許請求の範囲】
【請求項１】
全体座標系上における物体の、剛体としての変位である剛体変位と、弾性体としての振動である弾性振動とをコンピュータを用いて解析する方法であって、
前記全体座標系上を運動しつつ前記物体を局所的に観測するための局所観測座標系であって、その物体に固定されることなく物体とは独立しているが少なくともその物体にの近傍に位置して運動を行うものをその全体座標系と共に用いることにより、その物体の剛体変位と弾性振動とを、剛体変位は剛体モードにより、弾性振動は弾性モードによりそれぞれ表すとともに、それら剛体モードと弾性モードとを区別しないで一律に取り扱うモード解析により解析することを特徴とする物体変位・振動解析方法。
【請求項２】
前記局所観察座標系を用いて物体の剛体変位を表すことにより、局所観察座標上での変位が微小であるという条件を満足させる請求項１に記載の物体変位・振動解析方法。
【請求項３】
前記局所観察座標系を所定の周期で更新する請求項１または２に記載の物体変位・振動解析方法。
【請求項４】
前記局所観測座標系を時間の関数に従って運動させる請求項１ないし３に記載の物体変位・振動解析方法。
【請求項５】
前記局所観測座標系において、前記物体を局所的に観測することを周期的に繰り返すとともに、各観測周期のスタート時点においては前記局所観測座標系に物体の剛体変位と実質的に同じ運動を行わせ、スタート時点を除く期間においては物体の剛体変位とは独立に、時間の関数に従う運動を行わせ、前記物体の剛体変位と弾性振動との解析を、各観測周期ごとに行う請求項１に記載の物体変位・振動解析方法。
【請求項６】
前記局所観測座標系に、各観測周期のうちスタート時点を除く期間において、その局所観測座標系の変位と姿勢と速度と加速度と角速度と角加速度との中の少なくとも一つが、前記スタート時点における前記物体のそれと実質的に等しい大きさから変化しない運動を行わせるとともに、前記局所観測座標系を、各観測周期のエンド時点において前記物体の解析結果が得られるごとに、次回の観測に備えて、その物体の剛体変位に関する最新の解析結果と実質的に一致するように更新する請求項１ないし５のいずれかに記載の物体変位・振動解析方法。
【請求項７】
前記モード解析を、
(i)そのモード解析のための固有値λと、前記剛体モードを含むｎ個の正規直交化された複数のモードベクトルφj とを取得し、
(ii)前記局所観測座標系により定義されたベクトルを前記全体座標系により定義されるベクトルに変換するための方向余弦マトリクスＡを計算し、
(iii)最初の観測周期のスタート時点の前記全体座標系上における前記局所観測座標系の原点の位置ベクトルｒ，速度ベクトルｖおよび加速度ベクトルａを設定し、
(iv)各観測周期ごとに、前記物体に設定された複数の節点ｉの各々に作用する外力Ｆi を取得し、
(v)各観測周期ごとに、その観測周期のエンド時点における前記物体のモード座標系上におけるモード次数がｊであるモード変位ｑj ，モード速度ｑ.jおよびモード加速度ｑ:j を計算し、
(vi)前記各観測周期ごとに、次回の観測に備えて、その次回の観測周期のスタート時点の前記全体座標系上における前記局所観測座標系の原点の位置ベクトルｒdt(0) ，速度ベクトルｖdt(0) および加速度ベクトルａdt(0) を計算し、
(vii)前記各観測周期ごとに、次回の観測に備えて、前記物体の、モード次数がｊである剛体モードのモード変位ｑj ，モード速度ｑ.j およびモード加速度ｑ:j をいずれも０に初期化する
ことによって行う請求項１ないし６のいずれかに記載の物体変位・振動解析方法。
【請求項８】
当該物体変位・振動解析方法は、力の相互作用を有する複数の物体の剛体変位と弾性振動とを解析するものであり、
前記(iv)において、
各物体の節点ｉに作用する外力Ｆi を、関連する他の物体の節点ｉ’との間における力の相互作用に基づいて計算するとともに、その力の相互作用を、前記全体座標系上における前記複数の物体の、それらの弾性振動を考慮しない場合の各位置ベクトルｒi ，ｒi'と、各速度ベクトルｒ.i ，ｒ.i'と、各加速度ベクトルｒ:i ，ｒ:i'との少なくとも一つに基づいて計算する請求項７に記載の物体変位・振動解析方法。
【請求項９】
前記(ii)において、前記物体に予め設定された複数の節点のうち、前記局所座標系の原点に対応する節点を基点として、互いに線形独立な３つのベクトルｓ’1 ，ｓ’2 ，ｓ’3 を生成する３つの節点を選択し、その生成された３つのベクトルｓ’1 ，ｓ’2 ，ｓ’3 に基づいて前記方向余弦マトリクスＡを計算する請求項７または８に記載の物体変位・振動解析方法。
【請求項１０】
前記各観測周期ごとに、次回の観測に備えて、その次回の観測周期のスタート時点における前記方向余弦マトリクスＡdt(0) ，前記局所観測座標系の角速度ωの歪対称マトリクスω*dt (0) およびその局所観測座標系の角加速度αの歪対称マトリクスα*dt (0) を計算し、計算した方向余弦マトリクスＡdt(0) ，歪対称マトリクスω*dt (0) および歪対称マトリクスα*dt (0) を、前記ｎ個の正規直交化された複数のモードベクトルφj により定義されるモード座標系により前記モード解析を行うことに起因する誤差が、方向余弦マトリクスＡdt(0) ，歪対称マトリクスω*dt (0) および歪対称マトリクスα*dt (0) に生じないように補正する請求項７ないし９のいずれかに記載の物体変位・振動解析方法。
【請求項１１】
前記物体がエンジンであり、エンジンの振動を解析する請求項１ないし１０のいずれかに記載の物体変位・振動解析方法。
【請求項１２】
前記物体がエンジンであり、エンジンの車体に対する振動を解析する請求項１ないし１２のいずれかに記載の物体変位・振動解析方法。
【請求項１３】
コンピュータにより、運動する物体の、剛体としての変位である剛体変位と、弾性体としての振動である弾性振動とを解析するためのプログラムを記録した媒体であって、
前記プログラムは、
全体座標系を設定し、
この全体座標系上を運動しつつ、前記物体を局所的に観測する局所観測座標系であって、その物体に固定されることなく物体とは独立しているがその物体の少なくとも近傍に位置して運動を行うものを設定し、
２つの座標系を利用して、物体の剛体変位を剛体モードにより表し、物体の弾性振動を弾性モードにより表し、
これら剛体モードと弾性モードを一律に取り扱うモード解析により物体の剛体変位および弾性振動を解析することを特徴とする記録媒体。
【請求項１４】
運動する物体の、剛体としての変位である剛体変位と、弾性体としての振動である弾性振動とを解析するための装置であって、
全体座標系を設定する手段と、
この全体座標系上を運動しつつ、前記物体を局所的に観測する局所観測座標系であって、その物体に固定されることなく物体とは独立しているがその物体の少なくとも近傍に位置して運動を行うものを設定する手段と、
２つの座標系を利用して、物体の剛体変位を剛体モードで表し、物体の弾性振動を弾性モードにより表す手段と、
これら剛体モードと弾性モードを一律に取り扱うモード解析により物体の剛体変位および弾性振動を解析する手段と
を含み、それら全体座標系を設定する手段，局所観測座標系を設定する手段，物体の剛体変位および弾性振動を表す手段，物体の剛体変位および弾性振動を解析する手段がコンピュータにより構成されたことを特徴とする物体変位・振動解析装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、全体座標系上における物体の、剛体としての変位である剛体変位と弾性体としての振動である弾性振動とを解析する技術に関するものである。
【背景技術】
【０００２】
自動車のエンジン等、可動部品を有する機械においては、いかにして振動・騒音を低減させるかが重要な課題である。そのため、振動・騒音をより高い精度で解析できるより実用的な技術の確立が要望されている。機械要素単体の振動伝達特性に基づく振動解析については、汎用の構造解析ソフトウェア等を利用することにより、ほぼ実用レベルで実施できる段階に達していると言える。そして、近年では、複数の機械要素間の相互作用を考慮することにより、より高い精度で振動解析を行い得る技術の確立に重点が移行しつつある。
【０００３】
複数の機械要素間の相互作用を取り扱うためには、可動部品である物体の剛体変位と弾性振動との連成を解析することが重要な課題の一つである。そして、その課題を解決する一理論として機構解析による理論が既に知られている。この機構解析理論は、その当初は、複数の剛体が関節において互いに結合されたものを適用対象とするとともに、関節による結合を拘束条件式で表した微分代数方程式を用いる理論として構築された。その後、この機構解析理論は、物体の弾性振動を弾性モードの重ね合わせにより表現することにより、その適用対象を剛体から弾性体へと拡張された。この機構解析理論においては、物体の位置はデカルト座標、方向はオイラーパラメータ、弾性振動はモード座標によりそれぞれ表されている。モード座標は、弾性振動をモード解析する際に使用される、互いに直交する複数のモードベクトルにより定義されるものである。
【発明の開示】
【発明が解決しようとする課題，課題解決手段および発明の効果】
【０００４】
しかし、この機構解析理論に基づく従来の物体変位・振動解析方法では、物体の剛体変位と弾性振動とを記述する運動方程式において種類の異なる複数の座標系が混在する。このことは、剛性運動と弾性振動との連成を理由に、上記運動方程式において連成項の数を増加させたり、各連成項の構成を複雑にする可能性が高いことを意味する。そのため、この従来の物体変位・振動解析方法では、物体の剛体変位と弾性振動とを解析するための運動方程式の取り扱いが複雑で、計算に多くの工程や時間が必要であった。
【０００５】
さらに、この従来の物体変位・振動解析方法は、多関節の機構運動において低次の弾性振動の効果を考慮する場合には適しているが、弾性自由度が大きい弾性振動において剛体変位の効果を考慮する場合には十分に適しているとは言えなかった。
【０００６】
このような事情を背景とし、本発明は、物体の剛体変位と弾性振動とを解析する技術であって、弾性自由度が大きい弾性振動において剛体変位の効果を考慮する場合にも適しているものを提供することを課題としてなされたものであり、本発明によって下記各態様が得られる。各態様は、請求項と同様に、項に区分し、各項に番号を付し、必要に応じて他の項の番号を引用する形式で記載する。これは、本明細書に記載の技術的特徴およびそれらの組合せのいくつかの理解を容易にするためであり、本明細書に記載の技術的特徴やそれらの組合せが以下の態様に限定されると解釈されるべきではない。
【０００７】
（１）全体座標系上における物体の、剛体としての変位である剛体変位と、弾性体としての振動である弾性振動とをコンピュータを用いて解析する方法であって、
前記全体座標系上を運動しつつ前記物体を局所的に観測するための局所観測座標系であって、その物体に固定されることなく物体とは独立しているがその物体の少なくとも近傍に位置して運動を行うものをその全体座標系と共に用いることにより、その物体の剛体変位と弾性振動とを、剛体変位は剛体モードにより、弾性振動は弾性モードによりそれぞれ表すとともに、それら剛体モードと弾性モードとを区別しないで一律に取り扱うモード解析により解析することを特徴とする物体変位・振動解析方法（請求項１）。
この方法においては、全体座標系とは別に物体を観測する座標系として、その物体の少なくとも近傍に位置して運動を行う局所観測座標系が用いられている。この局所観測座標系は、物体に固定されてその物体の剛体変位と終始同じ運動をさせられる従来の物体固定座標系とは異なる。しかし、局所観測座標系は、物体の少なくとも近傍に位置してその物体を局所的に観測する。
したがって、この方法によれば、物体の剛体変位をモード解析する際の前提、すなわち、剛体変位がモード座標系上において微小であるという前提が容易に成立する。
よって、この方法によれば、小さな剛体変位を伴う物体のみならず大きな剛体変位を伴う物体についても容易に振動解析を行い得る。
ところで、例えば、自動車のエンジンの振動騒音を解析するためには、エンジンの回転中にそれのシリンダブロックに対してクランク軸が剛体変位すること，それらシリンダブロックの軸受とクランク軸のピンとの間に油膜が形成され、かつ、その油膜が非線形のばねとして作用すること等を考慮することが大切である。そして、この例においては、クランク軸がシリンダブロックに対して剛体変位する大きさは、そのクランク軸の弾性振動に比較すればかなり大きなものとなる。そして、本項に記載の方法は、この例の振動解析にも適用可能である。
また、本項に記載の方法は、後述のジャイロスコープの振動解析のように、単独の物体の振動解析に適用することが可能であり、また、上記の例、すなわち、シリンダブロックは車体に対して変位し、クランク軸はそのシリンダブロックに対して変位するエンジンの振動解析のように、相互作用を有する複数の物体の振動解析に適用することも可能である。
また、本項に記載の方法において「物体の少なくとも近傍に位置する」という要件は、それが成立すれば、物体の剛体変位をモード解析する際の前提、すなわち、剛体変位がモード座標系上において微小であるという前提が成立することになる要件であり、例えば、物体と局所観測座標系との相対変位（例えば、後述の個別変位ベクトルｕ’i ）が微小範囲にあるという要件であると定義したり、物体と局所観測座標系とが実質的に一致するという要件であると定義することが可能である。
また、本項に記載の方法において「物体の剛体変位」は、その物体の弾性振動の影響を受ける場合と、受けない場合とがある。受ける場合には、弾性振動の解析結果を考慮して剛体変位を決定することが必要であるが、受けない場合には、弾性振動の解析結果を待つことなくそれとは独立に剛体変位を決定することが可能である。よって、物体の剛体変位がその物体の弾性振動の影響を受ける場合には、物体の弾性振動を考慮して局所観測座標系の運動を定義することが必要であるが、受けない場合には、物体の弾性振動とは独立に局所観測座標系の運動を定義することが可能である。
（２）前記局所観察座標系を用いて物体の剛体変位を表すことにより、局所観察座標上での変位が微小であるという条件を満足させる (1)項に記載の物体変位・振動解析方法（請求項２）。
（３）前記局所観察座標系を所定の周期で更新する (1)項または (2)項に記載の物体変位・振動解析方法（請求項３）。
（４）前記局所観測座標系を時間の関数に従って運動させる (1)項ないし (3)項のいずれかに記載の物体変位・振動解析方法（請求項４）。
この方法によれば、局所観測座標系の運動とは独立にモード解析を行うことが可能となり、モード解析のための計算を容易に簡単化し得る。
（５）前記局所観測座標系において、前記物体を局所的に観測することを周期的に繰り返すとともに、各観測周期のスタート時点においては前記局所観測座標系に物体の剛体変位と実質的に同じ運動を行わせ、スタート時点を除く期間においては物体の剛体変位とは独立に、時間の関数に従う運動を行わせ、前記物体の剛体変位と弾性振動との解析を、各観測周期ごとに行う (1)項に記載の物体変位・振動解析方法（請求項５）。
この方法においては、局所観測座標系が、各観測周期のスタート時点においては物体の剛体変位と実質的に同じ運動を行い、スタート時点を除く期間においては物体の剛体変位とは独立に、時間の関数に従って運動を行わせられる。したがって、局所観測座標系は、各観測周期のスタート時点においては物体と実質的に同じ運動を行わせられるため、各観測周期の間、物体の近傍に位置してその物体を局所的に観測することが可能となる。
したがって、この方法によれば、物体の剛体変位をモード解析する際の前提、すなわち、剛体変位がモード座標系上において微小であるという前提が容易に成立する。
よって、この方法によれば、小さな剛体変位を伴う物体のみならず大きな剛体変位を伴う物体についても容易に振動解析を行い得る。
本項に記載の方法において「時間の関数」は、各観測周期のうちそれのスタート時点を除く期間においても、局所観測座標系が物体の近傍に位置し、前述の、モード解析における微小変位の前提が成立するように設定することが望ましい。
また、この方法において「周期的に」という用語は、各観測周期の長さが常に変化しないことを意味する場合のみならず、変化することを意味する場合をも包含する趣旨で用いられている。
（６）前記モード解析を、モード座標を一般座標として用いることによって前記物体の剛体変位と弾性振動とを記述する運動方程式を解くことにより行う (1)項ないし (5)項のいずれかに記載の物体変位・振動解析方法。
この方法においては、モード座標を一般座標として用いることによって物体を記述する運動方程式を解くことにより、物体の剛体変位と弾性振動とが解析される。一方、局所観測座標系は、各観測周期のスタート時点を除く期間においては物体の剛体変位とは独立に、時間の関数に従って運動を行わせられる。したがって、物体を記述する運動方程式において、局所観測座標系の運動に関連する変数をモード座標に対して定数として扱うことが可能となり、その運動方程式を用いて物体の剛体変位と弾性振動とを計算することが容易になる。
ここに「運動方程式」の一例はラグランジュの運動方程式である。また「モード座標」は一般に、互いに交差する複数のモードベクトルにより定義される座標であると定義される。
（７）前記モード解析を、互いに直交する複数のモードベクトルにより定義されるモード座標を用いて行う (1)項ないし (6)項のいずれかに記載の物体変位・振動解析方法。
この方法においては、互いに直交する複数のモードベクトルにより定義されるモード座標を用いてモード解析が行われる。したがって、この方法によれば、物体を記述する運動方程式において剛体モードと弾性モードとの間に生じる連成項の数を少なくしたり、連成項の構成を簡単にすることが容易になり、ひいては、運動方程式を用いた計算が容易になる。
ここに「モードベクトル」は、固有モードベクトルを選ぶことができるのはもちろんであるが、その他のモードベクトルを選ぶことも可能である。例えば、静的補正モードベクトルを選ぶことも可能であり、その静的補正モードベクトルには例えば、アタッチメントモードベクトル，拘束モードベクトル等がある。
（８）前記局所観測座標系に、各観測周期のうちスタート時点を除く期間において、その局所観測座標系の変位と姿勢と速度と加速度と角速度と角加速度との中の少なくとも一つが、前記スタート時点における前記物体のそれと実質的に等しい大きさから変化しない運動を行わせるとともに、前記局所観測座標系を、各観測周期のエンド時点において前記物体の解析結果が得られるごとに、次回の観測に備えて、その物体の剛体変位に関する最新の解析結果と実質的に一致するように更新する (1)項ないし (7)項のいずれかに記載の物体変位・振動解析方法（請求項６）。
この方法によれば、局所観測座標系の、各観測周期のスタート時点を除く期間における運動を比較的簡単に定義し得る。
（９）前記モード解析を、
(i)そのモード解析のための固有値λと、前記剛体モードを含むｎ個の正規直交化された複数のモードベクトルφj とを取得し、
(ii)前記局所観測座標系により定義されたベクトルを前記全体座標系により定義されるベクトルに変換するための方向余弦マトリクスＡを計算し、
(iii)最初の観測周期のスタート時点の前記全体座標系上における前記局所観測座標系の原点の位置ベクトルｒ，速度ベクトルｖおよび加速度ベクトルａを設定し、
(iv)各観測周期ごとに、前記物体に設定された複数の節点ｉの各々に作用する外力Ｆi を取得し、
(v)各観測周期ごとに、その観測周期のエンド時点における前記物体のモード座標系上におけるモード次数がｊであるモード変位ｑj ，モード速度ｑ.jおよびモード加速度ｑ:j を計算し、
(vi)前記各観測周期ごとに、次回の観測に備えて、その次回の観測周期のスタート時点の前記全体座標系上における前記局所観測座標系の原点の位置ベクトルｒdt(0) ，速度ベクトルｖdt(0) および加速度ベクトルａdt(0) を計算し、
(vii)前記各観測周期ごとに、次回の観測に備えて、前記物体の、モード次数がｊである剛体モードのモード変位ｑj ，モード速度ｑ. j およびモード加速度ｑ: j をいずれも０に初期化する
ことによって行う (1)項ないし (8)項のいずれかに記載の物体変位・振動解析方法（請求項７）。
（１０）当該物体変位・振動解析方法は、力の相互作用を有する複数の物体の剛体変位と弾性振動とを解析するものであり、
前記(iv)において、
各物体の節点ｉに作用する外力Ｆi を、関連する他の物体の節点ｉ’との間における力の相互作用に基づいて計算するとともに、その力の相互作用を、前記全体座標系上における前記複数の物体の、それらの弾性振動を考慮しない場合の各位置ベクトルｒi ，ｒi'と、各速度ベクトルｒ. i ，ｒ. i'と、各加速度ベクトルｒ: i ，ｒ: i'との少なくとも一つに基づいて計算する (9)項に記載の物体変位・振動解析方法（請求項８）。
この方法によれば、力の相互作用を有する複数の物体の剛体変位と弾性振動とを解析することができる。
（１１）前記(ii)において、前記物体に予め設定された複数の節点のうち、前記局所座標系の原点に対応する節点を基点として、互いに線形独立な３つのベクトルｓ’1 ，ｓ’2 ，ｓ’3 を生成する３つの節点を選択し、その生成された３つのベクトルｓ’1 ，ｓ’2 ，ｓ’3 に基づいて前記方向余弦マトリクスＡを計算する (9)項または(10)項に記載の物体変位・振動解析方法（請求項９）。
この方法によれば、方向余弦マトリクスＡを簡単に計算することができる。
（１２）前記各観測周期ごとに、次回の観測に備えて、その次回の観測周期のスタート時点における前記方向余弦マトリクスＡdt(0) ，前記局所観測座標系の角速度ωの歪対称マトリクスω*dt (0) およびその局所観測座標系の角加速度αの歪対称マトリクスα*dt (0) を計算し、計算した方向余弦マトリクスＡdt(0) ，歪対称マトリクスω*dt (0) および歪対称マトリクスα*dt (0) を、前記ｎ個の正規直交化された複数のモードベクトルφj により定義されるモード座標系により前記モード解析を行うことに起因する誤差が、方向余弦マトリクスＡdt(0) ，歪対称マトリクスω*dt (0) および歪対称マトリクスα*dt (0) に生じないように補正する (9)項ないし(11)項のいずれかに記載の物体変位・振動解析方法（請求項１０）。
（１３）前記物体がエンジンであり、エンジンの振動を解析する (1)項ないし(12)項のいずれかに記載の物体変位・振動解析方法（請求項１１）。
（１４）前記物体がエンジンであり、エンジンの車体に対する振動を解析する (1)項ないし(13)項のいずれかに記載の物体変位・振動解析方法（請求項１２）。
（１５） (1)項ないし(14)項のいずれかに記載の物体変位・振動解析方法を実施するためにコンピュータにより実行されるプログラムをコンピュータ読み取り可能に記録した記録媒体。
ここにおける「記録媒体」には例えば、フレキシブルディスク，磁気テープ，磁気ディスク，磁気ドラム，磁気カード，光ディスク，光磁気ディスク，ＲＯＭ，ＣＤ−ＲＯＭ，ＩＣカード，穿孔テープ等がある。
（１６）コンピュータにより、運動する物体の、剛体としての変位である剛体変位と、弾性体としての振動である弾性振動とを解析するためのプログラムを記録した媒体であって、
前記プログラムは、
全体座標系を設定し、
この全体座標系上を運動しつつ、前記物体を局所的に観測する局所観測座標系であって、その物体に固定されることなく物体とは独立しているがその物体の少なくとも近傍に位置して運動を行うものを設定し、
２つの座標系を利用して、物体の剛体変位を剛体モードにより表し、物体の弾性振動を弾性モードにより表し、
これら剛体モードと弾性モードを一律に取り扱うモード解析により物体の剛体変位および弾性振動を解析することを特徴とする記録媒体（請求項１３）。
（１７）運動する物体の、剛体としての変位である剛体変位と、弾性体としての振動である弾性振動とを解析するための装置であって、
全体座標系を設定する手段と、
この全体座標系上を運動しつつ、前記物体を局所的に観測する局所観測座標系であって、その物体に固定されることなく物体とは独立しているがその物体の少なくとも近傍に位置して運動を行うものを設定する手段と、
２つの座標系を利用して、物体の剛体変位を剛体モードで表し、物体の弾性振動を弾性モードにより表す手段と、
これら剛体モードと弾性モードを一律に取り扱うモード解析により物体の剛体変位および弾性振動を解析する手段と
を含み、それら全体座標系を設定する手段，局所観測座標系を設定する手段，物体の剛体変位および弾性振動を表す手段，物体の剛体変位および弾性振動を解析する手段がコンピュータにより構成されたを含むことを特徴とする物体変位・振動解析装置（請求項１４）。
【発明を実施するための最良の形態】
【０００８】
以下、本発明のさらに具体的な一実施形態を図面に基づいて詳細に説明する。
【０００９】
本実施形態は、物体の剛体変位と弾性振動とを解析する物体変位・振動解析方法である。この方法は、剛体自由度は小さいが弾性自由度は大きい機械の振動騒音の解析に適している。
【００１０】
本実施形態は、種々の新たな理論に基づいて構成されており、以下、その新たな理論を説明する。
【００１１】
１．局所観測座標系の導入
【００１２】
図１に示すように、全体座標系｛Ｘ，Ｙ，Ｚ｝に対し、局所観測座標系｛ｘ，ｙ，ｚ｝なる新しい概念の座標系を導入する。この座標系は、図２に示すように、従来の、物体に固定された座標系とは異なり、物体の近傍に位置しながら、物体の運動を局所的に観測するために設けられた座標系である。局所観測座標系ＬＯＣは、図３に示すように、ある時刻に物体Ｏと同一の運動をし、その状態から、物体Ｏの運動とは独立に、時間の関数に従って運動をする。
【００１３】
ここで、局所観測座標系と物体との関係を、位置，速度および加速度に関してそれぞれ説明する。
【００１４】
まず、図１に示すように、全体座標系から見た局所観測座標系の原点の位置ベクトルをｒ、方向余弦マトリクスをＡ、局所観測座標系から見た物体の節点ｉの位置ベクトルをｒ’i とすると、全体座標系から見た物体の節点ｉの位置ベクトルｒi は、
【００１５】
【数１】
2005091370.tif 000003
【００１６】
なる式で表される。
【００１７】
ここで補足説明すれば、上記方向余弦マトリクスＡは、座標系を局所観測座標系から全体座標系に変換する座標変換マトリクスの一種である。また、上記位置ベクトルｒ’i において「’」は、そのベクトルが局所観測座標系により定義されることを示す表記であり、この表記は他の変数を表記する際にも用いられる。また、観測すべき物体には、その物体を複数の要素に仮想的に分割するためにＮ個の節点が設定されており、上記節点ｉはそれらＮ個の節点のうち任意のものを示している。
【００１８】
物体の剛体変位および弾性振動を考慮しない場合の、局所観測座標系から見た節点ｉの位置ベクトルをｓ’i とし、それら剛体変位と弾性振動との和を個別変位ベクトルｕ’i とすると、
【００１９】
【数２】
2005091370.tif 000004
【００２０】
なる式が成立する。ここで、Ｎ個の節点の全変位ベクトルｕ’を、
【００２１】
【数３】
2005091370.tif 000005
【００２２】
なる式で定義する。ここで、全変位ベクトルｕ’が十分小さく、モード解析が前提とする微小変位の範囲内にあり、かつ、全変位ベクトルｕ’が、弾性モードのみならず剛体モードをも含むｎ個の正規直交化された固有モードベクトルφj （モード次数ｊ＝１，・・・，ｎ）の重ね合わせにより表すことができるとすると、全変位ベクトルｕ’が、モード変位ｑj を用いることにより、
【００２３】
【数４】
2005091370.tif 000006
【００２４】
なる式で表される。ここで、モードマトリクスΦを、固有モードベクトルφj を用いることにより、
【００２５】
【数５】
2005091370.tif 000007
【００２６】
なる式で定義すると、モードマトリクスΦは、固有モードの直交性を理由に、質量マトリクスＭと剛性マトリクスＫとの間において、
【００２７】
【数６】
2005091370.tif 000008
【数７】
2005091370.tif 000009
【００２８】
なる式で表される関係を有する。ここで、「Ｉ」は単位マトリクスであり、「λ」は、固有値λj を対角成分として有するマトリクス、すなわち、
【００２９】
【数８】
2005091370.tif 000010
【００３０】
である。ただし、モード次数ｊ＝１，・・・，ｎr が剛体モードを表すものであるとすると、固有値λi ＝０（ｊ＝１，・・・，ｎr ）となる。ここで、モード変位ベクトルｑを、マトリクスの転置を表す表記として「T 」を用いることにより、
【００３１】
【数９】
2005091370.tif 000011
【００３２】
なる式で定義すると、前記式(4) は、前記式(5) を用いることにより、
【００３３】
【数１０】
2005091370.tif 000012
【００３４】
に変形される。前記式(3) における全変位ベクトルｕ’と個別変位ベクトルｕ’i との関係と同様に、固有モードベクトルφj を、各節点ｉごとに分けてφijで表すことにより、
【００３５】
【数１１】
2005091370.tif 000013
【００３６】
なる式で定義すると、前記式(2) における個別変位ベクトルｕ’i は、
【００３７】
【数１２】
2005091370.tif 000014
【００３８】
なる式で表される。したがって、前記式(1) は、前記式(2) と上記式(12)とを用いることにより、
【００３９】
【数１３】
2005091370.tif 000015
【００４０】
に変形される。この式(13)における「ｒi 」は、図１に示すものと同じである。また、この式(13)は、局所観測座標系と物体との間の、位置に関する関係式を意味する。
【００４１】
上記式(13)を時間で微分すると、局所観測座標系と物体との間の、速度に関する関係式として、
【００４２】
【数１４】
2005091370.tif 000016
【００４３】
なる式が得られる。ここに、「Ａ. 」は、方向余弦マトリクスＡを時間で微分したものを意味し、「ｑ. j 」は、モード変位ｑj を時間で微分したものを意味する。
【００４４】
ここで、全体座標系に対する局所観測座標系の角速度ベクトルをω、その歪対称マトリクスをω* で表すと、
【００４５】
【数１５】
2005091370.tif 000017
【００４６】
なる式で表される関係が成立する。この式(15)を上記式(14)に代入し、全体座標系に対する局所観測座標系の速度ベクトルを改めてｖで表すと、速度の関係式として、
【００４７】
【数１６】
2005091370.tif 000018
【００４８】
なる式が得られる。
【００４９】
前記式(14)を時間で微分すると、局所観測座標系と物体との間の、加速度に関する関係式として、
【００５０】
【数１７】
2005091370.tif 000019
【００５１】
なる式が得られる。ここに、「Ａ: 」は、前記Ａ. を時間で微分したものを意味し、「ｑ: j 」は、前記ｑ. j を時間で微分したものを意味する。一方、前記式(15)を時間で微分すると、
【００５２】
【数１８】
2005091370.tif 000020
【００５３】
なる式が得られる。ここに、「ω.*」は、歪対称マトリクスω* を時間で微分したものを意味する。この式(18)に前記式(15)を代入するとともに、全体座標系に対する局所観測座標系の角加速度ベクトルをα、その歪対称マトリクスをα* としてそれぞれ表すと、
【００５４】
【数１９】
2005091370.tif 000021
【００５５】
なる式が得られる。この式(19)と前記式(15)とを前記式(17)に代入するとともに、全体座標系に対する局所観測座標系の加速度ベクトルを改めてａで表すと、
【００５６】
【数２０】
2005091370.tif 000022
【００５７】
なる式が得られる。
【００５８】
以上のようにして、局所観測座標系｛ｘ，ｙ，ｚ｝を導入することにより、前記式(13)，式(16)および式(20)が、弾性振動のみならず剛体変位をも伴う物体の節点ｉの位置，速度および加速度を、全体座標系と局所観測座標系とモード座標系とを用いて記述する式として得られた。
【００５９】
なお付言すれば、それら３式において全体座標系または局所観測座標系により定義される変数は、後述のように、物体の剛体変位とは独立に、時間の関数として定義される。すなわち、その変数は、時間依存変数と言うことができるのである。一方、本実施形態において採用された物体の運動方程式は、後述のように、モード座標を一般座標として用いる。したがって、その運動方程式において上記時間依存変数は、一般座標との関係においては定数として取り扱うことが可能となる。このような観点からすれば、前記３式は、局所観測座標系の導入により、弾性振動のみならず剛体変位をも伴う物体の節点ｉの位置，速度および加速度を、モード座標系のみを用いて上記運動方程式に組み込むことを可能にした式であると言うことができ、また、それら位置，速度および加速度を記述する座標系を、全体座標系および局所観測座標系からモード座標系に変換することを可能にした式であると言うこともできる。
【００６０】
２．モード座標系による運動方程式の導出
【００６１】
物体全体の運動エネルギーＴは、節点ｉの質量をｍi で表すと、
【００６２】
【数２１】
2005091370.tif 000023
【００６３】
なる式で表される。また、弾性変形による歪みエネルギーＵは、前記式(3) の全変位ベクトルｕ’を用いることにより、
【００６４】
【数２２】
2005091370.tif 000024
【００６５】
なる式で表される。この式(22)に前記式(10)を代入すると、
【００６６】
【数２３】
2005091370.tif 000025
【００６７】
なる式が得られ、さらに、前記式(7) により表される関係を利用すると、
【００６８】
【数２４】
2005091370.tif 000026
【００６９】
なる式が得られる。
【００７０】
一方、ラグランジュの運動方程式を、モード座標ｑk （モード次数ｋ＝１，・・・，ｎr ）を一般座標（generalized coordinates, general coordenates）として用いて表すと、
【００７１】
【数２５】
2005091370.tif 000027
【００７２】
なる式が得られる。この式(25)は、時間に依存しないため、局所観測座標系を、物体の運動とは独立に、時間の関数に従って運動させることによる影響を何ら受けない。したがって、この式(25)は、局所観測座標系の運動を考慮することなく、解くことができるのである。
【００７３】
上記式(25)の右辺における「Ｑk 」は、モード座標系での一般化力（generalized force ）を表している。
【００７４】
左辺の第１項の括弧内の変数は、前記式(21)を用いることにより、
【００７５】
【数２６】
2005091370.tif 000028
【００７６】
で表される。この式(26)は、前記式(16)を用いることにより、
【００７７】
【数２７】
2005091370.tif 000029
【００７８】
なる式に変形され、転置して整理すると、
【００７９】
【数２８】
2005091370.tif 000030
【００８０】
なる式が得られる。この式(28)を時間で微分し、前記式(15)を用いて整理すると、
【００８１】
【数２９】
2005091370.tif 000031
【００８２】
なる式が得られる。
【００８３】
一方、前記式(25)の左辺の第２項は、前記式(21)を用いることにより、
【００８４】
【数３０】
2005091370.tif 000032
【００８５】
で表される。この式(30)は、前記式(16)を用いることにより、
【００８６】
【数３１】
2005091370.tif 000033
【００８７】
に変形され、転置して整理すると、
【００８８】
【数３２】
2005091370.tif 000034
【００８９】
なる式が得られる。この式(32)の右辺の項は、前記式(29)の右辺の第２項と同一であるから、それら２式を引き算すれば、
【００９０】
【数３３】
2005091370.tif 000035
【００９１】
なる式が得られる。この式(33)に前記式(20)を代入すると、
【００９２】
【数３４】
2005091370.tif 000036
【００９３】
なる式が得られる。ここに、方向余弦マトリクスＡは正規直交マトリクスであるため、
【００９４】
【数３５】
2005091370.tif 000037
【００９５】
なる式が成立する。また、前記式(6) により表される、固有モードの直交性という条件は、節点ｉの質量ｍi を用いて表すと、
【００９６】
【数３６】
2005091370.tif 000038
【００９７】
となる。ただし、「δij」は、クロネッカーのデルタである。この式(36)と前記式(35)を用いると、前記式(34)は、
【００９８】
【数３７】
2005091370.tif 000039
【００９９】
なる式に変形される。
【０１００】
一方、前記式(25)の左辺の第３項は、前記式(24)と前記式(8) とを用いることにより、
【０１０１】
【数３８】
2005091370.tif 000040
【０１０２】
なる式で表される。
【０１０３】
また、前記式(25)の右辺の一般化力Ｑk については、まず、全体座標系から見た節点ｉに作用する外力をＦi として、仮想仕事δＷを考えるととにも、
【０１０４】
【数３９】
2005091370.tif 000041
【０１０５】
なる式を用いる。ここに、「δｒi 」は、前記位置ベクトルｒi の仮想的な微小変化を表す。この式(39)は、前記式(13)を用いることにより、
【０１０６】
【数４０】
2005091370.tif 000042
【０１０７】
なる式に変形される。ここに、「δｑj 」は、前記モード変位ｑj の仮想的な微小変化を表す。ここで、微小変化δｑj の係数をまとめることにより、
【０１０８】
【数４１】
2005091370.tif 000043
【０１０９】
なる式を定義すると、一般化力Ｑk は、
【０１１０】
【数４２】
2005091370.tif 000044
【０１１１】
なる式で表され、転置して整理すると、
【０１１２】
【数４３】
2005091370.tif 000045
【０１１３】
なる式が得られる。
【０１１４】
前記式(37)，式(38)および式(43)を前記式(25)、すなわち、ラグランジュの運動方程式に代入して整理すると、最終的に、モード座標系により記述された運動方程式として、
【０１１５】
【数４４】
2005091370.tif 000046
【０１１６】
なる式が得られる。
【０１１７】
この式(44)には、「φijｑj 」なる項と「φijｑ. j 」なる項とが存在し、それらは、異なる固有モード間で振動が連成することを表す連成項であるが、前述の従来の物体変位・振動解析方法を用いて同じ振動騒音問題を解析する場合に比較して、連成項の数が少なく、また、その構成も単純である。
【０１１８】
３．局所観測座標系の更新
【０１１９】
本実施形態においては、前述のように、局所観測座標系が導入されている。しかし、この座標系を同じ位置で使用し続けると、時間の経過に伴い、前記式(2) における個別変位ベクトルｕ’i が微小であるという仮定が成立しなくなってしまう。そのため、本実施形態においては、個別変位ベクトルｕ’i が微小である状態を維持するために局所観測座標系が適切に更新されるようになっている。以下、その更新の方法を説明する。
【０１２０】
(1) 各種記号の定義
【０１２１】
図４には、時刻０における局所観測座標系｛ｘ0 ，ｙ0 ，ｚ0 ｝と、時刻０からｄｔ秒が経過したとき（以下、「時刻ｄｔ」という）における局所観測座標系｛ｘdt，ｙdt，ｚdt｝とが示されている。時刻０においては、局所観測座標系の運動と物体の運動とは完全に一致し、その時刻を基点に、局所観測座標系は物体の運動とは独立に、時間の関数に従って運動をする。
【０１２２】
本実施形態においては、局所観測座標系が、時刻０から時刻ｄｔまでの間、すなわち、一回の観測周期のうちそのスタート時点を除く期間中、その局所観測座標系の、変位と姿勢とから成る組と、変位と姿勢と速度と角速度とから成る組と、変位と姿勢と速度と角速度と加速度と角加速度とから成る組とのいずれかが、時刻０における物体のそれと実質的に等しい大きさから実質的に変化しない運動を行わせられる。そして、例えば、局所観測座標系に等加速度・等角加速度運動をさせる場合には、局所観測座標系の運動を定義する時間の関数は例えば、局所観測座標系の速度については後述の式(86)、位置については式(87)、角速度については式(88)がそれぞれ用いられる。
【０１２３】
上記において「実質的に等しい」とは、局所観測座標系が物体に完全に一致する場合のみならず、変位ｕ' i が微小となる範囲に局所観測座標系が更新される場合も含む意味である。後者の場合、更新後の局所観測座標系上での剛体モードの変位，速度，加速度はゼロとならないので、これらを計算する過程が必要となる。
【０１２４】
一方、時刻ｄｔにおける局所観測座標系の加速度ベクトルをａ0 (dt)、速度ベクトルをｖ0 (dt)、位置ベクトルをｒ0 (dt)、角加速度ベクトルをα0 (dt)、角速度ベクトルをω0 (dt)、方向余弦マトリクスをＡ0 (dt)でそれぞれ表す。また、時刻ｄｔにおいて、局所観測座標系を更新し、その新しい局所観測座標系の加速度ベクトルをａdt(0) 、速度ベクトルをｖdt(0) 、位置ベクトルをｒdt(0) 、角加速度ベクトルをαdt(0) 、角速度ベクトルをωdt(0) 、方向余弦マトリクスをＡdt(0) でそれぞれ表すとともに、それらを求める。
【０１２５】
なお、以上説明した各種変数または関数の右肩添字は、その時刻に決定された変数または関数であることを示している。
【０１２６】
(2) 局所観測座標系の原点の位置ベクトルｒdt(0) に関する更新
【０１２７】
時刻ｄｔにおける新しい局所観測座標系の原点の位置ベクトルｒdt(0) は、モード次数ｊ＝１，・・・，ｎr が剛体モードを表すものであるとすると、前記式(13)を用い、かつ、原点（ｉ＝０）であることを考慮することにより、
【０１２８】
【数４５】
2005091370.tif 000047
【０１２９】
なる式で表されるように更新される。ここに、「φ0j」は、物体の複数の節点のうち、局所観測座標系の原点に対応する節点の固有モードベクトルを示している。なお、「原点に対応する節点」の意味については後述する。
【０１３０】
(3) 局所観測座標系の原点の速度ベクトルｖdt(0) に関する更新
【０１３１】
時刻ｄｔにおける新しい局所観測座標系の原点の速度ベクトルｖdt(0) は、前記式(16)を用いることにより、
【０１３２】
【数４６】
2005091370.tif 000048
【０１３３】
なる式で表されるように更新される。
【０１３４】
(4) 局所観測座標系の原点の加速度ベクトルａdt(0) に関する更新
【０１３５】
時刻ｄｔにおける新しい局所観測座標系の原点の加速度ベクトルａdt(0) は、前記式(20)を用いることにより、
【０１３６】
【数４７】
2005091370.tif 000049
【０１３７】
なる式で表されるように更新される。
【０１３８】
(5) 方向余弦マトリクスＡdt(0) の更新
【０１３９】
物体の剛体変位を考慮せず、かつ、時刻ｄｔである場合の、全体座標系から見た節点ｉの位置ベクトルｒdti は、前記式(13)を用いることにより、
【０１４０】
【数４８】
2005091370.tif 000050
【０１４１】
なる式で表される。また、時刻ｄｔに、前記位置ベクトルｓ’i を全体座標系から見た場合のベクトルｓdti は、
【０１４２】
【数４９】
2005091370.tif 000051
【０１４３】
なる式で表され、この式(49)に前記式(45)と式(48)とを代入すると、
【０１４４】
【数５０】
2005091370.tif 000052
【０１４５】
なる式が得られる。また、上記ベクトルｓdti は、時刻ｄｔにおける方向余弦マトリクスＡdt(0) を用いることにより、
【０１４６】
【数５１】
2005091370.tif 000053
【０１４７】
なる式で表される。この式(51)を方向余弦マトリクスＡdt(0) について解くために、図４に示すように、物体の複数の節点の中から、局所観測座標系の原点と共同して四面体を構成する３つの節点であって、かつ、原点を始点、各節点を終点とする３つのベクトルｓ’1 ，ｓ’2 およびｓ’3 が互いに線形独立であるものを選択する。この選択により、
【０１４８】
【数５２】
2005091370.tif 000054
【０１４９】
なる式が得られる。ここに、「ｓ’1 」，「ｓ’2 」および「ｓ’3 」はそれぞれ、同図において、原点を始点、節点１，２および３をそれぞれ終点とする、互いに線形独立な３つのベクトルを意味している。このように、ベクトルｓ’1 ，ｓ’2 およびｓ’3 は互いに線形独立であるため、それらベクトルｓ’1 ，ｓ’2 およびｓ’3 を成分とするマトリクスには逆マトリクスが存在する。よって、上記式(52)における方向余弦マトリクスＡdt(0) は、
【０１５０】
【数５３】
2005091370.tif 000055
【０１５１】
なる式で表されるように更新することができる。
【０１５２】
(6) 角速度の歪対称マトリクスω*dt (0) の更新
【０１５３】
時刻ｄｔにおける節点ｉの、物体の剛体変位に起因する速度成分ｒ.dt i は、前記式(16)を用いることにより、
【０１５４】
【数５４】
2005091370.tif 000056
【０１５５】
なる式で表される。ベクトルｓdti を時間で微分して得られるベクトルｓ.dt i は、全体座標系に対して並進方向の速度成分を持たないため、
【０１５６】
【数５５】
2005091370.tif 000057
【０１５７】
なる式で表される。この式(55)に前記式(46)と上記式(54)とを代入すると、
【０１５８】
【数５６】
2005091370.tif 000058
【０１５９】
なる式が得られる。さらに、ベクトルｓ.dt i は、前記式(53)と時刻ｄｔにおける歪対称マトリクスω*dt (0) とを用いることにより、
【０１６０】
【数５７】
2005091370.tif 000059
【０１６１】
なる式で表される。したがって、前記式(53)と同様に、前述の、互いに線形独立なベクトルｓ’1 ，ｓ’2 およびｓ’3 を用いることにより、
【０１６２】
【数５８】
2005091370.tif 000060
【０１６３】
なる式が求められる。よって、角速度の歪対称マトリクスω*dt (0) は、
【０１６４】
【数５９】
2005091370.tif 000061
【０１６５】
なる式で表されるように更新される。
【０１６６】
(7) 角加速度の歪対称マトリクスα*dt (0) の更新
【０１６７】
時刻ｄｔにおける節点ｉの、物体の剛体変位に起因する加速度成分ｒ:dt i は、前記式(20)を用いることにより、
【０１６８】
【数６０】
2005091370.tif 000062
【０１６９】
なる式で表される。ここに、「ｓ:dt i 」は、全体座標系に対して並進方向の加速度成分を持たないため、
【０１７０】
【数６１】
2005091370.tif 000063
【０１７１】
なる式で表される。この式(61)に前記式(47)と上記式(60)とを代入することにより、          
【０１７２】
【数６２】
2005091370.tif 000064
【０１７３】
なる式が得られる。さらに、ベクトルｓ:dt i は、前記式(53)と式(59)と時刻ｄｔにおける角加速度の歪対称マトリクスα*dt (0) とを用いることにより、
【０１７４】
【数６３】
2005091370.tif 000065
【０１７５】
なる式で表される。よって、前記式(53)および式(58)と同様に、前述の、互いに線形独立なベクトルｓ’1 ，ｓ’2 およびｓ’3 を用いることにより、
【０１７６】
【数６４】
2005091370.tif 000066
【０１７７】
なる式が求められる。よって、角加速度の歪対称マトリクスα*dt (0) は、
【０１７８】
【数６５】
2005091370.tif 000067
【０１７９】
なる式で表されるように更新される。
【０１８０】
(8) 方向余弦マトリクスＡdt (0)の補正
【０１８１】
本実施形態においては、物体の剛体変位を剛体モードの重ね合わせにより解析し、その際、複数の固有モードベクトルが互いに直交するモード座標系を用いて物体の剛体変位を記述するため、前記式(53)により求めた方向余弦マトリクスＡdt (0)には微小な誤差が存在する。以下、この誤差の存在を具体的に説明する。
【０１８２】
全体座標系から見た局所観測座標系の座標軸ベクトルをｘ< ，ｙ< ，ｚ< で表すこととすると、方向余弦マトリクスＡは、
【０１８３】
【数６６】
2005091370.tif 000068
【０１８４】
で表すことができる。各座標軸ベクトルは単位ベクトルであるため、
【０１８５】
【数６７】
2005091370.tif 000069
【０１８６】
なる式が成立する。さらに、３つの座標軸ベクトルは互いに直交するため、
【０１８７】
【数６８】
2005091370.tif 000070
【０１８８】
なる式が成立する。
【０１８９】
しかし、前記式(53)により更新された方向余弦マトリクスＡdt (0)の各列ベクトルは、小さい誤差を含んだｘ- ，ｙ- ，ｚ- であるため、上記式(67)により表される関係および上記式(68)により表される関係が厳密には成立しない。
【０１９０】
そこで、本実施形態においては、その誤差が解消されるように、方向余弦マトリクスＡdt (0)が適切に補正される。以下、その補正方法を説明する。
【０１９１】
まず、列ベクトルｘ- の方向は正しいと仮定として、列ベクトルｘ< を、
【０１９２】
【数６９】
2005091370.tif 000071
【０１９３】
なる式を用いて決定される。
【０１９４】
次に、求めるべき列ベクトルｙ< は、上記のようにして求められた列ベクトルｘ< と列ベクトルｙ- と同一平面上にあるとして、列ベクトルｚ< が、
【０１９５】
【数７０】
2005091370.tif 000072
【０１９６】
なる式を用いて決定される。
【０１９７】
最後に、列ベクトルｙ< が、
【０１９８】
【数７１】
2005091370.tif 000073
【０１９９】
なる式を用いて決定される。
【０２００】
(9) 角速度の歪対称マトリクスω*dt (0) の補正
【０２０１】
方向余弦マトリクスＡdt (0)の場合と同様に、更新された角速度ωの歪対称マトリクスω*dt (0) にも微小な誤差が含まれる。そこで、本実施形態においては、歪対称マトリクスω*dt (0) も適切に補正される。以下、その補正方法を説明する。
【０２０２】
前記式(59)で求めた歪対称マトリクスω*dt (0) は、一般にフルマトリクスとなるが、本来は、角速度ベクトルω、すなわち、
【０２０３】
【数７２】
2005091370.tif 000074
【０２０４】
の歪対称マトリクスであって、
【０２０５】
【数７３】
2005091370.tif 000075
【０２０６】
なる式で表される関係を満たすはずである。ここで、誤差を含んだ歪対称マトリクスω*dt (0) の各成分を、
【０２０７】
【数７４】
2005091370.tif 000076
【０２０８】
なる式で表す。歪対称マトリクスω*dt (0) の補正方法は種々考えられるが、本実施形態においては、次のように行われる。
【０２０９】
すなわち、まず、歪対称マトリクスω*dt (0) の対角成分が０でない場合にはすべて誤差であるとして、対角成分ω* ii＝０（ｉ＝１，２，３）とされる。次に、基本成分ωx ，ωy ，ωz が、
【０２１０】
【数７５】
2005091370.tif 000077
【数７６】
2005091370.tif 000078
【数７７】
2005091370.tif 000079
【０２１１】
のように補正される。
【０２１２】
(10)角加速度の歪対称マトリクスα*dt (0) の補正
【０２１３】
歪対称マトリクスα*dt (0) の補正方法は、歪対称マトリクスω*dt (0) の場合と基本的に同じである。以下、具体的に説明する。
【０２１４】
前記式(65)により更新された歪対称マトリクスα*dt (0) は、一般にフルマトリクスとなるが、本来は、角加速度ベクトル、すなわち、
【０２１５】
【数７８】
2005091370.tif 000080
【０２１６】
の歪対称マトリクスであるため、
【０２１７】
【数７９】
2005091370.tif 000081
【０２１８】
なる式で表される関係を満たすはずである。ここで、誤差を含んだ歪対称マトリクスα*dt (0) の各成分を、
【０２１９】
【数８０】
2005091370.tif 000082
【０２２０】
なる式で表す。そして、歪対称マトリクスα*dt (0) の補正においては、まず、歪対称マトリクスα*dt (0) の対角成分が０でない場合にはすべて誤差であるとして、対角成分ω* ii＝０（ｉ＝１，２，３）とされる。次に、基本成分αx ，αy ，αz が、
【０２２１】
【数８１】
2005091370.tif 000083
【数８２】
2005091370.tif 000084
【数８３】
2005091370.tif 000085
【０２２２】
のように補正される。
【０２２３】
(11)剛体モードの各種変数の初期化
【０２２４】
本実施形態においては、局所観測座標系｛ｘ，ｙ，ｚ｝が更新されるごとに、剛体モードにおけるモード変位ｑj （ｊ＝１，・・・，ｎr ），モード速度ｑ. j （ｊ＝１，・・・，ｎr ）およびモード加速度ｑ: j （ｊ＝１，・・・，ｎr ）がすべて０に初期化される。
【０２２５】
以上、本実施形態を構成するために採用された理論を詳細に説明したが、次に、その構成を詳細に説明する。
【０２２６】
図５には、本実施形態である物体変位・振動解析方法を実施するのに好適な物体変位・振動解析装置（以下、単に「解析装置」という）が示されている。
【０２２７】
解析装置は、入力装置１０とコンピュータ１２と出力装置１４と外部記憶装置１６とを備えている。入力装置１０はマウス，キーボード等を含むように構成される。コンピュータ１２はＣＰＵ等、プロセッサ２０と、ＲＯＭ，ＲＡＭ，ハードディスク等、メモリ２２と、それらプロセッサ２０とメモリ２２とを接続するバス２４とを含むように構成される。出力装置１４はディスプレイ，プリンタ，プロッタ等を含むように構成される。外部記憶装置１６は、ＣＤ−ＲＯＭ３０，書き込み可能なＦＤ３２（フレキシブルディスク）等、記録媒体が装填可能となっていて、装填状態においては、記録媒体に対するデータの読み取りおよび書き込みが必要に応じて行われる。
【０２２８】
本実施形態においては、物体の変位・振動解析（以下、単に「振動解析」という）に必要なプログラムがＣＤ−ＲＯＭ３０に記憶され、また、その振動解析に必要なデータがＦＤ３２に記憶されている。振動解析に必要なプログラムには、モード解析プログラムと振動解析プログラムとがある。そして、振動解析時には、それらＣＤ−ＲＯＭ３０およびＦＤ３２から必要なプログラムおよびデータが読み出されてコンピュータ１２のＲＡＭまたはハードディスクに転送され、その後、プロセッサ２０によりそのプログラムが実行される。
【０２２９】
上記モード解析プログラムは、有限要素法により、解析すべき物体の固有値λj と固有モードベクトルφj を計算する。このプログラムは、物体を仮想的に複数の要素に分割するためのデータと、物体の属性を規定するためのデータとに基づき、固有値λj と固有モードベクトルφj を計算する。そして、このプログラムを利用すれば、各節点の番号ｉ（ｉ＝１，・・・，Ｎ）と、各節点の座標ｓ’i と、固有値λj （ｊ＝１，・・・，ｎr ）と、固有モードベクトルφj （ｊ＝１，・・・，ｎr ）とが取得できる。このプログラムの一例は、日本ＭＳＣにより開発されて「ＭＳＣ−ＮＡＳＴＲＡＮ（バージョン６８，６９，７０）」という商品名で販売されたプログラムである。
【０２３０】
図６および図７には、上記振動解析プログラムがフローチャートで表されている。このプログラムは、本発明者らにより提案された前述の振動解析理論に基づいて設計されており、よって、各ステップの内容は既に詳しく説明されている。そのため、以下、そのフローチャートを参照しつつ本実施形態である物体変位・振動解析方法の時間的な流れを詳細に説明し、各ステップの内容については簡単に説明する。
【０２３１】
本プログラムにおいては、まず、ステップＳ１（以下、単に「Ｓ１」で表す。他のステップについても同じ）において、物体を解析するために設定された時間軸に対して、その解析を開始する開始時刻ｓＴと、終了する終了時刻ｅＴとが入力されるとともに、それら開始時刻ｓＴと終了時刻ｅＴとの間の期間を細かく分割する刻み時間ｄｔが入力される。これにより、開始時刻ｓＴと終了時刻ｅＴとの間において物体の観測および解析を行う観測周期が、刻み時間ｄｔごとに繰り返されることになる。
【０２３２】
次に、Ｓ２において、前記モード解析プログラムの実行により、節点番号ｉ（ｉ＝１，・・・，Ｎ）と、節点座標ｓ’i と、固有値λj （ｊ＝１，・・・，ｎr ）と、固有モードベクトルφj とが取得される。その後、Ｓ３において、物体の各節点が概念的に有するとされる質量ｍi が設定される。
【０２３３】
続いて、Ｓ４において、物体の複数の節点のうち、開始時刻ｓＴにおける局所観測座標系の原点に対応するものが対応節点として選択される。一致する節点が存在する場合には直ちにそれが対応節点として選択されるが、存在しない場合には、複数の節点のいずれかが原点に一致するように局所観測座標系上で物体が平行移動させられ、その平行移動後に原点に一致する節点が対応節点とされる。その後、Ｓ５において、物体の複数の節点のうち、局所観測座標系の原点と共同することにより、互いに線形独立な３つのベクトルｓ’1 ，ｓ’2 ，ｓ’3 を生成する３つの節点を選択することにより、それら３つのベクトルｓ’1 ，ｓ’2 ，ｓ’3 が生成される。
【０２３４】
続いて、Ｓ６において、局所観測座標系の加速度ベクトルａ，速度ベクトルｖ，位置ベクトルｒ，角加速度ベクトルα，角速度ベクトルωおよび方向余弦マトリクスＡのそれぞれの初期値が設定される。その後、Ｓ７において、物体の、モード次数がｊである弾性モードのモード変位ｑj ，モード速度ｑ. j およびモード加速度ｑ: j のそれぞれの初期値が設定される。
【０２３５】
続いて、Ｓ８において、物体の各節点ｉに作用する外力Ｆi が取得される。単独の物体の振動解析のためにその物体に対して本プログラムが実行される場合には、ユーザからの入力によって外力Ｆi が取得される。これに対して、力の相互作用を有する複数の物体の振動解析のためにそれら複数の物体の各々に対して本プログラムが実行される場合には、本プログラムの最新の解析結果を利用することにより外力Ｆi が逐次計算される。具体的には、ある物体の節点ｉに作用する外力Ｆi は、関連する他の物体の節点ｉ’との間における相対変位と相対変位速度との少なくとも一方に基づいて計算される。相対変位は、各物体の位置ベクトルｒi と、関連する他の物体の位置ベクトルｒi'との差により計算され、また、相対変位速度は、各物体の速度ベクトルｒ. i と、関連する他の物体の速度ベクトルｒ. i'との差により計算される。また、各位置ベクトルｒi ，ｒi'は前記式(13)により、各速度ベクトルｒ. i ，ｒ. i'は前記式(16)により、各加速度ベクトルｒ: i ，ｒ: i'は前記式(17)によりそれぞれ計算され得る。
【０２３６】
さらに具体的に説明すれば、例えば、力の相互作用を有する複数の物体がばねにより連結されているとみなし得る場合には、ばねの伸縮量に相当する相対変位に基づいて外力Ｆi が計算され、また、相互作用を有する複数の物体がばねとダンパとにより連結されているとみなし得る場合には、ばねの伸縮量に相当する相対変位と、ダンパの速度に相当する相対変位速度とに基づいて外力Ｆi が計算される。
【０２３７】
その後、Ｓ９において、前記式(44)を時間で積分することにより、時刻ｔ＝０から刻み時間ｄｔの経過時における物体の、モード次数がｊであるモード変位ｑj ，モード速度ｑ. j およびモード加速度ｑ: j が計算される。それらモード変位ｑj ，モード速度ｑ. j およびモード加速度ｑ: j は、剛体モードで生じるものと弾性モードで生じるものとの双方を含んでいる。前記式(44)を時間で積分するためにコンピュータにおいて使用される数値積分法の一例は、４次ルンゲ・クッタ法である。
【０２３８】
続いて、Ｓ１０において、局所観測座標系を更新するため、前述のようにして、その局所観測座標系の、次回の観測周期におけるスタート時点の位置ベクトルｒdt(0) ，速度ベクトルｖdt(0) および加速度ベクトルａdt(0) が計算される。その後、Ｓ１１において、局所観測座標系を更新するため、前述のようにして、その局所観測座標系の、次回の観測周期におけるスタート時点の方向余弦マトリクスＡdt(0) ，角速度ωの歪対称マトリクスω*dt (0) および角加速度αの歪対称マトリクスα*dt (0) が計算される。その後、Ｓ１３において、次回の観測に備えて、物体の、モード次数がｊである剛体モードのモード変位ｑj ，モード速度ｑ. j およびモード加速度ｑ: j のそれぞれの初期値が０に設定される。
【０２３９】
続いて、Ｓ１４において、現在の観測時刻ｔが、設定された終了時刻ｅＴに到達したか否かが判定される。今回は、到達してはいないと仮定すれば、判定がＮＯとなり、Ｓ１５において、観測時刻ｔが刻み時間ｄｔだけ増加させられ、その後、Ｓ８に戻る。そして、Ｓ８ないしＳ１３が前回と同様にして行われる。その後、それらＳ８ないしＳ１５が何回か繰り返された結果、現在の観測時刻ｔが終了時刻ｅＴに到達したと仮定すれば、Ｓ１４の判定がＹＥＳとなり、本プログラムの一回の実行が終了する。
【０２４０】
図８ないし図１０には、局所観測座標系ＬＯＣの運動（並進および回転）が物体Ｏの運動（剛体変位）との関係において概念的に示されている。
【０２４１】
まず、観測時刻ｔ＝０のときには、図８に示すように、局所観測座標系ＬＯＣが物体Ｏと同一の運動をさせられている。そのときから、刻み時間ｄｔが経過して観測時刻ｔ＝ｄｔとなるまでの間、すなわち、ある回の観測周期においては、局所観測座標系ＬＯＣは物体Ｏとは独立に運動させられる。したがって、観測時刻ｔ＝ｄｔのときに、局所観測座標系ＬＯＣの運動と物体Ｏの運動とは必ずしも一致しないが、その回の観測周期のスタート時点、すなわち、観測時刻ｔ＝０のときには、両者の運動が一致させられていたため、局所観測座標系ＬＯＣは常に物体Ｏの近傍に位置することとなり、観測時刻ｔ＝ｄｔのときに局所観測座標系ＬＯＣから物体Ｏを見ると、物体Ｏの剛体変位は微小となる。よって、前述の、モード解析において成立すべき前提が成立する。その後、図１０に示すように、次回の観測に備えて、局所観測座標系ＬＯＣの運動が、次回の観測周期のスタート時点、すなわち、今回の観測周期のエンド時点における物体Ｏと同じ運動をするように更新される。以上のような処理が、各観測周期ごとに行われる。
【０２４２】
以上、本実施形態において採用されている理論、ならびに本実施形態の構成を詳細に説明したが、次に、それら理論および構成の妥当性を、図１１に示すジャイロスコープを振動解析する場合を例にとり、説明する。
【０２４３】
１．解析モデルの説明
【０２４４】
上記ジャイロスコープは、質量ｍi （ｉ＝１，・・・，７）を有する７個の質点と、質量を無視できる長さｌ1 ，ｌ2 のはりと、質量を無視できる直径Ｄの円板とで構成されている。７個の質点は、正八面体の各頂点および重心に配置されている。このジャイロスコープは、質点７の位置において、全体座標系の原点に球面ジョイントを介して連結されており、また、−Ｚ方向に重量が作用している。このタイプのジャイロスコープは一般に、ラグランジュのこまと呼ばれている。
【０２４５】
２．局所観測座標系の定義および更新
【０２４６】
局所観測座標系は、初期状態において、図１１に示す座標系｛ｘ，ｙ，ｚ｝のように配置されており、かつ、等加速度および等角加速度で運動するようにされている。したがって、各観測周期（ｔ＝０〜ｄｔ）において、加速度ベクトルａおよび角加速度ベクトルαは、
【０２４７】
【数８４】
2005091370.tif 000086
【数８５】
2005091370.tif 000087
【０２４８】
なる２式でそれぞれ表される。したがって、速度ベクトルｖおよび位置ベクトルｒは、
【０２４９】
【数８６】
2005091370.tif 000088
【数８７】
2005091370.tif 000089
【０２５０】
なる２式でそれぞれ表される。また、同様に、角速度ベクトルωは、
【０２５１】
【数８８】
2005091370.tif 000090
【０２５２】
なる式で表される。
【０２５３】
方向余弦マトリクスＡは次のようにして計算される。
【０２５４】
まず、時刻０からｄｔまでの局所観測座標系の回転軸を表す単位ベクトルをｋ、回転角をθで表す。次に、その回転角θが微小であるとして、
【０２５５】
【数８９】
2005091370.tif 000091
【０２５６】
なる式が成立すると仮定する。この式における「ｋθ」が表す回転と等価な回転マトリクスをＲk-THで表し、ベクトルｋを、
【０２５７】
【数９０】
2005091370.tif 000092
【０２５８】
なる式で表すと、回転マトリクスＲk-THは、
【０２５９】
【数９１】
2005091370.tif 000093
【０２６０】
なる式で表される。ただし、ｃTH＝cos θ、ｓTH＝sin θ、ｖTH＝１−cos θである。そして、方向余弦マトリクスＡが、
【０２６１】
【数９２】
2005091370.tif 000094
【０２６２】
なる式として計算できる。
【０２６３】
なお、局所観測座標系は、前述の更新方法に従い、かつ、後述の数値積分におけると同じ刻み時間ｄｔごとに更新される。
【０２６４】
３．剛体変位が正しく計算できることの検証
【０２６５】
まず、本実施形態により剛体変位が正しく計算できることを確認するため、剛体モードのみ考慮して計算を行った。剛体モードのみ考慮すると、ジャイロスコープは質点７の位置において球面ジョイントを介して連結されているため、この運動系の自由度は３となる。そして、実際には、剛体モードとして、ｘ，ｙ，ｚ各軸まわりの３つの回転モードを選び、前記式(6) により正規化した固有モードベクトルを採用した。図１２には計算に必要な各種変数・定数の諸元、図１３には初期条件がそれぞれ示されている。なお、前記式(44)の数値積分法として４次ルンゲ・クッタ法を用いた。
【０２６６】
図１４において実線のグラフは、本実施形態により計算されたジャイロスコープの質点１の軌跡をＸＹ平面に正射影することによって得られたものである。また、同図において網かけ線のグラフは、比較のために、オイラー角を一般座標として用いたラグランジュの運動方程式を用いて計算された質点１の軌跡を同様にして示すものである。いずれの軌跡も、ジャイロスコーいずれの軌跡も、ジャイロスコープが、初期角速度と重力の影響により、章動を伴いながら歳差運動をする様子を示している。両者の計算結果は一致している。一方、オイラー角を一般座標として用いたラグランジュの運動方程式を用いた計算結果は、実際値と精度よく一致することが既に証明されている。このことから明らかなように、本実施形態は剛体変位を正しく解析できる物体変位・振動解析方法である。
【０２６７】
４．剛体変位と弾性振動が正しく計算できることの検証
【０２６８】
次に、上記ジャイロスコープにおける長さｌ1 ，ｌ2 のはりを弾性はりとして計算を行った。図１５には、弾性はりの諸元が示されている。また、前述のＭＳＣ−ＮＡＳＴＲＡＮにより固有値解析を行い、その結果に基づき、固有振動数が３４．７２７Ｈｚである円板の並進運動を伴う２つの曲げモードを採用した。さらに、初期条件として、図１３に示す条件に加えて、２つの弾性モードにそれぞれモード変位が０．００５である初期変位を与えた。
【０２６９】
図１６における実線のグラフは、本実施形態により計算した、弾性体ジャイロスコープの質点１の軌跡をＸＹ平面に正射影したものである。この軌跡は、ジャイロスコープが弾性振動と章動を伴いながら歳差運動をする様子を示している。同図における網かけ線のグラフは、剛体モードのみを採用した場合の軌跡を示しており、図１４における実線のグラフと同じものである。図１６において実線と網かけ線とを互いに比較してみると、弾性はりの弾性効果により、章動の大きさが若干変化していることが分かる。このことから明らかなように、本実施形態は剛体変位のみならず弾性振動をも正確に計算できる物体変位・振動解析方法である。
【０２７０】
５．全エネルギーが保存されているか否かの観点からの検証
【０２７１】
次に、本実施形態の妥当性をエネルギーの観点から検証するため、図１７には、ジャイロスコープの全エネルギーが時刻によって変化する履歴が示されている。同図における実線のグラフは、弾性体ジャイロスコープの全エネルギー、網かけ線のグラフは剛体ジャイロスコープの全エネルギーをそれぞれ示しており、図１６における実線のグラフと網かけ線のグラフとにそれぞれ対応している。図１７における両結果の差は、弾性ジャイロスコープに与えた初期変位の歪みエネルギー分に相当する。ジャイロスコープが弾性体である場合にも剛体である場合にも、全エネルギーが保存されていることが分かる。このことからも明らかなように、本実施形態は剛体変位のみならず弾性振動をも正確に計算できる物体変位・振動解析方法である。
【０２７２】
６．むすび
【０２７３】
以上の説明から明らかなように、本実施形態は、剛体変位および弾性振動の計算結果の観点ならびに全エネルギーの計算結果の観点からして、妥当であることが検証された。
【０２７４】
次に、本発明の別の実施形態である物体変位・振動解析方法を説明する。ただし、本実施形態は前記実施形態と共通する要素が多く、異なるのは局所観測座標系の運動を実現する要素についてのみであるため、共通する要素については詳細な説明を省略し、異なる要素についてのみ詳細に説明する。
【０２７５】
前記実施形態は、剛体としての動きが予め分かっていない物体について振動解析を行う場合の一例であり、この場合には、物体の剛体変位に合わせて局所観測座標系を逐次更新することが必要である。これに対して、物体の剛体としての動きが予め分かっている場合には、局所観測座標系にそのような動きを時間の関数として与えることが可能であり、このようにすれば、物体の剛体変位に合わせて局所観測座標系を更新することが不要となる。
【０２７６】
剛体としての動きが予め分かっている物体の一例は、図１８に示す自転車に装着されているチェーンである。このチェーンは、ペダルによって動くとともに、そのチェーンが巻き掛けられた一対の歯車対に沿って動くため、このチェーンの動きはある程度拘束される。したがって、そのような動きを、図１９に示すように、そのチェーンを構成する複数のリンクであって互いに回動可能に連結されたものの各々に局所観測座標系を設定するとともに、それら複数の局所観測座標系の各々に、対応するリンクと同じ動きを与える。これにより、局所観測座標系がチェーンの各リンクの少なくとも近傍に位置することが保証され、その結果、局所観測座標系を、各リンクの振動解析結果が得られるごとに更新することが不要となる。
【０２７７】
以上、本発明のいくつかの実施形態を図面に基づいて詳細に説明したが、これらは例示であり、本発明は、前記〔発明が解決しようとする課題，課題解決手段および発明の効果〕の項に記載された態様を始めとして、当業者の知識に基づいて種々の変形，改良を施した形態で実施することができる。
【図面の簡単な説明】
【０２７８】
【図１】本発明の一実施形態である物体変位・振動解析方法において使用されている全体座標系と局所観測座標系とを示す斜視図である。
【図２】従来の物体固定座標系の運動を説明するための正面図である。
【図３】図１における局所観測座標系の運動を説明するための正面図である。
【図４】図１に示す局所観測座標系が上記物体変位・振動解析方法により更新される様子を説明するための斜視図である。
【図５】上記物体変位・振動解析方法を実施するのに好適な物体変位・振動解析装置を示すブロック図である。
【図６】図５のＣＤ−ＲＯＭに記憶されている振動解析プログラムの前半部分を示すフローチャートである。
【図７】上記振動解析プログラムの後半部分を示すフローチャートである。
【図８】図６および図７の振動解析プログラムの実行内容を経時的に説明するための正面図である。
【図９】上記振動解析プログラムの実行内容を経時的に説明するための別の正面図である。
【図１０】上記振動解析プログラムの実行内容を経時的に説明するためのさらに別の正面図である。
【図１１】上記物体変位・振動解析方法の妥当性を検証するために用いた解析モデルを示す斜視図である。
【図１２】上記解析モデルの諸元を表形式で示す図である。
【図１３】上記解析モデルの初期条件を表形式で示す図である。
【図１４】上記解析モデルを剛体と仮定した場合の解析結果を上記物体変位・振動解析方法を用いた場合と、オイラー角を一般座標とする運動方程式を用いた場合とについて対比して示すグラフである。
【図１５】上記解析モデルを弾性体と仮定した場合の諸元を表形式で示す図である。
【図１６】上記物体変位・振動解析方法を用いた解析結果を、上記解析モデルを弾性体と仮定した場合と剛体と仮定した場合とについて対比して示すグラフである。
【図１７】上記物体変位・振動解析方法の妥当性を全エネルギーが保存されているか否かの観点から検証するためのグラフである。
【図１８】本発明の別の実施形態である物体変位・振動解析方法により振動解析が行われる物体としてのチェーンとそれを含む自転車とを示す正面図である。
【図１９】図１８におけるチェーンの一部を拡大するとともに概念化して示す正面図である。
【符号の説明】
【０２７９】
１２  コンピュータ
３０  ＣＤ−ＲＯＭ
ＬＯＣ  局所観測座標系
Ｏ  物体
【図１】
2005091370.tif 000095
【図２】
2005091370.tif 000096
【図３】
2005091370.tif 000097
【図４】
2005091370.tif 000098
【図５】
2005091370.tif 000099
【図６】
2005091370.tif 000100
【図７】
2005091370.tif 000101
【図８】
2005091370.tif 000102
【図９】
2005091370.tif 000103
【図１０】
2005091370.tif 000104
【図１１】
2005091370.tif 000105
【図１２】
2005091370.tif 000106
【図１３】
2005091370.tif 000107
【図１４】
2005091370.tif 000108
【図１５】
2005091370.tif 000109
【図１６】
2005091370.tif 000110
【図１７】
2005091370.tif 000111
【図１８】
2005091370.tif 000112
【図１９】
2005091370.tif 000113