JPWO2013042234A1 - 物体運動解析装置、物体運動解析方法、及び物体運動解析プログラム - Google Patents

Abstract

第１の物体を表現した複数の粒子を含む第１の粒子群の物理量に関するデータに基づき、前記第１の粒子群に含まれる前記複数の粒子の間を接続する接続要素を生成する生成部と、第２の物体に含まれる第２の粒子の物理量に関するデータを参照し、前記第２の粒子と前記生成された接続要素との最短距離を算出して、前記第１の物体と前記第２の粒子との間に作用する反発力を算出する算出部と、前記算出された反発力を用いて、前記第１の物体及び前記第２の物体の運動を解析する解析部と、を有することを特徴とする物体運動解析装置。

Description

本発明は、物体同士の反発力を算出して物体の運動を解析する物体運動解析装置、物体運動解析方法及び物体解析運動プログラムに関する。

従来、コンピュータ上で物体の運動を再現して挙動を解析する試みが、種々の分野で行われている。コンピュータによる物体の運動の再現結果は、例えば、クラブで打たれたゴルフボールの動き、津波を浴びた建築物の耐性、液体等を収容する容器に取り付けられたゴム膜の挙動等を調べるのに、広く利用することができる。

このような運動解析対象となる物体には、流体や弾性体が含まれる。従来、流体や弾性体等の連続体を解く数値計算手法として、格子をベースにして微分方程式の近似解を求解する有限差分法や有限要素法、有限体積法等が多く用いられてきた。近年では、数値計算をＣＡＥ（Computer Aided Engineering）などの応用分野で活用するため、これらの数値計算手法も発展し、流体と構造物が相互作用する問題が解かれるようになってきた。しかしながら、格子を用いる手法においては、自由表面などの界面の存在する問題や、流体-構造連成問題等の移動境界が発生する場合には、取り扱いが複雑なため、プログラム作成が困難となる場合が多い。

これに対して、物体を粒子の集合として解析する粒子法が知られている。代表的な粒子法としては、ＭＰＳ（Moving Particle Semi-implicit）法やＳＰＨ（Smoothed Particle Hydrodynamics）法等が知られている（例えば、非特許文献１、２参照）。粒子法では、移動境界の取り扱いに特別な処置を必要としないため、粒子法は流体−構造連成問題にも広く用いられるようになっている。

これに関連し、流体を表した第１の粒子毎に、及び流体に接触する壁を表した第２の粒子毎に、速度、位置、圧力に関するデータを記録しておき、接触角を加味して流体の挙動を解析する解析装置が知られている（例えば、特許文献１参照）。

ところで、粒子法の計算では、連続体を粒子分布に分けて、連続体の応力による合力を粒子間の相互作用の力の和として表現する。そして、粒子同士が互いに近づこうとする場合には、圧力が大きくなって反発力が発生し、粒子が余りに密集したり疎になったりしない仕組みになっている。

相互作用は、粒子の分布から連続場を構成するのに用いられるカーネル関数による重ね合わせにより表現される。例えばＳＰＨ法の場合、流体の運動を解く問題は、例えば、次式（１）〜（４）で表される常微分方程式系に帰着される。式（１）、（２）は、流体の運動量保存則を離散化したものであり、式（３）は、質量保存則を離散化したものである。また、式（４）は、流体の方程式であり、ここでは密度と圧力が線形の関係となる単純な例を挙げているが、式（４）とは異なる状態方程式を採用することも出来る。式中、添え字ｉは、粒子の番号を表している。ｉ番目の粒子の位置ベクトル、速度ベクトル、密度、圧力はそれぞれ#ｒ_ｉ、#ｖ_ｉ、ρ_ｉ、ｐ_ｉである。ここで、「#」は、後に続くアルファベットがベクトルであることを示すものとする。また、ｍｊはｊ番目の粒子の質量を示している。

上式（３）におけるＷには、カーネル関数として３次のスプライン関数などが用いられる。当該関数が次式（５）で表される場合、ｈはカーネル関数が値をもつ範囲を表す。

ｈは粒子間の影響半径であり、初期状態の平均粒子間隔の２倍から３倍程度が良く用いられる。βはカーネル関数の全空間積分量が１になるように調整された値であり、二次元の場合は０．７πｈ^２、三次元の場合はπｈ^３と決められる。

ＳＰＨ法では、式（３）より、粒子同士が近づこうとすると密度が大きくなろうとし、式（４）より、密度が大きくなると圧力が大きくなり、式（２）より、圧力が大きい方から小さい方へと力が働くため、粒子が反発力を得るという仕組みになっている。図１は、接近する粒子ｉ、ｊが互いに反発力ｆを得る様子を示す図である。

このように、従来の粒子法による計算手法では、境界に固定した粒子を配置することで、境界に近づこうとする粒子に反発力を発生させ、不透性条件（物体が他の物体の境界を突き抜けないこと）を表現している。流体-構造連成の場合には、構造物の境界粒子が受ける圧力に基づいて、流体と構造物の相互作用を表現している。図２は、流体を構成する粒子Ｄ１が、壁の粒子Ｅ１、Ｅ２、…から反発力ｆを受ける様子を示す図である。

特開２００８−１１１６７５号公報

S.Koshizuka, A.Nobe and Y.Oka, "Numerical Analysis of Breaking Waves Using The Moving Particle Semi-Implicit Method", International Journal for Numerical Method in Fluids, 26:751-769（1998） J.J.Monaghan, "Smoothed Particle Hydrodynamics", Annu. Rev. Astron. Astrophys. 30: 543-74（1992）

しかしながら、粒子法を用いた従来の手法では、流体の粒子が構造物の境界を越えない、という境界条件を保障することができない。図３は、粒子法を用いた従来の手法において、境界条件を満たさなくなる場面を示す図である。従来の手法では、壁側の物体が伸展、変形等して粒子間隔が大きくなり、接近する粒子Ｄ１の影響半径ｈの２倍を超えた場合、壁側の物体を構成する粒子から十分な反発力を受けず、境界を突き抜けてしまうという非物理的な計算結果が出力されてしまう。

１つの側面では、解析結果が非物理的な挙動を示してしまうことを抑制することが可能な物体運動解析装置等を提供することを目的とする。

