(57)【要約】【課題】  高精度なシミュレーションを実現する。【解決手段】  （ａ）第１の時刻において、第１の物体を構成する第１の粒子が、第２の物体の表面から離れており、第１の時刻の次の第２の時刻において、第１の粒子が第２の物体の表面よりも深い場所に位置するとき、第２の時刻における第１の粒子の位置を、凝着点として記憶する。（ｂ）第２の時刻の次の第３の時刻における第１の粒子の位置を数値計算により求める。（ｃ）第３の時刻における第１の粒子の位置と、凝着点の位置とに基づいて、第１の粒子に印加される摩擦力を算出する。（ｄ）摩擦力を加味して、第３の時刻の次の第４の時刻における第１の粒子の位置を数値計算により求める。【選択図】  図２
2010044710.tif 000002
【特許請求の範囲】
【請求項１】
物体を粒子の集まりで表現し、ある時間刻み幅で、各粒子の挙動を数値計算することにより、時間軸上で離散的に分布する時刻ごとの該物体の状態を解析する粒子法におけるシミュレーション方法であって、
（ａ）第１の時刻において、第１の物体を構成する第１の粒子が、第２の物体の表面から離れており、該第１の時刻の次の第２の時刻において、前記第１の粒子が前記第２の物体の表面よりも深い場所に位置するとき、該第２の時刻における前記第１の粒子の位置を、凝着点として記憶する工程と、
（ｂ）前記第２の時刻の次の第３の時刻における前記第１の粒子の位置を数値計算により求める工程と、
（ｃ）前記第３の時刻における前記第１の粒子の位置と、前記凝着点の位置とに基づいて、前記第１の粒子に印加される摩擦力を算出する工程と、
（ｄ）前記摩擦力を加味して、前記第３の時刻の次の第４の時刻における前記第１の粒子の位置を数値計算により求める工程と
を有するシミュレーション方法。
【請求項２】
前記工程（ｃ）において、前記第１の粒子に印加される摩擦力を、前記第３の時刻における前記第１の粒子の位置より前記第２の物体の表面に下ろした垂線の足から、前記凝着点より前記第２の物体の表面に下ろした垂線の足に向かうように算出する請求項１に記載のシミュレーション方法。
【請求項３】
前記工程（ｃ）において、前記第１の粒子に印加される摩擦力を、前記凝着点の前記第２の物体の表面からの深さと、前記第３の時刻における前記第１の粒子の位置の前記第２の物体の表面からの深さとのうち、小さくない方の値に比例した大きさをもつように算出する請求項１または２に記載のシミュレーション方法。
【請求項４】
前記工程（ｃ）が、前記第３の時刻における前記第１の粒子の位置より前記第２の物体の表面に下ろした垂線の足から、前記凝着点より前記第２の物体の表面に下ろした垂線の足までの距離を示す第１の値と、判定基準値との大小を判定する工程を含み、
前記第１の値が前記判定基準値よりも大きい場合には、前記第１の粒子に印加される摩擦力を算出せず、前記第１の値が前記判定基準値以下の場合に、前記第１の粒子に印加される摩擦力を算出する請求項１〜３のいずれか１項に記載のシミュレーション方法。
【請求項５】
前記第４の時刻以降の第ｍの時刻に至るまで、前記第１の粒子の位置より前記第２の物体の表面に下ろした垂線の足から、前記凝着点より前記第２の物体の表面に下ろした垂線の足までの距離が、前記判定基準値以下の場合には、前記第ｍの時刻における前記第１の粒子に印加される摩擦力であって、前記第ｍの時刻における前記第１の粒子の位置より前記第２の物体の表面に下ろした垂線の足から、前記凝着点より前記第２の物体の表面に下ろした垂線の足に向かい、前記第２〜第ｍの時刻における前記第１の粒子の前記第２の物体の表面からの深さの中で最も深い深さに比例した大きさをもつ摩擦力を算出し、算出された前記摩擦力を加味して、前記第ｍの時刻の次の第ｍ＋１の時刻における前記第１の粒子の位置を数値計算により求める請求項４に記載のシミュレーション方法。
【請求項６】
物体を粒子の集まりで表現し、ある時間刻み幅で、各粒子の挙動を数値計算することにより、時間軸上で離散的に分布する時刻ごとの該物体の状態を解析する粒子法におけるシミュレーションを行うためにコンピュータを、
第１の物体を構成する第１の粒子の初期状態を入力する手段、
第１の時刻において、前記第１の粒子が、第２の物体の表面から離れており、該第１の時刻の次の第２の時刻において、前記第１の粒子が前記第２の物体の表面よりも深い場所に位置するとき、該第２の時刻における前記第１の粒子の位置を、凝着点として画定する手段、
前記第２の時刻の次の第３の時刻における前記第１の粒子の位置を決定する手段、
前記第３の時刻における、前記第１の粒子に作用するポテンシャル力を得る手段、
前記第３の時刻における前記第１の粒子の位置と、前記凝着点の位置とに基づいて、前記第３の時刻に、前記第１の粒子に印加される摩擦力を得る手段、
前記第３の時刻におけるポテンシャル力と摩擦力とを加えて、前記第３の時刻の次の第４の時刻における前記第１の粒子の位置を決定する手段、
前記第４の時刻における前記第１の粒子を表示する表示手段
として機能させるためのプログラム。
【請求項７】
前記第１の粒子に印加される摩擦力を、前記第３の時刻における前記第１の粒子の位置より前記第２の物体の表面に下ろした垂線の足から、前記凝着点より前記第２の物体の表面に下ろした垂線の足に向かうように得る請求項６に記載のプログラム。
【請求項８】
前記第１の粒子に印加される摩擦力を、前記凝着点の前記第２の物体の表面からの深さと、前記第３の時刻における前記第１の粒子の位置の前記第２の物体の表面からの深さとのうち、小さくない方の値に比例した大きさをもつように得る請求項６または７に記載のプログラム。
【請求項９】
更に、コンピュータを、前記第３の時刻における前記第１の粒子の位置より前記第２の物体の表面に下ろした垂線の足から、前記凝着点より前記第２の物体の表面に下ろした垂線の足までの距離を示す第１の値と、判定基準値との大小を判定する手段として機能させ、
前記大小を判定する手段が、前記第１の値が前記判定基準値よりも大きいと判定した場合には、前記摩擦力を得る手段は、前記第３の時刻に前記第１の粒子に印加される摩擦力をゼロとする請求項６〜８のいずれか１項に記載のプログラム。
【請求項１０】
更に、コンピュータを、
前記第４の時刻以降の第ｍの時刻における、前記第１の粒子に作用するポテンシャル力を得る手段、
前記第ｍの時刻に至るまで、前記第１の粒子の位置より前記第２の物体の表面に下ろした垂線の足から、前記凝着点より前記第２の物体の表面に下ろした垂線の足までの距離が、前記判定基準値以下の場合には、前記第ｍの時刻における前記第１の粒子に印加される摩擦力であって、前記第ｍの時刻における前記第１の粒子の位置より前記第２の物体の表面に下ろした垂線の足から、前記凝着点より前記第２の物体の表面に下ろした垂線の足に向かい、前記第２〜第ｍの時刻における前記第１の粒子の前記第２の物体の表面からの深さの中で最も深い深さに比例した大きさをもつ摩擦力を得る手段、
前記第ｍの時刻におけるポテンシャル力と摩擦力とを加えて、前記第ｍの時刻の次の第ｍ＋１の時刻における前記第１の粒子の位置を決定する手段、
前記第ｍ＋１の時刻における前記第１の粒子を表示する表示手段
として機能させる請求項９に記載のプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、シミュレーション方法及びプログラムに関する。特に、物体を粒子の集まりで表現し、ある時間刻み幅で、各粒子の挙動を数値計算することにより、時間軸上で離散的に分布する時刻ごとの物体の状態を解析する粒子法におけるシミュレーション方法、及びそのシミュレーション方法をコンピュータで実行するプログラムに関する。
【背景技術】
【０００２】
コンピュータを用いて物質科学全般の現象を探求する方法として、分子動力学法に基づく分子シミュレーションが知られている。分子シミュレーションによって、分子のポテンシャルエネルギや最安定構造等物質の特性を、分子レベルで解明することが可能である。
【０００３】
分子動力学法に基づく分子シミュレーションに関して、種々の研究がなされている。たとえば少ない計算量で分子動力学計算を実行するシミュレーション方法の発明が開示されている（たとえば、特許文献１参照）。
【０００４】
粒子シミュレーションの動解析を行う場合、摩擦力のモデル化が問題となる。
【０００５】
現在提案されている代表的な粒子法には、ＭＰＳ(Moving Particle Semi-implicit)、ＳＰＨ(Smoothed Particle Hydrodynamics)、ＤＥＭ(Distinct Element Method)などがある。
【０００６】
これらの手法では、粒子の接触時にせん断力を決定し、その大きさが静止摩擦力を上回らないようにすることで摩擦力をモデル化している。また摩擦力は、接触粒子の運動を妨げる方向、すなわち接触時における粒子の相対速度の方向と反対方向に発生するとしている。
【０００７】
しかしながら粒子は、物体全体の重心運動に関わらず、常に高周波数で振動している。そのため上記の手法においては、粒子に働く摩擦力の方向は、物体全体の重心運動とは関係なく、原子の振動周期と同じオーダーで目まぐるしく変化することになる。ゆえに、これらの手法では、滑り摩擦時に要するエネルギを正しく計算できない可能性が高い。
【０００８】
図５（Ａ）及び（Ｂ）を参照して、上述の内容を詳述する。
【０００９】
図５（Ａ）に、物体全体の摺動の様子を示す。物体は摩擦力Ｆに逆らって、ｘ軸正方向に距離Ｌだけ移動する。物体の重心は常にｘ軸正方向に移動している。この場合、物体には移動中、常に摩擦力Ｆがｘ軸負方向に作用しており、物体に作用する摩擦力の大きさの平均もＦとなる。
【００１０】
図５（Ｂ）に物体の接触部粒子の挙動を示す。上述の粒子法によれば、摩擦力を受けるのは接触部の粒子のみであり、粒子の振動は物体の重心の運動方向と同じであるとは限らない。たとえば接触部の原子は、本図に示す軌跡を辿って、ｘ軸正方向に距離Ｌを移動する場合がある。
【００１１】
本図に示す例においては、接触部の原子は最初の時間ｄｔでｘ軸正方向に距離Ｌ１だけ移動する。続く時間ｄｔではｘ軸負方向に距離Ｌ２だけ移動する。最後の時間ｄｔではｘ軸正方向に距離Ｌ３だけ移動する。ここでＬ１−Ｌ２＋Ｌ３＝Ｌを満たす。このとき距離Ｌを摺動する間に、接触部の原子に働く摩擦力の大きさは、最初の時間ｄｔではｘ軸負方向にＦ、次の時間ｄｔではｘ軸正方向にＦ、最後の時間ｄｔではｘ軸負方向にＦであるから、平均して働く摩擦力は、ｘ軸負方向に、（Ｆ−Ｆ＋Ｆ）／３＝Ｆ／３  となり、図５（Ａ）を参照して考察した結果とは異なってしまう。
【００１２】
このことは時間刻みｄｔを十分小さくしても、摺動距離を長くしても起こりうることであり、防止策としては、粒子の大きさを十分大きくとる、もしくは十分な時間的、空間的なアンサンブルをとることなどが考えられる。
【００１３】
しかしながら、粒子を大きくすると複雑な形状が再現できない、時間アンサンブルを十分とるとフィードバックが必要となり、計算に余分な時間がかかる、空間アンサンブルが十分とれない場合は用いることができない、などの問題点があり実用的ではない。
【００１４】
従来の粒子法における摩擦のモデルでは滑り摩擦によるエネルギが正しく評価できず、滑り摩擦と転がり摩擦とを区別して評価する必要がある。しかし現実の摩擦は、滑り摩擦と転がり摩擦とが混在して生じているため、そこに人為的な仮定を含むことは望ましくない。
【００１５】
【特許文献１】特開２００６−２８５８６６号公報
【発明の開示】
【発明が解決しようとする課題】
【００１６】
本発明の目的は、高精度なシミュレーションを実現するシミュレーション方法を提供することである。
【００１７】
また、高精度なシミュレーションを実現するシミュレーション方法を実行させるプログラムを提供することである。
【課題を解決するための手段】
【００１８】
本発明の一観点によれば、物体を粒子の集まりで表現し、ある時間刻み幅で、各粒子の挙動を数値計算することにより、時間軸上で離散的に分布する時刻ごとの該物体の状態を解析する粒子法におけるシミュレーション方法であって、（ａ）第１の時刻において、第１の物体を構成する第１の粒子が、第２の物体の表面から離れており、該第１の時刻の次の第２の時刻において、前記第１の粒子が前記第２の物体の表面よりも深い場所に位置するとき、該第２の時刻における前記第１の粒子の位置を、凝着点として記憶する工程と、（ｂ）前記第２の時刻の次の第３の時刻における前記第１の粒子の位置を数値計算により求める工程と、（ｃ）前記第３の時刻における前記第１の粒子の位置と、前記凝着点の位置とに基づいて、前記第１の粒子に印加される摩擦力を算出する工程と、（ｄ）前記摩擦力を加味して、前記第３の時刻の次の第４の時刻における前記第１の粒子の位置を数値計算により求める工程とを有するシミュレーション方法が提供される。
【００１９】
また、本発明の他の観点によると、物体を粒子の集まりで表現し、ある時間刻み幅で、各粒子の挙動を数値計算することにより、時間軸上で離散的に分布する時刻ごとの該物体の状態を解析する粒子法におけるシミュレーションを行うためにコンピュータを、第１の物体を構成する第１の粒子の初期状態を入力する手段、第１の時刻において、前記第１の粒子が、第２の物体の表面から離れており、該第１の時刻の次の第２の時刻において、前記第１の粒子が前記第２の物体の表面よりも深い場所に位置するとき、該第２の時刻における前記第１の粒子の位置を、凝着点として画定する手段、前記第２の時刻の次の第３の時刻における前記第１の粒子の位置を決定する手段、前記第３の時刻における、前記第１の粒子に作用するポテンシャル力を得る手段、前記第３の時刻における前記第１の粒子の位置と、前記凝着点の位置とに基づいて、前記第３の時刻に、前記第１の粒子に印加される摩擦力を得る手段、前記第３の時刻におけるポテンシャル力と摩擦力とを加えて、前記第３の時刻の次の第４の時刻における前記第１の粒子の位置を決定する手段、前記第４の時刻における前記第１の粒子を表示する表示手段として機能させるためのプログラムが提供される。
【発明の効果】
【００２０】
本発明によれば、高精度なシミュレーションを実現するシミュレーション方法を提供することができる。
【００２１】
また、高精度なシミュレーションを実現するシミュレーション方法を実行させるプログラムを提供することができる。
【発明を実施するための最良の形態】
【００２２】
以下、物体を粒子の集まりで表現し、ある時間刻み幅Δｔで、各粒子の挙動を数値計算することにより、時間軸上で離散的に分布する時刻ごとの物体の状態を解析する粒子法におけるシミュレーション方法、及びそのシミュレーション方法をコンピュータで実行するプログラムの実施例について説明する。
【００２３】
図１を参照して、実施例によるシミュレーション方法に導入される摩擦力について説明する。
【００２４】
実施例においては、２つの物体ＡとＢとの間に働く摩擦力を、「摩擦の凝着説」に基づき以下のように導入する。「摩擦の凝着説」とは、摩擦力が、接触する２つの物体の凝着部をせん断するのに必要な力によって与えられる、と考える考え方である。
【００２５】
粒子ｋは物体Ａを構成する粒子の１つであり、粒子ｉ及びｊは物体Ｂを構成する粒子のうち、隣り合う２つの粒子である。粒子ｉの中心と粒子ｊの中心とを結んだ線分（線分ｉｊ）は、物体Ｂの表面の一部を構成する。
【００２６】
実施例によるシミュレーション方法においては、時刻ｔ１において物体Ｂの表面から離れていた物体Ａの粒子ｋが、時刻ｔ１から時間Δｔ後の時刻ｔ２において、物体Ｂの表面より深い場所に位置するとき、その位置を凝着点として設定、記憶する。
【００２７】
たとえば物体Ａと物体Ｂとの接触の際、時刻ｔ２において、物体Ａを構成する粒子ｋが、物体Ｂを構成する粒子ｉ、ｊ間の位置ｇに進入していた場合、位置ｇが凝着点として設定、記憶される。
【００２８】
粒子ｋは、時間の経過とともに、凝着点ｇから移動しうる。粒子ｋがある時刻ｔｎに物体Ｂから受ける垂直抗力Ｎは、当該時刻ｔｎにおいて、粒子ｋが物体Ｂに進入している深さ（粒子ｋの位置から物体Ｂの表面（図では直線ｉｊ）までの距離Δｒ）に比例する大きさを有し、物体Ｂの表面（直線ｉｊ）に垂直に、粒子ｋの進入深さ方向とは逆方向に作用すると考える。
【００２９】
そして粒子ｋが凝着点ｇから移動したとき、前記時刻ｔｎにおいて、粒子ｋは、凝着点ｇに向かう、物体Ｂの表面（直線ｉｊ）に平行な方向に摩擦力Ｆｋを受けると考える。摩擦力Ｆｋの大きさは、摩擦係数μを用い、μＮｍａｘで与えられるものとする。ここでＮｍａｘは、粒子ｋが物体Ｂに凝着した時刻ｔ２から時刻ｔｎまで（時刻ｔ２及び時刻ｔｎを含む。）に受けた垂直抗力Ｎの最大値である。粒子ｋに作用する摩擦力Ｆｋの大きさをμＮｍａｘとするのは、粒子ｋが、接触部で受けた垂直抗力に相当するエネルギで物体Ｂに凝着すると考えるためである。なお、μは、たとえば０．１である。
【００３０】
物体Ａを構成する粒子ｋが物体Ｂから受ける摩擦力Ｆｋの反作用として、物体Ｂを構成する粒子ｉ、ｊは、それぞれ物体Ａから摩擦力Ｆｉ、Ｆｊを受ける。粒子ｉ、ｊが受ける摩擦力Ｆｉ、Ｆｊの大きさを、それぞれ以下の式（１）、（２）で定める。
2010044710.tif 000003
・・（１）
2010044710.tif 000004
・・（２）
ここでＬｉは粒子ｋと粒子ｉの物体Ｂ表面（直線ｉｊ）方向に沿う距離を示し、Ｌｊは粒子ｋと粒子ｊの物体Ｂ表面（直線ｉｊ）方向に沿う距離を示す。摩擦力Ｆｉ、Ｆｊの作用する向きは、摩擦力Ｆｋの作用する向きと逆向きである。また、摩擦力Ｆｉの大きさと摩擦力Ｆｊの大きさとの和は、摩擦力Ｆｋの大きさに等しい。
【００３１】
なお、粒子ｋに作用する垂直抗力Ｎの反作用として、粒子ｉに力Ｎｉ、粒子ｊに力Ｎｊが働く。力Ｎｉ及び力Ｎｊの大きさはそれぞれ  Ｎｉ＝{Ｌｊ／（Ｌｉ＋Ｌｊ）}＊Ｎ、Ｎｊ＝{Ｌｉ／（Ｌｉ＋Ｌｊ）}＊Ｎ  で表される。力Ｎｉ及び力Ｎｊの作用する向きは、力Ｎの作用する向きと反対方向である。
【００３２】
次図を参照して詳述する実施例によるシミュレーション方法の一部概略は以下の通りである。
【００３３】
たとえば上記時刻ｔ２において凝着点ｇが画定されたとする。時刻ｔ２では、粒子ｋは凝着点ｇから移動していないため、粒子ｋに摩擦力は働かない。そこで粒子ｋに作用するたとえば粒子間ポテンシャル力から、時刻ｔ２から時間Δｔ後の時刻ｔ３における粒子ｋの位置を数値計算する。
【００３４】
時刻ｔ３において、粒子ｋが凝着点ｇから移動している場合、時刻ｔ３における粒子ｋには、たとえば粒子間ポテンシャル力の他、物体Ｂからの摩擦力が作用する。そこで実施例によるシミュレーション方法においては、この摩擦力を算出、加味して、時刻ｔ３から時間Δｔ後の時刻ｔ４における粒子ｋの位置を数値計算する。
【００３５】
なお、物体が移動すると凝着点も移動する。これを説明するために、カットオフ距離ｒｃを導入し、粒子ｋと凝着点ｇとの間の、物体Ｂ（直線ｉｊ）方向に沿う距離ｒｇｋがカットオフ距離ｒｃより離れると、凝着点ｇは消滅すると考える。カットオフ距離ｒｃは、滑り摩擦の場合、どのように設定してもよいが、転がり摩擦の場合、あまり長い距離に設定すると、シミュレーションにおいて、物体が不自然な挙動を示すことがある。カットオフ距離ｒｃは、接触される物体表面の原子間距離（図１においては粒子ｉ、ｊ間の距離）の１／２程度とすればよいであろう。
【００３６】
図２は、実施例によるシミュレーション方法を示すフローチャートである。
【００３７】
実施例によるシミュレーション方法では、まずステップＳ１０１において、粒子の位置や速度など、シミュレーションを行うに当たっての初期状態を決定する。
【００３８】
次に、ステップＳ１０２において、粒子間のポテンシャルから全ての粒子に働く力を導出する。物体を構成する粒子同士は、互いに粒子間のポテンシャルにより結びついている。ステップＳ１０２ではポテンシャルを基に、全ての粒子につき、粒子ごとに働く力を算出する。
【００３９】
粒子間のポテンシャルは、たとえばバネやＬｅｎｎａｒｄ−Ｊｏｎｅｓポテンシャルで表現される。どのようなポテンシャルを用いるかについては、ステップＳ１０１またはそれに先立って、シミュレーション実行者が任意に設定することができる。
【００４０】
ステップＳ１０３では、物体Ｂを構成する任意の粒子ｉ、ｊ間に、物体Ａを構成する任意の粒子ｋが接触しているか否か（任意の粒子ｋにつき物体Ｂ内部に進入しているかか否か）を判定する。
【００４１】
全ての粒子について接触していない場合には、物体Ａ及びＢを構成する粒子には、ポテンシャル場からの力のみが作用するとして、摩擦力の導入は行わず、ステップＳ１１３に進み、物体Ａ及びＢを構成する全ての粒子につき、運動方程式を解いて、時間Δｔ後の粒子の位置を数値計算で求め、粒子を移動させる。
【００４２】
ステップＳ１０３において、ある粒子ｉ、ｊ間に、ある粒子ｋが接触している場合、ステップＳ１０４に進み、その粒子ｋにつき、既に凝着点ｇが存在しているか否かを判定する。凝着点ｇが存在している場合には、ステップＳ１０６において、凝着点ｇと粒子ｋの物体Ｂ表面（直線ｉｊ）方向に沿う距離ｒｇｋとカットオフ距離ｒｃとの大小関係を判定する。凝着点ｇが不存在の場合には、ステップＳ１０５において凝着点ｇを決定した後ステップＳ１０６に進み、距離ｒｇｋと距離ｒｃとの大小関係を判定する。
【００４３】
カットオフ距離ｒｃの値については、ステップＳ１０１またはそれに先立って、シミュレーション実行者が任意に設定することができる。
【００４４】
ステップＳ１０６において、距離ｒｇｋが距離ｒｃよりも大きい（ｒｇｋ＞ｒｃ）場合、ステップＳ１１２において凝着点ｇを消去する。その粒子ｋについては摩擦力を導入することなく、ポテンシャルによる力のみが作用しているとして、ステップＳ１１３に進む。
【００４５】
ステップＳ１０６において、距離ｒｇｋが距離ｒｃよりも小さい（ｒｇｋ＜ｒｃ）場合、ステップＳ１０７において垂直抗力Ｎを導出し、ステップＳ１０８において、導出したＮとＮｍａｘとの大小を比較する。
【００４６】
Ｎ＜Ｎｍａｘの場合、ステップＳ１１０に進み、Ｎｍａｘを用いて粒子ｋ、ｉ、ｊに摩擦による力を導入する。Ｎ＞Ｎｍａｘの場合、ステップＳ１０９においてＮをＮｍａｘに置き換え（Ｎ＞ＮｍａｘであるＮを新たなＮｍａｘとし）た後、ステップＳ１１０に進み、Ｎｍａｘを用いて粒子ｋ、ｉ、ｊに摩擦による力を導入する。
【００４７】
ステップＳ１１０で粒子ｋ、ｉ、ｊに導入される摩擦力Ｆｋ、Ｆｉ、Ｆｊは、
2010044710.tif 000005
・・（３）
2010044710.tif 000006
・・（４）
2010044710.tif 000007
・・（５）
と表される。ここで、
2010044710.tif 000008
は、粒子ｋより物体Ｂ表面（直線ｉｊ）に下ろした垂線の足から、凝着点ｇより物体Ｂ表面（直線ｉｊ）に下ろした垂線の足に向かうベクトルを指す。
【００４８】
ステップＳ１１０で導入された摩擦力は、ステップＳ１１１において、ステップＳ１０２で導出されたポテンシャル場からの力に加えられ、ステップＳ１１３に進む。
【００４９】
ステップＳ１０３（粒子ｉ、ｊ間に粒子ｋが接触しているか否か（粒子ｋが物体Ｂ内部に進入しているか否か）の判定工程）で接触していると判定された全ての粒子ｉ、ｊ、ｋについて、ステップＳ１０４〜Ｓ１１２の工程は実施される。また、ステップＳ１０３で接触していないと判定された全ての粒子は、ステップＳ１０４〜Ｓ１１２の工程を経ることなく、ステップＳ１１３に進む。
【００５０】
そしてステップＳ１１３において、物体Ａ及びＢを構成する全ての粒子につき、運動方程式を解いて、時間Δｔ後の粒子の位置を数値計算で求め、粒子を移動させる。
【００５１】
ステップＳ１１３で全ての粒子についての移動が終わったら、ステップＳ１１４にでシミュレーションを終了するか否かを決定する。終了を選択した場合は、ステップＳ１１５に進み、シミュレーションは完了する。終了しない場合は再びステップＳ１０２に戻り、ステップＳ１０２〜Ｓ１１４の工程を繰り返す。
【００５２】
図３（Ａ）及び（Ｂ）を参照して実施例によるシミュレーション方法の効果について説明する。
【００５３】
図３（Ａ）に、実施例によるシミュレーション方法を用いて行ったシミュレーションにおける、シミュレーション対象（歯車）の一部を示す。
【００５４】
歯１とピンとの間の摩擦係数をμ１、歯２とピンとの間の摩擦係数をμ２とし、歯２を固定する条件で、歯１と歯２とをピンを介して噛み合わせながら運動させるシミュレーションを行った。
【００５５】
歯１が図中の矢印に示す軌跡を辿るとき、摩擦係数μ１による摩擦力と、摩擦係数μ２による摩擦力とはピンに対して逆向きに作用し、一方がピンの回転を促し、他方がピンの回転を妨げる方向に働く。
【００５６】
図３（Ｂ）は、これらの摩擦力がピンに及ぼすトルクの時間履歴を示すグラフである。グラフの横軸は時間を、縦軸はトルクを、それぞれ任意単位で示す。グラフにおいて、正のトルクはピンを反時計方向に回そうとするトルク、負のトルクはピンを時計方向に回そうとするトルクを表している。
【００５７】
図３（Ｂ）より、歯１の摩擦力によるトルクと、歯２の摩擦力によるトルクとは相互に釣り合っており、結果としてピンは回転せず滑ることがわかる。
【００５８】
図３（Ｂ）には、従来の摩擦モデル、すなわち図５（Ａ）及び（Ｂ）を参照して説明した、接触部の粒子の相対速度の方向と反対方向に摩擦力が発生すると考えるモデルを使ったシミュレーションにおいて発生するトルクを併せて示した。
【００５９】
従来の摩擦モデルを用いた場合、トルクがほとんど発生していないという結果が得られている。ピンは回転せずに滑るという結果が得られる点においては、実施例によるシミュレーション方法を用いた場合と同様であるが、従来の摩擦モデルを使用した場合には、滑り摩擦時に発生するエネルギを正しく評価することはできないことが理解される。
【００６０】
これに対して、実施例によるシミュレーション方法は、摩擦のメカニズムに基づいているため、これを用いてシミュレーションを行うと、転がり摩擦、滑り摩擦を区別することなく、摩擦を自然に再現することができ、たとえば滑り摩擦時に発生するエネルギを正しく評価することが可能となる。このため、高精度のシミュレーションを実現することができる。
【００６１】
実施例によるシミュレーション方法は、プログラムにより、コンピュータで実行することができる。
【００６２】
図４は、実施例によるシミュレーション方法を用いてシミュレーションを行うシミュレーション装置のシステム構成図である。本図に構成を示したシミュレーション装置を使用して、図２にフローチャートで示したシミュレーション方法を実施することができる。
【００６３】
まず、キーボードなどの入力装置から、図２のＳ１０１に対応して、粒子位置や粒子の速度などの粒子の初期状態が入力される。また、これと同時に、またはこれに先立って、キーボードなどの入力装置から、粒子間のポテンシャルを表現する関数等、及び、カットオフ距離ｒｃの値を設定することができる。
【００６４】
中央処理装置は、メインメモリ中の制御プログラムの指令を受け、メインメモリ中の粒子の移動位置決定プログラムを実行する。
【００６５】
粒子の移動位置決定プログラムは、粒子間のポテンシャルから全ての粒子について、作用するポテンシャル力を導出する部分（図２のＳ１０２で示される部分）、２つの物体が接触しているとき、凝着点を画定できる場合には画定する部分（図２のＳ１０３〜Ｓ１０６、及びＳ１１２で示される部分）、凝着点が画定された粒子、及びその間に凝着点が画定された相隣る２つの粒子については摩擦力を計算し、これをポテンシャル力に加算する部分（図２のＳ１０７〜Ｓ１１１で示される部分）、全ての粒子につき運動方程式を解いて移動位置を決定する部分（図２のＳ１１３で示される部分）に大別される。
【００６６】
入力装置から入力された情報に基づいて、図２を参照して説明したように、全ての粒子について、ポテンシャル場から受ける力が算出される。
【００６７】
次に、粒子の凝着点を画定できる場合には、凝着点の位置が計算され、決定される。
【００６８】
そして計算された凝着点の位置に基づいて、図２のＳ１１０に示される摩擦力を算出する。摩擦力は他の物体の粒子と接触している全ての粒子について算出されることになる。
【００６９】
摩擦力はポテンシャル場からの力に加算され、各々の粒子について作用する力が求められる。
【００７０】
全ての粒子につき運動方程式を解いて、粒子の移動位置を決定する。
【００７１】
粒子の移動位置が決定されたら、たとえば２次元的、または３次元的に定められた座標系において、決定された移動位置に粒子を表示する。表示結果は、出力装置、たとえばディスプレイに表示される。
【００７２】
以上実施例に沿って本発明を説明したが、本発明はこれらに制限されるものではない。種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。
【産業上の利用可能性】
【００７３】
摩擦を伴う種々の現象のシミュレーション、特に機構や弾性の動解析に好適に利用することができる。機構や弾性の動解析における、摩擦によるロス、及び発熱エネルギの見積もりに応用することができる。
【図面の簡単な説明】
【００７４】
【図１】実施例によるシミュレーション方法に導入される摩擦力について説明するための図である。
【図２】実施例によるシミュレーション方法を示すフローチャートである。
【図３】（Ａ）及び（Ｂ）は、実施例によるシミュレーション方法の効果について説明するための図である。
【図４】実施例によるシミュレーション方法を用いてシミュレーションを行うシミュレーション装置のシステム構成図である。
【図５】（Ａ）及び（Ｂ）は、摩擦力を導入する従来モデルにおける問題点を説明するための図である。
【符号の説明】
【００７５】
Ｓ１０１〜Ｓ１１５  ステップ
【図１】
2010044710.tif 000009
【図２】
2010044710.tif 000010
【図４】
2010044710.tif 000011
【図５】
2010044710.tif 000012
【図３】
2010044710.tif 000013