(57)【要約】【課題】  接触問題を扱う。【解決手段】  （ａ）第１の物体及び第２の物体を、それぞれ複数の粒子の集合で定義する。（ｂ）第１の物体を構成する粒子と第２の物体を構成する粒子との間の第１のポテンシャルを定める。（ｃ）第１の物体を構成する粒子と第２の物体を構成する粒子とが第１の距離だけ離れている場合に働く引力を、第１のポテンシャルにより求められる、第１の距離だけ離れた第１の物体の構成粒子と第２の物体の構成粒子との間に働く引力より小さくするという条件の下で、第１の物体を構成する粒子と第２の物体を構成する粒子との間に作用する力を加味して、第１及び第２の物体を構成する粒子の挙動を計算する。【選択図】          図４
2011013888.tif 000002
【特許請求の範囲】
【請求項１】
（ａ）第１の物体及び第２の物体を、それぞれ複数の粒子の集合で定義する工程と、
（ｂ）前記第１の物体を構成する粒子と前記第２の物体を構成する粒子との間の第１のポテンシャルを定める工程と、
（ｃ）前記第１の物体を構成する粒子と前記第２の物体を構成する粒子とが第１の距離だけ離れている場合に働く引力を、前記第１のポテンシャルにより求められる、前記第１の距離だけ離れた前記第１の物体の構成粒子と前記第２の物体の構成粒子との間に働く引力より小さくするという条件の下で、前記第１の物体を構成する粒子と前記第２の物体を構成する粒子との間に作用する力を加味して、前記第１及び第２の物体を構成する粒子の挙動を計算する工程と
を有するシミュレーション方法。
【請求項２】
前記工程（ｂ）において、前記第１のポテンシャルを、前記第１の物体を構成する物質のポテンシャルと、前記第２の物体を構成する物質のポテンシャルとに基づいて定める請求項１に記載のシミュレーション方法。
【請求項３】
前記第１の物体の表面を構成する粒子と前記第２の物体の表面を構成する粒子との相互作用を考慮しないという条件で、
前記工程（ｃ）において、前記第１の物体の表面から２層めを構成する粒子が、前記第２の物体を構成する粒子との関係で有する第２のポテンシャルによって、前記第１の物体を構成する粒子と前記第２の物体を構成する粒子との間に作用する力を加味する請求項１または２に記載のシミュレーション方法。
【請求項４】
前記第２のポテンシャルは、前記第２の物体を構成する粒子が近づいてきたときに斥力を受けはじめる距離であって、前記第１の物体の表面に垂直な方向に沿う第２の距離が、前記第２の物体を構成する粒子が近づいてきたときに斥力を受けはじめる距離であって、前記第１の物体の表面に平行な方向に沿う第３の距離よりも短くなるようなポテンシャルである請求項３に記載のシミュレーション方法。
【請求項５】
前記工程（ｃ）において、前記第１の物体を構成する粒子と前記第２の物体を構成する粒子との間に作用する力を加味するとともに、前記第２のポテンシャルによって前記第１の物体を構成する粒子が受ける作用を加味する請求項１〜４のいずれか１項に記載のシミュレーション方法。
【請求項６】
（ａ）第１の物体及び第２の物体を、それぞれ複数の粒子の集合で定義する工程と、
（ｂ）前記第１の物体を構成する粒子と前記第２の物体を構成する粒子との間に働く引力をゼロとする条件の下で、前記第１の物体を構成する粒子と前記第２の物体を構成する粒子との間に作用する力を加味して、前記第１及び第２の物体を構成する粒子の挙動を計算する工程と
を有するシミュレーション方法。
【請求項７】
前記第１の物体の表面を構成する粒子と前記第２の物体の表面を構成する粒子との相互作用を考慮しないという条件で、
前記工程（ｂ）において、前記第１の物体の表面から２層めを構成する粒子が、前記第２の物体を構成する粒子との関係で有するポテンシャルによって、前記第１の物体を構成する粒子と前記第２の物体を構成する粒子との間に作用する力を加味する請求項６に記載のシミュレーション方法。
【請求項８】
前記ポテンシャルは、前記第２の物体を構成する粒子が近づいてきたときに斥力を受けはじめる距離であって、前記第１の物体の表面に垂直な方向に沿う第１の距離が、前記第２の物体を構成する粒子が近づいてきたときに斥力を受けはじめる距離であって、前記第１の物体の表面に平行な方向に沿う第２の距離よりも短くなるようなポテンシャルである請求項７に記載のシミュレーション方法。
【請求項９】
前記工程（ｂ）において、前記第１の物体を構成する粒子と前記第２の物体を構成する粒子との間に作用する力を加味するとともに、前記ポテンシャルによって前記第１の物体を構成する粒子が受ける作用を加味する請求項６〜８のいずれか１項に記載のシミュレーション方法。
【請求項１０】
物体を粒子の集まりで表現し、粒子の挙動を数値計算する粒子法におけるシミュレーションを行うためにコンピュータを、
第１の物体を構成する粒子と第２の物体を構成する粒子の初期状態を入力する手段、
前記第１の物体を構成する粒子と前記第２の物体を構成する粒子とが第１の距離だけ離れている場合に働く引力を、前記第１の物体を構成する粒子と前記第２の物体を構成する粒子との間の第１のポテンシャルにより求められる、前記第１の距離だけ離れた前記第１の物体の構成粒子と前記第２の物体の構成粒子との間に働く引力より小さくするという条件の下で、前記第１の物体を構成する粒子と前記第２の物体を構成する粒子との間に作用する力を加味して、前記第１及び第２の物体を構成する粒子の挙動を計算する手段、
前記計算する手段によって粒子の挙動に関して表示を行う表示手段
として機能させるためのプログラム。
【請求項１１】
前記第１のポテンシャルは、前記第１の物体を構成する物質のポテンシャルと、前記第２の物体を構成する物質のポテンシャルとに基づいて定められる請求項１０に記載のプログラム。
【請求項１２】
前記第１の物体を構成する粒子と前記第２の物体を構成する粒子との間に作用する力を加味して、挙動を計算する手段は、
前記第１の物体の表面を構成する粒子と前記第２の物体の表面を構成する粒子との相互作用を考慮しないという条件で、
前記第１の物体の表面から２層めを構成する粒子が、前記第２の物体を構成する粒子との関係で有する第２のポテンシャルによって、前記第１の物体を構成する粒子と前記第２の物体を構成する粒子との間に作用する力を加味する請求項１０または１１に記載のプログラム。
【請求項１３】
前記第２のポテンシャルは、前記第２の物体を構成する粒子が近づいてきたときに斥力を受けはじめる距離であって、前記第１の物体の表面に垂直な方向に沿う第２の距離が、前記第２の物体を構成する粒子が近づいてきたときに斥力を受けはじめる距離であって、前記第１の物体の表面に平行な方向に沿う第３の距離よりも短くなるようなポテンシャルである請求項１２に記載のプログラム。
【請求項１４】
前記第１の物体を構成する粒子と前記第２の物体を構成する粒子との間に作用する力を加味して、挙動を計算する手段は、前記第１の物体を構成する粒子と前記第２の物体を構成する粒子との間に作用する力を加味するとともに、前記第２のポテンシャルによって前記第１の物体を構成する粒子が受ける作用を加味する請求項１０〜１３のいずれか１項に記載のプログラム。
【請求項１５】
物体を粒子の集まりで表現し、粒子の挙動を数値計算する粒子法におけるシミュレーションを行うためにコンピュータを、
第１の物体を構成する粒子と第２の物体を構成する粒子の初期状態を入力する手段、
前記第１の物体を構成する粒子と前記第２の物体を構成する粒子との間に働く引力をゼロとする条件の下で、前記第１の物体を構成する粒子と前記第２の物体を構成する粒子との間に作用する力を加味して、前記第１及び第２の物体を構成する粒子の挙動を計算する手段、
前記計算する手段によって粒子の挙動に関して表示を行う表示手段
として機能させるためのプログラム。
【請求項１６】
前記第１の物体を構成する粒子と前記第２の物体を構成する粒子との間に作用する力を加味して、挙動を計算する手段は、
前記第１の物体の表面を構成する粒子と前記第２の物体の表面を構成する粒子との相互作用を考慮しないという条件で、
前記第１の物体の表面から２層めを構成する粒子が、前記第２の物体を構成する粒子との関係で有するポテンシャルによって、前記第１の物体を構成する粒子と前記第２の物体を構成する粒子との間に作用する力を加味する請求項１５に記載のプログラム。
【請求項１７】
前記ポテンシャルは、前記第２の物体を構成する粒子が近づいてきたときに斥力を受けはじめる距離であって、前記第１の物体の表面に垂直な方向に沿う第１の距離が、前記第２の物体を構成する粒子が近づいてきたときに斥力を受けはじめる距離であって、前記第１の物体の表面に平行な方向に沿う第２の距離よりも短くなるようなポテンシャルである請求項１６に記載のプログラム。
【請求項１８】
前記第１の物体を構成する粒子と前記第２の物体を構成する粒子との間に作用する力を加味して、挙動を計算する手段は、前記第１の物体を構成する粒子と前記第２の物体を構成する粒子との間に作用する力を加味するとともに、前記ポテンシャルによって前記第１の物体を構成する粒子が受ける作用を加味する請求項１５〜１７のいずれか１項に記載のプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、シミュレーション方法及びプログラムに関する。特に、物体を粒子の集まりで表現し、粒子の挙動を数値計算する粒子法を用いたシミュレーション方法、及びそのシミュレーション方法をコンピュータで実行するプログラムに関する。
【背景技術】
【０００２】
解析対象を単純な形状をもつ小さな要素に分割して、各要素に関する計算を行い、これを組み合わせて、全体の解を求める有限要素法(Finite Elemental Method; FEM)が知られている。有限要素法は、たとえば構造解析の分野において広く利用される。離れている２物体が接触することによって発生する応力、変形を有限要素法を用いて解析するための接触解析方法に関する発明が開示されている（たとえば、特許文献１参照）。
【０００３】
有限要素法を利用する接触解析には、接触のアルゴリズムが導入される。このため、接触面の探索や接触判定に時間がかかる上、これをコンピュータで行う際にはバグの原因ともなっていた。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特許２６９２６６８号公報
【発明の概要】
【発明が解決しようとする課題】
【０００５】
本発明の目的は、接触問題を扱うことのできるシミュレーション方法及びプログラムを提供することである。
【０００６】
また、高い正確性で接触問題を扱うことのできるシミュレーション方法及びプログラムを提供することである。
【課題を解決するための手段】
【０００７】
本発明の一観点によれば、（ａ）第１の物体及び第２の物体を、それぞれ複数の粒子の集合で定義する工程と、（ｂ）前記第１の物体を構成する粒子と前記第２の物体を構成する粒子との間の第１のポテンシャルを定める工程と、（ｃ）前記第１の物体を構成する粒子と前記第２の物体を構成する粒子とが第１の距離だけ離れている場合に働く引力を、前記第１のポテンシャルにより求められる、前記第１の距離だけ離れた前記第１の物体の構成粒子と前記第２の物体の構成粒子との間に働く引力より小さくするという条件の下で、前記第１の物体を構成する粒子と前記第２の物体を構成する粒子との間に作用する力を加味して、前記第１及び第２の物体を構成する粒子の挙動を計算する工程とを有するシミュレーション方法が提供される。
【０００８】
また、本発明の他の観点によれば、（ａ）第１の物体及び第２の物体を、それぞれ複数の粒子の集合で定義する工程と、（ｂ）前記第１の物体を構成する粒子と前記第２の物体を構成する粒子との間に働く引力をゼロとする条件の下で、前記第１の物体を構成する粒子と前記第２の物体を構成する粒子との間に作用する力を加味して、前記第１及び第２の物体を構成する粒子の挙動を計算する工程とを有するシミュレーション方法が提供される。
【０００９】
更に、本発明の他の観点によれば、物体を粒子の集まりで表現し、粒子の挙動を数値計算する粒子法におけるシミュレーションを行うためにコンピュータを、第１の物体を構成する粒子と第２の物体を構成する粒子の初期状態を入力する手段、前記第１の物体を構成する粒子と前記第２の物体を構成する粒子とが第１の距離だけ離れている場合に働く引力を、前記第１の物体を構成する粒子と前記第２の物体を構成する粒子との間の第１のポテンシャルにより求められる、前記第１の距離だけ離れた前記第１の物体の構成粒子と前記第２の物体の構成粒子との間に働く引力より小さくするという条件の下で、前記第１の物体を構成する粒子と前記第２の物体を構成する粒子との間に作用する力を加味して、前記第１及び第２の物体を構成する粒子の挙動を計算する手段、前記計算する手段によって粒子の挙動に関して表示を行う表示手段として機能させるためのプログラムが提供される。
【００１０】
また、本発明の他の観点によれば、物体を粒子の集まりで表現し、粒子の挙動を数値計算する粒子法におけるシミュレーションを行うためにコンピュータを、第１の物体を構成する粒子と第２の物体を構成する粒子の初期状態を入力する手段、前記第１の物体を構成する粒子と前記第２の物体を構成する粒子との間に働く引力をゼロとする条件の下で、前記第１の物体を構成する粒子と前記第２の物体を構成する粒子との間に作用する力を加味して、前記第１及び第２の物体を構成する粒子の挙動を計算する手段、前記計算する手段によって粒子の挙動に関して表示を行う表示手段として機能させるためのプログラムが提供される。
【発明の効果】
【００１１】
本発明によれば、接触問題を扱うことの可能なシミュレーション方法及びプログラムを提供することができる。
【００１２】
また、高い正確性で接触問題を扱うことの可能なシミュレーション方法及びプログラムを提供することができる。
【図面の簡単な説明】
【００１３】
【図１】（Ａ）は、異なる２つの物体Ａ、Ｂ、及びそれらを構成する原子ａ〜ｅを示す概略図であり、（Ｂ）は、物体Ａを構成する原子のポテンシャルを示す概略的なグラフであり、（Ｃ）は、物体Ｂを構成する原子との関係における、物体Ａの構成原子のポテンシャルを示す概略的なグラフである。
【図２】角度依存型ポテンシャルの例を２次元的に示す。
【図３】物体Ａ及びＢを構成する原子の位置関係を３次元的に示す概略図である。
【図４】実施例によるシミュレーション方法を示すフローチャートである。
【図５】実施例によるシミュレーション方法を用いてシミュレーションを行うシミュレーション装置のシステム構成図である。
【図６】（Ａ）及び（Ｂ）は、球の衝突シミュレーションについて説明するための図である。
【発明を実施するための形態】
【００１４】
分子動力学法や粒子法において接触問題を扱う手法の実施例について説明する。
【００１５】
図１（Ａ）は、異なる２つの物体Ａ、Ｂ、及びそれらを構成する原子ａ〜ｅを示す概略図である。原子ａ及びｂは、物体Ａの表面を構成する原子であり、原子ａの存在位置と原子ｂの存在位置とを結ぶ線分は物体Ａの表面の一部を画定する。同様に、原子ｄ及びｅは、物体Ｂの表面を構成する原子であり、原子ｄの存在位置と原子ｅの存在位置とを結ぶ線分は物体Ｂの表面の一部を画定する。また、原子ｃは、物体Ａの最表面から２層めを構成する原子である。
【００１６】
図１（Ｂ）に、物体Ａを構成する原子のポテンシャルを示す。グラフの縦軸は、ポテンシャルφを示し、横軸は原子間距離ｒを示す。ポテンシャルφは、原子間距離が安定距離ｒ０のときに最小値をとる。物体Ａを構成する原子間においては、原子間距離がｒ０より小さいときに相互に斥力が働き、大きいときに相互に引力が働く。同様に、物体Ｂを構成する原子間においても、原子間距離が安定距離より小さいときには相互に斥力が働き、大きいときには相互に引力が働く。
【００１７】
図１（Ｃ）は、物体Ｂを構成する原子との関係における、物体Ａの構成原子のポテンシャルを示す概略的なグラフである。安定距離ｒ０より大きいｒの範囲ではポテンシャルφの値が一定となっている点で、図１（Ｂ）に示すグラフと異なる。すなわち、物体Ａを構成する原子は、物体Ｂを構成する原子に引力を及ぼさない。原子間距離がｒ０より小さい範囲において斥力を及ぼすのみである。
【００１８】
同様に、物体Ｂを構成する原子は、物体Ａを構成する原子に引力を及ぼさず、原子間距離が安定距離より小さい範囲において斥力を及ぼす。実施例においては、このように作用する力を導入する。
【００１９】
再び、図１（Ａ）を参照する。また、実施例においては、物体の表面を、その物体を構成する原子のポテンシャルを用いて表現する。この場合、物体の最表面を構成する原子のポテンシャルを用いると、実際の物体表面とポテンシャル面との懸隔が大きくなるため、実施例では、物体の最表面から２層めを構成する原子のポテンシャルを使って、物体の表面を表現することとし、たとえば図１（Ａ）に示す例においては、物体Ａの表面から２層めを構成する原子（たとえば原子ｃ）と、物体Ｂの表面を構成する原子（たとえば原子ｄ及びｅ）との相互作用を考える。なお、この場合、物体Ａの表面を構成する原子と、物体Ｂの表面を構成する原子とは相互作用を及ぼさないと考える。
【００２０】
原子のポテンシャルを使って物体の表面を表現すると、接触のアルゴリズムを導入する有限要素法の場合と異なり、分子動力学法の場合と同様に、物体の接触解析において接触面を判別しなくてもよい（接触面の判断を行うアルゴリズムを使用しなくてよい）という利点がある。
【００２１】
図１（Ａ）には、物体Ａの第２層め原子ｃのつくる等ポテンシャル面（半径ｒ０の球面）を点線で示した。物体Ｂの表面原子ｄ、ｅは、この等ポテンシャル面の内部にあるとき、原子ｃとの間に相互に斥力を及ぼしあう。しかし、この等ポテンシャル面の一部、たとえば原子ｄ、ｅ間のポテンシャル面を使用して物体Ａの表面の一部を表そうとすると、表された表面の一部と、実際の物体Ａ表面との間のずれ量（図１（Ａ）においては、たとえば矢印で図示した部分）が大きくなる。このためこのような近似で接触問題の解析を行うと、接触位置が不正確になることがある。
【００２２】
そこで本願発明者らは、等方向的ではないポテンシャルを導入することにした。
【００２３】
図２に、本願発明者らが導入する角度依存型ポテンシャルの例を２次元的に示す。本図に示す曲線αが、物体Ａの第２層め原子ｃのつくる等ポテンシャル面を示し、原子ｃと物体Ｂの（表面）原子ｄとは、曲線αの内部にあるとき相互に斥力を及ぼしあうと考える。曲線αは物体Ａ表面を表す線に接している。曲線αの一部、たとえば原子ａと原子ｂとの間にある一部を用いて、物体Ａの表面を表現する。
【００２４】
曲線αは、原子ｃの存在位置からの距離（安定距離）ｒ０´が、次式（１）で表される曲線である。
【００２５】
【数１】
2011013888.tif 000003
【００２６】
ここで、ｒ０は原子ｃと原子ａ及びｂとの間の距離を示す。また、θ０は、原子ａとｂとの間における物体Ａ表面の法線ベクトル（原子ｃから物体Ａ表面に向かうｎベクトル）と、原子ｃから原子ａまたはｂへ向かうベクトルとのなす角を示す。更に、θは、ｎベクトルと原子ｃから原子ｄ（物体Ｂの（表面）原子）へ向かうベクトルとのなす角を示す。
【００２７】
すなわち、本願発明者らが導入する角度依存型ポテンシャルにおいては、原子ｃの安定距離をθの関数とし、それを用いて物体Ａの表面を表現する。角度依存型ポテンシャルを使用することにより、図１（Ａ）に点線で示したような球状の二体ポテンシャルで物体Ａの表面を表現する場合よりも、表面の凹凸が少なくなり、表面のより正確な表現が可能になる。したがって接触問題において、接触位置を高い正確性で把握して解析を行うことができる。
【００２８】
なお、式（１）は、図２において、原子ｃを構成する電子のｐ軌道のうちの２つを２次元的に表す曲線β、曲線γ上の点と原子ｃの中心位置との間の距離を、それぞれｒ１、ｒ２とすると、
【００２９】
【数２】
2011013888.tif 000004
【００３０】
【数３】
2011013888.tif 000005
【００３１】
と表されることから、
【００３２】
【数４】
2011013888.tif 000006
【００３３】
として求められたものである。
【００３４】
また、式（１）の安定距離ｒ０´を用いて、原子ｃから距離ｒだけ離れた位置のポテンシャルφを、たとえば、
【００３５】
【数５】
2011013888.tif 000007
【００３６】
と表すことができる。
【００３７】
これは、原子間の二体ポテンシャルφをバネのポテンシャルを使って、
【００３８】
【数６】
2011013888.tif 000008
【００３９】
と表すとき、ｒ０＝ｒ０´とすると、
【００４０】
【数７】
2011013888.tif 000009
【００４１】
となり、ｒ＝０のときのφ（φｒ＝０）がｒ０´に依存してしまうため、この依存性を解消する目的で、式（７）におけるφｒ＝０を、
【００４２】
【数８】
2011013888.tif 000010
【００４３】
におけるφｒ＝０と等しくする、すなわち、
【００４４】
【数９】
2011013888.tif 000011
【００４５】
として導いた
【００４６】
【数１０】
2011013888.tif 000012
【００４７】
を式（７）に代入して得られた下式（１１）
【００４８】
【数１１】
2011013888.tif 000013
【００４９】
を変形して導出することができる。
【００５０】
図３は、物体Ａ及びＢを構成する原子の位置関係を３次元的に示す概略図である。本図においては、物体Ａの表面を構成する原子をｉｉ、ｊｊ、ｋｋ、第２層めを構成する原子をｊｍと表した。原子ｉｉ、ｊｊ、ｋｋの存在点を結んで得られる三角形は物体Ａの表面の一部を画定する。物体Ａの外部に向かう、原子ｉｉ、ｊｊ、ｋｋの張る平面の法線をｎベクトルと表示した。更に、物体Ａから見た外部原子をｉｍと表した。原子ｉｍは、たとえば物体Ｂの表面を構成する原子である。
【００５１】
物体Ａの表面原子ｉｉ、ｊｊ、ｋｋは、物体Ｂの表面原子ｉｍとの関係では、ポテンシャルを有していない。または、原子間距離で微分してもゼロとなる、フラットな形状のポテンシャルを有するのみである。物体Ｂの表面原子ｉｍも、物体Ａの表面原子との関係において同様である。物体Ａの第２層め原子ｊｍは、外部原子に対して式（５）で表されるポテンシャルを有している。式（５）におけるｒ０´は式（１）で示される。原子ｊｍと外部原子との間に引力は作用せず、斥力のみが作用しうる。本図を参照して、式（５）で表されるポテンシャルにより原子にかかる力を求める。原子にかかる力の向きは、斥力を正とする。
【００５２】
まず、原子ｉｍと原子ｊｍとの間に作用する力Ｆｉｍ−ｊｍを求める。力Ｆｉｍ−ｊｍは、以下の式（１２）で表される。
【００５３】
【数１２】
2011013888.tif 000014
【００５４】
ｒｉｍ−ｊｍは原子ｉｍと原子ｊｍとの間の距離を示す。また、これに矢印を付加した記号で、原子ｉｍの存在位置から原子ｊｍの存在位置に向かうベクトルを表す。以下、同種の記号で同様の意味を表す。
【００５５】
理解の便宜のために、ｘ成分のみを考える。力Ｆｉｍ−ｊｍのｘ成分Ｆｉｍ−ｊｍ（ｘ）は、次のように書ける。
【００５６】
【数１３】
2011013888.tif 000015
【００５７】
2011013888.tif 000016
【００５８】
2011013888.tif 000017
【００５９】
なお、ｘｉｍ−ｊｍは、ｒｉｍ−ｊｍのｘ成分を示す。
【００６０】
ここで、ｒ０´部分の微分は、
【００６１】
【数１４】
2011013888.tif 000018
【００６２】
2011013888.tif 000019
【００６３】
2011013888.tif 000020
【００６４】
2011013888.tif 000021
【００６５】
2011013888.tif 000022
【００６６】
となり、更にｃｏｓθは、ｎベクトルのｘ、ｙ、ｚ成分ｎｘ、ｎｙ、ｎｚを用いて、
【００６７】
【数１５】
2011013888.tif 000023
【００６８】
と書き換えられるので、ｃｏｓθの微分は、
【００６９】
【数１６】
2011013888.tif 000024
【００７０】
2011013888.tif 000025
【００７１】
2011013888.tif 000026
【００７２】
2011013888.tif 000027
【００７３】
2011013888.tif 000028
【００７４】
となる。以上をまとめると、原子ｉｍと原子ｊｍとの間に働く力のｘ成分Ｆｉｍ−ｊｍ（ｘ）は、
【００７５】
【数１７】
2011013888.tif 000029
【００７６】
2011013888.tif 000030
【００７７】
と表すことができる。これをベクトル成分で表示すると、力Ｆｉｍ−ｊｍ（ｘ）は、下式（１８）で示される。
【００７８】
【数１８】
2011013888.tif 000031
【００７９】
同様に、力Ｆｉｍ−ｊｍのｙ成分Ｆｉｍ−ｊｍ（ｙ）、及びｚ成分Ｆｉｍ−ｊｍ（ｚ）を求めることができる。こうして原子ｉｍと原子ｊｍとの間に作用する力Ｆｉｍ−ｊｍが求まる。力Ｆｉｍ−ｊｍとして引力が求まった場合には、これをゼロとする。
【００８０】
次に物体Ａの表面原子に働く力を求めるために、法線ベクトル成分を考える。法線ベクトルは、原子ｉｉと原子ｊｊとを結ぶベクトルと、原子ｉｉと原子ｋｋとを結ぶベクトルの外積から求めているため、ここから原子ｉｉと原子ｊｊとの間に働く力、及び原子ｉｉと原子ｋｋとの間に働く力が求まる。
【００８１】
原子ｉｍと原子ｊｍとの間に作用する力の計算の場合と同様に、ポテンシャルの微分を行う。ｘ成分のみを考えると、原子ｊｍが外部原子との関係で有するポテンシャルによって、原子ｉｉと原子ｊｊとの間に働く力のｘ成分Ｆｉｉ−ｊｊ（ｘ）は、
【００８２】
【数１９】
2011013888.tif 000032
【００８３】
2011013888.tif 000033
【００８４】
2011013888.tif 000034
【００８５】
となる。ｒ０´部分の微分は、
【００８６】
【数２０】
2011013888.tif 000035
【００８７】
2011013888.tif 000036
【００８８】
となり、ｃｏｓθの微分は、
【００８９】
【数２１】
2011013888.tif 000037
【００９０】
と書ける。ここで、
【００９１】
【数２２】
2011013888.tif 000038
【００９２】
とおくと、式（２１）は、
【００９３】
【数２３】
2011013888.tif 000039
【００９４】
2011013888.tif 000040
【００９５】
と表される。対称性に着目して、ｎｘ成分のみを考えると以下のようになる。
【００９６】
【数２４】
2011013888.tif 000041
【００９７】
2011013888.tif 000042
【００９８】
2011013888.tif 000043
【００９９】
2011013888.tif 000044
【０１００】
2011013888.tif 000045
【０１０１】
ここで、
【０１０２】
【数２５】
2011013888.tif 000046
【０１０３】
である。
【０１０４】
Ｎｘ、Ｎｙ、Ｎｚはそれぞれ、
【０１０５】
【数２６】
2011013888.tif 000047
【０１０６】
と表されるので、それらの微分は、
【０１０７】
【数２７】
2011013888.tif 000048
【０１０８】
となる。
【０１０９】
よって原子ｉｉ、原子ｊｊ間に作用する力は、たとえばｘ成分について、
【０１１０】
【数２８】
2011013888.tif 000049
【０１１１】
2011013888.tif 000050
【０１１２】
2011013888.tif 000051
【０１１３】
2011013888.tif 000052
【０１１４】
2011013888.tif 000053
【０１１５】
とまとめることができる。同様に、ｙｉｉ−ｊｊ、ｚｉｉ−ｊｊで微分し、原子ｉｉ、原子ｊｊ間に作用する力のｙ成分Ｆｉｉ−ｊｊ（ｙ）、及びｚ成分Ｆｉｉ−ｊｊ（ｚ）も求められる。また、たとえばｘｉｉ−ｋｋ、ｙｉｉ−ｋｋ、ｚｉｉ−ｋｋで微分することで、原子ｉｉ、原子ｋｋ間に作用する力のｘ成分、ｙ成分、ｚ成分を計算することができる。このようにして、物体Ａの第２層め原子ｊｍが外部原子との関係で有するポテンシャルによって物体Ａの表面原子に作用する力が求まる。同様に、物体Ｂの第２層め原子が外部原子との関係で有するポテンシャルによって物体Ｂの表面原子に作用する力を求めることができる。
【０１１６】
上述の例においては、物体の表面を、その物体を構成する粒子、特に、物体の最表面から２層めの粒子のポテンシャルを用いて表現し、その２層め粒子は、式（１）で表されるような、角度依存性をもつポテンシャルをつくると考えた。２層め粒子のつくるポテンシャルは式（１）で表されるものに限らない。
【０１１７】
物体Ｂを構成する粒子が近づいてきたときに斥力を受けはじめる距離であって、物体Ａの表面に垂直な方向に沿う距離（図２におけるＲ０）が、物体Ｂを構成する粒子が近づいてきたときに斥力を受けはじめる距離であって、物体Ａの表面に平行な方向に沿う距離よりも短くなるようなポテンシャルを用いることができる。
【０１１８】
また、上述の例においては、２つの異なる物体Ａ及びＢがある場合に、物体Ａを構成する粒子と物体Ｂを構成する粒子との間に作用する引力をゼロとし斥力のみが作用しうるとしたが、ゼロにするのではなく引力を弱くしてもよい。少なくとも、物体Ａの最表面を構成する粒子と、物体Ｂの最表面を構成する粒子との間に、そのような条件を設定する。そうすることによって、接触問題を扱うことが可能になる。
【０１１９】
以下、「物体Ａを構成する粒子と物体Ｂを構成する粒子との間に作用する引力を弱くする」ということについて説明する。両物体Ａ、Ｂ間に作用する引力は、物体Ａを構成する粒子と物体Ｂを構成する粒子との間のポテンシャルを粒子間距離で微分することで求めることができる。したがって、まず、当該ポテンシャルの設定について説明する。
【０１２０】
物体Ａ、物体Ｂが同種の物質で構成されている場合、（ｉ）物体Ａを構成する粒子間のポテンシャル、（ｉｉ）物体Ｂを構成する粒子間のポテンシャル、（ｉｉｉ）物体Ａを構成する粒子と物体Ｂを構成する粒子との間のポテンシャルとして、たとえば文献等に記載のある、物体Ａ及びＢを構成する物質のポテンシャルを採用する。たとえば物体Ａ及びＢが銅（Ｃｕ）で構成されている場合、（ｉ）〜（ｉｉｉ）のポテンシャルとして、すべて銅のポテンシャルを設定する。
【０１２１】
物体Ａ、物体Ｂが異種の物質で構成されている場合、（ｉ）物体Ａを構成する粒子間のポテンシャル、（ｉｉ）物体Ｂを構成する粒子間のポテンシャル、（ｉｉｉ）物体Ａを構成する粒子と物体Ｂを構成する粒子との間のポテンシャルとして、それぞれたとえば文献等に記載のある、物体Ａを構成する物質のポテンシャル、物体Ｂを構成する物質のポテンシャル、物体Ａを構成する物質と物体Ｂを構成する物質との間に働くポテンシャルを採用する。ここで、物体Ａを構成する物質と物体Ｂを構成する物質との間に働くポテンシャルは、物体Ａを構成する物質のポテンシャルと、物体Ｂを構成する物質のポテンシャルに基づいて求められるポテンシャルである。
【０１２２】
ただ、物体Ａを構成する物質と物体Ｂを構成する物質との間に働くポテンシャルは、文献等に記載されていない場合も少なくない。このようなときには、物体Ａを構成する物質のポテンシャルと、物体Ｂを構成する物質のポテンシャルをもとに、たとえばＭｉｘｉｎｇ  Ｒｕｌｅ（混合則）に従って、物体Ａを構成する物質と物体Ｂを構成する物質との間に働くポテンシャルを設定する。
【０１２３】
たとえば、物体Ａを構成する物質が銅（Ｃｕ）であり、物体Ｂを構成する物質がタングステン（Ｗ）である場合、（ｉ）のポテンシャルとしては銅のポテンシャルを採用し、（ｉｉ）のポテンシャルとしてはタングステンのポテンシャルを採用する。更に、（ｉｉｉ）のポテンシャルは、銅のポテンシャルとタングステンのポテンシャルとからＭｉｘｉｎｇ  Ｒｕｌｅによって設定する。銅やタングステンについては、通常、Ｍｏｒｓｅポテンシャルが使用される。このため、Ｍｏｒｓｅポテンシャルを例にとると、Ｍｉｘｉｎｇ  Ｒｕｌｅによって設定される（ｉｉｉ）のポテンシャルパラメータは以下のように示される。
【０１２４】
【数２９】
2011013888.tif 000054
【０１２５】
【数３０】
2011013888.tif 000055
【０１２６】
【数３１】
2011013888.tif 000056
【０１２７】
【数３２】
2011013888.tif 000057
【０１２８】
ここで、φはポテンシャル、Ｄはポテンシャルの深さ、Ａはポテンシャルの幅、ｒ０は安定距離（ポテンシャルの底となる原子間距離）を示す。また、添え字の「Ｃｕ」は銅のポテンシャルパラメータであることを表し、「Ｗ」はタングステンのそれであることを表す。更に、添え字の「Ｃｕ−Ｗ」は、銅とタングステンとの間のポテンシャルパラメータであることを表す。
【０１２９】
なお、物体Ａを構成する物質と物体Ｂを構成する物質とが等しく、たとえば銅である場合、（ｉｉｉ）物体Ａを構成する粒子と物体Ｂを構成する粒子との間のポテンシャルにＭｉｘｉｎｇ  Ｒｕｌｅを適用すると、ポテンシャルパラメータは以下のようになり、（ｉｉｉ）のポテンシャルとして銅のポテンシャルを設定したのと等しくなる。
【０１３０】
【数３３】
2011013888.tif 000058
【０１３１】
【数３４】
2011013888.tif 000059
【０１３２】
【数３５】
2011013888.tif 000060
【０１３３】
このように、物体Ａ及び物体Ｂが同種の物質で構成されている場合も含めて、物体Ａを構成する物質のポテンシャルと、物体Ｂを構成する物質のポテンシャルをもとに、（ｉｉｉ）物体Ａを構成する粒子と物体Ｂを構成する粒子との間に働くポテンシャルを設定する、ということができる。
【０１３４】
（ｉｉｉ）物体Ａを構成する粒子と物体Ｂを構成する粒子との間に働くポテンシャルが設定されたら、それを粒子間距離で微分して、物体Ａを構成する粒子と物体Ｂを構成する粒子との間に作用する引力を求めることができる。
【０１３５】
「物体Ａを構成する粒子と物体Ｂを構成する粒子との間に作用する引力を弱くする」とは、物体Ａを構成する粒子と物体Ｂを構成する粒子とが距離ｄだけ離れている場合に働く引力を、（ｉｉｉ）のポテンシャルの微分を通して求められる、距離ｄだけ離れた物体Ａの構成粒子と物体Ｂの構成粒子との間に働く引力より小さく設定することを意味する。そのまま引力を弱くしていき、引力がなくなった状態が斥力のみの状態である。
【０１３６】
図４は、実施例によるシミュレーション方法を示すフローチャートである。
【０１３７】
実施例によるシミュレーション方法では、まずステップＳ１０１において、第１の物体（たとえば物体Ａ）及び第２の物体（たとえば物体Ｂ）を、それぞれ複数の粒子の集合で定義する。次に、ステップＳ１０２において、第１の物体を構成する粒子間に作用する力を、第１のポテンシャル（たとえば図１（Ｂ）に示されるポテンシャル）を用いて求める。また、同様に、第２の物体を構成する粒子間に作用する力を、第２のポテンシャル（たとえば図１（Ｂ）に示されるポテンシャル、またはそれと類似するポテンシャル。第１のポテンシャルと等しい式で表される場合もあるし、そうでない場合もある。）を用いて求める。第１の物体の構成粒子間、及び第２の物体の構成粒子間に作用する力は、引力の場合もあれば斥力の場合もある。
【０１３８】
続いて、ステップＳ１０３において、第１の物体を構成する粒子と、第２の物体を構成する粒子との間に働く引力を弱くする、またはゼロとすることを前提として、第１の物体の表面から２層めを構成する粒子が、第２の物体を構成する粒子との関係で有する第３のポテンシャル（たとえば図１（Ｃ）に示されるポテンシャル、または式（１）及び式（５）で表されるポテンシャル）によって、第１の物体を構成する粒子と第２の物体を構成する粒子との間に作用する力、及び第１の物体を構成する粒子が受ける作用を求める。
【０１３９】
最後にステップＳ１０４において、Ｓ１０２で求められた力にＳ１０３で求められた力・作用を加えて運動方程式を解く等の計算を行う。すなわち、物体Ａ及びＢを構成する粒子の挙動を計算する。
【０１４０】
実施例によるシミュレーション方法は、プログラムにより、コンピュータで実行することができる。
【０１４１】
図５は、実施例によるシミュレーション方法を用いてシミュレーションを行うシミュレーション装置のシステム構成図である。本図に構成を示したシミュレーション装置を使用して、図４にフローチャートで示したシミュレーション方法を実施することができる。
【０１４２】
まず、キーボードなどの入力装置から、物体Ａ及びＢを構成する粒子の初期状態に関する情報が入力される。また、これと同時に、またはこれに先立って、キーボードなどの入力装置から、第１の物体を構成する粒子間、及び第２の物体を構成する粒子間に作用する力を求めるために必要な、たとえばポテンシャルなどの情報や、第１の物体の第２層めを構成する粒子が、第２の物体を構成する粒子との関係で有するポテンシャルに関する情報を入力し、設定することができる。
【０１４３】
中央処理装置は、メインメモリ中の制御プログラムの指令を受け、メインメモリ中のシミュレーション実行プログラムを実行する。
【０１４４】
シミュレーション実行プログラムは、ステップＳ１０１に対応して、第１の物体及び第２の物体を、それぞれ複数の粒子の集合で定義する部分、ステップＳ１０２に対応して、第１の物体を構成する粒子間に作用する力、及び第２の物体を構成する粒子間に作用する力を求める部分、ステップＳ１０３に対応して、第１の物体を構成する粒子と、第２の物体を構成する粒子との間に働く引力を弱くする、またはゼロとすることを前提として、第１の物体の表面から２層めを構成する粒子が、第２の物体を構成する粒子との関係で有するポテンシャルによって、第１の物体を構成する粒子と第２の物体を構成する粒子との間に作用する力、及び第１の物体を構成する粒子が受ける作用を求める部分、及び、ステップＳ１０４に対応して、Ｓ１０２で求められた力にＳ１０３で求められた力・作用を加え、物体Ａ及びＢを構成する粒子の挙動を計算して、位置、エネルギ等を把握する部分を含む。
【０１４５】
図４のステップＳ１０１を参照して説明したように、第１の物体及び第２の物体を、それぞれ複数の粒子の集合で定義する。
【０１４６】
次に、図４のステップＳ１０２を参照して説明したように、ある時刻における粒子の状態に関する情報に基づいて、第１の物体を構成する粒子間に作用する力、及び、第２の物体を構成する粒子間に作用する力を求める。
【０１４７】
続いて、図４のステップＳ１０３を参照して説明したように、第１の物体を構成する粒子と、第２の物体を構成する粒子との間に働く引力を弱くする、またはゼロとすることを前提として、第１の物体の表面から２層めを構成する粒子が、第２の物体を構成する粒子との関係で有するポテンシャルによって、第１の物体を構成する粒子と第２の物体を構成する粒子との間に作用する力、及び第１の物体を構成する粒子が受ける作用を求める。
【０１４８】
そして図４のステップＳ１０４を参照して説明したように、求められた力・作用を用いて、物体Ａ及びＢを構成する粒子の挙動を計算する。
【０１４９】
計算され、把握された粒子の位置、エネルギ等は、出力装置、たとえばディスプレイに表示される。
【０１５０】
本願発明者らは、実施例によるシミュレーション方法を用いて、球の衝突シミュレーションを行い、ステップＳ１０３において、第１の物体の表面から２層めを構成する粒子が、第２の物体を構成する粒子との関係で有する第３のポテンシャルを、図１（Ａ）に示される球状ポテンシャルとした場合と、図２に示される（式（１）及び（５）で表される）角度依存型ポテンシャルとした場合とについて、運動エネルギの変化を比較した。
【０１５１】
図６（Ａ）及び（Ｂ）を参照して、球の衝突シミュレーションについて説明する。図６（Ａ）に示すように、直径３０ｍｍと直径１５ｍｍの球（ヤング率２０８ＧＰａ）を、２球の中心を結ぶ直線に沿って２００ｍ／ｓの速度で衝突させるシミュレーションを行った。
【０１５２】
図６（Ｂ）は、運動エネルギの変化を示すグラフである。グラフの横軸は、時刻を単位「１×１０−６秒」で示し、縦軸は、２球の運動エネルギの合計を「任意単位」で示す。運動エネルギが減少しはじめている時点が、接触の開始時刻である。
【０１５３】
球状ポテンシャルを用いたシミュレーションにおける接触開始時刻と、接触のアルゴリズムを用いた比較シミュレーションにおけるそれとの間にはずれが生じている。これは球状ポテンシャルの斥力の働きはじめる位置が、物体表面を正確に表現できていないためである。他方、角度依存型ポテンシャルを用いたシミュレーションにおいては、接触開始時刻が、接触のアルゴリズムを用いた比較シミュレーションの場合とほぼ等しい。実施例によるシミュレーション方法において、角度依存型ポテンシャルを使用した場合、斥力の働きはじめる位置で物体表面をうまく表現できるため、一層正確性の高いシミュレーションを行うことができる。
【０１５４】
以上、実施例に沿って本発明を説明したが、本発明はこれらに限定されるものではない。種々の変更、改良、組み合わせ等が可能なことは当業者には自明であろう。
【産業上の利用可能性】
【０１５５】
減速機等の接触が起きる機構、弾性の連成解析に応用可能である。
【符号の説明】
【０１５６】
Ｓ１０１〜Ｓ１０４  ステップ
【図１】
2011013888.tif 000061
【図２】
2011013888.tif 000062
【図３】
2011013888.tif 000063
【図４】
2011013888.tif 000064
【図５】
2011013888.tif 000065
【図６】
2011013888.tif 000066