特許 5897992 高分子の粘弾性計算装置、その方法及びプログラム

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

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

(11)【特許番号】5897992

(24)【登録日】2016年3月11日

(45)【発行日】2016年4月6日

(54)【発明の名称】高分子の粘弾性計算装置、その方法及びプログラム

(51)【国際特許分類】

   G01N 3/00 20060101AFI20160324BHJP

   G01N 19/00 20060101ALI20160324BHJP

【ＦＩ】

   G01N3/00 K

   G01N19/00 D

【請求項の数】6

【全頁数】18

(21)【出願番号】特願2012-127323(P2012-127323)

(22)【出願日】2012年6月4日

(65)【公開番号】特開2013-250242(P2013-250242A)

(43)【公開日】2013年12月12日

【審査請求日】2015年4月28日

(73)【特許権者】

【識別番号】000003148

【氏名又は名称】東洋ゴム工業株式会社

(74)【代理人】

【識別番号】100059225

【弁理士】

【氏名又は名称】蔦田 璋子

(74)【代理人】

【識別番号】100076314

【弁理士】

【氏名又は名称】蔦田 正人

(74)【代理人】

【識別番号】100112612

【弁理士】

【氏名又は名称】中村 哲士

(74)【代理人】

【識別番号】100112623

【弁理士】

【氏名又は名称】富田 克幸

(74)【代理人】

【識別番号】100124707

【弁理士】

【氏名又は名称】夫 世進

(74)【代理人】

【識別番号】100163393

【弁理士】

【氏名又は名称】有近 康臣

(72)【発明者】

【氏名】日野 理

【審査官】 福田 裕司

(56)【参考文献】

【文献】 国際公開第９６／０２１０７（ＷＯ，Ａ１）

【文献】 特開２０１１−０６９７８１（ＪＰ，Ａ）

【文献】 米国特許第５５５３００４（ＵＳ，Ａ）

【文献】 特開２００２−１８４９６９（ＪＰ，Ａ）

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

Ｇ０１Ｎ ３／００

Ｇ０１Ｎ １９／００

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

【請求項1】

高分子をモデル化したシステムを作成するシステム作成部と、
前記システムを所定温度で平衡化させる平衡化部と、
前記システムに与える歪みの周波数を決定する周波数決定部と、
前記システムに前記周波数で所定振幅の歪みを与えながら分子動力学シミュレーションにより熱雑音を含む応力を計算する応力計算部と、
前記熱雑音を含む応力と前記歪みから時間相関関数として歪み−応力相互相関関数を計算する相互相関関数計算部と、
前記歪み−応力相互相関関数から前記熱雑音を除去した応力の振幅および該応力と歪みとの位相差を計算して前記高分子の粘弾性に関する指標値を得る第１指標値計算部と、
を有することを特徴とする高分子の粘弾性計算装置。

【請求項2】

前記第１指標値計算部が、前記歪み−応力相互相関関数を用いて歪みと熱雑音の相互相関関数を除去した正弦曲線にフィッティングすることで前記応力の振幅および前記位相差を計算することを特徴とする請求項１記載の粘弾性計算装置。

【請求項3】

予め設定された複数の周波数のうち全ての周波数で前記粘弾性に関する指標値が計算されたか否かを判定し、全ての周波数で計算されるまで、平衡化部、周波数決定部、応力計算部、相互相関関数計算部および第１指標値計算部の処理を繰り返すように制御する判定部を備えたことを特徴とする請求項１又は２記載の粘弾性計算装置。

【請求項4】

前記応力計算部により得られた前記熱雑音を含む応力から時間相関関数として応力自己相関関数を計算する自己相関関数計算部と、
前記応力自己相関関数から前記熱雑音を除去した応力の振幅を計算して前記高分子の粘弾性に関する指標値を得る第２指標値計算部と、
前記第１指標値計算部により得られた指標値を、前記第２指標値計算部により得られた指標値と対比して、精度を評価する評価部と、
を備えたことを特徴とする請求項１〜３のいずれか１項に記載の粘弾性計算装置。

【請求項5】

高分子をモデル化したシステムを作成するシステム作成ステップと、
前記システムを所定温度で平衡化させる平衡化ステップと、
前記システムに与える歪みの周波数を決定する周波数決定ステップと、
前記システムに前記周波数で所定振幅の歪みを与えながら分子動力学シミュレーションにより熱雑音を含む応力を計算する応力計算ステップと、
前記熱雑音を含む応力と前記歪みから時間相関関数として歪み−応力相互相関関数を計算する相互相関関数計算ステップと、
前記歪み−応力相互相関関数から前記熱雑音を除去した応力の振幅および該応力と歪みとの位相差を計算して前記高分子の粘弾性に関する指標値を得る第１指標値計算ステップと、
を有することを特徴とする高分子の粘弾性計算方法。

【請求項6】

