図面 (/)

技術 シミュレーション方法及びその装置

出願人 キヤノン株式会社
発明者 野田智之
出願日 1996年7月15日 (24年5ヶ月経過) 出願番号 1996-184624
公開日 1998年2月3日 (22年11ヶ月経過) 公開番号 1998-031659
状態 未査定
技術分野 特定用途計算機 複合演算
主要キーワード 主慣性モーメント 常微分方程式 反対称 角速度ベクトル 時間刻み 中間変数 数値解 慣性主軸
関連する未来課題
重要な関連分野

この項目の情報は公開日時点(1998年2月3日)のものです。
また、この項目は機械的に抽出しているため、正しく解析できていない場合があります

図面 (5)

課題

計算機内部の記憶容量の消費量が少なく、計算時間が短く、誤差が小さい剛体物理的な運動シミュレーション方法及びその装置を提供する。

解決手段

剛体の主慣性モーメントと、或る時刻の剛体の傾き12に対応する行列Rと、或る時刻の剛体の傾きの時間微分13と、或る時刻の剛体に働く力のモーメント11とから、次の時刻の剛体の傾き14(及び傾きの時間微分15)を求めるために、行列Rに基づくオイラーラグランジェ方程式から得られる行列Rの2階常微分方程式数値解を、Gearの予測子修正子法等により求める(S4)。

概要

背景

仮想現実感を取り扱う分野では、仮想空間内仮想物体シミュレーションを行なうので、剛体運動方程式数値的に解く数値解法が必要となる。また、分子動力学では分子を剛体として取り扱って分子の集合体の集団的運動を追跡するので、同様に剛体の運動方程式の数値解法が必要となる。剛体は、回転運動並進運動とからなる6つの自由度を持つが、並進運動の運動方程式は単純な形であるので数値的に解くことは容易である。一方、回転の運動方程式は複雑であり、色々な数値解法が研究されてきた。

図4に示すように、従来の数値解法では、力のモーメント角運動量ベクトル時間変化との間に成り立つ微分方程式と、角運動量ベクトルと角速度ベクトル関係式と、角度ベクトルと傾きの時間変化との間に成り立つ微分方程式とを用いて、剛体の回転運動の数値解を得る。そのため、或る時刻の角運動量ベクトル42と力のモーメント41とから次の時刻の角運動量ベクトル45を求める手段S1と、或る時刻の角運動量ベクトル42を或る時刻の角速度ベクトル44に変換する手段S2と、或る時刻の角速度ベクトル44と傾き43とから次の時刻の傾き46を求める手段S3とを備える必要がある。

ところが、剛体の傾きをオイラーの角で表現する方法(コンピュータシミュレーション、上田顯著、書店、P146)では、或る時刻の角速度ベクトル44と傾き43とから次の時刻の傾き46を求める上記手段S3において、角速度ベクトルとオイラーの角の時間変化との間に成り立つ微分方程式が特異点を持つため、誤差が大きくなる場合がある。この問題を回避するため、上記手段S3において傾きを4元数というもので表現する方法がある。

概要

計算機内部の記憶容量の消費量が少なく、計算時間が短く、誤差が小さい剛体の物理的な運動のシミュレーション方法及びその装置を提供する。

剛体の主慣性モーメントと、或る時刻の剛体の傾き12に対応する行列Rと、或る時刻の剛体の傾きの時間微分13と、或る時刻の剛体に働く力のモーメント11とから、次の時刻の剛体の傾き14(及び傾きの時間微分15)を求めるために、行列Rに基づくオイラー・ラグランジェ方程式から得られる行列Rの2階常微分方程式の数値解を、Gearの予測子修正子法等により求める(S4)。

目的

本発明の目的は、計算機内部の記憶容量の消費量が少なく、計算時間が短く、誤差が小さい剛体の物理的な運動のシミュレーション方法及びその装置を提供することにある。

効果

実績

技術文献被引用数
0件
牽制数
0件

この技術が所属する分野

ライセンス契約や譲渡などの可能性がある特許掲載中! 開放特許随時追加・更新中 詳しくはこちら

請求項1

剛体運動方程式数値解を得ることにより剛体の物理的な運動シミュレーションするシミュレーション方法であって、剛体の主慣性モーメントと、或る時刻の剛体の傾きに対応する行列Rと、或る時刻の剛体の傾きの時間微分と、或る時刻の剛体に働く力のモーメントとから、次の時刻の剛体の傾きを求める工程を含み、該工程を繰り返すことで剛体の物理的な運動をシミュレーションすることを特徴とするシミュレーション方法。

請求項2

