特許 6091402 解析装置および解析方法

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】6091402

(24)【登録日】2017年2月17日

(45)【発行日】2017年3月8日

(54)【発明の名称】解析装置および解析方法

(51)【国際特許分類】

   G06F 19/00 20110101AFI20170227BHJP

【ＦＩ】

   G06F19/00 110

【請求項の数】4

【全頁数】11

(21)【出願番号】特願2013-232776(P2013-232776)

(22)【出願日】2013年11月11日

(65)【公開番号】特開2015-95008(P2015-95008A)

(43)【公開日】2015年5月18日

【審査請求日】2016年1月18日

(73)【特許権者】

【識別番号】000002107

【氏名又は名称】住友重機械工業株式会社

(74)【代理人】

【識別番号】100105924

【弁理士】

【氏名又は名称】森下 賢樹

(74)【代理人】

【識別番号】100109047

【弁理士】

【氏名又は名称】村田 雄祐

(74)【代理人】

【識別番号】100109081

【弁理士】

【氏名又は名称】三木 友由

(74)【代理人】

【識別番号】100116274

【弁理士】

【氏名又は名称】富所 輝観夫

(72)【発明者】

【氏名】広瀬 良太

(72)【発明者】

【氏名】市嶋 大路

【審査官】 宮地 匡人

(56)【参考文献】

【文献】 特開２０１０−１４６３６９（ＪＰ，Ａ）

【文献】 特開平０８−１７４０３５（ＪＰ，Ａ）

【文献】 米国特許出願公開第２０１１／０３０１９２４（ＵＳ，Ａ１）

【文献】 宮内 勇馬，軸受け摩耗予測計算と薄層放射化法による予測精度検証，マツダ技報，２００８年 ４月，No.26，pp.147-152

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ｇ０６Ｆ １９／００

ＪＳＴＰｌｕｓ（ＪＤｒｅａｍＩＩＩ）

ＪＳＴ７５８０（ＪＤｒｅａｍＩＩＩ）

(57)【特許請求の範囲】

【請求項1】

潤滑剤が介在する２つの物体の相対運動を数値解析によりシミュレートする解析装置であって、
一方の物体を表す粒子系に含まれる粒子と他方の物体を表す粒子系に含まれる粒子とからなる粒子対について、２つの物体の接触に関与する粒子対を特定する粒子対特定部と、
２つの物体間の摩擦状態であって、流体潤滑状態、混合潤滑状態および境界潤滑状態の各摩擦状態に応じた摩擦係数を用いて算出される摩擦力を、前記粒子対特定部によって特定された粒子対に含まれる粒子のそれぞれに付与する摩擦力付与部と、
前記摩擦力付与部によって摩擦力が付与された粒子を含む各粒子系の各粒子の運動を支配する支配方程式に基づき各粒子の状態を更新する粒子状態更新部と、を備えることを特徴とする解析装置。

【請求項2】

前記摩擦力付与部は、潤滑剤の粘度、２つの物体間の相対速度、および、２つの物体間に作用する面圧から導出されるパラメータを用いて算出される摩擦係数を用いて算出される摩擦力を付与することを特徴とする請求項１に記載の解析装置。

【請求項3】

前記摩擦力付与部は、前記パラメータの大きさに応じてその大きさが連続的に変化する摩擦係数を用いて算出される摩擦力を付与することを特徴とする請求項２に記載の解析装置。

【請求項4】

潤滑剤が介在する２つの物体の相対運動を数値解析によりシミュレートする解析方法であって、
一方の物体を表す粒子系に含まれる粒子と他方の物体を表す粒子系に含まれる粒子とからなる粒子対について、２つの物体の接触に関与する粒子対を特定するステップと、
２つの物体間の摩擦状態であって、流体潤滑状態、混合潤滑状態および境界潤滑状態の各摩擦状態に応じた摩擦係数を用いて算出される摩擦力を、特定された粒子対に含まれる粒子のそれぞれに付与するステップと、
摩擦力が付与された粒子を含む各粒子系の各粒子の運動を支配する支配方程式に基づき各粒子の状態を更新するステップと、を含むことを特徴とする解析方法。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、解析装置および解析方法に関する。

【背景技術】

【0002】