上記目的を達成するための本発明の一態様は、第１の物体を表現した複数の粒子を含む第１の粒子群の物理量に関するデータに基づき、前記第１の粒子群に含まれる前記複数の粒子の間を接続する接続要素を生成する生成部と、第２の物体に含まれる第２の粒子の物理量に関するデータを参照し、前記第２の粒子と前記生成された接続要素との最短距離を算出して、前記第１の物体と前記第２の粒子との間に作用する反発力を算出する算出部と、前記算出された反発力を用いて、前記第１の物体及び前記第２の物体の運動を解析する解析部と、を有することを特徴とする物体運動解析装置である。

１つの側面では、解析結果が非物理的な挙動を示してしまうことを抑制することが可能な、物体運動解析装置、物体運動解析方法、及び物体運動解析プログラムを提供することができる。

接近する粒子ｉ、ｊが互いに反発力ｆを得る様子を示す図である。 流体を構成する粒子Ｄ１が、壁の粒子Ｅ１、Ｅ２、…から反発力ｆを受ける様子を示す図である。 粒子法を用いた従来の手法において、境界条件を満たさなくなる場面を示す図である。 本発明の第１実施例に係る物体運動解析装置１のハードウエア構成例である。 本発明の第１実施例に係る物体運動解析装置１の機能構成例である。 物体Ｅの外縁に位置する粒子からパネル要素ｉを設定する様子を示す図である。 Ｄ_#ｐ，ｉが０から１の間にある場合の位置関係及び最短距離を示す図である。 Ｄ_#ｐ，ｉが０から１の間に無い場合の位置関係及び最短距離を示す図である。 ポテンシャル関数ψの性質を示す図である。 粒子Ｄ１と物体Ｅとの間に作用する反発力ｆ_Ｄ１Ｅを模式的に示す図である。 本発明の第２実施例に係る物体運動解析装置２の機能構成例である。 物体Ｅの外縁に位置する粒子からパネル要素ｉ（ｊ_１２３、ｊ_２３４、…）を設定する様子を示す図である。 点#ｐから三角形パネル要素ｉを含む平面への射影点を求める様子を示す図である。 点#ｐから三角形パネル要素ｉを含む平面への射影点が存在する領域を場合分けした図である。 第２実施例の反発力算出部４４により実行される処理の流れを示すフローチャートである。 本発明の第３実施例に係る物体運動解析装置３の機能構成例である。 物体運動解析装置３が扱う二次元の流体−構造連成問題を模式的に示す図である。 ゴム膜を分割して分点（粒子）を定義し、これに基づき分点（粒子）を接続するパネル要素を生成する様子を示す図である。 本実施例の物体運動解析装置３により実行される処理の流れを示すフローチャートである。 本適用例において説明した手法を用いて実際にシミュレーションを行った際の各過程の状態を示す図である。

以下、本発明を実施するための形態について、添付図面を参照しながら実施例を挙げて説明する。第１実施例では二次元平面上での運動解析について説明し、第２実施例では三次元空間上での運動解析について説明し、第３実施例では、二次元の流体−構造連成問題における運動解析について説明する。

なお、図面及び数式内ではベクトル表記を採用するが、文章中では、「#」が、後に続くアルファベットがベクトルであることを示すものとする。

＜第１実施例＞
以下、本発明の第１実施例に係る物体運動解析装置、物体運動解析方法及び物体運動解析プログラムについて説明する。

［ハードウエア構成］
図４は、本発明の第１実施例に係る物体運動解析装置１のハードウエア構成例である。物体運動解析装置１は、例えば、ＣＰＵ（Central Processing Unit）１０と、ドライブ装置１２と、補助記憶装置１６と、メモリ装置１８と、インターフェース装置２０と、入力装置２２と、ディスプレイ装置２４と、を備える情報処理装置である。これらの構成要素は、バスやシリアル回線等を介して接続されている。

ＣＰＵ１０は、例えば、プログラムカウンタや命令デコーダ、各種演算器、ＬＳＵ（Load Store Unit）、汎用レジスタ等を有するプロセッサである。

ドライブ装置１２は、記憶媒体１４からプログラムやデータを読み込み可能な装置である。プログラムを記録した記憶媒体１４がドライブ装置１２に装着されると、プログラムが記憶媒体１４からドライブ装置１２を介して補助記憶装置１６にインストールされる。記憶媒体１４は、例えば、ＣＤ−ＲＯＭ（Compact Disk Read Only Memory）、ＤＶＤ（Digital Versatile Disk）、ＵＳＢ（Universal Serial Bus）メモリ等の可搬型の記憶媒体である。また、補助記憶装置１６は、例えば、ＨＤＤ（Hard Disk Drive）やフラッシュメモリである。

プログラムのインストールは、上記のように記憶媒体１４を用いる他、インターフェース装置２０がネットワークを介して他のコンピュータよりダウンロードし、補助記憶装置１６にインストールすることによって行うこともできる。ネットワークは、インターネット、ＬＡＮ（Local Area Network）、無線ネットワーク等である。また、プログラムは、情報処理装置の出荷時に、予め補助記憶装置１６やＲＯＭ（Read Only Memory）等に格納されていてもよい。

このようにしてインストール又は予め格納されたプログラムをＣＰＵ１０が実行することにより、図４に示す態様の情報処理装置が、本実施例の物体運動解析装置１として機能することができる。

メモリ装置１８は、例えば、ＲＡＭ（Random Access Memory）やＥＥＰＲＯＭ（Electrically Erasable and Programmable Read Only Memory）である。インターフェース装置２０は、上記ネットワークとの接続等を制御する。

入力装置２２は、例えば、キーボードやマウス、タッチパッド、タッチパネル、マイク等である。また、ディスプレイ装置２４は、例えば、ＬＣＤ（Liquid Crystal Display）やＣＲＴ（Cathode Ray Tube）等の表示装置である。物体運動解析装置１は、ディスプレイ装置２４の他、プリンタ、スピーカ等の他の種類の出力装置を備えてもよい。

［機能構成］
図５は、本発明の第１実施例に係る物体運動解析装置１の機能構成例である。物体運動解析装置１は、パラメータ入力部３０と、パネル要素生成部３２と、反発力算出部３４と、運動解析部３６と、を備える。これらの機能ブロックは、補助記憶装置１６等に格納されたプログラム・ソフトウエアをＣＰＵ１０が実行することにより機能する。なお、これらの機能ブロックが明確に分離したプログラムによって実現される必要はなく、サブルーチンや関数として他のプログラムによって呼び出されるものであってもよい。また、機能ブロックの一部が、ＩＣ（Integrated Circuit）やＦＰＧＡ（Field Programmable Gate Array）等のハードウエア手段であっても構わない。

