特許 7117973 シミュレーション装置、シミュレーション方法、及びプログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】

(24)【登録日】2022-08-04

(45)【発行日】2022-08-15

(54)【発明の名称】シミュレーション装置、シミュレーション方法、及びプログラム

(51)【国際特許分類】

   G16Z 99/00 20190101AFI20220805BHJP

   G16C 10/00 20190101ALI20220805BHJP

【ＦＩ】

G16Z99/00

G16C10/00

【請求項の数】 5

(21)【出願番号】P 2018204194

(22)【出願日】2018-10-30

(65)【公開番号】P2020071602

(43)【公開日】2020-05-07

【審査請求日】2021-06-16

(73)【特許権者】

【識別番号】000002107

【氏名又は名称】住友重機械工業株式会社

(74)【代理人】

【識別番号】100105887

【弁理士】

【氏名又は名称】来山 幹雄

(72)【発明者】

【氏名】大西 良孝

(72)【発明者】

【氏名】宮崎 修司

【審査官】宮地 匡人

(56)【参考文献】

【文献】特開２０１７－１８２２３１（ＪＰ，Ａ）

【文献】特開２００６－０１５８６６（ＪＰ，Ａ）

【文献】特開２００９－０３７３３４（ＪＰ，Ａ）

【文献】米国特許出願公開第２００２／０１５７４７８（ＵＳ，Ａ１）

【文献】大山 尚武，力学系はいかにして熱平衡に達するか－非線形格子振動の計算機実験－，応用物理，Vol.38 No.7，1969年，pp.689-695

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ｇ１６Ｚ ９９／００

Ｇ１６Ｃ １０／００

ＪＳＴＰｌｕｓ（ＪＤｒｅａｍＩＩＩ）

ＪＳＴ７５８０（ＪＤｒｅａｍＩＩＩ）

(57)【特許請求の範囲】

【請求項1】

シミュレーション対象の材料に応じて決められている粒子間の相互作用ポテンシャルの線形項を定義するパラメータの値、非線形項を定義するパラメータの値、及び粒子の配置の初期条件を取得する入力部と、
前記入力部で取得された前記パラメータの値で定義される相互作用ポテンシャルを用い、前記入力部で取得された初期条件に基づいて分子動力学法によって粒子の挙動を解析する処理部と
を有し、
前記入力部は、前記パラメータの初期値、及びシミュレーション対象の材料の応力と歪みとの実測値から得られた第１応力歪み関係データを取得し、
前記処理部は、
前記パラメータの値を初期値から更新しながら、前記パラメータの値の組み合わせごとに分子動力学法によって応力と歪みとの関係を表す第２応力歪み関係データを求め、
前記第１応力歪み関係データと前記第２応力歪み関係データとの比較結果に基づいて、前記パラメータの値の最適値を決定するシミュレーション装置。

【請求項2】

前記処理部は、さらに、
前記第１応力歪み関係データからヤング率と歪みとの関係を表す第１ヤング率歪み関係データを求め、
前記第２応力歪み関係データからヤング率と歪みとの関係を表す第２ヤング率歪み関係データを求め、
前記第１ヤング率歪み関係データと前記第２ヤング率歪み関係データとの比較結果に基づいて、前記パラメータの最適値を決定する請求項１に記載のシミュレーション装置。

【請求項3】

シミュレーション条件として、シミュレーション対象の材料に応じて決められている粒子間の相互作用ポテンシャルの線形項を定義するパラメータの値、及び非線形項を定義するパラメータの値で定義される相互作用ポテンシャルを用いて分子動力学法によって粒子の挙動を解析するシミュレーション方法であって、
前記粒子の挙動を解析する前に、さらに、前記パラメータの値を初期値から更新しながら、前記パラメータの値の組み合わせごとに分子動力学法によって応力と歪みとの関係を表す第２応力歪み関係データを求め、
シミュレーション対象の材料の応力と歪みとの実測値から得られた第１応力歪み関係データと、前記第２応力歪み関係データとの比較結果に基づいて、前記パラメータの最適値を決定し、
前記パラメータの最適値を用いて、前記粒子の挙動を解析するシミュレーション方法。

【請求項4】

さらに、
前記第１応力歪み関係データからヤング率と歪みとの関係を表す第１ヤング率歪み関係データを求め、
前記第２応力歪み関係データからヤング率と歪みとの関係を表す第２ヤング率歪み関係データを求め、
前記第１ヤング率歪み関係データと前記第２ヤング率歪み関係データとの比較結果に基づいて、前記パラメータの最適値を決定する請求項３に記載のシミュレーション方法。

