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

Abstract

物体を表現する複数の粒子を表すデータを記憶する記憶部と、複数の粒子に含まれる着目粒子を中心とした所定範囲内に存在する他の粒子の間に所定の偏在関係が存在する場合、着目粒子を分割した複数の分割粒子を表すデータを記憶部に記憶する分割処理部と、記憶部が記憶する複数の分割粒子を表すデータを用いて、物体の運動を物体の種類に応じて解析する運動解析部と、を有することを特徴とする運動解析装置。

Description

本発明は、運動解析装置、運動解析方法及び運動解析プログラムに関する。

従来、コンピュータ上で物体の運動を再現したシミュレーション結果を出力する装置が、幅広い分野で用いられている。これに関連し、物体を粒子の集合として認識し、各粒子の運動を解析することにより粒子が集合した物体の運動を解析するための粒子法という手法が知られている。

代表的な粒子法としては、ＳＰＨ（Smoothed Particle Hydrodynamics）法やＭＰＳ（Moving Particle Semi-implicit）法等が知られている（例えば、非特許文献１、２参照）。これらの手法では、ある粒子から影響半径と呼ばれる距離ｈ以内に存在する粒子を近傍粒子として扱い、影響半径ｈ以内の粒子間には反発力が作用するものとして粒子の運動を解析する。

次式（１）は、ＳＰＨ法における運動方程式を離散化した式の一例である。式中、下付添え字ａ、ｂは、粒子を識別するための識別子である。また、#ｒ_ａ、#ｖ_ａ、ρ_ａ、Ｐ_ａ、ｍ_ａ、П_ａｂはそれぞれ粒子ａの位置ベクトル、速度ベクトル、密度、圧力、質量、粘性応力テンソルである。なお、「#」は、後に続くアルファベットがベクトルであることを示している。

上式（１）において、Ｗはカーネル関数である。カーネル関数Ｗは、粒子の分布から連続場を構成するのに用いられる。カーネル関数Ｗは、具体的には、例えば次式（２）で示す三次のスプライン関数が用いられる。式中、ｈは粒子間の影響半径であり、初期状態の平均粒子間隔の２倍から３倍程度が良く用いられる。βはカーネル関数の全空間積分量が１になるように調整された値であり、二次元の場合は０．７πｈ^２、三次元の場合はπｈ^３と決められる。

また、上式（１）は、流体の粘性を無視して完全流体と扱うことができるような場合には、粘性応力テンソルП_ａｂの項を無視することにより、次式（３）のように簡易に表すこともできる。

粒子法を用いた流動解析を行う方法等についての発明が開示されている（例えば、特許文献１参照）。この方法では、計算機毎に異なる粒子径の粒子を取り扱っており、各計算機は、担当する解析領域内に存在する粒子の粒子径が上限値を上回っている場合に、当該粒子を複数個の粒子に分割する。

特開２００７−１４０９９３号公報

J.J.Monaghan, "Smoothed Particle Hydrodynamics", Annu. Rev. Astron. Astrophys. 30: 543-74（1992） 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）

ところで、粒子法を用いて流体や弾性体の運動を解析する際には、流体や弾性体が破砕やせん断等により大きく変形すると、物体を構成する粒子が移動して疎になる箇所が部分的に出現し得る。この結果、粒子の周囲に十分な数の他の粒子が存在せず、粒子に作用すべき力（反発力やせん断力等）が正確に算出されないため、運動解析結果が現実のものから解離してしまう場合がある。

１つの側面では、本発明は、物体の運動の解析の正確性の向上を図る目的とする。

上記目的を達成するための本発明の一態様は、
物体を表現する複数の粒子を表すデータを記憶する記憶部と、
前記複数の粒子に含まれる着目粒子を中心とした所定範囲内に存在する他の粒子の間に所定の偏在関係が存在する場合、前記着目粒子を分割した複数の分割粒子を表すデータを前記記憶部に記憶する分割処理部と、
前記記憶部が記憶する前記複数の分割粒子を表すデータを用いて、前記物体の運動を前記物体の種類に応じて解析する運動解析部と、
を有することを特徴とする運動解析装置である。