パラメータ入力部３０は、解析対象の物体を構成する粒子の座標、速度、結合関係、質量等に関する物理量のデータの入力を受け付け、メモリ装置１８に格納する。これらのデータは、予め補助記憶装置１６に格納されていたものであってもよいし、ユーザが入力装置２２を用いて入力したものであってもよい。また、これらのデータは、記憶媒体１４から補助記憶装置１６にインストールされたものであってもよいし、インターフェース装置２０によりネットワークから取得されたものであってもよい。

入力される解析対象の物体の態様は如何なる態様であってもよい。本実施例では、特に、ある平面上に存在する物体Ｄに含まれる粒子Ｄ１と、同一平面上に存在する他の物体Ｅに含まれる複数の粒子（Ｅ１、Ｅ２、Ｅ３、…）と、の間に作用する反発力を算出する手法について説明する。係る反発力以外の物理量（重力や速度、回転速度等）については、他の一般的な物理的手法を用いて算出すればよい。

なお、以下では物体Ｄが単一粒子の物体であるかのように説明するが、例えば第３実施例で説明するように、物体Ｄに含まれる複数の粒子が、それぞれ反発力を物体Ｅから受けるものとして運動を解析することも可能である。

パネル要素生成部３２は、メモリ装置１８に格納されたデータを読み込み、物体Ｅの外縁（境界線）に位置する粒子（Ｅ１、Ｅ２、Ｅ３、…）を抽出してメモリ装置１８に格納する。そして、パネル要素生成部３２は、抽出した粒子を直線で連結し、粒子間を接続する各線分を接続要素である各パネル要素ｉと定義してメモリ装置１８に格納する。図６は、物体Ｅの外縁に位置する粒子からパネル要素ｉ（ｉ_１２、ｉ_２３、ｉ_３４、…）を設定する様子を示す図である。

以下、パネル要素ｉは、始点の座標#ｒ_ｓ，ｉと、終点の座標#ｒ_ｅ，ｉで特定されるものとする。

反発力算出部３４は、まず、粒子Ｄ１と、上記設定した各パネル要素との最短距離を算出する。ここでは、粒子Ｄ１の位置を点#ｐとする。反発力算出部３４は、例えば次式（６）〜（８）で表される量を算出し、粒子Ｄ１とパネル要素ｉとの位置関係を定義する。

式（８）におけるＤ_#ｐ，ｉが０＜Ｄ#ｐ，ｉ＜１の範囲にある場合は、点#ｐからパネル要素ｉへの垂線が、パネル要素ｉと交わるような位置関係にある。この場合、反発力算出部３４は、パネル要素ｉに落とした垂線の長さを、点#ｐとパネル要素ｉとの最短距離として採用する。一方、反発力算出部３４は、Ｄ_#ｐ，ｉがＤ#ｐ，ｉ≦０又は１≦Ｄ#ｐ，ｉの範囲にある場合は、パネル要素ｉの両端である点#ｒ_ｓ，ｉと点#ｒ_ｅ，ｉのうち点#ｐから近い方の点と点#ｐとの距離を、点#ｐとパネル要素ｉとの最短距離として採用する。

図７は、Ｄ_#ｐ，ｉが０＜Ｄ#ｐ，ｉ＜１の範囲にある場合の位置関係及び最短距離を示す図であり、図８は、Ｄ_#ｐ，ｉがＤ#ｐ，ｉ≦０又は１≦Ｄ#ｐ，ｉの範囲にある場合の位置関係及び最短距離を示す図である。

点#ｐとパネル要素ｉとの最短距離ｄ_#ｐ，２（#ｒ_ｓ，ｉ，#ｒ_ｅ，ｉ）は、次式（９）、（１０）によって求められる。

反発力算出部３４は、各パネル要素（ｉ＝１〜ｎ）について、点#ｐとの最短距離ｄ_#ｐ，２（#ｒ_ｓ，ｉ，#ｒ_ｅ，ｉ）を計算する。そして、反発力算出部３４は、次式（１１）により、各パネル要素（ｉ＝１〜ｎ）についての点＃ｐとの最短距離ｄ_#ｐ，２（#ｒ_ｓ，ｉ，#ｒ_ｅ，ｉ）のうち最も値の小さい距離ｄ_{#ｐ，ｍｉｎ，２}を抽出し、距離ｄ_{#ｐ，ｍｉｎ，２}に基づき次式（１２）によりポテンシャル関数ψを定義する。図９は、ポテンシャル関数ψの性質を示す図である。図示するようにポテンシャル関数ψは、距離ｄ_{#ｐ，ｍｉｎ，２}がゼロに近づくと値及び傾斜が大きくなり、距離ｄ_{#ｐ，ｍｉｎ，２}がｒ_０に近づくとゼロに漸近する。更に、反発力算出部３４は、ポテンシャル関数ψを微分することによって、粒子Ｄ１と物体Ｅとの間に作用する反発力ｆ_Ｄ１Ｅを算出する。図９に示すポテンシャル関数ψの性質から、反発力ｆ_Ｄ１Ｅは、距離ｄ_{#ｐ，ｍｉｎ，２}が減少してゼロに近づくと値が大きくなり、距離ｄ_{#ｐ，ｍｉｎ，２}が増加してｒ_０に近づくとゼロに近づく傾向を有している。図１０は、粒子Ｄ１と物体Ｅとの間に作用する反発力ｆ_Ｄ１Ｅを模式的に示す図である。図中、ｉ_４５が、粒子Ｄ１との最短距離が最も小さい（すなわち粒子Ｄ１に最も近い）パネル要素である。

運動解析部３６は、反発力算出部３４によって算出された反発力ｆ_Ｄ１Ｅ、及びその他の物理量（重力や速度、回転速度等）を用いて、物体Ｄ及びＥの運動を解析する。物体に作用する力に基づく運動解析については、既に種々の手法が公知となっているため、説明を省略する。

ここで、従来の手法との比較について説明する。前述したように、粒子法を用いた従来の手法では、壁側の物体が伸展、変形等して粒子間隔が大きくなると、接近する粒子が壁側の物体を構成する粒子から十分な反発力を受けず、境界を突き抜けてしまうという非物理的な計算結果が出力されてしまう場合がある。また、これを避けるために、構造物の境界に多くの粒子を配置すると、計算量が大幅に増加してしまう。

