"物体を表現する複数の粒子を表すデータを記憶する記憶部と、前記複数の粒子に含まれる着目粒子を中心とした所定範囲内に存在する他の粒子の間に所定の偏在関係が存在する場合、前記着目粒子を分割した複数の分割粒子を表すデータを前記記憶部に記憶する分割処理部と、前記記憶部が記憶する前記複数の分割粒子を表すデータを用いて、前記物体の運動を前記物体の種類に応じて解析する運動解析部と、を有することを特徴とする運動解析装置。前記運動解析装置において、前記所定の偏在関係は、前記着目粒子を中心とした所定範囲内に存在する他の粒子が所定個数未満の関係であることを特徴とする請求項1記載の運動解析装置。前記運動解析装置において、前記所定の偏在関係は、前記着目粒子と前記他の粒子が略同一平面上又は略同一直線上に存在する関係であることを特徴とする請求項1記載の運動解析装置。前記運動解析装置において、前記分割処理部は、前記着目粒子と前記他の粒子との分布を表す関数の値を成分とする再標準化行列の行列式の値が所定値以上である場合、前記着目粒子を分割した複数の分割粒子を表すデータを前記記憶部に記憶することを特徴とする請求項1記載の運動解析装置。前記運動解析装置において、前記再標準化行列は、であり、前記数式中、Wはカーネル関数、hは前記着目粒子の影響半径、rベクトル表記は前記着目粒子の位置ベクトル、x,y,zは前記着目粒子のベクトル成分、rbベクトル表記は前記他の粒子の位置ベクトル、xb,yb,zbは前記他の粒子のベクトル成分であることを特徴とする請求項4記載の運動解析装置。前記運動解析装置において、前記分割処理部は、前記着目粒子と前記他の粒子の間に所定の偏在関係が存在する場合、前記着目粒子と前記他の粒子の偏在方向と直交する方向に、前記着目粒子を分割することを特徴とする請求項1記載の運動解析装置。前記運動解析装置において、前記分割処理部は、前記複数の分割粒子の位置を決定した後、前記複数の分割粒子の分割前後の密度関数の差分が小さくなるように、前記複数の分割粒子の質量を決定することを特徴とする請求項1記載の運動解析装置。物体を表現する複数の粒子の運動を解析する運動解析方法において、コンピュータが、前記複数の粒子に含まれる着目粒子を中心とした所定範囲内に存在する他の粒子の間に所定の偏在関係が存在する場合、前記着目粒子を表すデータを複数の分割粒子を表すデータに分割し、前記複数の分割粒子を表すデータを用いて、前記物体の運動を前記物体の種類に応じて解析することを特徴とする運動解析方法。物体を表現する複数の粒子の運動を解析する運動解析プログラムにおいて、コンピュータに、前記複数の粒子に含まれる着目粒子を中心とした所定範囲内に存在する他の粒子の間に所定の偏在関係が存在する場合、前記着目粒子を表すデータを複数の分割粒子を表すデータに分割させ、前記分割粒子を表すデータを用いて、前記物体の運動を前記物体の種類に応じて解析させることを特徴とする運動解析プログラム。ところで、粒子法を用いて流体や弾性体の運動を解析する際には、流体や弾性体が破砕やせん断等により大きく変形すると、物体を構成する粒子が移動して疎になる箇所が部分的に出現し得る。この結果、粒子の周囲に十分な数の他の粒子が存在せず、粒子に作用すべき力反発力やせん断力等が正確に算出されないため、運動解析結果が現実のものから解離してしまう場合がある。1つの側面では、本発明は、物体の運動の解析の正確性の向上を図る目的とする。上記目的を達成するための本発明の一態様は、物体を表現する複数の粒子を表すデータを記憶する記憶部と、前記複数の粒子に含まれる着目粒子を中心とした所定範囲内に存在する他の粒子の間に所定の偏在関係が存在する場合、前記着目粒子を分割した複数の分割粒子を表すデータを前記記憶部に記憶する分割処理部と、前記記憶部が記憶する前記複数の分割粒子を表すデータを用いて、前記物体の運動を前記物体の種類に応じて解析する運動解析部と、を有することを特徴とする運動解析装置である。1実施態様によれば、物体の運動の解析の正確性の向上を図ることができる。以下、本発明を実施するための形態について、添付図面を参照しながら実施例を挙げて説明する。なお、図面及び数式内では原則としてベクトル表記を採用するが、文章中では、が、後に続くアルファベットがベクトルであることを示すものとする。以下、本発明の一実施例に係る運動解析装置1及び運動解析方法、並びに運動解析プログラムについて説明する。運動解析装置1は、仮想的な三次元空間又は二次元平面上で、流体や弾性体、剛体等の物体の運動を粒子法に基づき解析する装置である。[ハードウエア構成]図1は、本発明の第1実施例に係る運動解析装置1のハードウエア構成例である。運動解析装置1は、例えば、CPUCentralProcessingUnit10と、ドライブ装置12と、補助記憶装置16と、メモリ装置18と、インターフェース装置20と、入力装置22と、ディスプレイ装置24と、を備える情報処理装置である。これらの構成要素は、バスやシリアル回線等を介して接続されている。CPU10は、例えば、プログラムカウンタや命令デコーダ、各種演算器、LSULoadStoreUnit、汎用レジスタ等を有するプロセッサ等の演算処理装置である。ドライブ装置12は、記憶媒体14からプログラムやデータを読み込み可能な装置である。プログラムを記憶した記憶媒体14がドライブ装置12に装着されると、プログラムが記憶媒体14からドライブ装置12を介して補助記憶装置16にインストールされる。記憶媒体14は、例えば、CD-ROM、DVDディスク、USBメモリ等の可搬型の記憶媒体である。また、補助記憶装置16は、例えば、HDDHardDiskDriveやフラッシュメモリである。プログラムのインストールは、上記のように記憶媒体14を用いる他、インターフェース装置20がネットワークを介して他のコンピュータよりダウンロードし、補助記憶装置16にインストールすることによって行うこともできる。ネットワークは、インターネット、LANLocalAreaNetwork、無線ネットワーク等である。また、プログラムは、情報処理装置の出荷時に、予め補助記憶装置16やROMReadOnlyMemory等に格納されていてもよい。このようにしてインストール又は予め格納されたプログラムをCPU10が実行することにより、図1に示す態様の情報処理装置が、本実施例の運動解析装置1として機能することができる。メモリ装置18は、例えば、RAMRandomAccessMemoryやEEPROMElectricallyErasableandProgrammableReadOnlyMemoryである。インターフェース装置20は、上記ネットワークとの接続等を制御する。入力装置22は、例えば、キーボード、マウス、タッチパッド、タッチパネルやマイク等である。また、ディスプレイ装置24は、例えば、LCDLiquidCrystalDisplayやCRTCathodeRayTube等の表示装置である。運動解析装置1は、ディスプレイ装置24の他、プリンタ、スピーカ等、他の種類の出力装置を備えてもよい。[機能構成]図2は、本発明の一実施例に係る運動解析装置1の機能構成例である。運動解析装置1は、データ入力部30と、分割処理部32と、運動解析部34と、を備える。これらの機能ブロックは、補助記憶装置16等に格納されたプログラムソフトウエアをCPU10が実行することにより機能する。なお、これらの機能ブロックが明確に分離したプログラムによって実現される必要はなく、サブルーチンや関数として他のプログラムから呼び出されるものであってもよい。また、機能ブロックの一部が、ICIntegratedCircuitやFPGAFieldProgrammableGateArray等のハードウエア手段であっても構わない。データ入力データ入力部30は、物体を構成する粒子に関するデータの入力を受け付け、メモリ装置18に格納する。入力データは、流体や弾性体等をモデル化した粒子の集合データであり、具体的には、粒子の中心座標、速度、影響半径、密度、質量、変形勾配テンソル、材料物性、温度等を含む。これらのデータは、予め補助記憶装置16に格納されていたものであってもよいし、ユーザが入力装置22を用いて入力したものであってもよい。また、これらのデータは、記憶媒体14から補助記憶装置16にインストールされたものであってもよいし、インターフェース装置20によりネットワークから取得されたものであってもよい。ここで、変形勾配テンソルFijは、一般的には次式4で表される。式中、Xjは粒子の位置間の距離のj成分、uiは粒子の変位のi成分、Iは恒等テンソルである。i成分及びj成分は、三次元であれば三次元座標を表すx、y、zの3成分を表している。粒子法における粒子aの変形勾配テンソルは、粒子aを中心として粒子aとの距離が粒子aの影響半径以内の領域に存在する粒子bからの影響を重ね合わせて、例えば次式5により表される。式中、duは粒子aの変位を表している。以下、粒子aを中心として粒子aとの距離が粒子aの影響半径以内の領域に存在することを、単に粒子aの周囲に存在すると表記する。なお、入力データが粒子を用いて表現されていない場合に備え、データ入力部30は、物体のデータから粒子を用いたデータを設定する機能を有してもよい。分割処理部32は、メモリ装置18に格納された各粒子の位置、影響半径等の各データを考慮し、粒子の分割の判定及び分割処理を行う。分割の判定粒子の分割処理は、ある粒子aが、当該粒子aの周囲に存在する粒子から、本来受けるべき反発力が正確に作用しなくなる状況で行われる。係る状況は、粒子aと他の粒子bの間に所定の偏在関係が存在する状況である。所定の偏在関係とは、例えばA粒子aの周囲に存在する他の粒子bが所定個数例えば三次元ならば3個程度、二次元ならば2個程度未満である状況、又はB粒子aと粒子aの周囲に存在する全ての他の粒子bが、略同一平面上又は略同一直線上に存在する関係である。図3は、上記Aの状況を示す図であり粒子bが2個、図4、5は、上記Bの状況を示す図である。図3中、fは粒子aが粒子bから受ける反発力を示している。また、図4において平面をSと、図5において直線をLとそれぞれ表記した。なお、運動解析装置1の扱う空間が二次元の場合、粒子aと他の粒子bが略同一平面上に存在する状況は所定の偏在関係から除かれる。ここで、略同一平面上に存在するとは、例えば粒子aと、粒子aの周辺に存在する任意の2個の粒子bを含む平面を基準平面Sとし、全ての粒子bについて基準平面Sからの距離D1が所定値Dref1以内であることと定義することができる。また、これに限らず、略同一平面上に存在することを、粒子aから全ての粒子bへのベクトルと基準平面Sのなす角度θ1が所定角度θref1以内であることと定義してもよい。図6、7は、略同一平面上に存在すると定義される粒子aと粒子bの位置関係を例示した図である。また、略同一直線上に存在するとは、例えば粒子aと、粒子aの周辺に存在する任意の1個の粒子bを含む直線を基準直線Lとし、全ての粒子bについて基準直線Lからの距離D2が所定値Dref2以内であることと定義することができる。また、これに限らず、略同一直線上に存在することを、粒子aから全ての粒子bへのベクトルと基準直線Lのなす角度θ2が所定角度θref2以内であることと定義してもよい。図8、9は、略同一直線上に存在すると定義される粒子aと粒子bの位置関係を例示した図である。上記Aの状況では、粒子aは周囲の粒子bから十分な反発力を受けないため、比較的自由に移動できる。また、上記Bの状況では、粒子aは基準平面S又は基準直線Lに直交する方向に関して比較的自由に移動できる。いずれの状況も流体や弾性体を構成する粒子の現実のふるまいとしては適切でないため、分割処理部32は、粒子aについて上記A、Bの状況が生じると、粒子aを複数の粒子に分割して分割後の新たな粒子の位置及び質量等のデータをメモリ装置18に格納する。分割処理部32は、例えば次式6で表される再標準化行列Lrを計算し、行列式detLrを求め、行列式detLrを閾値Refと比較することにより、粒子の分割の判定を行う。具体的には、分割処理部32は、行列式detLrが閾値Refを超える場合には粒子の分割が必要であると判定し、行列式detLrが閾値Ref以下である場合には粒子の分割が不要であると判定する。なお、次式6は、運動解析装置1の扱う空間が三次元の場合の再標準化行列Lrを示している。rは粒子aの位置ベクトルであり、x,y,zを成分とする。運動解析装置1の扱う空間が二次元の場合、再標準化行列Lは2×2の行列となり、次式7で表される。この場合、rはx,yを成分とする。上式6、7のいずれの場合にも、粒子aの周囲に存在する粒子bが均等に分布しているならば、detLrは1に近い値となる。一方、前述したA、Bの状況下では、detLrは増加し、無限大までの値をとり得ることになる。従って、分割処理部32は、行列式detLrが閾値Ref以上である場合に、粒子aの分割処理を行う。閾値Refは、上記のような行列式detLrの性質から、1を超え、2を超えない程度の範囲内で任意に設定される。なお、閾値Refは、略同一平面上や略同一直線上をどの程度まで含めるかを決定付けるパラメータである。閾値が大きい値であればあるほど、分割処理を行う可能性がより低くなるため、略同一平面上や略同一直線上はより狭い範囲を含むことになる。一方、閾値が小さい値であればあるほど、分割処理を行う可能性がより高くなるため、略同一平面上や略同一直線上はより広い範囲を含むことになる。分割処理分割処理部32は、上記のような手法により粒子aの分割処理が必要であると判定すると、以下のような手順で粒子aの分割処理を実行する。まず、粒子aを分割する個数iを決定する。個数iについては任意に定めてもよいが、以下のように、状況に応じた分割ルールを予め定めておくと好適である。粒子aの周囲に存在する他の粒子bが所定個数未満である場合には、例えば粒子aをx軸方向、y軸方向、z軸方向にそれぞれ分割して例えば6個の新たな粒子とする。図10は、粒子aを6個の粒子a_1~a_6に分割する様子を示す図である。また、粒子aと粒子aの周囲に存在する全ての他の粒子bが略同一平面上に存在する場合には、例えば粒子aを前述した基準平面Sに直交する方向に分割して例えば2個の新たな粒子とする。図11は、粒子aを2個の粒子a_1、a_2に分割する様子を示す図である。また、粒子aと粒子aの周囲に存在する全ての他の粒子bが略同一直線上に存在する場合には、例えば粒子aを前述した基準直線Lに直交する2方向に分割して例えば4個の新たな粒子とする。図12は、粒子aを4個の粒子a_1~a_4に分割する様子を示す図である。このように、粒子aと粒子aの周囲に存在する粒子bが、偏在関係すなわち略同一平面上や略同一直線状に存在するという関係を有する場合は、偏在方向すなわち平面に平行な方向や直線の方向に、直交する方向に粒子aを分割すると好適である。分割処理部32は、以下のように、分割された粒子間の距離、及び分割された粒子の影響範囲を決定する。分割処理部32は、分割前の粒子の位置をxとすると、当該粒子からの距離が影響半径以下の領域内で、前述した分割ルールに従って分割後の粒子の位置xiを決定するi=1、2、...。分割後の粒子の影響半径は、分割前の粒子の影響半径と同じものを用いてよい。なお、粒子の分割処理は、各粒子を独立に分割させるのではなく、2つの粒子を合計3つの粒子に分割するなどしても良い。ここで、粒子法に基づく解析の性質上、分割後の粒子群によって表される密度分布は、分割前の粒子の密度分布にできるだけ近づけることが好ましい。分割前の粒子xの密度ρは、次式8で表される。式中、Mは分割前の粒子xの質量であり、Wは前述のカーネル関数である。図13は、カーネル関数Wによって表現される密度ρを概念的に示す図である。一方、分割後の粒子群があらわす密度ρは、粒子法の重ねあわせの式により、次式9で表される。式中、mii=1、2、...は分割後の各粒子の質量であり、制約条件としてM=Σimiを満たす必要がある。分割処理部32は、上記制約条件下で密度ρを密度ρに近づけるための質量miの最適解を求めるために、ラグランジュの未定乗数法LagrangeMultiplierを用いる。まず、分割処理部32は、次式10で表されるエネルギーJmiを定義する。そして、分割処理部32は、定義したエネルギーJmiが極値となる次式11の解miを、分割後の各粒子の最適な質量とする。具体的には、分割処理部32は、次式11が任意のδmi、δλに対して成立し、上式10がmi、λに対して二次以下の式であるため、式11よりmi、λに関する線形連立方程式を得る。そして、分割処理部32は、この線形連立方程式を解くことにより、分割前の粒子により表される密度と、分割後の粒子群により表される密度の誤差を最小にする最適値を求める。以上説明した手法は、まず分割後の粒子群の位置を決定し、次に分割された各粒子の最適な質量を求めるものであるが、以下に説明するように、まず分割後の粒子群の質量を決定し、次に分割された各粒子の最適な位置を求めることもできる。この場合、分割処理部32は、例えば粒子の分割前の位置からの変位の合計値が0になるという、次式12で表される制約条件を設定する。更に、分割処理部32は、次式13で表されるエネルギーJxiを定義する。そして、分割処理部32は、定義したエネルギーJxiが極値となる次式14の解xiを、分割後の各粒子の最適な位置とする。具体的には、分割処理部32は、次式14が任意のδxi、δλに対して成立し、上式13がxi、λに対して二次以下の式であるため、式14よりxi、λに関する線形連立方程式を得る。そして、分割処理部32は、この線形連立方程式を解くことにより、分割前の粒子により表される密度と、分割後の粒子群により表される密度の誤差を最小にする最適値を求める。なお、位置を決定してから最適な質量を求める手法と、質量を決定してから最適な位置を求める手法を例示したが、分割前後の質量を一定にしつつ分割前後の密度の変化等を小さくする手法であれば、如何なる手法を用いて分割処理を行ってもよい。運動解析運動解析部34は、流体、弾性体、剛体等の物体の種類毎に予め設定されている物理モデルを適用して、粒子が構成する物体の運動を解析する。以下、運動解析部34によって用いられる物理モデルを例示する。例えば、運動解析部34は、物体が流体である場合は、次式15~17の非圧縮性流体の運動方程式、及びナヴィエストークス方程式を離散化して求められる上式1の運動方程式を用いることにより、加速度や温度の時間微分を求める。式中、uは内部エネルギーであり、Пは粘性応力テンソルであり、κは熱伝導係数であり、Tは温度である。運動解析部34は、物体が弾性体である場合は、上式4により求められる変形勾配テンソルから、次式18~20を解くことにより、コーシー応力σを算出する。式中、CはCauchyGreen変形テンソルであり、EはGreenLagrange歪テンソルに相当し、jは体積変化率であり、Qは弾性ポテンシャルである。運動解析部34は、上記弾性ポテンシャルQとして、次式21で表される線形バネのポテンシャルや、次式22~24で表されるMooneyRivlinポテンシャルを用いる。式中、kはバネ定数であり、trは行列演算のトレースである。運動解析部34は、応力テンソルであるコーシー応力σを算出すると、次式25によって運動方程式を得る。なお、運動解析部34は、これらの物体の種類に応じた物理モデルに加えて、周知のニュートンの運動方程式等を適用して物体の運動を解析することができる。[処理フロー]以下、運動解析装置1によって実行される処理の流れをフローチャートに則して説明する。図14は、運動解析装置1によって実行される全体的な処理の流れを示すフローチャートである。まず、データ入力部30が、物体を構成する粒子に関するデータの入力を受け付け、メモリ装置18に格納するS100。入力されるデータは、前述したデータに加え、上限ループ数Nが含まれる。上限ループ数Nは、S110~S170のループ処理を何回行うかを指定する値である。ループ処理の回数は、ループ処理を行なう毎に処理ループ数nとしてカウントされる。処理ループ数nの初期値は0である。分割処理部32は、メモリ装置18に格納された粒子の各データを、例えば格納順に一つ取り出しS110、分割処理が必要か否かを判定するS120。図15は、分割処理が必要か否かを判定する際の処理の流れを示すフローチャートである。分割処理部32は、まず、再標準化行列Lrを生成しS121、行列式detLrを算出するS122。次に、分割処理部32は、行列式detLrを閾値Refと比較することにより、粒子の分割の判定を行うS123。分割処理部32は、行列式detLrが閾値Refを超える場合には粒子の分割が必要であると判定しS124、行列式detLrが閾値Ref以下である場合には粒子の分割は不要であると判定するS125。分割処理部32は、上記S120~S125において分割処理が必要と判定された場合は、当該粒子を分割して分割後の粒子の位置及び質量等のデータをメモリ装置18に格納するS130、S140。一方、分割処理部32は、S120~S125において分割処理が不要と判定された場合は、分割処理を実行しないS130。図16は、分割処理の流れを示すフローチャートである。分割処理部32は、分割対象となる粒子の周囲に存在する他の粒子の分布に基づき分割後の粒子の個数を決定するS141。次に、分割処理部32は、分割後の各粒子の位置及び影響半径を決定しS142、分割後の各粒子の最適な質量を決定するS143。そして、分割処理部32は、決定した分割後の粒子の位置及び質量等のデータをメモリ装置18に格納するS144。次に、分割処理部32は、解析対象の全ての粒子についてS110~S144の処理を実行したか否かを判定するS150。分割処理部32が解析対象の全ての粒子についてS110~S144の処理を実行していない場合は、S110に戻る。分割処理部32が解析対象の全ての粒子についてS110~S144の処理を実行すると、S160に進む。運動解析部34は、メモリ装置18に格納された粒子の各データに対して、物体の種類に応じた物理モデルを適用して、粒子の加速度を算出するS160。そして、運動解析部34は、加速度を時間積分して規定時間Δt経過後の粒子の位置、及び速度を算出しS170、必要に応じて算出結果を外部ファイルとして出力するS180。また、運動解析部34は、各ルーチンの算出結果をディスプレイ装置24に表示させてもよい。次に、運動解析部34は、処理ループ数nを1インクリメントしS190、処理ループ数nがNと等しいか否かを判定するS200。処理ループ数nがNと異なる場合はS110に戻り、処理ループ数nがNと等しい場合は本フローを終了する。[まとめ]本実施例の運動解析装置1は、解析対象の物体である流体や弾性体が破砕等により大きく変形し、物体を構成する粒子が移動して疎になる箇所が部分的に出現したとしても、粒子を分割して粒子数を増やすことができる。部分的に疎になる箇所とは、前述のように粒子の周囲に存在する他の粒子が所定個数未満である箇所や、粒子と他の粒子の間に所定の偏在関係が存在する箇所である。この結果、運動解析装置1は、粒子の周囲に十分な数の他の粒子が存在する環境を維持することができ、粒子に作用すべき力反発力やせん断力等を正確に算出することができる。従って、運動解析装置1は、より正確に物体の運動を解析することができる。なお、流体や弾性体が大きく変形しても粒子の分布が疎になる箇所が出現しないように、多くの粒子を最初から配置することも考えられる。しかしながら、この場合、演算対象となる粒子の数が増加することにより、演算負荷や演算時間が過大となってしまう場合がある。例えば、前述のように粒子をx軸方向、y軸方向、z軸方向にそれぞれ分割して例えば6個の新たな粒子とする処理を、最初から各粒子に対して実行したとすると、粒子の分布が密となり、部分的に疎になる箇所は生じにくくなる。反面、各粒子について加速度等を演算する処理が6倍に増加するため、演算負荷や演算時間が過大となる可能性がある。また、流体や弾性体の変形は予測が困難な場合があるため、どれだけ密に粒子を配置したとしても、粒子の分布が疎になる箇所が出現しないとは言い切れない。これに対し、本実施例の運動解析装置1は、必要な箇所について分割処理を行うため、多くの粒子を最初から配置する場合に比して、トータルの演算負荷を軽減することができる。この結果、演算速度を向上させることができる。更に、必要に応じて粒子を分割することができるため、流体や弾性体が変形し続けたとしても、粒子の分布が疎になるのを持続的に解消することができる。以上説明した本実施例の運動解析装置1によれば、より正確に物体の運動を解析することができる。以上、本発明を実施するための最良の形態について実施例を用いて説明したが、本発明はこうした実施例に何等限定されるものではなく、本発明の要旨を逸脱しない範囲内において種々の変形及び置換を加えることができる。例えば、分割処理部32は、再標準化行列Lrの行列式detLrを閾値Refと比較することにより粒子の分割の判定を行うことを例示したが、粒子の分割の判定手法は、これに限られない。分割処理部32は、粒子の配置から、当該粒子を中心とする所定範囲内に存在する他の粒子が所定個数未満である場合や、当該粒子と他の粒子が略同一平面上又は略同一直線上に存在する場合をベクトル演算等によって判別し、粒子の分割の判定を行ってもよい。"