JP4814846B2 - 対象物衝突模擬装置、その方法、そのプログラム、及びその記録媒体 - Google Patents

Description

この発明は、対象物同士が衝突をする際、衝突による変形や変位をリアルに模擬する対象物衝突模擬装置、その方法、そのプログラム、及びその記録媒体に関する。

衝突は日常生活にてよく発生する物理現象である。対象物同士が衝突を起こすと、当該対象物同士は運動量交換を行う。当該運動量交換による各対象物における運動量の増減は、衝撃量(Impulse)としてそれぞれの対象物に作用する。衝撃量からなる衝撃力(Impulse Force)は、各対象物の運動状態を変動させ、衝突前後にて各対象物の加速度は不連続に変わる。

図１Ａに、剛体同士の衝突の様子を示し、図１Ｂに弾性体同士の衝突の様子を示す。弾性体とは、人間や動物の皮膚組織、木材、ゴム、など弾性を有する連続体である。以下では、連続弾性体という場合もある。対象物を剛体と仮定した場合、図１Ａに示すように、剛体同士が衝突すると、衝突によって加速度の不連続転移のみが発生する。対象物を弾性体と仮定した場合、図１Ｂに示すように、弾性体同士が衝突すると、剛体衝突のような衝突前後における加速度の不連続転移に加え、衝突による弾性変形ももたらされる。衝突による弾性変形は、対象物が柔らかいほど、すなわち剛性が低いほど著しくなる。また、衝突前後における加速度の不連続転移は、衝突直前の速度及び対象物の材質によって決まる。この不連続転移における衝突前速度や対象物の材質の影響は運動力学にて定義され、いわば反発係数(Coefficient of Restitution)として定量化し、表現できる。

連続弾性体の動的挙動を計算力学的手法で模擬する際、連続弾性体を複数の離散質点に網状の弾性要素を張って構成した離散要素近似力学系としてモデル化し、解析を行う場合がある。離散要素近似力学系は、例えば有限要素法(FEM)などと比べて、モデル化が容易で、計算時間が短いなどの利点がある。
弾性体の衝突の模擬方法については、コンピュータグラフィックスの分野から、非特許文献１〜３に示すように、いくつかの手法が提案されている。
R. Keiser, M. Muller, B. Heidelberger, M. Teschner, and M. Gross, "Contact Handling for Deformable Point-Based Objects," Proc. of Vision, Modeling, Visualization VMV'04, pp. 339-347, Stanford, USA, Nov. 2004 J. Spillmann and M. Teschner, "Contact Surface Computation for Coarsely Sampled Deformable Objects," Proc. of Vision, Modeling, Visualization VMV'05, pp. 289-296, Erlangen, Germany, Nov. 2005 G. Hirota, S. Fisher, and M. C. Lin "Simulation of non-penetrating Elastic Bodies Using Distance Fields," Technical Report TR00-018, University of North Carolina at Chapel Hill, 2000

非特許文献１〜３記載の対象物の衝突の模擬方法では、対象物同士の衝突をリアルに模擬できないという問題があった。本発明は、複数の対象物の衝突の挙動を衝突の際に行われる対象物間の運動量交換を考慮することで、自然にかつリアルに模擬する対象物衝突模擬装置を提供することを目的とする。

この発明は、複数の対象物をそれぞれ離散質点（以下、「質点」という。）及び粘弾性要素による離散系モデルに近似して、当該対象物同士の衝突を模擬する対象物衝突模擬装置である。当該対象物衝突模擬装置は、衝突確認部と、衝突後等価速度計算部と、全系衝突衝撃力計算部と、弾性衝撃力計算部と、剛性衝撃力計算部と、衝突衝撃力計算部と、適用部と、を備える。衝突確認部は、対象物同士が衝突状態かを確認する。衝突後等価速度計算部は、対象物の衝突後の等価速度ベクトルｖ_ｅｑ−Ｉ（ただし、Ｉは対象物を示す値）を求める。全系衝突衝撃力計算部は、等価速度ベクトルｖ_ｅ−Ｉを用いて、全系衝突衝撃力ベクトルＦ_Ｉを求める。弾性衝撃力計算部は、質点の質量ｍ_Ｉ，ｉ（ただし、ｉは質点を示す値）と、前記全系衝突衝撃力ベクトルのスカラー量│Ｆ_Ｉ│と、質点の法線ベクトルｎ_Ｉ，ｉを用いて、質点ごとの弾性衝撃力ベクトルｆ_ＩＤ，ｉを求める。剛性衝撃力計算部は、全系衝突衝撃力ベクトルのスカラー量│Ｆ_Ｉ│と、衝突後の対象物の等価加速度ベクトルａ_ｅｑ−Ｉを用いて、質点ごとの剛性衝撃力ベクトルｆ_ＩＲ，ｉを求める。衝突衝撃力計算部は、弾性衝撃力ベクトルｆ_ＩＤ，ｉと、剛性衝撃力ベクトルｆ_ＩＲ，ｉを加算して、質点ごとの衝突衝撃力ベクトルｆ_Ｉ，ｉを求める。適用部は、衝突衝撃力ベクトルｆ_ＩＤ，ｉから、あらかじめ定めた時間Δｔ後の質点の速度ベクトルｖ_Ｉ，ｉと位置ベクトルｘ_Ｉ，ｉとを求める。

また、衝突後速度計算部は、反発係数を用いて、衝突後の速度ベクトルを計算してもよい。
更に、質点が貫入増加状態にある場合には、貫入防止部を備えても良い。貫入防止部は、適用部で求めた質点の速度ベクトルｖ_Ｉ，ｉと位置ベクトルｘ_Ｉ，ｉとを、貫入解消状態となるように調整する。

当該貫入防止部は、貫入増加状態質点探索手段と、貫入解消速度計算手段と、置き換え手段と、を備えても良い。貫入増加状態質点探索手段は、質点の速度ベクトルｖ_Ｉ，ｉ（ｐ）と、質点の法線ベクトルｎ_ｉを用いて、全ての前記衝突した質点の中から、貫入増加状態にある質点を探索する。貫入解消速度計算手段は、貫入増加状態にある質点に対して、速度ベクトルｖ_Ｉ，ｉ（ｐ）との内積が負となる貫入解消速度ベクトルｖ’’_Ｉ，ｉ（ｐ）を求める。置き換え手段は、適用部で求めた質点の速度ベクトルｖ_Ｉ，ｉ（ｐ）を、前記貫入解消速度ベクトルｖ’’_Ｉ，ｉ（ｐ）に置き換える。