高分子の粘弾性を計算するプログラムであって、
コンピュータに、
高分子をモデル化したシステムを作成するシステム作成機能と、
前記システムを所定温度で平衡化させる平衡化機能と、
前記システムに与える歪みの周波数を決定する周波数決定機能と、
前記システムに前記周波数で所定振幅の歪みを与えながら分子動力学シミュレーションにより熱雑音を含む応力を計算する応力計算機能と、
前記熱雑音を含む応力と前記歪みから時間相関関数として歪み−応力相互相関関数を計算する相互相関関数計算機能と、
前記歪み−応力相互相関関数から前記熱雑音を除去した応力の振幅および該応力と歪みとの位相差を計算して前記高分子の粘弾性に関する指標値を得る第１指標値計算機能と、
を実現するための高分子の粘弾性計算プログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、高分子の粘弾性を計算するための計算装置、その方法及びそのプログラムに関するものである。

【背景技術】

【0002】

ゴムなどの高分子の材料開発を行う場合、粘弾性特性などの物性を予測することが求められ、計算機シミュレーションにより高分子材料の物性を予測することが行われている。

【0003】

例えば、特許文献１には、高分子材料の分子動力学シミュレーションから得られた位相空間情報と、線形応答理論を用いて高分子材料の各種物性（粘弾性特性を含む）を求めることが開示されている。しかしながら、線形応答理論は物質に与えられた微小変位に対する応答に関する物理量の計算理論であり、変位量依存した物理量の変化を扱うことができない。すなわち、線形応答理論では、微小外力が仮定されており、線形領域の感受率しか計算できない（例えば、R. Kubo and K. Tomita, Journal of Physical Society of Japan, vol.9, pp. 888-919 (1954)、Statistical Mechanics, Donald A. McQuarrie, University Science Books (2000) 第22章を参照）。

【先行技術文献】

【特許文献】

【0004】

【特許文献1】ＷＯ９７／０６４９６号

【発明の概要】

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

【0005】

本発明は、以上の点に鑑み、線形応答理論を用いることなく、分子動力学シミュレーションにより粘弾性特性を計算することができる高分子の粘弾性計算装置、その方法及びそのプログラムを提供することを目的とする。

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

【0006】

本発明に係る高分子の粘弾性計算装置は、高分子をモデル化したシステムを作成するシステム作成部と、前記システムを所定温度で平衡化させる平衡化部と、前記システムに与える歪みの周波数を決定する周波数決定部と、前記システムに前記周波数で所定振幅の歪みを与えながら分子動力学シミュレーションにより熱雑音を含む応力を計算する応力計算部と、前記熱雑音を含む応力と前記歪みから時間相関関数として歪み−応力相互相関関数を計算する相互相関関数計算部と、前記歪み−応力相互相関関数から前記熱雑音を除去した応力の振幅および該応力と歪みとの位相差を計算して前記高分子の粘弾性に関する指標値を得る第１指標値計算部と、を有するものである。

【発明の効果】

【0007】

本発明によれば、線形応答理論ではなく、システムに歪みを与えて応力応答を直接計算するので、任意の歪みに対して粘弾性を計算することができる。また、分子動力学シミュレーションにおいて応力は熱雑音の影響を強く受けるが、歪み−応力相互相関関数を計算することにより熱雑音の影響を除去することができるので、粘弾性を精度よく計算することができる。

【図面の簡単な説明】

【0008】

【図1】第１実施形態に係る計算装置のブロック図である。

【図2】第１実施形態に係る計算装置のフローチャートである。

【図3】第２実施形態に係る計算装置のブロック図である。

【図4】第２実施形態に係る計算装置のフローチャートである。

【図5】高分子モデルの一例を示す図である。

【図6】実施例１における歪みの時間変化を示すグラフである。

【図7】実施例１における応力の時間変化を示すグラフである。

【図8】実施例１における歪み−応力相関の時間変化を示すグラフである。

【図9】実施例１における歪み−応力相関のフィッティング結果を示すグラフである。

【図10】実施例１における応力自己相関の時間変化を示すグラフである。

【図11】実施例１における応力自己相関のフィッティング結果を示すグラフである。

【図12】実施例１における応力のフィッティング結果を示すグラフである。

【図13】実施例２における歪みの時間変化を示すグラフである。

【図14】実施例２における応力の時間変化を示すグラフである。

【図15】実施例２における歪み−応力相関の時間変化を示すグラフである。

【図16】実施例２における歪み−応力相関のフィッティング結果を示すグラフである。

【図17】実施例２における応力自己相関の時間変化を示すグラフである。

【図18】実施例２における応力自己相関のフィッティング結果を示すグラフである。

【図19】実施例２における応力のフィッティング結果を示すグラフである。

【図20】粘弾性計算結果として周波数と貯蔵弾性率の関係を示すグラフである。

【図21】粘弾性計算結果として周波数とｔａｎδの関係を示すグラフである。

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

【0009】

以下、本発明の実施形態について図面に基づいて説明する。

【0010】

［第１実施形態］
第１実施形態に係る高分子の粘弾性の計算装置１０は、図１に示すように、入力部１２、システム作成部１４、平衡化部１６、周波数決定部１８、応力計算部２０、相互相関関数計算部２２、第１指標値計算部２４、判定部２６、及び出力部２８を有する。

【0011】