１実施態様によれば、物体の運動の解析の正確性の向上を図ることができる。

本発明の第１実施例に係る運動解析装置１のハードウエア構成例である。 本発明の一実施例に係る運動解析装置１の機能構成例である。 粒子ａの周囲に存在する他の粒子ｂが所定個数未満である状況を示す図である。 粒子ａと他の粒子ｂが略同一平面上に存在する状況を示す図である。 粒子ａと他の粒子ｂが略同一直線上に存在する状況を示す図である。 略同一平面上に存在すると定義される粒子ａと粒子ｂの位置関係を例示した図である。 略同一平面上に存在すると定義される粒子ａと粒子ｂの位置関係を例示した図である。 略同一直線上に存在すると定義される粒子ａと粒子ｂの位置関係を例示した図である。 略同一直線上に存在すると定義される粒子ａと粒子ｂの位置関係を例示した図である。 粒子ａを６個の粒子ａ＿１〜ａ＿６に分割する様子を示す図である。 粒子ａを２個の粒子ａ＿１、ａ＿２に分割する様子を示す図である。 粒子ａを４個の粒子ａ＿１〜ａ＿４に分割する様子を示す図である。 カーネル関数Ｗによって表現される密度ρ*を概念的に示す図である。 運動解析装置１によって実行される全体的な処理の流れを示すフローチャートである。 分割処理が必要か否かを判定する際の処理の流れを示すフローチャートである。 分割処理の流れを示すフローチャートである。

以下、本発明を実施するための形態について、添付図面を参照しながら実施例を挙げて説明する。

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

以下、本発明の一実施例に係る運動解析装置１及び運動解析方法、並びに運動解析プログラムについて説明する。運動解析装置１は、仮想的な三次元空間又は二次元平面上で、流体や弾性体、剛体等の物体の運動を粒子法に基づき解析する装置である。

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

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

ドライブ装置１２は、記憶媒体１４からプログラムやデータを読み込み可能な装置である。プログラムを記憶した記憶媒体１４がドライブ装置１２に装着されると、プログラムが記憶媒体１４からドライブ装置１２を介して補助記憶装置１６にインストールされる。記憶媒体１４は、例えば、ＣＤ−ＲＯＭ、ＤＶＤディスク、ＵＳＢメモリ等の可搬型の記憶媒体である。また、補助記憶装置１６は、例えば、ＨＤＤ（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）等のハードウエア手段であっても構わない。

〔データ入力〕
データ入力部３０は、物体を構成する粒子に関するデータの入力を受け付け、メモリ装置１８に格納する。入力データは、流体や弾性体等をモデル化した粒子の集合データであり、具体的には、粒子の中心座標、速度、影響半径、密度、質量、変形勾配テンソル、材料物性、温度等を含む。

これらのデータは、予め補助記憶装置１６に格納されていたものであってもよいし、ユーザが入力装置２２を用いて入力したものであってもよい。また、これらのデータは、記憶媒体１４から補助記憶装置１６にインストールされたものであってもよいし、インターフェース装置２０によりネットワークから取得されたものであってもよい。

ここで、変形勾配テンソルＦ_ｉｊは、一般的には次式（４）で表される。式中、Ｘ_ｊは粒子の位置間の距離のｊ成分、ｕ_ｉは粒子の変位のｉ成分、Ｉは恒等テンソルである。ｉ成分及びｊ成分は、三次元であれば三次元座標を表すｘ、ｙ、ｚの３成分を表している。

粒子法における粒子ａの変形勾配テンソルは、粒子ａを中心として粒子ａとの距離が粒子ａの影響半径以内の領域に存在する粒子ｂからの影響を重ね合わせて、例えば次式（５）により表される。式中、ｄｕは粒子ａの変位を表している。以下、「粒子ａを中心として粒子ａとの距離が粒子ａの影響半径以内の領域に存在する」ことを、単に「粒子ａの周囲に存在する」と表記する。