古典力学や量子力学等を基に計算機を用いて物質科学全般の現象を探るための方法として、分子動力学法（Molecular Dynamics Method、以下ＭＤ法と称す）や、量子分子動力学法（Quantum Molecular Dynamics Method）や、ＭＤ法をマクロスケールの系を扱えるように発展させた繰り込み群分子動力学法（Renormalized Molecular Dynamics、以下ＲＭＤ法と称す）に基づくシミュレーションが知られている。

【0003】

例えば本出願人は特許文献１において、物体Ａ、Ｂを構成する粒子ａ、ｂがある場合、粒子ａ、ｂの相対速度から動摩擦を計算し、計算された動摩擦を粒子の運動方程式に取り入れるＲＭＤ法を提案している。この技術によると、動摩擦を考慮したシミュレーションを実現できる。

【先行技術文献】

【特許文献】

【0004】

【特許文献1】特開２０１０−１４６３６９号公報

【発明の概要】

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

【0005】

特許文献１に記載のシミュレーション技術では、摩擦係数を一定としている。したがって、このシミュレーション技術の適用範囲は動摩擦が支配的となっている摩擦現象のシミュレーションなどであり、比較的限られている。

【0006】

本発明はこうした課題に鑑みてなされたものであり、その目的は、摩擦を含む現象のシミュレーションの適用範囲をより広げることができる解析技術の提供にある。

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

【0007】

本発明のある態様は、解析装置に関する。この解析装置は、潤滑剤が介在する２つの物体の相対運動を数値解析によりシミュレートする解析装置であって、一方の物体を表す粒子系に含まれる粒子と他方の物体を表す粒子系に含まれる粒子とからなる粒子対について、２つの物体の接触に関与する粒子対を特定する粒子対特定部と、２つの物体間の摩擦状態であって、流体潤滑状態、混合潤滑状態および境界潤滑状態の各摩擦状態に応じた摩擦係数を用いて算出される摩擦力を、粒子対特定部によって特定された粒子対に含まれる粒子のそれぞれに付与する摩擦力付与部と、摩擦力付与部によって摩擦力が付与された粒子を含む各粒子系の各粒子の運動を支配する支配方程式に基づき各粒子の状態を更新する粒子状態更新部と、を備える。

【0008】

なお、以上の構成要素の任意の組み合わせや、本発明の構成要素や表現を装置、方法、システム、コンピュータプログラム、コンピュータプログラムを格納した記録媒体などの間で相互に置換したものもまた、本発明の態様として有効である。

【発明の効果】

【0009】

本発明によれば、摩擦を含む現象のシミュレーションの適用範囲をより広げることができる。

【図面の簡単な説明】

【0010】

【図1】ストライベック曲線を示す図である。

【図2】実施の形態に係る解析装置の機能および構成を示すブロック図である。

【図3】図２の粒子データ保持部の一例を示すデータ構造図である。

【図4】図２の解析装置における一連の処理の一例を示すフローチャートである。

【図5】検証モデルを示す模式図である。

【図6】摩擦係数と、粘度、相対速度および面圧との関係を示す図である。

【発明を実施するための形態】

【0011】

以下、各図面に示される同一または同等の構成要素、部材、処理には、同一の符号を付するものとし、適宜重複した説明は省略する。

【0012】

実施の形態に係る解析装置における解析手法の原理を説明する。
図１は、潤滑剤（流体）が介在する物体間の摩擦係数と、粘度、相対速度および面圧との関係を示したストライベック曲線を示す。図１に示すように、物体間の摩擦状態には、流体潤滑状態、混合潤滑状態および境界潤滑領域の３つの状態がある。特許文献１に記載されている技術では、潤滑剤が非常に薄くなり物体同士の接触が生じている境界潤滑状態の再現は可能であるが、流体潤滑状態、混合潤滑状態および境界潤滑領域のすべてを再現したシミュレーションを実施することはできない。また、２つの物体の接触面の状態をより忠実に再現することで３つの状態を原理的に導くことも考えられるが、そのためには２つの物体の接触面により多くの粒子を配置し、かつ、潤滑剤の粒子系を導入する必要があり、粒子総数が増大して計算負荷が増える可能性がある。そこで、本実施の形態では、特許文献１に記載されている解析手法における摩擦係数部分を粘度、相対速度および面圧の関数に変更する。より具体的には、ストライベック曲線を再現する関数に変更する。これにより、粒子数や計算負荷をそれほど増やすことなく、特に潤滑剤の粒子系を導入することなく、３つの状態を再現したシミュレーションを実現する。