前記工程は、行列Rに基づくオイラーラグランジェ方程式から得られる行列Rの2階常微分方程式の数値解を求める工程を含むことを特徴とする請求項1記載のシミュレーション方法。

請求項3

剛体の運動方程式の数値解を得ることによって剛体の物理的な運動をシミュレーションするシミュレーション装置であって、剛体の主慣性モーメントと、或る時刻の剛体の傾きに対応する行列Rと、或る時刻の剛体の傾きの時間微分と、剛体に働く力のモーメントとを記憶する記憶手段と、前記剛体の主慣性モーメントと、或る時刻の剛体の傾きに対応する行列Rと、或る時刻の剛体の傾きの時間微分と、剛体に働く力のモーメントとから、次の時刻の剛体の傾きを求めて出力する数値解獲得手段と、得られた前記次の時刻の剛体の傾きにより前記記憶手段の内容を更新する更新手段とを備えることを特徴とするシミュレーション装置。

請求項4

剛体の傾き及び傾きの時間微分の初期値と、剛体に働く力のモーメントの初期値と、剛体の主慣性モーメントの初期値と、時間刻みの値とを入力する入力手段を更に備えることを特徴とする請求項3記載のシミュレーション装置。

請求項5

剛体の運動方程式の数値解を得ることにより剛体の物理的な運動をシミュレーションするためのコンピュータ読み出し可能なプログラムを格納するコンピュータ可読メモリであって、少なくとも、剛体の主慣性モーメントと、或る時刻の剛体の傾きに対応する行列Rと、或る時刻の剛体の傾きの時間微分と、或る時刻の剛体に働く力のモーメントとから、次の時刻の剛体の傾きを求めるために、行列Rに基づくオイラー・ラグランジェ方程式から得られる行列Rの2階常微分方程式の数値解を求める数値解モジュールを含むことを特徴とするコンピュータ可読メモリ。

技術分野

0001

本発明は、剛体の回転の運動方程式数値解法及び剛体の物理的な運動シミュレーション方法及びその装置に関するものである。

背景技術

0002

仮想現実感を取り扱う分野では、仮想空間内仮想物体シミュレーションを行なうので、剛体の運動方程式を数値的に解く数値解法が必要となる。また、分子動力学では分子を剛体として取り扱って分子の集合体の集団的運動を追跡するので、同様に剛体の運動方程式の数値解法が必要となる。剛体は、回転運動並進運動とからなる6つの自由度を持つが、並進運動の運動方程式は単純な形であるので数値的に解くことは容易である。一方、回転の運動方程式は複雑であり、色々な数値解法が研究されてきた。

0003

図4に示すように、従来の数値解法では、力のモーメント角運動量ベクトル時間変化との間に成り立つ微分方程式と、角運動量ベクトルと角速度ベクトル関係式と、角度ベクトルと傾きの時間変化との間に成り立つ微分方程式とを用いて、剛体の回転運動の数値解を得る。そのため、或る時刻の角運動量ベクトル42と力のモーメント41とから次の時刻の角運動量ベクトル45を求める手段S1と、或る時刻の角運動量ベクトル42を或る時刻の角速度ベクトル44に変換する手段S2と、或る時刻の角速度ベクトル44と傾き43とから次の時刻の傾き46を求める手段S3とを備える必要がある。

0004

ところが、剛体の傾きをオイラーの角で表現する方法(コンピュータシミュレーション、上田顯著、書店、P146)では、或る時刻の角速度ベクトル44と傾き43とから次の時刻の傾き46を求める上記手段S3において、角速度ベクトルとオイラーの角の時間変化との間に成り立つ微分方程式が特異点を持つため、誤差が大きくなる場合がある。この問題を回避するため、上記手段S3において傾きを4元数というもので表現する方法がある。

発明が解決しようとする課題

0005

しかしながら、オイラーの角の方法並びに4元数の方法などの上記従来例には、ある時刻での剛体の傾きを求めるために、角運動量ベクトルと角速度ベクトルという中間変数を用いるために、計算機内部の記憶容量を余計に消費するという問題点、並びに、2つの微分方程式を解くために計算時間がかかり、誤差が大きくなるという問題点があった。

0006

本発明の目的は、計算機内部の記憶容量の消費量が少なく、計算時間が短く、誤差が小さい剛体の物理的な運動のシミュレーション方法及びその装置を提供することにある。

課題を解決するための手段

0007

