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

Description

本発明は、シミュレーション方法及びプログラムに関する。

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

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

予め設定された複雑な形状の物性値または物理量を求める分子シミュレーションを行う場合、適当に原子を配置したとしても、配置された原子位置に基づいてポテンシャルを修正することにより、原子間相互作用のパラメータを定めることが可能でる。

しかしながら、原子間距離があまりに異なる領域が存在すると、強度（ヤング率）に異方性が生じてしまう。具体的には、Ｖｏｒｏｎｏｉ解析を行ったときの要素（メッシュ）の面積（２次元の場合）、もしくは体積（３次元の場合）の分散が大きいときに不具合が生じやすい。

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

本発明の目的は、高効率でシミュレーションを行うシミュレーション方法を提供することである。

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

本発明の一観点によれば、
（ａ）シミュレーション対象の１つの物体を密な領域と疎な領域とに分ける工程と、
（ｂ）前記密な領域を、第１の粒子間距離以下の粒子間距離で分布する複数の粒子で表し、前記疎な領域を、前記第１の粒子間距離よりも大きい粒子間距離で分布する複数の粒子で表す工程と、
（ｃ）前記密な領域、前記疎な領域の各々の中で、粒子位置に基づいてポテンシャルを修正することにより粒子の力学的相互作用を定義する工程と、
（ｄ）前記疎な領域との境界における前記密な領域の表面に存在する粒子と、前記密な領域との境界における前記疎な領域の表面に存在する粒子との間の力学的相互作用を、前記密な領域と前記疎な領域との境界に存在する粒子からなる系内で力とトルクのバランスが保たれるように定義する工程と、
（ｅ）前記工程（ｃ）及び（ｄ）において定義された力学的相互作用に基づいて、前記密な領域及び前記疎な領域内の粒子の挙動をシミュレーションする工程と
を有するシミュレーション方法が提供される。

また、本発明の他の観点によると、
種々のシミュレーションを行うためにコンピュータを、
シミュレーション対象の１つの物体を密な領域と疎な領域とに分けるための情報を入力する手段、
前記密な領域を、第１の粒子間距離以下の粒子間距離で分布する複数の粒子で表し、前記疎な領域を、前記第１の粒子間距離よりも大きい粒子間距離で分布する複数の粒子で表すための演算を行う手段、
前記密な領域、前記疎な領域の各々の中で、粒子位置に基づいてポテンシャルを修正することにより粒子の力学的相互作用を定義する第１の定義手段、
前記疎な領域との境界における前記密な領域の表面に存在する粒子と、前記密な領域との境界における前記疎な領域の表面に存在する粒子との間の力学的相互作用を、前記密な領域と前記疎な領域との境界に存在する粒子からなる系内で力とトルクのバランスが保たれるように定義する第２の定義手段、
前記第１及び第２の定義手段で定義された力学的相互作用に基づいて、前記密な領域及び前記疎な領域内の粒子の挙動をシミュレーションする手段
として機能させるためのプログラムが提供される。

本発明によれば、高効率でシミュレーションを行うシミュレーション方法を提供することができる。

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

２次元の場合と３次元の場合とに分けて、原子（粒子）が細かく配置された領域と、粗く配置された領域とを、両領域間で物理量が連続的に変化するように接着（Ｇｌｕｅ）する手法の実施例について説明する。

まず、図１〜図７を参照し、２次元の場合について、実施例による粗密領域接着手法を説明する。

図１（Ａ）に原子が密に配置された２次元領域と、原子が疎に配置された２次元領域を示す。

図１（Ｂ）に示すように、これら２つの領域を接着して、１つの領域を形成する。この際、密な領域の内部と疎な領域の内部は、それぞれポテンシャル修正法で相互作用のパラメータを決定しておく。

図２（Ａ）に、密に配置された領域の原子と、疎に配置された領域の原子の接着前の状態を示す。

図の上方の原子ｇ_０〜ｇ_４は密に配置された領域に属する原子である。原子ｇ_０とｇ_１とはバネ定数ｋ_０１のバネで結合しているとする。同様に、原子ｇ_１とｇ_２とはバネ定数ｋ_１２のバネで、原子ｇ_２とｇ_３とはバネ定数ｋ_２３のバネで、原子ｇ_３とｇ_４とはバネ定数ｋ_３４のバネで結合していると考える。