なお、この計算装置１０は、例えば、マウスとキーボードを有する汎用のコンピュータを基本ハードウェアとして用いることでも実現することが可能である。すなわち、入力部１２、システム作成部１４、平衡化部１６、周波数決定部１８、応力計算部２０、相互相関関数計算部２２、第１指標値計算部２４、判定部２６、及び出力部２８は、上記のコンピュータに搭載されたプロセッサにプログラムを実行させることにより実現することができる。このとき、計算装置１０は、上記のプログラムをコンピュータに予めインストールすることで実現してもよいし、ＣＤ−ＲＯＭ等の記憶媒体に記憶して、又はネットワークを介して上記のプログラムを配布して、このプログラムをコンピュータに適宜インストールすることで実現してもよい。

【0012】

以下、上記各部の構成と機能について順番に説明する。

【0013】

［１］入力部１２
入力部１２は、計算対象となる高分子の構造に関する情報を取得する。取得する情報としては、システム作成部１４で高分子モデルを作成するのに必要な各種データが含まれ、例えば、高分子を構成するモノマーなどの繰り返し単位の構造、複数種の単位からなる場合にはその構成比、重合度、及び分子鎖の数などが挙げられる。また、入力部１２において入力されるデータとしては、分子動力学シミュレーションにおいて必要なシステムの温度及び圧力、並びにシステムに付与する歪みの振幅、周波数の範囲（計算を行う複数の周波数）などが含まれる。

【0014】

計算対象となる高分子としては、結晶性高分子でも、無定形（アモルファス）高分子でもよく、特に限定されないが、好ましくは無定形高分子である。無定形高分子としては、例えば、ブタジエン、イソプレン、クロロプレンなどのジエン系モノマーの単独重合体又は共重合体、アクリル酸やメタクリル酸又はそれらのエステル、酢酸ビニル、塩化ビニル、スチレン、アクリロニトリルなどのビニル系モノマーの単独重合体又は共重合体、更には、これらジエン系モノマーとビニル系モノマーとの共重合体などが挙げられ、ジエン系ゴム材料を対象とすることがより好ましい。

【0015】

［２］システム作成部１４
システム作成部１４は、入力部１２で取得した情報を用いて、高分子をモデル化したシステムを作成する。

【0016】

モデル化する手法としては、分子動力学シミュレーションを行うことが可能な各種のモデルを挙げることができ、特に限定されない。例えば、全原子モデルの他、ビーズスプリングモデルやユナイテッドアトムモデルなど、高分子を粗視化したモデルを用いることができる。ビーズスプリングモデルは、ｋｕｈｎ長に相当する程度のいくつかのモノマーユニットを１つのビーズ（セグメント）とみなし、仮想的なバネで結合してモデル化するものである。ユナイテッドアトムモデルは、水素を重原子（例えば炭素）に含めて１つの原子（質点）として取り扱う仮想原子モデルである。

【0017】

図５は、高分子のビーズスプリングモデルの一例を示した図であり、参考文献１：Dynamics of entangled linear polymer melts: A molecular-dynamics simulation, Kurt Kremer, Gray S. Grest, Journal of chemical physics, v.92, pp.5057-pp.5086 (1990)を用いて作成することができる。

【0018】

なお、モデル化する高分子が無定形高分子の場合、ガラス転移温度よりも高温にてモデル系、即ちシステムを作成する。

【0019】

［３］平衡化部１６
平衡化部１６は、上記システムを所定の温度で平衡化させる。所定の温度としては、粘弾性特性を求めたい温度に設定され、無定形高分子の場合、ガラス転移温度よりも高い温度が設定される。

【0020】

例えば、粒子数、圧力及び温度が一定の系であるＮＰＴアンサンブルの場合、温度を設定すれば、時間発展によりシステムは平衡化するので、温度設定後、所定時間刻みにて分子動力学計算を行い、ある一定時間を経過した段階で平衡化したと判断すればよい。

【0021】

［４］周波数決定部１８
周波数決定部１８は、上記システムに与える歪みの周波数を決定する。周波数は、粘弾性特性を求めたい所定の周波数が決定され、入力部１２で入力された複数の周波数の中から選択される。例えば、複数の周波数の低周波数側から順次に選択してもよく、高周波数側から順次に選択してもよい。

【0022】

［５］応力計算部２０
応力計算部２０は、上記システムに前記周波数で所定振幅の歪みを与えながら、分子動力学シミュレーションを行い、各時刻の応力を計算する。

【0023】

分子動力学シミュレーションは、例えば参考文献２（岡崎進、吉井範行著「コンピュータシミュレーションの基礎」、第２版、化学同人発行、２０１１年）に記載されたように公知である。例えば、高分子を原子モデルで構成する場合、原子の間にいくつかのポテンシャル関数を設定し、ニュートンの運動方程式に従って原子位置を徐々にずらしていき、このようにして時々刻々と変化する原子位置が統計的な平衡状態に達する状況をシミュレートするのが分子動力学シミュレーションである。

【0024】

詳細には、上記システムの温度Ｔを一定に保ちながら、下記式（１）で表される振幅Ｘ_０および周波数ωの正弦波の歪みをシステムに与える。かかる周期的歪みを与える周期の繰り返し数は特に限定されないが、通常は５〜５０回程度である。

【数1】

【0025】

かかる条件の下で分子動力学シミュレーションを行い、時間に依存する応力応答Ｂ（ｔ）を計算する。分子動力学計算で得られる応力値は熱雑音を含むものであり、すなわち、強い熱雑音Ｎ（ｔ）と真の応力応答値σ（ｔ）の和になっているので、応力応答Ｂ（ｔ）は、下記式（２）で表される。