この点、本実施例の物体運動解析装置１では、物体Ｄから、物体Ｅの粒子間を接続するパネル要素までの最短距離に基づいて反発力を計算するため、仮に物体Ｅが大きく変形した場合でも、物体が境界を突き抜けるような計算結果は出力されない。また、物体Ｅの外縁（境界線）に多くの粒子を設定しなくとも物体Ｅの外縁をパネル要素で隙間無く覆うことができるため、計算コストの増加を抑制することができる。従って、非物理的な挙動を防止しつつ迅速に結果を出力することができる。

以上説明した本実施例の物体運動解析装置、物体運動解析方法及び物体運動解析プログラムによれば、非物理的な挙動を防止しつつ迅速に結果を出力することができる。

＜第２実施例＞
以下、本発明の第２実施例に係る物体運動解析装置、物体運動解析方法及び物体運動解析プログラムについて説明する。

［ハードウエア構成］については、第１実施例と同様であるため、説明を省略する。

［機能構成］
図１１は、本発明の第２実施例に係る物体運動解析装置２の機能構成例である。物体運動解析装置２は、パラメータ入力部４０と、パネル要素生成部４２と、反発力算出部４４と、運動解析部４６と、を備える。これらの機能ブロックは、補助記憶装置１６等に格納されたプログラム・ソフトウエアをＣＰＵ１０が実行することにより機能する。なお、これらの機能ブロックが明確に分離したプログラムによって実現される必要はなく、サブルーチンや関数として他のプログラムによって呼び出されるものであってもよい。また、機能ブロックの一部が、ＩＣ（Integrated Circuit）やＦＰＧＡ（Field Programmable Gate Array）等のハードウエア手段であっても構わない。

パラメータ入力部４０は、物体を構成する粒子の座標、速度、結合関係、質量等に関するデータの入力を受け付け、メモリ装置１８に格納する。本実施例のパラメータ入力部４０が受け付ける物理量のデータは、予め補助記憶装置１６に格納されていたものであってもよいし、ユーザが入力装置２２を用いて入力したものであってもよい。また、これらのデータは、記憶媒体１４から補助記憶装置１６にインストールされたものであってもよい。入力される物体の態様は如何なる態様をもとり得るが、本実施例では、特に、ある空間に存在する物体Ｄに含まれる粒子Ｄ１と、同一空間に存在する他の物体Ｅに含まれる複数の粒子（Ｅ１、Ｅ２、Ｅ３、…）と、の間に作用する反発力を算出する手法について説明する。係る反発力以外の物理量（重力や速度、回転速度等）については、他の一般的な物理的手法を用いて算出すればよい。

パネル要素生成部４２は、メモリ装置１８に格納されたデータを読み込み、物体Ｅの外縁（境界面）に位置する粒子（Ｅ１、Ｅ２、Ｅ３、…）を抽出してメモリ装置１８に格納する。そして、パネル要素生成部３２は、抽出した粒子をメッシュ状に連結し、三粒子間を接続する線分によって区画される三角形の面を、接続要素であるパネル要素ｉと定義してメモリ装置１８に格納する。図１２は、物体Ｅの外縁に位置する粒子からパネル要素ｉ（ｉ_１２３、ｉ_２３４、…）を設定する様子を示す図である。なお、三次元で解析を行う場合のパネル要素ｉは、三角形に限らず、四角形その他の形状を有してもよい。

以下、パネル要素ｉは、三点の座標#ｒ_１，ｉ，#ｒ_２，ｉ，#ｒ_３，ｉで特定されるものとする。

反発力算出部４４は、まず、粒子Ｄ１と、上記設定された各パネル要素との最短距離を算出する。以下、第１実施例と同様、粒子Ｄ１の位置を点#ｐとする。反発力算出部４４は、以下の手順で、点#ｐと、三点で構成される閉領域である三角形パネル要素ｉとの最短距離を算出する。

まず、反発力算出部４４は、点#ｐから三角形パネル要素ｉを含む平面への射影点#ｈ_ｉを求める。図１３は、点#ｐから三角形パネル要素ｉを含む平面への射影点#ｈ_ｉを求める様子を示す図である。

反発力算出部４４は、射影点#ｈ_ｉが三角形パネル要素ｉに含まれる場合は、点#ｐから射影点#ｈ_ｉまでの長さを最短距離とする。

一方、反発力算出部４４は、射影点#ｈ_ｉが三角形パネル要素ｉに含まれない場合は、点#ｐから、三角形パネル要素ｉの外縁線において射影点#ｈ_ｉに最も近い点までの距離を最短距離とする。

以下、上記の手順を数式により表す。まず、点#ｐから三角形パネル要素ｉを含む平面への射影点#ｈ_ｉは、次式（１３）により表される。

ここで、ｄ#ｒ_２，ｉ＝ｒ_２，ｉ−ｒ_１，ｉ、ｄ#ｒ_３，ｉ＝ｒ_３，ｉ−ｒ_１，ｉとすると、#ｐ−＃ｈ_ｉは、ｄ#ｒ_２，ｉとｄ#ｒ_３，ｉに対して垂直であるから、次式（１４）、（１５）が成立する。この両式を満たすようにパラメータｓ_２、ｓ_３を決定することで、射影点#ｈ_ｉを求めることができる。

反発力算出部４４は、射影点#ｈ_ｉを求めると、以下のように場合分けして、各点#ｐと三角形パネル要素ｉとの最短距離を算出する。

（１）ｓ_２≧０，ｓ_３≧０，且つｓ_２＋ｓ_３≧１の場合は、射影点#ｈ_ｉが三角形パネル要素ｉに含まれるため、点#ｐと各三角形パネル要素ｉとの最短距離ｄ_#ｐ，３（#ｒ_１，ｉ，#ｒ_２，ｉ，#ｒ_３，ｉ）は、次式（１６）で表される。