図の下方の原子ｉ及びｊは疎に配置された領域に属する原子である。原子ｉとｊとはバネ定数ｋ_ｉｊのバネで結合しているとする。

疎な領域の原子ｉ、ｊ上に、密な領域の原子ｇ_０〜ｇ_４を接着する。原子ｇ_０〜ｇ_４は、疎な領域の表面を構成するライン（原子ｉと原子ｊとを結ぶ直線）上に存在する。

図２（Ｂ）を参照する。接着に当たって、原子ｇ_０、ｇ_４の位置が、それぞれ原子ｉ、ｊの位置と一致するときには、密な領域の原子ｇ_０、ｇ_４を消去し、代わりに疎な領域の原子ｉ、ｊを用いる。原子ｇ_０、ｇ_４の位置が、それぞれ原子ｉ、ｊの位置と極めて近い場合も同じく処理する。

また、密な領域の原子同士を結合していたバネ、及び疎な領域の原子同士を結合していたバネをすべて消去する。代わりに、隣り合う原子同士がすべてバネ定数ｋ_ｉｊのバネで結合すると考える。具体的には、原子ｉと原子ｇ_１、原子ｇ_１と原子ｇ_２、原子ｇ_２と原子ｇ_３、及び原子ｇ_３と原子ｊとはすべてバネ定数ｋ_ｉｊのバネで結合すると考える。

図２（Ｃ）を参照して、接着後、原子ｉと原子ｊとを結ぶ方向に直交する方向に作用する力について説明する。

原子ｉと原子ｊとを結ぶ直線から、当該直線と直交する方向にΔｒだけ変位している原子には、変位している方向と反対方向に変位Δｒに比例する力が加わると考える。たとえば、図２（Ｃ）において、原子ｉと原子ｊとを結ぶ直線と直交する方向にΔｒだけ変位している原子ｇ_１には、

・・（１）

の力が作用するとする。ここで負の符号は力Ｆの作用する向きが変位Δｒと反対方向であることを示す。力Ｆは原子ｉと原子ｊとを結ぶ直線に直交する方向に働く。

また、式（１）の力Ｆの反作用として、原子ｉに作用する力Ｆ_ｉ、原子ｊに作用する力Ｆ_ｊをそれぞれ以下の式（２）、（３）で導入する。

・・（２）

・・（３）

ここでｌ_ｉ、ｌ_ｊはそれぞれ原子ｇ_１の存する点から直線ｉｊに下ろした垂線の足と、原子ｉ、原子ｊとの間の距離を示す。力Ｆ_ｉ、Ｆ_ｊの働く向きは、力Ｆの働く向きと反対方向である。

力Ｆ、Ｆ_ｉ、Ｆ_ｊを式（１）〜（３）で示すように定めることで、図示の系内で力とトルクのバランスが保たれる。

このようにして、原子が密に配置された２次元領域と、疎に配置された２次元領域とを接着する。

本願発明者らは、上述のように原子密度の異なる２つの２次元領域を接着して、精度確認のためのシミュレーションを行った。

図３〜図７を参照して、シミュレーションについて説明する。

図３（Ａ）及び（Ｂ）に、精度確認のための計算系を示す。

図３（Ａ）を参照する。半径６ｍｍの半円筒と、半径２０ｍｍの半円筒の円筒面を、お互いの頂点部分が接触するように接触させ、半径６ｍｍの半円筒側から接触点に垂直に、５００Ｎ／ｍｍの荷重をかけるシミュレーションを行った。

シミュレーションに当たっては、２つの半円筒を形成する材料にＳＵＳを想定し、ヤング率を２０８ＧＰａ、密度を７８００ｋｇ／ｍ^３、ポアソン比を０．２７に設定した。

図３（Ｂ）は、２つの半円筒の接触部分、及び粗密両領域が接着された部分の拡大写真である。シミュレーションにおいては、２つの半円筒の接触部分の近傍にのみ、密領域を導入し、２つの半円筒の接近距離、及び接近距離と理論値との間の誤差を求めた。なお、上記条件で接触された２つの円筒が接近する距離の理論値は１６．８４μｍである。