【0013】

特許文献１の式（８）、式（１８）をそれぞれ式（１）、式（２）として以下に示す。

【数1】

ここでＭは、それぞれ接触している物体Ａ、Ｂを表す粒子系に属する粒子ｉ、ｊの換算質量であり、ｖ_ｒはそれらの粒子の相対速度である。μは摩擦係数であり、Ｎは垂直抗力である。式（１）、式（２）より、与えられている摩擦力ｆ_ｆｒｉｃは、

【数2】

となる。

【0014】

本解析手法では、ストライベック曲線を再現するため摩擦係数μを以下に示すように潤滑剤の粘度、２つの物体間の相対速度、および、２つの物体間に作用する面圧の関数とする。

【数3】

よって摩擦力は、

【数4】

となる。ここでηは粘度、ｖ_ｒは相対速度、Ｐは面圧であり、ａ、ｂ、ｃ、ｄ、ｆ、ｎ、ｍは係数である。

【0015】

式（４）に示すように、摩擦係数μは、粘度η、相対速度ｖ_ｒおよび面圧Ｐの大きさの応じて連続的にすなわち滑らかに変化する。なお、式（４）は、例えば階段関数などの不連続関数で表される関係であってもよい。

【0016】

式（５）によると、摩擦力ｆ_ｆｒｉｃは相対速度ｖ_ｒの向きの単位ベクトルに比例するベクトルであるから、摩擦力ｆ_ｆｒｉｃの向きは相対速度ｖ_ｒの大きさに依らず相対速度ｖ_ｒの向きである。摩擦力ｆ_ｆｒｉｃの大きさは相対速度ｖ_ｒの大きさに応じて変化する。

【0017】

図２は、実施の形態に係る解析装置１００の機能および構成を示すブロック図である。ここに示す各ブロックは、ハードウエア的には、コンピュータのＣＰＵ（central processing unit）をはじめとする素子や機械装置で実現でき、ソフトウエア的にはコンピュータプログラム等によって実現されるが、ここでは、それらの連携によって実現される機能ブロックを描いている。したがって、これらの機能ブロックはハードウエア、ソフトウエアの組合せによっていろいろなかたちで実現できることは、本明細書に触れた当業者には理解されるところである。

【0018】

本実施の形態ではＲＭＤ法に倣って粒子系を解析する場合について説明するが、ＤＥＭ（Distinct Element Method）やＳＰＨ（Smoothed Particle Hydrodynamics）やＭＰＳ（Moving Particle Semi-implicit）などの粒子法や繰り込みを行わないＭＤ法に倣って粒子系を解析する場合にも、本実施の形態に係る技術的思想を適用できることは本明細書に触れた当業者には明らかである。

【0019】

解析装置１００は入力装置１０２およびディスプレイ１０４と接続される。入力装置１０２は、解析装置１００で実行される処理に関係するユーザの入力を受けるためのキーボード、マウスなどであってもよい。入力装置１０２は、インターネットなどのネットワークやＣＤ、ＤＶＤなどの記録媒体から入力を受けるよう構成されていてもよい。

【0020】

解析装置１００は、粒子系取得部１０８と、数値演算部１２０と、表示制御部１１８と、粒子データ保持部１１４と、を備える。

【0021】

粒子系取得部１０８は、入力装置１０２を介してユーザから取得する入力情報に基づき、１、２または３次元の仮想空間内に定義されるＮ（Ｎは自然数）個の粒子からなる粒子系のデータを取得する。粒子系はＲＭＤ法を使用して繰り込まれた粒子系である。

【0022】

粒子系取得部１０８は、入力情報に基づき仮想空間内にＮ個の粒子を配置し、配置されたそれぞれの粒子に速度を付与する。すなわち、粒子系取得部１０８は、粒子系に初期位置および初期速度を付与する。粒子系取得部１０８は、配置された粒子を特定する粒子ＩＤと、その粒子の位置と、その粒子の速度と、を対応付けて粒子データ保持部１１４に登録する。

【0023】