【数2】

【0026】

ここで、十分に長い時間周期的な歪みを与えると、真の応力応答値σ（ｔ）は周波数ωの振動を行うようになるので、応力応答Ｂ（ｔ）は、応力の振幅をσ_０とし、歪みと応力の位相差をδとして、下記式（３）により表される。

【数3】

【0027】

なお、システムに付与する歪みの振幅は、特に限定されず、例えば、システムのＸ軸方向においてその平衡化状態の長さの０．１〜０．６倍とすることができる。すなわち、平衡化状態の長さを１００としたとき、振幅が１０〜６０である歪みを付与することができる。

【0028】

［６］相互相関関数計算部２２
相互相関関数計算部２２は、上記で得られた熱雑音を含む応力と歪みから、時間相関関数として歪み−応力相互相関関数を計算する。

【0029】

高分子材料モデルの粘弾性は、貯蔵弾性率と位相差δとによって定まる。ここで、貯蔵弾性率は、歪みの振幅と応力の振幅とにより下記式（４）で表される。

【数4】

【0030】

ここで、σ_０とｔａｎδは、上記式（３）を用いて正弦曲線にフィッティングすることにより求めることができる。しかしながら、分子動力学シミュレーションにおいて応力は熱雑音の影響を強く受けるので、式（３）の応力値Ｂ（ｔ）を直接用いると、熱雑音Ｎ（ｔ）により精度が悪くなる。そこで、本実施形態では、歪みと応力応答の時間相関関数、即ち、歪み−応力相互相関関数を計算して熱雑音を除去する。歪み−応力相互相関関数は、遅れ時間をＴとして、下記式（５）により与えられる。

【0031】

【数5】

【0032】

ここで、Ｃ（ｔ）は歪みと熱雑音の相互相関関数であり、通常絶対値は小さくなる。

【0033】

［７］第１指標値計算部２４
第１指標値計算部２４は、上記歪み−応力相互相関関数から、熱雑音を除去した応力の振幅σ_０および該応力と歪みとの位相差δを計算して、高分子の粘弾性に関する指標値を算出する。本実施形態では、該指標値として、上記式（４）で表される貯蔵弾性率Ｇ’（ω）と、ｔａｎδを求める。

【0034】

上記応力の振幅σ_０と位相差δを計算するに際しては、上記式（５）を用いたフィッティングを行うことが好ましい。すなわち、第１指標値計算部２４は、歪み−応力相互相関関数を用いて、歪みと熱雑音の相互相関関数Ｃ（ｔ）を除去した正弦曲線にフィッティングすることで、応力の振幅σ_０と位相差δを計算する。ここでいう正弦曲線には余弦曲線も含まれる。より詳細には、本実施形態では、下記式（６）を用いて、応力の時間変化を最小平均二乗法によりフィッティングする。

【数6】

【0035】

これにより、応力の振幅σ_０と位相差δが得られるので、更に歪みの振幅Ｘ_０を用いて、上記式（４）により貯蔵弾性率Ｇ’（ω）と、ｔａｎδを算出することができる。

【0036】

［８］判定部２６
判定部２６は、入力部１２で設定された複数の周波数のうち全ての周波数で粘弾性に関する指標値が計算されたか否かを判定し、全ての周波数で計算されるまで、平衡化部１６、周波数決定部１８、応力計算部２０、相互相関関数計算部２２および第１指標値計算部２４の処理を繰り返すように制御する。すなわち、全ての周波数で計算されていなければ、上記処理を繰り返すように制御し、全ての周波数で計算されていれば、計算を終了する。

【0037】

［９］出力部２８
出力部２８は、上記により得られた高分子の粘弾性特性を出力する。粘弾性特性の出力は、ディスプレイによって表示したり、プリンタによって印刷したりすることにより行うことができる。

【0038】

次に、本実施形態に係る計算装置１０の動作状態について、図２のフローチャートに基づいて説明する。

【0039】

ステップＳ１において、入力部１２が、計算対象となる高分子の分子構造を含む情報を取得する。

【0040】

次いで、ステップＳ２において、システム作成部１４が、取得した情報を用いて高分子をモデル化したシステムを作成する。そしてステップＳ３に進む。

【0041】

ステップＳ３において、平衡化部１６が、粘弾性特性を求めたい温度で上記システムを平衡化する。そしてステップＳ４に進む。

【0042】

ステップＳ４において、周波数決定部１８が、上記システムに与える歪みの周波数を決定する。歪みの周波数は、ステップＳ１で入力された複数の周波数の中から選択され、例えば低周波数側から順次に計算する場合、まず最初に、最も低い周波数が選択される。そしてステップＳ５に進む。

【0043】

ステップＳ５において、応力計算部２０が、上記システムに対し、ステップＳ４で決定した周波数で、所定振幅の周期的歪みを所定周期にわたり与えながら、分子動力学シミュレーションを行い、各時刻の応力応答を計算する。これにより得られる応力は、上記のように熱雑音を含むものである。

【0044】

次いで、ステップＳ６において、相互相関関数計算部２２が、上記で得られた熱雑音を含む応力と歪みから、時間相関関数として、上記式（５）で表される歪み−応力相互相関関数を計算する。そしてステップＳ７に進む。