図４（Ａ）及び（Ｂ）を参照する。まず、比較例として２つの半円筒に粗密領域を設けず、メッシュ幅を一様としてシミュレーションを行った。メッシュ幅は半円筒を構成する原子の原子間距離に相当する。

図４（Ａ）を参照して、第１の比較例について説明する。第１の比較例においては、半径６ｍｍの半円筒のメッシュ幅を０．３ｍｍとし、半径２０ｍｍの半円筒のメッシュ幅を０．５ｍｍとしてシミュレーションを行った。

第１の比較例において、２つの半円筒の接近距離は１７．７８μｍであり、理論値との誤差は、５．５８％であった。

図４（Ｂ）を参照して、第２の比較例について説明する。第２の比較例においては、半径６ｍｍの半円筒のメッシュ幅を０．０８ｍｍとし、半径２０ｍｍの半円筒のメッシュ幅を０．１ｍｍとしてシミュレーションを行った。

第２の比較例において、２つの半円筒の接近距離は１７．２０μｍであり、理論値との誤差は、２．１４％であった。

第１の比較例と第２の比較例とを比べたとき、メッシュ幅の小さい第２の比較例の方が理論値との誤差が小さいことがわかる。

図５を参照して、基本シミュレーションについて説明する。基本シミュレーションにおいては、疎領域については、第１の比較例と等しいメッシュ幅を用い、密領域については、第２の比較例と等しいメッシュ幅を用いた。すなわち、半径６ｍｍの半円筒については、疎領域のメッシュ幅を０．３ｍｍ、密領域のメッシュ幅を０．０８ｍｍとし、半径２０ｍｍの半円筒については、疎領域のメッシュ幅を０．５ｍｍ、密領域のメッシュ幅を０．１ｍｍとして、２つの半円筒の接近距離を計算した。その結果、２つの半円筒の接近距離は１７．４８μｍ、理論値との誤差は３．２１％であった。

全体のメッシュ幅を一様に細かくした第２の比較例には及ばないが、メッシュ幅が一様に粗い第１の比較例に比べると、理論値との誤差が小さくなっていることがわかる。

次に本願発明者らは、上述の基本シミュレーションよりも、疎領域のメッシュ幅を粗くしてシミュレーションを行った。

図６（Ａ）を参照する。本図に示す例では、半径６ｍｍの半円筒の疎領域のメッシュ幅を０．６ｍｍ、半径２０ｍｍの半円筒の疎領域のメッシュ幅を１．０ｍｍとした。その結果、２つの半円筒の接近距離は１７．５０μｍ、理論値との誤差は３．９２％であった。

図６（Ｂ）を参照する。本図に示すのは、図６（Ａ）に示す例より更に疎領域のメッシュ幅を粗くした例である。本例においては、半径６ｍｍの半円筒の疎領域のメッシュ幅を１．２ｍｍ、半径２０ｍｍの半円筒の疎領域のメッシュ幅を２．０ｍｍとした。２つの半円筒の接近距離は１７．３３μｍ、理論値との誤差は２．９１％という結果が得られた。

図６（Ａ）及び（Ｂ）を参照して説明したシミュレーションと、基本シミュレーションとの比較から、誤差は、疎領域のメッシュ幅の大小にかかわらず、ほぼ等しくなることがわかる。

最後に、本願発明者らは、理論値との誤差とアスペクト比との関係について調査した。アスペクト比は、疎領域のメッシュ幅を密領域のメッシュ幅で除したものとして定義される。

調査に当たっては、図３（Ａ）及び（Ｂ）を参照して説明した計算系において、２つの半円筒の疎領域のメッシュ幅をともに０．５ｍｍに固定し、密領域のメッシュ幅を様々に変化させることでアスペクト比を変化させた。そして変化させたアスペクト比を与えるメッシュ幅ごとに、２つの半円筒の接近距離、及び接近距離と理論値との間の誤差を求めた。

図７は、アスペクト比と誤差との関係を示すグラフである。グラフの横軸はアスペクト比を任意単位で表し、縦軸は、誤差を単位「％」で表す。

アスペクト比が大きくなるほど誤差は小さくなっている。これはアスペクト比を大きくするために、密領域のメッシュ幅を小さくしたためであると考えられる。

グラフから、実施例による接着手法は、アスペクト比を大きくすることによる不都合は生じないものと考えられる。