（２）ｓ_２＜０，ｓ_３＜０，且つｓ_２＋ｓ_３≧１の場合は、射影点#ｈ_ｉが図１４における領域Ｑ１に存在することになる。図１４は、粒子Ｅ１の位置から各三角形パネル要素ｉを含む平面への射影点が存在する領域を場合分けした図である。この場合、点#ｐと各三角形パネル要素ｉとの最短距離ｄ_#ｐ，３（#ｒ_１，ｉ，#ｒ_２，ｉ，#ｒ_３，ｉ）を求める問題は、点#ｐと線分の最短距離を求める問題となる。従って、反発力算出部４４は、第１実施例の手法を援用し、#ｒ_ｓ，ｉ＝#ｒ_２，ｉ、#ｒ_ｅ，ｉ＝#ｒ_３，ｉと代入して求めた次式（１７）の値を、点#ｐと接続要素である三角形パネル要素ｉとの最短距離ｄ_#ｐ，３（#ｒ_１，ｉ，#ｒ_２，ｉ，#ｒ_３，ｉ）とする。

（３）ｓ_２≧０，ｓ_３＜０，且つｓ_２＋ｓ_３≧１の場合は、射影点#ｈ_ｉが図１４における領域Ｑ２に存在することになる。この場合も、点#ｐと各三角形パネル要素ｉとの最短距離ｄ_#ｐ，３（#ｒ_１，ｉ，#ｒ_２，ｉ，#ｒ_３，ｉ）を求める問題は、点#ｐと線分の最短距離を求める問題となる。従って、反発力算出部４４は、第１実施例の手法を援用し、#ｒ_ｓ，ｉ＝#ｒ_１，ｉ、#ｒ_ｅ，ｉ＝#ｒ_２，ｉと代入して求めた次式（１８）の値を、点#ｐと接続要素である三角形パネル要素ｉとの最短距離ｄ_#ｐ，３（#ｒ_１，ｉ，#ｒ_２，ｉ，#ｒ_３，ｉ）とする。

（４）ｓ_２＜０，ｓ_３≧０，且つｓ_２＋ｓ_３≧１の場合は、射影点#ｈ_ｉが図１４における領域Ｑ３に存在することになる。この場合も、点#ｐと各三角形パネル要素ｉとの最短距離ｄ_#ｐ，３（#ｒ_１，ｉ，#ｒ_２，ｉ，#ｒ_３，ｉ）を求める問題は、点#ｐと線分の最短距離を求める問題となる。従って、反発力算出部４４は、第１実施例の手法を援用し、#ｒ_ｓ，ｉ＝#ｒ_３，ｉ、#ｒ_ｅ，ｉ＝#ｒ_１，ｉと代入して求めた次式（１９）の値を、点#ｐと接続要素である三角形パネル要素ｉとの最短距離ｄ_#ｐ，３（#ｒ_１，ｉ，#ｒ_２，ｉ，#ｒ_３，ｉ）とする。

（５）上記（１）〜（４）以外の場合は、射影点#ｈ_ｉと各三角形パネル要素ｉとの最短距離ｄ_#ｐ，３（#ｒ_１，ｉ，#ｒ_２，ｉ，#ｒ_３，ｉ）は、射影点#ｈ_ｉと各三角形パネル要素ｉのいずれかの頂点との距離となる。従って、反発力算出部４４は、次式（２０）に基づき点#ｐと各三角形パネル要素ｉとの最短距離ｄ_#ｐ，３（#ｒ_１，ｉ，#ｒ_２，ｉ，#ｒ_３，ｉ）を算出する。

反発力算出部４４は、以上の手順によって、各三角形パネル要素ｉについて、点#ｐとの最短距離ｄ_#ｐ，３（#ｒ_１，ｉ，#ｒ_２，ｉ，#ｒ_３，ｉ）を算出する。そして、次式（２１）により、各三角形パネル要素ｉについての点＃ｐとの最短距離ｄ_#ｐ，３（#ｒ_１，ｉ，#ｒ_２，ｉ，#ｒ_３，ｉ）のうち、最も値の小さい距離ｄ_{#ｐ，ｍｉｎ，３}を算出する。更に、反発力算出部４４は、距離ｄ_{#ｐ，ｍｉｎ，３}に基づき、次式（２２）によりポテンシャル関数ψを定義し、ポテンシャル関数ψを微分することによって、粒子Ｄ１と物体Ｅとの間に作用する反発力ｆ_Ｄ１Ｅを算出する。

運動解析部４６は、反発力算出部４４によって算出された反発力ｆ_Ｄ１Ｅ、及びその他の物理量（重力や速度、回転速度等）を考慮し、物体Ｄ及びＥの運動を解析する。物体に作用する力に基づく運動解析については、既に種々の手法が公知となっているため、説明を省略する。

以下、上記説明した反発力算出部４４の処理を、フローチャートの形式で説明する。図１５は、第２実施例の反発力算出部４４により実行される処理の流れを示すフローチャートである。本フローでは、物体Ｄが複数の粒子（Ｄ１、Ｄ２、Ｄ３、…）を含むものとする。

反発力算出部４４は、Ｓ１００〜Ｓ１１２の処理を、物体Ｄを構成する各粒子Ｄｊ（ｊ＝１〜ｋ）について実行する。また、反発力算出部４４は、Ｓ１００〜Ｓ１０８の処理を、各三角形パネル要素ｉ（ｉ＝１〜ｎ）について実行する。

反発力算出部４４は、粒子Ｄｊから三角形パネル要素ｉを含む平面への射影点＃ｈ_ｉを求める（Ｓ１００）。

次に、反発力算出部４４は、射影点が三角形パネル要素ｉに含まれるか否かを判定する（Ｓ１０２）。

反発力算出部４４は、射影点が三角形パネル要素ｉに含まれる場合には、粒子Ｄｊから射影点までの距離を、粒子Ｄｊから三角形パネル要素ｉまでの最短距離とする（Ｓ１０４）。

反発力算出部４４は、射影点が三角形パネル要素ｉに含まれない場合には、三角形パネル要素ｉの外縁線（頂点を含む）上で、粒子Ｄｊから最も近い点を抽出する（Ｓ１０６）。次に、反発力算出部４４は、粒子Ｄｊから抽出した点までの距離を、粒子Ｄｊから三角形パネル要素ｉまでの最短距離とする（Ｓ１０８）。

反発力算出部４４は、Ｓ１００〜Ｓ１０８の処理を全ての三角形パネル要素ｉについて実行すると、各三角形パネル要素ｉについて算出した最短距離の中で、最も値が小さいものを抽出する（Ｓ１１０）。

そして、反発力算出部４４は、抽出した値を用いて、粒子Ｄｊと物体Ｅと間の反発力を表すポテンシャル関数を算出する（Ｓ１１２）。