以下、解析装置１００が、潤滑剤が介在する２つの物体Ａ、Ｂの相対運動を数値解析によりシミュレートする場合について説明する。Ｎ個の粒子からなる上記粒子系は、一方の物体Ａを表す第１粒子系と、他方の物体Ｂを表す第２粒子系と、を含む。なお、上述したように、２つの物体Ａ、Ｂの間に介在する潤滑剤の粒子系は導入しない。

【0024】

以下では粒子系の粒子は全て同質または同等なものとして設定され、かつ、ポテンシャルエネルギ関数は２体のポテンシャルであって粒子によらずに同じ形を有するものとして設定される場合について説明する。しかしながら、他の場合にも本実施の形態に係る技術的思想を適用できることは、本明細書に触れた当業者には明らかである。

【0025】

数値演算部１２０は、粒子データ保持部１１４によって保持されるデータが表す第１粒子系、第２粒子系それぞれの各粒子の運動を支配する支配方程式を数値的に演算する。特に数値演算部１２０は、離散化された粒子の運動方程式にしたがった繰り返し演算を行う。
数値演算部１２０は、接触粒子対特定部１１０と、摩擦力付与部１３２と、力演算部１２２と、粒子状態演算部１２４と、状態更新部１２６と、終了条件判定部１２８と、を含む。

【0026】

接触粒子対特定部１１０は、第１粒子系に含まれる粒子と第２粒子系に含まれる粒子とからなる粒子対について、２つの物体Ａ、Ｂの接触に関与する接触粒子対を特定する。接触粒子対特定部１１０における接触粒子対の特定は、例えば特許文献１の主に図１を参照して開示される技術を使用して実現されてもよく、または他の接触判定技術を使用して実現されてもよい。

【0027】

摩擦力付与部１３２は、２つの物体間の摩擦状態であって、流体潤滑状態、混合潤滑状態および境界潤滑状態の各摩擦状態に応じた摩擦係数を用いて算出される摩擦力を、接触粒子対特定部１１０によって特定された接触粒子対に含まれる粒子のそれぞれに付与する。

【0028】

摩擦力付与部１３２は、面圧を求めるのに必要となる２つの物体間の接触面積を計算する。続いて、摩擦力付与部１３２は、２つの物体間の摩擦状態に応じた摩擦係数を式（４）に基づき算出し、算出した摩擦係数と式（５）に基づき摩擦力を算出する。摩擦力付与部１３２は、接触粒子対に含まれる２つの粒子に付与すべき摩擦力の向きを、相対速度に沿った方向で互いに逆向きとする。摩擦力付与部１３２は、接触粒子対に含まれる各粒子について決定された付与すべき摩擦力を粒子データ保持部１１４に登録する。

【0029】

力演算部１２２は粒子データ保持部１１４によって保持される第１粒子系、第２粒子系のデータを参照し、各粒子系の各粒子について、粒子間の距離に基づきその粒子に働く摩擦力以外の力を演算する。力演算部１２２は、第１粒子系の演算対象の粒子について、その演算対象の粒子との距離が所定のカットオフ距離よりも小さな粒子（以下、近接粒子と称す）を決定する。

【0030】

力演算部１２２は、各近接粒子について、その近接粒子と演算対象の粒子との間のポテンシャルエネルギ関数およびその近接粒子と演算対象の粒子との距離に基づいて、その近接粒子が演算対象の粒子に及ぼす力を演算する。特に力演算部１２２は、その近接粒子と演算対象の粒子との距離の値におけるポテンシャルエネルギ関数のグラジエント（Gradient）の値から力を算出する。力演算部１２２は、近接粒子が演算対象の粒子に及ぼす力を全ての近接粒子について足し合わせることによって、演算対象の粒子に働く力を算出する。
力演算部１２２は、第２粒子系の粒子に働く力も同様に算出する。

【0031】

粒子状態演算部１２４は粒子データ保持部１１４に保持される第１粒子系、第２粒子系のデータを参照し、各粒子系の各粒子について、離散化された粒子の運動方程式に力演算部１２２によって演算された力および粒子データ保持部１１４によって保持される摩擦力を適用することによって粒子の位置および速度のうちの少なくともひとつを演算する。本実施の形態では、粒子状態演算部１２４は粒子の位置および速度の両方を演算する。

【0032】

