節点を境界とする有限数の要素に物体を分割して行われる有限要素法による前記物体の変位及び変形の解析の結果から同物体の剛体運動及び弾性変形を解析する方法であって、固定された全体座標系における各節点の初期の座標値及び各要素の質量行列を取得する第1のステップと、前記物体の初期の重心位置及び慣性主軸を原点及び軸ベクトルとする局所座標系を設定するとともに、その局所座標系における前記各節点の初期の座標値を計算する第2のステップと、前記質量行列及び前記局所座標系における前記各節点の初期の座標値から、前記各節点の剛体モードの固有ベクトルを計算する第3のステップと、を行い、剛体運動及び弾性変形による前記各節点の変位を同各節点の基準変位として設定した上で、前記固有ベクトル、前記質量行列及び前記基準変位から、前記局所座標系における前記物体の剛体モード変位を計算する第4のステップと、前記剛体モード変位及び前記固有ベクトルから、前記局所座標系における前記剛体運動及び弾性変形後の前記物体の前記重心位置の座標値及び慣性主軸の回転行列をそれぞれ計算する第5のステップと、前記剛体運動及び弾性変形後の前記物体の重心位置及び慣性主軸を原点及び軸ベクトルとする新たな局所座標系を設定するとともに、前記局所座標系における前記各節点の初期の座標値と前記新たな局所座標系における前記各節点の剛体運動及び弾性変形後の座標値との偏差を前記各節点の弾性変形による変位として計算する第6のステップと、を行い、以後、前記全体座標系における前記局所座標系と前記新たな局所座標系との偏差が既定値以下となるまで、前記新たな局所座標系を前記局所座標系として、前記各節点の弾性変形による変位を同各節点の前記基準変位としてそれぞれ設定し直した上で、前記第4~第6のステップを繰り返し行い、前記偏差が前記既定値以下となった時点の前記全体座標系における前記新たな局所座標系の原点の座標値及び軸ベクトルの回転行列から前記物体の剛体運動を、同時点の前記各節点の弾性変形による変位から前記物体の弾性変形をそれぞれ求めるようにしたことを特徴とする剛体運動及び弾性変形の解析方法。こうした従来の解析方法では、解析対象となる物体の剛体運動による節点変位以下、剛体変位と記載すると弾性変形による節点変位とを分離して、それぞれ個別に表示することが可能である。しかしながら、こうした従来の解析方法では、解析対象となる物体は、既知の弾性モードで表現された物体である必要があり、剛性行列が時々刻々と変化することから弾性モードが未知の有限要素で表現された物体の非線形過渡応答計算結果の解析には、適用できないものとなっている。本発明は、こうした実情に鑑みてなされたものであり、その解決しようとする課題は、剛性行列が時々刻々と変化することから弾性モードが未知の有限要素で表現された物体の場合にも、その剛体変位と弾性変形とを分離して解析することのできる剛体運動及び弾性変形の解析方法を提供することにある。上記課題を解決するため、剛体運動及び弾性変形の解析方法としての請求項1に記載の発明は、節点を境界とする有限数の要素に物体を分割して行われる有限要素法による前記物体の変位及び変形の解析の結果から同物体の剛体運動及び弾性変形を解析する方法であって、1固定された全体座標系における各節点の初期の座標値及び各要素の質量行列を取得する第1のステップ、2前記物体の初期の重心位置及び慣性主軸を原点及び軸ベクトルとする局所座標系を設定するとともに、その局所座標系における前記各節点の初期の座標値を計算する第2のステップ、3前記質量行列及び前記局所座標系における前記各節点の初期の座標値から、前記各節点の剛体モードの固有ベクトルを計算する第3のステップ、を行い、剛体運動及び弾性変形による前記各節点の変位を同各節点の基準変位として設定した上で、4前記固有ベクトル、前記質量行列及び前記基準変位から、前記局所座標系における前記物体の剛体モード変位を計算する第4のステップ、5前記剛体モード変位及び前記固有ベクトルから、前記局所座標系における前記剛体運動及び弾性変形後の前記物体の前記重心位置の座標値及び慣性主軸の回転行列をそれぞれ計算する第5のステップ、6前記剛体運動及び弾性変形後の前記物体の重心位置及び慣性主軸を原点及び軸ベクトルとする新たな局所座標系を設定するとともに、前記局所座標系における前記各節点の初期の座標値と前記新たな局所座標系における前記各節点の剛体運動及び弾性変形後の座標値との偏差を前記各節点の弾性変形による変位として計算する第6のステップ、を行い、以後、前記全体座標系における前記局所座標系と前記新たな局所座標系との偏差が既定値以下となるまで、前記新たな局所座標系を前記局所座標系として、前記各節点の弾性変形による変位を同各節点の前記基準変位としてそれぞれ設定し直した上で、前記第4~第6のステップを繰り返し行い、前記偏差が前記既定値以下となった時点の前記全体座標系における前記新たな局所座標系の原点の座標値及び軸ベクトルの回転行列から前記物体の剛体運動を、同時点の前記各節点の弾性変形による変位から前記物体の弾性変形をそれぞれ求めるようにしている。こうした本解析方法では、非線形FEMの解析結果が得られると、まず第1のステップとして、固定された全体座標系における各節点の初期の座標値及び各要素の質量行列が取得される。そして第2のステップとして、物体の初期の重心位置及び慣性主軸を原点及び軸ベクトルとする局所座標系が設定され、その局所座標系における各節点の初期の座標値が計算される。続いて第3のステップとして、質量行列及び局所座標系における各節点の初期の座標値から、各節点の剛体モードの固有ベクトルが計算される。次に、剛体運動及び弾性変形による各節点の変位が非線形有限要素法の解析結果から取得され、その変位が各節点の基準変位として設定される。そして下記の第4~第6のステップにより構成された収束ループの処理に移行される。収束ループの処理が開始されると、まず第4のステップとして、固有ベクトル、質量行列及び基準変位から、局所座標系における物体の剛体モード変位が計算される。また第5のステップとして、剛体モード変位及び固有ベクトルから、局所座標系における剛体運動及び弾性変形後の物体の重心位置の座標値及び慣性主軸の回転行列がそれぞれ計算される。そして第6のステップとして、剛体運動及び弾性変形後の物体の重心位置及び慣性主軸を原点及び軸ベクトルとする新たな局所座標系が設定され、第2のステップで設定された局所座標系における各節点の初期の座標値とその新たな局所座標系における各節点の剛体運動及び弾性変形後の座標値との偏差が各節点の弾性変形による変位として計算される。この時点で物体の剛体モード変位を微小≒0とすることができれば、物体の剛体運動は、その重心位置の並進及び回転として上記新たな局所座標系からは除かれており、第6のステップで計算された各節点の弾性変形による変位には、剛体運動による変位は一切含まれていないことになる。しかしながら、物体の局所座標の位置及び回転行列が既知でない以上、剛体モード変位は、微小≒0である保証はなく、この時点で計算された各節点の弾性変形による変位には、剛体運動分が未だ若干残っている。そこで、本解析方法では、新たな局所座標系を局所座標系に、先に計算した各節点の弾性変形による変位を各節点の基準変位にそれぞれ設定し直した上で、上記収束ループの処理を再実行している。収束ループの処理が再実行されると、その処理において更新された各節点の基準変位からは、物体の剛体運動分が更に取り除かれるようになる。こうした収束ループの処理の再実行は、全体座標系における局所座標系と新たな局所座標系との偏差ずれが既定値以下となるまで繰り返される。したがって、収束ループを抜ける直前に計算された各節点の弾性変形による変位は、物体の剛体運動分が十分に取り除かれたものとなる。すなわち、このときの各節点の弾性変形による変位の計算値は、各節点の弾性変形による変位をほぼ正確に示すものとなる。またこのときの全体座標系における新たな局所座標系の原点の座標値及びその軸ベクトルの回転行列は、物体の剛体運動の並進分及び回転分をそれぞれ示すものとなる。したがって、上記解析方法によれば、剛性行列が時々刻々と変化することから弾性モードが未知の有限要素で表現された物体の場合にも、その剛体変位と弾性変形とを分離して解析することができる。以下、本発明の剛体変位及び弾性変形の解析方法を具体化した一実施の形態を、図1~図7を参照して詳細に説明する。なお、本実施の形態の解析方法は、例えばコネクティングロッドのような可動部品の動作中の剛体変位及び弾性変形の解析に用いることができる。本実施形態の解析方法は、非線形FEMによる物体の変位及び変形の解析結果からその物体の剛体変位と弾性変形とを分離して提示するものとなっている。非線形FEMによる物体の変形の解析は、既存の非線形FEMソフトウェアを用いて行われる。非線形FEMによる物体の変形の解析では、まず解析対象となる物体を、節点を境界とする有限数の要素に分割し、各節点の初期の座標値や拘束条件、弾性率やポアソン比といった各要素の物理量、各節点に加わる荷重などを入力する。非線形FEMソフトウェアは、この入力されたデータをもとに、各々の要素の挙動をモデルで近似して形成された連立方程式を解くことで、各節点の座標値や各要素に作用する応力などを、時間ステップ毎に求めている。なお、非線形FEMソフトウェアの解析結果は、ポストプロセッサーによって処理されることで、変形図やミーゼス応力図などの視覚化されたデータとして提示される。本実施の形態の解析方法では、こうした非線形FEMの解析結果から解析対象となる物体の剛体変位と弾性変形とが分離して求められる。次に、簡単なモデルを用いて、その基本的な原理について説明する。ここでは、図1に示すような、4つの節点を有する単一の要素からなる物体Mの剛体変位と弾性変形とを求める場合を考える。同図には、全体座標系における物体Mの初期の位置及び形状と、同じく全体座標系における物体Mの剛体運動及び弾性変形後の位置及び形状とが示されている。なお、全体座標系における各節点の初期の座標値は既知であり、全体座標系における各節点の剛体運動及び弾性変形後の座標値は、非線形FEMの解析結果から求められている。ここで、節点の一つ節点Pについて考える。ここでは、初期の節点Pの位置を、全体座標系の原点を始点とし、節点Pの初期の位置を終点とするベクトルp0で表す。また剛体運動及び弾性変形による節点Pの変位を、節点Pの初期の位置を始点とし、節点Pの剛体運動及び弾性変形後の位置を終点とするベクトルu0で表す。このときの変形後の節点Pの位置は、ベクトルp0とベクトルu0の和として表される。また、初期における物体Mの重心位置を基準とした節点Pの位置を、物体Mの初期の重心位置Gの位置を始点とし、節点Pの初期の位置を終点とするベクトルpで、剛体運動及び弾性変形による物体Mの重心位置Gの変位の並進分及び回転分をベクトルr及び回転行列Aで、節点Pの弾性変形による変位をベクトルuでそれぞれ表す。このときの全体座標系における節点Pの剛体運動及び弾性変形後の位置、すなわちベクトルp0とベクトルu0の和は、下式のように表される。そして、節点Pの剛体運動はrApで、その弾性変形はAuで表されることになる。p0u0=rApu=rApAuそのため、全体座標系における各節点の剛体運動及び弾性変形後の位置{p0}{u0}は、下式のように表される。なお、下式において、列ベクトル{p0}は、各節点の初期の位置を、列ベクトル{u0}は、各節点の変位を、列ベクトル{r}及び回転行列のマトリックス[A]は、剛体運動及び弾性変形による物体Mの重心位置Gの変位の並進分及び回転分を、列ベクトル{u}は、各節点の弾性変形による変位をそれぞれ表している。{p0}{u0}={r}[A]{p}[A]{u}上式において、{r}[A]{p}は各節点の剛体運動を、[A]{u}は各節点の弾性変形をそれぞれ示している。したがって、初期、並びに剛体運動及び弾性変形後の各節点の位置、物体Mの重心位置Gの変位の並進分及び回転分、及び初期の物体Mの重心位置Gを基準とした各節点の位置さえ解れば、弾性モードが未知であり、各節点の弾性変形による変位が不明であっても、変形後の各節点の剛体運動及び弾性変形を求めることができる。すなわち、図2に示すように、物体Mの重心位置Gを原点とし、その慣性主軸を軸ベクトルとする局所座標系を設定すると、物体Mの剛体運動は、全体座標系上での局所座標系の並進並びに回転として求めることができる。また、物体Mの弾性変形は、局所座標系上での各節点の変位として求めることができる。次に、本実施形態での剛体変位及び弾性変形の解析の具体的な手順について説明する。この解析は、図3の剛体運動及び弾性変形解析ルーチンに示される手順で行われる。非線形FEMによる解析結果が得られると、まずステップS100において、非線形FEMの解析結果から、固定された全体座標系における各節点の初期の座標値{p0}が取得される。またステップS101において、非線形FEMの解析結果から、物体Mの質量行列[M]が取得される。なお、質量行列[M]は、物体Mを構成する各要素の質量を行列式で表したものである。次にステップS102において、物体Mの形状と質量行列[M]とから、全体座標系における物体Mの重心位置の初期の座標値{g}、主慣性モーメント及び慣性主軸の回転行列[I]がそれぞれ計算される。その後、ステップS103において、全体座標系における各節点の初期の座標値{p0}から、物体Mの初期の重心位置を基準とした各節点の初期の座標値{p}が計算される。この座標値{p}は、全体座標系における物体Mの初期の重心位置及び慣性主軸を原点及び軸ベクトルとする局所座標系における各節点の座標値を指している。次に、ステップS104において、上記各節点の初期の座標値{p}と質量行列[M]とから、各節点の剛体モードの固有ベクトル[Φr]が計算される。以上の処理が終了すると、非線形FEMの解析結果の時間ステップ毎に、以下のステップS105~S107の処理が実行され、各時間ステップの物体Mの剛体運動及び弾性変形が求められる。そして剛体運動及び弾性変形の解析は、最後の時間ステップの処理が完了した時点でS108:YES終了される。すなわち、ステップS105においては、物体Mの初期の重心位置及び慣性主軸を原点及び軸ベクトルとする座標系が局所座標系として設定される。またステップS106においては、該当時間ステップでの剛体運動及び弾性変形による各節点の変位{u0}が、非線形FEMの解析結果から取得され、その取得された各節点の変位{u0}が基準変位{u}として設定される。そしてステップS107において、当該時間ステップにおける物体Mの剛体運動及び弾性変形を計算するための剛体運動及び弾性変形計算処理が実行される。剛体運動及び弾性変形計算処理は、図4に示す剛体運動及び弾性変形計算ルーチンの処理を通じて行われる。本ルーチンの処理が開始されると、まずステップS200において、局所座標系における物体Mの剛体モード変位{qr}が下式により計算される。なお、下式においてt[Φr]は、固有ベクトル[Φr]の転置行列を示している。{qr}=t[Φr][M]{u}次に、ステップS201において、各節点の剛体モードの固有ベクトル[Φr]と剛体モード変位{qr}とから、局所座標系における物体Mの重心位置の変位局所座標系の原点からのずれの並進分及び回転分が計算される。そして、ステップS202において、その計算された並進分及び回転分から、局所座標系における、今回の時間ステップにおける物体Mの重心位置の座標値{ri}及びその慣性主軸の回転行列[Ai]が計算される。続いて、ステップS203において、今回の時間ステップにおける物体Mの重心位置及び慣性主軸を原点及び軸ベクトルとする新たな局所座標系が設定される。そしてステップS204において、各節点の弾性変形による変位{u’}が下式により計算される。なお、下式においてt[Ai]は、回転行列[Ai]の転置行列を示している。ここで計算される各節点の弾性変形による変位{u’}は、厳密には、局所座標系における各節点の初期の座標値{p}と、新たな局所座標系における各節点の剛体運動及び弾性変形後の座標値との偏差であり、本ルーチンの処理が最終段階に進むまでは、剛体運動による変位分が含まれている。{u’}=t[Ai]{p}{u}-{ri}-{p}続いて、ステップS205において、ステップS203において設定された新たな局所座標系が局所座標系に設定される。すなわち、全体座標系における局所座標系の原点の座標値{r}及び軸ベクトルの回転行列[A]が下式に従い更新される。またステップS205においては、ステップS204で計算された各節点の弾性変形による変位{u’}が基準変位{u}に設定される。{r}={r}[A]{ri}[A]=[A][Ai]続くステップS206では、更新前後の局所座標系の偏差が、すなわちステップS205での更新前の局所座標系とステップS203で設定された新たな局所座標系との偏差が既定値以下であるか否かが確認される。この確認は、更新前の局所座標系の原点の座標値{ri-1}と新たな局所座標系の原点の座標値{ri}との偏差が既定値α以下であり、かつ更新前の局所座標系の軸ベクトルの回転行列[Ai-1]と新たな局所座標系の軸ベクトルの回転行列[Ai]との偏差が既定値β以下であるか否かによって行われている。ここで、更新前後の局所座標系の偏差が既定値以下でなければS206:NO、処理がステップS200に戻される。そして更新された局所座標系及び基準変位{u}を用いて、ステップS200~ステップS206の処理が再度実行される。なお、以下では、ステップS200からステップS206までの処理のループを収束ループと記載する。一方、更新前後の局所座標系の偏差が既定値以下であればS206:YES、収束ループを抜けて処理がステップS207に進められる。そしてそのステップS207において、全体座標系における現状の局所座標系の原点の座標値{r}及びその軸ベクトルの回転行列[A]が、今回の時間ステップにおける物体Mの剛体運動として保存される。またステップS208において、現状の各節点の基準変位{u}が、今回の時間ステップにおける物体Mの弾性変形として保存される。そしてその後、今回の本ルーチンの処理が終了され、剛体運動及び弾性変形解析ルーチンのステップS108に処理が進められる。次に、エンジンのピストンとクランクシャフトとを繋ぐコネクティングロッドを解析対象とした場合の剛体運動及び弾性変形の解析態様を説明する。図5は、初期のコネクティングロッドの位置及び形状を一点鎖線で、剛体運動及び弾性変形後のコネクティングロッドの位置及び形状を実線で、それぞれ全体座標系上に示したものとなっている。非線形FEMによる解析は、節点を境界とする有限数の要素にコネクティングロッドを分割して行われている。本実施の形態では、こうした非線形FEMの解析結果が得られた後、以下の態様で各時間ステップにおけるコネクティングロッドの剛体運動及び弾性変形が求められる。すなわち、非線形FEMの解析結果が得られると、まず第1のステップとして、全体座標系における各節点の初期の座標値{p0}及び各要素の質量行列[M]が取得される。そして第2のステップとして、コネクティングロッドの初期の重心位置及び慣性主軸を原点及び軸ベクトルとする局所座標系が設定され、その局所座標系における各節点の初期の座標値{p}が計算される。続いて第3のステップとして、質量行列[M]及び局所座標系における各節点の初期の座標値{p}から、各節点の剛体モードの固有ベクトル[Φr]が計算される。次に、コネクティングロッドの各節点の剛体運動及び弾性変形による変位{u0}が非線形FEMの解析結果から取得され、その変位{u0}が各節点の基準変位{u}として設定される。そして下記の第4~第6のステップにより構成された収束ループの処理に移行される。収束ループの処理が開始されると、まず第4のステップとして、固有ベクトル[Φr]、質量行列[M]及び基準変位{u}から、局所座標系におけるコネクティングロッドの剛体モード変位{qr}が計算される。また第5のステップとして、剛体モード変位{qr}及び固有ベクトル[Φr]から、局所座標系における剛体運動及び弾性変形後の物体の重心位置の座標値{ri}及び慣性主軸の回転行列[Ai]がそれぞれ計算される。そして第6のステップとして、剛体運動及び弾性変形後のコネクティングロッドの重心位置及び慣性主軸を原点及び軸ベクトルとする新たな局所座標系が設定される。また、第2のステップで設定された局所座標系における各節点の初期の座標値{p}とその新たな局所座標系における各節点の剛体運動及び弾性変形後の座標値との偏差が、各節点の弾性変形による変位{u’}として計算される。図6aは、この時点の、すなわち1回目の収束ループの処理が終了した時点の、新たな局所座標系におけるコネクティングロッドの位置及び形状を示したものである。この時点でコネクティングロッドの物体の剛体モード変位{qr}を微小≒0とすることができれば、コネクティングロッドの剛体運動は、その重心位置の並進及び回転として上記新たな局所座標系からは除かれており、第6のステップで計算された各節点の弾性変形による変位{u’}には、剛体運動による変位は一切含まれていないことになる。しかしながら、コネクティングロッドの局所座標系の位置及び回転行列は未知であるため、剛体モード変位は微小≒0である保証はなく、この時点で計算された各節点の弾性変形による変位{u’}には、剛体運動分が未だ若干残っている。図6bは、この時点の、新たな局所座標系におけるコネクティングロッドの位置及び形状を、変形を100倍に拡大して示したものである。同図に示すように、この時点では、剛体運動の除去は、未だ不十分となっている。そこで、本実施の形態では、新たな局所座標系を局所座標系に、先に計算した各節点の弾性変形による変位{u’}を各節点の基準変位{u}にそれぞれ設定し直した上で、上記収束ループの処理を再実行している。収束ループの処理が再実行されると、その処理において更新された各節点の基準変位{u}からは、コネクティングロッドの剛体運動分が更に取り除かれる。こうした収束ループの処理の再実行は、更新前後の局所座標系の偏差が既定値以下となるまで、より具体的には更新前後の局所座標系の原点の座標値{r}及び軸ベクトルの回転行列[A]の偏差が既定値以下となるまで繰り返される。したがって、収束ループを抜ける直前に計算された各節点の弾性変形による変位{u’}からは、物体の剛体運動分がほぼ完全に取り除かれるようになる。したがって、収束ループを抜ける直前の各節点の弾性変形による変位{u’}の計算値、すなわちステップS205で更新された基準変位{u}は、各節点の弾性変形による変位をほぼ正確に示すものとなる。また収束ループを抜ける直前に設定された新たな局所座標系の、すなわちステップS205で更新された局所座標系の、全体座標系における原点の座標値{r}及びその軸ベクトルの回転行列[A]は、コネクティングロッドの剛体運動の並進成分及び回転成分をそれぞれ示すものとなる。図7aは、この時点の、すなわち収束ループを抜けた時点の、局所座標系におけるコネクティングロッドの位置及び形状を示したものである。また図7bは、この時点の局所座標系におけるコネクティングロッドの位置及び形状を、変形を10000倍に拡大して示したものである。このように、この時点の局所座標系におけるコネクティングロッドの位置及び形状からは、剛体運動による変位がほぼ完全に除かれている。そのため、このときの局所座標系におけるコネクティングロッドの位置及び形状を見ることで、その弾性変形の様相を容易かつ的確に把握することができる。以上説明した本実施の形態の剛体運動及び弾性変形の解析方法によれば、以下の効果を奏することができる。1剛性行列が時々刻々と変化することから弾性モードが未知の有限要素で表現された物体の場合にも、その剛体変位と弾性変形とを分離して解析することができる。2剛体運動及び弾性変形の解析の前提となるFEMによる物体の変位及び変形の解析に、既存のFEMソフトウェアをそのまま利用することができる。なお、上記実施の形態は、以下のように変更して実施することもできる。上記実施の形態では、物体Mの剛体運動及び弾性変形の時間ステップ毎の過渡応答計算を行うようにしていたが、最終的な物体Mの剛体運動及び弾性変形さえ解れば良い場合には、最終的な物体Mの剛体運動及び弾性変形のみ求めるようにしても良い。すなわち、この場合には、図3の剛体変位及び弾性変形解析ルーチンのステップS106において、最終的な各節点の変位を取得するようにするとともに、ステップS107の処理を完了した時点で同ルーチンの処理を終了するようにすることになる。上記実施の形態では、非線形FEMの解析結果から物体の剛体運動及び弾性変形を求めていたが、線形FEMの解析結果からの物体の剛体運動及び弾性変形の解析にも本発明は同様に適用することができる。