【0045】

ステップＳ７において、第１指標値計算部２４が、上記歪み−応力相互相関関数から、熱雑音を除去した応力の振幅σ_０と位相差δを計算して、高分子の粘弾性に関する指標値として貯蔵弾性率とｔａｎδを算出する。そしてステップＳ８に進む。

【0046】

ステップＳ８において、判定部２６が、ステップＳ１で入力された複数の周波数のうち全ての周波数で粘弾性に関する指標値が計算されたか否かを判定し、計算されていなければ、ステップＳ３、Ｓ４、Ｓ５、Ｓ６及びＳ７の処理を繰り返すように制御し、全ての周波数で計算されていれば、計算を終了して、ステップＳ９に進む。

【0047】

ステップＳ９において、出力部２８が、上記により得られた高分子の粘弾性特性を出力する。

【0048】

以上よりなる本実施形態によれば、線形応答理論ではなく、システムに歪みを与えて応力応答を直接計算するので、任意の歪みに対して粘弾性を計算することができ、粘弾性の歪み依存性をシミュレートすることができる。また、分子動力学シミュレーションにおいて応力は熱雑音の影響を強く受けるが、歪み−応力相互相関関数を計算することにより熱雑音の影響を除去することができるので、与えられた歪みに対する応力応答を精密に計算することができ、正確な粘弾性特性を得ることができる。

【0049】

［第２実施形態］
第２実施形態に係る高分子の粘弾性の計算装置１０Ａは、図３に示すように、入力部１２、システム作成部１４、平衡化部１６、周波数決定部１８、応力計算部２０、相互相関関数計算部２２、第１指標値計算部２４、自己相関関数計算部３０、第２指標値計算部３２、評価部３４、判定部２６、及び出力部２８を有する。すなわち、第２実施形態に係る計算装置１０Ａは、第１実施形態に係る計算装置１０に対して、自己相関関数計算部３０、第２指標値計算部３２及び評価部３４を追加したものであり、その他の構成は共通している。以下、追加した構成部分について詳細に説明する。

【0050】

［１０］自己相関関数計算部３０
自己相関関数計算部３０は、応力計算部２０により得られた熱雑音を含む応力から時間相関関数として応力自己相関関数を計算する。応力自己相関関数は、遅れ時間をＴとして、下記式（７）により与えられる。

【0051】

【数7】

【0052】

ここで、Ｄ（ｔ）は熱雑音の自己相関関数であり、原点付近に強いピークを持ち急速に減衰する。

【0053】

［１１］第２指標値計算部３２
第２指標値計算部３２は、上記応力自己相関関数から、熱雑音を除去した応力の振幅σ_０を計算して、高分子の粘弾性に関する指標値を算出する。本実施形態では、該指標値として、上記式（４）で表される貯蔵弾性率Ｇ’（ω）を求める。

【0054】

該応力の振幅σ_０を計算するに際しては、上記式（７）を用いたフィッティングを行うことが好ましい。すなわち、第２指標値計算部３２は、応力自己相互相関関数を用いて、熱雑音の自己相互相関関数Ｄ（ｔ）を除去した正弦曲線にフィッティングすることで、応力の振幅σ_０を計算する。ここでいう正弦曲線には余弦曲線も含まれる。より詳細には、本実施形態では、下記式（８）を用いて、応力の時間変化を最小平均二乗法によりフィッティングする。

【数8】

【0055】

上記のようにＤ（ｔ）は原点付近に強いピークを持ち急速に減衰するので、フィッティングはＤ（ｔ）が十分にゼロに近づいた時刻から行うにより、真の応力応答の振幅σ_０を精度良く計算することができる。このように応力の振幅σ_０が得られるので、歪みの振幅Ｘ_０を用いて、上記式（４）により貯蔵弾性率Ｇ’（ω）を算出することができる。

【0056】

［１２］評価部３４
評価部３４は、第１指標値計算部２４により得られた指標値（以下、第１指標値）を、第２指標値計算部３２により得られた指標値（以下、第２指標値）と対比して、その精度を評価する。対比する指標値としては、応力の振幅σ_０でもよく、あるいはまた貯蔵弾性率Ｇ’（ω）でもよい。

【0057】

より詳細には、評価部３４は、第１指標値と第２指標値を対比し、両者の差が所定範囲内であれば、第１指標値計算部２４により算出される粘弾性特性の精度が確保されていると判断して、判定ステップＳ８に進むように制御する。一方、両者の差が所定範囲内でなければ、第１指標値計算部２４により算出される粘弾性特性の精度が確保されていないと判断して、応力計算ステップＳ５に進み、分子動力学シミュレーションにおけるシミュレーション条件を変更した上で再度応力を計算するように制御する。すなわち、シミュレーションの精度を高めるために、例えば、システムに付与する歪み周期の繰り返し数を増加させ、両者の差が所定範囲内になるまで、応力計算部２０、相互相関関数計算部２２、第１指標値計算部２４、自己相関関数計算部３０及び第２指標値計算部３２の処理を繰り返すように制御する。

【0058】

ここで、第１指標値と第２指標値の差に関する上記所定範囲としては、特に限定しないが、例えば、両者の差が、第１指標値の±１０％の範囲内にあるかどうかを基準とすることができる。

