(57)【要約】【課題】  高精度なシミュレーションを実現する。【解決手段】  物体を粒子の集まりで表現し、ある時間刻み幅で、各粒子の挙動を数値計算することにより、時間軸上で離散的に分布する時刻ごとの物体の状態を解析する粒子法におけるシミュレーション方法であって、（ａ）第１の物体を構成する第１の粒子と、第２の物体を構成する第２の粒子との、第１の時刻における相対速度に基づいて、第１の粒子及び第２の粒子に印加される摩擦力を算出する工程と、（ｂ）工程（ａ）で算出された摩擦力を加味して、第１の時刻の次の第２の時刻における、第１の粒子、及び第２の粒子の位置を求める工程とを有するシミュレーション方法を提供する。【選択図】  図２
2010146369.tif 000002
【特許請求の範囲】
【請求項１】
物体を粒子の集まりで表現し、ある時間刻み幅で、各粒子の挙動を数値計算することにより、時間軸上で離散的に分布する時刻ごとの該物体の状態を解析する粒子法におけるシミュレーション方法であって、
（ａ）第１の物体を構成する第１の粒子と、第２の物体を構成する第２の粒子との、第１の時刻における相対速度に基づいて、前記第１の粒子及び前記第２の粒子に印加される摩擦力を算出する工程と、
（ｂ）前記工程（ａ）で算出された摩擦力を加味して、前記第１の時刻の次の第２の時刻における、前記第１の粒子、及び前記第２の粒子の位置を求める工程と
を有するシミュレーション方法。
【請求項２】
前記工程（ａ）において、前記摩擦力を、前記第１の粒子と前記第２の粒子との相対運動で散逸されるエネルギが、前記第１の粒子と前記第２の粒子との摩擦によって散逸されるエネルギと等しいとして算出する請求項１に記載のシミュレーション方法。
【請求項３】
前記第１の粒子が、前記第２の物体を構成する隣接する３つの粒子で形成される三角形の周または内部に接触し、前記第２の粒子が前記３つの粒子のうち、前記第１の粒子に最も近い粒子である請求項１または２に記載のシミュレーション方法。
【請求項４】
前記第１の粒子が、前記第２の粒子が属するボロノイ多面体内にある請求項１または２に記載のシミュレーション方法。
【請求項５】
物体を粒子の集まりで表現し、ある時間刻み幅で、各粒子の挙動を数値計算することにより、時間軸上で離散的に分布する時刻ごとの該物体の状態を解析する粒子法におけるシミュレーションを行うためにコンピュータを、
第１の物体を構成する第１の粒子、及び、第２の物体を構成する第２の粒子の初期状態を入力する手段、
前記第１の粒子と前記第２の粒子との、第１の時刻における相対速度に基づいて、前記第１の粒子及び前記第２の粒子に印加される摩擦力を算出する手段、
前記摩擦力を加味して、前記第１の時刻の次の第２の時刻における、前記第１の粒子、及び前記第２の粒子の位置を求める手段、
前記第２の時刻における前記第１の粒子及び前記第２の粒子を表示する表示手段
として機能させるためのプログラム。
【請求項６】
前記摩擦力を、前記第１の粒子と前記第２の粒子との相対運動で散逸されるエネルギが、前記第１の粒子と前記第２の粒子との摩擦によって散逸されるエネルギと等しいとして算出する請求項５に記載のプログラム。
【請求項７】
前記第１の粒子が、前記第２の物体を構成する隣接する３つの粒子で形成される三角形の周または内部に接触し、前記第２の粒子が前記３つの粒子のうち、前記第１の粒子に最も近い粒子である請求項５または６に記載のプログラム。
【請求項８】
前記第１の粒子が、前記第２の粒子が属するボロノイ多面体内にある請求項５または６に記載のプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、シミュレーション方法及びプログラムに関する。特に、物体を粒子の集まりで表現し、ある時間刻み幅で、各粒子の挙動を数値計算することにより、時間軸上で離散的に分布する時刻ごとの物体の状態を解析する粒子法におけるシミュレーション方法、及びそのシミュレーション方法をコンピュータで実行するプログラムに関する。
【背景技術】
【０００２】
コンピュータを用いて物質科学全般の現象を探求する方法として、分子動力学法に基づく分子シミュレーションが知られている。分子シミュレーションによって、分子のポテンシャルエネルギや最安定構造等物質の特性を、分子レベルで解明することが可能である。
【０００３】
分子動力学法に基づく分子シミュレーションに関して、種々の研究がなされている。たとえば少ない計算量で分子動力学計算を実行するシミュレーション方法の発明が開示されている（たとえば、特許文献１参照）。
【０００４】
粒子シミュレーションの動解析を行う場合、摩擦力のモデル化が問題となる。
【０００５】
現在提案されている代表的な粒子法には、ＭＰＳ(Moving Particle Semi-implicit)、ＳＰＨ(Smoothed Particle Hydrodynamics)、ＤＥＭ(Distinct Element Method)などがある。
【０００６】
これらの手法では、粒子の接触時にせん断力を決定し、その大きさが静止摩擦力を上回らないようにすることで摩擦力をモデル化している。また摩擦力は、接触粒子の運動を妨げる方向、すなわち接触時における粒子の相対速度の方向と反対方向に発生するとしている。
【０００７】
しかしながら粒子は、物体全体の重心運動に関わらず、常に高周波数で振動している。そのため上記の手法においては、粒子に働く摩擦力の方向は、物体全体の重心運動とは関係なく、原子の振動周期と同じオーダーで目まぐるしく変化することになる。ゆえに、これらの手法では、滑り摩擦時に要するエネルギを正しく計算できない可能性が高い。
【０００８】
図６（Ａ）及び（Ｂ）を参照して、上述の内容を詳述する。
【０００９】
図６（Ａ）に、物体全体の摺動の様子を示す。物体は摩擦力Ｆに逆らって、ｘ軸正方向に距離Ｌだけ移動する。物体の重心は常にｘ軸正方向に移動している。この場合、物体には移動中、常に摩擦力Ｆがｘ軸負方向に作用しており、物体に作用する摩擦力の大きさの平均もＦとなる。
【００１０】
図６（Ｂ）に物体の接触部粒子の挙動を示す。上述の粒子法によれば、摩擦力を受けるのは接触部の粒子のみであり、粒子の振動は物体の重心の運動方向と同じであるとは限らない。たとえば接触部の原子は、本図に示す軌跡を辿って、ｘ軸正方向に距離Ｌを移動する場合がある。
【００１１】
本図に示す例においては、接触部の原子は最初の時間ｄｔでｘ軸正方向に距離Ｌ１だけ移動する。続く時間ｄｔではｘ軸負方向に距離Ｌ２だけ移動する。最後の時間ｄｔではｘ軸正方向に距離Ｌ３だけ移動する。ここでＬ１−Ｌ２＋Ｌ３＝Ｌを満たす。このとき距離Ｌを摺動する間に、接触部の原子に働く摩擦力の大きさは、最初の時間ｄｔではｘ軸負方向にＦ、次の時間ｄｔではｘ軸正方向にＦ、最後の時間ｄｔではｘ軸負方向にＦであるから、平均して働く摩擦力は、ｘ軸負方向に、（Ｆ−Ｆ＋Ｆ）／３＝Ｆ／３  となり、図６（Ａ）を参照して考察した結果とは異なってしまう。
【００１２】
このことは時間刻みｄｔを十分小さくしても、摺動距離を長くしても起こりうることであり、防止策としては、粒子の大きさを十分大きくとる、もしくは十分な時間的、空間的なアンサンブルをとることなどが考えられる。
【００１３】
しかしながら、粒子を大きくすると複雑な形状が再現できない、時間アンサンブルを十分とるとフィードバックが必要となり、計算に余分な時間がかかる、空間アンサンブルが十分とれない場合は用いることができない、などの問題点があり実用的ではない。
【００１４】
従来の粒子法における摩擦のモデルでは滑り摩擦によるエネルギ散逸が正しく評価できず、また、滑り摩擦と転がり摩擦とを区別して評価する必要がある。しかし現実の摩擦は、滑り摩擦と転がり摩擦とが混在して生じているため、そこに人為的な仮定を含むことは望ましくない。
【００１５】
【特許文献１】特開２００６−２８５８６６号公報
【発明の開示】
【発明が解決しようとする課題】
【００１６】
本発明の目的は、高精度なシミュレーションを実現するシミュレーション方法を提供することである。
【００１７】
また、高精度なシミュレーションを実現するシミュレーション方法を実行させるプログラムを提供することである。
【課題を解決するための手段】
【００１８】
本発明の一観点によれば、物体を粒子の集まりで表現し、ある時間刻み幅で、各粒子の挙動を数値計算することにより、時間軸上で離散的に分布する時刻ごとの該物体の状態を解析する粒子法におけるシミュレーション方法であって、（ａ）第１の物体を構成する第１の粒子と、第２の物体を構成する第２の粒子との、第１の時刻における相対速度に基づいて、前記第１の粒子及び前記第２の粒子に印加される摩擦力を算出する工程と、（ｂ）前記工程（ａ）で算出された摩擦力を加味して、前記第１の時刻の次の第２の時刻における、前記第１の粒子、及び前記第２の粒子の位置を求める工程とを有するシミュレーション方法が提供される。
【００１９】
また、本発明の他の観点によると、物体を粒子の集まりで表現し、ある時間刻み幅で、各粒子の挙動を数値計算することにより、時間軸上で離散的に分布する時刻ごとの該物体の状態を解析する粒子法におけるシミュレーションを行うためにコンピュータを、第１の物体を構成する第１の粒子、及び、第２の物体を構成する第２の粒子の初期状態を入力する手段、前記第１の粒子と前記第２の粒子との、第１の時刻における相対速度に基づいて、前記第１の粒子及び前記第２の粒子に印加される摩擦力を算出する手段、前記摩擦力を加味して、前記第１の時刻の次の第２の時刻における、前記第１の粒子、及び前記第２の粒子の位置を求める手段、前記第２の時刻における前記第１の粒子及び前記第２の粒子を表示する表示手段として機能させるためのプログラムが提供される。
【発明の効果】
【００２０】
本発明によれば、高精度なシミュレーションを実現するシミュレーション方法を提供することができる。
【００２１】
また、高精度なシミュレーションを実現するシミュレーション方法を実行させるプログラムを提供することができる。
【発明を実施するための最良の形態】
【００２２】
以下、物体を粒子の集まりで表現し、ある時間刻み幅Δｔで、各粒子の挙動を数値計算することにより、時間軸上で離散的に分布する時刻ごとの物体の状態を解析する粒子法におけるシミュレーション方法、及びそのシミュレーション方法をコンピュータで実行するプログラムの実施例について説明する。
【００２３】
質量ｍｉの粒子ｉと質量ｍｊの粒子ｊが摩擦した場合を考える。２つの粒子についての散逸方程式は、拡散係数γ（ｔ）を用いて、下の式（１）及び（２）で表される。
2010146369.tif 000003
【００２４】
・・（１）
2010146369.tif 000004
【００２５】
・・（２）
ここでｖｉ、ｖｊは、それぞれ粒子ｉ、ｊの速度を示す。
【００２６】
また、粒子ｊに対する粒子ｉの相対速度をｖｒ、粒子ｉと粒子ｊの重心の速度をｖｇとすると、ｖｒとｖｇは、それぞれ下の（３）、（４）式で表される。
2010146369.tif 000005
【００２７】
・・（３）
2010146369.tif 000006
【００２８】
・・（４）
（３）、（４）式を、ｖｉ、ｖｊについて解くと、下の（５）、（６）式を得る。
2010146369.tif 000007
【００２９】
・・（５）
2010146369.tif 000008
【００３０】
・・（６）
式（３）、（５）を式（１）に代入して得られた式と、式（３）、（６）を式（２）に代入して得られた式とを連立することで、以下の式（７）及び（８）が導かれる。
2010146369.tif 000009
【００３１】
・・（７）
2010146369.tif 000010
【００３２】
・・（８）
ここで、
2010146369.tif 000011
【００３３】
・・（９）
である。
式（８）を解くと、
2010146369.tif 000012
【００３４】
・・（１０）
が得られる。
【００３５】
式（７）、（１０）より、散逸過程においては、重心の運動量は保存され、相対運動のエネルギが散逸されることがわかる。
【００３６】
以下、数値積分を実行する。式（８）を離散化すると、ある時刻ｎについて、
2010146369.tif 000013
【００３７】
・・（１１）
を得る。
【００３８】
式（１１）を、ｖｒ［ｎ＋１／２］及びｖｒ［ｎ−１／２］についてまとめると、
2010146369.tif 000014
【００３９】
・・（１２）
となり、更にこれは、
2010146369.tif 000015
【００４０】
・・（１３）
と変形することができる。
【００４１】
一方、粒子ｉについての離散化された散逸方程式は、式（１）、（８）、（１１）より、
2010146369.tif 000016
【００４２】
・・（１４）
と書ける。式（１３）を式（１４）に代入すると、
2010146369.tif 000017
【００４３】
・・（１５）
を得、これを整理すると、
2010146369.tif 000018
【００４４】
・・（１６）
となる。
【００４５】
ここで本願発明者らは、従来のように粒子に阻止力を加える摩擦モデルではなく、２粒子の相対運動で散逸されるエネルギを、摩擦によって散逸されるエネルギと関連づけて考える、たとえば両者を同等と考える摩擦モデルを導入することにした。
【００４６】
拡散係数γとクーロン摩擦力とを関係づけるために、２粒子の相対運動による散逸パワー（単位時間当たりの散逸エネルギ）と、摩擦による散逸パワーとを等しいとして、
2010146369.tif 000019
【００４７】
・・（１７）
すなわち、
2010146369.tif 000020
【００４８】
・・（１８）
という関係式を導入する。ここで、μは摩擦係数、Ｎは垂直抗力である。
【００４９】
式（１６）に式（１８）を代入し、代入した式の右辺に、粒子ｉに働く力のうち、粒子ｊとの間の摩擦力を除く力をｆｉとして加えると、粒子ｉについての離散化された運動方程式（１９）が得られる。
2010146369.tif 000021
【００５０】
・・（１９）
同様に、粒子ｊについても、粒子ｊに働く力のうち、粒子ｉとの間の摩擦力を除く力ｆｊを用いて、離散化された運動方程式（２０）を得ることができる。
2010146369.tif 000022
【００５１】
・・（２０）
力ｆｉ、ｆｊには、たとえば、重力、電磁力、原子間力等が含まれる。粒子ｉ、ｊの運動は、運動方程式（１９）、（２０）を解くことによって把握される。式（１９）、（２０）から明らかなように、本願発明者らの提案する摩擦モデルによれば、粒子ｉが粒子ｊから受ける摩擦力Ｆｉは、
2010146369.tif 000023
【００５２】
・・（２１）
で表され、粒子ｊが粒子ｉから受ける摩擦力Ｆｊは、
2010146369.tif 000024
【００５３】
・・（２２）
で表される。
【００５４】
図１及び図２を参照して、実施例によるシミュレーション方法について説明する。
【００５５】
まず、シミュレーション方法全体の説明に先立って、図１を参照し、相互間に摩擦力を生じる粒子の決定方法について説明する。
【００５６】
たとえば物体Ａを構成する粒子の１つである粒子ｐが、物体Ｂを構成する隣接粒子である粒子ｉ、ｊ、及びｋの中心の張る平面において、粒子ｉ、ｊ、ｋの定める三角形（三角要素）に接触したとする。ここで隣接粒子とは、多数の粒子で構成される物体Ｂをボロノイ分割したときに、隣り合う粒子をいう。また、以下粒子の位置が問題となるときには、粒子の中心の位置を粒子の位置と考える。
【００５７】
粒子ｉ、ｊ、ｋの定める三角形の外心をＣ、外心Ｃから線分ｉｊに下ろした垂線の足をＨｉｊ、線分ｊｋに下ろした垂線の足をＨｊｋ、線分ｋｉに下ろした垂線の足をＨｋｉとし、四角形ｉＨｉｊＣＨｋｉの周及び内部から、線分ＣＨｉｊ及び線分ＣＨｋｉを除いた領域を、粒子ｉの支配領域と定める。同様に、四角形ｊＨｊｋＣＨｉｊの周及び内部から、線分ＣＨｊｋ及び線分ＣＨｉｊを除いた領域を粒子ｊの支配領域、四角形ｋＨｋｉＣＨｊｋの周及び内部から、線分ＣＨｋｉ及び線分ＣＨｊｋを除いた領域を粒子ｋの支配領域と定める。
【００５８】
また、線分ＣＨｉｊ上の点は、粒子ｉの支配領域、粒子ｊの支配領域のいずれかに含める。同様に、線分ＣＨｊｋ上の点は、粒子ｊ、ｋのいずれかの支配領域に含め、線分ＣＨｋｉ上の点は、粒子ｋ、ｉのいずれかの支配領域に含める。
【００５９】
そして粒子ｐが、粒子ｉ、ｊ、ｋのいずれの粒子の支配領域に存在するかを確認し、粒子ｐが存在する領域を支配する粒子との間に摩擦力が発生すると考える。たとえば粒子ｐが粒子ｉの支配領域に存在する場合、粒子ｐと粒子ｉとの間に、式（２１）、（２２）で表される摩擦力を導入する。そして両粒子ｐ、ｉは、式（１９）、（２０）に則った相互作用をすると考える。また、粒子ｐが粒子ｉの支配領域に存在する場合、粒子ｐと粒子ｊ、ｋ間の摩擦力は０とする。
【００６０】
なお、上記摩擦力導入の考え方によれば、三角形ｉｊｋの周または内部に接触した粒子ｐは、粒子ｉ、ｊ、ｋのうち最も近い位置にある粒子との間に摩擦力を導入されることになる。
【００６１】
次に、図２に示すフローチャートを用いて、実施例によるシミュレーション方法を説明する。
【００６２】
実施例によるシミュレーション方法では、まずステップＳ１０１において、粒子の位置や速度など、シミュレーションを行うに当たっての初期状態を決定する。
【００６３】
次に、ステップＳ１０２において、粒子間のポテンシャル、及び外場が及ぼす力ｆを全ての粒子につき、粒子ごとに導出する。ｆは粒子によって異なりうる力である。
【００６４】
物体を構成する粒子同士は、互いに粒子間のポテンシャルにより結びついている。粒子間のポテンシャルは、たとえばバネやＬｅｎｎａｒｄ−Ｊｏｎｅｓポテンシャルで表現される。どのようなポテンシャルを用いるかについては、ステップＳ１０１またはそれに先立って、シミュレーション実行者が任意に設定することができる。また、外場が及ぼす力とは、たとえば重力場や電磁場が各粒子に与える力を指す。外場の設定も、ステップＳ１０１またはそれに先立って、シミュレーション実行者が任意に行うことができる。
【００６５】
ステップＳ１０３では、たとえば物体Ｂを構成する任意の隣接粒子ｉ、ｊ、ｋの張る平面（粒子ｉ、ｊ、ｋの定める三角要素の周及び内部）に、物体Ａを構成する任意の粒子ｐが接触しているか否かを判定する。
【００６６】
全ての粒子について接触していない場合には、物体Ａ及びＢを構成する粒子には、ポテンシャル場、及び外場からの力のみが作用するとして、摩擦力の導入は行わず、ステップＳ１０７に進み、物体Ａ及びＢを構成する全ての粒子につき、運動方程式を解いて、時間Δｔ後の粒子の位置を数値計算で求め、粒子を移動させる。
【００６７】
ステップＳ１０３において、粒子ｉ、ｊ、ｋの定める三角要素の周及び内部に、ある粒子ｐが接触している場合、ステップＳ１０４に進み、粒子ｐの接触点を含む領域を支配する粒子ｍを、粒子ｉ、ｊ、及びｋのうちから決定する。決定方法は図１を参照し、前述した通りである。
【００６８】
ステップＳ１０５においては、ステップＳ１０４で決定された粒子ｍと、粒子ｐとの間に摩擦力を導入する。粒子ｐが粒子ｍから受ける摩擦力をｆｐ、粒子ｍが粒子ｐから受ける摩擦力をｆｍとすると、ｆｐ、ｆｍは、式（２１）、（２２）より、たとえば、
2010146369.tif 000025
【００６９】
・・（２３）
2010146369.tif 000026
【００７０】
・・（２４）
の両式を用いて計算される。
【００７１】
ステップＳ１０５における摩擦力の導入は、ステップＳ１０３で接触していると判定された全ての粒子、及び、ステップＳ１０４で決定された、それらの粒子の接触点を含む領域を支配する粒子について行われる。なお、摩擦力ｆｐ、ｆｍは、粒子によって異なりうる力である。
【００７２】
ステップＳ１０６において、ステップＳ１０２で導出された力ｆに、ステップＳ１０５で導入された摩擦力を加える。これにより、摩擦関係の設定された粒子に働く力は、ｆ＋ｆｐ、または、ｆ＋ｆｍ（粒子によって異なりうる。）となる。また、他の粒子と摩擦関係にない粒子に働く力はｆ（粒子によって異なりうる。）となる。粒子ごとに働く力は、たとえば物体Ａ及び物体Ｂを構成する全ての粒子につき決定される。
【００７３】
そしてステップＳ１０７において、物体Ａ及びＢを構成する全ての粒子につき、ステップＳ１０６で決定された力を用いて、運動方程式を解き、時間Δｔ後の各粒子の位置を数値計算で求め、各粒子を移動させる。
【００７４】
ステップＳ１０７で全ての粒子についての移動が終わったら、ステップＳ１０８にでシミュレーションを終了するか否かを決定する。終了を選択した場合は、ステップＳ１０９に進み、シミュレーションは完了する。終了しない場合は再びステップＳ１０２に戻り、ステップＳ１０２〜Ｓ１０８の工程を繰り返す。
【００７５】
実施例によるシミュレーション方法によれば、２粒子の相対運動で散逸されるエネルギを、摩擦によって散逸されるエネルギと関連づけて、たとえば両者を同等と考えて、摩擦モデルを導入することで、高精度のシミュレーションを実行することができる。
【００７６】
本願発明者らは、実施例によるシミュレーション方法の妥当性を、滑り摩擦を伴う運動、転がり摩擦を伴う運動の双方について検証した。
【００７７】
図３（Ａ）及び（Ｂ）を参照して、滑り摩擦を伴う運動のシミュレーションの妥当性について説明する。
【００７８】
図３（Ａ）に検証モデルを示す。大きな板の上に小さな板を載置し、小さな板の上面に１０（Ｎ）の荷重を加えながら、側面を１０（Ｎ）の力で引っ張る。上面への荷重は、小さな板の最表面を構成する粒子に均等に与えられ、側面を引く力は、当該側面の最表面を構成する粒子に均等に与えられるものとした。小さな板は、大きな板の上面を、引っ張り方向に等加速度直線運動を行う。
【００７９】
図３（Ｂ）に、摩擦係数を、０．０、０．１及び０．５としたときの加速度の理論値、及び、実施例によるシミュレーション方法を用いて行ったシミュレーションでの加速度値を示す。
【００８０】
摩擦係数が０．０の場合、加速度の理論値は、１．２１７×１０５（ｍ／ｓ２）であり、実施例によるシミュレーションにおける加速度は、１．２２０×１０５（ｍ／ｓ２）である。また、摩擦係数が０．１、０．５の場合、それぞれ加速度の理論値は、１．０９５×１０４（ｍ／ｓ２）、６．０８３×１０４（ｍ／ｓ２）であり、実施例によるシミュレーションにおける加速度は、１．０９８×１０４（ｍ／ｓ２）、６．１３４×１０４（ｍ／ｓ２）である。
【００８１】
実施例によるシミュレーション方法を用いて行ったシミュレーションでの加速度値と理論値との間の誤差は、すべての摩擦係数について１．０％未満である。このことから、実施例によるシミュレーション方法を用いて、高精度にシミュレーションを行えることがわかる。
【００８２】
次に、図４（Ａ）及び（Ｂ）を参照して、転がり摩擦を伴う運動のシミュレーションの妥当性について説明する。転がり摩擦に関しては、円筒の側面を２枚の板で挟み、下側の板を固定し、上側の板を摺動する力学系を使って検証した。上側の板と円筒側面の固定点をマーカー（指標）とし、実施例によるシミュレーション方法でシミュレーションを行ったときのマーカーの動きを確認した。
【００８３】
図４（Ａ）及び（Ｂ）に、シミュレーションにおけるマーカーの位置を示す。両図においては、矢印の方向に時間が進行している。また、図４（Ａ）の右端写真に示される状態と、図４（Ｂ）の左端写真に示される状態との間にも、時間の経過がある。
【００８４】
時間の経過とともに、上側の板に付けられたマーカーと、円筒の側面に付けられたマーカーとが、前者は直線的に、後者は円周に沿って、ほぼ等しい距離だけ移動していることが認められる。実施例によるシミュレーション方法によれば、板の摺動に伴う円筒の転がり運動を適切にシミュレートすることが可能である。
【００８５】
このように、実施例によるシミュレーション方法を用いることによって、滑り摩擦、転がり摩擦を区別することなく、これらを伴う運動を自然に再現することができる。また、摩擦時に発生するエネルギを正しく評価することも可能である。このため、高精度のシミュレーションを行うことができる。
【００８６】
実施例によるシミュレーション方法は、プログラムにより、コンピュータで実行することができる。
【００８７】
図５は、実施例によるシミュレーション方法を用いてシミュレーションを行うシミュレーション装置のシステム構成図である。本図に構成を示したシミュレーション装置を使用して、図２にフローチャートで示したシミュレーション方法を実施することができる。
【００８８】
まず、キーボードなどの入力装置から、図２のＳ１０１に対応して、粒子位置や粒子の速度などの粒子の初期状態が入力される。また、これと同時に、またはこれに先立って、キーボードなどの入力装置から、粒子間のポテンシャルや外場を表現する関数を設定することができる。
【００８９】
中央処理装置は、メインメモリ中の制御プログラムの指令を受け、メインメモリ中の粒子の移動位置決定プログラムを実行する。
【００９０】
粒子の移動位置決定プログラムは、粒子間のポテンシャル及び外場から作用する力を、全ての粒子について導出する部分（図２のＳ１０２で示される部分）、物体Ｂを構成する隣接粒子ｉ、ｊ、及びｋの定める三角要素の周及び内部に、物体Ａを構成する粒子ｐが接触しているとき、粒子ｐの接触点を含む領域を支配する粒子（∈ｉ、ｊ、ｋ）を決定する部分（図２のＳ１０３及びＳ１０４で示される部分）、粒子ｐと、粒子ｐの接触点を含む領域を支配する粒子（∈ｉ、ｊ、ｋ）とが相互に及ぼし合う摩擦力を計算し、これを粒子間ポテンシャル及び外場が与える力に加算する部分（図２のＳ１０５及びＳ１０６で示される部分）、全ての粒子につき運動方程式を解いて移動位置を決定する部分（図２のＳ１０７で示される部分）に大別される。
【００９１】
図２を参照して説明したように、まず、入力装置から入力された情報に基づいて、全ての粒子について、粒子ポテンシャル場及び外場から受ける力が算出される。
【００９２】
次に、隣接粒子ｉ、ｊ、及びｋの定める三角要素の周及び内部に、粒子ｐが接触しているか否かが判定され、接触している場合には、粒子ｐの接触点を含む領域を支配する粒子（∈ｉ、ｊ、ｋ）が決定される。
【００９３】
そして粒子ｐ、及び、粒子ｐの接触点を含む領域を支配する粒子（∈ｉ、ｊ、ｋ）に作用する摩擦力が計算される。計算は、たとえば式（２３）、（２４）に基づいて行われる。摩擦力の算出は、物体Ｂを構成する任意の隣接３粒子の定める三角要素の周及び内部に接触している物体Ａの全ての粒子、及びそれらの粒子の接触点を含む領域を支配する物体Ｂの粒子の全てについて行われる。
【００９４】
摩擦力は粒子間ポテンシャル場及び外場が与える力に加算され、各々の粒子について作用する力が求められる。
【００９５】
それから全ての粒子につき運動方程式を解いて、粒子の移動位置を決定する。
【００９６】
粒子の移動位置が決定されたら、たとえば２次元的、または３次元的に定められた座標系において、決定された移動位置に粒子を表示する。表示結果は、出力装置、たとえばディスプレイに表示される。
【００９７】
以上実施例に沿って本発明を説明したが、本発明はこれらに制限されるものではない。
【００９８】
たとえば、実施例においては、図１を参照して支配領域を２次元的に説明したが、粒子の支配領域を３次元的に把握することも可能である。たとえば、粒子ｉ、ｊ、ｋ等の支配領域を各粒子の含まれるボロノイ多面体の内部とし、物体Ａの構成粒子ｐが、物体Ｂの構成粒子ｉ、ｊ、ｋ等のいずれかの支配領域に属するようになったとき、粒子ｐと、粒子ｐが属するボロノイ多面体内の粒子との間に、式（２１）、（２２）で表される摩擦力を導入する。なお、ボロノイ面上の粒子ｐについては、当該ボロノイ面をはさむいずれの粒子の支配領域下においてもよい。
【００９９】
また、実施例においては、粒子ｐと、粒子ｉ、ｊ、ｋのいずれか１つの間に摩擦関係を導入した。粒子ｐと、複数粒子間に摩擦力を導入することもできる。この場合、粒子ｐと、各粒子との間に、重み付けをして摩擦力を設定することが可能である。
【０１００】
その他、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。
【産業上の利用可能性】
【０１０１】
摩擦を伴う種々の現象のシミュレーション、特に機構や弾性の動解析に好適に利用することができる。機構や弾性の動解析における、摩擦によるロス、及び発熱エネルギの見積もりに応用することができる。
【図面の簡単な説明】
【０１０２】
【図１】相互間に摩擦力を生じる粒子の決定方法について説明する図である。
【図２】実施例によるシミュレーション方法を説明するフローチャートである。
【図３】（Ａ）及び（Ｂ）は、滑り摩擦を伴う運動のシミュレーションの妥当性について説明する図である。
【図４】（Ａ）及び（Ｂ）は、転がり摩擦を伴う運動のシミュレーションの妥当性について説明する図である。
【図５】実施例によるシミュレーション方法を用いてシミュレーションを行うシミュレーション装置のシステム構成図である。
【図６】（Ａ）及び（Ｂ）は、摩擦力を導入する従来モデルにおける問題点を説明するための図である。
【符号の説明】
【０１０３】
Ｓ１０１〜Ｓ１０９  ステップ
【図１】
2010146369.tif 000027
【図２】
2010146369.tif 000028
【図５】
2010146369.tif 000029
【図６】
2010146369.tif 000030
【図３】
2010146369.tif 000031
【図４】
2010146369.tif 000032