JP5241469B2 - シミュレーション方法及びプログラム - Google Patents

Description

本発明は、シミュレーション方法及びプログラムに関する。特に、物体を粒子の集まりで表現し、ある時間刻み幅で、各粒子の挙動を数値計算することにより、時間軸上で離散的に分布する時刻ごとの物体の状態を解析する粒子法におけるシミュレーション方法、及びそのシミュレーション方法をコンピュータで実行するプログラムに関する。

コンピュータを用いて物質科学全般の現象を探求する方法として、分子動力学法に基づく分子シミュレーションが知られている。分子シミュレーションによって、分子のポテンシャルエネルギや最安定構造等物質の特性を、分子レベルで解明することが可能である。

分子動力学法に基づく分子シミュレーションに関して、種々の研究がなされている。たとえば少ない計算量で分子動力学計算を実行するシミュレーション方法の発明が開示されている（たとえば、特許文献１参照）。

粒子シミュレーションの動解析を行う場合、摩擦力のモデル化が問題となる。

現在提案されている代表的な粒子法には、ＭＰＳ(Moving Particle Semi-implicit)、ＳＰＨ(Smoothed Particle Hydrodynamics)、ＤＥＭ(Distinct Element Method)などがある。

これらの手法では、粒子の接触時にせん断力を決定し、その大きさが静止摩擦力を上回らないようにすることで摩擦力をモデル化している。また摩擦力は、接触粒子の運動を妨げる方向、すなわち接触時における粒子の相対速度の方向と反対方向に発生するとしている。

しかしながら粒子は、物体全体の重心運動に関わらず、常に高周波数で振動している。そのため上記の手法においては、粒子に働く摩擦力の方向は、物体全体の重心運動とは関係なく、原子の振動周期と同じオーダーで目まぐるしく変化することになる。ゆえに、これらの手法では、滑り摩擦時に要するエネルギを正しく計算できない可能性が高い。

図６（Ａ）及び（Ｂ）を参照して、上述の内容を詳述する。

図６（Ａ）に、物体全体の摺動の様子を示す。物体は摩擦力Ｆに逆らって、ｘ軸正方向に距離Ｌだけ移動する。物体の重心は常にｘ軸正方向に移動している。この場合、物体には移動中、常に摩擦力Ｆがｘ軸負方向に作用しており、物体に作用する摩擦力の大きさの平均もＦとなる。

図６（Ｂ）に物体の接触部粒子の挙動を示す。上述の粒子法によれば、摩擦力を受けるのは接触部の粒子のみであり、粒子の振動は物体の重心の運動方向と同じであるとは限らない。たとえば接触部の原子は、本図に示す軌跡を辿って、ｘ軸正方向に距離Ｌを移動する場合がある。

本図に示す例においては、接触部の原子は最初の時間ｄｔでｘ軸正方向に距離Ｌ_１だけ移動する。続く時間ｄｔではｘ軸負方向に距離Ｌ_２だけ移動する。最後の時間ｄｔではｘ軸正方向に距離Ｌ_３だけ移動する。ここでＬ_１−Ｌ_２＋Ｌ_３＝Ｌを満たす。このとき距離Ｌを摺動する間に、接触部の原子に働く摩擦力の大きさは、最初の時間ｄｔではｘ軸負方向にＦ、次の時間ｄｔではｘ軸正方向にＦ、最後の時間ｄｔではｘ軸負方向にＦであるから、平均して働く摩擦力は、ｘ軸負方向に、（Ｆ−Ｆ＋Ｆ）／３＝Ｆ／３ となり、図６（Ａ）を参照して考察した結果とは異なってしまう。

このことは時間刻みｄｔを十分小さくしても、摺動距離を長くしても起こりうることであり、防止策としては、粒子の大きさを十分大きくとる、もしくは十分な時間的、空間的なアンサンブルをとることなどが考えられる。

しかしながら、粒子を大きくすると複雑な形状が再現できない、時間アンサンブルを十分とるとフィードバックが必要となり、計算に余分な時間がかかる、空間アンサンブルが十分とれない場合は用いることができない、などの問題点があり実用的ではない。