次に、図８〜図１３を参照し、３次元の場合について、実施例による粗密領域接着手法を説明する。

図８（Ａ）に原子が密に配置された３次元領域と、原子が疎に配置された３次元領域を示す。

図８（Ｂ）に示すように、これら２つの領域を接着して、１つの領域を形成する。この際、密な領域の内部と疎な領域の内部は、それぞれポテンシャル修正法で相互作用のパラメータを決定しておく。

図９（Ａ）に、密に配置された領域の原子と、疎に配置された領域の原子の接着前の状態を示す。

図の左側の原子ｇ_０〜ｇ_５は密に配置された領域に属する原子である。原子ｇ_０、ｇ_２、ｇ_５が形成する三角形の原子ｇ_０、ｇ_２を結ぶ辺の中点に原子ｇ_１が配置されている。また、原子ｇ_２、ｇ_５を結ぶ辺の中点に原子ｇ_４が、原子ｇ_５、ｇ_０を結ぶ辺の中点に原子ｇ_３が配置されている。

図の右側の原子ｉ、ｊ、及びｋは疎に配置された領域に属する原子である。原子ｉ、ｊ、ｋが形成する三角形は、密領域の原子ｇ_０、ｇ_２、ｇ_５が形成する三角形と合同である。

原子ｉとｊとはバネ定数ｋ_ｉｊのバネで結合しているとする。同様に、原子ｊとｋとはバネ定数ｋ_ｊｋのバネで、原子ｋとｉとはバネ定数ｋ_ｋｉのバネで結合しているとする。

疎な領域の原子ｉ、ｊ、ｋのつくる三角形上またはその内部に、密な領域の原子ｇ_０〜ｇ_５を接着する。

図９（Ｂ）を参照する。接着に当たって、密な領域の原子位置と疎な領域の原子位置とが一致した場合、密な領域の原子を取り除き、代わりに疎な領域の原子を用いる。

図９（Ｂ）においては、それぞれ疎な領域の原子ｉ、ｊ、ｋと位置が一致する密な領域の原子ｇ_０、ｇ_２、ｇ_５を取り除き、その位置には原子ｉ、ｊ、ｋを残す。原子ｇ_０、ｇ_２、ｇ_５の位置が、それぞれ原子ｉ、ｊ、ｋの位置と極めて近い場合も同じく処理する。

図９（Ｃ）を参照する。疎な領域の原子同士を結合していたバネをすべて消去し、代わりに、隣り合う原子同士を結ぶバネを想定する。原子ｉ、ｊを結ぶ辺に平行な方向に隣り合う原子間を結ぶバネのバネ定数はすべてｋ_ｉｊとする。同様に、原子ｊ、ｋを結ぶ辺に平行な方向に隣り合う原子間を結ぶバネのバネ定数はすべてｋ_ｊｋとし、原子ｋ、ｉを結ぶ辺に平行な方向に隣り合う原子間を結ぶバネのバネ定数はすべてｋ_ｋｉとする。

三次元形状にメッシュを切る場合、図９（Ａ）〜（Ｃ）に示すように、接着箇所の粗密両領域の原子配置を相似形とすることが望ましい。しかしながら実際には、そうすることは困難であり、全く異なった形状に配置された粗密領域の原子同士を接着する場合が多い。そこでそのような場合には、以下のルールで結合させる。以下のルールを用いることにより、疎密両領域のメッシュを相似にしなくても接着が可能となる。

図１０を参照する。

原子ｉ、ｊ、ｋを疎な領域に属する原子とする。原子ｉｊ間の距離をｒ_ｉｊ、原子ｊｋ間の距離をｒ_ｊｋ、原子ｋｉ間の距離をｒ_ｋｉとする。また、原子ｉとｊとはバネ定数ｋ_ｉｊのバネで結合しているとする。同様に、原子ｊとｋとはバネ定数ｋ_ｊｋのバネで、原子ｋとｉとはバネ定数ｋ_ｋｉのバネで結合しているとする。

疎な領域の原子ｉ、ｊ、ｋで形成される三角形上またはその内部（三角形要素上）に、密な領域の原子ｇを接着する。原子ｉ、ｊ、ｋで形成される三角形は、原子ｇを囲む、Ｄｅｌａｕｎａｙ三角形分割によって得られた三角形である。