このような処理によって奏される効果は、第１実施例と同様である。従って、従来技術との比較等については、第１実施例の説明を三次元に拡張すればよいため、説明を省略する。

以上説明した本実施例の物体運動解析装置、物体運動解析方法、及び物体運動解析プログラムによれば、非物理的な挙動を防止しつつ迅速に結果を出力することができる。

＜第３実施例＞
以下、本発明の第３実施例に係る物体運動解析装置、物体運動解析方法、及び物体運動解析プログラムについて説明する。第３実施例の物体運動解析装置３は、第１実施例の物体運動解析装置１が有する機能を基本とし、二次元の流体−構造連成問題に適用可能となるように機能を拡張したものである。

ハードウエア構成については、第１実施例や第２実施例と同様であるため、説明を省略する。

［機能構成］
図１６は、本発明の第３実施例に係る物体運動解析装置３の機能構成例である。物体運動解析装置３は、パラメータ入力部５０と、粒子点設定部５１と、パネル要素生成部５２と、反発力算出部５４と、運動解析部５６と、を備える。これらの機能ブロックは、補助記憶装置１６等に格納されたプログラム・ソフトウエアをＣＰＵ１０が実行することにより機能する。なお、これらの機能ブロックが明確に分離したプログラムによって実現される必要はなく、サブルーチンや関数として他のプログラムによって呼び出されるものであってもよい。また、機能ブロックの一部が、ＩＣ（Integrated Circuit）やＦＰＧＡ（Field Programmable Gate Array）等のハードウエア手段であっても構わない。なお、各機能ブロックの機能については、後述する図１９のフローに則して説明する。

物体運動解析装置３は、二次元のゴム膜（第１実施例における物体Ｅの一例）と流体（第１実施例における物体Ｄの一例）が相互作用する支配方程式系の問題を取り扱うことができる。図１７は、物体運動解析装置３が扱う二次元の流体−構造連成問題を模式的に示す図である。図示するように、物体運動解析装置３は、ゴム膜Ｅの上に純水等の流体Ｄが置かれた系における、ゴム膜Ｅや流体Ｄの運動を解析する。以下、この系における高さをｙ座標で、横方向の位置をｘ座標で表すものとする。

まず、物体運動解析装置３が扱う支配方程式を、次式（２３）〜（２６）で定義する。式中、ｕ（ｘ，ｔ）はゴム膜を構成する各粒子の時間ｔにおける位置ｘに対する高さを表す関数であり、ρ（ｘ，ｙ，ｔ）、ｐ（ｘ，ｙ，ｔ）、#ｖ（ｘ，ｙ，ｔ）はそれぞれ、流体の密度場、圧力場、速度場である。#Ｆは外力であり、ここでは重力を表すものとする。式（２３）はゴム膜の運動量保存則であり、式（２４）は流体の質量保存則であり、式（２５）は流体の運動量保存則であり、式（２６）は流体の状態方程式である。

また、Ｐｕは、流体がゴム膜を押す力を表す。この項は以下のように導出される。（ｘ−ｄｘ／２）と（ｘ＋ｄｘ／２）の間にある膜の微小要素に作用する力は、−ｐ（ｘ，ｕ（ｘ））ｄ#Ｓ・#ｙ＾で表される。ここで、ｄ#Ｓは上記微小要素の面積を有し、向きが外向き法線方向の面積ベクトルであり、#ｙ＾はｙ方向の単位ベクトルである。また、ｄ#Ｓ＝（−ｕ_ｘ，１）ｄｘであるため、Ｐｕは、単位長さあたりの力として表現すると、次式（２７）で表される。

Ｐｕ＝−ｐ（ｘ，ｕ（ｘ）） …（２７）。

ここで、∂ｕ／∂ｔ＝ηとし、ｕを、空間についてはｄｘでｍ等分すると共に、時間についてはｄｔで一様に分割し、ｕ（ｉ×ｄｘ，ｎ×ｄｔ）＝ｕ_ｉ ^ｎという表記を用いると、上式（２３）を、次式（２８）〜（３０）のように離散化することができる。式中、#ｕ*＝（ｕ_０*，ｕ_１*，…，ｕ_ｉ*，…，ｕ_ｍ*）は、各分点における値ｕ_ｉ*を並べたベクトルである。また、ψ_ｂは、流体を構成する任意の粒子ｂと、パネル要素との距離から得られるポテンシャル関数である。また、上付添え字「*」は、時間に関する中間値を表している。なお、式（２９）の第２項は、張力係数γが乗算されたゴム膜の内部力を示している。

図１８は、ゴム膜を分割して分点（粒子）を定義し、これに基づき分点（粒子）を接続するパネル要素を生成する様子を示す図である。図１８に示すように、高さがｕ（ｘ，ｔ）で表されるゴム膜は、例えばｄｘ幅を持つようにｍ等分され、分割された切片の中心座標が仮想的な粒子の座標（ｕ_１、ｕ_２、…ｕ_ｉ、ｕ_ｉ＋１、…）と認識される。そして、分割認識された仮想的な粒子を連結する各線分が、接続要素である各パネル要素ｉとして認識及び処理される。
点（ｉ×ｄｘ，ｕ_ｉ*）と点（（ｉ＋１）×ｄｘ，ｕ_ｉ＋１*）で作られる各線分要素をｉ番目のパネル要素とすると、ゴム膜は、ｍ−１個のパネル要素で構成される。

として実施例１の手順でポテンシャル関数ψ_ｂを定義すると、ψ_ｂは#ｒ_ｂ*と#ｕ*のみの関数となるため、ポテンシャル関数ψ_ｂをψ_ｂ（#ｒ_ｂ*，#ｕ*）と表記した。

ここで、ポテンシャル関数ψ_ｂの微分∂ψ_ｂ（#ｒ_ｂ*，#ｕ*）／∂ｕ_ｉ*を解析的に求めるのは困難なため、本実施例では、次式（３１）に示す差分近似により求める。式中、ｕ_ｉ，＋*＝（ｕ_０*，ｕ_１*，…，ｕ_ｉ*＋δ，…，ｕ_ｍ*）、ｕ_ｉ，−*＝（ｕ_０*，ｕ_１*，…，ｕ_ｉ*−δ，…，ｕ_ｍ*）である。また、δは、例えば１０^−８等の定数である。次式（３１）は、パネル要素が流体から受ける反発力を示している。