従来の粒子法における摩擦のモデルでは滑り摩擦によるエネルギ散逸が正しく評価できず、また、滑り摩擦と転がり摩擦とを区別して評価する必要がある。しかし現実の摩擦は、滑り摩擦と転がり摩擦とが混在して生じているため、そこに人為的な仮定を含むことは望ましくない。

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

本発明の目的は、高精度なシミュレーションを実現するシミュレーション方法を提供することである。

また、高精度なシミュレーションを実現するシミュレーション方法を実行させるプログラムを提供することである。

本発明の一観点によれば、
物体を粒子の集まりで表現し、ある時間刻み幅で、各粒子の挙動を数値計算することにより、時間軸上で離散的に分布する時刻ごとの該物体の状態を解析する粒子法におけるシミュレーション方法であって、
（ａ）第１の物体を構成する第１の粒子と、第２の物体を構成する第２の粒子との、第１の時刻における相対速度に基づいて、前記第１の粒子と前記第２の粒子との相対運動で散逸されるエネルギが、前記第１の粒子と前記第２の粒子との摩擦によって散逸されるエネルギと等しいとして、前記第１の粒子及び前記第２の粒子に印加される摩擦力を算出する工程と、
（ｂ）前記第１の粒子及び前記第２の粒子に働く力として、粒子間ポテンシャルに基づく力、外場から受ける力、及び前記工程（ａ）で算出された摩擦力を用いて、運動方程式を解くことにより、前記第１の時刻の次の第２の時刻における、前記第１の粒子、及び前記第２の粒子の位置を求める工程と
を有するシミュレーション方法が提供される。

本発明の他の観点によると、
物体を粒子の集まりで表現し、ある時間刻み幅で、各粒子の挙動を数値計算することにより、時間軸上で離散的に分布する時刻ごとの該物体の状態を解析する粒子法におけるシミュレーションを行うためにコンピュータを、
第１の物体を構成する第１の粒子、及び、第２の物体を構成する第２の粒子の初期状態を入力する手段、
前記第１の粒子と前記第２の粒子との、第１の時刻における相対速度に基づいて、前記
第１の粒子と前記第２の粒子との相対運動で散逸されるエネルギが、前記第１の粒子と前記第２の粒子との摩擦によって散逸されるエネルギと等しいとして、前記第１の粒子及び前記第２の粒子に印加される摩擦力を算出する手段、
前記第１の粒子及び前記第２の粒子に働く力として、粒子間ポテンシャルに基づく力、外場から受ける力、及び前記摩擦力を用いて、運動方程式を解くことにより、前記第１の時刻の次の第２の時刻における、前記第１の粒子、及び前記第２の粒子の位置を求める手段、
前記第２の時刻における前記第１の粒子及び前記第２の粒子を表示する表示手段
として機能させるためのプログラムが提供される。

本発明によれば、高精度なシミュレーションを実現するシミュレーション方法を提供することができる。

また、高精度なシミュレーションを実現するシミュレーション方法を実行させるプログラムを提供することができる。

以下、物体を粒子の集まりで表現し、ある時間刻み幅Δｔで、各粒子の挙動を数値計算することにより、時間軸上で離散的に分布する時刻ごとの物体の状態を解析する粒子法におけるシミュレーション方法、及びそのシミュレーション方法をコンピュータで実行するプログラムの実施例について説明する。

質量ｍ_ｉの粒子ｉと質量ｍ_ｊの粒子ｊが摩擦した場合を考える。２つの粒子についての散逸方程式は、拡散係数γ（ｔ）を用いて、下の式（１）及び（２）で表される。

・・（１）

・・（２）

ここでｖ_ｉ、ｖ_ｊは、それぞれ粒子ｉ、ｊの速度を示す。

また、粒子ｊに対する粒子ｉの相対速度をｖ_ｒ、粒子ｉと粒子ｊの重心の速度をｖ_ｇとすると、ｖ_ｒとｖ_ｇは、それぞれ下の（３）、（４）式で表される。

・・（３）

・・（４）

（３）、（４）式を、ｖ_ｉ、ｖ_ｊについて解くと、下の（５）、（６）式を得る。

・・（５）

・・（６）

式（３）、（５）を式（１）に代入して得られた式と、式（３）、（６）を式（２）に代入して得られた式とを連立することで、以下の式（７）及び（８）が導かれる。