なお、入力データが粒子を用いて表現されていない場合に備え、データ入力部３０は、物体のデータから粒子を用いたデータを設定する機能を有してもよい。

分割処理部３２は、メモリ装置１８に格納された各粒子の位置、影響半径等の各データを考慮し、粒子の分割の判定及び分割処理を行う。

〔分割の判定〕
粒子の分割処理は、ある粒子ａが、当該粒子ａの周囲に存在する粒子から、本来受けるべき反発力が正確に作用しなくなる状況で行われる。

係る状況は、粒子ａと他の粒子ｂの間に所定の偏在関係が存在する状況である。「所定の偏在関係」とは、例えば（Ａ）粒子ａの周囲に存在する他の粒子ｂが所定個数（例えば三次元ならば３個程度、二次元ならば２個程度）未満である状況、又は（Ｂ）粒子ａと粒子ａの周囲に存在する全ての他の粒子ｂが、略同一平面上又は略同一直線上に存在する関係である。

図３は、上記（Ａ）の状況を示す図であり（粒子ｂが２個）、図４、５は、上記（Ｂ）の状況を示す図である。図３中、ｆは粒子ａが粒子ｂから受ける反発力を示している。また、図４において平面をＳと、図５において直線をＬとそれぞれ表記した。なお、運動解析装置１の扱う空間が二次元の場合、粒子ａと他の粒子ｂが略同一平面上に存在する状況は「所定の偏在関係」から除かれる。

ここで、「略同一平面上に存在する」とは、例えば「粒子ａと、粒子ａの周辺に存在する任意の２個の粒子ｂを含む平面を基準平面Ｓとし、全ての粒子ｂについて基準平面Ｓからの距離Ｄ１が所定値Ｄｒｅｆ１以内である」ことと定義することができる。また、これに限らず、「略同一平面上に存在する」ことを、「粒子ａから全ての粒子ｂへのベクトルと基準平面Ｓのなす角度θ１が所定角度θｒｅｆ１以内である」ことと定義してもよい。図６、７は、略同一平面上に存在すると定義される粒子ａと粒子ｂの位置関係を例示した図である。

また、「略同一直線上に存在する」とは、例えば「粒子ａと、粒子ａの周辺に存在する任意の１個の粒子ｂを含む直線を基準直線Ｌとし、全ての粒子ｂについて基準直線Ｌからの距離Ｄ２が所定値Ｄｒｅｆ２以内である」ことと定義することができる。また、これに限らず、「略同一直線上に存在する」ことを、「粒子ａから全ての粒子ｂへのベクトルと基準直線Ｌのなす角度θ２が所定角度θｒｅｆ２以内である」ことと定義してもよい。図８、９は、略同一直線上に存在すると定義される粒子ａと粒子ｂの位置関係を例示した図である。

上記（Ａ）の状況では、粒子ａは周囲の粒子ｂから十分な反発力を受けないため、比較的自由に移動できる。また、上記（Ｂ）の状況では、粒子ａは基準平面Ｓ又は基準直線Ｌに直交する方向に関して比較的自由に移動できる。いずれの状況も流体や弾性体を構成する粒子の現実のふるまいとしては適切でないため、分割処理部３２は、粒子ａについて上記（Ａ）、（Ｂ）の状況が生じると、粒子ａを複数の粒子に分割して分割後の新たな粒子の位置及び質量等のデータをメモリ装置１８に格納する。