【請求項5】

シミュレーション条件として、シミュレーション対象の材料に応じて決められている粒子間の相互作用ポテンシャルの線形項を定義するパラメータの値及び非線形項を定義するパラメータの値で定義される相互作用ポテンシャルを用いて分子動力学法によって粒子の挙動を解析する機能、
前記粒子の挙動を解析する前に、さらに、前記パラメータの値を初期値から更新しながら、前記パラメータの値の組み合わせごとに分子動力学法によって応力と歪みとの関係を表す第２応力歪み関係データを求める機能、
シミュレーション対象の材料の応力と歪みとの実測値から得られた第１応力歪み関係データと、前記第２応力歪み関係データとの比較結果に基づいて、前記パラメータの最適値を決定する機能、及び
前記粒子の挙動を解析する際に、前記パラメータの最適値を用いて解析する機能
をコンピュータに実現させるプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

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

【背景技術】

【0002】

任意の形状の材料に外力が加わったときの変形を、分子動力学法またはくりこみ群分子動力学法により解析するシミュレーション方法が公知である（例えば、特許文献１）。以下、本明細書において分子動力学法及びくりこみ群分子動力学法を合わせて、単に分子動力学法という。

【0003】

分子動力学法によるシミュレーション方法では、シミュレーション対象の材料を複数の粒子の集合体で表し、粒子同士をバネで結合したモデルを用いて弾性領域の解析を行うことができる。粒子同士を結合するバネのバネ定数は、複数の粒子の位置を節点とするテトラメッシュの情報から導出されるボロノイ面積に基づいて、粒子の配置によらず当該材料のヤング率を再現するように定められる。

【先行技術文献】

【特許文献】

【0004】

【文献】特開２００９－３７３３４号公報

【発明の概要】

【発明が解決しようとする課題】

【0005】

線形弾性材料の解析では、公知の方法で、粒子を結合するバネのバネ定数を決定することができる。この方法でバネ定数を決定する従来のシミュレーション方法では、非線形弾性材料の変形等のシミュレーションを行うことができない。

【0006】

本発明の目的は、非線形弾性材料の変形等のシミュレーションを行うことができるシミュレーション装置、シミュレーション方法、及びプログラムを提供することである。

【課題を解決するための手段】

【0007】

本発明の一観点によると、
シミュレーション対象の材料に応じて決められている粒子間の相互作用ポテンシャルの線形項を定義するパラメータの値、非線形項を定義するパラメータの値、及び粒子の配置の初期条件を取得する入力部と、
前記入力部で取得された前記パラメータの値で定義される相互作用ポテンシャルを用い、前記入力部で取得された初期条件に基づいて分子動力学法によって粒子の挙動を解析する処理部と
を有し、
前記入力部は、前記パラメータの初期値、及びシミュレーション対象の材料の応力と歪みとの実測値から得られた第１応力歪み関係データを取得し、
前記処理部は、
前記パラメータの値を初期値から更新しながら、前記パラメータの値の組み合わせごとに分子動力学法によって応力と歪みとの関係を表す第２応力歪み関係データを求め、
前記第１応力歪み関係データと前記第２応力歪み関係データとの比較結果に基づいて、前記パラメータの値の最適値を決定するシミュレーション装置が提供される。

【0008】

本発明の他の観点によると、
シミュレーション条件として、シミュレーション対象の材料に応じて決められている粒子間の相互作用ポテンシャルの線形項を定義するパラメータの値、及び非線形項を定義するパラメータの値で定義される相互作用ポテンシャルを用いて分子動力学法によって粒子の挙動を解析するシミュレーション方法であって、
前記粒子の挙動を解析する前に、さらに、前記パラメータの値を初期値から更新しながら、前記パラメータの値の組み合わせごとに分子動力学法によって応力と歪みとの関係を表す第２応力歪み関係データを求め、
シミュレーション対象の材料の応力と歪みとの実測値から得られた第１応力歪み関係データと、前記第２応力歪み関係データとの比較結果に基づいて、前記パラメータの最適値を決定し、
前記パラメータの最適値を用いて、前記粒子の挙動を解析するシミュレーション方法が提供される。

【0009】