・・（７）

・・（８）

ここで、

・・（９）

である。

式（８）を解くと、

・・（１０）

が得られる。

式（７）、（１０）より、散逸過程においては、重心の運動量は保存され、相対運動のエネルギが散逸されることがわかる。

以下、数値積分を実行する。式（８）を離散化すると、ある時刻ｎについて、

・・（１１）

を得る。

式（１１）を、ｖ_ｒ ^{［ｎ＋１／２］}及びｖ_ｒ ^{［ｎ−１／２］}についてまとめると、

・・（１２）

となり、更にこれは、

・・（１３）

と変形することができる。

一方、粒子ｉについての離散化された散逸方程式は、式（１）、（８）、（１１）より、

・・（１４）

と書ける。式（１３）を式（１４）に代入すると、

・・（１５）

を得、これを整理すると、

・・（１６）

となる。

ここで本願発明者らは、従来のように粒子に阻止力を加える摩擦モデルではなく、２粒子の相対運動で散逸されるエネルギを、摩擦によって散逸されるエネルギと関連づけて考える、たとえば両者を同等と考える摩擦モデルを導入することにした。

拡散係数γとクーロン摩擦力とを関係づけるために、２粒子の相対運動による散逸パワー（単位時間当たりの散逸エネルギ）と、摩擦による散逸パワーとを等しいとして、

・・（１７）

すなわち、

・・（１８）

という関係式を導入する。ここで、μは摩擦係数、Ｎは垂直抗力である。

式（１６）に式（１８）を代入し、代入した式の右辺に、粒子ｉに働く力のうち、粒子ｊとの間の摩擦力を除く力をｆ_ｉとして加えると、粒子ｉについての離散化された運動方程式（１９）が得られる。

・・（１９）

同様に、粒子ｊについても、粒子ｊに働く力のうち、粒子ｉとの間の摩擦力を除く力ｆ_ｊを用いて、離散化された運動方程式（２０）を得ることができる。

・・（２０）

力ｆ_ｉ、ｆ_ｊには、たとえば、重力、電磁力、原子間力等が含まれる。粒子ｉ、ｊの運動は、運動方程式（１９）、（２０）を解くことによって把握される。式（１９）、（２０）から明らかなように、本願発明者らの提案する摩擦モデルによれば、粒子ｉが粒子ｊから受ける摩擦力Ｆ_ｉは、

・・（２１）

で表され、粒子ｊが粒子ｉから受ける摩擦力Ｆ_ｊは、

・・（２２）

で表される。

図１及び図２を参照して、実施例によるシミュレーション方法について説明する。

まず、シミュレーション方法全体の説明に先立って、図１を参照し、相互間に摩擦力を生じる粒子の決定方法について説明する。

たとえば物体Ａを構成する粒子の１つである粒子ｐが、物体Ｂを構成する隣接粒子である粒子ｉ、ｊ、及びｋの中心の張る平面において、粒子ｉ、ｊ、ｋの定める三角形（三角要素）に接触したとする。ここで隣接粒子とは、多数の粒子で構成される物体Ｂをボロノイ分割したときに、隣り合う粒子をいう。また、以下粒子の位置が問題となるときには、粒子の中心の位置を粒子の位置と考える。

粒子ｉ、ｊ、ｋの定める三角形の外心をＣ、外心Ｃから線分ｉｊに下ろした垂線の足をＨ_ｉｊ、線分ｊｋに下ろした垂線の足をＨ_ｊｋ、線分ｋｉに下ろした垂線の足をＨ_ｋｉとし、四角形ｉＨ_ｉｊＣＨ_ｋｉの周及び内部から、線分ＣＨ_ｉｊ及び線分ＣＨ_ｋｉを除いた領域を、粒子ｉの支配領域と定める。同様に、四角形ｊＨ_ｊｋＣＨ_ｉｊの周及び内部から、線分ＣＨ_ｊｋ及び線分ＣＨ_ｉｊを除いた領域を粒子ｊの支配領域、四角形ｋＨ_ｋｉＣＨ_ｊｋの周及び内部から、線分ＣＨ_ｋｉ及び線分ＣＨ_ｊｋを除いた領域を粒子ｋの支配領域と定める。