分割処理部３２は、例えば次式（６）で表される再標準化行列Ｌ（#ｒ）を計算し、行列式ｄｅｔＬ（#ｒ）を求め、行列式ｄｅｔＬ（#ｒ）を閾値Ｒｅｆと比較することにより、粒子の分割の判定を行う。具体的には、分割処理部３２は、行列式ｄｅｔＬ（#ｒ）が閾値Ｒｅｆを超える場合には粒子の分割が必要であると判定し、行列式ｄｅｔＬ（#ｒ）が閾値Ｒｅｆ以下である場合には粒子の分割が不要であると判定する。なお、次式（６）は、運動解析装置１の扱う空間が三次元の場合の再標準化行列Ｌ（#ｒ）を示している。#ｒは粒子ａの位置ベクトルであり、（ｘ，ｙ，ｚ）を成分とする。

運動解析装置１の扱う空間が二次元の場合、再標準化行列Ｌは２×２の行列となり、次式（７）で表される。この場合、#ｒは（ｘ，ｙ）を成分とする。

上式（６）、（７）のいずれの場合にも、粒子ａの周囲に存在する粒子ｂが均等に分布しているならば、ｄｅｔＬ（#ｒ）は１に近い値となる。

一方、前述した（Ａ）、（Ｂ）の状況下では、ｄｅｔＬ（#ｒ）は増加し、無限大までの値をとり得ることになる。従って、分割処理部３２は、行列式ｄｅｔＬ（#ｒ）が閾値Ｒｅｆ以上である場合に、粒子ａの分割処理を行う。閾値Ｒｅｆは、上記のような行列式ｄｅｔＬ（#ｒ）の性質から、１を超え、２を超えない程度の範囲内で任意に設定される。

なお、閾値Ｒｅｆは、「略同一平面上」や「略同一直線上」をどの程度まで含めるかを決定付けるパラメータである。閾値が大きい値であればあるほど、分割処理を行う可能性がより低くなるため、「略同一平面上」や「略同一直線上」はより狭い範囲を含むことになる。一方、閾値が小さい値であればあるほど、分割処理を行う可能性がより高くなるため、「略同一平面上」や「略同一直線上」はより広い範囲を含むことになる。

〔分割処理〕
分割処理部３２は、上記のような手法により粒子ａの分割処理が必要であると判定すると、以下のような手順で粒子ａの分割処理を実行する。

まず、粒子ａを分割する個数ｉを決定する。個数ｉについては任意に定めてもよいが、以下のように、状況に応じた分割ルールを予め定めておくと好適である。

粒子ａの周囲に存在する他の粒子ｂが所定個数未満である場合には、例えば粒子ａをｘ軸方向、ｙ軸方向、ｚ軸方向にそれぞれ分割して例えば６個の新たな粒子とする。図１０は、粒子ａを６個の粒子ａ＿１〜ａ＿６に分割する様子を示す図である。

また、粒子ａと粒子ａの周囲に存在する全ての他の粒子ｂが略同一平面上に存在する場合には、例えば粒子ａを前述した基準平面Ｓに直交する方向に分割して例えば２個の新たな粒子とする。図１１は、粒子ａを２個の粒子ａ＿１、ａ＿２に分割する様子を示す図である。

また、粒子ａと粒子ａの周囲に存在する全ての他の粒子ｂが略同一直線上に存在する場合には、例えば粒子ａを前述した基準直線Ｌに直交する２方向に分割して例えば４個の新たな粒子とする。図１２は、粒子ａを４個の粒子ａ＿１〜ａ＿４に分割する様子を示す図である。

このように、粒子ａと粒子ａの周囲に存在する粒子ｂが、偏在関係すなわち略同一平面上や略同一直線状に存在するという関係を有する場合は、偏在方向すなわち平面に平行な方向や直線の方向に、直交する方向に粒子ａを分割すると好適である。

分割処理部３２は、以下のように、分割された粒子間の距離、及び分割された粒子の影響範囲を決定する。分割処理部３２は、分割前の粒子の位置を#ｘ*とすると、当該粒子からの距離が影響半径以下の領域内で、前述した分割ルールに従って分割後の粒子の位置#ｘ^ｉを決定する（ｉ＝１、２、…）。分割後の粒子の影響半径は、分割前の粒子の影響半径と同じものを用いてよい。