【0059】

第２実施形態に係る計算装置１０Ａの動作状態は、図４のフローチャートに示した通りであり、第１実施形態とは、ステップＳ１０、Ｓ１１及びＳ１２が追加された点で異なる。以下、追加したステップについて説明する。

【0060】

ステップＳ７で第１指標値計算部２４が第１指標値としての貯蔵弾性率Ｇ’（ω）とｔａｎδを算出した後、ステップＳ１０において、自己相関関数計算部３０が、ステップＳ５で得られた熱雑音を含む応力から、時間相関関数として、上記式（７）で表される応力自己相互相関関数を計算する。そしてステップＳ１１に進む。

【0061】

ステップＳ１１において、第２指標値計算部３２が、上記応力自己相互相関関数から、熱雑音を除去した応力の振幅σ_０を計算して、第２指標値としての貯蔵弾性率Ｇ’（ω）を算出する。そしてステップＳ１２に進む。

【0062】

ステップＳ１２において、評価部３４が、ステップＳ７で得られた第１指標値としての貯蔵弾性率を、ステップＳ１１で得られた第２指標値としての貯蔵弾性率と対比して、両者の差が所定範囲内であれば、ステップＳ７で得られた第１指標値の精度が確保されていると判断してステップＳ８に進み、所定範囲内でなければ、精度が確保されていないと判断してステップＳ５に進む。

【0063】

このようにしてステップＳ１２から戻ってきたステップＳ５では、分子動力学シミュレーションにおけるシミュレーション条件を変更した上で再度応力を計算し、ステップＳ６に進む。

【0064】

なお、ステップＳ８において、判定部２６は、ステップＳ７で得られた第１指標値が全ての周波数で計算されているか否かを判定する。また、ステップＳ９において、出力部２８は、上記により得られた第１指標値を高分子の粘弾性特性として出力する。

【0065】

第２実施形態によれば、歪み−応力相互相関関数に加えて、応力自己相関関数を計算し、応力自己相関関数から得られた第２指標値を用いて、歪み−応力相互相関関数から得られた第１指標値の精度を評価するので、より正確な粘弾性特性を得ることができる。

【0066】

第２実施形態において、ステップＳ４における周波数の決定は、ステップＳ１で入力された複数の周波数について低周波数側から順次に選択することが好ましい。分子動力学シミュレーションにより計算される粘弾性特性の精度は、一般に低周波数側ほどノイズが大きく、精度に劣る。そのため、低周波数側で計算精度が確保されていれば、より高周波数側でも計算精度が確保されている可能性が高い。そのため、低周波数側から計算することにより、シミュレーション条件を変更することによるそれまでの計算工数の無駄をできるだけ小さくすることができる。

【0067】

なお、図４に示す例では、ステップＳ１２において、評価部３４が、精度が確保されていないと判断したときに、ステップＳ５に進むようにしているが、ステップＳ２に進むように制御してもよい。すなわち、評価部３４は、第１指標値と第２指標値の差が所定範囲内でないとき、システム作成ステップＳ２に進み、システムを再構成した上で再度応力を計算するように制御してもよい。例えば、システムを構成する高分子モデルの分子鎖数を多くすることで精度を高めることができるので、このようにして分子鎖数を増加させながら、両者の差が所定範囲内になるまで、システム作成部１４、平衡化部１６、周波数決定部１８、応力計算部２０、相互相関関数計算部２２、第１指標値計算部２４、自己相関関数計算部３０及び第２指標値計算部３２の処理を繰り返すように制御してもよい。

【0068】

上記第１及び第２実施形態では、複数の周波数についての粘弾性特性を求めているが、他の実施形態として、単一の周波数についての粘弾性特性を求めるようにしてもよい。その場合、周波数決定部１８は、入力部１２で取得した周波数を、システムに与える歪みの周波数として決定する。

【0069】

また、上記第１及び第２実施形態では、単一の振幅についての粘弾性特性を求めているが、他の実施形態として、複数の振幅についての粘弾性特性を求めるようにしてもよい。その場合、システムに与える歪みの振幅を設定する振幅設定部を設けておき、応力計算部１８では、該振幅設定部で設定した振幅の歪みをシステムに与えながら分子動力学シミュレーションにより応力を計算する。そして、予め設定された複数の振幅のうち全ての振幅で粘弾性に関する指標値が計算されたか否かを判定し、全ての振幅で計算されるまで、指標値を取得する計算を繰り返すように制御する第２判定部を設けておけばよい。

【0070】

また、上記第１及び第２実施形態では、単一の温度についての粘弾性特性を求めているが、他の実施形態として、複数の温度についての粘弾性特性を求めるようにしてもよい。その場合、システムを平衡化する温度を設定する温度設定部を設けておき、平衡化部１６では該温度設定部で設定した温度でシステムを平衡化させ、当該温度で分子動力学シミュレーションにより粘弾性に関する指標値を計算する。そして、予め設定された複数の温度のうち全ての粘弾性に関する指標値が計算されたか否かを判定し、全ての温度で計算されるまで、指標値を取得する計算を繰り返すように制御する第３判定部を設けておけばよい。

【0071】