なお、Ｄｅｌａｕｎａｙ三角形分割とは、隣接するＶｏｒｏｎｏｉ点同士をつないで領域を分割することをいう。Ｖｏｒｏｎｏｉ点とは、Ｖｏｒｏｎｏｉ多面体（２次元ならば垂直二等分線、３次元ならば垂直二等分面で囲まれた領域）を構成する点である。

このとき原子ｇは、原子ｉ、ｊ、ｋとそれぞれ下の式（４）〜（６）で表されるバネ定数ｋ_ｇｉ、ｋ_ｇｊ、ｋ_ｇｋのバネで結合されると考える。

・・（４）

・・（５）

・・（６）

ここで、ｒ_ｇｉ、ｒ_ｇｊ、ｒ_ｇｋはそれぞれ、原子ｇと原子ｉ、原子ｇと原子ｊ、原子ｇと原子ｋ間の距離を表す。またＮは、原子ｉ、ｊ、ｋで形成される三角形上またはその内部に存在する密領域の原子（疎領域に接着される原子）の数である。更に、ｒｒ、ｋｋはそれぞれ下の式（７）、（８）で求められる値である。

・・（７）

・・（８）

すなわちｒｒは、原子ｉ、ｊ、ｋの原子間距離ｒ_ｉｊ、ｒ_ｊｋ、ｒ_ｋｉの平均値であり、ｋｋは原子ｉ、ｊ、ｋ間を結ぶバネのバネ定数ｋ_ｉｊ、ｋ_ｊｋ、ｋ_ｋｉの平均値である。

次に、図１１を参照して、接着後、原子ｉ、ｊ、ｋで定められる平面に垂直な方向に作用する力について説明する。疎領域原子ｉ、ｊ、ｋで形成される三角形上またはその内部に接着された密領域の原子ｇが、原子ｉ、ｊ、ｋで定められる平面からΔｒだけ離れた場合、原子ｇ、ｉ、ｊ、ｋに働く、平面に垂直な方向の力を、それぞれＦ_ｇ、Ｆ_ｉ、Ｆ_ｊ、Ｆ_ｋとすると、それらはそれぞれ下式（９）〜（１２）で表されると考える。

・・（９）

・・（１０）

・・（１１）

・・（１２）

ここで、Ｋは比例定数であり、（９）式において負符号は変位Δｒの方向とは反対向きに力Ｆ_ｇが作用することを示す。また、ｌ_１、ｌ_２、ｌ_３はそれぞれ、原子ｇから、原子ｉ、ｊ、ｋで定められる平面に下ろした垂線の足Ｈから、原子ｊ、ｋ、ｉまでの距離である。更に、θ_１は、原子ｉと垂線の足Ｈを結ぶ直線と、原子ｊと垂線の足Ｈを結ぶ直線とのなす角度であり、θ_２は、原子ｉと垂線の足Ｈを結ぶ直線と、原子ｋと垂線の足Ｈを結ぶ直線とのなす角度である。また、δは下式（１３）で計算される値である。

・・（１３）

力Ｆ_ｇ、Ｆ_ｉ、Ｆ_ｊ、Ｆ_ｋを式（９）〜（１２）で示すように定めることで、図示の系内で力とトルクのバランスが保たれる。

このようにして、原子が密に配置された３次元領域と、疎に配置された３次元領域とを接着する。

本願発明者らは、２次元の粗密領域の接着の場合と同様に、原子密度の異なる２つの３次元領域を接着して、精度確認のためのシミュレーションを行った。

図１２及び図１３を参照して、シミュレーションについて説明する。

図１２（Ａ）及び（Ｂ）に、精度確認のための計算系を示す。

図１２（Ａ）を参照する。直径３０ｍｍの球と直径１３ｍｍの球とを、２球の中心を結ぶ直線に沿って、相対速度１００ｍ／ｓで衝突させるシミュレーションを行い、２球の接触時間τからヤング率を導出し、設定値との比較を行った。シミュレーションにおいて、ヤング率は２０８ＧＰａに設定した。なお、接触時間τとヤング率Ｅとの関係は、下式（１４）〜（１６）を用いて表される。

・・（１４）

・・（１５）

・・（１６）

式（１４）においてμは換算質量、νは２球の相対速度を示す。また、式（１５）において、Ｒ及びＲ´は、２球それぞれの半径を示す。更に、式（１６）におけるσはポアソン比である。