なお、粒子の分割処理は、各粒子を独立に分割させるのではなく、２つの粒子を合計３つの粒子に分割するなどしても良い。

ここで、粒子法に基づく解析の性質上、分割後の粒子群によって表される密度分布は、分割前の粒子の密度分布にできるだけ近づけることが好ましい。分割前の粒子#ｘ*の密度ρ*は、次式（８）で表される。式中、Ｍは分割前の粒子#ｘ*の質量であり、Ｗは前述のカーネル関数である。図１３は、カーネル関数Ｗによって表現される密度ρ*を概念的に示す図である。

一方、分割後の粒子群があらわす密度ρは、粒子法の重ねあわせの式により、次式（９）で表される。式中、ｍ_ｉ（ｉ＝１、２、…）は分割後の各粒子の質量であり、制約条件としてＭ＝Σ_ｉｍ_ｉを満たす必要がある。

分割処理部３２は、上記制約条件下で密度ρを密度ρ*に近づけるための質量ｍ_ｉの最適解を求めるために、ラグランジュの未定乗数法（Lagrange Multiplier）を用いる。まず、分割処理部３２は、次式（１０）で表されるエネルギーＪ（ｍ_ｉ）を定義する。

そして、分割処理部３２は、定義したエネルギーＪ（ｍ_ｉ）が極値となる次式（１１）の解ｍ_ｉを、分割後の各粒子の最適な質量とする。具体的には、分割処理部３２は、次式（１１）が任意のδｍ_ｉ、δλに対して成立し、上式（１０）がｍ_ｉ、λに対して二次以下の式であるため、式（１１）よりｍ_ｉ、λに関する線形連立方程式を得る。そして、分割処理部３２は、この線形連立方程式を解くことにより、分割前の粒子により表される密度と、分割後の粒子群により表される密度の誤差を最小にする最適値を求める。

以上説明した手法は、まず分割後の粒子群の位置を決定し、次に分割された各粒子の最適な質量を求めるものであるが、以下に説明するように、まず分割後の粒子群の質量を決定し、次に分割された各粒子の最適な位置を求めることもできる。

この場合、分割処理部３２は、例えば粒子の分割前の位置からの変位の合計値が０になるという、次式（１２）で表される制約条件を設定する。

更に、分割処理部３２は、次式（１３）で表されるエネルギーＪ（#ｘ_ｉ）を定義する。

そして、分割処理部３２は、定義したエネルギーＪ（#ｘ_ｉ）が極値となる次式（１４）の解#ｘ_ｉを、分割後の各粒子の最適な位置とする。具体的には、分割処理部３２は、次式（１４）が任意のδｘ_ｉ、δλに対して成立し、上式（１３）が#ｘ_ｉ、λに対して二次以下の式であるため、式（１４）より#ｘ_ｉ、λに関する線形連立方程式を得る。そして、分割処理部３２は、この線形連立方程式を解くことにより、分割前の粒子により表される密度と、分割後の粒子群により表される密度の誤差を最小にする最適値を求める。

なお、位置を決定してから最適な質量を求める手法と、質量を決定してから最適な位置を求める手法を例示したが、分割前後の質量を一定にしつつ分割前後の密度の変化等を小さくする手法であれば、如何なる手法を用いて分割処理を行ってもよい。

〔運動解析〕
運動解析部３４は、流体、弾性体、剛体等の物体の種類毎に予め設定されている物理モデルを適用して、粒子が構成する物体の運動を解析する。以下、運動解析部３４によって用いられる物理モデルを例示する。

例えば、運動解析部３４は、物体が流体である場合は、次式（１５）〜（１７）の非圧縮性流体の運動方程式、及びナヴィエ・ストークス方程式を離散化して求められる上式（１）の運動方程式を用いることにより、加速度や温度の時間微分を求める。式中、ｕは内部エネルギーであり、Пは粘性応力テンソルであり、κは熱伝導係数であり、Ｔは温度である。