以上のように、本実施形態によれば、高分子材料の粘弾性特性を分子動力学シミュレーションによって予測することができる。一般にゴムに代表されるエラストマー材料では、特定の周波数の振動に対する損失を制御することが重要であり、この損失物性は粘弾性特性によって決定されるので、本発明はゴムなどの高分子の材料開発を行う場合に有益である。

【実施例】

【0072】

［実施例１］
上記実施形態により高分子の粘弾性特性が正確に計算できることを、実施例１により示す。分子動力学シミュレーションには、公開プログラムのLAMMPS［Large-Scale Atomic/Molecular Massively Parallel Simulator：米国サンディア国立研究所］を用いた。

【0073】

上記参考文献１に基づき、高分子のビーズスプリングモデルシステムとして、５０ビーズの高分子鎖４０本からなるアモルファス高分子モデルを作成した。該システムにおいて、結合ポテンシャルと非結合ポテンシャルは、参考文献１に従い、下記式の通りとした。

【0074】

【数9】

【0075】

ここで、Ｕ_FENE（ｒ）は伸びきりのある結合ポテンシャルであり、ｋとＲ_０は、それぞれ、結合の強さおよび伸びきり距離（これ以上結合した粒子間の距離が大きくならない）を示している。ｒは粒子間距離を表わしている。Ｕ_LJ（ｒ）はレナードジョーンズ非結合ポテンシャルで、Ａは相互作用の強さ、Ｂは粒子の大きさを表わすパラメータである。式（9b）ではカットオフ距離２^1/6Ｂを設定し、それよりも大きな距離にある粒子間の非結合相互作用をゼロとしている。

【0076】

次いで、該システムを、圧力４．８５で一定の下、温度１．０で平衡化した。ここで、温度や圧力などの各パラメータの単位は、レナードジョーンズ単位である。

【0077】

周波数を４．０×１０^−４［１／τ］とし、図６に示すように該周波数で振幅０．２の歪みを１０周期にわたり与えながら、分子動力学シミュレーションを行い、各時刻の応力を計算した。歪みは、シミュレーション中の座標系においてｘ軸方向に与えた。分子動力学シミュレーションにより得られた応力の時間変化は図７に示す通りである。応力の単位はレナードジョーンズパラメータのσ（長さ）とε（エネルギー）で与えられており、この計算では、両者とも１．０、即ちσ＝１．０、ε＝１．０とした。図７に示すように、分子動力学シミュレーションにより得られた応力の時間変化は熱雑音の大きいものであった。

【0078】

次に、この熱雑音を含む応力から、上記式（５）で表される歪み−応力相互相関関数を計算した。結果は図８に示す通りであり、歪み−応力相互相関（以下、歪み−応力相関ということがある。）の時間変化では、熱雑音の影響を大幅に除去することができた（即ち、式（５）においてＣ（ｔ）で表される歪みと熱雑音の相互相関関数は小さい）。なお、歪み−応力相関の次元は、歪み×応力となるので、単位は応力と同じである。

【0079】

次に、上記式（６）を用いて、応力の時間変化を最小平均二乗法により正弦曲線にフィッティングしたところ、図９に示す通りとなり、以下の結果が得られた。歪みの振幅Ｘ_０＝０．２、応力の振幅σ_０＝０．０１０８８［ε／σ^３］、周波数（角振動数）ω＝２π×４．０×１０^−４［１／τ］、位相差δ＝１．９８０４８［ｒａｄ］、残差平均の平方根／振幅＝０．０２７２４８。ここから計算した貯蔵弾性率は０．０５４２３２［ε／σ^３］であり、ｔａｎδは２．３０２７２８であった。

【0080】

また、上記熱雑音を含む応力から、上記式（７）で表される応力自己相互相関関数を計算した。結果は図１０に示す通りであり、応力自己相互相関の時間変化では、図７に示す応力自体の時間変化に比べて熱雑音の影響を除去することができた。次いで、上記式（８）を用いて、応力の時間変化を最小平均二乗法により正弦曲線にフィッティングしたところ、図１１に示す通りとなり、以下の結果が得られた。歪みの振幅Ｘ_０＝０．２、応力の振幅σ_０＝０．０１０８６［ε／σ^３］、周波数（角振動数）ω＝２π×４．０×１０^−４［１／τ］、残差平均の平方根／振幅＝１．２１７９５７。ここから計算した貯蔵弾性率は０．０５４３１８［ε／σ^３］であった。この値は、歪み−応力相関から得られた貯蔵弾性率＝０．０５４２３２とよく一致していた。

【0081】

一方、上記実施例による効果を示すために、比較例として、図７に示す応力自体の時間変化について、相関関数を計算することなく、そのままフィッティングにより粘弾性特性を計算した。すなわち、下記式（１０）を用いて、応力の時間変化を最小平均二乗法によりフィッティングしたところ、図１２に示す白抜きの曲線が得られた。結果は、歪みの振幅Ｘ_０＝０．２、応力の振幅σ_０＝０．２４０８６［ε／σ^３］、周波数（角振動数）ω＝２π×４．０×１０^−４［１／τ］、位相差δ＝２．０６２９７［ｒａｄ］、残差平均の平方根／振幅＝３０．７９８８であった。この結果は、残差がフィッティング対象の値に比べて非常に大きく、信頼性に乏しい。位相差δは、歪み−応力相互相関関数から得られたものに近いが、貯蔵弾性率は１０倍以上の違いがある。これは、応答応力に比べて非常に大きな雑音がフィッティング精度に悪影響を与えているためと考えられる。

