(57)【要約】【課題】  高効率でシミュレーションを行う。【解決手段】  （ａ）１つの物体を密な領域と疎な領域とに分ける。（ｂ）密な領域を、第１の粒子間距離以下の粒子間距離で分布する複数の粒子で表し、疎な領域を、第１の粒子間距離よりも大きい粒子間距離で分布する複数の粒子で表す。（ｃ）密な領域、疎な領域の各々の中で、粒子の力学的相互作用を定義する。（ｄ）疎な領域との境界における密な領域の表面に存在する粒子と、密な領域との境界における疎な領域の表面に存在する粒子との間の力学的相互作用を定義する。（ｅ）工程（ｃ）及び（ｄ）において定義された力学的相互作用に基づいて、密な領域及び疎な領域内の粒子の挙動をシミュレーションする。【選択図】  図１４
2010039623.tif 000002
【特許請求の範囲】
【請求項１】
（ａ）１つの物体を密な領域と疎な領域とに分ける工程と、
（ｂ）前記密な領域を、第１の粒子間距離以下の粒子間距離で分布する複数の粒子で表し、前記疎な領域を、前記第１の粒子間距離よりも大きい粒子間距離で分布する複数の粒子で表す工程と、
（ｃ）前記密な領域、前記疎な領域の各々の中で、粒子の力学的相互作用を定義する工程と、
（ｄ）前記疎な領域との境界における前記密な領域の表面に存在する粒子と、前記密な領域との境界における前記疎な領域の表面に存在する粒子との間の力学的相互作用を定義する工程と、
（ｅ）前記工程（ｃ）及び（ｄ）において定義された力学的相互作用に基づいて、前記密な領域及び前記疎な領域内の粒子の挙動をシミュレーションする工程と
を有するシミュレーション方法。
【請求項２】
前記工程（ｄ）において、前記疎な領域の表面に接触した前記密な領域の粒子と、該粒子が接触した前記疎な領域の表面位置を囲むＤｅｌａｕｎａｙ三角形分割によって得られた三角形の頂点に位置する、前記疎な領域の表面の３つの粒子との間に、力学的相互作用を定義する請求項１に記載のシミュレーション方法。
【請求項３】
前記工程（ｂ）において、前記疎な領域との境界における前記密な領域の粒子の位置と、前記密な領域との境界における前記疎な領域の粒子の位置とが一致するように粒子を分布させ、
前記工程（ｄ）において、前記密な領域の粒子と前記疎な領域の粒子の位置が一致しているとき、前記密な領域の粒子を消滅させて、該位置に前記疎な領域の粒子のみを配置する請求項１または２に記載のシミュレーション方法。
【請求項４】
種々のシミュレーションを行うためにコンピュータを、
１つの物体を密な領域と疎な領域とに分けるための情報を入力する手段、
前記密な領域を、第１の粒子間距離以下の粒子間距離で分布する複数の粒子で表し、前記疎な領域を、前記第１の粒子間距離よりも大きい粒子間距離で分布する複数の粒子で表すための演算を行う手段、
前記密な領域、前記疎な領域の各々の中で、粒子の力学的相互作用を定義する第１の定義手段、
前記疎な領域との境界における前記密な領域の表面に存在する粒子と、前記密な領域との境界における前記疎な領域の表面に存在する粒子との間の力学的相互作用を定義する第２の定義手段、
前記第１及び第２の定義手段で定義された力学的相互作用に基づいて、前記密な領域及び前記疎な領域内の粒子の挙動をシミュレーションする手段
として機能させるためのプログラム。
【請求項５】
前記第２の定義手段は、前記疎な領域の表面に接触した前記密な領域の粒子と、該粒子が接触した前記疎な領域の表面位置を囲むＤｅｌａｕｎａｙ三角形分割によって得られた三角形の頂点に位置する、前記疎な領域の表面の３つの粒子との間に、力学的相互作用を定義する請求項４に記載のプログラム。
【請求項６】
前記第２の定義手段は、前記密な領域の粒子と前記疎な領域の粒子の位置が一致しているとき、前記密な領域の粒子を消滅させて、該位置に前記疎な領域の粒子のみを配置する請求項４または５に記載のプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、シミュレーション方法及びプログラムに関する。
【背景技術】
【０００２】
コンピュータを用いて物質科学全般の現象を探求する方法として、分子動力学法に基づく分子シミュレーションが知られている。分子シミュレーションによって、分子のポテンシャルエネルギや最安定構造等物質の特性を、分子レベルで解明することが可能である。
【０００３】
分子動力学法に基づく分子シミュレーションに関して、種々の研究がなされている。たとえば少ない計算量で分子動力学計算を実行するシミュレーション方法の発明が開示されている（たとえば、特許文献１参照）。
【０００４】
予め設定された複雑な形状の物性値または物理量を求める分子シミュレーションを行う場合、適当に原子を配置したとしても、配置された原子位置に基づいてポテンシャルを修正することにより、原子間相互作用のパラメータを定めることが可能でる。
【０００５】
しかしながら、原子間距離があまりに異なる領域が存在すると、強度（ヤング率）に異方性が生じてしまう。具体的には、Ｖｏｒｏｎｏｉ解析を行ったときの要素（メッシュ）の面積（２次元の場合）、もしくは体積（３次元の場合）の分散が大きいときに不具合が生じやすい。
【０００６】
【特許文献１】特開２００６−２８５８６６号公報
【発明の開示】
【発明が解決しようとする課題】
【０００７】
本発明の目的は、高効率でシミュレーションを行うシミュレーション方法を提供することである。
【０００８】
また、高効率でシミュレーションを行うシミュレーション方法を実行させるプログラムを提供することである。
【課題を解決するための手段】
【０００９】
本発明の一観点によれば、（ａ）１つの物体を密な領域と疎な領域とに分ける工程と、（ｂ）前記密な領域を、第１の粒子間距離以下の粒子間距離で分布する複数の粒子で表し、前記疎な領域を、前記第１の粒子間距離よりも大きい粒子間距離で分布する複数の粒子で表す工程と、（ｃ）前記密な領域、前記疎な領域の各々の中で、粒子の力学的相互作用を定義する工程と、（ｄ）前記疎な領域との境界における前記密な領域の表面に存在する粒子と、前記密な領域との境界における前記疎な領域の表面に存在する粒子との間の力学的相互作用を定義する工程と、（ｅ）前記工程（ｃ）及び（ｄ）において定義された力学的相互作用に基づいて、前記密な領域及び前記疎な領域内の粒子の挙動をシミュレーションする工程とを有するシミュレーション方法が提供される。
【００１０】
また、本発明の他の観点によると、種々のシミュレーションを行うためにコンピュータを、１つの物体を密な領域と疎な領域とに分けるための情報を入力する手段、前記密な領域を、第１の粒子間距離以下の粒子間距離で分布する複数の粒子で表し、前記疎な領域を、前記第１の粒子間距離よりも大きい粒子間距離で分布する複数の粒子で表すための演算を行う手段、前記密な領域、前記疎な領域の各々の中で、粒子の力学的相互作用を定義する第１の定義手段、前記疎な領域との境界における前記密な領域の表面に存在する粒子と、前記密な領域との境界における前記疎な領域の表面に存在する粒子との間の力学的相互作用を定義する第２の定義手段、前記第１及び第２の定義手段で定義された力学的相互作用に基づいて、前記密な領域及び前記疎な領域内の粒子の挙動をシミュレーションする手段として機能させるためのプログラムが提供される。
【発明の効果】
【００１１】
本発明によれば、高効率でシミュレーションを行うシミュレーション方法を提供することができる。
【００１２】
また、高効率でシミュレーションを行うシミュレーション方法を実行させるプログラムを提供することができる。
【発明を実施するための最良の形態】
【００１３】
２次元の場合と３次元の場合とに分けて、原子（粒子）が細かく配置された領域と、粗く配置された領域とを、両領域間で物理量が連続的に変化するように接着（Ｇｌｕｅ）する手法の実施例について説明する。
【００１４】
まず、図１〜図７を参照し、２次元の場合について、実施例による粗密領域接着手法を説明する。
【００１５】
図１（Ａ）に原子が密に配置された２次元領域と、原子が疎に配置された２次元領域を示す。
【００１６】
図１（Ｂ）に示すように、これら２つの領域を接着して、１つの領域を形成する。この際、密な領域の内部と疎な領域の内部は、それぞれポテンシャル修正法で相互作用のパラメータを決定しておく。
【００１７】
図２（Ａ）に、密に配置された領域の原子と、疎に配置された領域の原子の接着前の状態を示す。
【００１８】
図の上方の原子ｇ０〜ｇ４は密に配置された領域に属する原子である。原子ｇ０とｇ１とはバネ定数ｋ０１のバネで結合しているとする。同様に、原子ｇ１とｇ２とはバネ定数ｋ１２のバネで、原子ｇ２とｇ３とはバネ定数ｋ２３のバネで、原子ｇ３とｇ４とはバネ定数ｋ３４のバネで結合していると考える。
【００１９】
図の下方の原子ｉ及びｊは疎に配置された領域に属する原子である。原子ｉとｊとはバネ定数ｋｉｊのバネで結合しているとする。
【００２０】
疎な領域の原子ｉ、ｊ上に、密な領域の原子ｇ０〜ｇ４を接着する。原子ｇ０〜ｇ４は、疎な領域の表面を構成するライン（原子ｉと原子ｊとを結ぶ直線）上に存在する。
【００２１】
図２（Ｂ）を参照する。接着に当たって、原子ｇ０、ｇ４の位置が、それぞれ原子ｉ、ｊの位置と一致するときには、密な領域の原子ｇ０、ｇ４を消去し、代わりに疎な領域の原子ｉ、ｊを用いる。原子ｇ０、ｇ４の位置が、それぞれ原子ｉ、ｊの位置と極めて近い場合も同じく処理する。
【００２２】
また、密な領域の原子同士を結合していたバネ、及び疎な領域の原子同士を結合していたバネをすべて消去する。代わりに、隣り合う原子同士がすべてバネ定数ｋｉｊのバネで結合すると考える。具体的には、原子ｉと原子ｇ１、原子ｇ１と原子ｇ２、原子ｇ２と原子ｇ３、及び原子ｇ３と原子ｊとはすべてバネ定数ｋｉｊのバネで結合すると考える。
【００２３】
図２（Ｃ）を参照して、接着後、原子ｉと原子ｊとを結ぶ方向に直交する方向に作用する力について説明する。
【００２４】
原子ｉと原子ｊとを結ぶ直線から、当該直線と直交する方向にΔｒだけ変位している原子には、変位している方向と反対方向に変位Δｒに比例する力が加わると考える。たとえば、図２（Ｃ）において、原子ｉと原子ｊとを結ぶ直線と直交する方向にΔｒだけ変位している原子ｇ１には、
2010039623.tif 000003
【００２５】
・・（１）
の力が作用するとする。ここで負の符号は力Ｆの作用する向きが変位Δｒと反対方向であることを示す。力Ｆは原子ｉと原子ｊとを結ぶ直線に直交する方向に働く。
【００２６】
また、式（１）の力Ｆの反作用として、原子ｉに作用する力Ｆｉ、原子ｊに作用する力Ｆｊをそれぞれ以下の式（２）、（３）で導入する。
2010039623.tif 000004
【００２７】
・・（２）
2010039623.tif 000005
【００２８】
・・（３）
ここでｌｉ、ｌｊはそれぞれ原子ｇ１の存する点から直線ｉｊに下ろした垂線の足と、原子ｉ、原子ｊとの間の距離を示す。力Ｆｉ、Ｆｊの働く向きは、力Ｆの働く向きと反対方向である。
【００２９】
力Ｆ、Ｆｉ、Ｆｊを式（１）〜（３）で示すように定めることで、図示の系内で力とトルクのバランスが保たれる。
【００３０】
このようにして、原子が密に配置された２次元領域と、疎に配置された２次元領域とを接着する。
【００３１】
本願発明者らは、上述のように原子密度の異なる２つの２次元領域を接着して、精度確認のためのシミュレーションを行った。
【００３２】
図３〜図７を参照して、シミュレーションについて説明する。
【００３３】
図３（Ａ）及び（Ｂ）に、精度確認のための計算系を示す。
【００３４】
図３（Ａ）を参照する。半径６ｍｍの半円筒と、半径２０ｍｍの半円筒の円筒面を、お互いの頂点部分が接触するように接触させ、半径６ｍｍの半円筒側から接触点に垂直に、５００Ｎ／ｍｍの荷重をかけるシミュレーションを行った。
【００３５】
シミュレーションに当たっては、２つの半円筒を形成する材料にＳＵＳを想定し、ヤング率を２０８ＧＰａ、密度を７８００ｋｇ／ｍ３、ポアソン比を０．２７に設定した。
【００３６】
図３（Ｂ）は、２つの半円筒の接触部分、及び粗密両領域が接着された部分の拡大写真である。シミュレーションにおいては、２つの半円筒の接触部分の近傍にのみ、密領域を導入し、２つの半円筒の接近距離、及び接近距離と理論値との間の誤差を求めた。なお、上記条件で接触された２つの円筒が接近する距離の理論値は１６．８４μｍである。
【００３７】
図４（Ａ）及び（Ｂ）を参照する。まず、比較例として２つの半円筒に粗密領域を設けず、メッシュ幅を一様としてシミュレーションを行った。メッシュ幅は半円筒を構成する原子の原子間距離に相当する。
【００３８】
図４（Ａ）を参照して、第１の比較例について説明する。第１の比較例においては、半径６ｍｍの半円筒のメッシュ幅を０．３ｍｍとし、半径２０ｍｍの半円筒のメッシュ幅を０．５ｍｍとしてシミュレーションを行った。
【００３９】
第１の比較例において、２つの半円筒の接近距離は１７．７８μｍであり、理論値との誤差は、５．５８％であった。
【００４０】
図４（Ｂ）を参照して、第２の比較例について説明する。第２の比較例においては、半径６ｍｍの半円筒のメッシュ幅を０．０８ｍｍとし、半径２０ｍｍの半円筒のメッシュ幅を０．１ｍｍとしてシミュレーションを行った。
【００４１】
第２の比較例において、２つの半円筒の接近距離は１７．２０μｍであり、理論値との誤差は、２．１４％であった。
【００４２】
第１の比較例と第２の比較例とを比べたとき、メッシュ幅の小さい第２の比較例の方が理論値との誤差が小さいことがわかる。
【００４３】
図５を参照して、基本シミュレーションについて説明する。基本シミュレーションにおいては、疎領域については、第１の比較例と等しいメッシュ幅を用い、密領域については、第２の比較例と等しいメッシュ幅を用いた。すなわち、半径６ｍｍの半円筒については、疎領域のメッシュ幅を０．３ｍｍ、密領域のメッシュ幅を０．０８ｍｍとし、半径２０ｍｍの半円筒については、疎領域のメッシュ幅を０．５ｍｍ、密領域のメッシュ幅を０．１ｍｍとして、２つの半円筒の接近距離を計算した。その結果、２つの半円筒の接近距離は１７．４８μｍ、理論値との誤差は３．２１％であった。
【００４４】
全体のメッシュ幅を一様に細かくした第２の比較例には及ばないが、メッシュ幅が一様に粗い第１の比較例に比べると、理論値との誤差が小さくなっていることがわかる。
【００４５】
次に本願発明者らは、上述の基本シミュレーションよりも、疎領域のメッシュ幅を粗くしてシミュレーションを行った。
【００４６】
図６（Ａ）を参照する。本図に示す例では、半径６ｍｍの半円筒の疎領域のメッシュ幅を０．６ｍｍ、半径２０ｍｍの半円筒の疎領域のメッシュ幅を１．０ｍｍとした。その結果、２つの半円筒の接近距離は１７．５０μｍ、理論値との誤差は３．９２％であった。
【００４７】
図６（Ｂ）を参照する。本図に示すのは、図６（Ａ）に示す例より更に疎領域のメッシュ幅を粗くした例である。本例においては、半径６ｍｍの半円筒の疎領域のメッシュ幅を１．２ｍｍ、半径２０ｍｍの半円筒の疎領域のメッシュ幅を２．０ｍｍとした。２つの半円筒の接近距離は１７．３３μｍ、理論値との誤差は２．９１％という結果が得られた。
【００４８】
図６（Ａ）及び（Ｂ）を参照して説明したシミュレーションと、基本シミュレーションとの比較から、誤差は、疎領域のメッシュ幅の大小にかかわらず、ほぼ等しくなることがわかる。
【００４９】
最後に、本願発明者らは、理論値との誤差とアスペクト比との関係について調査した。アスペクト比は、疎領域のメッシュ幅を密領域のメッシュ幅で除したものとして定義される。
【００５０】
調査に当たっては、図３（Ａ）及び（Ｂ）を参照して説明した計算系において、２つの半円筒の疎領域のメッシュ幅をともに０．５ｍｍに固定し、密領域のメッシュ幅を様々に変化させることでアスペクト比を変化させた。そして変化させたアスペクト比を与えるメッシュ幅ごとに、２つの半円筒の接近距離、及び接近距離と理論値との間の誤差を求めた。
【００５１】
図７は、アスペクト比と誤差との関係を示すグラフである。グラフの横軸はアスペクト比を任意単位で表し、縦軸は、誤差を単位「％」で表す。
【００５２】
アスペクト比が大きくなるほど誤差は小さくなっている。これはアスペクト比を大きくするために、密領域のメッシュ幅を小さくしたためであると考えられる。
【００５３】
グラフから、実施例による接着手法は、アスペクト比を大きくすることによる不都合は生じないものと考えられる。
【００５４】
次に、図８〜図１３を参照し、３次元の場合について、実施例による粗密領域接着手法を説明する。
【００５５】
図８（Ａ）に原子が密に配置された３次元領域と、原子が疎に配置された３次元領域を示す。
【００５６】
図８（Ｂ）に示すように、これら２つの領域を接着して、１つの領域を形成する。この際、密な領域の内部と疎な領域の内部は、それぞれポテンシャル修正法で相互作用のパラメータを決定しておく。
【００５７】
図９（Ａ）に、密に配置された領域の原子と、疎に配置された領域の原子の接着前の状態を示す。
【００５８】
図の左側の原子ｇ０〜ｇ５は密に配置された領域に属する原子である。原子ｇ０、ｇ２、ｇ５が形成する三角形の原子ｇ０、ｇ２を結ぶ辺の中点に原子ｇ１が配置されている。また、原子ｇ２、ｇ５を結ぶ辺の中点に原子ｇ４が、原子ｇ５、ｇ０を結ぶ辺の中点に原子ｇ３が配置されている。
【００５９】
図の右側の原子ｉ、ｊ、及びｋは疎に配置された領域に属する原子である。原子ｉ、ｊ、ｋが形成する三角形は、密領域の原子ｇ０、ｇ２、ｇ５が形成する三角形と合同である。
【００６０】
原子ｉとｊとはバネ定数ｋｉｊのバネで結合しているとする。同様に、原子ｊとｋとはバネ定数ｋｊｋのバネで、原子ｋとｉとはバネ定数ｋｋｉのバネで結合しているとする。
【００６１】
疎な領域の原子ｉ、ｊ、ｋのつくる三角形上またはその内部に、密な領域の原子ｇ０〜ｇ５を接着する。
【００６２】
図９（Ｂ）を参照する。接着に当たって、密な領域の原子位置と疎な領域の原子位置とが一致した場合、密な領域の原子を取り除き、代わりに疎な領域の原子を用いる。
【００６３】
図９（Ｂ）においては、それぞれ疎な領域の原子ｉ、ｊ、ｋと位置が一致する密な領域の原子ｇ０、ｇ２、ｇ５を取り除き、その位置には原子ｉ、ｊ、ｋを残す。原子ｇ０、ｇ２、ｇ５の位置が、それぞれ原子ｉ、ｊ、ｋの位置と極めて近い場合も同じく処理する。
【００６４】
図９（Ｃ）を参照する。疎な領域の原子同士を結合していたバネをすべて消去し、代わりに、隣り合う原子同士を結ぶバネを想定する。原子ｉ、ｊを結ぶ辺に平行な方向に隣り合う原子間を結ぶバネのバネ定数はすべてｋｉｊとする。同様に、原子ｊ、ｋを結ぶ辺に平行な方向に隣り合う原子間を結ぶバネのバネ定数はすべてｋｊｋとし、原子ｋ、ｉを結ぶ辺に平行な方向に隣り合う原子間を結ぶバネのバネ定数はすべてｋｋｉとする。
【００６５】
三次元形状にメッシュを切る場合、図９（Ａ）〜（Ｃ）に示すように、接着箇所の粗密両領域の原子配置を相似形とすることが望ましい。しかしながら実際には、そうすることは困難であり、全く異なった形状に配置された粗密領域の原子同士を接着する場合が多い。そこでそのような場合には、以下のルールで結合させる。以下のルールを用いることにより、疎密両領域のメッシュを相似にしなくても接着が可能となる。
【００６６】
図１０を参照する。
【００６７】
原子ｉ、ｊ、ｋを疎な領域に属する原子とする。原子ｉｊ間の距離をｒｉｊ、原子ｊｋ間の距離をｒｊｋ、原子ｋｉ間の距離をｒｋｉとする。また、原子ｉとｊとはバネ定数ｋｉｊのバネで結合しているとする。同様に、原子ｊとｋとはバネ定数ｋｊｋのバネで、原子ｋとｉとはバネ定数ｋｋｉのバネで結合しているとする。
【００６８】
疎な領域の原子ｉ、ｊ、ｋで形成される三角形上またはその内部（三角形要素上）に、密な領域の原子ｇを接着する。原子ｉ、ｊ、ｋで形成される三角形は、原子ｇを囲む、Ｄｅｌａｕｎａｙ三角形分割によって得られた三角形である。
【００６９】
なお、Ｄｅｌａｕｎａｙ三角形分割とは、隣接するＶｏｒｏｎｏｉ点同士をつないで領域を分割することをいう。Ｖｏｒｏｎｏｉ点とは、Ｖｏｒｏｎｏｉ多面体（２次元ならば垂直二等分線、３次元ならば垂直二等分面で囲まれた領域）を構成する点である。
【００７０】
このとき原子ｇは、原子ｉ、ｊ、ｋとそれぞれ下の式（４）〜（６）で表されるバネ定数ｋｇｉ、ｋｇｊ、ｋｇｋのバネで結合されると考える。
2010039623.tif 000006
【００７１】
・・（４）
2010039623.tif 000007
【００７２】
・・（５）
2010039623.tif 000008
【００７３】
・・（６）
ここで、ｒｇｉ、ｒｇｊ、ｒｇｋはそれぞれ、原子ｇと原子ｉ、原子ｇと原子ｊ、原子ｇと原子ｋ間の距離を表す。またＮは、原子ｉ、ｊ、ｋで形成される三角形上またはその内部に存在する密領域の原子（疎領域に接着される原子）の数である。更に、ｒｒ、ｋｋはそれぞれ下の式（７）、（８）で求められる値である。
2010039623.tif 000009
【００７４】
・・（７）
2010039623.tif 000010
【００７５】
・・（８）
すなわちｒｒは、原子ｉ、ｊ、ｋの原子間距離ｒｉｊ、ｒｊｋ、ｒｋｉの平均値であり、ｋｋは原子ｉ、ｊ、ｋ間を結ぶバネのバネ定数ｋｉｊ、ｋｊｋ、ｋｋｉの平均値である。
【００７６】
次に、図１１を参照して、接着後、原子ｉ、ｊ、ｋで定められる平面に垂直な方向に作用する力について説明する。疎領域原子ｉ、ｊ、ｋで形成される三角形上またはその内部に接着された密領域の原子ｇが、原子ｉ、ｊ、ｋで定められる平面からΔｒだけ離れた場合、原子ｇ、ｉ、ｊ、ｋに働く、平面に垂直な方向の力を、それぞれＦｇ、Ｆｉ、Ｆｊ、Ｆｋとすると、それらはそれぞれ下式（９）〜（１２）で表されると考える。
2010039623.tif 000011
【００７７】
・・（９）
2010039623.tif 000012
【００７８】
・・（１０）
2010039623.tif 000013
【００７９】
・・（１１）
2010039623.tif 000014
【００８０】
・・（１２）
ここで、Ｋは比例定数であり、（９）式において負符号は変位Δｒの方向とは反対向きに力Ｆｇが作用することを示す。また、ｌ１、ｌ２、ｌ３はそれぞれ、原子ｇから、原子ｉ、ｊ、ｋで定められる平面に下ろした垂線の足Ｈから、原子ｊ、ｋ、ｉまでの距離である。更に、θ１は、原子ｉと垂線の足Ｈを結ぶ直線と、原子ｊと垂線の足Ｈを結ぶ直線とのなす角度であり、θ２は、原子ｉと垂線の足Ｈを結ぶ直線と、原子ｋと垂線の足Ｈを結ぶ直線とのなす角度である。また、δは下式（１３）で計算される値である。
2010039623.tif 000015
【００８１】
・・（１３）
力Ｆｇ、Ｆｉ、Ｆｊ、Ｆｋを式（９）〜（１２）で示すように定めることで、図示の系内で力とトルクのバランスが保たれる。
【００８２】
このようにして、原子が密に配置された３次元領域と、疎に配置された３次元領域とを接着する。
【００８３】
本願発明者らは、２次元の粗密領域の接着の場合と同様に、原子密度の異なる２つの３次元領域を接着して、精度確認のためのシミュレーションを行った。
【００８４】
図１２及び図１３を参照して、シミュレーションについて説明する。
【００８５】
図１２（Ａ）及び（Ｂ）に、精度確認のための計算系を示す。
【００８６】
図１２（Ａ）を参照する。直径３０ｍｍの球と直径１３ｍｍの球とを、２球の中心を結ぶ直線に沿って、相対速度１００ｍ／ｓで衝突させるシミュレーションを行い、２球の接触時間τからヤング率を導出し、設定値との比較を行った。シミュレーションにおいて、ヤング率は２０８ＧＰａに設定した。なお、接触時間τとヤング率Ｅとの関係は、下式（１４）〜（１６）を用いて表される。
2010039623.tif 000016
【００８７】
・・（１４）
2010039623.tif 000017
【００８８】
・・（１５）
2010039623.tif 000018
【００８９】
・・（１６）
式（１４）においてμは換算質量、νは２球の相対速度を示す。また、式（１５）において、Ｒ及びＲ´は、２球それぞれの半径を示す。更に、式（１６）におけるσはポアソン比である。
【００９０】
シミュレーションに当たっては、２つの球を形成する材料にＳＵＳを想定し、密度を７８００ｋｇ／ｍ３、ポアソン比を０．２７に設定した。また、先に述べたように、ヤング率を２０８ＧＰａとした。このとき式（１４）〜（１６）を用いて、接触時間τの理論値は、約２１．６１μｍと計算される。
【００９１】
図１２（Ｂ）に、２つの球の中心を通る断面を示す。シミュレーションにおいては、２つの球の表面近傍に密領域を導入し、それ以外の球の内部領域（中心部）に疎領域を画定して、上述のルールで接着した。
【００９２】
シミュレーションは、疎密両領域のアスペクト比が２、３、４、５、１０の各場合について行った。なお、比較例としてアスペクト比１の場合（領域の密度が一様とされ、疎密領域の接着が行われない場合）のシミュレーションも実施した。
【００９３】
図１３（Ａ）及び（Ｂ）にシミュレーションの結果を示す。図１３（Ａ）は、アスペクト比と接触時間τとの関係を表し、図１３（Ｂ）は、アスペクト比とヤング率Ｅとの関係を表す。両図ともに、横軸はアスペクト比を任意単位で示す。また、図１３（Ａ）における縦軸は、接触時間τを単位「μｓ」で示し、図１３（Ｂ）における縦軸は、ヤング率Ｅを単位「ＧＰａ」で示す。
【００９４】
アスペクト比が大きくなるにつれ、接触時間τは短く（図１３（Ａ））、ヤング率Ｅは大きくなる（図１３（Ｂ））ことがわかる。しかしアスペクト比が１０以下であるならば、たとえばヤング率Ｅは設定値から５％以内の誤差で求められることが確認される。
【００９５】
図１４は、実施例によるシミュレーション方法を示すフローチャートである。
【００９６】
まず、ステップＳ１０１において、１つの物体を密な領域と疎な領域とに分ける。
【００９７】
そして、ステップＳ１０２において、密な領域、疎な領域のそれぞれに粒子を配置する。すなわち、疎な領域を、疎に分布する複数の粒子で表し、密な領域を、密に分布する複数の粒子で表す。
【００９８】
次に、ステップＳ１０３において、疎密それぞれの領域に対してポテンシャル修正法を用い、相互作用のパラメータを決定する。すなわち、疎な領域、密な領域の各々の粒子の力学的相互作用を定義する。
【００９９】
続いて、ステップＳ１０４において、疎な領域の表面に存在する粒子と密な領域の表面に存在する粒子との間の力学的相互作用を定義する。
【０１００】
具体的には、ステップＳ１０４ａにおいて、疎密それぞれの領域が２次元領域である場合、疎な領域の表面を構成するライン上に存在する、密な領域の粒子を探索する。疎密両域が３次元領域である場合は、疎な領域の表面を構成する三角形要素上に存在する、密な領域の粒子を探索する。
【０１０１】
そして、ステップＳ１０４ｂにおいて、２次元の場合はライン上、３次元の場合は三角形要素上にある密な領域の粒子と、疎な領域の粒子とを、図１〜図１３を参照して説明した接着手法を用いて接着（Ｇｌｕｅ）し、相互作用のパラメータを決定する。
【０１０２】
最後に、ステップＳ１０５において、分子動力学計算（粒子法によるシミュレーション）を実行する。ステップＳ１０３及びＳ１０４において定義された力学的相互作用に基づいて、密な領域及び疎な領域内の粒子の挙動をシミュレーションする。
【０１０３】
実施例によるシミュレーション方法を用いると、たとえば接触問題を含む分子シミュレーションを行うに当たって、物体同士が接触する箇所には密に粒子を配置し、それ以外の箇所には疎に粒子を配置し、これらを接着（Ｇｌｕｅ）することで、全体の粒子数を少なくし、シミュレーションの効率を高めることができる。
【０１０４】
従来のポテンシャルを修正する手法だけでは、モデルに配置する粒子間距離（メッシュの大きさ）は、一番小さい領域に律速してしまい、特に３次元解析を行う場合、極端に粒子数が多くなってしまって、実質的には計算を行うことが不可能であった。しかしながら、実施例によるシミュレーション方法によれば、接着（Ｇｌｕｅ）の考え方を用いることにより、２次元解析であれば１０分の１、３次元解析であれば１００分の１程度の粒子数で、接着（Ｇｌｕｅ）をしない場合と同等の結果を得ることが可能である。
【０１０５】
図１４にフローチャートを示したシミュレーション方法は、プログラムにより、コンピュータで実行することができる。
【０１０６】
図１５は、実施例によるシミュレーション方法を実行するシミュレーション実行装置のシステム構成図である。
【０１０７】
まず、キーボードなどの入力装置から、シミュレーション対象（物体）に関する情報、たとえば大きさ、形状、材料などの情報が入力される。そしてそのシミュレーション対象（物体）を、どのように疎な領域と密な領域に分けるかを決定する情報、たとえば画定される疎密両領域が２次元領域か、３次元領域か、それらの領域の形状や大きさ（範囲）、更に、それぞれの領域における粒子の粒子間距離等の、疎領域及び密領域に関する情報が入力される。
【０１０８】
中央処理装置は、メインメモリ中の制御プログラムの指令を受け、メインメモリに保存されている相互作用パラメータ決定プログラムを実行する。
【０１０９】
相互作用パラメータ決定プログラムは、入力された情報を基に、前図Ｓ１０１に示したように、１つの物体を密な領域と疎な領域とに分ける。そしてＳ１０２に示したように、密な領域、疎な領域のそれぞれに粒子を配置する位置を決定する。すなわち、疎な領域を、疎に分布する複数の粒子で表し、密な領域を、密に分布する複数の粒子で表すための演算を行う。
【０１１０】
そしてＳ１０３に示したように、疎密それぞれの領域に対してポテンシャル修正法を用い、相互作用のパラメータを決定する。すなわち、疎な領域、密な領域の各々の中での粒子の力学的相互作用を定義する。
【０１１１】
更に、ステップＳ１０４に示したように、疎な領域の表面に存在する粒子と密な領域の表面に存在する粒子との間の力学的相互作用を定義する。具体的には、ステップＳ１０４ａに示したように、疎密それぞれの領域が２次元領域である場合、疎な領域の表面を構成するライン上に存在する、密な領域の粒子を探索する計算を行う。疎密両域が３次元領域である場合は、疎な領域の表面を構成する三角形要素上に存在する、密な領域の粒子を探索する計算を行う。そして、Ｓ１０４ｂに示したように、２次元の場合はライン上、３次元の場合は三角形要素上にある密な領域の粒子と、疎な領域の粒子とを、図１〜図１３を参照して説明した接着手法を用いて接着（Ｇｌｕｅ）し、相互作用のパラメータを決定する。
【０１１２】
その後、中央処理装置は、メインメモリ中の制御プログラムの指令を受け、メインメモリに保存されているシミュレーション実行プログラムを実行し、Ｓ１０５に示したように、粒子法によるシミュレーションを実行する。
【０１１３】
シミュレーションの結果は、ディスプレイやプリンタなどの出力装置を介して表示される。
【０１１４】
以上実施例に沿って本発明を説明したが、本発明はこれらに制限されるものではない。種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。
【産業上の利用可能性】
【０１１５】
原理的には、量子的な現象を除くすべての古典現象を解析するシミュレーションに利用することができる。特に、機構、弾性の動解析における摩擦による発熱、破壊の計算に好適に用いることができる。
【図面の簡単な説明】
【０１１６】
【図１】（Ａ）及び（Ｂ）は、２次元疎密領域の接着について示す図である。
【図２】（Ａ）は、密に配置された領域の原子と、疎に配置された領域の原子の接着前の状態を示し、（Ｂ）は、両領域の原子の接着後の状態を示し、（Ｃ）は、接着後、原子ｉと原子ｊとを結ぶ方向に直交する方向に作用する力を示す図である。
【図３】（Ａ）及び（Ｂ）は、精度確認のための計算系を示す図である。
【図４】（Ａ）及び（Ｂ）は、それぞれ第１、第２の比較例について説明するための図である。
【図５】基本シミュレーションについて説明するための図である。
【図６】（Ａ）及び（Ｂ）は、基本シミュレーションよりも、疎領域のメッシュ幅を粗くして行ったシミュレーションについて説明するための図である。
【図７】アスペクト比と誤差との関係を示すグラフである。
【図８】（Ａ）及び（Ｂ）は、３次元疎密領域の接着について示す図である。
【図９】（Ａ）は、密に配置された領域の原子と、疎に配置された領域の原子の接着前の状態を示し、（Ｂ）及び（Ｃ）は、両領域の原子の接着後の状態を示す。
【図１０】密に配置された領域の原子と、疎に配置された領域の原子の接着後の状態を示す図である。
【図１１】疎密両領域の接着後、疎領域原子ｉ、ｊ、ｋで定められる平面に垂直な方向に作用する力について説明するための図である。
【図１２】（Ａ）及び（Ｂ）は、精度確認のための計算系を示す図である。
【図１３】（Ａ）は、アスペクト比と接触時間τとの関係を表す図であり、（Ｂ）は、アスペクト比とヤング率Ｅとの関係を表す図である。
【図１４】実施例によるシミュレーション方法を示すフローチャートである。
【図１５】実施例によるシミュレーション方法を実行するシミュレーション実行装置のシステム構成図である。
【符号の説明】
【０１１７】
Ｓ１０１〜Ｓ１０５、Ｓ１０４ａ、Ｓ１０４ｂ  ステップ
【図２】
2010039623.tif 000019
【図７】
2010039623.tif 000020
【図９】
2010039623.tif 000021
【図１０】
2010039623.tif 000022
【図１１】
2010039623.tif 000023
【図１３】
2010039623.tif 000024
【図１４】
2010039623.tif 000025
【図１５】
2010039623.tif 000026
【図１】
2010039623.tif 000027
【図３】
2010039623.tif 000028
【図４】
2010039623.tif 000029
【図５】
2010039623.tif 000030
【図６】
2010039623.tif 000031
【図８】
2010039623.tif 000032
【図１２】
2010039623.tif 000033