運動解析部３４は、物体が弾性体である場合は、上式（４）により求められる変形勾配テンソルから、次式（１８）〜（２０）を解くことにより、コーシー応力σを算出する。式中、ＣはCauchy-Green変形テンソルであり、ＥはGreen-Lagrange歪テンソルに相当し、ｊは体積変化率であり、Ｑは弾性ポテンシャルである。

運動解析部３４は、上記弾性ポテンシャルＱとして、次式（２１）で表される線形バネのポテンシャルや、次式（２２）〜（２４）で表されるMooney-Rivlinポテンシャルを用いる。式中、ｋはバネ定数であり、ｔｒは行列演算のトレースである。

運動解析部３４は、応力テンソルであるコーシー応力σを算出すると、次式（２５）によって運動方程式を得る。

なお、運動解析部３４は、これらの物体の種類に応じた物理モデルに加えて、周知のニュートンの運動方程式等を適用して物体の運動を解析することができる。

［処理フロー］
以下、運動解析装置１によって実行される処理の流れをフローチャートに則して説明する。図１４は、運動解析装置１によって実行される全体的な処理の流れを示すフローチャートである。

まず、データ入力部３０が、物体を構成する粒子に関するデータの入力を受け付け、メモリ装置１８に格納する（Ｓ１００）。入力されるデータは、前述したデータに加え、上限ループ数Ｎが含まれる。上限ループ数Ｎは、Ｓ１１０〜Ｓ１７０のループ処理を何回行うかを指定する値である。ループ処理の回数は、ループ処理を行なう毎に処理ループ数ｎとしてカウントされる。処理ループ数ｎの初期値は０である。

分割処理部３２は、メモリ装置１８に格納された粒子の各データを、例えば格納順に一つ取り出し（Ｓ１１０）、分割処理が必要か否かを判定する（Ｓ１２０）。

図１５は、分割処理が必要か否かを判定する際の処理の流れを示すフローチャートである。分割処理部３２は、まず、再標準化行列Ｌ（#ｒ）を生成し（Ｓ１２１）、行列式ｄｅｔＬ（#ｒ）を算出する（Ｓ１２２）。次に、分割処理部３２は、行列式ｄｅｔＬ（#ｒ）を閾値Ｒｅｆと比較することにより、粒子の分割の判定を行う（Ｓ１２３）。分割処理部３２は、行列式ｄｅｔＬ（#ｒ）が閾値Ｒｅｆを超える場合には粒子の分割が必要であると判定し（Ｓ１２４）、行列式ｄｅｔＬ（#ｒ）が閾値Ｒｅｆ以下である場合には粒子の分割は不要であると判定する（Ｓ１２５）。

分割処理部３２は、上記Ｓ１２０〜Ｓ１２５において分割処理が必要と判定された場合は、当該粒子を分割して分割後の粒子の位置及び質量等のデータをメモリ装置１８に格納する（Ｓ１３０、Ｓ１４０）。一方、分割処理部３２は、Ｓ１２０〜Ｓ１２５において分割処理が不要と判定された場合は、分割処理を実行しない（Ｓ１３０）。

図１６は、分割処理の流れを示すフローチャートである。分割処理部３２は、分割対象となる粒子の周囲に存在する他の粒子の分布に基づき分割後の粒子の個数を決定する（Ｓ１４１）。次に、分割処理部３２は、分割後の各粒子の位置及び影響半径を決定し（Ｓ１４２）、分割後の各粒子の最適な質量を決定する（Ｓ１４３）。そして、分割処理部３２は、決定した分割後の粒子の位置及び質量等のデータをメモリ装置１８に格納する（Ｓ１４４）。

次に、分割処理部３２は、解析対象の全ての粒子についてＳ１１０〜Ｓ１４４の処理を実行したか否かを判定する（Ｓ１５０）。分割処理部３２が解析対象の全ての粒子についてＳ１１０〜Ｓ１４４の処理を実行していない場合は、Ｓ１１０に戻る。分割処理部３２が解析対象の全ての粒子についてＳ１１０〜Ｓ１４４の処理を実行すると、Ｓ１６０に進む。

