JP3799437B2 - 物体変位・振動解析方法，装置および記録媒体 - Google Patents
物体変位・振動解析方法，装置および記録媒体 Download PDFInfo
- Publication number
- JP3799437B2 JP3799437B2 JP11267499A JP11267499A JP3799437B2 JP 3799437 B2 JP3799437 B2 JP 3799437B2 JP 11267499 A JP11267499 A JP 11267499A JP 11267499 A JP11267499 A JP 11267499A JP 3799437 B2 JP3799437 B2 JP 3799437B2
- Authority
- JP
- Japan
- Prior art keywords
- coordinate system
- mode
- displacement
- rigid body
- movable part
- 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 - Lifetime
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01H—MEASUREMENT OF MECHANICAL VIBRATIONS OR ULTRASONIC, SONIC OR INFRASONIC WAVES
- G01H1/00—Measuring characteristics of vibrations in solids by using direct conduction to the detector
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- Complex Calculations (AREA)
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
- Testing Of Engines (AREA)
Description
【０００１】
【発明の属する技術分野】
本発明は、全体座標系上における物体の、剛体としての変位である剛体変位と弾性体としての振動である弾性振動とを解析する技術に関するものである。
【０００２】
【従来の技術】
自動車のエンジン等、可動部品を有する機械においては、いかにして振動・騒音を低減させるかが重要な課題である。そのため、振動・騒音をより高い精度で解析できるより実用的な技術の確立が要望されている。機械要素単体の振動伝達特性に基づく振動解析については、汎用の構造解析ソフトウェア等を利用することにより、ほぼ実用レベルで実施できる段階に達していると言える。そして、近年では、複数の機械要素間の相互作用を考慮することにより、より高い精度で振動解析を行い得る技術の確立に重点が移行しつつある。
【０００３】
複数の機械要素間の相互作用を取り扱うためには、可動部品である物体の剛体変位と弾性振動との連成を解析することが重要な課題の一つである。そして、その課題を解決する一理論として機構解析による理論が既に知られている。この機構解析理論は、その当初は、複数の剛体が関節において互いに結合されたものを適用対象とするとともに、関節による結合を拘束条件式で表した微分代数方程式を用いる理論として構築された。その後、この機構解析理論は、物体の弾性振動を弾性モードの重ね合わせにより表現することにより、その適用対象を剛体から弾性体へと拡張された。この機構解析理論においては、物体の位置はデカルト座標、方向はオイラーパラメータ、弾性振動はモード座標によりそれぞれ表されている。モード座標は、弾性振動をモード解析する際に使用される、互いに直交する複数のモードベクトルにより定義されるものである。
【０００４】
【発明が解決しようとする課題，課題解決手段および発明の効果】
しかし、この機構解析理論に基づく従来の物体変位・振動解析方法では、物体の剛体変位と弾性振動とを記述する運動方程式において種類の異なる複数の座標系が混在する。このことは、剛性運動と弾性振動との連成を理由に、上記運動方程式において連成項の数を増加させたり、各連成項の構成を複雑にする可能性が高いことを意味する。そのため、この従来の物体変位・振動解析方法では、物体の剛体変位と弾性振動とを解析するための運動方程式の取り扱いが複雑で、計算に多くの工程や時間が必要であった。
【０００５】
さらに、この従来の物体変位・振動解析方法は、多関節の機構運動において低次の弾性振動の効果を考慮する場合には適しているが、弾性自由度が大きい弾性振動において剛体変位の効果を考慮する場合には十分に適しているとは言えなかった。
【０００６】
このような事情を背景とし、本発明は、エンジンの構成部材である可動部品の剛体変位と弾性振動とを解析する技術であって、弾性自由度が大きい弾性振動において剛体変位の効果を考慮する場合にも適しているものを提供することを課題としてなされたものであり、本発明によって、下記各態様のエンジン部品変位・振動解析方法，エンジン部品変位・振動解析用記録媒体およびエンジン部品変位・振動解析装置が得られる。各態様は、請求項と同様に、項に区分し、各項に番号を付し、必要に応じて他の項の番号を引用する形式で記載する。これは、本明細書に記載の技術的特徴およびそれらの組合わせのいくつかの理解を容易にするためであり、本明細書に記載の技術的特徴やそれらの組合わせが以下の態様に限定されると解釈されるべきではない。
【０００７】
（１）エンジンの構成部材である可動部品の、全体座標系上における剛体としての変位である剛体変位と、弾性体としての振動である弾性振動とをコンピュータにより解析する方法であって、
前記可動部品を有限要素法により解析すべく、その可動部品を複数の仮想的な要素に分割した有限要素モデルを前記可動部品の諸元に基づき作成するステップと、
前記有限要素モデルによりモード解析を行うために、前記全体座標系と、その全体座標系上を運動しつつ前記有限要素モデルを局所的に観測するための局所観測座標系であって、その有限要素モデルに固定されることなく有限要素モデルとは独立しているが少なくともその有限要素モデルの近傍に位置して運動を行うものとを設定するステップと、
それら設定した全体座標系および局所観測座標系を用いて、前記有限要素モデルの剛体変位と弾性振動とを、剛体変位は剛体モードにより、弾性振動は弾性モードによりそれぞれ表すとともに、それら剛体モードと弾性モードとを区別しないで一律に取り扱うモード解析により解析するステップと
を含み、かつ、前記全体座標系と前記局所観測座標系とを設定するステップが、前記局所観測座標系において前記可動部品を局所的に観測することを周期的に繰り返す観測周期の各々のスタート時点を除く期間において、局所観測座標系に、その局所観測座標系の変位と姿勢と速度と角速度と加速度と角加速度との中の少なくとも一つが、前記スタート時点における前記可動部品のそれと実質的に等しい大きさから変化しない運動を行わせるとともに、各観測周期のエンド時点において前記可動部品の解析結果が得られるごとに、局所観測座標系を、次回の観測に備えて、その可動部品の剛体変位に関する最新の解析結果と実質的に一致するように更新するステップを含むことを特徴とするエンジン可動部品の変位・振動解析方法（請求項１）。
この方法においては、全体座標系とは別に可動部品を観測する座標系として、その可動部品の少なくとも近傍に位置して運動を行う局所観測座標系が用いられる。この局所観測座標系は、可動部品に固定されてその可動部品の剛体変位と終始同じ運動をさせられる従来の物体固定座標系とは異なり、可動部品の有限要素モデルに固定されることなく有限要素モデルとは独立しているが、少なくとも有限要素モデルの近傍に位置して運動を行うものとされ、有限要素モデルを局所的に観測する。この局所的な観測が周期的に繰り返され、局所観測座標系は、各観測周期のスタート時点を除く期間においては、スタート時点における有限要素モデルの変位と姿勢と速度と角速度と加速度と角加速度との中の少なくとも一つが実質的に等しい大きさから変化しない運動を行わせられるとともに、各観測周期のエンド時点において有限要素モデルの解析結果が得られるごとに、局所観測座標系が、次回の観測に備えて、有限要素モデルの剛体変位に関する最新の解析結果と実質的に一致するように更新される。
この方法によれば、可動部品の剛体変位をモード解析する際の前提、すなわち、剛体変位がモード座標系上において微小であるという前提が容易に成立する。しかも、局所観測座標系の、各観測周期のスタート時点を除く期間における運動を比較的簡単に定義し得る。
よって、この方法によれば、小さな剛体変位を伴う可動部品のみならず大きな剛体変位を伴う可動部品についても容易に振動解析を行い得る。
ところで、例えば、自動車のエンジンの振動騒音を解析するためには、エンジンの回転中にそれのシリンダブロックに対してクランク軸が剛体変位すること，それらシリンダブロックの軸受とクランク軸のピンとの間に油膜が形成され、かつ、その油膜が非線形のばねとして作用すること等を考慮することが大切である。そして、この例においては、クランク軸がシリンダブロックに対して剛体変位する大きさは、そのクランク軸の弾性振動に比較すればかなり大きなものとなる。そして、本項に記載の方法は、この例の振動解析にも適用可能である。
また、本項に記載の方法は、単独の可動部品の振動解析に適用することが可能であり、また、上記の例、すなわち、シリンダブロックは車体に対して変位し、クランク軸はそのシリンダブロックに対して変位するエンジンの振動解析のように、相互作用を有する複数の可動部品の振動解析に適用することも可能である。
また、本項に記載の方法において「可動部品の少なくとも近傍に位置する」という要件は、それが成立すれば、可動部品の剛体変位をモード解析する際の前提、すなわち、剛体変位がモード座標系上において微小であるという前提が成立することになる要件であり、例えば、可動部品と局所観測座標系との相対変位（例えば、後述の個別変位ベクトルｕ’i ）が微小範囲にあるという要件であると定義したり、可動部品と局所観測座標系とが実質的に一致するという要件であると定義することが可能である。
また、本項に記載の方法において「可動部品の剛体変位」は、その可動部品の弾性振動の影響を受ける場合と、受けない場合とがある。受ける場合には、弾性振動の解析結果を考慮して剛体変位を決定することが必要であるが、受けない場合には、弾性振動の解析結果を待つことなくそれとは独立に剛体変位を決定することが可能である。よって、可動部品の剛体変位がその可動部品の弾性振動の影響を受ける場合には、可動部品の弾性振動を考慮して局所観測座標系の運動を定義することが必要であるが、受けない場合には、可動部品の弾性振動とは独立に局所観測座標系の運動を定義することが可能である。
（２）前記モード解析を、モード座標を一般座標として用いることによって前記可動部品の剛体変位と弾性振動とを記述する運動方程式を解くことにより行う(1)項に記載のエンジン可動部品の変位・振動解析方法。
この方法においては、モード座標を一般座標として用いることによって可動部品を記述する運動方程式を解くことにより、可動部品の剛体変位と弾性振動とが解析される。一方、局所観測座標系は、各観測周期のスタート時点を除く期間においては可動部品の剛体変位とは独立に、時間の関数に従って運動を行わせられる。したがって、可動部品を記述する運動方程式において、局所観測座標系の運動に関連する変数をモード座標に対して定数として扱うことが可能となり、その運動方程式を用いて可動部品の剛体変位と弾性振動とを計算することが容易になる。
ここに「運動方程式」の一例はラグランジュの運動方程式である。また「モード座標」は一般に、互いに交差する複数のモードベクトルにより定義される座標であると定義される。
（３）前記モード解析を、互いに直交する複数のモードベクトルにより定義されるモード座標を用いて行う(1)項または(2)項に記載のエンジン可動部品の変位・振動解析方法。
この方法においては、互いに直交する複数のモードベクトルにより定義されるモード座標を用いてモード解析が行われる。したがって、この方法によれば、可動部品を記述する運動方程式において剛体モードと弾性モードとの間に生じる連成項の数を少なくしたり、連成項の構成を簡単にすることが容易になり、ひいては、運動方程式を用いた計算が容易になる。
ここに「モードベクトル」は、固有モードベクトルを選ぶことができるのはもちろんであるが、その他のモードベクトルを選ぶことも可能である。例えば、静的補正モードベクトルを選ぶことも可能であり、その静的補正モードベクトルには例えば、アタッチメントモードベクトル，拘束モードベクトル等がある。
（４）前記モード解析を行うステップが、
(i)そのモード解析のための固有値λと、前記剛体モードを含むｎ個の正規直交化された複数のモードベクトルφj とを取得するステップと、
(ii)前記局所観測座標系により定義されたベクトルを前記全体座標系により定義されるベクトルに変換するための方向余弦マトリクスＡを計算するステップと、
(iii)最初の観測周期のスタート時点の前記全体座標系上における前記局所観測座標系の原点の位置ベクトルｒ，速度ベクトルｖおよび加速度ベクトルａを設定するステップと、(iv)各観測周期ごとに、前記可動部品に設定された複数の節点ｉの各々に作用する外力Ｆi を取得するステップと、
(v)各観測周期ごとに、その観測周期のエンド時点における前記可動部品のモード座標系上におけるモード次数がｊであるモード変位ｑj ，モード速度ｑ. jおよびモード加速度ｑ: j を計算するステップと、
(vi)前記各観測周期ごとに、次回の観測に備えて、その次回の観測周期のスタート時点の前記全体座標系上における前記局所観測座標系の原点の位置ベクトルｒdt(0) ，速度ベクトルｖdt(0) および加速度ベクトルａdt(0) を計算するステップと、
(vii)前記各観測周期ごとに、次回の観測に備えて、前記可動部品の、モード次数がｊである剛体モードのモード変位ｑj ，モード速度ｑ. j およびモード加速度ｑ: j をいずれも０に初期化するステップと
を含む(1)項ないし(3)項のいずれかに記載のエンジン可動部品の変位・振動解析方法（請求項２）。
（５）当該エンジン可動部品の変位・振動解析方法は、力の相互作用を有する複数の可動部品の剛体変位と弾性振動とを解析するものであり、
前記(iv)のステップが、各可動部品の節点ｉに作用する外力Ｆi を、関連する他の可動部品の節点ｉ’との間における力の相互作用に基づいて計算するとともに、その力の相互作用を、前記全体座標系上における前記複数の可動部品の、それらの弾性振動を考慮しない場合の各位置ベクトルｒi ，ｒi'と、各速度ベクトルｒ. i ，ｒ. i'と、各加速度ベクトルｒ: i ，ｒ: i'との少なくとも一つに基づいて計算するステップを含む(4)項に記載のエンジン可動部品の変位・振動解析方法（請求項３）。
この方法によれば、力の相互作用を有するエンジンの複数の可動部品の剛体変位と弾性振動とを解析することができる。
（６）前記(ii)のステップにおいて、前記可動部品に予め設定された複数の節点のうち、前記局所座標系の原点に対応する節点を基点として、互いに線形独立な３つのベクトルｓ’1 ，ｓ’2 ，ｓ’3 を生成する３つの節点を選択し、その生成された３つのベクトルｓ’1 ，ｓ’2 ，ｓ’3 に基づいて前記方向余弦マトリクスＡを計算する(4)項または(5)項に記載のエンジン可動部品の変位・振動解析方法（請求項４）。
この方法によれば、方向余弦マトリクスＡを簡単に計算することができる。
（７）前記各観測周期ごとに、次回の観測に備えて、その次回の観測周期のスタート時点における前記方向余弦マトリクスＡdt(0) ，前記局所観測座標系の角速度ωの歪対称マトリクスω*dt (0) およびその局所観測座標系の角加速度αの歪対称マトリクスα*dt (0) を計算し、計算した方向余弦マトリクスＡdt(0) ，歪対称マトリクスω*dt (0) および歪対称マトリクスα*dt (0) を、前記ｎ個の正規直交化された複数のモードベクトルφj により定義されるモード座標系により前記モード解析を行うことに起因する誤差が、方向余弦マトリクスＡdt(0) ，歪対称マトリクスω*dt (0) および歪対称マトリクスα*dt (0) に生じないように補正する(4)項ないし(6)項のいずれかに記載のエンジン可動部品の変位・振動解析方法（請求項５）。
（８）(1)項ないし(7)項のいずれかに記載のエンジン可動部品の変位・振動解析方法あるいは物体変位・振動解析方法を実施するためにコンピュータにより実行されるプログラムをコンピュータ読み取り可能に記録した記録媒体。
ここにおける「記録媒体」には例えば、フレキシブルディスク，磁気テープ，磁気ディスク，磁気ドラム，磁気カード，光ディスク，光磁気ディスク，ＲＯＭ，ＣＤ−ＲＯＭ，ＩＣカード，穿孔テープ等がある。
（９）コンピュータにより、エンジンの構成部材である可動部品の、全体座標系上における剛体としての変位である剛体変位と、弾性体としての振動である弾性振動とを解析するためのプログラムを記録した媒体であって、
前記プログラムは、
前記可動部品を有限要素法により解析すべく、その可動部品を複数の仮想的な要素に分割した有限要素モデルを前記可動部品の諸元に基づき作成するステップと、
前記有限要素モデルによりモード解析を行うために、前記全体座標系と、その全体座標系上を運動しつつ前記有限要素モデルを局所的に観測するための局所観測座標系であって、その有限要素モデルに固定されることなく有限要素モデルとは独立しているが少なくともその有限要素モデルの近傍に位置して運動を行うものとを設定するステップと、
それら全体座標系と局所観測座標系とを用いて、前記有限要素モデルの剛体変位と弾性振動とを、剛体変位は剛体モードにより、弾性振動は弾性モードによりそれぞれ表すとともに、それら剛体モードと弾性モードとを区別しないで一律に取り扱うモード解析により解析するステップと
を含み、かつ、前記全体座標系と前記局所観測座標系とを設定するステップが、前記局所観測座標系において前記可動部品を局所的に観測することを周期的に繰り返す観測周期の各々のスタート時点を除く期間において、局所観測座標系に、各その局所観測座標系の変位と姿勢と速度と角速度と加速度と角加速度との中の少なくとも一つが、前記スタート時点における前記可動部品のそれと実質的に等しい大きさから変化しない運動を行わせるとともに、各観測周期のエンド時点において前記可動部品の解析結果が得られるごとに、局所観測座標系を、次回の観測に備えて、その可動部品の剛体変位に関する最新の解析結果と実質的に一致するように更新するステップを含むことを特徴とする記録媒体（請求項６）。
前記 (2) 項ないし (7) 項の各々に記載の特徴は、本項の記録媒体にも適用可能である。
（１０）エンジンの構成部材である可動部品の、剛体としての変位である剛体変位と、弾性体としての振動である弾性振動とを解析するための装置であって、
前記可動部品を有限要素法により解析すべく、その可動部品を複数の仮想的な要素に分割した有限要素モデルを前記可動部品の諸元に基づき作成する手段と、
前記有限要素モデルによりモード解析を行うために、前記全体座標系と、その全体座標系上を運動しつつ前記有限要素モデルを局所的に観測するための局所観測座標系であって、その有限要素モデルに固定されることなく有限要素モデルとは独立しているが少なくともその有限要素モデルの近傍に位置して運動を行うものとを設定する手段と、
それら設定した全体座標系および局所観測座標系を用いて、前記有限要素モデルの剛体変位と弾性振動とを、剛体変位は剛体モードにより、弾性振動は弾性モードによりそれぞれ表すとともに、それら剛体モードと弾性モードとを区別しないで一律に取り扱うモード解析により解析する手段と
を含み、かつ、前記全体座標系と前記局所観測座標系とを設定する手段が、前記局所観測座標系において前記可動部品を局所的に観測することを周期的に繰り返す観測周期の各々のスタート時点を除く期間において、局所観測座標系に、その局所観測座標系の変位と姿勢と速度と角速度と加速度と角加速度との中の少なくとも一つが、前記スタート時点における前記可動部品のそれと実質的に等しい大きさから変化しない運動を行わせるとともに、各観測周期のエンド時点において前記可動部品の解析結果が得られるごとに、局所観測座標系を、次回の観測に備えて、その可動部品の剛体変位に関する最新の解析結果と実質的に一致するように更新する手段を含むことを特徴とするエンジン可動部品の変位・振動解析装置（請求項７）。
前記 (2)項ないし(7)項の各々に記載の特徴は、本項のエンジン可動部品の剛体変位・振動解析装置にも適用可能である。
【０００８】
【発明の実施の形態】
以下、本発明のさらに具体的な一実施形態を図面に基づいて詳細に説明する。
【０００９】
本実施形態は、物体の剛体変位と弾性振動とを解析する物体変位・振動解析方法である。この方法は、剛体自由度は小さいが弾性自由度は大きい機械の振動騒音の解析に適している。
【００１０】
本実施形態は、種々の新たな理論に基づいて構成されており、以下、その新たな理論を説明する。
【００１１】
１．局所観測座標系の導入
【００１２】
図１に示すように、全体座標系｛Ｘ，Ｙ，Ｚ｝に対し、局所観測座標系｛ｘ，ｙ，ｚ｝なる新しい概念の座標系を導入する。この座標系は、図２に示すように、従来の、物体に固定された座標系とは異なり、物体の近傍に位置しながら、物体の運動を局所的に観測するために設けられた座標系である。局所観測座標系ＬＯＣは、図３に示すように、ある時刻に物体Ｏと同一の運動をし、その状態から、物体Ｏの運動とは独立に、時間の関数に従って運動をする。
【００１３】
ここで、局所観測座標系と物体との関係を、位置，速度および加速度に関してそれぞれ説明する。
【００１４】
まず、図１に示すように、全体座標系から見た局所観測座標系の原点の位置ベクトルをｒ、方向余弦マトリクスをＡ、局所観測座標系から見た物体の節点ｉの位置ベクトルをｒ’i とすると、全体座標系から見た物体の節点ｉの位置ベクトルｒi は、
【００１５】
【数１】
【００１６】
なる式で表される。
【００１７】
ここで補足説明すれば、上記方向余弦マトリクスＡは、座標系を局所観測座標系から全体座標系に変換する座標変換マトリクスの一種である。また、上記位置ベクトルｒ’i において「’」は、そのベクトルが局所観測座標系により定義されることを示す表記であり、この表記は他の変数を表記する際にも用いられる。また、観測すべき物体には、その物体を複数の要素に仮想的に分割するためにＮ個の節点が設定されており、上記節点ｉはそれらＮ個の節点のうち任意のものを示している。
【００１８】
物体の剛体変位および弾性振動を考慮しない場合の、局所観測座標系から見た節点ｉの位置ベクトルをｓ’i とし、それら剛体変位と弾性振動との和を個別変位ベクトルｕ’i とすると、
【００１９】
【数２】
【００２０】
なる式が成立する。ここで、Ｎ個の節点の全変位ベクトルｕ’を、
【００２１】
【数３】
【００２２】
なる式で定義する。ここで、全変位ベクトルｕ’が十分小さく、モード解析が前提とする微小変位の範囲内にあり、かつ、全変位ベクトルｕ’が、弾性モードのみならず剛体モードをも含むｎ個の正規直交化された固有モードベクトルφj （モード次数ｊ＝１，・・・，ｎ）の重ね合わせにより表すことができるとすると、全変位ベクトルｕ’が、モード変位ｑj を用いることにより、
【００２３】
【数４】
【００２４】
なる式で表される。ここで、モードマトリクスΦを、固有モードベクトルφj を用いることにより、
【００２５】
【数５】
【００２６】
なる式で定義すると、モードマトリクスΦは、固有モードの直交性を理由に、質量マトリクスＭと剛性マトリクスＫとの間において、
【００２７】
【数６】
【数７】
【００２８】
なる式で表される関係を有する。ここで、「Ｉ」は単位マトリクスであり、「λ」は、固有値λj を対角成分として有するマトリクス、すなわち、
【００２９】
【数８】
【００３０】
である。ただし、モード次数ｊ＝１，・・・，ｎr が剛体モードを表すものであるとすると、固有値λi ＝０（ｊ＝１，・・・，ｎr ）となる。ここで、モード変位ベクトルｑを、マトリクスの転置を表す表記として「T 」を用いることにより、
【００３１】
【数９】
【００３２】
なる式で定義すると、前記式(4) は、前記式(5) を用いることにより、
【００３３】
【数１０】
【００３４】
に変形される。前記式(3) における全変位ベクトルｕ’と個別変位ベクトルｕ’i との関係と同様に、固有モードベクトルφj を、各節点ｉごとに分けてφijで表すことにより、
【００３５】
【数１１】
【００３６】
なる式で定義すると、前記式(2) における個別変位ベクトルｕ’i は、
【００３７】
【数１２】
【００３８】
なる式で表される。したがって、前記式(1) は、前記式(2) と上記式(12)とを用いることにより、
【００３９】
【数１３】
【００４０】
に変形される。この式(13)における「ｒi 」は、図１に示すものと同じである。また、この式(13)は、局所観測座標系と物体との間の、位置に関する関係式を意味する。
【００４１】
上記式(13)を時間で微分すると、局所観測座標系と物体との間の、速度に関する関係式として、
【００４２】
【数１４】
【００４３】
なる式が得られる。ここに、「Ａ. 」は、方向余弦マトリクスＡを時間で微分したものを意味し、「ｑ. j 」は、モード変位ｑj を時間で微分したものを意味する。
【００４４】
ここで、全体座標系に対する局所観測座標系の角速度ベクトルをω、その歪対称マトリクスをω* で表すと、
【００４５】
【数１５】
【００４６】
なる式で表される関係が成立する。この式(15)を上記式(14)に代入し、全体座標系に対する局所観測座標系の速度ベクトルを改めてｖで表すと、速度の関係式として、
【００４７】
【数１６】
【００４８】
なる式が得られる。
【００４９】
前記式(14)を時間で微分すると、局所観測座標系と物体との間の、加速度に関する関係式として、
【００５０】
【数１７】
【００５１】
なる式が得られる。ここに、「Ａ: 」は、前記Ａ. を時間で微分したものを意味し、「ｑ: j 」は、前記ｑ. j を時間で微分したものを意味する。一方、前記式(15)を時間で微分すると、
【００５２】
【数１８】
【００５３】
なる式が得られる。ここに、「ω.*」は、歪対称マトリクスω* を時間で微分したものを意味する。この式(18)に前記式(15)を代入するとともに、全体座標系に対する局所観測座標系の角加速度ベクトルをα、その歪対称マトリクスをα* としてそれぞれ表すと、
【００５４】
【数１９】
【００５５】
なる式が得られる。この式(19)と前記式(15)とを前記式(17)に代入するとともに、全体座標系に対する局所観測座標系の加速度ベクトルを改めてａで表すと、
【００５６】
【数２０】
【００５７】
なる式が得られる。
【００５８】
以上のようにして、局所観測座標系｛ｘ，ｙ，ｚ｝を導入することにより、前記式(13)，式(16)および式(20)が、弾性振動のみならず剛体変位をも伴う物体の節点ｉの位置，速度および加速度を、全体座標系と局所観測座標系とモード座標系とを用いて記述する式として得られた。
【００５９】
なお付言すれば、それら３式において全体座標系または局所観測座標系により定義される変数は、後述のように、物体の剛体変位とは独立に、時間の関数として定義される。すなわち、その変数は、時間依存変数と言うことができるのである。一方、本実施形態において採用された物体の運動方程式は、後述のように、モード座標を一般座標として用いる。したがって、その運動方程式において上記時間依存変数は、一般座標との関係においては定数として取り扱うことが可能となる。このような観点からすれば、前記３式は、局所観測座標系の導入により、弾性振動のみならず剛体変位をも伴う物体の節点ｉの位置，速度および加速度を、モード座標系のみを用いて上記運動方程式に組み込むことを可能にした式であると言うことができ、また、それら位置，速度および加速度を記述する座標系を、全体座標系および局所観測座標系からモード座標系に変換することを可能にした式であると言うこともできる。
【００６０】
２．モード座標系による運動方程式の導出
【００６１】
物体全体の運動エネルギーＴは、節点ｉの質量をｍi で表すと、
【００６２】
【数２１】
【００６３】
なる式で表される。また、弾性変形による歪みエネルギーＵは、前記式(3) の全変位ベクトルｕ’を用いることにより、
【００６４】
【数２２】
【００６５】
なる式で表される。この式(22)に前記式(10)を代入すると、
【００６６】
【数２３】
【００６７】
なる式が得られ、さらに、前記式(7) により表される関係を利用すると、
【００６８】
【数２４】
【００６９】
なる式が得られる。
【００７０】
一方、ラグランジュの運動方程式を、モード座標ｑk （モード次数ｋ＝１，・・・，ｎr ）を一般座標（generalized coordinates, general coordenates）として用いて表すと、
【００７１】
【数２５】
【００７２】
なる式が得られる。この式(25)は、時間に依存しないため、局所観測座標系を、物体の運動とは独立に、時間の関数に従って運動させることによる影響を何ら受けない。したがって、この式(25)は、局所観測座標系の運動を考慮することなく、解くことができるのである。
【００７３】
上記式(25)の右辺における「Ｑk 」は、モード座標系での一般化力（generalized force ）を表している。
【００７４】
左辺の第１項の括弧内の変数は、前記式(21)を用いることにより、
【００７５】
【数２６】
【００７６】
で表される。この式(26)は、前記式(16)を用いることにより、
【００７７】
【数２７】
【００７８】
なる式に変形され、転置して整理すると、
【００７９】
【数２８】
【００８０】
なる式が得られる。この式(28)を時間で微分し、前記式(15)を用いて整理すると、
【００８１】
【数２９】
【００８２】
なる式が得られる。
【００８３】
一方、前記式(25)の左辺の第２項は、前記式(21)を用いることにより、
【００８４】
【数３０】
【００８５】
で表される。この式(30)は、前記式(16)を用いることにより、
【００８６】
【数３１】
【００８７】
に変形され、転置して整理すると、
【００８８】
【数３２】
【００８９】
なる式が得られる。この式(32)の右辺の項は、前記式(29)の右辺の第２項と同一であるから、それら２式を引き算すれば、
【００９０】
【数３３】
【００９１】
なる式が得られる。この式(33)に前記式(20)を代入すると、
【００９２】
【数３４】
【００９３】
なる式が得られる。ここに、方向余弦マトリクスＡは正規直交マトリクスであるため、
【００９４】
【数３５】
【００９５】
なる式が成立する。また、前記式(6) により表される、固有モードの直交性という条件は、節点ｉの質量ｍi を用いて表すと、
【００９６】
【数３６】
【００９７】
となる。ただし、「δij」は、クロネッカーのデルタである。この式(36)と前記式(35)を用いると、前記式(34)は、
【００９８】
【数３７】
【００９９】
なる式に変形される。
【０１００】
一方、前記式(25)の左辺の第３項は、前記式(24)と前記式(8) とを用いることにより、
【０１０１】
【数３８】
【０１０２】
なる式で表される。
【０１０３】
また、前記式(25)の右辺の一般化力Ｑk については、まず、全体座標系から見た節点ｉに作用する外力をＦi として、仮想仕事δＷを考えるととにも、
【０１０４】
【数３９】
【０１０５】
なる式を用いる。ここに、「δｒi 」は、前記位置ベクトルｒi の仮想的な微小変化を表す。この式(39)は、前記式(13)を用いることにより、
【０１０６】
【数４０】
【０１０７】
なる式に変形される。ここに、「δｑ j 」は、前記モード変位ｑ j の仮想的な微小変化を表す。ここで、微小変化δｑ j の係数をまとめることにより、
【０１０８】
【数４１】
【０１０９】
なる式を定義すると、一般化力Ｑk は、
【０１１０】
【数４２】
【０１１１】
なる式で表され、転置して整理すると、
【０１１２】
【数４３】
【０１１３】
なる式が得られる。
【０１１４】
前記式(37)，式(38)および式(43)を前記式(25)、すなわち、ラグランジュの運動方程式に代入して整理すると、最終的に、モード座標系により記述された運動方程式として、
【０１１５】
【数４４】
【０１１６】
なる式が得られる。
【０１１７】
この式(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)を用い、かつ、原点（ｉ＝０）であることを考慮することにより、
【０１２８】
【数４５】
【０１２９】
なる式で表されるように更新される。ここに、「φ0j」は、物体の複数の節点のうち、局所観測座標系の原点に対応する節点の固有モードベクトルを示している。なお、「原点に対応する節点」の意味については後述する。
【０１３０】
(3) 局所観測座標系の原点の速度ベクトルｖdt(0) に関する更新
【０１３１】
時刻ｄｔにおける新しい局所観測座標系の原点の速度ベクトルｖdt(0) は、前記式(16)を用いることにより、
【０１３２】
【数４６】
【０１３３】
なる式で表されるように更新される。
【０１３４】
(4) 局所観測座標系の原点の加速度ベクトルａdt(0) に関する更新
【０１３５】
時刻ｄｔにおける新しい局所観測座標系の原点の加速度ベクトルａdt(0) は、前記式(20)を用いることにより、
【０１３６】
【数４７】
【０１３７】
なる式で表されるように更新される。
【０１３８】
(5) 方向余弦マトリクスＡdt(0) の更新
【０１３９】
物体の剛体変位を考慮せず、かつ、時刻ｄｔである場合の、全体座標系から見た節点ｉの位置ベクトルｒdt i は、前記式(13)を用いることにより、
【０１４０】
【数４８】
【０１４１】
なる式で表される。また、時刻ｄｔに、前記位置ベクトルｓ’i を全体座標系から見た場合のベクトルｓdt i は、
【０１４２】
【数４９】
【０１４３】
なる式で表され、この式(49)に前記式(45)と式(48)とを代入すると、
【０１４４】
【数５０】
【０１４５】
なる式が得られる。また、上記ベクトルｓdt i は、時刻ｄｔにおける方向余弦マトリクスＡdt(0) を用いることにより、
【０１４６】
【数５１】
【０１４７】
なる式で表される。この式(51)を方向余弦マトリクスＡdt(0) について解くために、図４に示すように、物体の複数の節点の中から、局所観測座標系の原点と共同して四面体を構成する３つの節点であって、かつ、原点を始点、各節点を終点とする３つのベクトルｓ’1 ，ｓ’2 およびｓ’3 が互いに線形独立であるものを選択する。この選択により、
【０１４８】
【数５２】
【０１４９】
なる式が得られる。ここに、「ｓ’1 」，「ｓ’2 」および「ｓ’3 」はそれぞれ、同図において、原点を始点、節点１，２および３をそれぞれ終点とする、互いに線形独立な３つのベクトルを意味している。このように、ベクトルｓ’1 ，ｓ’2 およびｓ’3 は互いに線形独立であるため、それらベクトルｓ’1 ，ｓ’2 およびｓ’3 を成分とするマトリクスには逆マトリクスが存在する。よって、上記式(52)における方向余弦マトリクスＡdt(0) は、
【０１５０】
【数５３】
【０１５１】
なる式で表されるように更新することができる。
【０１５２】
(6) 角速度の歪対称マトリクスω*dt (0) の更新
【０１５３】
時刻ｄｔにおける節点ｉの、物体の剛体変位に起因する速度成分ｒ.dt i は、前記式(16)を用いることにより、
【０１５４】
【数５４】
【０１５５】
なる式で表される。ベクトルｓdt i を時間で微分して得られるベクトルｓ.dt i は、全体座標系に対して並進方向の速度成分を持たないため、
【０１５６】
【数５５】
【０１５７】
なる式で表される。この式(55)に前記式(46)と上記式(54)とを代入すると、
【０１５８】
【数５６】
【０１５９】
なる式が得られる。さらに、ベクトルｓ.dt i は、前記式(53)と時刻ｄｔにおける歪対称マトリクスω*dt (0) とを用いることにより、
【０１６０】
【数５７】
【０１６１】
なる式で表される。したがって、前記式(53)と同様に、前述の、互いに線形独立なベクトルｓ’1 ，ｓ’2 およびｓ’3 を用いることにより、
【０１６２】
【数５８】
【０１６３】
なる式が求められる。よって、角速度の歪対称マトリクスω*dt (0) は、
【０１６４】
【数５９】
【０１６５】
なる式で表されるように更新される。
【０１６６】
(7) 角加速度の歪対称マトリクスα*dt (0) の更新
【０１６７】
時刻ｄｔにおける節点ｉの、物体の剛体変位に起因する加速度成分ｒ:dt i は、前記式(20)を用いることにより、
【０１６８】
【数６０】
【０１６９】
なる式で表される。ここに、「ｓ:dt i 」は、全体座標系に対して並進方向の加速度成分を持たないため、
【０１７０】
【数６１】
【０１７１】
なる式で表される。この式(61)に前記式(47)と上記式(60)とを代入することにより、
【０１７２】
【数６２】
【０１７３】
なる式が得られる。さらに、ベクトルｓ:dt i は、前記式(53)と式(59)と時刻ｄｔにおける角加速度の歪対称マトリクスα*dt (0) とを用いることにより、
【０１７４】
【数６３】
【０１７５】
なる式で表される。よって、前記式(53)および式(58)と同様に、前述の、互いに線形独立なベクトルｓ’1 ，ｓ’2 およびｓ’3 を用いることにより、
【０１７６】
【数６４】
【０１７７】
なる式が求められる。よって、角加速度の歪対称マトリクスα*dt (0) は、
【０１７８】
【数６５】
【０１７９】
なる式で表されるように更新される。
【０１８０】
(8) 方向余弦マトリクスＡdt (0)の補正
【０１８１】
本実施形態においては、物体の剛体変位を剛体モードの重ね合わせにより解析し、その際、複数の固有モードベクトルが互いに直交するモード座標系を用いて物体の剛体変位を記述するため、前記式(53)により求めた方向余弦マトリクスＡdt (0)には微小な誤差が存在する。以下、この誤差の存在を具体的に説明する。
【０１８２】
全体座標系から見た局所観測座標系の座標軸ベクトルをｘ< ，ｙ< ，ｚ< で表すこととすると、方向余弦マトリクスＡは、
【０１８３】
【数６６】
【０１８４】
で表すことができる。各座標軸ベクトルは単位ベクトルであるため、
【０１８５】
【数６７】
【０１８６】
なる式が成立する。さらに、３つの座標軸ベクトルは互いに直交するため、
【０１８７】
【数６８】
【０１８８】
なる式が成立する。
【０１８９】
しかし、前記式(53)により更新された方向余弦マトリクスＡdt (0)の各列ベクトルは、小さい誤差を含んだｘ- ，ｙ- ，ｚ- であるため、上記式(67)により表される関係および上記式(68)により表される関係が厳密には成立しない。
【０１９０】
そこで、本実施形態においては、その誤差が解消されるように、方向余弦マトリクスＡdt (0)が適切に補正される。以下、その補正方法を説明する。
【０１９１】
まず、列ベクトルｘ- の方向は正しいと仮定として、列ベクトルｘ< を、
【０１９２】
【数６９】
【０１９３】
なる式を用いて決定される。
【０１９４】
次に、求めるべき列ベクトルｙ< は、上記のようにして求められた列ベクトルｘ< と列ベクトルｙ- と同一平面上にあるとして、列ベクトルｚ< が、
【０１９５】
【数７０】
【０１９６】
なる式を用いて決定される。
【０１９７】
最後に、列ベクトルｙ< が、
【０１９８】
【数７１】
【０１９９】
なる式を用いて決定される。
【０２００】
(9) 角速度の歪対称マトリクスω*dt (0) の補正
【０２０１】
方向余弦マトリクスＡdt (0)の場合と同様に、更新された角速度ωの歪対称マトリクスω*dt (0) にも微小な誤差が含まれる。そこで、本実施形態においては、歪対称マトリクスω*dt (0) も適切に補正される。以下、その補正方法を説明する。
【０２０２】
前記式(59)で求めた歪対称マトリクスω*dt (0) は、一般にフルマトリクスとなるが、本来は、角速度ベクトルω、すなわち、
【０２０３】
【数７２】
【０２０４】
の歪対称マトリクスであって、
【０２０５】
【数７３】
【０２０６】
なる式で表される関係を満たすはずである。ここで、誤差を含んだ歪対称マトリクスω*dt (0) の各成分を、
【０２０７】
【数７４】
【０２０８】
なる式で表す。歪対称マトリクスω*dt (0) の補正方法は種々考えられるが、本実施形態においては、次のように行われる。
【０２０９】
すなわち、まず、歪対称マトリクスω*dt (0) の対角成分が０でない場合にはすべて誤差であるとして、対角成分ω * ii ＝０（ｉ＝１，２，３）とされる。次に、基本成分ωx ，ωy ，ωz が、
【０２１０】
【数７５】
【数７６】
【数７７】
【０２１１】
のように補正される。
【０２１２】
(10)角加速度の歪対称マトリクスα*dt (0) の補正
【０２１３】
歪対称マトリクスα*dt (0) の補正方法は、歪対称マトリクスω*dt (0) の場合と基本的に同じである。以下、具体的に説明する。
【０２１４】
前記式(65)により更新された歪対称マトリクスα*dt (0) は、一般にフルマトリクスとなるが、本来は、角加速度ベクトル、すなわち、
【０２１５】
【数７８】
【０２１６】
の歪対称マトリクスであるため、
【０２１７】
【数７９】
【０２１８】
なる式で表される関係を満たすはずである。ここで、誤差を含んだ歪対称マトリクスα*dt (0) の各成分を、
【０２１９】
【数８０】
【０２２０】
なる式で表す。そして、歪対称マトリクスα*dt (0) の補正においては、まず、歪対称マトリクスα*dt (0) の対角成分が０でない場合にはすべて誤差であるとして、対角成分α * ii ＝０（ｉ＝１，２，３）とされる。次に、基本成分αx ，αy ，αz が、
【０２２１】
【数８１】
【数８２】
【数８３】
【０２２２】
のように補正される。
【０２２３】
(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 のはりと、質量を無視できる直径Ｄの円板とで構成されている。７個の質点は、正八面体の各頂点および重心に配置されている。このジャイロスコープは、質点７の位置において、全体座標系の原点に球面ジョイントを介して連結されており、また、−Ｚ方向に重量が作用している。このタイプのジャイロスコープは一般に、ラグランジュのこまと呼ばれている。
【０２４５】
２．局所観測座標系の定義および更新
【０２４６】
局所観測座標系は、初期状態において、図１１に示す座標系｛ｘ，ｙ，ｚ｝のように配置されており、かつ、等加速度および等角加速度で運動するようにされている。したがって、各観測周期（ｔ＝０〜ｄｔ）において、加速度ベクトルａおよび角加速度ベクトルαは、
【０２４７】
【数８４】
【数８５】
【０２４８】
なる２式でそれぞれ表される。したがって、速度ベクトルｖおよび位置ベクトルｒは、
【０２４９】
【数８６】
【数８７】
【０２５０】
なる２式でそれぞれ表される。また、同様に、角速度ベクトルωは、
【０２５１】
【数８８】
【０２５２】
なる式で表される。
【０２５３】
方向余弦マトリクスＡは次のようにして計算される。
【０２５４】
まず、時刻０からｄｔまでの局所観測座標系の回転軸を表す単位ベクトルをｋ、回転角をθで表す。次に、その回転角θが微小であるとして、
【０２５５】
【数８９】
【０２５６】
なる式が成立すると仮定する。この式における「ｋθ」が表す回転と等価な回転マトリクスをＲk-THで表し、ベクトルｋを、
【０２５７】
【数９０】
【０２５８】
なる式で表すと、回転マトリクスＲk-THは、
【０２５９】
【数９１】
【０２６０】
なる式で表される。ただし、ｃTH＝cos θ、ｓTH＝sin θ、ｖTH＝１−cos θである。そして、方向余弦マトリクスＡが、
【０２６１】
【数９２】
【０２６２】
なる式として計算できる。
【０２６３】
なお、局所観測座標系は、前述の更新方法に従い、かつ、後述の数値積分におけると同じ刻み時間ｄｔごとに更新される。
【０２６４】
３．剛体変位が正しく計算できることの検証
【０２６５】
まず、本実施形態により剛体変位が正しく計算できることを確認するため、剛体モードのみ考慮して計算を行った。剛体モードのみ考慮すると、ジャイロスコープは質点７の位置において球面ジョイントを介して連結されているため、この運動系の自由度は３となる。そして、実際には、剛体モードとして、ｘ，ｙ，ｚ各軸まわりの３つの回転モードを選び、前記式(6) により正規化した固有モードベクトルを採用した。図１２には計算に必要な各種変数・定数の諸元、図１３には初期条件がそれぞれ示されている。なお、前記式(44)の数値積分法として４次ルンゲ・クッタ法を用いた。
【０２６６】
図１４において実線のグラフは、本実施形態により計算されたジャイロスコープの質点１の軌跡をＸＹ平面に正射影することによって得られたものである。また、同図において網かけ線のグラフは、比較のために、オイラー角を一般座標として用いたラグランジュの運動方程式を用いて計算された質点１の軌跡を同様にして示すものである。いずれの軌跡も、ジャイロスコーいずれの軌跡も、ジャイロスコープが、初期角速度と重力の影響により、章動を伴いながら歳差運動をする様子を示している。両者の計算結果は一致している。一方、オイラー角を一般座標として用いたラグランジュの運動方程式を用いた計算結果は、実際値と精度よく一致することが既に証明されている。このことから明らかなように、本実施形態は剛体変位を正しく解析できる物体変位・振動解析方法である。
【０２６７】
４．剛体変位と弾性振動が正しく計算できることの検証
【０２６８】
次に、上記ジャイロスコープにおける長さｌ1 ，ｌ2 のはりを弾性はりとして計算を行った。図１５には、弾性はりの諸元が示されている。また、前述のＭＳＣ−ＮＡＳＴＲＡＮにより固有値解析を行い、その結果に基づき、固有振動数が３４．７２７Ｈｚである円板の並進運動を伴う２つの曲げモードを採用した。さらに、初期条件として、図１３に示す条件に加えて、２つの弾性モードにそれぞれモード変位が０．００５である初期変位を与えた。
【０２６９】
図１６における実線のグラフは、本実施形態により計算した、弾性体ジャイロスコープの質点１の軌跡をＸＹ平面に正射影したものである。この軌跡は、ジャイロスコープが弾性振動と章動を伴いながら歳差運動をする様子を示している。同図における網かけ線のグラフは、剛体モードのみを採用した場合の軌跡を示しており、図１４における実線のグラフと同じものである。図１６において実線と網かけ線とを互いに比較してみると、弾性はりの弾性効果により、章動の大きさが若干変化していることが分かる。このことから明らかなように、本実施形態は剛体変位のみならず弾性振動をも正確に計算できる物体変位・振動解析方法である。
【０２７０】
５．全エネルギーが保存されているか否かの観点からの検証
【０２７１】
次に、本実施形態の妥当性をエネルギーの観点から検証するため、図１７には、ジャイロスコープの全エネルギーが時刻によって変化する履歴が示されている。同図における実線のグラフは、弾性体ジャイロスコープの全エネルギー、網かけ線のグラフは剛体ジャイロスコープの全エネルギーをそれぞれ示しており、図１６における実線のグラフと網かけ線のグラフとにそれぞれ対応している。図１７における両結果の差は、弾性ジャイロスコープに与えた初期変位の歪みエネルギー分に相当する。ジャイロスコープが弾性体である場合にも剛体である場合にも、全エネルギーが保存されていることが分かる。このことからも明らかなように、本実施形態は剛体変位のみならず弾性振動をも正確に計算できる物体変位・振動解析方法である。
【０２７２】
６．むすび
【０２７３】
以上の説明から明らかなように、本実施形態は、剛体変位および弾性振動の計算結果の観点ならびに全エネルギーの計算結果の観点からして、妥当であることが検証された。
【０２７７】
以上、本発明の一実施形態を図面に基づいて詳細に説明したが、これは文字通り例示であり、本発明は、前記〔発明が解決しようとする課題，課題解決手段および発明の効果〕の項に記載された態様を始めとして、当業者の知識に基づいて種々の変形，改良を施した形態で実施することができる。
【図面の簡単な説明】
【図１】本発明の一実施形態である物体変位・振動解析方法において使用されている全体座標系と局所観測座標系とを示す斜視図である。
【図２】従来の物体固定座標系の運動を説明するための正面図である。
【図３】図１における局所観測座標系の運動を説明するための正面図である。
【図４】図１に示す局所観測座標系が上記物体変位・振動解析方法により更新される様子を説明するための斜視図である。
【図５】上記物体変位・振動解析方法を実施するのに好適な物体変位・振動解析装置を示すブロック図である。
【図６】図５のＣＤ−ＲＯＭに記憶されている振動解析プログラムの前半部分を示すフローチャートである。
【図７】上記振動解析プログラムの後半部分を示すフローチャートである。
【図８】図６および図７の振動解析プログラムの実行内容を経時的に説明するための正面図である。
【図９】上記振動解析プログラムの実行内容を経時的に説明するための別の正面図である。
【図１０】上記振動解析プログラムの実行内容を経時的に説明するためのさらに別の正面図である。
【図１１】上記物体変位・振動解析方法の妥当性を検証するために用いた解析モデルを示す斜視図である。
【図１２】上記解析モデルの諸元を表形式で示す図である。
【図１３】上記解析モデルの初期条件を表形式で示す図である。
【図１４】上記解析モデルを剛体と仮定した場合の解析結果を上記物体変位・振動解析方法を用いた場合と、オイラー角を一般座標とする運動方程式を用いた場合とについて対比して示すグラフである。
【図１５】上記解析モデルを弾性体と仮定した場合の諸元を表形式で示す図である。
【図１６】上記物体変位・振動解析方法を用いた解析結果を、上記解析モデルを弾性体と仮定した場合と剛体と仮定した場合とについて対比して示すグラフである。
【図１７】上記物体変位・振動解析方法の妥当性を全エネルギーが保存されているか否かの観点から検証するためのグラフである。
【符号の説明】
１２ コンピュータ
３０ ＣＤ−ＲＯＭ
ＬＯＣ 局所観測座標系
Ｏ 物体
【発明の属する技術分野】
本発明は、全体座標系上における物体の、剛体としての変位である剛体変位と弾性体としての振動である弾性振動とを解析する技術に関するものである。
【０００２】
【従来の技術】
自動車のエンジン等、可動部品を有する機械においては、いかにして振動・騒音を低減させるかが重要な課題である。そのため、振動・騒音をより高い精度で解析できるより実用的な技術の確立が要望されている。機械要素単体の振動伝達特性に基づく振動解析については、汎用の構造解析ソフトウェア等を利用することにより、ほぼ実用レベルで実施できる段階に達していると言える。そして、近年では、複数の機械要素間の相互作用を考慮することにより、より高い精度で振動解析を行い得る技術の確立に重点が移行しつつある。
【０００３】
複数の機械要素間の相互作用を取り扱うためには、可動部品である物体の剛体変位と弾性振動との連成を解析することが重要な課題の一つである。そして、その課題を解決する一理論として機構解析による理論が既に知られている。この機構解析理論は、その当初は、複数の剛体が関節において互いに結合されたものを適用対象とするとともに、関節による結合を拘束条件式で表した微分代数方程式を用いる理論として構築された。その後、この機構解析理論は、物体の弾性振動を弾性モードの重ね合わせにより表現することにより、その適用対象を剛体から弾性体へと拡張された。この機構解析理論においては、物体の位置はデカルト座標、方向はオイラーパラメータ、弾性振動はモード座標によりそれぞれ表されている。モード座標は、弾性振動をモード解析する際に使用される、互いに直交する複数のモードベクトルにより定義されるものである。
【０００４】
【発明が解決しようとする課題，課題解決手段および発明の効果】
しかし、この機構解析理論に基づく従来の物体変位・振動解析方法では、物体の剛体変位と弾性振動とを記述する運動方程式において種類の異なる複数の座標系が混在する。このことは、剛性運動と弾性振動との連成を理由に、上記運動方程式において連成項の数を増加させたり、各連成項の構成を複雑にする可能性が高いことを意味する。そのため、この従来の物体変位・振動解析方法では、物体の剛体変位と弾性振動とを解析するための運動方程式の取り扱いが複雑で、計算に多くの工程や時間が必要であった。
【０００５】
さらに、この従来の物体変位・振動解析方法は、多関節の機構運動において低次の弾性振動の効果を考慮する場合には適しているが、弾性自由度が大きい弾性振動において剛体変位の効果を考慮する場合には十分に適しているとは言えなかった。
【０００６】
このような事情を背景とし、本発明は、エンジンの構成部材である可動部品の剛体変位と弾性振動とを解析する技術であって、弾性自由度が大きい弾性振動において剛体変位の効果を考慮する場合にも適しているものを提供することを課題としてなされたものであり、本発明によって、下記各態様のエンジン部品変位・振動解析方法，エンジン部品変位・振動解析用記録媒体およびエンジン部品変位・振動解析装置が得られる。各態様は、請求項と同様に、項に区分し、各項に番号を付し、必要に応じて他の項の番号を引用する形式で記載する。これは、本明細書に記載の技術的特徴およびそれらの組合わせのいくつかの理解を容易にするためであり、本明細書に記載の技術的特徴やそれらの組合わせが以下の態様に限定されると解釈されるべきではない。
【０００７】
（１）エンジンの構成部材である可動部品の、全体座標系上における剛体としての変位である剛体変位と、弾性体としての振動である弾性振動とをコンピュータにより解析する方法であって、
前記可動部品を有限要素法により解析すべく、その可動部品を複数の仮想的な要素に分割した有限要素モデルを前記可動部品の諸元に基づき作成するステップと、
前記有限要素モデルによりモード解析を行うために、前記全体座標系と、その全体座標系上を運動しつつ前記有限要素モデルを局所的に観測するための局所観測座標系であって、その有限要素モデルに固定されることなく有限要素モデルとは独立しているが少なくともその有限要素モデルの近傍に位置して運動を行うものとを設定するステップと、
それら設定した全体座標系および局所観測座標系を用いて、前記有限要素モデルの剛体変位と弾性振動とを、剛体変位は剛体モードにより、弾性振動は弾性モードによりそれぞれ表すとともに、それら剛体モードと弾性モードとを区別しないで一律に取り扱うモード解析により解析するステップと
を含み、かつ、前記全体座標系と前記局所観測座標系とを設定するステップが、前記局所観測座標系において前記可動部品を局所的に観測することを周期的に繰り返す観測周期の各々のスタート時点を除く期間において、局所観測座標系に、その局所観測座標系の変位と姿勢と速度と角速度と加速度と角加速度との中の少なくとも一つが、前記スタート時点における前記可動部品のそれと実質的に等しい大きさから変化しない運動を行わせるとともに、各観測周期のエンド時点において前記可動部品の解析結果が得られるごとに、局所観測座標系を、次回の観測に備えて、その可動部品の剛体変位に関する最新の解析結果と実質的に一致するように更新するステップを含むことを特徴とするエンジン可動部品の変位・振動解析方法（請求項１）。
この方法においては、全体座標系とは別に可動部品を観測する座標系として、その可動部品の少なくとも近傍に位置して運動を行う局所観測座標系が用いられる。この局所観測座標系は、可動部品に固定されてその可動部品の剛体変位と終始同じ運動をさせられる従来の物体固定座標系とは異なり、可動部品の有限要素モデルに固定されることなく有限要素モデルとは独立しているが、少なくとも有限要素モデルの近傍に位置して運動を行うものとされ、有限要素モデルを局所的に観測する。この局所的な観測が周期的に繰り返され、局所観測座標系は、各観測周期のスタート時点を除く期間においては、スタート時点における有限要素モデルの変位と姿勢と速度と角速度と加速度と角加速度との中の少なくとも一つが実質的に等しい大きさから変化しない運動を行わせられるとともに、各観測周期のエンド時点において有限要素モデルの解析結果が得られるごとに、局所観測座標系が、次回の観測に備えて、有限要素モデルの剛体変位に関する最新の解析結果と実質的に一致するように更新される。
この方法によれば、可動部品の剛体変位をモード解析する際の前提、すなわち、剛体変位がモード座標系上において微小であるという前提が容易に成立する。しかも、局所観測座標系の、各観測周期のスタート時点を除く期間における運動を比較的簡単に定義し得る。
よって、この方法によれば、小さな剛体変位を伴う可動部品のみならず大きな剛体変位を伴う可動部品についても容易に振動解析を行い得る。
ところで、例えば、自動車のエンジンの振動騒音を解析するためには、エンジンの回転中にそれのシリンダブロックに対してクランク軸が剛体変位すること，それらシリンダブロックの軸受とクランク軸のピンとの間に油膜が形成され、かつ、その油膜が非線形のばねとして作用すること等を考慮することが大切である。そして、この例においては、クランク軸がシリンダブロックに対して剛体変位する大きさは、そのクランク軸の弾性振動に比較すればかなり大きなものとなる。そして、本項に記載の方法は、この例の振動解析にも適用可能である。
また、本項に記載の方法は、単独の可動部品の振動解析に適用することが可能であり、また、上記の例、すなわち、シリンダブロックは車体に対して変位し、クランク軸はそのシリンダブロックに対して変位するエンジンの振動解析のように、相互作用を有する複数の可動部品の振動解析に適用することも可能である。
また、本項に記載の方法において「可動部品の少なくとも近傍に位置する」という要件は、それが成立すれば、可動部品の剛体変位をモード解析する際の前提、すなわち、剛体変位がモード座標系上において微小であるという前提が成立することになる要件であり、例えば、可動部品と局所観測座標系との相対変位（例えば、後述の個別変位ベクトルｕ’i ）が微小範囲にあるという要件であると定義したり、可動部品と局所観測座標系とが実質的に一致するという要件であると定義することが可能である。
また、本項に記載の方法において「可動部品の剛体変位」は、その可動部品の弾性振動の影響を受ける場合と、受けない場合とがある。受ける場合には、弾性振動の解析結果を考慮して剛体変位を決定することが必要であるが、受けない場合には、弾性振動の解析結果を待つことなくそれとは独立に剛体変位を決定することが可能である。よって、可動部品の剛体変位がその可動部品の弾性振動の影響を受ける場合には、可動部品の弾性振動を考慮して局所観測座標系の運動を定義することが必要であるが、受けない場合には、可動部品の弾性振動とは独立に局所観測座標系の運動を定義することが可能である。
（２）前記モード解析を、モード座標を一般座標として用いることによって前記可動部品の剛体変位と弾性振動とを記述する運動方程式を解くことにより行う(1)項に記載のエンジン可動部品の変位・振動解析方法。
この方法においては、モード座標を一般座標として用いることによって可動部品を記述する運動方程式を解くことにより、可動部品の剛体変位と弾性振動とが解析される。一方、局所観測座標系は、各観測周期のスタート時点を除く期間においては可動部品の剛体変位とは独立に、時間の関数に従って運動を行わせられる。したがって、可動部品を記述する運動方程式において、局所観測座標系の運動に関連する変数をモード座標に対して定数として扱うことが可能となり、その運動方程式を用いて可動部品の剛体変位と弾性振動とを計算することが容易になる。
ここに「運動方程式」の一例はラグランジュの運動方程式である。また「モード座標」は一般に、互いに交差する複数のモードベクトルにより定義される座標であると定義される。
（３）前記モード解析を、互いに直交する複数のモードベクトルにより定義されるモード座標を用いて行う(1)項または(2)項に記載のエンジン可動部品の変位・振動解析方法。
この方法においては、互いに直交する複数のモードベクトルにより定義されるモード座標を用いてモード解析が行われる。したがって、この方法によれば、可動部品を記述する運動方程式において剛体モードと弾性モードとの間に生じる連成項の数を少なくしたり、連成項の構成を簡単にすることが容易になり、ひいては、運動方程式を用いた計算が容易になる。
ここに「モードベクトル」は、固有モードベクトルを選ぶことができるのはもちろんであるが、その他のモードベクトルを選ぶことも可能である。例えば、静的補正モードベクトルを選ぶことも可能であり、その静的補正モードベクトルには例えば、アタッチメントモードベクトル，拘束モードベクトル等がある。
（４）前記モード解析を行うステップが、
(i)そのモード解析のための固有値λと、前記剛体モードを含むｎ個の正規直交化された複数のモードベクトルφj とを取得するステップと、
(ii)前記局所観測座標系により定義されたベクトルを前記全体座標系により定義されるベクトルに変換するための方向余弦マトリクスＡを計算するステップと、
(iii)最初の観測周期のスタート時点の前記全体座標系上における前記局所観測座標系の原点の位置ベクトルｒ，速度ベクトルｖおよび加速度ベクトルａを設定するステップと、(iv)各観測周期ごとに、前記可動部品に設定された複数の節点ｉの各々に作用する外力Ｆi を取得するステップと、
(v)各観測周期ごとに、その観測周期のエンド時点における前記可動部品のモード座標系上におけるモード次数がｊであるモード変位ｑj ，モード速度ｑ. jおよびモード加速度ｑ: j を計算するステップと、
(vi)前記各観測周期ごとに、次回の観測に備えて、その次回の観測周期のスタート時点の前記全体座標系上における前記局所観測座標系の原点の位置ベクトルｒdt(0) ，速度ベクトルｖdt(0) および加速度ベクトルａdt(0) を計算するステップと、
(vii)前記各観測周期ごとに、次回の観測に備えて、前記可動部品の、モード次数がｊである剛体モードのモード変位ｑj ，モード速度ｑ. j およびモード加速度ｑ: j をいずれも０に初期化するステップと
を含む(1)項ないし(3)項のいずれかに記載のエンジン可動部品の変位・振動解析方法（請求項２）。
（５）当該エンジン可動部品の変位・振動解析方法は、力の相互作用を有する複数の可動部品の剛体変位と弾性振動とを解析するものであり、
前記(iv)のステップが、各可動部品の節点ｉに作用する外力Ｆi を、関連する他の可動部品の節点ｉ’との間における力の相互作用に基づいて計算するとともに、その力の相互作用を、前記全体座標系上における前記複数の可動部品の、それらの弾性振動を考慮しない場合の各位置ベクトルｒi ，ｒi'と、各速度ベクトルｒ. i ，ｒ. i'と、各加速度ベクトルｒ: i ，ｒ: i'との少なくとも一つに基づいて計算するステップを含む(4)項に記載のエンジン可動部品の変位・振動解析方法（請求項３）。
この方法によれば、力の相互作用を有するエンジンの複数の可動部品の剛体変位と弾性振動とを解析することができる。
（６）前記(ii)のステップにおいて、前記可動部品に予め設定された複数の節点のうち、前記局所座標系の原点に対応する節点を基点として、互いに線形独立な３つのベクトルｓ’1 ，ｓ’2 ，ｓ’3 を生成する３つの節点を選択し、その生成された３つのベクトルｓ’1 ，ｓ’2 ，ｓ’3 に基づいて前記方向余弦マトリクスＡを計算する(4)項または(5)項に記載のエンジン可動部品の変位・振動解析方法（請求項４）。
この方法によれば、方向余弦マトリクスＡを簡単に計算することができる。
（７）前記各観測周期ごとに、次回の観測に備えて、その次回の観測周期のスタート時点における前記方向余弦マトリクスＡdt(0) ，前記局所観測座標系の角速度ωの歪対称マトリクスω*dt (0) およびその局所観測座標系の角加速度αの歪対称マトリクスα*dt (0) を計算し、計算した方向余弦マトリクスＡdt(0) ，歪対称マトリクスω*dt (0) および歪対称マトリクスα*dt (0) を、前記ｎ個の正規直交化された複数のモードベクトルφj により定義されるモード座標系により前記モード解析を行うことに起因する誤差が、方向余弦マトリクスＡdt(0) ，歪対称マトリクスω*dt (0) および歪対称マトリクスα*dt (0) に生じないように補正する(4)項ないし(6)項のいずれかに記載のエンジン可動部品の変位・振動解析方法（請求項５）。
（８）(1)項ないし(7)項のいずれかに記載のエンジン可動部品の変位・振動解析方法あるいは物体変位・振動解析方法を実施するためにコンピュータにより実行されるプログラムをコンピュータ読み取り可能に記録した記録媒体。
ここにおける「記録媒体」には例えば、フレキシブルディスク，磁気テープ，磁気ディスク，磁気ドラム，磁気カード，光ディスク，光磁気ディスク，ＲＯＭ，ＣＤ−ＲＯＭ，ＩＣカード，穿孔テープ等がある。
（９）コンピュータにより、エンジンの構成部材である可動部品の、全体座標系上における剛体としての変位である剛体変位と、弾性体としての振動である弾性振動とを解析するためのプログラムを記録した媒体であって、
前記プログラムは、
前記可動部品を有限要素法により解析すべく、その可動部品を複数の仮想的な要素に分割した有限要素モデルを前記可動部品の諸元に基づき作成するステップと、
前記有限要素モデルによりモード解析を行うために、前記全体座標系と、その全体座標系上を運動しつつ前記有限要素モデルを局所的に観測するための局所観測座標系であって、その有限要素モデルに固定されることなく有限要素モデルとは独立しているが少なくともその有限要素モデルの近傍に位置して運動を行うものとを設定するステップと、
それら全体座標系と局所観測座標系とを用いて、前記有限要素モデルの剛体変位と弾性振動とを、剛体変位は剛体モードにより、弾性振動は弾性モードによりそれぞれ表すとともに、それら剛体モードと弾性モードとを区別しないで一律に取り扱うモード解析により解析するステップと
を含み、かつ、前記全体座標系と前記局所観測座標系とを設定するステップが、前記局所観測座標系において前記可動部品を局所的に観測することを周期的に繰り返す観測周期の各々のスタート時点を除く期間において、局所観測座標系に、各その局所観測座標系の変位と姿勢と速度と角速度と加速度と角加速度との中の少なくとも一つが、前記スタート時点における前記可動部品のそれと実質的に等しい大きさから変化しない運動を行わせるとともに、各観測周期のエンド時点において前記可動部品の解析結果が得られるごとに、局所観測座標系を、次回の観測に備えて、その可動部品の剛体変位に関する最新の解析結果と実質的に一致するように更新するステップを含むことを特徴とする記録媒体（請求項６）。
前記 (2) 項ないし (7) 項の各々に記載の特徴は、本項の記録媒体にも適用可能である。
（１０）エンジンの構成部材である可動部品の、剛体としての変位である剛体変位と、弾性体としての振動である弾性振動とを解析するための装置であって、
前記可動部品を有限要素法により解析すべく、その可動部品を複数の仮想的な要素に分割した有限要素モデルを前記可動部品の諸元に基づき作成する手段と、
前記有限要素モデルによりモード解析を行うために、前記全体座標系と、その全体座標系上を運動しつつ前記有限要素モデルを局所的に観測するための局所観測座標系であって、その有限要素モデルに固定されることなく有限要素モデルとは独立しているが少なくともその有限要素モデルの近傍に位置して運動を行うものとを設定する手段と、
それら設定した全体座標系および局所観測座標系を用いて、前記有限要素モデルの剛体変位と弾性振動とを、剛体変位は剛体モードにより、弾性振動は弾性モードによりそれぞれ表すとともに、それら剛体モードと弾性モードとを区別しないで一律に取り扱うモード解析により解析する手段と
を含み、かつ、前記全体座標系と前記局所観測座標系とを設定する手段が、前記局所観測座標系において前記可動部品を局所的に観測することを周期的に繰り返す観測周期の各々のスタート時点を除く期間において、局所観測座標系に、その局所観測座標系の変位と姿勢と速度と角速度と加速度と角加速度との中の少なくとも一つが、前記スタート時点における前記可動部品のそれと実質的に等しい大きさから変化しない運動を行わせるとともに、各観測周期のエンド時点において前記可動部品の解析結果が得られるごとに、局所観測座標系を、次回の観測に備えて、その可動部品の剛体変位に関する最新の解析結果と実質的に一致するように更新する手段を含むことを特徴とするエンジン可動部品の変位・振動解析装置（請求項７）。
前記 (2)項ないし(7)項の各々に記載の特徴は、本項のエンジン可動部品の剛体変位・振動解析装置にも適用可能である。
【０００８】
【発明の実施の形態】
以下、本発明のさらに具体的な一実施形態を図面に基づいて詳細に説明する。
【０００９】
本実施形態は、物体の剛体変位と弾性振動とを解析する物体変位・振動解析方法である。この方法は、剛体自由度は小さいが弾性自由度は大きい機械の振動騒音の解析に適している。
【００１０】
本実施形態は、種々の新たな理論に基づいて構成されており、以下、その新たな理論を説明する。
【００１１】
１．局所観測座標系の導入
【００１２】
図１に示すように、全体座標系｛Ｘ，Ｙ，Ｚ｝に対し、局所観測座標系｛ｘ，ｙ，ｚ｝なる新しい概念の座標系を導入する。この座標系は、図２に示すように、従来の、物体に固定された座標系とは異なり、物体の近傍に位置しながら、物体の運動を局所的に観測するために設けられた座標系である。局所観測座標系ＬＯＣは、図３に示すように、ある時刻に物体Ｏと同一の運動をし、その状態から、物体Ｏの運動とは独立に、時間の関数に従って運動をする。
【００１３】
ここで、局所観測座標系と物体との関係を、位置，速度および加速度に関してそれぞれ説明する。
【００１４】
まず、図１に示すように、全体座標系から見た局所観測座標系の原点の位置ベクトルをｒ、方向余弦マトリクスをＡ、局所観測座標系から見た物体の節点ｉの位置ベクトルをｒ’i とすると、全体座標系から見た物体の節点ｉの位置ベクトルｒi は、
【００１５】
【数１】
なる式で表される。
【００１７】
ここで補足説明すれば、上記方向余弦マトリクスＡは、座標系を局所観測座標系から全体座標系に変換する座標変換マトリクスの一種である。また、上記位置ベクトルｒ’i において「’」は、そのベクトルが局所観測座標系により定義されることを示す表記であり、この表記は他の変数を表記する際にも用いられる。また、観測すべき物体には、その物体を複数の要素に仮想的に分割するためにＮ個の節点が設定されており、上記節点ｉはそれらＮ個の節点のうち任意のものを示している。
【００１８】
物体の剛体変位および弾性振動を考慮しない場合の、局所観測座標系から見た節点ｉの位置ベクトルをｓ’i とし、それら剛体変位と弾性振動との和を個別変位ベクトルｕ’i とすると、
【００１９】
【数２】
なる式が成立する。ここで、Ｎ個の節点の全変位ベクトルｕ’を、
【００２１】
【数３】
なる式で定義する。ここで、全変位ベクトルｕ’が十分小さく、モード解析が前提とする微小変位の範囲内にあり、かつ、全変位ベクトルｕ’が、弾性モードのみならず剛体モードをも含むｎ個の正規直交化された固有モードベクトルφj （モード次数ｊ＝１，・・・，ｎ）の重ね合わせにより表すことができるとすると、全変位ベクトルｕ’が、モード変位ｑj を用いることにより、
【００２３】
【数４】
なる式で表される。ここで、モードマトリクスΦを、固有モードベクトルφj を用いることにより、
【００２５】
【数５】
なる式で定義すると、モードマトリクスΦは、固有モードの直交性を理由に、質量マトリクスＭと剛性マトリクスＫとの間において、
【００２７】
【数６】
なる式で表される関係を有する。ここで、「Ｉ」は単位マトリクスであり、「λ」は、固有値λj を対角成分として有するマトリクス、すなわち、
【００２９】
【数８】
である。ただし、モード次数ｊ＝１，・・・，ｎr が剛体モードを表すものであるとすると、固有値λi ＝０（ｊ＝１，・・・，ｎr ）となる。ここで、モード変位ベクトルｑを、マトリクスの転置を表す表記として「T 」を用いることにより、
【００３１】
【数９】
なる式で定義すると、前記式(4) は、前記式(5) を用いることにより、
【００３３】
【数１０】
に変形される。前記式(3) における全変位ベクトルｕ’と個別変位ベクトルｕ’i との関係と同様に、固有モードベクトルφj を、各節点ｉごとに分けてφijで表すことにより、
【００３５】
【数１１】
なる式で定義すると、前記式(2) における個別変位ベクトルｕ’i は、
【００３７】
【数１２】
なる式で表される。したがって、前記式(1) は、前記式(2) と上記式(12)とを用いることにより、
【００３９】
【数１３】
に変形される。この式(13)における「ｒi 」は、図１に示すものと同じである。また、この式(13)は、局所観測座標系と物体との間の、位置に関する関係式を意味する。
【００４１】
上記式(13)を時間で微分すると、局所観測座標系と物体との間の、速度に関する関係式として、
【００４２】
【数１４】
なる式が得られる。ここに、「Ａ. 」は、方向余弦マトリクスＡを時間で微分したものを意味し、「ｑ. j 」は、モード変位ｑj を時間で微分したものを意味する。
【００４４】
ここで、全体座標系に対する局所観測座標系の角速度ベクトルをω、その歪対称マトリクスをω* で表すと、
【００４５】
【数１５】
なる式で表される関係が成立する。この式(15)を上記式(14)に代入し、全体座標系に対する局所観測座標系の速度ベクトルを改めてｖで表すと、速度の関係式として、
【００４７】
【数１６】
なる式が得られる。
【００４９】
前記式(14)を時間で微分すると、局所観測座標系と物体との間の、加速度に関する関係式として、
【００５０】
【数１７】
なる式が得られる。ここに、「Ａ: 」は、前記Ａ. を時間で微分したものを意味し、「ｑ: j 」は、前記ｑ. j を時間で微分したものを意味する。一方、前記式(15)を時間で微分すると、
【００５２】
【数１８】
なる式が得られる。ここに、「ω.*」は、歪対称マトリクスω* を時間で微分したものを意味する。この式(18)に前記式(15)を代入するとともに、全体座標系に対する局所観測座標系の角加速度ベクトルをα、その歪対称マトリクスをα* としてそれぞれ表すと、
【００５４】
【数１９】
なる式が得られる。この式(19)と前記式(15)とを前記式(17)に代入するとともに、全体座標系に対する局所観測座標系の加速度ベクトルを改めてａで表すと、
【００５６】
【数２０】
なる式が得られる。
【００５８】
以上のようにして、局所観測座標系｛ｘ，ｙ，ｚ｝を導入することにより、前記式(13)，式(16)および式(20)が、弾性振動のみならず剛体変位をも伴う物体の節点ｉの位置，速度および加速度を、全体座標系と局所観測座標系とモード座標系とを用いて記述する式として得られた。
【００５９】
なお付言すれば、それら３式において全体座標系または局所観測座標系により定義される変数は、後述のように、物体の剛体変位とは独立に、時間の関数として定義される。すなわち、その変数は、時間依存変数と言うことができるのである。一方、本実施形態において採用された物体の運動方程式は、後述のように、モード座標を一般座標として用いる。したがって、その運動方程式において上記時間依存変数は、一般座標との関係においては定数として取り扱うことが可能となる。このような観点からすれば、前記３式は、局所観測座標系の導入により、弾性振動のみならず剛体変位をも伴う物体の節点ｉの位置，速度および加速度を、モード座標系のみを用いて上記運動方程式に組み込むことを可能にした式であると言うことができ、また、それら位置，速度および加速度を記述する座標系を、全体座標系および局所観測座標系からモード座標系に変換することを可能にした式であると言うこともできる。
【００６０】
２．モード座標系による運動方程式の導出
【００６１】
物体全体の運動エネルギーＴは、節点ｉの質量をｍi で表すと、
【００６２】
【数２１】
なる式で表される。また、弾性変形による歪みエネルギーＵは、前記式(3) の全変位ベクトルｕ’を用いることにより、
【００６４】
【数２２】
なる式で表される。この式(22)に前記式(10)を代入すると、
【００６６】
【数２３】
なる式が得られ、さらに、前記式(7) により表される関係を利用すると、
【００６８】
【数２４】
なる式が得られる。
【００７０】
一方、ラグランジュの運動方程式を、モード座標ｑk （モード次数ｋ＝１，・・・，ｎr ）を一般座標（generalized coordinates, general coordenates）として用いて表すと、
【００７１】
【数２５】
なる式が得られる。この式(25)は、時間に依存しないため、局所観測座標系を、物体の運動とは独立に、時間の関数に従って運動させることによる影響を何ら受けない。したがって、この式(25)は、局所観測座標系の運動を考慮することなく、解くことができるのである。
【００７３】
上記式(25)の右辺における「Ｑk 」は、モード座標系での一般化力（generalized force ）を表している。
【００７４】
左辺の第１項の括弧内の変数は、前記式(21)を用いることにより、
【００７５】
【数２６】
で表される。この式(26)は、前記式(16)を用いることにより、
【００７７】
【数２７】
なる式に変形され、転置して整理すると、
【００７９】
【数２８】
なる式が得られる。この式(28)を時間で微分し、前記式(15)を用いて整理すると、
【００８１】
【数２９】
なる式が得られる。
【００８３】
一方、前記式(25)の左辺の第２項は、前記式(21)を用いることにより、
【００８４】
【数３０】
で表される。この式(30)は、前記式(16)を用いることにより、
【００８６】
【数３１】
に変形され、転置して整理すると、
【００８８】
【数３２】
なる式が得られる。この式(32)の右辺の項は、前記式(29)の右辺の第２項と同一であるから、それら２式を引き算すれば、
【００９０】
【数３３】
なる式が得られる。この式(33)に前記式(20)を代入すると、
【００９２】
【数３４】
なる式が得られる。ここに、方向余弦マトリクスＡは正規直交マトリクスであるため、
【００９４】
【数３５】
なる式が成立する。また、前記式(6) により表される、固有モードの直交性という条件は、節点ｉの質量ｍi を用いて表すと、
【００９６】
【数３６】
となる。ただし、「δij」は、クロネッカーのデルタである。この式(36)と前記式(35)を用いると、前記式(34)は、
【００９８】
【数３７】
なる式に変形される。
【０１００】
一方、前記式(25)の左辺の第３項は、前記式(24)と前記式(8) とを用いることにより、
【０１０１】
【数３８】
なる式で表される。
【０１０３】
また、前記式(25)の右辺の一般化力Ｑk については、まず、全体座標系から見た節点ｉに作用する外力をＦi として、仮想仕事δＷを考えるととにも、
【０１０４】
【数３９】
なる式を用いる。ここに、「δｒi 」は、前記位置ベクトルｒi の仮想的な微小変化を表す。この式(39)は、前記式(13)を用いることにより、
【０１０６】
【数４０】
なる式に変形される。ここに、「δｑ j 」は、前記モード変位ｑ j の仮想的な微小変化を表す。ここで、微小変化δｑ j の係数をまとめることにより、
【０１０８】
【数４１】
なる式を定義すると、一般化力Ｑk は、
【０１１０】
【数４２】
なる式で表され、転置して整理すると、
【０１１２】
【数４３】
なる式が得られる。
【０１１４】
前記式(37)，式(38)および式(43)を前記式(25)、すなわち、ラグランジュの運動方程式に代入して整理すると、最終的に、モード座標系により記述された運動方程式として、
【０１１５】
【数４４】
なる式が得られる。
【０１１７】
この式(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)を用い、かつ、原点（ｉ＝０）であることを考慮することにより、
【０１２８】
【数４５】
なる式で表されるように更新される。ここに、「φ0j」は、物体の複数の節点のうち、局所観測座標系の原点に対応する節点の固有モードベクトルを示している。なお、「原点に対応する節点」の意味については後述する。
【０１３０】
(3) 局所観測座標系の原点の速度ベクトルｖdt(0) に関する更新
【０１３１】
時刻ｄｔにおける新しい局所観測座標系の原点の速度ベクトルｖdt(0) は、前記式(16)を用いることにより、
【０１３２】
【数４６】
なる式で表されるように更新される。
【０１３４】
(4) 局所観測座標系の原点の加速度ベクトルａdt(0) に関する更新
【０１３５】
時刻ｄｔにおける新しい局所観測座標系の原点の加速度ベクトルａdt(0) は、前記式(20)を用いることにより、
【０１３６】
【数４７】
なる式で表されるように更新される。
【０１３８】
(5) 方向余弦マトリクスＡdt(0) の更新
【０１３９】
物体の剛体変位を考慮せず、かつ、時刻ｄｔである場合の、全体座標系から見た節点ｉの位置ベクトルｒdt i は、前記式(13)を用いることにより、
【０１４０】
【数４８】
なる式で表される。また、時刻ｄｔに、前記位置ベクトルｓ’i を全体座標系から見た場合のベクトルｓdt i は、
【０１４２】
【数４９】
なる式で表され、この式(49)に前記式(45)と式(48)とを代入すると、
【０１４４】
【数５０】
なる式が得られる。また、上記ベクトルｓdt i は、時刻ｄｔにおける方向余弦マトリクスＡdt(0) を用いることにより、
【０１４６】
【数５１】
なる式で表される。この式(51)を方向余弦マトリクスＡdt(0) について解くために、図４に示すように、物体の複数の節点の中から、局所観測座標系の原点と共同して四面体を構成する３つの節点であって、かつ、原点を始点、各節点を終点とする３つのベクトルｓ’1 ，ｓ’2 およびｓ’3 が互いに線形独立であるものを選択する。この選択により、
【０１４８】
【数５２】
なる式が得られる。ここに、「ｓ’1 」，「ｓ’2 」および「ｓ’3 」はそれぞれ、同図において、原点を始点、節点１，２および３をそれぞれ終点とする、互いに線形独立な３つのベクトルを意味している。このように、ベクトルｓ’1 ，ｓ’2 およびｓ’3 は互いに線形独立であるため、それらベクトルｓ’1 ，ｓ’2 およびｓ’3 を成分とするマトリクスには逆マトリクスが存在する。よって、上記式(52)における方向余弦マトリクスＡdt(0) は、
【０１５０】
【数５３】
なる式で表されるように更新することができる。
【０１５２】
(6) 角速度の歪対称マトリクスω*dt (0) の更新
【０１５３】
時刻ｄｔにおける節点ｉの、物体の剛体変位に起因する速度成分ｒ.dt i は、前記式(16)を用いることにより、
【０１５４】
【数５４】
なる式で表される。ベクトルｓdt i を時間で微分して得られるベクトルｓ.dt i は、全体座標系に対して並進方向の速度成分を持たないため、
【０１５６】
【数５５】
なる式で表される。この式(55)に前記式(46)と上記式(54)とを代入すると、
【０１５８】
【数５６】
なる式が得られる。さらに、ベクトルｓ.dt i は、前記式(53)と時刻ｄｔにおける歪対称マトリクスω*dt (0) とを用いることにより、
【０１６０】
【数５７】
なる式で表される。したがって、前記式(53)と同様に、前述の、互いに線形独立なベクトルｓ’1 ，ｓ’2 およびｓ’3 を用いることにより、
【０１６２】
【数５８】
なる式が求められる。よって、角速度の歪対称マトリクスω*dt (0) は、
【０１６４】
【数５９】
なる式で表されるように更新される。
【０１６６】
(7) 角加速度の歪対称マトリクスα*dt (0) の更新
【０１６７】
時刻ｄｔにおける節点ｉの、物体の剛体変位に起因する加速度成分ｒ:dt i は、前記式(20)を用いることにより、
【０１６８】
【数６０】
なる式で表される。ここに、「ｓ:dt i 」は、全体座標系に対して並進方向の加速度成分を持たないため、
【０１７０】
【数６１】
なる式で表される。この式(61)に前記式(47)と上記式(60)とを代入することにより、
【０１７２】
【数６２】
なる式が得られる。さらに、ベクトルｓ:dt i は、前記式(53)と式(59)と時刻ｄｔにおける角加速度の歪対称マトリクスα*dt (0) とを用いることにより、
【０１７４】
【数６３】
なる式で表される。よって、前記式(53)および式(58)と同様に、前述の、互いに線形独立なベクトルｓ’1 ，ｓ’2 およびｓ’3 を用いることにより、
【０１７６】
【数６４】
なる式が求められる。よって、角加速度の歪対称マトリクスα*dt (0) は、
【０１７８】
【数６５】
なる式で表されるように更新される。
【０１８０】
(8) 方向余弦マトリクスＡdt (0)の補正
【０１８１】
本実施形態においては、物体の剛体変位を剛体モードの重ね合わせにより解析し、その際、複数の固有モードベクトルが互いに直交するモード座標系を用いて物体の剛体変位を記述するため、前記式(53)により求めた方向余弦マトリクスＡdt (0)には微小な誤差が存在する。以下、この誤差の存在を具体的に説明する。
【０１８２】
全体座標系から見た局所観測座標系の座標軸ベクトルをｘ< ，ｙ< ，ｚ< で表すこととすると、方向余弦マトリクスＡは、
【０１８３】
【数６６】
で表すことができる。各座標軸ベクトルは単位ベクトルであるため、
【０１８５】
【数６７】
なる式が成立する。さらに、３つの座標軸ベクトルは互いに直交するため、
【０１８７】
【数６８】
なる式が成立する。
【０１８９】
しかし、前記式(53)により更新された方向余弦マトリクスＡdt (0)の各列ベクトルは、小さい誤差を含んだｘ- ，ｙ- ，ｚ- であるため、上記式(67)により表される関係および上記式(68)により表される関係が厳密には成立しない。
【０１９０】
そこで、本実施形態においては、その誤差が解消されるように、方向余弦マトリクスＡdt (0)が適切に補正される。以下、その補正方法を説明する。
【０１９１】
まず、列ベクトルｘ- の方向は正しいと仮定として、列ベクトルｘ< を、
【０１９２】
【数６９】
なる式を用いて決定される。
【０１９４】
次に、求めるべき列ベクトルｙ< は、上記のようにして求められた列ベクトルｘ< と列ベクトルｙ- と同一平面上にあるとして、列ベクトルｚ< が、
【０１９５】
【数７０】
なる式を用いて決定される。
【０１９７】
最後に、列ベクトルｙ< が、
【０１９８】
【数７１】
なる式を用いて決定される。
【０２００】
(9) 角速度の歪対称マトリクスω*dt (0) の補正
【０２０１】
方向余弦マトリクスＡdt (0)の場合と同様に、更新された角速度ωの歪対称マトリクスω*dt (0) にも微小な誤差が含まれる。そこで、本実施形態においては、歪対称マトリクスω*dt (0) も適切に補正される。以下、その補正方法を説明する。
【０２０２】
前記式(59)で求めた歪対称マトリクスω*dt (0) は、一般にフルマトリクスとなるが、本来は、角速度ベクトルω、すなわち、
【０２０３】
【数７２】
の歪対称マトリクスであって、
【０２０５】
【数７３】
なる式で表される関係を満たすはずである。ここで、誤差を含んだ歪対称マトリクスω*dt (0) の各成分を、
【０２０７】
【数７４】
なる式で表す。歪対称マトリクスω*dt (0) の補正方法は種々考えられるが、本実施形態においては、次のように行われる。
【０２０９】
すなわち、まず、歪対称マトリクスω*dt (0) の対角成分が０でない場合にはすべて誤差であるとして、対角成分ω * ii ＝０（ｉ＝１，２，３）とされる。次に、基本成分ωx ，ωy ，ωz が、
【０２１０】
【数７５】
のように補正される。
【０２１２】
(10)角加速度の歪対称マトリクスα*dt (0) の補正
【０２１３】
歪対称マトリクスα*dt (0) の補正方法は、歪対称マトリクスω*dt (0) の場合と基本的に同じである。以下、具体的に説明する。
【０２１４】
前記式(65)により更新された歪対称マトリクスα*dt (0) は、一般にフルマトリクスとなるが、本来は、角加速度ベクトル、すなわち、
【０２１５】
【数７８】
の歪対称マトリクスであるため、
【０２１７】
【数７９】
なる式で表される関係を満たすはずである。ここで、誤差を含んだ歪対称マトリクスα*dt (0) の各成分を、
【０２１９】
【数８０】
なる式で表す。そして、歪対称マトリクスα*dt (0) の補正においては、まず、歪対称マトリクスα*dt (0) の対角成分が０でない場合にはすべて誤差であるとして、対角成分α * ii ＝０（ｉ＝１，２，３）とされる。次に、基本成分αx ，αy ，αz が、
【０２２１】
【数８１】
のように補正される。
【０２２３】
(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 のはりと、質量を無視できる直径Ｄの円板とで構成されている。７個の質点は、正八面体の各頂点および重心に配置されている。このジャイロスコープは、質点７の位置において、全体座標系の原点に球面ジョイントを介して連結されており、また、−Ｚ方向に重量が作用している。このタイプのジャイロスコープは一般に、ラグランジュのこまと呼ばれている。
【０２４５】
２．局所観測座標系の定義および更新
【０２４６】
局所観測座標系は、初期状態において、図１１に示す座標系｛ｘ，ｙ，ｚ｝のように配置されており、かつ、等加速度および等角加速度で運動するようにされている。したがって、各観測周期（ｔ＝０〜ｄｔ）において、加速度ベクトルａおよび角加速度ベクトルαは、
【０２４７】
【数８４】
なる２式でそれぞれ表される。したがって、速度ベクトルｖおよび位置ベクトルｒは、
【０２４９】
【数８６】
なる２式でそれぞれ表される。また、同様に、角速度ベクトルωは、
【０２５１】
【数８８】
なる式で表される。
【０２５３】
方向余弦マトリクスＡは次のようにして計算される。
【０２５４】
まず、時刻０からｄｔまでの局所観測座標系の回転軸を表す単位ベクトルをｋ、回転角をθで表す。次に、その回転角θが微小であるとして、
【０２５５】
【数８９】
なる式が成立すると仮定する。この式における「ｋθ」が表す回転と等価な回転マトリクスをＲk-THで表し、ベクトルｋを、
【０２５７】
【数９０】
なる式で表すと、回転マトリクスＲk-THは、
【０２５９】
【数９１】
なる式で表される。ただし、ｃTH＝cos θ、ｓTH＝sin θ、ｖTH＝１−cos θである。そして、方向余弦マトリクスＡが、
【０２６１】
【数９２】
なる式として計算できる。
【０２６３】
なお、局所観測座標系は、前述の更新方法に従い、かつ、後述の数値積分におけると同じ刻み時間ｄｔごとに更新される。
【０２６４】
３．剛体変位が正しく計算できることの検証
【０２６５】
まず、本実施形態により剛体変位が正しく計算できることを確認するため、剛体モードのみ考慮して計算を行った。剛体モードのみ考慮すると、ジャイロスコープは質点７の位置において球面ジョイントを介して連結されているため、この運動系の自由度は３となる。そして、実際には、剛体モードとして、ｘ，ｙ，ｚ各軸まわりの３つの回転モードを選び、前記式(6) により正規化した固有モードベクトルを採用した。図１２には計算に必要な各種変数・定数の諸元、図１３には初期条件がそれぞれ示されている。なお、前記式(44)の数値積分法として４次ルンゲ・クッタ法を用いた。
【０２６６】
図１４において実線のグラフは、本実施形態により計算されたジャイロスコープの質点１の軌跡をＸＹ平面に正射影することによって得られたものである。また、同図において網かけ線のグラフは、比較のために、オイラー角を一般座標として用いたラグランジュの運動方程式を用いて計算された質点１の軌跡を同様にして示すものである。いずれの軌跡も、ジャイロスコーいずれの軌跡も、ジャイロスコープが、初期角速度と重力の影響により、章動を伴いながら歳差運動をする様子を示している。両者の計算結果は一致している。一方、オイラー角を一般座標として用いたラグランジュの運動方程式を用いた計算結果は、実際値と精度よく一致することが既に証明されている。このことから明らかなように、本実施形態は剛体変位を正しく解析できる物体変位・振動解析方法である。
【０２６７】
４．剛体変位と弾性振動が正しく計算できることの検証
【０２６８】
次に、上記ジャイロスコープにおける長さｌ1 ，ｌ2 のはりを弾性はりとして計算を行った。図１５には、弾性はりの諸元が示されている。また、前述のＭＳＣ−ＮＡＳＴＲＡＮにより固有値解析を行い、その結果に基づき、固有振動数が３４．７２７Ｈｚである円板の並進運動を伴う２つの曲げモードを採用した。さらに、初期条件として、図１３に示す条件に加えて、２つの弾性モードにそれぞれモード変位が０．００５である初期変位を与えた。
【０２６９】
図１６における実線のグラフは、本実施形態により計算した、弾性体ジャイロスコープの質点１の軌跡をＸＹ平面に正射影したものである。この軌跡は、ジャイロスコープが弾性振動と章動を伴いながら歳差運動をする様子を示している。同図における網かけ線のグラフは、剛体モードのみを採用した場合の軌跡を示しており、図１４における実線のグラフと同じものである。図１６において実線と網かけ線とを互いに比較してみると、弾性はりの弾性効果により、章動の大きさが若干変化していることが分かる。このことから明らかなように、本実施形態は剛体変位のみならず弾性振動をも正確に計算できる物体変位・振動解析方法である。
【０２７０】
５．全エネルギーが保存されているか否かの観点からの検証
【０２７１】
次に、本実施形態の妥当性をエネルギーの観点から検証するため、図１７には、ジャイロスコープの全エネルギーが時刻によって変化する履歴が示されている。同図における実線のグラフは、弾性体ジャイロスコープの全エネルギー、網かけ線のグラフは剛体ジャイロスコープの全エネルギーをそれぞれ示しており、図１６における実線のグラフと網かけ線のグラフとにそれぞれ対応している。図１７における両結果の差は、弾性ジャイロスコープに与えた初期変位の歪みエネルギー分に相当する。ジャイロスコープが弾性体である場合にも剛体である場合にも、全エネルギーが保存されていることが分かる。このことからも明らかなように、本実施形態は剛体変位のみならず弾性振動をも正確に計算できる物体変位・振動解析方法である。
【０２７２】
６．むすび
【０２７３】
以上の説明から明らかなように、本実施形態は、剛体変位および弾性振動の計算結果の観点ならびに全エネルギーの計算結果の観点からして、妥当であることが検証された。
【０２７７】
以上、本発明の一実施形態を図面に基づいて詳細に説明したが、これは文字通り例示であり、本発明は、前記〔発明が解決しようとする課題，課題解決手段および発明の効果〕の項に記載された態様を始めとして、当業者の知識に基づいて種々の変形，改良を施した形態で実施することができる。
【図面の簡単な説明】
【図１】本発明の一実施形態である物体変位・振動解析方法において使用されている全体座標系と局所観測座標系とを示す斜視図である。
【図２】従来の物体固定座標系の運動を説明するための正面図である。
【図３】図１における局所観測座標系の運動を説明するための正面図である。
【図４】図１に示す局所観測座標系が上記物体変位・振動解析方法により更新される様子を説明するための斜視図である。
【図５】上記物体変位・振動解析方法を実施するのに好適な物体変位・振動解析装置を示すブロック図である。
【図６】図５のＣＤ−ＲＯＭに記憶されている振動解析プログラムの前半部分を示すフローチャートである。
【図７】上記振動解析プログラムの後半部分を示すフローチャートである。
【図８】図６および図７の振動解析プログラムの実行内容を経時的に説明するための正面図である。
【図９】上記振動解析プログラムの実行内容を経時的に説明するための別の正面図である。
【図１０】上記振動解析プログラムの実行内容を経時的に説明するためのさらに別の正面図である。
【図１１】上記物体変位・振動解析方法の妥当性を検証するために用いた解析モデルを示す斜視図である。
【図１２】上記解析モデルの諸元を表形式で示す図である。
【図１３】上記解析モデルの初期条件を表形式で示す図である。
【図１４】上記解析モデルを剛体と仮定した場合の解析結果を上記物体変位・振動解析方法を用いた場合と、オイラー角を一般座標とする運動方程式を用いた場合とについて対比して示すグラフである。
【図１５】上記解析モデルを弾性体と仮定した場合の諸元を表形式で示す図である。
【図１６】上記物体変位・振動解析方法を用いた解析結果を、上記解析モデルを弾性体と仮定した場合と剛体と仮定した場合とについて対比して示すグラフである。
【図１７】上記物体変位・振動解析方法の妥当性を全エネルギーが保存されているか否かの観点から検証するためのグラフである。
【符号の説明】
１２ コンピュータ
３０ ＣＤ−ＲＯＭ
ＬＯＣ 局所観測座標系
Ｏ 物体
Claims (7)
- エンジンの構成部材である可動部品の、全体座標系上における剛体としての変位である剛体変位と、弾性体としての振動である弾性振動とをコンピュータにより解析する方法であって、
前記可動部品を有限要素法により解析すべく、その可動部品を複数の仮想的な要素に分割した有限要素モデルを前記可動部品の諸元に基づき作成するステップと、
前記有限要素モデルによりモード解析を行うために、前記全体座標系と、その全体座標系上を運動しつつ前記有限要素モデルを局所的に観測するための局所観測座標系であって、その有限要素モデルに固定されることなく有限要素モデルとは独立しているが少なくともその有限要素モデルの近傍に位置して運動を行うものとを設定するステップと、
それら設定した全体座標系および局所観測座標系を用いて、前記有限要素モデルの剛体変位と弾性振動とを、剛体変位は剛体モードにより、弾性振動は弾性モードによりそれぞれ表すとともに、それら剛体モードと弾性モードとを区別しないで一律に取り扱うモード解析により解析するステップと
を含み、かつ、前記全体座標系と前記局所観測座標系とを設定するステップが、前記局所観測座標系において前記可動部品を局所的に観測することを周期的に繰り返す観測周期の各々のスタート時点を除く期間において、局所観測座標系に、その局所観測座標系の変位と姿勢と速度と角速度と加速度と角加速度との中の少なくとも一つが、前記スタート時点における前記可動部品のそれと実質的に等しい大きさから変化しない運動を行わせるとともに、各観測周期のエンド時点において前記可動部品の解析結果が得られるごとに、局所観測座標系を、次回の観測に備えて、その可動部品の剛体変位に関する最新の解析結果と実質的に一致するように更新するステップを含むことを特徴とするエンジン可動部品の変位・振動解析方法。 - 前記モード解析を行うステップが、
(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) に生じないように補正する請求項２ないし４のいずれかに記載のエンジン可動部品の変位・振動解析方法。
- コンピュータにより、エンジンの構成部材である可動部品の、全体座標系上における剛体としての変位である剛体変位と、弾性体としての振動である弾性振動とを解析するためのプログラムを記録した媒体であって、
前記プログラムは、
前記可動部品を有限要素法により解析すべく、その可動部品を複数の仮想的な要素に分割した有限要素モデルを前記可動部品の諸元に基づき作成するステップと、
前記有限要素モデルによりモード解析を行うために、前記全体座標系と、その全体座標系上を運動しつつ前記有限要素モデルを局所的に観測するための局所観測座標系であって、その有限要素モデルに固定されることなく有限要素モデルとは独立しているが少なくともその有限要素モデルの近傍に位置して運動を行うものとを設定するステップと、
それら全体座標系と局所観測座標系とを用いて、前記有限要素モデルの剛体変位と弾性振動とを、剛体変位は剛体モードにより、弾性振動は弾性モードによりそれぞれ表すとともに、それら剛体モードと弾性モードとを区別しないで一律に取り扱うモード解析により解析するステップと
を含み、かつ、前記全体座標系と前記局所観測座標系とを設定するステップが、前記局所観測座標系において前記可動部品を局所的に観測することを周期的に繰り返す観測周期の各々のスタート時点を除く期間において、局所観測座標系に、各その局所観測座標系の変位と姿勢と速度と角速度と加速度と角加速度との中の少なくとも一つが、前記スタート時点における前記可動部品のそれと実質的に等しい大きさから変化しない運動を行わせるとともに、各観測周期のエンド時点において前記可動部品の解析結果が得られるごとに、局所観測座標系を、次回の観測に備えて、その可動部品の剛体変位に関する最新の解析結果と実質的に一致するように更新するステップを含むことを特徴とする記録媒体。 - エンジンの構成部材である可動部品の、剛体としての変位である剛体変位と、弾性体としての振動である弾性振動とを解析するための装置であって、
前記可動部品を有限要素法により解析すべく、その可動部品を複数の仮想的な要素に分割した有限要素モデルを前記可動部品の諸元に基づき作成する手段と、
前記有限要素モデルによりモード解析を行うために、前記全体座標系と、その全体座標系上を運動しつつ前記有限要素モデルを局所的に観測するための局所観測座標系であって、その有限要素モデルに固定されることなく有限要素モデルとは独立しているが少なくともその有限要素モデルの近傍に位置して運動を行うものとを設定する手段と、
それら設定した全体座標系および局所観測座標系を用いて、前記有限要素モデルの剛体変位と弾性振動とを、剛体変位は剛体モードにより、弾性振動は弾性モードによりそれぞれ表すとともに、それら剛体モードと弾性モードとを区別しないで一律に取り扱うモード解析により解析する手段と
を含み、かつ、前記全体座標系と前記局所観測座標系とを設定する手段が、前記局所観測座標系において前記可動部品を局所的に観測することを周期的に繰り返す観測周期の各々 のスタート時点を除く期間において、局所観測座標系に、その局所観測座標系の変位と姿勢と速度と角速度と加速度と角加速度との中の少なくとも一つが、前記スタート時点における前記可動部品のそれと実質的に等しい大きさから変化しない運動を行わせるとともに、各観測周期のエンド時点において前記可動部品の解析結果が得られるごとに、局所観測座標系を、次回の観測に備えて、その可動部品の剛体変位に関する最新の解析結果と実質的に一致するように更新する手段を含むことを特徴とするエンジン可動部品の変位・振動解析装置。
Priority Applications (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP11267499A JP3799437B2 (ja) | 1999-04-20 | 1999-04-20 | 物体変位・振動解析方法，装置および記録媒体 |
| US09/552,728 US6684168B1 (en) | 1999-04-20 | 2000-04-19 | Body displacement and vibration analysis method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP11267499A JP3799437B2 (ja) | 1999-04-20 | 1999-04-20 | 物体変位・振動解析方法，装置および記録媒体 |
Related Child Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2004280712A Division JP2005091370A (ja) | 2004-09-27 | 2004-09-27 | 物体変位・振動解析方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JP2000305922A JP2000305922A (ja) | 2000-11-02 |
| JP3799437B2 true JP3799437B2 (ja) | 2006-07-19 |
Family
ID=14592648
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP11267499A Expired - Lifetime JP3799437B2 (ja) | 1999-04-20 | 1999-04-20 | 物体変位・振動解析方法，装置および記録媒体 |
Country Status (2)
| Country | Link |
|---|---|
| US (1) | US6684168B1 (ja) |
| JP (1) | JP3799437B2 (ja) |
Families Citing this family (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP4599723B2 (ja) * | 2001-02-13 | 2010-12-15 | 株式会社豊田中央研究所 | 数値解析方法、それを実施するためのプログラムおよびそれを記録した記録媒体 |
| US7630869B2 (en) * | 2003-05-27 | 2009-12-08 | University Of Washington | Method for predicting vibrational characteristics of rotating structures |
| US20050250086A1 (en) * | 2004-05-10 | 2005-11-10 | Johnson Lyle A | Method of modeling gravimetric flow and plotting lagrange points |
| US8403479B2 (en) * | 2009-12-17 | 2013-03-26 | Johnson & Johnson Vision Care, Inc. | Contact lens eye model |
| CN103530272B (zh) * | 2013-09-26 | 2019-02-05 | 广东工业大学 | 一种用于界定机构运动高速区域的判别方法 |
| JP6127945B2 (ja) * | 2013-12-04 | 2017-05-17 | マツダ株式会社 | エンジンのコネクティングロッドの設計支援方法及びその設計支援装置 |
| US10386339B2 (en) | 2017-08-04 | 2019-08-20 | Crystal Instruments Corporation | Modal vibration analysis system |
| CN111680368B (zh) * | 2019-02-25 | 2023-07-25 | 中国石油天然气集团有限公司 | 底部框架支撑式塔器结构的获取方法及装置 |
| CN110362855B (zh) * | 2019-05-30 | 2023-04-25 | 中国海洋大学 | 一种基于交叉模态能敏感度的模态优选及模型修正方法 |
| WO2023276803A1 (ja) * | 2021-07-01 | 2023-01-05 | 京セラ株式会社 | ジャイロ式振動計 |
Family Cites Families (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US3597598A (en) * | 1964-12-14 | 1971-08-03 | North American Rockwell | Method and means for transforming the acceleration signals generated by accelerometers in a first coordinate system into acceleration signals in a second coordinate system |
| US4592034A (en) * | 1982-11-15 | 1986-05-27 | Cornell Research Foundation, Inc. | Acoustic emission source location on plate-like structures using a small array of transducers |
| US4870588A (en) * | 1985-10-21 | 1989-09-26 | Sundstrand Data Control, Inc. | Signal processor for inertial measurement using coriolis force sensing accelerometer arrangements |
| JP2730195B2 (ja) * | 1989-06-30 | 1998-03-25 | 三菱電機株式会社 | 結合振動特性解析装置 |
| JP2536958B2 (ja) * | 1990-09-04 | 1996-09-25 | 日産自動車株式会社 | 車両の振動制御装置 |
| US5533399A (en) * | 1992-09-30 | 1996-07-09 | Wayne State University | Method and apparatus for non-destructive measurement of elastic properties of structural materials |
| US5342017A (en) * | 1992-09-30 | 1994-08-30 | Cincinnati Milacron Inc. | Damped cantilevered support arm assembly |
| JP3320138B2 (ja) * | 1993-05-07 | 2002-09-03 | キヤノン株式会社 | 座標入力装置及び方法 |
| US5625575A (en) * | 1993-08-03 | 1997-04-29 | Lucent Technologies Inc. | Apparatus for modelling interaction of rigid bodies |
| US5400033A (en) * | 1994-02-07 | 1995-03-21 | Rockwell International Corporation | Tracking system for tracking targets with a spacecraft |
| JPH0882571A (ja) * | 1994-09-13 | 1996-03-26 | Mazda Motor Corp | 振動測定システム |
| US5915278A (en) * | 1995-02-27 | 1999-06-22 | Mallick; Brian C. | System for the measurement of rotation and translation for modal analysis |
| US5831164A (en) * | 1997-01-21 | 1998-11-03 | Conrad Technologies, Inc. | Linear and rotational accelerometer |
-
1999
- 1999-04-20 JP JP11267499A patent/JP3799437B2/ja not_active Expired - Lifetime
-
2000
- 2000-04-19 US US09/552,728 patent/US6684168B1/en not_active Expired - Fee Related
Also Published As
| Publication number | Publication date |
|---|---|
| US6684168B1 (en) | 2004-01-27 |
| JP2000305922A (ja) | 2000-11-02 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Felis | RBDL: an efficient rigid-body dynamics library using recursive algorithms | |
| Brüls et al. | Lie group generalized-α time integration of constrained flexible multibody systems | |
| JP2730195B2 (ja) | 結合振動特性解析装置 | |
| Wallrapp et al. | Simulation of deployment of a flexible solar array | |
| JP3799437B2 (ja) | 物体変位・振動解析方法，装置および記録媒体 | |
| Pappalardo et al. | Stability analysis of rigid multibody mechanical systems with holonomic and nonholonomic constraints | |
| CN111368466A (zh) | 一种基于频响函数参数修正的机械振动预测方法 | |
| RU2308764C2 (ru) | Перемещение виртуального сочлененного объекта в виртуальном пространстве с предотвращением столкновений сочлененного объекта с элементами окружающего пространства | |
| Psiaki et al. | Numerical solution of a generalized wahba problem for a spinning spacecraft | |
| Fisette et al. | Fully symbolic generation of complex multibody models | |
| Ashrafiuon et al. | Analysis and optimal design of spatial mechanical systems | |
| Buffinton | Kane’s Method in Robotics | |
| Wasfy | Modeling spatial rigid multibody systems using an explicit-time integration finite element solver and a penalty formulation | |
| Brank et al. | Dynamics and time-stepping schemes for elastic shells undergoing finite rotations | |
| JP4612989B2 (ja) | メカトロニクスシステムのシミュレーション方法 | |
| JP3383282B2 (ja) | ６自由度シミュレーション方法 | |
| JP2018015436A (ja) | 打具の挙動の解析装置 | |
| Hsu et al. | Recursive sensitivity analysis for constrained multi-rigid-body dynamic systems design optimization | |
| JP2005091370A (ja) | 物体変位・振動解析方法 | |
| Mishra et al. | Dynamic modelling and control of flexible link manipulators: Methods and scope-Part-1 | |
| Lugrís et al. | Efficient calculation of the inertia terms in floating frame of reference formulations for flexible multibody dynamics | |
| Dwarshuis et al. | A multinode superelement in the generalized strain formulation | |
| Zou et al. | An algorithm for rigid-body angular velocity and attitude estimation based on isotropic accelerometer strapdowns | |
| JP7202567B2 (ja) | 多体系の最小力学パラメータ同定装置、方法及びプログラム | |
| Cuadrado et al. | A hybrid global-topological real-time formulation for multibody systems |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20040727 |
|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20040927 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20050809 |
|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20051011 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20060322 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20060404 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20090512 Year of fee payment: 3 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100512 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110512 Year of fee payment: 5 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110512 Year of fee payment: 5 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120512 Year of fee payment: 6 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120512 Year of fee payment: 6 |
|
| S531 | Written request for registration of change of domicile |
Free format text: JAPANESE INTERMEDIATE CODE: R313532 |
|
| R350 | Written notification of registration of transfer |
Free format text: JAPANESE INTERMEDIATE CODE: R350 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20130512 Year of fee payment: 7 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20130512 Year of fee payment: 7 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20140512 Year of fee payment: 8 |
|
| EXPY | Cancellation because of completion of term |