上記の構成により、衝突の際に各対象物に加わる全系衝突衝撃力を、弾性衝撃力と剛性衝撃力とに分けて考える。弾性衝撃力とは、対象物を弾性体と仮定した場合に弾性変形に寄与する成分である。剛性衝撃力とは、対象物を剛性体と仮定した場合に、衝突前後にわたる加速度の不連続転移を引き起こす成分である。このように、各対象物に加わる衝撃力を弾性衝撃力と剛性衝撃力とに分けて考えることで、リアルで自然に、対象物同士の衝突を模擬できる。

以下に、発明を実施するための最良の形態を示す。なお、同じ機能を持つ構成部や同じ処理を行う過程には同じ番号を付し、重複説明を省略する。

図２は、実施例１の対象物衝突模擬装置２などの機能構成例を示すブロック図であり、図３は、主な処理の流れを示すフローチャートである。対象物衝突模擬装置２は衝突確認部８、衝突後等価速度計算部１０、全系衝突衝撃力計算部１２、弾性衝撃力計算部１４、剛性衝撃力計算部１６、衝突衝撃力計算部１８、適用部２０、とで構成されている。この実施例１では、２つの対象物の衝突について考える。当該２つの対象物を人間の上唇、下唇と仮定する。上唇を対象物Ａとし、下唇を対象物Ｂとする。上唇、下唇は、人間の皮膚組織の一部であり、ゲル状のものとして一定の比重をもち、弾性、粘性を有する連続弾性体である。対象物は剛性体でもよい。

まず対象物入力部４に対象物が電子的に変換されて入力される（ステップＳ１）。当該入力は、例えば、ステレオカメラなどの、光学計測装置などで、対象物を読み取る手法や、図面がある場合、図面から直接読み取る手法などで入力される。図４は上唇と下唇を示した図である。図４には、２つの対象物Ａ、Ｂとしての上唇、下唇、これらの境界線と共に示されている。電子的に入力された対象物Ａ、Ｂは、対象物近似部６に入力される。

対象物近似部６は入力された対象物を複数の質点と粘弾性要素とに近似する（ステップＳ２）。図５に近似された上唇、下唇を示す。図５の記載では、黒丸が質点であり、質点同士を結ぶ線が粘弾性要素であり、質点、粘弾性要素で構成したいわゆる離散要素近似力学系モデルである。次に、離散時刻ｐの初期値を例えば、「０」と設定する。ｐは整数である。また離散時刻ｐの時の対象物の質点の速度ｖ_Ｉ、ｉ（ｐ）の初期値、対象物の質点の位置ｘ_Ｉ、ｉ（ｐ）の初期値を設定する（ステップＳ３）。ここで、左側の添え字は対象物を示す記号であり、右側の添え字は質点を示す記号である。この実施例の説明では、対象物が２つ（上唇と下唇）であるので、Ｉの値は「１」もしくは「２」とする。このことは以下の説明でも同様である。

衝突確認部８は対象物Ａ、Ｂ同士が衝突状態か否かを確認する（ステップＳ４）。以下に、確認処理の一例を説明する。ステップＳ４は、「対象物Ａ、Ｂの共有領域Ｒがあるか否かを判定する過程（ステップＳ５）」と「対象物Ａ、Ｂが衝突しているか否かを判定する過程（ステップＳ６）」とを有する。