上記目的を達成するために、本発明のシミュレーション方法は、剛体の運動方程式の数値解を得ることにより剛体の物理的な運動をシミュレーションするシミュレーション方法であって、剛体の主慣性モーメントと、或る時刻の剛体の傾きに対応する行列Rと、或る時刻の剛体の傾きの時間微分と、或る時刻の剛体に働く力のモーメントとから、次の時刻の剛体の傾きを求める工程を含み、該工程を繰り返すことで剛体の物理的な運動をシミュレーションすることを特徴とする。ここで、前記工程は、行列Rに基づくオイラー・ラグランジェ方程式から得られる行列Rの2階常微分方程式の数値解を求める工程を含む。

0008

又、本発明のシミュレーション装置は、剛体の運動方程式の数値解を得ることによって剛体の物理的な運動をシミュレーションするシミュレーション装置であって、剛体の主慣性モーメントと、或る時刻の剛体の傾きに対応する行列Rと、或る時刻の剛体の傾きの時間微分と、剛体に働く力のモーメントとを記憶する記憶手段と、前記剛体の主慣性モーメントと、或る時刻の剛体の傾きに対応する行列Rと、或る時刻の剛体の傾きの時間微分と、剛体に働く力のモーメントとから、次の時刻の剛体の傾きを求めて出力する数値解獲得手段と、得られた前記次の時刻の剛体の傾きにより前記記憶手段の内容を更新する更新手段とを備えることを特徴とする。剛体の傾き及び傾きの時間微分の初期値と、剛体に働く力のモーメントの初期値と、剛体の主慣性モーメントの初期値と、時間刻みの値とを入力する入力手段を更に備える。

0009

上記構成において、中間変数を用いずに1つの工程で剛体の傾きを求めることが出来るため、計算機内部の記憶容量の消費量を少なく、計算時間を短く、誤差を小さくすることが出来る。

発明を実施するための最良の形態

0010

以下、本発明の実施の形態について説明する。
<本シミュレーション方法の概略>図1は、本実施の形態における剛体の傾きを求める方法を説明するものである。

0011

図1の手段S4は、角運動量ベクトルと角速度ベクトルという中間変数を用いず、力のモーメントと剛体の傾きの関係を記した剛体の傾きの運動方程式であるところの2階常微分方程式を数値的に解く。すなわち、或る時刻の力のモーメント11と或る時刻の傾き12及び傾きの時間微分13とから、2階常微分方程式を数値的に解くことにより、次の時刻の傾き314及び傾きの時間微分15を求める。

0012

(行列Rの2階常微分方程式の導出)最初に、角運動量ベクトルと角速度ベクトルという中間変数を用いず、力のモーメントと剛体の傾きの関係を記した剛体の傾きの運動方程式を立てる方法の概略を述べる。まず、剛体の傾きに対応する行列Rを用いてラグランジェアンを書き下し、そのラグランジェアンからオイラー=ラグランジェ方程式を導く。このオイラー=ラグランジェ方程式を行列Rの時間による2回微分について整理すると、行列Rの2階常微分方程式が得られ、本実施の形態のシミュレーション装置は、この行列Rの2階常微分方程式の解を求める手段S4を備える。

0013

以下、角運動量ベクトルと角速度ベクトルという中間変数を用いず、力のモーメントと剛体の傾きの関係を記した剛体の傾きの運動方程式を立てる方法を、更に具体的に述べる。剛体の主慣性モーメントをi1 ,i2 ,i3 とし、剛体の傾きを表す行列Rを、

0014

0015

とする。ただし、慣性主軸座標軸に平行であるとき傾きを表す行列Rが単位行列に等しくなるように、傾きを表す行列Rを定める。また、剛体に固定した座標系でrと表されるベクトルが、空間座標系ではRrと表されるように定める。

0016

簡単な考察から、角速度ベクトルΩ=(Ωx ,Ωy ,Ωz )によって、行列

0017

0018

であることが分かる。行列Iを、

0019

0020

とすると、回転のエネルギーは、

0021

0022

であることが簡単な考察から分かる。さて、剛体に働く力のモーメントKは、剛体の一点rに働く力をFとすると、
K=r×F (5)
で表される。

0023

剛体がx軸の周り微小角Δθ1 回転したときに外部から受けた仕事Δφ1 は、

0024

0025

但し、Eは3×3の単位行列である。同様に、剛体がy軸の周りに微小角Δθ2 回転したときに外部から受けた仕事Δφ2 は、
Δφ2 =(r×F)y Δθ2 (7)
剛体がz軸の周りに微小角Δθ3 回転したときに外部から受けた仕事Δφ3 は、
Δθ3 =(r×F)zΔθ3 (8)
従って、(5)式より力のモーメントKは、

0026

0027

と表される。剛体の回転の位置エネルギーφの定義ができるとき、

0028

0029

同様に、

0030

0031

従って、

0032

0033