一方、流体の方程式（２４）〜（２６）は、ＳＰＨ法に基づき、次式（３２）〜（３５）により離散化される。式中、βは定数である。また、Ｌ_２（#ｒ_ａ*）は二次元の再標準化行列であり、次式（３６）で表される。

ここでは、

が成立する。ｐ_ａｂ ^{ｎ＋１／２}、ｖ_ａｂ ^{ｎ＋１／２}は、流体を構成する任意の粒子ａｂ間での一次元リーマン問題を解いて求めた時空間の中間値である。ｐ_ａｂ ^{ｎ＋１／２}、ｖ_ａｂ ^{ｎ＋１／２}は、具体的には、以下のように決定される。

まず流体に含まれる粒子ａ、ｂに関して、次式（３７）〜（４０）で表される特性関数を定める。

更に、勾配を次式（４１）〜（４６）により計算する。式中、速度の上付き添え字は、各座標の成分を表している。

これらを用いて、ｐ_ａｂ ^{ｎ＋１／２}、ｖ_ａｂ ^{ｎ＋１／２}を次式（４７）〜（５１）により決定することができ、上式（３３）に示すポテンシャル関数ψａの微分を、次式（５２）の差分近似によって求めることができる。

ここでは、

が成立する。次式（５２）は、流体がパネル要素から受ける反発力を示している。

以上のような計算により、第１実施例で説明した運動解析原理を、二次元の流体−構造連成問題に適用することができる。また、同様に、第２実施例で説明した運動解析原理についても、流体や構造体としてデータが与えられたモデルに対して適用することが可能である。

図１９は、本実施例の物体運動解析装置３により実行される処理の流れを示すフローチャートである。なお、本フローチャートは、構造物と流体を処理対象とするが、構造物同士、流体同士等についても適用可能である。

まず、パラメータ入力部５０が、構造物（例えばゴム膜）及び流体についてのデータを取得する（Ｓ２００）。パラメータ入力部５０により取得されるデータは、構造物の座標、速度、質量、及び結合関係（張力係数等を含む）、並びに流体を構成する粒子の座標、速度、質量等である。

次に、粒子点設定部５１が、上記のように構造物を分割して仮想的な複数の粒子を設定する（Ｓ２０２）。本ステップの処理の概要は、図１８に則して説明した通りである。このように、粒子点設定部５１は、構造物の特に外縁（境界線）を所望の間隔で分割し、仮想的な粒子を設定する。これによって、入力データが粒子として与えられていない場合であっても、パネル要素生成部５２にパネル要素を生成させることができる。この結果、入力データが粒子データとして与えられていない場合でも、第１実施例や第２実施例と同様に、物体が境界を突き抜けるような非物理的な挙動を防止しつつ迅速に結果を出力することができる。

次に、パネル要素生成部５２が、粒子点設定部５１が設定した仮想的な粒子を接続してパネル要素を生成する（Ｓ２０４）。

次に、運動解析部５６が、Ｓ２００において入力された座標及び速度、並びにＳ２０２及びＳ２０４において生成されたパネル要素のデータを用いて、時刻をｄｔ／２進め、構造物及び流体の粒子及びパネル要素等の座標を更新する（Ｓ２０６）。

次に、反発力算出部５４が、Ｓ２０８〜Ｓ２１０の処理を、流体を構成する各粒子（以下、流体粒子と称する）について実行する。

まず、流体粒子を順に一つ取り出し、現時点における各パネル要素との位置関係に基づきポテンシャル関数ψを算出する（Ｓ２０８）。そして、ポテンシャル関数の微分を計算して相互作用項（流体粒子とパネル要素との間に作用する反発力）を算出する（Ｓ２１０）。

運動解析部５６は、各流体粒子とパネル要素との間に作用する反発力がＳ２０８〜Ｓ２１０で算出されると、算出された反発力を構造物の物理モデルに適用して構造物の加速度を算出し（Ｓ２１２）、算出した加速度を用いて構造物の速度を更新する（Ｓ２１４）。

次に、運動解析部５６は、上記算出された反発力を流体の物理モデルに適用して流体粒子の加速度を算出し（Ｓ２１６）、算出した加速度を用いて流体粒子の速度を更新する（Ｓ２１８）。

そして、運動解析部５６は、上記更新した構造物の速度、及び流体粒子の速度を用いて、時刻をｄｔ進め、構造物及び流体粒子の座標を更新する（Ｓ２２０）。

更に、運動解析部５６は、上記更新した構造物及び流体粒子の座標を、ディスプレイ装置２４等を用いて出力し（Ｓ２２２）、Ｓ２０８に戻る。なお、物体運動解析装置３は、Ｓ２０８〜Ｓ２２２のループを所定回数繰り返すと本フローを終了する等と、予め設定されていてもよい。

図２０は、本適用例において説明した手法を用いて実際にシミュレーションを行った際の各過程の状態を示す図である。物体運動解析装置１に与えたパラメータは、Ｌ＝0.292［ｍ］、ｃ＝100［ｍ／ｓ］、Ｍ＝80、ｄｘ＝2Ｌ／Ｍ［ｍ］、ｄｔ＝ｄｘ／（5ｃ）［ｓ］、σ＝10［ｋｇ／ｍ］、α＝3［ｋｇｍ］、γ＝250［Ｎ／ｓ］である。

また、係るシミュレーションでは、流体の初期密度はρ_ａ ^０＝1000［ｋｇ／ｍ^２］で一定にし、流体の粒子の質量は全てρ_ａ ^０／（ｄｘ）^２とし、ｈ＝2ｄｘとし、粒子数を2054個とした。

図示するように、ゴム膜Ｅが大きく伸びた場合でも、流体Ｄがゴム膜Ｅを突き抜けるといった非物理的挙動が起きずに計算を行うことが出来ることが示された。

以上説明した本実施例の物体運動解析装置、物体運動解析方法、及び物体運動解析プログラムによれば、非物理的な挙動を防止しつつ迅速に結果を出力することができる。

また、構造物を分割して仮想的な複数の粒子を設定する粒子点設定部５１を備えることにより、入力データが粒子として与えられていない場合であっても、パネル要素生成部５２にパネル要素を生成させることができる。この結果、入力データが粒子として与えられていない場合でも、第１実施例や第２実施例と同様に、物体が境界を突き抜けるような非物理的な挙動を防止しつつ迅速に結果を出力することができる。