シミュレーションに当たっては、２つの球を形成する材料にＳＵＳを想定し、密度を７８００ｋｇ／ｍ^３、ポアソン比を０．２７に設定した。また、先に述べたように、ヤング率を２０８ＧＰａとした。このとき式（１４）〜（１６）を用いて、接触時間τの理論値は、約２１．６１μｍと計算される。

図１２（Ｂ）に、２つの球の中心を通る断面を示す。シミュレーションにおいては、２つの球の表面近傍に密領域を導入し、それ以外の球の内部領域（中心部）に疎領域を画定して、上述のルールで接着した。

シミュレーションは、疎密両領域のアスペクト比が２、３、４、５、１０の各場合について行った。なお、比較例としてアスペクト比１の場合（領域の密度が一様とされ、疎密領域の接着が行われない場合）のシミュレーションも実施した。

図１３（Ａ）及び（Ｂ）にシミュレーションの結果を示す。図１３（Ａ）は、アスペクト比と接触時間τとの関係を表し、図１３（Ｂ）は、アスペクト比とヤング率Ｅとの関係を表す。両図ともに、横軸はアスペクト比を任意単位で示す。また、図１３（Ａ）における縦軸は、接触時間τを単位「μｓ」で示し、図１３（Ｂ）における縦軸は、ヤング率Ｅを単位「ＧＰａ」で示す。

アスペクト比が大きくなるにつれ、接触時間τは短く（図１３（Ａ））、ヤング率Ｅは大きくなる（図１３（Ｂ））ことがわかる。しかしアスペクト比が１０以下であるならば、たとえばヤング率Ｅは設定値から５％以内の誤差で求められることが確認される。

図１４は、実施例によるシミュレーション方法を示すフローチャートである。

まず、ステップＳ１０１において、１つの物体を密な領域と疎な領域とに分ける。

そして、ステップＳ１０２において、密な領域、疎な領域のそれぞれに粒子を配置する。すなわち、疎な領域を、疎に分布する複数の粒子で表し、密な領域を、密に分布する複数の粒子で表す。

次に、ステップＳ１０３において、疎密それぞれの領域に対してポテンシャル修正法を用い、相互作用のパラメータを決定する。すなわち、疎な領域、密な領域の各々の粒子の力学的相互作用を定義する。

続いて、ステップＳ１０４において、疎な領域の表面に存在する粒子と密な領域の表面に存在する粒子との間の力学的相互作用を定義する。

具体的には、ステップＳ１０４ａにおいて、疎密それぞれの領域が２次元領域である場合、疎な領域の表面を構成するライン上に存在する、密な領域の粒子を探索する。疎密両域が３次元領域である場合は、疎な領域の表面を構成する三角形要素上に存在する、密な領域の粒子を探索する。

そして、ステップＳ１０４ｂにおいて、２次元の場合はライン上、３次元の場合は三角形要素上にある密な領域の粒子と、疎な領域の粒子とを、図１〜図１３を参照して説明した接着手法を用いて接着（Ｇｌｕｅ）し、相互作用のパラメータを決定する。

最後に、ステップＳ１０５において、分子動力学計算（粒子法によるシミュレーション）を実行する。ステップＳ１０３及びＳ１０４において定義された力学的相互作用に基づいて、密な領域及び疎な領域内の粒子の挙動をシミュレーションする。

実施例によるシミュレーション方法を用いると、たとえば接触問題を含む分子シミュレーションを行うに当たって、物体同士が接触する箇所には密に粒子を配置し、それ以外の箇所には疎に粒子を配置し、これらを接着（Ｇｌｕｅ）することで、全体の粒子数を少なくし、シミュレーションの効率を高めることができる。

従来のポテンシャルを修正する手法だけでは、モデルに配置する粒子間距離（メッシュの大きさ）は、一番小さい領域に律速してしまい、特に３次元解析を行う場合、極端に粒子数が多くなってしまって、実質的には計算を行うことが不可能であった。しかしながら、実施例によるシミュレーション方法によれば、接着（Ｇｌｕｅ）の考え方を用いることにより、２次元解析であれば１０分の１、３次元解析であれば１００分の１程度の粒子数で、接着（Ｇｌｕｅ）をしない場合と同等の結果を得ることが可能である。