粒子状態演算部１２４は、接触粒子対特定部１１０によって特定された接触粒子対に含まれる粒子については、力演算部１２２によって演算された力および粒子データ保持部１１４によって保持される摩擦力を含む離散化された粒子の運動方程式から粒子の速度を演算する。粒子状態演算部１２４は、それ以外の粒子については、力演算部１２２によって演算された力を含む離散化された粒子の運動方程式から粒子の速度を演算する。

【0033】

粒子状態演算部１２４は、接触粒子対に含まれる第１粒子系の粒子について、蛙跳び法やオイラー法などの所定の数値解析の手法に基づき所定の微小な時間刻みΔｔを使用して離散化された粒子の運動方程式に、力演算部１２２によって演算された力および粒子データ保持部１１４によって保持される摩擦力を代入することによって、粒子の速度を演算する。この演算には以前の繰り返し演算のサイクルで演算された粒子の速度が使用される。

【0034】

粒子状態演算部１２４は、演算された粒子の速度に基づいて粒子の位置を算出する。粒子状態演算部１２４は、接触粒子対に含まれる第１粒子系の粒子について、所定の数値解析の手法に基づき時間刻みΔｔを使用して離散化された粒子の位置と速度の関係式に、演算された粒子の速度を適用することによって、粒子の位置を演算する。この演算には以前の繰り返し演算のサイクルで演算された粒子の位置が使用される。

【0035】

粒子状態演算部１２４における、接触粒子対に含まれない第１粒子系の粒子についての演算は摩擦力を考慮しない点を除いて上記と同様である。また、粒子状態演算部１２４は、第２粒子系の粒子についても同様に速度、位置を演算する。

【0036】

状態更新部１２６は、粒子状態演算部１２４における演算結果に基づき各粒子系の各粒子の状態を更新する。状態更新部１２６は、粒子データ保持部１１４に保持される各粒子系の各粒子の位置および速度のそれぞれを、粒子状態演算部１２４によって演算された位置および速度で更新する。

【0037】

終了条件判定部１２８は、数値演算部１２０における繰り返し演算を終了すべきか否かを判定する。繰り返し演算を終了すべき終了条件は、例えば繰り返し演算が所定の回数行われたことや、外部から終了の指示を受け付けたことや、粒子系が定常状態に達したことである。終了条件判定部１２８は、終了条件が満たされる場合、数値演算部１２０における繰り返し演算を終了させる。終了条件判定部１２８は、終了条件が満たされない場合、処理を接触粒子対特定部１１０に戻す。すると接触粒子対特定部１１０は、状態更新部１２６によって更新された粒子の位置で再び接触粒子対を特定する。

【0038】

表示制御部１１８は、粒子データ保持部１１４に保持されるデータが表す各粒子系の各粒子の位置、速度に基づき、ディスプレイ１０４に粒子系の時間発展の様子やある時刻における状態を表示させる。この表示は、静止画または動画の形式で行われてもよい。

【0039】

図３は、粒子データ保持部１１４の一例を示すデータ構造図である。粒子データ保持部１１４は、粒子ＩＤと、粒子の位置と、粒子の速度と、粒子に摩擦力が付与されている場合はその摩擦力と、を対応付けて保持する。

【0040】

上述の実施の形態において、保持部の例は、ハードディスクやメモリである。また、本明細書の記載に基づき、各部を、図示しないＣＰＵや、インストールされたアプリケーションプログラムのモジュールや、システムプログラムのモジュールや、ハードディスクから読み出したデータの内容を一時的に記憶するメモリなどにより実現できることは本明細書に触れた当業者には理解されるところである。

【0041】

以上の構成による解析装置１００の動作を説明する。
図４は、解析装置１００における一連の処理の一例を示すフローチャートである。粒子系取得部１０８は、ＲＭＤ法に倣って繰り込まれた第１粒子系および第２粒子系を取得する（Ｓ２０２）。接触粒子対特定部１１０は、特定されていない接触粒子対の有無を判定する（Ｓ２０４）。特定されていない接触粒子対が存在すると判定された場合（Ｓ２０４のＹ）、接触粒子対特定部１１０はこれまで特定されていない接触粒子対を特定する（Ｓ２０６）。摩擦力付与部１３２は、２つの物体間の摩擦状態に応じた摩擦係数を算出し、各粒子対に付与すべき摩擦力を算出する（Ｓ２０８）。摩擦力付与部１３２は、接触粒子対に含まれる各粒子に算出した摩擦力を付与（Ｓ２１０）。その後、処理はステップＳ２０４に戻る。