運動解析部３４は、メモリ装置１８に格納された粒子の各データに対して、物体の種類に応じた物理モデルを適用して、粒子の加速度を算出する（Ｓ１６０）。

そして、運動解析部３４は、加速度を時間積分して規定時間Δｔ経過後の粒子の位置、及び速度を算出し（Ｓ１７０）、必要に応じて算出結果を外部ファイルとして出力する（Ｓ１８０）。また、運動解析部３４は、各ルーチンの算出結果をディスプレイ装置２４に表示させてもよい。

次に、運動解析部３４は、処理ループ数ｎを１インクリメントし（Ｓ１９０）、処理ループ数ｎがＮと等しいか否かを判定する（Ｓ２００）。処理ループ数ｎがＮと異なる場合はＳ１１０に戻り、処理ループ数ｎがＮと等しい場合は本フローを終了する。

［まとめ］
本実施例の運動解析装置１は、解析対象の物体である流体や弾性体が破砕等により大きく変形し、物体を構成する粒子が移動して疎になる箇所が部分的に出現したとしても、粒子を分割して粒子数を増やすことができる。部分的に疎になる箇所とは、前述のように粒子の周囲に存在する他の粒子が所定個数未満である箇所や、粒子と他の粒子の間に所定の偏在関係が存在する箇所である。

この結果、運動解析装置１は、粒子の周囲に十分な数の他の粒子が存在する環境を維持することができ、粒子に作用すべき力（反発力やせん断力等）を正確に算出することができる。従って、運動解析装置１は、より正確に物体の運動を解析することができる。

なお、流体や弾性体が大きく変形しても粒子の分布が疎になる箇所が出現しないように、多くの粒子を最初から配置することも考えられる。しかしながら、この場合、演算対象となる粒子の数が増加することにより、演算負荷や演算時間が過大となってしまう場合がある。例えば、前述のように粒子をｘ軸方向、ｙ軸方向、ｚ軸方向にそれぞれ分割して例えば６個の新たな粒子とする処理を、最初から各粒子に対して実行したとすると、粒子の分布が密となり、部分的に疎になる箇所は生じにくくなる。反面、各粒子について加速度等を演算する処理が６倍に増加するため、演算負荷や演算時間が過大となる可能性がある。また、流体や弾性体の変形は予測が困難な場合があるため、どれだけ密に粒子を配置したとしても、粒子の分布が疎になる箇所が出現しないとは言い切れない。

これに対し、本実施例の運動解析装置１は、必要な箇所について分割処理を行うため、多くの粒子を最初から配置する場合に比して、トータルの演算負荷を軽減することができる。この結果、演算速度を向上させることができる。更に、必要に応じて粒子を分割することができるため、流体や弾性体が変形し続けたとしても、粒子の分布が疎になるのを持続的に解消することができる。

以上説明した本実施例の運動解析装置１によれば、より正確に物体の運動を解析することができる。

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

例えば、分割処理部３２は、再標準化行列Ｌ（#ｒ）の行列式ｄｅｔＬ（#ｒ）を閾値Ｒｅｆと比較することにより粒子の分割の判定を行うことを例示したが、粒子の分割の判定手法は、これに限られない。分割処理部３２は、粒子の配置から、当該粒子を中心とする所定範囲内に存在する他の粒子が所定個数未満である場合や、当該粒子と他の粒子が略同一平面上又は略同一直線上に存在する場合をベクトル演算等によって判別し、粒子の分割の判定を行ってもよい。

本発明は、物体の運動を粒子法により解析することを通じて、建築業や製造業、コンピュータソフトウエア産業等に利用することができる。例えば、本発明を河川や海の水の流体解析に対して適用することで、防災の計画立案に役立てることができる。また、本発明により鋳造過程を解析することで金属製品等の製品設計に用いることが出来る。また、本発明を弾性体に対して適用することで、製品設計の際に封止ゲルの形状などを適切に決定することが出来る。また、本発明を利用したプログラム自体を販売することが出来る。