図１４にフローチャートを示したシミュレーション方法は、プログラムにより、コンピュータで実行することができる。

図１５は、実施例によるシミュレーション方法を実行するシミュレーション実行装置のシステム構成図である。

まず、キーボードなどの入力装置から、シミュレーション対象（物体）に関する情報、たとえば大きさ、形状、材料などの情報が入力される。そしてそのシミュレーション対象（物体）を、どのように疎な領域と密な領域に分けるかを決定する情報、たとえば画定される疎密両領域が２次元領域か、３次元領域か、それらの領域の形状や大きさ（範囲）、更に、それぞれの領域における粒子の粒子間距離等の、疎領域及び密領域に関する情報が入力される。

中央処理装置は、メインメモリ中の制御プログラムの指令を受け、メインメモリに保存されている相互作用パラメータ決定プログラムを実行する。

相互作用パラメータ決定プログラムは、入力された情報を基に、前図Ｓ１０１に示したように、１つの物体を密な領域と疎な領域とに分ける。そしてＳ１０２に示したように、密な領域、疎な領域のそれぞれに粒子を配置する位置を決定する。すなわち、疎な領域を、疎に分布する複数の粒子で表し、密な領域を、密に分布する複数の粒子で表すための演算を行う。

そしてＳ１０３に示したように、疎密それぞれの領域に対してポテンシャル修正法を用い、相互作用のパラメータを決定する。すなわち、疎な領域、密な領域の各々の中での粒子の力学的相互作用を定義する。

更に、ステップＳ１０４に示したように、疎な領域の表面に存在する粒子と密な領域の表面に存在する粒子との間の力学的相互作用を定義する。具体的には、ステップＳ１０４ａに示したように、疎密それぞれの領域が２次元領域である場合、疎な領域の表面を構成するライン上に存在する、密な領域の粒子を探索する計算を行う。疎密両域が３次元領域である場合は、疎な領域の表面を構成する三角形要素上に存在する、密な領域の粒子を探索する計算を行う。そして、Ｓ１０４ｂに示したように、２次元の場合はライン上、３次元の場合は三角形要素上にある密な領域の粒子と、疎な領域の粒子とを、図１〜図１３を参照して説明した接着手法を用いて接着（Ｇｌｕｅ）し、相互作用のパラメータを決定する。

その後、中央処理装置は、メインメモリ中の制御プログラムの指令を受け、メインメモリに保存されているシミュレーション実行プログラムを実行し、Ｓ１０５に示したように、粒子法によるシミュレーションを実行する。

シミュレーションの結果は、ディスプレイやプリンタなどの出力装置を介して表示される。

以上実施例に沿って本発明を説明したが、本発明はこれらに制限されるものではない。種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。

原理的には、量子的な現象を除くすべての古典現象を解析するシミュレーションに利用することができる。特に、機構、弾性の動解析における摩擦による発熱、破壊の計算に好適に用いることができる。

（Ａ）及び（Ｂ）は、２次元疎密領域の接着について示す図である。 （Ａ）は、密に配置された領域の原子と、疎に配置された領域の原子の接着前の状態を示し、（Ｂ）は、両領域の原子の接着後の状態を示し、（Ｃ）は、接着後、原子ｉと原子ｊとを結ぶ方向に直交する方向に作用する力を示す図である。 （Ａ）及び（Ｂ）は、精度確認のための計算系を示す図である。 （Ａ）及び（Ｂ）は、それぞれ第１、第２の比較例について説明するための図である。 基本シミュレーションについて説明するための図である。 （Ａ）及び（Ｂ）は、基本シミュレーションよりも、疎領域のメッシュ幅を粗くして行ったシミュレーションについて説明するための図である。 アスペクト比と誤差との関係を示すグラフである。 （Ａ）及び（Ｂ）は、３次元疎密領域の接着について示す図である。 （Ａ）は、密に配置された領域の原子と、疎に配置された領域の原子の接着前の状態を示し、（Ｂ）及び（Ｃ）は、両領域の原子の接着後の状態を示す。 密に配置された領域の原子と、疎に配置された領域の原子の接着後の状態を示す図である。 疎密両領域の接着後、疎領域原子ｉ、ｊ、ｋで定められる平面に垂直な方向に作用する力について説明するための図である。 （Ａ）及び（Ｂ）は、精度確認のための計算系を示す図である。 （Ａ）は、アスペクト比と接触時間τとの関係を表す図であり、（Ｂ）は、アスペクト比とヤング率Ｅとの関係を表す図である。 実施例によるシミュレーション方法を示すフローチャートである。 実施例によるシミュレーション方法を実行するシミュレーション実行装置のシステム構成図である。