また、線分ＣＨ_ｉｊ上の点は、粒子ｉの支配領域、粒子ｊの支配領域のいずれかに含める。同様に、線分ＣＨ_ｊｋ上の点は、粒子ｊ、ｋのいずれかの支配領域に含め、線分ＣＨ_ｋｉ上の点は、粒子ｋ、ｉのいずれかの支配領域に含める。

そして粒子ｐが、粒子ｉ、ｊ、ｋのいずれの粒子の支配領域に存在するかを確認し、粒子ｐが存在する領域を支配する粒子との間に摩擦力が発生すると考える。たとえば粒子ｐが粒子ｉの支配領域に存在する場合、粒子ｐと粒子ｉとの間に、式（２１）、（２２）で表される摩擦力を導入する。そして両粒子ｐ、ｉは、式（１９）、（２０）に則った相互作用をすると考える。また、粒子ｐが粒子ｉの支配領域に存在する場合、粒子ｐと粒子ｊ、ｋ間の摩擦力は０とする。

なお、上記摩擦力導入の考え方によれば、三角形ｉｊｋの周または内部に接触した粒子ｐは、粒子ｉ、ｊ、ｋのうち最も近い位置にある粒子との間に摩擦力を導入されることになる。

次に、図２に示すフローチャートを用いて、実施例によるシミュレーション方法を説明する。

実施例によるシミュレーション方法では、まずステップＳ１０１において、粒子の位置や速度など、シミュレーションを行うに当たっての初期状態を決定する。

次に、ステップＳ１０２において、粒子間のポテンシャル、及び外場が及ぼす力ｆを全ての粒子につき、粒子ごとに導出する。ｆは粒子によって異なりうる力である。

物体を構成する粒子同士は、互いに粒子間のポテンシャルにより結びついている。粒子間のポテンシャルは、たとえばバネやＬｅｎｎａｒｄ−Ｊｏｎｅｓポテンシャルで表現される。どのようなポテンシャルを用いるかについては、ステップＳ１０１またはそれに先立って、シミュレーション実行者が任意に設定することができる。また、外場が及ぼす力とは、たとえば重力場や電磁場が各粒子に与える力を指す。外場の設定も、ステップＳ１０１またはそれに先立って、シミュレーション実行者が任意に行うことができる。

ステップＳ１０３では、たとえば物体Ｂを構成する任意の隣接粒子ｉ、ｊ、ｋの張る平面（粒子ｉ、ｊ、ｋの定める三角要素の周及び内部）に、物体Ａを構成する任意の粒子ｐが接触しているか否かを判定する。

全ての粒子について接触していない場合には、物体Ａ及びＢを構成する粒子には、ポテンシャル場、及び外場からの力のみが作用するとして、摩擦力の導入は行わず、ステップＳ１０７に進み、物体Ａ及びＢを構成する全ての粒子につき、運動方程式を解いて、時間Δｔ後の粒子の位置を数値計算で求め、粒子を移動させる。

ステップＳ１０３において、粒子ｉ、ｊ、ｋの定める三角要素の周及び内部に、ある粒子ｐが接触している場合、ステップＳ１０４に進み、粒子ｐの接触点を含む領域を支配する粒子ｍを、粒子ｉ、ｊ、及びｋのうちから決定する。決定方法は図１を参照し、前述した通りである。

ステップＳ１０５においては、ステップＳ１０４で決定された粒子ｍと、粒子ｐとの間に摩擦力を導入する。粒子ｐが粒子ｍから受ける摩擦力をｆ_ｐ、粒子ｍが粒子ｐから受ける摩擦力をｆ_ｍとすると、ｆ_ｐ、ｆ_ｍは、式（２１）、（２２）より、たとえば、

・・（２３）

・・（２４）

の両式を用いて計算される。

ステップＳ１０５における摩擦力の導入は、ステップＳ１０３で接触していると判定された全ての粒子、及び、ステップＳ１０４で決定された、それらの粒子の接触点を含む領域を支配する粒子について行われる。なお、摩擦力ｆ_ｐ、ｆ_ｍは、粒子によって異なりうる力である。