本発明のさらに他の観点によると、
シミュレーション条件として、シミュレーション対象の材料に応じて決められている粒子間の相互作用ポテンシャルの線形項を定義するパラメータの値及び非線形項を定義するパラメータの値で定義される相互作用ポテンシャルを用いて分子動力学法によって粒子の挙動を解析する機能、
前記粒子の挙動を解析する前に、さらに、前記パラメータの値を初期値から更新しながら、前記パラメータの値の組み合わせごとに分子動力学法によって応力と歪みとの関係を表す第２応力歪み関係データを求める機能、
シミュレーション対象の材料の応力と歪みとの実測値から得られた第１応力歪み関係データと、前記第２応力歪み関係データとの比較結果に基づいて、前記パラメータの最適値を決定する機能、及び
前記粒子の挙動を解析する際に、前記パラメータの最適値を用いて解析する機能
をコンピュータに実現させるプログラムが提供される。

【発明の効果】

【0010】

非線形項を定義するパラメータを含む相互作用ポテンシャルを用いることにより、非線形弾性材料の変形等のシミュレーションを行うことが可能になる。

【図面の簡単な説明】

【0011】

【図1】図１は、実施例によるシミュレーション装置のブロック図である。

【図2】図２は、粒子ｉ、粒子ｊ、及び両者の近傍のボロノイ多面体の一部を示す模式図である。

【図3】図３は、相互作用ポテンシャルの形状を示すグラフである。

【図4】図４Ａはシミュレーションモデルを示す斜視図であり、図４Ｂは、印加した応力と、シミュレーションにより得られた歪みとの関係を示すグラフである。

【図5】図５は、寸法の異なる３つのシミュレーションモデルの斜視図である。

【図6】図６Ａ及び図６Ｂは、非線形パラメータであるａ_３＝ａ_４＝０とした相互作用ポテンシャルを用いて解析を行ったときの各モデルの歪みと応力との関係を示すグラフであり、図６Ｃは、ａ_３＝－５、ａ_４＝－１とした場合のシミュレーション結果を示すグラフであり、図６Ｄは、ａ_３＝－１０、ａ_４＝２０とした場合のシミュレーション結果を示すグラフである。

【図7】図７は、パラメータＡ、ａ_３、ａ_４の最適値を見つけ出す手順のフローチャートである。

【図8】図８Ａは、第１応力歪み関係データの一例を示すグラフであり、図８Ｂは、第１応力歪み関係データから求められるヤング率と歪みεとの関係を示すグラフである。

【図9】図９は、ヤング率歪み関係の決定係数を横軸とし、応力歪み関係の決定係数を縦軸とする散布図である。

【図10】図１０Ａは、パラメータＡ、ａ_３、ａ_４を最適値に設定して分子動力学法により求めた第２応力歪み関係データと、実測値から求めた第１応力歪み関係データとを示すグラフであり、図１０Ｂは、パラメータＡ、ａ_３、ａ_４を最適値に設定して分子動力学法により求めた第２ヤング率歪み関係データと、実測値から求めた第１ヤング率歪み関係データとを示すグラフである。

【図11】図１１は実施例によるシミュレーション方法のフローチャートである。

【発明を実施するための形態】

【0012】

図１～図１１を参照して、実施例によるシミュレーション装置及びシミュレーション方法について説明する。

【0013】

図１は、実施例によるシミュレーション装置のブロック図である。実施例によるシミュレーション装置は、入力部２０、処理部２１、出力部２２、及び記憶部２３を含む。入力部２０から処理部２１にシミュレーション条件等が入力される。さらに、オペレータから入力部２０に各種指令（コマンド）等が入力される。入力部２０は、例えば通信装置、リムーバブルメディア読取装置、キーボード等で構成される。

【0014】

処理部２１は、入力されたシミュレーション条件に基づいて分子動力学法によるシミュレーションを行い、シミュレーション結果を出力部２２に出力する。シミュレーション結果には、シミュレーション対象を複数の粒子の集合体で表したときの粒子の挙動を表す情報が含まれる。処理部２１は、例えばコンピュータを含み、分子動力学法によるシミュレーションをコンピュータに実行させるためのプログラムが記憶部２３に記憶されている。出力部２２は、通信装置、リムーバブルメディア書込み装置、ディスプレイ等を含む。

【0015】

本実施例では、非線形弾性材料の変形をシミュレーションするために、分子動力学法で用いる粒子間の相互作用ポテンシャルに非線形項を導入する。以下、相互作用ポテンシャルの非線形項の導入方法について説明する。