以上、本発明を実施するための最良の形態について実施例を用いて説明したが、本発明はこうした実施例に何等限定されるものではなく、本発明の要旨を逸脱しない範囲内において種々の変形及び置換を加えることができる。

例えば、第１〜第３実施例において、粒子とパネル要素の最短距離のうち最も値の小さい距離を抽出し、これに基づきポテンシャル関数ψを定義するものとした。すなわち、粒子に最も近いパネル要素のみをポテンシャル関数ψの定義に反映させているが、これに限らず、ある粒子と全てのパネル要素についてポテンシャル関数ψを定義して反発力を算出してもよい。

本発明は、種々の構造物や弾性体、流体の製造産業、製造支援産業、並びにコンピュータソフトウエア産業等に広く利用可能である。

１、２、３ 物体運動解析装置
３０、４０、５０ パラメータ入力部
３２、４２、５２ パネル要素生成部
３４、４４、５４ 反発力算出部
３６、４６、５６ 運動解析部
５１ 粒子点設定部

Claims (15)

1. 第１の物体を表現した複数の粒子を含む第１の粒子群の物理量に関するデータに基づき、前記第１の粒子群に含まれる前記複数の粒子の間を接続する接続要素を生成する生成部と、
   第２の物体に含まれる第２の粒子の物理量に関するデータを参照し、前記第２の粒子と前記生成された接続要素との最短距離を算出して、前記第１の物体と前記第２の粒子との間に作用する反発力を算出する算出部と、
   前記算出された反発力を用いて、前記第１の物体及び前記第２の物体の運動を解析する解析部と、
   を有することを特徴とする物体運動解析装置。
2. 前記物体運動解析装置において、
   前記算出部は、前記粒子と前記接続要素との最短距離が減少してゼロに近づくにつれて前記反発力が大きくなり、前記第２の粒子と前記接続要素との最短距離が増加して所定距離に近づくにつれて前記反発力がゼロに近づくように、前記反発力を算出することを特徴とする、
   請求項１記載の物体運動解析装置。
3. 前記物体運動解析装置において、
   構造物を表現するデータを分割することにより、前記複数の粒子に関するデータを生成する分割部を備えることを特徴とする、
   請求項１記載の物体運動解析装置。
4. 前記物体運動解析装置において、
   前記第１の物体及び前記第２の物体は、同一平面上に存在する物体であり、
   前記接続要素は、前記複数の第１の粒子間を接続する線分要素であることを特徴とする、
   請求項１記載の物体運動解析装置。
5. 前記物体運動解析装置において、
   前記第１の物体及び前記第２の物体は、三次元空間に存在する物体であり、
   前記接続要素は、前記複数の第１の粒子から選択された三個の第１の粒子を含む閉領域であることを特徴とする、
   請求項１記載の物体運動解析装置。
6. 物体の運動を解析する物体運動解析方法において、
   コンピュータが、
   第１の物体に含まれる複数の第１の粒子の物理量に関するデータに基づき、前記複数の第１の粒子間を接続する接続要素を生成し、
   第２の物体に含まれる第２の粒子の物理量に関するデータを参照し、前記第２の粒子と前記生成された接続要素との最短距離を算出して、前記第１の物体と前記第２の粒子との間に作用する反発力を算出し、
   前記算出された反発力を用いて、前記第１の物体及び前記第２の物体の運動を解析する、
   ことを特徴とする物体運動解析方法。
7. 前記物体運動解析方法において、
   前記最短距離が減少してゼロに近づくと前記反発力が大きくなり、前記第２の粒子と前記接続要素との最短距離が増加して所定距離に近づくと前記反発力がゼロに近づくように、前記反発力を算出することを特徴とする、
   請求項６記載の物体運動解析方法。
8. 前記物体運動解析方法はさらに、
   前記コンピュータが、
   構造物を表現するデータを分割することにより、前記複数の粒子に関するデータを生成することを特徴とする、
   請求項６記載の物体運動解析方法。
9. 前記物体運動解析方法において、
   前記第１の物体及び前記第２の物体は、同一平面上に存在する物体であり、
   前記接続要素は、前記複数の第１の粒子間を接続する線分要素であることを特徴とする、
   請求項６記載の物体運動解析方法。
10. 前記物体運動解析方法において、
    前記第１の物体及び前記第２の物体は、三次元空間に存在する物体であり、
    前記接続要素は、前記複数の第１の粒子から選択された三個の第１の粒子を含む閉領域であることを特徴とする、
    請求項６記載の物体運動解析方法。
11. 物体の運動を解析する物体運動解析プログラムにおいて、
    コンピュータに、
    第１の物体に含まれる複数の第１の粒子の物理量に関するデータに基づき、前記複数の第１の粒子間を接続する接続要素を生成させ、
    第２の物体に含まれる第２の粒子の物理量に関するデータを参照し、前記第２の粒子と前記生成された接続要素との最短距離を算出して、前記第１の物体と前記第２の粒子との間に作用する反発力を算出させ、
    前記算出された反発力を用いて、前記第１の物体及び前記第２の物体の運動を解析させる、
    ことを特徴とする物体運動解析プログラム。
12. 前記物体運動解析プログラムにおいて、
    前記最短距離が減少してゼロに近づくと前記反発力が大きくなり、前記第２の粒子と前記接続要素との最短距離が増加して所定距離に近づくと前記反発力がゼロに近づくように、前記反発力を算出させることを特徴とする、
    請求項１１記載の物体運動解析プログラム。
13. 前記物体運動解析プログラムはさらに、
    前記コンピュータに、
    構造物を表現するデータを分割することにより、前記複数の粒子に関するデータを生成させることを特徴とする、
    請求項１１記載の物体運動解析プログラム。
14. 前記物体運動解析プログラムにおいて、
    前記第１の物体及び前記第２の物体は、同一平面上に存在する物体であり、
    前記接続要素は、前記複数の第１の粒子間を接続する線分要素であることを特徴とする、
    請求項１１記載の物体運動解析プログラム。
15. 前記物体運動解析プログラムにおいて、
    前記第１の物体及び前記第２の物体は、三次元空間に存在する物体であり、
    前記接続要素は、前記複数の第１の粒子から選択された三個の第１の粒子を含む閉領域であることを特徴とする、
    請求項１１記載の物体運動解析プログラム。

Images