対象物Ａ、Ｂ同士が衝突状態か否かを確認するためには、例えば、上記非特許文献１等に記載されているＡＡＢＢ（Axis-Aligned Bounding Boxes (AABB:座標軸と平行な辺を持つ直方体）を用いる。ＡＡＢＢにより、衝突を起こす可能性のある衝突候補の質点群を検知する。図６に、対象物Ａ、対象物Ｂ、それぞれのＡＡＢＢ_Ｉを示す。ただしＩは対象物を示す値（この説明では、Ｉ＝１、２）である。図６に示すように、ＡＡＢＢ_１、ＡＡＢＢ_２の共有領域Ｒにある質点が、衝突候補の質点である。衝突候補の質点が存在するか、つまり、対象物Ａ、Ｂの共有領域Ｒがあるか否かを判定する（ステップＳ５）。

次に、対象物Ａ、Ｂが衝突しているか否かを判定する（ステップＳ６）。図７Ａ、Ｂに共有領域Ｒ内にある質点、つまり、衝突候補の質点を示す。図７Ａは対象物Ａ、Ｂが衝突していない状態を示し、図７Ｂは対象物Ａ、Ｂが衝突している状態を示す。白丸が対象物Ａの質点であり、黒丸が対象物Ｂの質点である。対象物Ａの共有領域Ｒ内にある任意の質点を質点ｉとする。質点ｉから最も近い対象物Ｂの質点ｊを質点ｉの投影点と定義する。質点ｉ、質点ｊの位置ベクトルをそれぞれｘ_ｉ、ｘ_ｊとする。質点ｉと質点ｊの位置ベクトルの差を差ベクトルｚ_ｊｉ（＝ｘ_ｉ−ｘ_ｊ）とする。質点における法線ベクトルは、対象物の領域から外側向きのものを用いる。そして、以下の式（１）を満たす質点ｉと質点ｊとを衝突質点とし、衝突質点が存在すれば、対象物Ａと対象物Ｂとは衝突していると判断する（図７Ｂ参照）。
（ｘ_ｉ−ｘ_ｊ）・ｎ_ｊ＝ｚ_ｊｉ・ｎ_ｊ＜０ （１）
ただし、ｚ_ｊｉ・ｎ_ｊはベクトルｚ_ｊｉとベクトルｎ_ｊとの内積を示す。式（１）を満たすということは、差ベクトルｚ_ｊｉと投影点の法線ベクトルｎ_ｊとが９０度より大きい内角をなしているということである。

説明を図２、図３に戻す。衝突確認部８が、対象物Ａ、Ｂが衝突していると確認すると、衝突後等価速度計算部１０は、対象物Ａ、Ｂの衝突後の等価速度ベクトルｖ_ｅｑ−Ｉを求める（ステップＳ８）。以下、求め方の一例を示す。
まず、対象物Ａ、Ｂの運動量Ｌ_１、Ｌ_２を求める。対象物Ａ、Ｂの運動量Ｌ_１、Ｌ_２は、対象物Ａ、Ｂを構成している全質点の運動量の和とできる。つまり、以下の式（２ａ）（２ｂ）により求める。
Ｌ_１＝Σ_ｉ＝１ ^Ｐｍ_１、ｉｖ_１、ｉ（ｐ） （２ａ）
Ｌ_２＝Σ_ｉ＝１ ^Ｑｍ_２、ｉｖ_２、ｉ（ｐ） （２ｂ）

ただし、Ｐ、Ｑはそれぞれ、対象物Ａ、Ｂの質点の総数であり、ｍ_１、ｉ、ｍ_２、ｉはそれぞれ、対象物Ａの質点ｉの質量、対象物Ｂの質点ｉの質量であり、ｖ_１、ｉ（ｐ）、ｖ_２、ｉ（ｐ）はそれぞれ、離散時刻ｐ（ｐ番目）の対象物Ａの質点ｉの速度ベクトル、離散時刻ｐの対象物Ｂの質点ｉの速度ベクトルである。また、対象物Ａ、Ｂの運動量は、以下の式（３ａ）（３ｂ）を用いても、求めることができる。
Ｌ_１＝Σ_ｉ＝１ ^Ｐｍ_１、ｉｖ_Ｉ、ｉ（ｐ）＝ｖ_ｅｑ−１（ｐ）Σ_ｉ＝１ ^Ｐｍ_１、ｉ
＝Ｍ_１・ｖ_ｅｑ−１（ｐ） （３ａ）
Ｌ_２＝Σ_ｉ＝１ ^Ｑｍ_２、ｉｖ_２、ｉ（ｐ）＝ｖ_ｅｑ−２（ｐ）Σ_ｊ＝１ ^Ｑｍ_２、ｉ
＝Ｍ_２・ｖ_ｅｑ−２（ｐ） （３ｂ）

ただし、ｖ_ｅｑ−１（ｐ）、ｖ_ｅｑ−２（ｐ）はそれぞれ、離散時刻ｐの対象物Ａ、Ｂの等価速度ベクトルで、対象物内の全質点は、同一の等価速度を持つ。Ｍ_１、Ｍ_２はそれぞれ対象物Ａ、Ｂの総質量である。従って、対象物Ａ、Ｂの等価速度ｖ_ｅｑ−１（ｐ）、ｖ_ｅｑ−２（ｐ）は、以下の式（４ａ）（４ｂ）で求めることができる。
ｖ_ｅｑ−１（ｐ）＝Ｌ_１／Ｍ_１ （４ａ）
ｖ_ｅｑ−２（ｐ）＝Ｌ_２／Ｍ_２ （４ｂ）
式（３ａ）（３ｂ）、式（４ａ）（４ｂ）において、対象物Ａ、Ｂの等価速度ｖ_ｅｑ−１（ｐ）、ｖ_ｅｑ−２（ｐ）は、連続弾性体と同じ質量、同じ運動量を持つ剛体の移動速度、すなわち、連続弾性体を剛体に置き換えた場合、両者が力学的に等価の状態にあるための必要速度として解釈できる。

また、運動量保存の法則より、衝突前後の対象物Ａ、Ｂの運動量を以下の式（５）で表現できる。
Ｍ_１ｖ_ｅｑ−１（ｐ）＋Ｍ_２ｖ_ｅｑ−２（ｐ）＝Ｍ_１ｖ’_ｅｑ−１（ｐ）＋Ｍ_２ｖ’_ｅｑ−２（ｐ） （５）
ただし、ｖ’_ｅｑ−１（ｐ）、ｖ’_ｅｑ−２（ｐ）はそれぞれ、対象物Ａ、Ｂの衝突後の等価速度である（以下、衝突後等価速度という）。左辺が衝突前の対象物Ａ、Ｂの運動量の和であり、右辺が衝突後の対象物Ａ、Ｂの運動量の和である。対象物Ａ、Ｂの衝突を考慮する際、反発係数を考慮することが好ましい。当該反発係数ｅは以下の式（６）で表すことができる。
ｅ＝（ｖ_ｅｑ−１（ｐ）−ｖ_ｅｑ−２（ｐ））／（ｖ’_ｅｑ−２（ｐ）−ｖ’_ｅｑ−１（ｐ）） （６）

対象物の材質を考慮して反発係数ｅの値を与え、式（５）、（６）を連立して解けば、対象物Ａ、Ｂの離散時刻ｐにおける衝突後等価速度ｖ’_ｅｑ−１（ｐ）、ｖ’_ｅｑ−２（ｐ）を以下の式（７ａ）、（７ｂ）で表現できる。
ｖ’_ｅｑ−１（ｐ）＝
｛Ｍ_１ｖ_ｅｑ−１（ｐ）＋Ｍ_２ｖ_ｅｑ−２（ｐ）＋Ｍ_２ｅ（ｖ_ｅｑ−２（ｐ）−ｖ_ｅｑ−１（ｐ））｝／（Ｍ_１＋Ｍ_２） （７ａ）
ｖ’_ｅｑ−２（ｐ）＝
｛Ｍ_１ｖ_ｅｑ−１（ｐ）＋Ｍ_２ｖ_ｅｑ−２（ｐ）＋Ｍ_１ｅ（ｖ_ｅｑ−１（ｐ）−ｖ_ｅｑ−２（ｐ））｝／（Ｍ_１＋Ｍ_２） （７ｂ）
全系衝突衝撃力計算部１２は、衝突後等価速度ベクトルｖ’_ｅｑ−１（ｐ）、ｖ’_ｅｑ−２（ｐ）を用いて、全系衝突衝撃力ベクトルＦ_Ｉを求める（ステップＳ１０）。運動量保存の法則により、以下の式（８ａ）、（８ｂ）が成り立つ。
Ｍ_１ｖ_ｅｑ−１（ｐ）＋Ｆ_１Δｔ＝Ｍ_１ｖ’_ｅｑ−１（ｐ） （８ａ）
Ｍ_２ｖ_ｅｑ−２（ｐ）＋Ｆ_２Δｔ＝Ｍ_２ｖ’_ｅｑ−２（ｐ） （８ｂ）

ただし、Ｆ_１は対象物Ａが対象物Ｂから与えられる全系衝突衝撃力であり、Ｆ_２は対象物Ｂが対象物Ａから与えられる全系衝突衝撃力である。また、Δｔは予め定められた時間であるが、式（８ａ）、（８ｂ）ではΔｔを衝突による運動量変換が行われる時間（衝突時間）に近似して用いている。式（８ａ）、（８ｂ）より、全系衝突衝撃力Ｆ_１、Ｆ_２は、以下の式（９ａ）、（９ｂ）で表現できる。
Ｆ_１＝Ｍ_１（ｖ’_ｅｑ−１（ｐ）−ｖ_ｅｑ−１（ｐ））／Δｔ （９ａ）
Ｆ_２＝Ｍ_２（ｖ’_ｅｑ−１（ｐ）−ｖ_ｅｑ−１（ｐ））／Δｔ （９ｂ）

全系衝突衝撃力Ｆ_１、Ｆ_２の質点ごとの衝突衝撃力をｆ_１、ｉ、ｆ_２、ｉ、とすれば、ｆ_Ｉ、ｉを質点ごとの弾性衝撃力ｆ_ＩＤ，ｉと、質点ごとの剛体衝撃力ｆ_ＩＲ，ｉとに分けて考える。弾性衝撃力ｆ_ＩＤ，ｉとは、対象物を弾性体と仮定した場合に弾性変形に寄与する成分である。剛性衝撃力ｆ_ＩＲ，ｉとは、対象物を剛性体と仮定した場合に、衝突前後にわたる加速度の不連続転移を引き起こす成分である。

弾性衝撃力計算部１４は、質点の質量ｍ_Ｉ，ｉと、全系衝突衝撃力ベクトルのスカラー量│Ｆ_Ｉ│と、質点の法線ベクトルｎ_Ｉ，ｉを用いて、質点ごとの弾性衝撃力ベクトルｆ_ＩＤ，ｉを求める（ステップＳ１２）。人間の唇のように高い柔軟性を持つゲル状のものは、流体と同様に圧力変動に対して、境界面変形がもたらされる。こういう柔軟性の高い連続弾性体同士が衝突を起こした場合、衝突による運動量変換から発生する衝突衝撃力は、衝突境界面における衝突体の移動速度を急激に減少させ、衝突体の運動エネルギーもその分、減少してしまう。運動エネルギーの減少は流体におけるエネルギー保存則から内部圧力の上昇を誘発し、この内部圧力の上昇が、境界面変形を引き起こす原因である。以上のように柔軟性の高い連続弾性体の衝突にて、全系衝突衝撃力Ｆ_Ｉは軟組織の内部圧力に変換され、それによる境界面変形がもたらされている。従って、この実施例では、全系衝突衝撃力Ｆ_Ｉを上唇（対象物Ａ）及び下唇（対象物Ｂ）の境界面を構成する質点に対して、法線方向、すなわち境界面に垂直した方向に作用させ、衝突による弾性変形を実現する。従って、上唇（対象物Ａ）の境界面質点ｉの弾性衝撃力ベクトルｆ_１Ｄ，ｉ、下唇（対象物Ｂ）の境界面質点ｉの弾性衝撃力ベクトルｆ_２Ｄ，ｉはそれぞれ以下の式（１０ａ）、（１０ｂ）で求めることができる。
ｆ_１Ｄ，ｉ＝ｋ_Ｄｍ_１，ｉ│Ｆ_１│ｎ_１，ｉ （１０ａ）
ｆ_２Ｄ，ｉ＝ｋ_Ｄｍ_２，ｉ│Ｆ_２│ｎ_２，ｉ （１０ｂ）

ただし、ｎ_１，ｉは対象物Ａ側から見て、物体領域の外側方向に向く、質点ｉの法線ベクトル（図７Ａ参照）であり、ｎ_２，ｉは、対象物２側から見て物体領域の外側方向に向く、質点ｉの法線ベクトルであり、ｋ_Ｄは比例定数であり、│Ｆ_Ｉ│はベクトルＦ_Ｉのスカラー値である。

一方、剛性衝撃力計算部１６は、全系衝突衝撃力ベクトルのスカラー量│Ｆ_Ｉ│と、衝突後の対象物の等価加速度ベクトルａ_ｅｑ−Ｉを用いて、質点ごとの剛性衝撃力ベクトルｆ_ＩＲ，ｉを求める（ステップＳ１４）。対象物Ａ、Ｂが衝突した場合の運動量保存の法則の適用の際、対象物Ａ、Ｂの運動学的状態を等価速度ｖ_ｅｑ−１、ｖ_ｅｑ−２を用いていることから、衝突が対象物Ａ、Ｂにもたらす加速度も以下の式（１１ａ）、（１１ｂ）の様に等価加速度ベクトルａ_ｅｑ−１、ａ_ｅｑ−２として、表すことができる。
ａ_ｅｑ−１＝（ｖ’_ｅｑ−１（ｐ）−ｖ_ｅｑ−１（ｐ））／Δｔ （１１ａ）
ａ_ｅｑ−２＝（ｖ’_ｅｑ−２（ｐ）−ｖ_ｅｑ−２（ｐ））／Δｔ （１１ｂ）

剛体衝撃力は、上述のように、連続弾性体の衝突の問題にて、加速度の不連続転移といった剛体衝突特性の寄与分を意味している。従って、剛体衝撃力ベクトルｆ_１Ｒ，ｉ、ｆ_２Ｒ，ｉは、その方向は加速度ベクトルａ_ｅｑ−１、ａ_ｅｑ−２の方向に従わせ、以下の式（１２ａ）（１２ｂ）で求めることができる。
ｆ_１Ｒ，ｉ＝ｋ_Ｒ│Ｆ_１│ｕ_ａ１ （１２ａ）
ｆ_２Ｒ，ｉ＝ｋ_Ｒ│Ｆ_２│ｕ_ａ２ （１２ｂ）
ただし、ｕ_ａ１、ｕ_ａ２はそれぞれ、ａ_ｅｑ−１、ａ_ｅｑ−２の単位ベクトルであり、ｋ_Ｒは比例定数である。

衝突衝撃力計算部１８は、弾性衝撃力ベクトルｆ_ＩＤ，ｉと、剛性衝撃力ベクトルｆ_ＩＲ，ｉを加算して、質点ごとの衝突衝撃力ベクトルｆ_Ｉ，ｉを求める（ステップＳ１６）。つまり、対象物Ａの質点ｉに作用させる衝突衝撃力ベクトルｆ_１，ｉ、対象物２の質点ｉに作用させる衝突衝撃力ベクトルｆ_２，ｉを以下の式（１３ａ）、（１３ｂ）により求める。
ｆ_１，ｉ＝ｆ_１Ｄ，ｉ＋ｆ_１Ｒ，ｉ （１３ａ）
ｆ_２，ｉ＝ｆ_２Ｄ，ｉ＋ｆ_２Ｒ，ｉ （１３ｂ）

適用部２０は、衝突衝撃力ベクトルｆ_Ｉ，ｉから、あらかじめ定めた時間Δｔ後の質点の速度ベクトルｖ_Ｉ，ｉと位置ベクトルｘ_Ｉ，ｉとを求める（ステップＳ１７）。まず、対象物Ａ、Ｂの質点ｉにそれぞれｆ_１，ｉ、ｆ_２，ｉを加えた場合、離散時刻ｐの次の離散時刻ｐ＋１の速度ｖ_１、ｉ（ｐ＋１）、ｖ_２、ｉ（ｐ＋１）を以下の式（１４ａ）、（１４ｂ）により求める。
ｍ_１，ｉｖ_１、ｉ（ｐ）＋ｆ_１，ｉΔｔ＝ｍ_１，ｉｖ_１、ｉ（ｐ＋１） （１４ａ）
ｍ_２，ｉｖ_２、ｉ（ｐ）＋ｆ_２，ｉΔｔ＝ｍ_２，ｉｖ_２、ｉ（ｐ＋１） （１４ｂ）

そして、全質点の位置ｘ_Ｉ、ｉ（ｐ＋１）について、対象物Ａの質点の速度ベクトルｖ_１、ｉ（ｐ）をｖ_１、ｉ（ｐ＋１）に更新し、対象物２の質点の速度ベクトルｖ_２、ｉ（ｐ）をｖ_２、ｉ（ｐ＋１）に更新する（ステップＳ１８）。現在の離散時刻ｐが「１」進み、ｐ＋１になる（ステップＳ２０）。次に、離散時刻ｐ＋１の対象物Ａ、Ｂの質点の位置ベクトルｘ_１，ｉ（ｐ＋１）、ｘ_２，ｉ（ｐ＋１）を以下の式（１５ａ）（１５ｂ）により求め、更新する（ステップＳ２２）。
ｘ_１，ｉ（ｐ＋１）＝ｘ_１，ｉ（ｐ）＋ｖ_１、ｉ（ｐ）Δｔ （１５ａ）
ｘ_２，ｉ（ｐ＋１）＝ｘ_２，ｉ（ｐ）＋ｖ_２、ｉ（ｐ）Δｔ （１５ｂ）

そして、更新されたｖ_１、ｉ（ｐ）、ｖ_２、ｉ（ｐ）が衝突確認部８に入力され、衝突確認部８、衝突後等価速度計算部１０、全系衝突衝撃力計算部１２、弾性衝撃力計算部１４、剛性衝撃力計算部１６、衝突衝撃力計算部１８、適用部２０の処理が繰り返される。また、更新された位置ベクトルｘ_１，ｉ（ｐ）、ｘ_２，ｉ（ｐ）に対応する情報が表示部２６に表示される。表示部２６は例えばモニタである。このようにして、表示部２６に質点の挙動が表示され、総合的に対象物Ａ、Ｂの挙動が表示される。

ステップＳ５において対象物同士の共有領域Ｒがないと判断された場合、および、ステップＳ６において対象物同士が衝突していないと判断された場合、対象物Ａ、Ｂについて何らの処理を行わず、対象物Ａ、Ｂの質点の位置（挙動）を表示部２６に表示する。
この実施例１の対象物模擬装置２のように、全系衝突衝撃力Ｆ_Ｉを、弾性衝撃力ｆ_ＩＤ，ｉと剛性衝撃力ｆ_ＩＲ，ｉとに分けて考えることで、リアルで自然な対象物同士の衝突を模擬できる。

実施例１の対象物衝突模擬装置２の構成で、対象物の衝突を模擬すると、対象物Ａ、Ｂの質点が貫入するという現象が生じる。貫入現象を簡単に説明する。例えば、対象物Ａと対象物２とが衝突を起こした場合、この衝突に伴って発生する境界面の付近の質点を衝突境界面付近質点とする。貫入現象とは、当該衝突により、例えば、対象物Ａの衝突境界面付近質点が対象物２の領域に入り込む現象である（図７Ｂ参照）。上記非特許文献１〜３には、この貫入現象を防ぐ手法が記載されている。その手法とは、交差が生じた対象物の境界面同士の間に、擬似的に新たな擬似的境界面を１つ生成する。そして、２つの対象物が擬似的境界面で互いに接触していると仮定して、境界面付近質点を強制的に擬似的境界面に移動させて模擬する。

しかし、この手法では、境界面貫入の防止には有効であるが、衝突による時間経過に伴う弾性変形を模擬できない。従って、対象物が高い柔軟性をもつ弾性体である場合、自然な衝突を模擬できない。また、この手法では、擬似擬的境界面の設定では、その位置を直接指定して衝突境界面の変位を形成している。しかし、衝突境界面付近質点から離れたところに位置する質点に時間領域にて運動方程式から逐次に運動状態が発生し、不自然な衝突しか模擬できない。これらの欠点については非特許文献３にも記載されている。

この実施例２の対象物衝突模擬装置３０により、貫入防止を達成し、かつ、対象物が弾性体であっても、衝突を自然に模擬できる。図２を用いて、対象物衝突模擬装置３０（参照番号を括弧書きで示す）を説明し、図３を用いて主な処理の流れを説明する。対象物衝突模擬装置３０は、破線で示す貫入防止部２２を有する点で対象物衝突模擬装置２と異なる。また、図８に貫入防止部２２の機能構成例を示し、図９に貫入防止部２２の主な処理の流れを示す。貫入防止部２２は、貫入増加状態にある質点に対して、適用部２０で求めた質点の速度ベクトルｖ_Ｉ，ｉと位置ベクトルｘ_Ｉ，ｉとを、貫入が解消するように調整する。貫入防止部２２は貫入増加状態質点探索手段２２２、貫入解消速度計算手段２２４、置き換え手段２２６により構成される。以下、貫入を解消させる一例を示す。衝突している対象物の領域内に貫入している質点は、当該質点の速度ベクトルと質点位置における境界面の法線ベクトルとの相対方位から、以下の２つに分類できる。
貫入増加状態・・・衝突質点が更に、衝突している対象物の領域内部に進み、より深い貫入へ移行しようとしている状態
貫入解消状態・・・衝突質点が衝突している対象物の内部から離れ、貫入を解消する方向へと進もうとしている状態

ステップＳ１８において、離散時刻ｐ＋１の全ての質点の速度ベクトルｖ_Ｉ、ｉ（ｐ＋１）を求めた後に、貫入増加状態質点探索手段２２２は、質点の速度ベクトルｖ_Ｉ，ｉ（ｐ＋１）と、質点の法線ベクトルｎ_ｉを用いて、全ての衝突した質点の中から、貫入増加状態にある質点（以下、「貫入増加状態質点」という。）を探索する（ステップＳ２６２）。具体的には、質点の速度ベクトルｖ_Ｉ，ｉ（ｐ＋１）と、質点の法線ベクトルｎ_ｉの内積を用いて探索する。図１０は、対象物Ａと対象物Ｂとが衝突を起こした場合の、対象物Ａと対象物Ｂとの境界線付近の質点を示した図である。質点ｉにおける速度ベクトルｖ_Ｉ，ｉ（ｐ＋１）と法線ベクトルｎ_ｉの内積が正（ｖ_Ｉ，ｉ（ｐ＋１）・ｎ_ｉ＞０）、つまり、ｖ_Ｉ，ｉ（ｐ＋１）とｎ_ｉのなす内角θが９０度未満の場合、質点ｉは貫入増加の状態にあると判断する。また、速度ベクトルｖ_Ｉ，ｉ（ｐ＋１）と法線ベクトルｎ_ｉの内積が０以下（ｖ_Ｉ，ｉ（ｐ＋１）・ｎ_ｉ≦０）、つまり、θが９０度以上の場合、質点ｉは貫入解消の状態にあると判断する。例えば、図１０に示すように、質点１や質点２は貫入増加状態にあり、質点３や質点４は貫入解消状態にある。

貫入解消速度計算手段２２４は、貫入増加状態質点に対して、速度ベクトルｖ_Ｉ，ｉ（ｐ＋１）との内積が負となる貫入解消速度ベクトルｖ’’_Ｉ，ｉ（ｐ＋１）を求める（ステップＳ２６３）。ステップＳ２６３は、速度ベクトルの向きを変更する過程（ステップＳ２６４）と、速度ベクトルのゲインを調整する過程（ステップＳ２６６）とを有する。まず、貫入解消状態になるように、速度ベクトルの向きを変更させる。具体的には、以下の式（１６）が成り立つ速度ベクトルｖ’_Ｉ，ｉ（ｐ＋１）を求める。
ｖ’_Ｉ，ｉ（ｐ＋１）・ｖ_Ｉ，ｉ（ｐ＋１）≦０ （１６）
ここで、計算量を削減するため、ｖ’_Ｉ，ｉ（ｐ＋１）＝−ｖ_Ｉ，ｉ（ｐ＋１）とすることが好ましい。

次に、速度ベクトルｖ’_Ｉ，ｉ（ｐ＋１）のゲインを調整して、貫入解消速度ベクトルｖ’’_Ｉ，ｉ（ｐ＋１）を求める（ステップＳ２６６）。具体的には、以下の式（１７）により調整する。
ｖ’’_Ｉ，ｉ（ｐ＋１）＝ｃｖ_Ｉ，ｉ’（ｐ＋１） （１７）
ただし、ｃは比例定数である。ｃの値は、計算結果が不自然にならない限り１でよい。不自然になる場合は、適切な挙動が得られるように、ｃを調整する。ｃが大きいほど、すばやい貫入解消が実現する。しかし、極端に大きなｃは不自然な衝突挙動を誘発する。そこで、ｃの値は０．５＜ｃ＜１．５にするのがよい。

置き換え手段２２６は、適用部２０で求めた質点の速度ベクトルｖ_Ｉ，ｉ（ｐ＋１）を、貫入解消速度ベクトルｖ’’_Ｉ，ｉ（ｐ＋１）に置き換える（ステップＳ２６８）。この置き換え処理により、貫入増加状態質点は貫入解消状態に向かう。ステップＳ２６８の後はステップＳ２０（図３参照）に進む。
ステップＳ２６２において、貫入増加状態質点がない場合は、速度ベクトルｖ_Ｉ，ｉ（ｐ）の向きを変更する過程（ステップＳ２６４）に進まず、速度ベクトルｖ_Ｉ，ｉ（ｐ＋１）のゲインを調整する過程（ステップＳ２６８）に進む。

そして、貫入解消速度ベクトルｖ’’_Ｉ，ｉ（ｐ＋１）は、衝突確認部８に入力され、衝突確認部８、衝突後等価速度計算部１０、全系衝突衝撃力計算部１２、弾性衝撃力計算部１４、剛性衝撃力計算部１６、衝突衝撃力計算部１８、適用部２０の処理が繰り返される。
この実施例２の対象物衝突模擬装置３０のように、貫入防止部２２を設けることで、貫入防止を達成し、かつ、対象物が弾性体であっても、自然な衝突を模擬できる。実施例１、２では、対象物が２つの場合を説明したが、対象物が３つ以上の場合であっても、対象物衝突模擬装置２、３０を適用できる。

［実験結果］
図１１Ａ〜Ｄは対象物衝突模擬装置２、３０を用いていない離散要素近似力学系モデルを用いた、上唇と下唇との閉口運動シミュレーションの結果である。図１２Ａ〜Ｄは対象物衝突模擬装置３０を用いた上唇と下唇との閉口運動シミュレーションの結果である。図１３Ａ〜Ｄはそれぞれ、図１１Ａ〜Ｄの中央縦断面を示した図であり、図１４Ａ〜Ｄはそれぞれ、図１２Ａ〜Ｄの中央縦断面を示した図である。また、対象物衝突模擬装置３０において、反発係数ｅ＝０とする。図１１Ａ、図１２Ａはｔ＝０．００秒、図１１Ｂ、図１２Ｂはｔ＝０．０４秒、図１１Ｃ、図１２Ｃはｔ＝０．１６秒、図１１Ｄ、図１２Ｄはｔ＝０．４２秒、の口唇の挙動を示している。ｔ＝０．００秒、ｔ＝０．０４秒、ｔ＝０．１６秒、ｔ＝０．４２秒、はそれぞれ初期状態、衝突開始、衝突展開中及びシミュレーション終了の時刻に対応する。このことは以下の図１５〜図２０についても同様である。

図１１Ａ〜Ｄ、図１３Ａ〜Ｄより理解されるように、口唇間衝突のモデルが設けられていない各口唇にて衝突相手、すなわち反対側口唇の存在の検知がされておらず、上唇と下唇が相手口唇の領域内に、ひどく貫入していることが分かる。
一方、図１２Ａ〜Ｄ、図１４Ａ〜Ｄより理解されるように、対象物衝突模擬装置３０を用いると、衝突検知及び衝突挙動のモデリングを設けてから行ったシミュレーションによる衝突挙動で、有効な衝突検知はもちろん、衝突後挙動では、境界面における貫入が防止され、更に衝突に伴う弾性変形も、自然に実現されていることが分かる。

図１５〜図１７は、それぞれ対象物衝突模擬装置３０を用いて、反発係数ｅ＝０．０、ｅ＝０．５、ｅ＝１．０、に設定した場合の、口唇間の衝突挙動シミュレーション結果である。図１８〜図２０は、図１５〜図１７の中央縦断面を示した図である。
ｅ＝０．０の反発係数を持つモデル（図１５、図１８）にて最後まで上唇と下唇の相互接触状態が維持している。これに対して、ｅ＝０．５の反発係数を持つモデル（図１６、図１９）においては、衝突後、上唇と下唇が跳ね返るような挙動でお互いに遠ざかり、最後は開口の状態に戻っていることが分かる。また、ｅ＝１．０の反発係数を持つモデル、つまり、図１６、図１９に示すように、反発係数の増加につれ更に著しくなる。これは反発係数が大きいほど速い速度で跳ね返ってくるという、反発係数の力学的定義とも合致する結果である。

Ａは剛体同士の衝突を示した図であり、Ｂは弾性体同士の衝突を示した図である。 対象物衝突模擬装置２、３０の機能構成例を示したブロック図である。 対象物衝突模擬装置２、３０の主な処理の流れを示した図である。 上唇、下唇を示した図である。 上唇、下唇を離散質点及び粘弾性要素で構成される離散系モデルに近似した図である。 対象物Ａ、ＢのＡＡＢＢを示した図である。 Ａは対象物Ａ、Ｂが衝突していない場合の共有領域Ｒ中の質点を示した図であり、Ｂは対象物Ａ、Ｂが衝突している場合の共有領域Ｒ中の質点を示した図である。 貫入防止部２２の機能構成例を示した図である。 貫入防止部２２の主な処理の流れを示した図である。 貫入増加状態質点を示した図である。 従来の技術を用いた上唇と下唇との衝突模擬の実験結果（以下、衝突模擬実験結果という）であり、Ａはｔ＝０．００秒の時を示し、Ｂはｔ＝０．０４秒の時を示し、Ｃはｔ＝０．１６秒の時を示し、Ｃはｔ＝０．４２秒の時を示す。なお、Ａ〜Ｄの時刻については図１２〜図２０についても同様である。 対象物衝突模擬装置３０を用いた場合の衝突模擬実験結果を示す図。 図１１の中央縦断面を示した図である。 図１２の中央縦断面を示した図である。 対象物衝突模擬装置３０を用い、反発係数ｅ＝０．０の場合の衝突模擬実験結果を示す図である。 対象物衝突模擬装置３０を用い、反発係数ｅ＝０．５の場合の衝突模擬実験結果を示す図である。 対象物衝突模擬装置３０を用い、反発係数ｅ＝１．０の場合の衝突模擬実験結果を示す図である。 図１５の中央縦断面を示した図である。 図１６の中央縦断面を示した図である。 図１７の中央縦断面を示した図である。

Claims (10)

1. 複数の対象物をそれぞれ離散質点（以下、「質点」という。）及び粘弾性要素による離散系モデルで近似して、当該対象物同士の衝突を模擬する対象物衝突模擬装置であって、
   対象物同士が衝突状態かを確認する衝突確認部と、
   前記対象物の衝突後の等価速度ベクトルｖ_ｅｑ−Ｉ（ただし、Ｉは対象物を示す値）を運動量保存の法則により求める衝突後等価速度計算部と、
   前記等価速度ベクトルｖ_ｅｑ−Ｉを用いて、全系衝突衝撃力ベクトルＦ_Ｉを運動量保存の法則により求める全系衝突衝撃力計算部と、
   質点の質量ｍ_Ｉ，ｉ（ただし、ｉは質点を示す値）と、前記全系衝突衝撃力ベクトルのスカラー量│Ｆ_Ｉ│と、質点の法線ベクトルｎ_Ｉ，ｉを用いて、運動量保存の法則と弾性体の衝突を同時に考慮した質点ごとの弾性衝撃力ベクトルｆ_ＩＤ，ｉを求める弾性衝撃力計算部と、
   前記全系衝突衝撃力ベクトルのスカラー量│Ｆ_Ｉ│と、衝突後の対象物の等価加速度ベクトルａ_ｅｑ−Ｉを用いて、質点ごとの剛性衝撃力ベクトルｆ_ＩＲ，ｉ ＝ｋ _Ｒ ｜Ｆ _Ｉ ｜ｕ _ａＩ （但し、ｋ _Ｒ は比例定数であり、ｕ _ａＩ は等価加速度ベクトルａ _ｅｑ−Ｉ の単位ベクトル）を求める剛性衝撃力計算部と、
   前記弾性衝撃力ベクトルｆ_ＩＤ，ｉと、前記剛性衝撃力ベクトルｆ_ＩＲ，ｉを加算して、質点ごとの衝突衝撃力ベクトルｆ_Ｉ，ｉを求める衝突衝撃力計算部と、
   前記衝突衝撃力ベクトルｆ _Ｉ，ｉから、あらかじめ定めた時間Δｔ後の質点の速度ベクトルｖ_Ｉ，ｉと位置ベクトルｘ_Ｉ，ｉとを求める適用部と、
   を備えることを特徴とする対象物衝突模擬装置。
2. 請求項１記載の対象物衝突模擬装置であって、
   前記衝突後等価速度計算部は、反発係数を用いて、衝突後の等価速度ベクトルを計算することを特徴とする対象物衝突模擬装置。
3. 請求項１または２記載の対象物衝突模擬装置であって、
   衝突質点が更に、衝突している対象物の領域内部に進み、より深い貫入へ移行しようとしている貫入増加状態にある場合には、
   前記適用部で求めた質点の速度ベクトルｖ_Ｉ，ｉと位置ベクトルｘ_Ｉ，ｉとを、
   衝突質点が衝突している対象物の内部から離れ、貫入を解消する方向へと進もうとしている貫入解消状態となるように調整する貫入防止部も備えることを特徴とする対象物衝突模擬装置。
4. 請求項３記載の対象物衝突模擬装置であって、
   前記貫入防止部は、
   質点の速度ベクトルｖ_Ｉ，ｉと、質点の法線ベクトルｎ_ｉを用いて、全ての前記衝突した質点の中から、貫入増加状態にある質点を探索する貫入増加状態質点探索手段と、
   前記貫入増加状態にある質点に対して、速度ベクトルｖ_Ｉ，ｉとの内積が負となる貫入解消速度ベクトルｖ’’_Ｉ，ｉを求める貫入解消速度計算手段と、
   前記適用部で求めた質点の速度ベクトルｖ_Ｉ，ｉを、前記貫入解消速度ベクトルｖ’’_Ｉ，ｉに置き換える置き換え手段と、
   を備えることを特徴とする対象物衝突模擬装置。
5. 複数の対象物をそれぞれ離散質点（以下、「質点」という。）及び粘弾性要素による離散系モデルで近似して、当該対象物同士の衝突を模擬する対象物衝突模擬方法であって、
   衝突確認部が、対象物同士が衝突状態かを確認する衝突確認過程と、
   衝突後等価速度計算部が、前記対象物の衝突後の等価速度ベクトルｖ_ｅｑ−Ｉ（ただし、Ｉは対象物を示す値）を運動量保存の法則により求める衝突後等価速度計算過程と、
   全系衝突衝撃力計算部が、前記等価速度ベクトルｖ_ｅｑ−Ｉを用いて、全系衝突衝撃力ベクトルＦ_Ｉを運動量保存の法則により求める全系衝突衝撃力計算過程と、
   弾性衝撃力計算部が、質点の質量ｍ_Ｉ，ｉ（ただし、ｉは質点を示す値）と、前記全系衝突衝撃力ベクトルのスカラー量│Ｆ_Ｉ│と、質点の法線ベクトルｎ_Ｉ，ｉを用いて、運動量保存の法則と弾性体の衝突を同時に考慮した質点ごとの弾性衝撃力ベクトルｆ_ＩＤ，ｉを求める弾性衝撃力計算過程と、
   剛性衝撃力計算部が、前記全系衝突衝撃力ベクトルのスカラー量│Ｆ_Ｉ│と、衝突後の対象物の等価加速度ベクトルａ_ｅｑ−Ｉを用いて、質点ごとの剛性衝撃力ベクトルｆ_ＩＲ，ｉ ＝ｋ _Ｒ ｜Ｆ _Ｉ ｜ｕ _ａＩ （但し、ｋ _Ｒ は比例定数であり、ｕ _ａＩ は等価加速度ベクトルａ _ｅｑ−Ｉ の単位ベクトル）を求める剛性衝撃力計算過程と、
   衝突衝撃力計算部が、前記弾性衝撃力ベクトルｆ_ＩＤ，ｉと、前記剛性衝撃力ベクトルｆ_ＩＲ，ｉを加算して、質点ごとの衝突衝撃力ベクトルｆ_Ｉ，ｉを求める衝突衝撃力計算過程と、
   適用部が、前記衝突衝撃力ベクトルｆ _Ｉ，ｉから、あらかじめ定めた時間Δｔ後の質点の速度ベクトルｖ_Ｉ，ｉと位置ベクトルｘ_Ｉ，ｉとを求める適用過程と、
   を有することを特徴とする対象物衝突模擬方法。
6. 請求項５記載の対象物衝突模擬方法であって、
   前記衝突後等価速度計算過程は、反発係数を用いて、衝突後の等価速度ベクトルを計算することを特徴とする対象物衝突模擬方法。
7. 請求項５または６記載の対象物衝突模擬方法であって、
   衝突質点が更に、衝突している対象物の領域内部に進み、より深い貫入へ移行しようとしている貫入増加状態にある場合には、
   前記適用過程で求めた質点の速度ベクトルｖ_Ｉ，ｉと位置ベクトルｘ_Ｉ，ｉとを、衝突質点が衝突している対象物の内部から離れ、貫入を解消する方向へと進もうとしている貫入解消状態となるように調整する貫入防止過程も有することを特徴とする対象物衝突模擬方法。
8. 請求項７記載の対象物衝突模擬方法であって、
   前記貫入防止過程は、
   貫入増加状態質点探索手段が、質点の速度ベクトルｖ_Ｉ，ｉと、質点の法線ベクトルｎ_ｉを用いて、全ての前記衝突した質点の中から、貫入増加状態にある質点を探索する貫入増加状態質点探索ステップと、
   貫入解消速度計算手段が、前記貫入増加状態にある質点に対して、速度ベクトルｖ_Ｉ，ｉとの内積が負となる貫入解消速度ベクトルｖ’’_Ｉ，ｉを求める貫入解消速度計算ステップと、
   置き換え手段が、前記適用部で求めた質点の速度ベクトルｖ_Ｉ，ｉを、前記貫入解消速度ベクトルｖ’’_Ｉ，ｉに置き換える置き換えステップと、
   を有することを特徴とする対象物衝突模擬方法。
9. 請求項５〜８何れかに記載の対象物衝突模擬方法の各過程をコンピュータに実行させるための対象物衝突模擬プログラム。
10. 請求項９記載の対象物衝突模擬プログラムを記録したコンピュータ読み取り可能な記録媒体。