【0016】

粒子ｉと粒子ｊとの間の相互作用ポテンシャルφ（ｒ_ｉｊ）は、以下の式で表される。

【数1】

ここで、ｋ_ｉｊは、粒子ｉと粒子ｊとを結合するバネのバネ定数であり、ｒ_ｉｊは、粒子ｉと粒子ｊとの距離であり、ｒ_ｉｊ０は、バネが自然長のときの粒子ｉと粒子ｊとの距離である。

【0017】

次に、図２を参照してバネ定数ｋ_ｉｊの決定方法について説明する。バネ定数ｋ_ｉｊの決定方法については上記特許文献１に詳しく説明されているため、ここでは簡単に説明する。まず、シミュレーション対象の三次元形状に基づいてテトラメッシュを生成する。このテトラメッシュの節点に粒子を配置し、ボロノイ多面体解析を行う。

【0018】

図２は、粒子ｉ、粒子ｊ、及び両者の近傍のボロノイ多面体の一部を示す模式図である。ボロノイ多面体を構成する複数の界面のうち粒子ｉと粒子ｊとを両端とする線分Ｌ_ｉｊを横切る界面の面積をＳ_ｉｊで表す。線形弾性材料の場合には、材料のヤング率と面積Ｓ_ｉｊとからバネ定数ｋ_ｉｊを決定することができる。

【0019】

本実施例では、式（１）の相互作用ポテンシャルに非線形項を加えて相互作用ポテンシャルを以下の式（２）で定義する。

【数2】

式（２）では、非線形項として３次及び４次の項を追加しているが、さらに高次の項まで追加してもよい。本実施例では、相互作用ポテンシャルの非線形項として４次の項までを考える。

【0020】

式（２）の３次の項及び４次の項の係数ａ_ｉｊ３及びａ_ｉｊ４を、以下の式（３）で定義する。

【数3】

【0021】

扱う材料が線形弾性材料である場合には、図２を参照して説明したように材料のヤング率からバネ定数ｋ_ｉｊを決定することができる。ところが、非線形弾性材料を扱う場合には、ヤング率を明確に定義することができないため、ヤング率からバネ定数を決定することもできない。そこで、新たに調整可能なパラメータＡを加えて、バネ定数をＡ・ｋ_ｉｊと定義することとする。

【0022】

さらに、非線形弾性材料を圧縮したときの変形は、ほぼ線形であると仮定し、非線形弾性材料に引っ張り応力を加えたときに非線形性が現れると仮定する。
非線形弾性材料に縮み歪を発生させたとき、すなわちｒ_ｉｊ≦ｒ_ｉｊ０のとき、相互作用ポテンシャルは以下の式（４）で表すことができる。

【数4】

非線形弾性材料に延び歪を発生させたとき、すなわちｒ_ｉｊ＞ｒ_ｉｊ０のとき、相互作用ポテンシャルは以下の式（５）で表すことができる。

【数5】

【0023】

Ａ＝１、ａ_３＝ａ_４＝０のとき、式（５）の相互作用ポテンシャルは、線形弾性材料を扱うときの相互作用ポテンシャルと同じ形状になる。

【0024】

図３は、相互作用ポテンシャルの形状を示すグラフである。横軸はバネの自然長ｒ_ｉｊ０で正規化した粒子間距離ｒ_ｉｊ／ｒ_ｉｊ０を表し、縦軸は相互作用ポテンシャルの大きさを表す。図３のグラフに示した太い実線はａ_３＝ａ_４＝０の場合、細い実線はａ_３＝－５、ａ_４＝－１の場合、破線はａ_３＝－１０、ａ_４＝２０の場合の相互作用ポテンシャルを示す。材料に延び歪みが発生している領域において非線形性が現れていることがわかる。

【0025】

式（４）及び式（５）の相互作用ポテンシャルを用いて引張試験を行ったときの材料の変形をシミュレーションにより求めた。以下、このシミュレーション結果について説明する。

【0026】

図４Ａはシミュレーションモデルを示す斜視図である。シミュレーションモデルは、縦、横、高さをそれぞれ１０ｍｍ、１０ｍｍ、２０ｍｍの直方体とした。この直方体の底面を固定し、上面に引っ張り応力を印加するという条件で、分子動力学法により定常状態になるまで解析し、歪みを求めた。