であることが分かる。さて、ラグランジェアンLは、(4)式より、

0034

0035

である。但し、傾きを表す行列Rは、RRt =Eという関係を常に満たす。従って、拘束条件RRt =Eに対応するLagrangeの未定係数項をラグランジェアンに付け加えると、

0036

0037

但し、未定係数行列λは対称行列である。このラグランジェアンLから、オイラー=ラグランジェ方程式を導き出す。

0038

0039

上式に右からR-1をかけ、λを消去すると、

0040

0041

右辺

0042

0043

は、(13)式により力のモーメントに結びつけられる。右辺が反対称、−RIRt+1/2・Tr(I)が対称であるので、右辺を、

0044

0045

−RIRt +1/2・Tr(I)を、

0046

0047

と置く。(2)式より、

0048

0049

であるので、結局、(17)式を

0050

0051

について解くと、

0052

0053

となる。ところで、

0054

0055

であるので、(20),(21),(22)式より、行列Rの2階常微分方程式、

0056

0057

とすると,

0058

0059

が得られる。
(Gearの予測子修正子法)以下に、Gearの予測子修正子法(Gear,C.W.,1966,Report ANL7126,Argonne National Laboratry.)を例に用いて、この行列Rの2階常微分方程式の数値解を求める工程について延べる。

0060

まず、Gearの予測子修正子法そのものについて説明する。予測子修正子法は、次の時間の状態を予測し、その時間の状態の微分方程式からのずれを用いて、その時間の状態を修正する方法である。Gearの方法は、ある時間の座標と、座標の時間による1〜n階微分値とを保持し、それから次の時間の座標と、座標の時間による1〜n階微分値とを予測し、予測した値を微分方程式に代入して求めた誤差の値を用いて、予測した値を修正する。Gearの予測子修正子法は、4元数の方法でも同様に用いられる。

0061

x”=f(x,x’) (25)
という微分方程式を解く場合のこと考える。

0062

0063

と置く。但し、Δtは時間刻みである。予測子段階は、

0064

0065

と書ける。この行列Aは、テーラー定理に従って決められ、nにのみ依存する。修正子段階では、

0066

0067

このベクトルaは、解が安定になるように決められ、nにのみ依存する。結局、予測子段階と修正子段階とを行なうことで、次の時間のx,x’が求められる。

0068

(行列Rの2階常微分方程式の数値解)次に、Gearの予測子修正子法を例に用いて、行列Rの2階常微分方程式の数値解を求める工程について、更に具体的に述べる。

0069

0070

と置く。ただし、1≦i,j≦3である。予測子段階では、

0071

0072

修正子段階では、

0073

0074

但し、Kp は予測子段階終了時の力のモーメントである。以上の予測子段階と修正子段階とをもって、行列Rの2階常微分方程式の数値解を求める工程とする。尚、本実施例では、数値解を求める工程でGearの方法を用いたが、他の2階常微分方程式を解く方法を用いてもよい。
<本シミュレーション装置の構成例>本実施の形態であるシミュレーション装置は、前記の工程及びその手段を備える。

0075

図2は本実施の形態の方法を処理するための解析装置である。解析装置は、CPU21とRAM22から構成されるコンピュータ本体20と、表示装置であるCRT24、入力装置であるキーボード25,マウス26、Gearの予測子修正子法を例に用いて、行列Rの2階常微分方程式の数値解を求める数値解モジュール271を有する外部記憶装置27、出力装置であるプリンタ28等とが、バス29で接続されて構成されている。

0076

RAM22には、剛体の主慣性モーメントI、剛体の傾き及び傾きの1〜n階時間微分Rtc、剛体に働く力のモーメントK、時間刻みΔtからなるデータ221を格納するワーク領域221と、行列Rの2階常微分方程式を記憶する領域222と、プログラムを前記外部記憶装置からロードするプログラムロード領域223とが含まれる。。

0077

キーボード及びマウスは、剛体の傾き及び傾きの1〜n階時間微分Rtcの初期値、剛体に働く力のモーメントKの初期値、剛体の主慣性モーメントI、時間刻みΔtを入力するのに用いられる。これらのデータは一旦外部記憶装置27に格納しても良いが、計算処理速度のためにRAM22に常駐させる。CRT24やプリンタ28は、シミュレーション結果を出力するために使用される。

0078

