(57)【特許請求の範囲】
【請求項１】
分子動力学法を使用して多数の粒子を含む系について粒子の運動方程式を数値的に演算することによりその系を解析する解析方法であって、
粒子に働く非保存力であって前記系のエネルギを散逸させる前提で得られる非保存力を取得するステップと、
取得された非保存力を含む粒子の運動方程式を数値的に演算するステップと、
演算結果に基づいて、前記系が定常状態に達したか否かを判定するステップと、を含み、
前記非保存力を取得するステップは、粒子の速度に依存した非保存力を取得するステップを含むことを特徴とする解析方法。
【請求項２】
前記非保存力を取得するステップは、粒子に働く非保存力であって前記系の内部の振動を減衰させる前提で得られる非保存力を取得するステップを含むことを特徴とする請求項１に記載の解析方法。
【請求項３】
前記非保存力を取得するステップは、粒子に働く非保存力であって前記系の並進運動を減衰させる前提で得られる非保存力を取得するステップを含むことを特徴とする請求項１または２に記載の解析方法。
【請求項４】
分子動力学法を使用して多数の粒子を含む系について粒子の運動方程式を数値的に演算することによりその系を解析する解析装置であって、
粒子に働く非保存力であって前記系のエネルギを散逸させる前提で得られる非保存力を取得する非保存力取得部と、
前記非保存力取得部によって取得された非保存力を含む粒子の運動方程式を数値的に演算する演算部と、
前記演算部における演算結果に基づいて、前記系が定常状態に達したか否かを判定する定常状態判定部と、を備え、
前記非保存力取得部は、粒子の速度に依存した非保存力を取得することを特徴とする解析装置。
【請求項５】
分子動力学法を使用して多数の粒子を含む系について粒子の運動方程式を数値的に演算することによりその系を解析する機能をコンピュータに実現させるコンピュータプログラムであって、
粒子に働く非保存力であって前記系のエネルギを散逸させる前提で得られる非保存力を取得する機能と、
取得された非保存力を含む粒子の運動方程式を数値的に演算する機能と、
演算結果に基づいて、前記系が定常状態に達したか否かを判定する機能と、を前記コンピュータに実現させ、
前記取得する機能は、粒子の速度に依存した非保存力を取得することを特徴とするコンピュータプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、系を解析する解析方法および解析装置に関する。
【背景技術】
【０００２】
一定の力やトルクを入力した系がどのような状態でつりあい位置に収束するかを知るには、有限要素法による陰解法が広く用いられている。しかしながら、特に多数の接触を伴う場合等は極めて収束が悪くなり、特に物体が剛体運動を伴う場合、計算が破綻することも珍しくない。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】越塚  誠一、「弾性体に対する陰的な粒子計算モデルの開発」、計算工学講演会論文集、１９９９、４、ｐ．３７−４０
【発明の概要】
【発明が解決しようとする課題】
【０００４】
そこで近年、粒子法（例えば、非特許文献１参照）の陽解法による解析が行われ始めているが、そのままではエネルギ保存則より振動が残留してしまう。
【０００５】
本発明はこうした課題に鑑みてなされたものであり、その目的は、陽解法において効率よくエネルギを散逸し、系の定常状態を素早く得るための手法の提供にある。
【課題を解決するための手段】
【０００６】
本発明のある態様は解析方法に関する。この解析方法は、多数の粒子を含む系について粒子の運動方程式を数値的に演算することによりその系を解析する解析方法であって、粒子に働く非保存力であって系のエネルギを散逸させる前提で得られる非保存力を取得するステップと、取得された非保存力を含む粒子の運動方程式を数値的に演算するステップと、演算結果に基づいて、系が定常状態に達したか否かを判定するステップと、を含む。
【０００７】
この態様によると、系のエネルギを散逸させる前提で得られる非保存力を考慮することで、系のエネルギを効率よく散逸できる。
【０００８】
なお、以上の構成要素の任意の組み合わせや、本発明の構成要素や表現を装置、方法、システム、コンピュータプログラム、コンピュータプログラムを格納した記録媒体などの間で相互に置換したものもまた、本発明の態様として有効である。
【発明の効果】
【０００９】
本発明によれば、陽解法において効率よくエネルギを散逸できる。
【図面の簡単な説明】
【００１０】
【図１】実施の形態に係る解析処理を説明するためのフローチャートである。
【図２】解析装置の機能および構成を示すブロック図である。
【図３】図３（ａ）、（ｂ）は、実施の形態に係る解析方法の機構解析への応用例を示す説明図である。
【図４】各種エネルギ散逸過程を取り入れた場合とそうでない場合とで運動エネルギの減衰を示すグラフである。
【発明を実施するための形態】
【００１１】
粒子法による陽解法による解析では振動が残留するので、定常状態を得るためには何らかの方法でエネルギを除去する必要がある。このとき、例えばある一定の間隔で運動エネルギをゼロにするなどの方法も考えられるが、問題によっては適切でない場合もある。
【００１２】
一般的な有限要素法では、陰解法が主流であり、陽解法は主に非定常な動解析で用いられる場合が多い。最近では静的陽解法という手法も提案されているが、これは力の釣り合い位置を探索しながら時間発展させる方法である。これらのことから、有限要素法ではエネルギ散逸をさせるという発想そのものがないと言える。
【００１３】
粒子法は主に陽解法で解かれるが，解析の目的が破壊や摩擦などの動的なものが多いため、やはりエネルギを散逸させるような発想は無い。片持ち梁の曲げ位置を求める解析を行った例もある（例えば、非特許文献１参照）が、この場合は陰的なアルゴリズムを用いている。
【００１４】
ＤＥＭ(Distinct Element Method)では衝突時の減衰係数を経験的に決定しているが、これは反発係数を模擬しているのであり、定常状態を得る目的で用いられているのではない。
【００１５】
以下、本発明を好適な実施の形態をもとに図面を参照しながら説明する。各図面に示される同一または同等の構成要素、部材、処理には、同一の符号を付するものとし、適宜重複した説明は省略する。
【００１６】
実施の形態に係る解析方法では、粒子法の陽解法において、解析対象の系のエネルギの散逸を考慮に入れてその系の定常状態を求める。これにより、より素早く系の定常状態を得ることができる。
【００１７】
図１は、実施の形態に係る解析処理を説明するためのフローチャートである。図１に示される処理は、実施の形態に係る解析装置１０（図２参照）により実行される。解析装置１０はまず、解析モデルを生成する（Ｓ１０２）。解析装置１０は、ユーザからの入力情報に基づいて解析モデルを作成してもよいし、既に作成されている解析モデルを読み込むことによりモデルを取得してもよい。この解析モデルは、分子動力学法に適する粒子モデルであり、Ｎ個の粒子からなる粒子系Ｓとして構成される。粒子は原子または分子であってもよい。
粒子系ＳではＮ個の粒子のそれぞれが初期配置され、またＮ個の粒子のそれぞれに初期速度が与えられる。
【００１８】
解析装置１０は、粒子系Ｓの各粒子に働く保存力を取得する（Ｓ１０４）。解析装置１０は、ユーザからの入力情報に基づいて粒子に働く保存力を取得してもよい。あるいはまた、解析装置１０は、粒子系Ｓにおいて粒子間の相互作用を表すポテンシャルエネルギ関数（以下では適宜、ポテンシャル関数または単にポテンシャルとも称する）に基づいて、粒子に働く保存力を演算して取得してもよい。
【００１９】
解析装置１０は、粒子に働く非保存力であって粒子系Ｓのエネルギを散逸させる前提で得られる非保存力を取得する（Ｓ１０６）。非保存力を取得するステップＳ１０６では、解析装置１０は、粒子に働く非保存力のうち、粒子系Ｓの内部の振動を減衰させる前提で得られる散逸力と、粒子系Ｓの並進運動を減衰させる前提で得られる粘性力と、を取得する（Ｓ１０８、Ｓ１１０）。粒子系Ｓ内部の粒子ｉ、ｊがバネ定数ｋで結合されていると仮定する場合の、粒子ｉ、ｊの一方が他方に及ぼす散逸力は以下の式１から計算され、粒子系Ｓ内部の粒子ｉに働く粘性力は以下の式２から計算される。
【００２０】
【数１】
0005546339.tif 000002
但し、ｍｉは粒子ｉの質量、ｍｊは粒子ｊの質量、ｖｉは粒子ｉの速度、ｖｊは粒子ｊの速度、ｔは時間、ｖｒは粒子ｉ、ｊの相対速度、
【数２】
0005546339.tif 000003
はｎを任意の整数として（ｎ−１／２）ｄｔにおける粒子ｉ、ｊの相対速度、Ｍは粒子ｉ、ｊの相対運動における換算質量、である。
【００２１】
【数３】
0005546339.tif 000004
但し、
【数４】
0005546339.tif 000005
はｎを任意の整数として（ｎ−１／２）ｄｔにおける粒子ｉの速度、τは粘性抵抗の大きさによって定まる時定数、である。
式１および式２からも理解されるとおり、散逸力や粘性力は粒子の速度に依存した力である。
【００２２】
解析装置１０は、各粒子について散逸力と粘性力と保存力とを合算し、その粒子に働く合計の力（以下、合算力と称す）を得る（Ｓ１１２）。
解析装置１０は、得られた合算力を含む粒子の運動方程式を数値的に解くことにより、その粒子の位置および速度のうちの少なくともひとつを演算する。解析装置１０はその演算結果に基づき粒子を移動させる（Ｓ１１４）。
【００２３】
解析装置１０は、粒子の位置および速度のうちの少なくともひとつの演算結果に基づき、粒子系Ｓが定常状態に達したか否かを判定する（Ｓ１１６）。解析装置１０は、例えば粒子系Ｓの全運動エネルギが所定の下限値より小さい場合、粒子系Ｓが定常状態に達したと判定する。解析装置１０は、粒子系Ｓが定常状態に達していないと判定された場合（Ｓ１１６のＮ）、ステップＳ１０４に処理を戻す。解析装置１０は、粒子系Ｓが定常状態に達していると判定された場合（Ｓ１１６のＹ）、分子動力学計算処理を終了し、定常状態における粒子系Ｓの物理量をユーザに提供する。
【００２４】
図２は、解析装置１０の機能および構成を示すブロック図である。ここに示す各ブロックは、ハードウエア的には、コンピュータのＣＰＵ（central processing unit）をはじめとする素子や機械装置で実現でき、ソフトウエア的にはコンピュータプログラム等によって実現されるが、ここでは、それらの連携によって実現される機能ブロックを描いている。したがって、これらの機能ブロックはハードウエア、ソフトウエアの組合せによっていろいろなかたちで実現できることは、本明細書に触れた当業者には理解されるところである。
【００２５】
解析装置１０は、演算部１２と、記憶部２４と、入力部２６と、出力部２８と、を備える。演算部１２は、モデル生成部１４と、保存力取得部１６と、非保存力取得部１８と、合算部２０と、分子動力学演算部２２と、定常状態判定部３０と、を含む。
【００２６】
モデル生成部１４は、入力部２６を介したユーザからの入力情報または記憶部２４に既に記憶されている情報に基づき、粒子系ＳのＮ個の粒子の初期位置および初期速度を設定する。
保存力取得部１６は粒子系Ｓの各粒子に働く保存力を取得する。保存力取得部１６は、粒子系Ｓにおいて粒子間の相互作用を表すポテンシャルのグラジエント（Gradient）を演算し、演算されたグラジエントに応じた保存力を取得する。
【００２７】
非保存力取得部１８は、式１に基づき散逸力を演算して取得する。非保存力取得部１８は、式２に基づき粘性力を演算して取得する。以下、式１および式２の導出過程を説明する。
【００２８】
（散逸力、式１）
散逸力は、粒子系Ｓの内部の振動を減衰させる前提で導出される。粒子系Ｓの粒子ｉ、ｊが，バネ定数kで結合されている場合を考える．散逸方程式は以下の式３で与えられる。
【数５】
0005546339.tif 000006
但し、γ（ｔ）は散逸の強さを示す係数である。
【００２９】
ここで、以下の式４の通りに相対速度ｖｒと重心速度ｖｇを定義し、重心の運動と相対運動に分ける。
【数６】
0005546339.tif 000007
【００３０】
このとき、式３の散逸方程式は、
【数７】
0005546339.tif 000008
となる。相対運動の解は、
【数８】
0005546339.tif 000009
となる。
【００３１】
式５、式６から、散逸過程においても，重心の運動量は保存され，相対運動のエネルギが散逸されることが分かる。以下、これらに対し数値積分を実行する。
式６を離散化すると、
【数９】
0005546339.tif 000010
となる。但し、
【数１０】
0005546339.tif 000011
はｎを任意の整数として（ｎ＋１／２）ｄｔにおける粒子ｉ、ｊの相対速度、γｎ−１／２はｎを任意の整数として（ｎ−１／２）ｄｔにおけるγの値、である。
式８をまとめると、
【数１１】
0005546339.tif 000012
となる。
【００３２】
式３から、粒子ｉに対する離散化された散逸方程式は
【数１２】
0005546339.tif 000013
となるので、式９を式１０に代入すると、
【数１３】
0005546339.tif 000014
が得られる。
【００３３】
ここで、粒子系Ｓを振動させず過減衰させるためには、
【数１４】
0005546339.tif 000015
とすればよい。式１２を式１１に代入することで、式１のうちｖｉについての式が得られる。式１のうちｖｊについての式も同様にして得られる。
【００３４】
この離散化は相対速度の値として前ステップのものを用いているため、部分的に陰解法となっており（過去の情報のみ使用している）、極端にバネ定数ｋが大きい場合などを除いて破綻しにくい特徴がある。
【００３５】
（粘性力、式２）
粘性力は、粒子系Ｓの並進運動を減衰させる前提で導出される。この並進運動の減衰は粘性抵抗による減衰として表される。
粒子系Ｓの内部の粒子ｉに働く粘性抵抗に係る方程式（以下、粘性方程式と称す）は、以下の式１３で与えられる。
【数１５】
0005546339.tif 000016
但し、γ（ｔ）は粘性の強さを示す係数である。
式１３を解くことにより、
【数１６】
0005546339.tif 000017
を得る。
【００３６】
式１４より、時定数τは、
【数１７】
0005546339.tif 000018
で与えられる。
このτの値によって粘性抵抗の大きさが決まる。本発明者の当業者としての経験から、τは１０００ｄｔ〜１００００ｄｔ程度の大きさが適している。
【００３７】
式８、式９、式１０と同様の手順により、
【数１８】
0005546339.tif 000019
が得られる。したがって、式１５を式１６に代入することで式２が得られる。式１および式２を使用することで、振動・並進のエネルギを効率良く散逸することが可能である。
【００３８】
合算部２０は、粒子ｉについて、他の粒子ｊ（ｊ＝１、２、…、Ｎ（ｊ≠ｉ））との相互作用を表すポテンシャルのグラジエントに応じた保存力ＦＣｉｊを保存力取得部１６から取得してｊについて合算する。粒子ｉについて、このように合算された保存力を合算保存力ΣＦＣｉと称す。
合算部２０は、粒子ｉについて、他の粒子（ｊ＝１、２、…、Ｎ（ｊ≠ｉ））から及ぼされる散逸力ＦＤｉｊを非保存力取得部１８から取得してｊについて合算する。粒子ｉについて、このように合算された散逸力を合算散逸力ΣＦＤｉと称す。
合算部２０は、粒子ｉについて、粘性力ＦＶｉを非保存力取得部１８から取得する。
合算部２０は、粒子ｉについて、合算保存力ΣＦＣｉと合算散逸力ΣＦＤｉと粘性力ＦＶｉとを合算して、粒子ｉに働く合算力ＦＴｉ（＝ΣＦＣｉ＋ΣＦＤｉ＋ＦＶｉ）を得る。
【００３９】
分子動力学演算部２２は、合算部２０で得られた合算力ＦＴｉを含む粒子ｉの運動方程式を数値積分する。なお、数値積分には、Ｖｅｒｌｅｔ法が用いられることが多い。Ｖｅｒｌｅｔ法については、例えば、Ｊ．Ｍ．Ｔｈｉｊｓｓｅｎ，“Ｃｏｍｐｕｔａｔｉｏｎａｌ  Ｐｈｙｓｉｃｓ”，Ｃａｍｂｒｉｄｇｅ  Ｕｎｉｖｅｒｓｉｔｙ  Ｐｒｅｓｓ（１９９９）のｐ．１７５に説明されている。
【００４０】
分子動力学演算部２２は、粒子ｉの位置および速度を、運動方程式を数値積分することによって得られる粒子ｉの新たな位置および速度で更新する。なお、位置および速度の両方を更新することは必須ではなく、いずれか一方のみを更新するようにしてもよい。
【００４１】
定常状態判定部３０は、粒子系Ｓの全運動エネルギが下限値より小さい場合、粒子系Ｓが定常状態に達したと判定する。定常状態判定部３０は、粒子系Ｓが定常状態に達していないと判定した場合、処理を保存力取得部１６に戻し、分子動力学演算部２２で更新された各粒子の位置および速度を使用した新たな分子動力学演算を開始させる。定常状態判定部３０は、粒子系Ｓが定常状態に達していると判定した場合、分子動力学演算を終了する。
【００４２】
記憶部２４は、実施の形態に係る解析方法を実施するための分子動力学計算プログラムや粒子モデル生成プログラム等を格納する。また、記憶部２４は、これらプログラムの実行に要する各種のデータを格納する。演算部１２は記憶部２４から必要なプログラム及びデータを読み込んで処理を実行する。
【００４３】
入力部２６は、演算部１２で実行される処理に関係するユーザの入力を受けるためのキーボード、マウス等の入力装置を含む。入力部２６は、フレキシブルディスク、ＣＤ、ＤＶＤ等の記録媒体から入力を受けるよう構成されていてもよい。
出力部２８は、ディスプレイ等の表示機器やプリンタ等を含む出力装置である。ユーザは、出力部２８を介して、定常状態に達した粒子系Ｓの応力などの各種物理量を確認する。出力部２８は、記録媒体に結果を出力するよう構成されていてもよい。
【００４４】
図３（ａ）、（ｂ）は、実施の形態に係る解析方法の機構解析への応用例を示す説明図である。図３（ａ）は、トルクを与えて回転する機構５０と固定剛体ピン６０との初期状態を示す説明図である。機構５０は、軸部５２と歯部５４とを有する。本応用例では、軸部５２に図３（ａ）のようにトルクが付加され、同時に機構５０全体が固定剛体ピン６０に向けて動かされる。機構５０と固定剛体ピン６０とが十分近づくと機構５０の歯部５４と固定剛体ピン６０とが接触する。本応用例では、定常状態、すなわち接触後の釣り合い位置を決定する。
【００４５】
このとき、エネルギを散逸させないと、全エネルギが保存することから、機構５０は固定剛体ピン６０との間で衝突を繰り返し続けることになる。これを避けるために、定期的に運動エネルギを強制的にゼロにする場合を考える。
【００４６】
この場合、機構５０は固定剛体ピン６０に当たるまでの距離を移動する必要があるため、一旦運動エネルギをゼロにしてしまうとまた加速からやり直さなければならず、効率が悪い。また，モデルによってはエネルギ散逸のタイミングと振動の周期が共振してしまい、振動が大きくなる場合もある。
【００４７】
そこで、本実施の形態に係る解析方法のように、系のエネルギを散逸させる前提で得られる非保存力を導入する。この場合で、非保存力として散逸力を導入する場合、機構５０は速やかに固定剛体ピン６０に当たるので、シミュレーションにかかる時間を短縮できる。また、散逸力を導入することで定常状態に至るまでの時間を十分に低減できる場合は、粘性力と散逸力との両方を考慮する場合と比べて計算負荷を低減できるので好適である。
ただし、散逸力の場合では、機構５０が固定剛体ピン６０に当たった後、幾度か固定剛体ピン６０との衝突を繰り返すことになりうる。その際、次第に衝突間隔は短くなりいずれは釣り合い位置に達すると考えられる。しかしながらその釣り合い位置を得るためには、衝突が終わるのを待つ必要がある。
【００４８】
非保存力として粘性力を導入する場合、やはり機構５０は速やかに固定剛体ピン６０に当たるので、シミュレーションにかかる時間を短縮できる。また、粘性力を導入することで定常状態に至るまでの時間を十分に低減できる場合は、粘性力と散逸力との両方を考慮する場合と比べて計算負荷を低減できるので好適である。
ただし、固定剛体ピン６０との衝突後，その衝撃で内部に発生した振動はなかなか散逸されないと予想される。仮に内部振動がすぐに散逸するほど粘性抵抗が大きいようだと、運動エネルギを定期的にゼロとする場合と同様に固定剛体ピン６０に衝突するまでの時間が無視できなくなる可能性がある。
【００４９】
散逸力と粘性力の両方を導入する場合、衝突の繰り返しを抑えつつ衝突の衝撃で内部に発生した振動も効果的に減衰できるので、機構５０のエネルギを効率的に散逸して素早く釣り合い位置を得ることができる。
【００５０】
図３（ｂ）は、機構５０と固定剛体ピン６０とが釣り合い位置となったときのＭｉｓｅｓ応力の分布を示す説明図である。本応用例では、実施の形態に係る解析方法に、図３（ａ）の初期状態を与えて図３（ｂ）の定常状態を素早く得ることができ、シミュレーションにかかる時間を短縮できる。
【００５１】
図４は、各種エネルギ散逸過程を取り入れた場合とそうでない場合とで運動エネルギの減衰を示すグラフである。図４の縦軸は機構５０の運動エネルギを任意の単位で示し、横軸は時間を任意の単位で示す。図４の符号７０で示されるデータは、図３に示される応用例において、散逸力も粘性力も取り入れない場合に対応する。図４の符号７２で示されるデータは、図３に示される応用例において、散逸力を取り入れた場合に対応する。図４の符号７４で示されるデータは、図３に示される応用例において、散逸力および粘性力を取り入れた場合に対応する。図４から明らかなとおり、本実施の形態に係る解析方法によると、散逸力および粘性力を分子動力学演算に導入することで、好適に系のエネルギを散逸させてより速く定常状態を得ることができる。
【００５２】
上述の実施の形態において、記憶装置の例は、ハードディスクやメモリである。また、本明細書の記載に基づき、各部を、図示しないＣＰＵや、インストールされたアプリケーションプログラムのモジュールや、システムプログラムのモジュールや、ハードディスクから読み出したデータの内容を一時的に記憶するメモリなどにより実現できることは本明細書に触れた当業者には理解されるところである。
【００５３】
以上、実施の形態に係る解析方法および解析装置１０の構成と動作について説明した。これらの実施の形態は例示であり、その各構成要素や各処理の組み合わせにいろいろな変形例が可能なこと、またそうした変形例も本発明の範囲にあることは当業者に理解されるところである。
【００５４】
図３および図４では、機構５０と固定剛体ピン６０からなる系に本実施の形態に係る解析方法を適用する場合について説明したが、これに限られず、例えば負荷印加時の歯車の接触解析にも応用できる。特に、歯車などの接触箇所が多い場合の構造・機構解析について、接触箇所が多いと陰的な解析では収束解が得られにくいが、本実施の形態に係る解析方法は陽的な解法であり必ず解が見つかるので好適である。
【００５５】
実施の形態では、非保存力として散逸力や粘性力を採用する場合について説明したが、これに限られず、例えば摩擦力を採用してもよい。
【符号の説明】
【００５６】
１０  解析装置、  １２  演算部、  １４  モデル生成部、  １６  保存力取得部、  １８  非保存力取得部、  ２０  合算部、  ２２  分子動力学演算部、  ２４  記憶部、  ２６  入力部、  ２８  出力部、  ３０  定常状態判定部。
【図１】
0005546339.tif 000020
【図２】
0005546339.tif 000021
【図４】
0005546339.tif 000022
【図３】
0005546339.tif 000023