【0027】

図４Ｂは、印加した応力と、シミュレーションにより得られた歪みとの関係を示すグラフである。横軸は歪みを表し、縦軸は応力を単位「ＧＰａ」で表す。図４Ｂのグラフに示した太い実線は式（５）においてａ_３＝ａ_４＝０とした場合、細い実線はａ_３＝－５、ａ_４＝－１とした場合、破線はａ_３＝－１０、ａ_４＝２０とした場合の応力と歪みとの関係を示す。ａ_３＝ａ_４＝０のときには、歪みと応力との関係が線形である。その他の場合には、歪みと応力との間に非線形の関係が現れていることが確認された。

【0028】

式（５）に示した相互作用ポテンシャルを用いてシミュレーションを行う場合、シミュレーションモデルの寸法が異なっても、求められる歪みと応力との関係が同一でなければならない。以下、寸法の異なるシミュレーションモデルの歪みと応力との関係を求めた結果について説明する。

【0029】

図５は、寸法の異なる３つのシミュレーションモデルの斜視図である。いずれのシミュレーションモデルも直方体とした。モデル１の縦、横、高さは、それぞれ１０ｍｍ、１０ｍｍ、２０ｍｍである。モデル２の縦、横、高さは、それぞれ１５ｍｍ、１０ｍｍ、４０ｍｍである。モデル３の縦、横、高さは、それぞれ３０ｍｍ、２０ｍｍ、５０ｍｍである。これらのシミュレーションモデルの底面を固定し、上面に引っ張り応力を印加して、歪みと応力との関係を求めた。式（５）のＡ・ｋ_ｉｊは、ヤング率が２０８ＧＰａの材料を再現するように設定した。

【0030】

図６Ａ及び図６Ｂは、式（５）においてａ_３＝ａ_４＝０とした相互作用ポテンシャルを用いて解析を行ったときの各モデルの歪みと応力との関係を示すグラフである。図６Ｂは図６Ａの一部分を拡大したグラフである。参考のために、ヤング率が２０８ＧＰａの材料の理論値も一緒に示している。いずれのモデルにおいても、ほぼ理論値に近いシミュレーション結果が得られていることが確認された。

【0031】

図６Ｃは、式（５）においてａ_３＝－５、ａ_４＝－１とした場合のシミュレーション結果を示し、図６Ｄは、式（５）においてａ_３＝－１０、ａ_４＝２０とした場合のシミュレーション結果を示すグラフである。参考のために、線形弾性材料に引っ張り応力を加えたときの理論値を示している。いずれのモデルにおいても、線形弾性材料の理論値に対して非線形性が現れていることがわかる。また、モデル１、２、３のシミュレーション結果は、ほとんど区別がつかない程度に一致していることがわかる。これにより、式（５）の相互作用ポテンシャルはシミュレーションモデルの寸法に依存しないことが確認された。

【0032】

式（５）に示した相互作用ポテンシャルを用いてシミュレーションを行う場合、シミュレーション対象の材料の物性を再現するようにパラメータＡ、ａ_３、ａ_４を決定しなければならない。ところが、これらの非線形項のパラメータが応力と歪みとの関係にどのように作用するのか不明であるため、材料の物性値から直ちにこれらのパラメータの値を決定することはできない。これらのパラメータの値は、例えば試行錯誤を繰り返すことにより見つけ出すことができる。以下、これらのパラメータの最適値を見つけ出す手順の一例について説明する。なお、ここで「最適値」とは、すべての値の組み合わせのうち最適な値であることを意味しているわけではなく、十分高い精度でシミュレーションを行うことができる好ましい値を意味している。

【0033】

図７は、シミュレーション対象の材料に応じてパラメータＡ、ａ_３、ａ_４の最適値を見つけ出す手順のフローチャートである。この手順は、図１に示したシミュレーション装置の処理部２１により実行される。なお、シミュレーション装置とは異なるコンピュータでこの手順を実行するようにしてもよい。

【0034】

まず、パラメータＡ、ａ_３、ａ_４の初期値を設定する（ステップＳＡ１）。この初期値は、例えば、シミュレーション対象である非線形弾性材料の線形変形領域におけるヤング率等に基づいて、経験則から決定するとよい。パラメータＡ、ａ_３、ａ_４の設定された値に基づいて式（５）の相互作用ポテンシャルを用い、分子動力学法により応力と歪みとの関係を求める（ステップＳＡ２）。この関係を、第２応力歪み関係データということとする。