【数10】

【0082】

［実施例２］
本実施形態により粘弾性の歪み依存性を計算できることを示すために、実施例２では歪み振幅を０．５として、その他は実施例１と同様の計算を行った。

【0083】

詳細には、圧力、温度及び周波数を実施例１と同じに設定し、図１３に示す振幅０．５の歪みを１０周期にわたり与えながら、分子動力学シミュレーションを行い、各時刻の応力を計算した。分子動力学シミュレーションにより得られた応力の時間変化は図１４に示す通りであった。

【0084】

次に、この熱雑音を含む応力から、上記式（５）で表される歪み−応力相互相関関数を計算した。結果は図１５に示す通りであり、歪み−応力相互相関の時間変化では、熱雑音の影響を大幅に除去することができた。次いで、上記式（６）を用いて、応力の時間変化を最小平均二乗法により正弦曲線にフィッティングしたところ、図１６に示す通りとなり、以下の結果が得られた。歪みの振幅Ｘ_０＝０．５、応力の振幅σ_０＝０．１３００２５［ε／σ^３］、周波数（角振動数）ω＝２π×４．０×１０^−４［１／τ］、位相差δ＝１．７３４０４５［ｒａｄ］、残差平均の平方根／振幅＝０．０１３９４９。ここから計算した貯蔵弾性率は０．２６００５［ε／σ^３］であり、ｔａｎδは６．０７１１２６であった。

【0085】

また、上記熱雑音を含む応力から、上記式（７）で表される応力自己相互相関関数を計算した。結果は図１７に示す通りであった。また、上記式（８）を用いて、応力の時間変化を最小平均二乗法により正弦曲線にフィッティングしたところ、図１８に示す通りとなり、以下の結果が得られた。歪みの振幅Ｘ_０＝０．５、応力の振幅σ_０＝０．１２３８７２［ε／σ^３］、周波数（角振動数）ω＝２π×４．０×１０^−４［１／τ］、残差平均の平方根／振幅＝０．１０１７８４。ここから計算した貯蔵弾性率は０．２４７７４５［ε／σ^３］であった。この値は、歪み−応力相関から得られた貯蔵弾性率＝０．２６００５とよく一致していた。

【0086】

一方、上記実施例による効果を示すために、比較例として、図１４に示す応力自体の時間変化について、相関関数を計算することなく、そのままフィッティングにより粘弾性特性を計算した。すなわち、上記式（１０）を用いて、応力の時間変化を最小平均二乗法によりフィッティングしたところ、図１９に示す白抜きの曲線が得られた。結果は、歪みの振幅Ｘ_０＝０．５、応力の振幅σ_０＝０．５４７９４１［ε／σ^３］、周波数（角振動数）ω＝２π×４．０×１０^−４［１／τ］、位相差δ＝１．７２４０９［ｒａｄ］、残差平均の平方根／振幅＝０．６０５１９であった。歪み０．２の場合に比べて残差は小さくなっているが、それでも歪み−応力相関に比べて数十倍の値となっており、信頼性に乏しい。位相差δは、歪み−応力相互相関関数から得られたものに近いが、貯蔵弾性率は５倍程度の違いを生じている。これは、応答応力に比べて非常に大きな雑音がフィッティング精度に悪影響を与えているためと考えられる。

【0087】

以上の実施例１及び２によれば、下記表１に示すように、システムに与える歪みの大きさにより粘弾性が変化しており、すなわち歪みに依存する粘弾性データが得られることが示された。

【0088】

【表1】

【0089】

［実施例３］
上記参考文献１に基づき作成した１００ビーズの高分子鎖５０本からなるアモルファス高分子モデルについて、温度を１．０とし、周波数を１．０×１０^−１〜３．９１×１０^−４［１／τ］として、ｘ軸方向に振幅０．２の歪みを１０周期にわたり与えながら、実施例１と同様の計算を行った。該アモルファス高分子モデルについての結合ポテンシャルと非結合ポテンシャルは、上記式（9a）及び（9b）の通りとした。

【0090】

結果は、下記表２、図２０及び図２１に示す通りであり、各周波数に対して貯蔵弾性率とｔａｎδの計算結果が得られた。

【0091】

【表2】

【0092】

上記では本発明の一実施形態を説明したが、この実施形態は、例として提示したものであり、発明の範囲を限定することは意図していない。これら新規な実施形態は、その他の様々な形態で実施されることが可能であり、発明の主旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。これら実施形態やその変形は、発明の範囲や要旨に含まれるとともに、特許請求の範囲に記載された発明とその均等の範囲に含まれる。

【符号の説明】

【0093】

１０…計算装置 １２…入力部 １４…システム作成部
１６…平衡化部 １８…周波数決定部 ２０…応力計算部
２２…相互相関関数計算部 ２４…第１指標値計算部 ２６…判定部
３０…自己相関関数計算部 ３２…第２指標値計算部 ３４…評価部

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】

【図13】

【図14】

【図15】

【図16】

【図17】

【図18】

【図19】

【図20】

【図21】