ステップＳ１０６において、ステップＳ１０２で導出された力ｆに、ステップＳ１０５で導入された摩擦力を加える。これにより、摩擦関係の設定された粒子に働く力は、ｆ＋ｆ_ｐ、または、ｆ＋ｆ_ｍ（粒子によって異なりうる。）となる。また、他の粒子と摩擦関係にない粒子に働く力はｆ（粒子によって異なりうる。）となる。粒子ごとに働く力は、たとえば物体Ａ及び物体Ｂを構成する全ての粒子につき決定される。

そしてステップＳ１０７において、物体Ａ及びＢを構成する全ての粒子につき、ステップＳ１０６で決定された力を用いて、運動方程式を解き、時間Δｔ後の各粒子の位置を数値計算で求め、各粒子を移動させる。

ステップＳ１０７で全ての粒子についての移動が終わったら、ステップＳ１０８にでシミュレーションを終了するか否かを決定する。終了を選択した場合は、ステップＳ１０９に進み、シミュレーションは完了する。終了しない場合は再びステップＳ１０２に戻り、ステップＳ１０２〜Ｓ１０８の工程を繰り返す。

実施例によるシミュレーション方法によれば、２粒子の相対運動で散逸されるエネルギを、摩擦によって散逸されるエネルギと関連づけて、たとえば両者を同等と考えて、摩擦モデルを導入することで、高精度のシミュレーションを実行することができる。

本願発明者らは、実施例によるシミュレーション方法の妥当性を、滑り摩擦を伴う運動、転がり摩擦を伴う運動の双方について検証した。

図３（Ａ）及び（Ｂ）を参照して、滑り摩擦を伴う運動のシミュレーションの妥当性について説明する。

図３（Ａ）に検証モデルを示す。大きな板の上に小さな板を載置し、小さな板の上面に１０（Ｎ）の荷重を加えながら、側面を１０（Ｎ）の力で引っ張る。上面への荷重は、小さな板の最表面を構成する粒子に均等に与えられ、側面を引く力は、当該側面の最表面を構成する粒子に均等に与えられるものとした。小さな板は、大きな板の上面を、引っ張り方向に等加速度直線運動を行う。

図３（Ｂ）に、摩擦係数を、０．０、０．１及び０．５としたときの加速度の理論値、及び、実施例によるシミュレーション方法を用いて行ったシミュレーションでの加速度値を示す。

摩擦係数が０．０の場合、加速度の理論値は、１．２１７×１０^５（ｍ／ｓ^２）であり、実施例によるシミュレーションにおける加速度は、１．２２０×１０^５（ｍ／ｓ^２）である。また、摩擦係数が０．１、０．５の場合、それぞれ加速度の理論値は、１．０９５×１０^４（ｍ／ｓ^２）、６．０８３×１０^４（ｍ／ｓ^２）であり、実施例によるシミュレーションにおける加速度は、１．０９８×１０^４（ｍ／ｓ^２）、６．１３４×１０^４（ｍ／ｓ^２）である。

実施例によるシミュレーション方法を用いて行ったシミュレーションでの加速度値と理論値との間の誤差は、すべての摩擦係数について１．０％未満である。このことから、実施例によるシミュレーション方法を用いて、高精度にシミュレーションを行えることがわかる。

次に、図４（Ａ）及び（Ｂ）を参照して、転がり摩擦を伴う運動のシミュレーションの妥当性について説明する。転がり摩擦に関しては、円筒の側面を２枚の板で挟み、下側の板を固定し、上側の板を摺動する力学系を使って検証した。上側の板と円筒側面の固定点をマーカー（指標）とし、実施例によるシミュレーション方法でシミュレーションを行ったときのマーカーの動きを確認した。

図４（Ａ）及び（Ｂ）に、シミュレーションにおけるマーカーの位置を示す。両図においては、矢印の方向に時間が進行している。また、図４（Ａ）の右端写真に示される状態と、図４（Ｂ）の左端写真に示される状態との間にも、時間の経過がある。

時間の経過とともに、上側の板に付けられたマーカーと、円筒の側面に付けられたマーカーとが、前者は直線的に、後者は円周に沿って、ほぼ等しい距離だけ移動していることが認められる。実施例によるシミュレーション方法によれば、板の摺動に伴う円筒の転がり運動を適切にシミュレートすることが可能である。