【0035】

ステップＳＡ２の計算を行っている途中で計算が破綻した場合（ステップＳＡ３）には、計算を終了し、パラメータＡ、ａ_３、ａ_４の値を更新して（ステップＳＡ６）、ステップＳＡ２の処理を繰り返す。ここで、計算が破綻する場合とは、例えばゼロで除算する計算を行う場合等を意味する。パラメータＡ、ａ_３、ａ_４の値の最大値、最小値、及び更新の際の刻み幅は、予め決められている。

【0036】

ステップＳＡ２の計算が破綻することなく第２応力歪み関係データが求められた場合には、実測による応力と歪みとの関係を表す第１応力歪み関係データと、計算によって求められた第２応力歪関係データとを比較する。両者を比較することにより、第１応力歪み関係データに対する第２応力歪関係データの誤差を求める（ステップＳＡ４）。

【0037】

図８Ａは、第１応力歪み関係データの一例を示すグラフである。横軸は歪みεを表し、縦軸は応力を単位「ＧＰａ」で表す。図８Ａに示したグラフにおいて、実測値を丸記号で示しており、実測値を二次関数で近似した二次曲線を実線で示している。応力を歪みεの関数としてｆ（ε）と表記することとする。第１応力歪み関係データは、関数ｆ（ε）で定義することができる。

【0038】

図８Ｂは、第１応力歪み関係データから求められるヤング率と歪みεとの関係を示すグラフである。横軸は歪みεを表し、縦軸はヤング率を単位「ＧＰａ」で表す。ヤング率は、関数ｆ（ε）を歪みεで微分することにより得られる。第１応力歪み関係データから求められるヤング率と歪みεとの関係を第１ヤング率歪み関係データということとする。

【0039】

第１応力歪み関係データに対する第２応力歪関係データの誤差として、例えば実測値から求まる第１応力歪み関係データを表す関数ｆ（ε）を回帰方程式とし、計算結果から求まる第２応力歪み関係データを標本値とする決定係数（Ｒ^２値）を採用するとよい。以下、この決定係数を「応力歪み関係の決定係数」という。

【0040】

さらに、計算によって求められた第２応力歪み関係データから、第２ヤング率歪み関係データを求める。第２ヤング率歪み関係データは、例えば、離散的な第２応力歪み関係データの隣り合う２つの点の、歪みεの増分に対する応力の増分の比で表すことができる。ステップＳＡ４において、第１応力歪み関係データに対する第２応力歪関係データの誤差とともに、第１ヤング率歪み関係データに対する第２ヤング率歪み関係データの誤差を求める。この誤差として、実測値から求まる第１ヤング率歪み関係データを表す関数ｄｆ（ε）／ｄε（図８Ｂ）を回帰方程式とし、計算結果から求まる第２ヤング率歪み関係データを標本値とする決定係数（Ｒ^２値）を採用するとよい。以下、この決定係数を、「ヤング率歪み関係の決定係数」という。

【0041】

ステップＳＡ２（図７）からステップＳＡ４（図７）までの処理、及びステップＳＡ６（図７）の処理を、規定計算回数に達するまで繰り返す（ステップＳＡ５）。ステップＳＡ２からステップＳＡ４までの処理、及びステップＳＡ６の処理が規定計算回数に達すると、第１応力歪み関係データと第２応力歪み関係データとの比較結果、及び第１ヤング率歪み関係データと第２ヤング率歪み関係データとの比較結果に基づいて、パラメータＡ、ａ_３、ａ_４の最適値を決定し、出力する（ステップＳＡ７）。

【0042】

以下、図９を参照してパラメータＡ、ａ_３、ａ_４の最適値を決定する方法について説明する。

【0043】

図９は、ヤング率歪み関係の決定係数を横軸とし、応力歪み関係の決定係数を縦軸とする散布図である。図９に示した散布図は、パラメータＡの最小値、最大値、刻み幅をそれぞれ０．０１、０．２、０．０２とし、パラメータａ_３の最小値、最大値、刻み幅をそれぞれ－５００、５００、１．０とし、パラメータａ_４の最小値、最大値、刻み幅をそれぞれ－５００、５００、１．０として計算した結果を示している。

【0044】