１ 運動解析装置
１０ ＣＰＵ
１２ ドライブ装置
１４ 記憶媒体
１６ 補助記憶装置
１８ メモリ装置
２０ インターフェース装置
２２ 入力装置
２４ ディスプレイ装置
３０ データ入力部
３２ 分割処理部
３４ 運動解析部

Claims (9)

1. 物体を表現する複数の粒子を表すデータを記憶する記憶部と、
   前記複数の粒子に含まれる着目粒子を中心とした所定範囲内に存在する他の粒子の間に所定の偏在関係が存在する場合、前記着目粒子を分割した複数の分割粒子を表すデータを前記記憶部に記憶する分割処理部と、
   前記記憶部が記憶する前記複数の分割粒子を表すデータを用いて、前記物体の運動を前記物体の種類に応じて解析する運動解析部と、
   を有することを特徴とする運動解析装置。
2. 前記運動解析装置において、
   前記所定の偏在関係は、
   前記着目粒子を中心とした所定範囲内に存在する他の粒子が所定個数未満の関係であることを特徴とする請求項１記載の運動解析装置。
3. 前記運動解析装置において、
   前記所定の偏在関係は、
   前記着目粒子と前記他の粒子が略同一平面上又は略同一直線上に存在する関係であることを特徴とする請求項１記載の運動解析装置。
4. 前記運動解析装置において、
   前記分割処理部は、
   前記着目粒子と前記他の粒子との分布を表す関数の値を成分とする再標準化行列の行列式の値が所定値以上である場合、前記着目粒子を分割した複数の分割粒子を表すデータを前記記憶部に記憶することを特徴とする請求項１記載の運動解析装置。
5. 前記運動解析装置において、
   前記再標準化行列は、
   

   

   であり、
   前記数式中、Ｗはカーネル関数、ｈは前記着目粒子の影響半径、ｒ（ベクトル表記）は前記着目粒子の位置ベクトル、（ｘ，ｙ，ｚ）は前記着目粒子のベクトル成分、ｒ_ｂ（ベクトル表記）は前記他の粒子の位置ベクトル、（ｘ_ｂ，ｙ_ｂ，ｚ_ｂ）は前記他の粒子のベクトル成分であることを特徴とする請求項４記載の運動解析装置。
6. 前記運動解析装置において、
   前記分割処理部は、
   前記着目粒子と前記他の粒子の間に所定の偏在関係が存在する場合、前記着目粒子と前記他の粒子の偏在方向と直交する方向に、前記着目粒子を分割することを特徴とする請求項１記載の運動解析装置。
7. 前記運動解析装置において、
   前記分割処理部は、
   前記複数の分割粒子の位置を決定した後、前記複数の分割粒子の分割前後の密度関数の差分が小さくなるように、前記複数の分割粒子の質量を決定することを特徴とする請求項１記載の運動解析装置。
8. 物体を表現する複数の粒子の運動を解析する運動解析方法において、
   コンピュータが、
   前記複数の粒子に含まれる着目粒子を中心とした所定範囲内に存在する他の粒子の間に所定の偏在関係が存在する場合、前記着目粒子を表すデータを複数の分割粒子を表すデータに分割し、
   前記複数の分割粒子を表すデータを用いて、前記物体の運動を前記物体の種類に応じて解析することを特徴とする運動解析方法。
9. 物体を表現する複数の粒子の運動を解析する運動解析プログラムにおいて、
   コンピュータに、
   前記複数の粒子に含まれる着目粒子を中心とした所定範囲内に存在する他の粒子の間に所定の偏在関係が存在する場合、前記着目粒子を表すデータを複数の分割粒子を表すデータに分割させ、
   前記分割粒子を表すデータを用いて、前記物体の運動を前記物体の種類に応じて解析させることを特徴とする運動解析プログラム。

Images