このように、実施例によるシミュレーション方法を用いることによって、滑り摩擦、転がり摩擦を区別することなく、これらを伴う運動を自然に再現することができる。また、摩擦時に発生するエネルギを正しく評価することも可能である。このため、高精度のシミュレーションを行うことができる。

実施例によるシミュレーション方法は、プログラムにより、コンピュータで実行することができる。

図５は、実施例によるシミュレーション方法を用いてシミュレーションを行うシミュレーション装置のシステム構成図である。本図に構成を示したシミュレーション装置を使用して、図２にフローチャートで示したシミュレーション方法を実施することができる。

まず、キーボードなどの入力装置から、図２のＳ１０１に対応して、粒子位置や粒子の速度などの粒子の初期状態が入力される。また、これと同時に、またはこれに先立って、キーボードなどの入力装置から、粒子間のポテンシャルや外場を表現する関数を設定することができる。

中央処理装置は、メインメモリ中の制御プログラムの指令を受け、メインメモリ中の粒子の移動位置決定プログラムを実行する。

粒子の移動位置決定プログラムは、粒子間のポテンシャル及び外場から作用する力を、全ての粒子について導出する部分（図２のＳ１０２で示される部分）、物体Ｂを構成する隣接粒子ｉ、ｊ、及びｋの定める三角要素の周及び内部に、物体Ａを構成する粒子ｐが接触しているとき、粒子ｐの接触点を含む領域を支配する粒子（∈ｉ、ｊ、ｋ）を決定する部分（図２のＳ１０３及びＳ１０４で示される部分）、粒子ｐと、粒子ｐの接触点を含む領域を支配する粒子（∈ｉ、ｊ、ｋ）とが相互に及ぼし合う摩擦力を計算し、これを粒子間ポテンシャル及び外場が与える力に加算する部分（図２のＳ１０５及びＳ１０６で示される部分）、全ての粒子につき運動方程式を解いて移動位置を決定する部分（図２のＳ１０７で示される部分）に大別される。

図２を参照して説明したように、まず、入力装置から入力された情報に基づいて、全ての粒子について、粒子ポテンシャル場及び外場から受ける力が算出される。

次に、隣接粒子ｉ、ｊ、及びｋの定める三角要素の周及び内部に、粒子ｐが接触しているか否かが判定され、接触している場合には、粒子ｐの接触点を含む領域を支配する粒子（∈ｉ、ｊ、ｋ）が決定される。

そして粒子ｐ、及び、粒子ｐの接触点を含む領域を支配する粒子（∈ｉ、ｊ、ｋ）に作用する摩擦力が計算される。計算は、たとえば式（２３）、（２４）に基づいて行われる。摩擦力の算出は、物体Ｂを構成する任意の隣接３粒子の定める三角要素の周及び内部に接触している物体Ａの全ての粒子、及びそれらの粒子の接触点を含む領域を支配する物体Ｂの粒子の全てについて行われる。

摩擦力は粒子間ポテンシャル場及び外場が与える力に加算され、各々の粒子について作用する力が求められる。

それから全ての粒子につき運動方程式を解いて、粒子の移動位置を決定する。

粒子の移動位置が決定されたら、たとえば２次元的、または３次元的に定められた座標系において、決定された移動位置に粒子を表示する。表示結果は、出力装置、たとえばディスプレイに表示される。

以上実施例に沿って本発明を説明したが、本発明はこれらに制限されるものではない。

たとえば、実施例においては、図１を参照して支配領域を２次元的に説明したが、粒子の支配領域を３次元的に把握することも可能である。たとえば、粒子ｉ、ｊ、ｋ等の支配領域を各粒子の含まれるボロノイ多面体の内部とし、物体Ａの構成粒子ｐが、物体Ｂの構成粒子ｉ、ｊ、ｋ等のいずれかの支配領域に属するようになったとき、粒子ｐと、粒子ｐが属するボロノイ多面体内の粒子との間に、式（２１）、（２２）で表される摩擦力を導入する。なお、ボロノイ面上の粒子ｐについては、当該ボロノイ面をはさむいずれの粒子の支配領域下においてもよい。