符号の説明

Ｓ１０１〜Ｓ１０５、Ｓ１０４ａ、Ｓ１０４ｂ ステップ

Claims (6)

1. （ａ）シミュレーション対象の１つの物体を密な領域と疎な領域とに分ける工程と、
   （ｂ）前記密な領域を、第１の粒子間距離以下の粒子間距離で分布する複数の粒子で表し、前記疎な領域を、前記第１の粒子間距離よりも大きい粒子間距離で分布する複数の粒子で表す工程と、
   （ｃ）前記密な領域、前記疎な領域の各々の中で、粒子位置に基づいてポテンシャルを修正することにより粒子の力学的相互作用を定義する工程と、
   （ｄ）前記疎な領域との境界における前記密な領域の表面に存在する粒子と、前記密な領域との境界における前記疎な領域の表面に存在する粒子との間の力学的相互作用を、前記密な領域と前記疎な領域との境界に存在する粒子からなる系内で力とトルクのバランスが保たれるように定義する工程と、
   （ｅ）前記工程（ｃ）及び（ｄ）において定義された力学的相互作用に基づいて、前記密な領域及び前記疎な領域内の粒子の挙動をシミュレーションする工程と
   を有するシミュレーション方法。
2. 前記工程（ｄ）において、前記疎な領域の表面に接触した前記密な領域の粒子と、該粒子が接触した前記疎な領域の表面位置を囲むＤｅｌａｕｎａｙ三角形分割によって得られた三角形の頂点に位置する、前記疎な領域の表面の３つの粒子との間に、力学的相互作用を定義する請求項１に記載のシミュレーション方法。
3. 前記工程（ｂ）において、前記疎な領域との境界における前記密な領域の粒子の位置と、前記密な領域との境界における前記疎な領域の粒子の位置とが一致するように粒子を分布させ、
   前記工程（ｄ）において、前記密な領域の粒子と前記疎な領域の粒子の位置が一致しているとき、前記密な領域の粒子を消滅させて、該位置に前記疎な領域の粒子のみを配置する請求項１または２に記載のシミュレーション方法。
4. 種々のシミュレーションを行うためにコンピュータを、
   シミュレーション対象の１つの物体を密な領域と疎な領域とに分けるための情報を入力
   する手段、
   前記密な領域を、第１の粒子間距離以下の粒子間距離で分布する複数の粒子で表し、前記疎な領域を、前記第１の粒子間距離よりも大きい粒子間距離で分布する複数の粒子で表すための演算を行う手段、
   前記密な領域、前記疎な領域の各々の中で、粒子位置に基づいてポテンシャルを修正することにより粒子の力学的相互作用を定義する第１の定義手段、
   前記疎な領域との境界における前記密な領域の表面に存在する粒子と、前記密な領域との境界における前記疎な領域の表面に存在する粒子との間の力学的相互作用を、前記密な領域と前記疎な領域との境界に存在する粒子からなる系内で力とトルクのバランスが保たれるように定義する第２の定義手段、
   前記第１及び第２の定義手段で定義された力学的相互作用に基づいて、前記密な領域及び前記疎な領域内の粒子の挙動をシミュレーションする手段
   として機能させるためのプログラム。
5. 前記第２の定義手段は、前記疎な領域の表面に接触した前記密な領域の粒子と、該粒子が接触した前記疎な領域の表面位置を囲むＤｅｌａｕｎａｙ三角形分割によって得られた三角形の頂点に位置する、前記疎な領域の表面の３つの粒子との間に、力学的相互作用を定義する請求項４に記載のプログラム。
6. 前記第２の定義手段は、前記密な領域の粒子と前記疎な領域の粒子の位置が一致しているとき、前記密な領域の粒子を消滅させて、該位置に前記疎な領域の粒子のみを配置する請求項４または５に記載のプログラム。