ヤング率歪み関係の決定係数、及び応力歪み関係の決定係数が最も大きくなるときのパラメータＡ、ａ_３、ａ_４の値を最適値として採用する。図９に示した例では、パラメータＡ、ａ_３、ａ_４の最適値を与える点のヤング率歪み関係の決定係数が０．９８３、応力歪み関係の決定係数が０．９９９７であり、計算結果と実測結果とのよい一致が得られた。

【0045】

図１０Ａは、パラメータＡ、ａ_３、ａ_４を最適値に設定して分子動力学法により求めた第２応力歪み関係データと、実測値から求めた第１応力歪み関係データとを示すグラフである。三角記号は、計算値から求めた第２応力歪み関係データを示し、実線は、実測値による第１応力歪み関係データを表す関数ｆ（ε）を示す。

【0046】

図１０Ｂは、パラメータＡ、ａ_３、ａ_４を最適値に設定して分子動力学法により求めた第２ヤング率歪み関係データと、実測値から求めた第１ヤング率歪み関係データとを示すグラフである。三角記号は、計算値から求めた第２ヤング率歪み関係データを示し、実線は、第１ヤング率歪み関係データを表す関数ｄｆ（ε）／ｄεを示す。

【0047】

図１０Ａ及び図１０Ｂに示すように、計算値と実測値とがよく一致していることがわかる。このように、図７に示した方法で、高い精度でシミュレーションを行うことが可能なパラメータＡ、ａ_３、ａ_４の最適値を決めることができることが確認された。

【0048】

図１１は実施例によるシミュレーション方法のフローチャートである。図１１に示した処理は、本実施例のシミュレーション装置の処理部２１によって実行される。

【0049】

処理部２１は、シミュレーション対象の材料に応じて決められている粒子間の相互作用ポテンシャルの線形項を定義するパラメータＡの値、バネ定数ｋ_ｉｊの値、非線形項を定義するパラメータａ_３、ａ_４の値、及び粒子の配置の初期条件等のシミュレーション条件を取得する（ステップＳＢ１）。なお、パラメータＡとバネ定数ｋ_ｉｊとを、Ａ・ｋ_ｉｊとして１つの線形パラメータとして扱ってもよい。取得したシミュレーション条件に基づいて粒子を配置してシミュレーションモデルを生成する（ステップＳＢ２）。

【0050】

パラメータＡ、バネ定数ｋ_ｉｊ、パラメータａ_３、ａ_４の入力された値に基づく式（４）及び式（５）の相互作用ポテンシャルを用いて、分子動力学法により粒子の挙動を解析する（ステップＳＢ３）。解析後、解析結果を出力部２２に出力する（ステップＳＢ４）。出力部２２には、例えば粒子の位置の変化またはテトラメッシュの形状の変化を時系列で図形として表示するとよい。

【0051】

次に、実施例の優れた効果について説明する。
本実施例により、非線形弾性材料の非線形領域を考慮した機構解析を行うことができる。図７に示したパラメータ最適化方法を用いて、種々の非線形弾性材料に応じて非線形パラメータの最適値を決定することができる。これにより、種々の非線形弾性材材料を用いた機構の解析を行うことができる。

【0052】

次に、上記実施例の変形例について説明する。
上記実施例では、図７のステップＳＡ７において、ヤング率歪み関係の決定係数及び応力歪み関係の決定係数を用いてパラメータＡ、ａ_３、ａ_４の最適値を決定したが、いずれか一方の決定係数を用いてパラメータＡ、ａ_３、ａ_４の最適値を決定してもよい。

【0053】

なお、図９に示したように、応力歪み関係の決定係数が１に近い場合でも、ヤング率歪み関係の決定係数が１より大幅に小さい場合がある。応力歪み関係の決定係数のみを用いてパラメータＡ、ａ_３、ａ_４の最適値を決定すると、ヤング率歪み関係の決定係数が小さい点を、最適値を与える点と誤認してしまう可能性が高まる。ヤング率歪み関係の決定係数と、応力歪み関係の決定係数とを併用することにより、このような誤認を回避することができる。

【0054】

上記実施例では、式（５）の相互作用ポテンシャルとして４次の非線形項までを考慮したが、５次以上の高次の非線形項を考慮してもよい。この場合には、非線形パラメータの個数が多くなる。

【0055】

上述の実施例は例示であり、本発明は上述の実施例に制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。

【符号の説明】

【0056】

２０ 入力部
２１ 処理部
２２ 出力部
２３ 記憶部

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】