また、実施例においては、粒子ｐと、粒子ｉ、ｊ、ｋのいずれか１つの間に摩擦関係を導入した。粒子ｐと、複数粒子間に摩擦力を導入することもできる。この場合、粒子ｐと、各粒子との間に、重み付けをして摩擦力を設定することが可能である。

その他、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。

摩擦を伴う種々の現象のシミュレーション、特に機構や弾性の動解析に好適に利用することができる。機構や弾性の動解析における、摩擦によるロス、及び発熱エネルギの見積もりに応用することができる。

相互間に摩擦力を生じる粒子の決定方法について説明する図である。 実施例によるシミュレーション方法を説明するフローチャートである。 （Ａ）及び（Ｂ）は、滑り摩擦を伴う運動のシミュレーションの妥当性について説明する図である。 （Ａ）及び（Ｂ）は、転がり摩擦を伴う運動のシミュレーションの妥当性について説明する図である。 実施例によるシミュレーション方法を用いてシミュレーションを行うシミュレーション装置のシステム構成図である。 （Ａ）及び（Ｂ）は、摩擦力を導入する従来モデルにおける問題点を説明するための図である。

符号の説明

Ｓ１０１〜Ｓ１０９ ステップ

Claims (6)

1. 物体を粒子の集まりで表現し、ある時間刻み幅で、各粒子の挙動を数値計算することにより、時間軸上で離散的に分布する時刻ごとの該物体の状態を解析する粒子法におけるシミュレーション方法であって、
   （ａ）第１の物体を構成する第１の粒子と、第２の物体を構成する第２の粒子との、第１の時刻における相対速度に基づいて、前記第１の粒子と前記第２の粒子との相対運動で散逸されるエネルギが、前記第１の粒子と前記第２の粒子との摩擦によって散逸されるエネルギと等しいとして、前記第１の粒子及び前記第２の粒子に印加される摩擦力を算出する工程と、
   （ｂ）前記第１の粒子及び前記第２の粒子に働く力として、粒子間ポテンシャルに基づく力、外場から受ける力、及び前記工程（ａ）で算出された摩擦力を用いて、運動方程式を解くことにより、前記第１の時刻の次の第２の時刻における、前記第１の粒子、及び前記第２の粒子の位置を求める工程と
   を有するシミュレーション方法。
2. 前記第１の粒子が、前記第２の物体を構成する隣接する３つの粒子で形成される三角形の周または内部に接触し、前記第２の粒子が前記３つの粒子のうち、前記第１の粒子に最も近い粒子である請求項１に記載のシミュレーション方法。
3. 前記第１の粒子が、前記第２の粒子が属するボロノイ多面体内にある請求項１に記載のシミュレーション方法。
4. 物体を粒子の集まりで表現し、ある時間刻み幅で、各粒子の挙動を数値計算することにより、時間軸上で離散的に分布する時刻ごとの該物体の状態を解析する粒子法におけるシミュレーションを行うためにコンピュータを、
   第１の物体を構成する第１の粒子、及び、第２の物体を構成する第２の粒子の初期状態を入力する手段、
   前記第１の粒子と前記第２の粒子との、第１の時刻における相対速度に基づいて、前記第１の粒子と前記第２の粒子との相対運動で散逸されるエネルギが、前記第１の粒子と前
   記第２の粒子との摩擦によって散逸されるエネルギと等しいとして、前記第１の粒子及び前記第２の粒子に印加される摩擦力を算出する手段、
   前記第１の粒子及び前記第２の粒子に働く力として、粒子間ポテンシャルに基づく力、外場から受ける力、及び前記摩擦力を用いて、運動方程式を解くことにより、前記第１の時刻の次の第２の時刻における、前記第１の粒子、及び前記第２の粒子の位置を求める手段、
   前記第２の時刻における前記第１の粒子及び前記第２の粒子を表示する表示手段
   として機能させるためのプログラム。
5. 前記第１の粒子が、前記第２の物体を構成する隣接する３つの粒子で形成される三角形の周または内部に接触し、前記第２の粒子が前記３つの粒子のうち、前記第１の粒子に最も近い粒子である請求項４に記載のプログラム。
6. 前記第１の粒子が、前記第２の粒子が属するボロノイ多面体内にある請求項４に記載のプログラム。