<本シミュレーションの手順例>図3は、本シミュレーションの手順例を示すフローチャートである。まず、ステップS10で、剛体の傾き及び傾きの1〜n階時間微分Rtcの初期値、剛体に働く力のモーメントKの初期値、剛体の主慣性モーメントI、時間刻みΔtが入力される。ステップS30で、数値解モジュール271に基づき、傾きの行列Rに基づくオイラー・ラグランジェ方程式から得られる行列Rの2階常微分方程式の数値解を、Gearの予測子修正子法により求める。ステップS40で、求められた次の時刻の傾き及び傾きの時間微分により値を更新する。ステップS50で、その結果を数値で、あるいは画像でCRTやプリンタに出力する。ステップS60で、シミュレーションの終了か否かを判断し、終了でなければステップS30に戻って、時間刻みΔtだけ進んだ時刻の次の時刻の傾き及び傾きの時間微分を繰り返し求める。

0079

なお、本発明は、複数の機器(例えばホストコンピュータインタフェイス機器リーダ,プリンタなど)から構成されるシステムに適用しても、1つの機器からなる装置に適用してもよい。また、本発明の目的は、前述した実施形態の機能を実現するソフトウェアプログラムコードを記録した記憶媒体を、システムあるいは装置に供給し、そのシステムあるいは装置のコンピュータ(またはCPUやMPU)が記憶媒体に格納されたプログラムコードを読出し実行することによっても、達成されることは言うまでもない。この場合、記憶媒体から読出されたプログラムコード自体が前述した実施形態の機能を実現することになり、そのプログラムコードを記憶した記憶媒体は本発明を構成することになる。

0080

プログラムコードを供給するための記憶媒体としては、例えば、フロッピディスクハードディスク光ディスク光磁気ディスクCD−ROM,CD−R,磁気テープ不揮発性メモリカード,ROMなどを用いることができる。また、コンピュータが読出したプログラムコードを実行することにより、前述した実施形態の機能が実現されるだけでなく、そのプログラムコードの指示に基づき、コンピュータ上で稼働しているOS(オペレーティングシステム)などが実際の処理の一部または全部を行い、その処理によって前述した実施形態の機能が実現される場合も含まれることは言うまでもない。

0081

さらに、記憶媒体から読出されたプログラムコードが、コンピュータに挿入された機能拡張ボードやコンピュータに接続された機能拡張ユニットに備わるメモリ書込まれた後、そのプログラムコードの指示に基づき、その機能拡張ボードや機能拡張ユニットに備わるCPUなどが実際の処理の一部または全部を行い、その処理によって前述した実施形態の機能が実現される場合も含まれることは言うまでもない。

0082

本発明を上記記憶媒体に適用する場合、その記憶媒体には、先に説明したフローチャートに対応するプログラムコードを格納することになるが、簡単に説明すると、図2メモリマップ例に示す各モジュールを記憶媒体に格納することになる。すなわち、少なくとも数値解モジュールのプログラムコードを記憶媒体に格納すればよい。

発明の効果

0083

以上に詳説したように、本発明は、従来の方法と異なり、角運動量ベクトルと角速度ベクトルという中間変数を用いないため、中間変数に計算機内部の記憶容量を割り当てないので、計算機内部の記憶容量の消費量を従来より少なくすることができる。また、従来の数値解法では複数の工程を要していたが、本発明では1つの工程で済むため計算速度計算精度が向上する。

図面の簡単な説明

0084

図1本実施の形態に係る計算処理手順を示す図である。
図2本実施の形態に係るシミュレーション装置の構成例を示す図である。
図3本実施の形態に係るシミュレーション装置の動作手順例を示すフローチャートである。
図4従来例に係る計算処理手順を示す図である。

ページトップへ

この技術を出願した法人

この技術を発明した人物

ページトップへ

関連する挑戦したい社会課題

関連する公募課題

ページトップへ

技術視点だけで見ていませんか?

この技術の活用可能性がある分野

分野別動向を把握したい方- 事業化視点で見る -

ページトップへ

おススメ サービス

おススメ astavisionコンテンツ

新着 最近 公開された関連が強い技術

この 技術と関連性が強い人物

関連性が強い人物一覧

この 技術と関連する社会課題

関連する挑戦したい社会課題一覧

この 技術と関連する公募課題

関連する公募課題一覧

astavision 新着記事

サイト情報について

本サービスは、国が公開している情報(公開特許公報、特許整理標準化データ等)を元に構成されています。出典元のデータには一部間違いやノイズがあり、情報の正確さについては保証致しかねます。また一時的に、各データの収録範囲や更新周期によって、一部の情報が正しく表示されないことがございます。当サイトの情報を元にした諸問題、不利益等について当方は何ら責任を負いかねることを予めご承知おきのほど宜しくお願い申し上げます。

主たる情報の出典

特許情報…特許整理標準化データ(XML編)、公開特許公報、特許公報、審決公報、Patent Map Guidance System データ