【0042】

ステップＳ２０４において存在しないと判定された場合（Ｓ２０４のＮ）、力演算部１２２は、粒子間の距離から粒子に働く摩擦力以外の力を演算する（Ｓ２１２）。粒子状態演算部１２４は、粒子の運動方程式から粒子の速度と位置を演算する（Ｓ２１４）。ここで、接触粒子対に含まれる粒子の運動方程式には、摩擦力が導入される。状態更新部１２６は、粒子データ保持部１１４に保持される粒子の位置、速度を演算された位置、速度で更新する（Ｓ２１６）。終了条件判定部１２８は、終了条件が満たされるか否かを判定する（Ｓ２１８）。終了条件が満たされない場合（Ｓ２１８のＮ）、処理はＳ２０４に戻される。終了条件が満たされる場合（Ｓ２１８のＹ）、処理は終了する。

【0043】

本実施の形態に係る解析装置１００によると、２つの物体の接触面に配置する粒子が特許文献１と同程度で、かつ、潤滑剤の粒子系を導入することなく、すなわち粒子数や計算負荷をそれほど増やすことなく、流体潤滑状態、混合潤滑状態および境界潤滑状態の３つの状態を再現したシミュレーションを実現できる。その結果、摩擦が関与する現象をより正確にシミュレートすることができる。

【0044】

本発明者は、本実施の形態に係る解析装置１００を使用して以下の検証を行った。
図５は、検証モデルを示す模式図である。大きな板３０２の上に小さな板３０４を置く。小さな板３０４には上から加重Ｆｎが加えられている。また、大きな板３０２と小さな板３０４とには、粘度ηの潤滑剤が介在している。この状態で、バネ３０６およびダンパー３０８を介して小さな板３０４とつながれた物体３１０を一定速度Ｖで引っ張る。ダンパー３０８は小さな板３０４の速度に対する減衰項を与える。また、摩擦係数μがこの検証モデルにおけるストライベック曲線に合うよう式（５）の各係数を、ａ＝３．０×１０^−１、ｂ＝１．０×１０^９、ｃ＝３．８、ｄ＝７．０×１０^５、ｆ＝０．０、ｎ＝１、ｍ＝１とした。

【0045】

計算条件はＶ＝０．０１〜０．１［ｍ／ｓ］、Ｆｎ＝１００［Ｎ］、η＝９．８×１０^−４［Ｐａ・ｓ］とし、このときの摩擦力（図中のＦｆ）を測定した。つまり、物体３１０を引っ張る速度を変化させながら、角速度での摩擦力Ｆｆを測定した。

【0046】

図６は、検証モデルにおけるストライベック曲線と、本実施の形態に係る手法を使用した計算結果から得られた摩擦係数と、粘度、相対速度および面圧との関係を示す。図６より、ストライベック曲線（薄い点線）と、測定した摩擦力Ｆｆから算出した摩擦係数と、粘度、相対速度および面圧との関係（濃い点線）とが比較的一致していることが分かる。これより、流体潤滑状態、混合潤滑状態および境界潤滑状態を考慮したシミュレーションが実現できることが分かる。

【0047】

以上、実施の形態に係る解析装置１００の構成と動作について説明した。この実施の形態は例示であり、その各構成要素や各処理の組み合わせにいろいろな変形例が可能なこと、またそうした変形例も本発明の範囲にあることは当業者に理解されるところである。

【0048】

実施の形態では、数値演算部１２０において粒子の位置と速度の両方を演算する場合について説明したが、これに限られない。例えば、数値解析の手法にはＶｅｒｌｅｔ法のように、粒子の位置を演算する際に粒子に働く力から粒子の位置を直接演算し、粒子の速度は陽に計算しなくてもよい手法もあり、本実施の形態に係る技術的思想をそのような手法に適用してもよい。

【符号の説明】

【0049】

１００ 解析装置、 １０２ 入力装置、 １０４ ディスプレイ、 １０８ 粒子系取得部、 １１０ 接触粒子対特定部、 １１４ 粒子データ保持部、 １１８ 表示制御部、 １２０ 数値演算部、 １２２ 力演算部、 １２４ 粒子状態演算部、 １２６ 状態更新部、 １２８ 終了条件判定部、 １３２ 摩擦力付与部。

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】