JP2009271719A - 高分子材料のシミュレーション方法 - Google Patents

Abstract

【課題】 能率良く高分子材料のエネルギーロス発生箇所を特定する。
【解決手段】 分子動力学計算を用いて高分子材料のシミュレーションを行う方法であって、解析対象の高分子材料から、その分子構造に基づいて、原子又はその集合体を表す有限個の粒子と、この粒子間を結合して各粒子の相対位置を定義する結合鎖とを含む三次元構造の分子鎖モデルを設定するモデル設定ステップＳ１、前記分子鎖モデルを分子動力学計算に基づいて安定化させる構造緩和のシミュレーションを行うステップＳ２、前記構造緩和後の分子鎖モデルを分子鎖部分に分解しかつそれぞれの分子鎖部分に生じる応力を計算する応力計算ステップＳ３、及び前記応力計算ステップで計算された各分子鎖部分の応力に基づいてエネルギーロスの主要な発生箇所を決定する決定ステップＳ４を含むことを特徴とする高分子材料のシミュレーション方法。
【選択図】 図１

Description

本発明は、能率良く高分子材料のエネルギーロス発生箇所を分子レベルで特定することができるシミュレーション方法に関する。

プラスチック又はゴムのような高分子材料において、そのエネルギーロスは製品の様々な特性に影響を及ぼす重要な物理量である。例えば、ゴム製品であるタイヤにおいて、エネルギーロスは燃費性能やグリップ性能に密接に係わっており、適切な制御が必要と考えられる。従来、目的とするゴム材料を得るために、ゴム材料の配合設計、試作及びエネルギーロスの測定試験といった工程を含む開発サイクルが繰り返される。

しかしながら、従来の開発サイクルでは、ゴムのどの部分、具体的にはどの分子でエネルギーロスが多く発生しているか等の情報を知ることはできない。従って、例えば、エネルギー伝達効率の良いゴム材料を開発するためには、これまでの知見や経験などに依存した試行錯誤的な作業が多く効率的ではない。

このような実情に鑑み、本件出願人は、分子動力学計算を用いて、高分子材料のエネルギーロス発生箇所を分子レベルで決定しうるコンピュータシミュレーション方法を既に提案した（下記特許文献１参照）。この方法では、下記のイないしニのステップを経て行われる。

（イ）高分子材料の分子構造に基づいた粒子・結合鎖を具える分子鎖モデルを含む分子構造モデルの設定
（ロ）人為的影響を取り除くための分子構造モデルの構造緩和処理
（ハ）構造緩和後、分子構造モデルに任意の軸方向についての負荷変形及び除荷変形を与える変形シミュレーションの実行
（ニ）変形シミュレーションの後、特定の歪状態における負荷変形時における前記軸方向に対する結合鎖の角度αi と同一の歪状態かつ除荷変形時の結合鎖の角度α'iとの差（α'i−αi ）が大きい箇所をエネルギーロスの発生箇所として特定

しかしながら、上述のような変形シミュレーション（ステップ（ハ））は、原子や分子の動きを追跡しなければならないため、計算の時間増分がフェムト秒オーダのように非常に小さくならざるを得ない。この結果、変形シミュレーションの計算には、非常に多くの時間が必要であった。また、変形シミュレーションに要する計算時間を短縮させるために、変形速度を高速にする必要があるなど、種々の制約が存在する。このような制約は次のような不具合をもたらす。即ち、高分子材料は変形速度によって特性や変形挙動が大きく異なる。例えば、タイヤにおいては、そのゴム部分や走行速度によって変形速度が種々異なる。従って、これらを評価するためには、種々の変形速度でのシミュレーションを行う必要がある。しかしながら、これまでの分子動力学による変形シミュレーションでは、非常に高速な変形速度の特性、変形挙動しか得られず十分ではない。

特開２００７−１０７９６８号公報

発明者らは、変形シミュレーションに非常に多くの時間を要するという事実に鑑み、種々の研究を重ねた。その結果、詳細な理由は後述するが、高分子材料の分子鎖モデルの中で主要なエネルギーロスを発生する箇所は、意外にも分子構造モデルを構造緩和した後に特定の圧縮応力状態の分子鎖部分であるという新規な知見を得た。従って、構造緩和後の分子鎖モデルの中で、このような特定の応力状態の分子鎖部分を知ることができれば、従来のような変形シミュレーションを行うことなしに分子鎖モデルの中のエネルギーロスの主たる発生箇所を決定でき、大幅に開発能率を高めることができる。

以上のように、本発明は、高分子材料の主要なエネルギーロスの発生箇所を、変形シミュレーションを行うことなく分子レベルで突き止めることにより、計算時間の大幅な短縮を可能とし、ひいては高分子材料の設計ないし開発の能率化を図るのに役立つシミュレーション方法を提供することを主たる目的としている。

本発明のうち請求項１記載の発明は、分子動力学計算を用いて高分子材料のシミュレーションを行う方法であって、解析対象の高分子材料の分子構造に基づき、任意の空間を有するセル内に、原子又はその集合体を表す有限個の粒子と、この粒子間を結合して各粒子の相対位置を定義する結合鎖とからなる三次元構造の分子鎖モデルを配置した分子構造モデルを設定するステップ、前記分子構造モデルを分子動力学計算に基づいて安定化させる構造緩和のシミュレーションを行うステップ、前記構造緩和後の分子鎖モデルを分子鎖部分に分解しかつそれぞれの分子鎖部分に生じる応力を計算する応力計算ステップ、及び前記応力計算ステップで計算された各分子鎖部分の応力に基づいてエネルギーロスの主要な発生箇所を決定する決定ステップを含むことを特徴とする。

また請求項２記載の発明は、前記決定ステップは、圧縮応力が大きい分子鎖部分をエネルギーロスの主要な発生箇所として決定する請求項１記載の高分子材料のシミュレーション方法である。

また請求項３記載の発明は、前記応力計算ステップは、結合鎖を介して連続する４粒子が構成する立体配座に基づいて前記分子鎖モデルを分子鎖部分に分解する請求項１又は２記載の高分子材料のシミュレーション方法である。

また請求項４記載の発明は、前記応力計算ステップは、前記分子鎖モデルを２粒子の結合長に基づいて分解しかつそれぞれの結合長を中心軸とする立体配座毎に区別して応力を計算する第１の計算ステップと、この第１の計算ステップの結果に基づいて、大きな圧縮応力が生じる結合長−立体配座の組み合わせを特定するステップと、前記特定された結合長−立体配座の組み合わせについて、その両端に結合される粒子をさらに考慮して分子鎖モデルを３つの連続する立体配座の並び別にグループ化しかつ各グループ毎の応力を計算する第２の計算ステップとを含むとともに、前記決定ステップは、前記第２の計算ステップの結果に基づいて、圧縮応力が大きいグループを、エネルギーロスの主要な発生箇所として特定する請求項１記載の高分子材料のシミュレーション方法である。

また請求項５記載の発明は、前記エネルギーロスの主要な発生箇所として決定された部分を他の部分と差別化して表示するステップをさらに含む請求項１乃至４のいずれかに記載の高分子材料のシミュレーション方法である。

本発明のシミュレーション方法では、構造緩和後の分子鎖モデルを分子鎖部分に分解しかつそれぞれの分子鎖部分に生じる応力を計算し、その各分子鎖部分の応力に基づいてエネルギーロスの主要な発生箇所を決定することができる。従って、高分子材料の分子鎖モデルの具体的な変形シミュレーションを行うことなく、エネルギーロスの発生箇所を原子又は分子レベルで特定できる。よって、高分子材料の開発を短時間でかつ能率的に行うことができる。

以下、本発明の実施の一形態が図面に基づき説明される。
本発明では、分子動力学計算を用いた高分子材料のシミュレーションがコンピュータ装置（図示せず）によって行われる。

前記分子動力学（ Molecular Dynamics : ＭＤ）計算は、コンピュータ装置を用いたシミュレーションの一種である。分子動力学では、先ず、解析の対象となる高分子材料の分子構造に基づいて多数の原子又はその集合体（分子を含む）からなる粒子を配置したモデル（系）が設定される。そして、このモデルは、配置された全ての粒子が古典力学に従うものとされ、ニュートンの運動方程式を適用して各時刻における全ての粒子の動きが追跡される。このような分子動力学計算によれば、粒子の微視的な運動をできる限り正確に追跡することができ、実験結果などに頼ることなく、物質の性質や運動を明らかにすることができる。また、構造緩和や追跡時間等の調節により、粒子の初期配置といった人為的操作に依存しない正確なシミュレーション結果を得ることができる。

前記高分子材料は、高分子を含む材料であれば特に限定されないが、例えば天然ゴム、合成ゴム及び合成樹脂を少なくとも含む有機高分子材料が好適である。本実施形態では、高分子材料としてcis-１，４ポリブタジエンを用いて以後の説明が行われる。

また、前記コンピュータ装置は、特に制限されないが、好ましくは計算処理能力の高いスーパーコンピュータ等が好適に用いられる。

図１には、本実施形態のシミュレーション方法の具体的な処理手順が示される。
本実施形態のシミュレーション方法では、先ず解析対象となる高分子材料であるcis-１，４ポリブタジエンから、その分子構造に基づいて、シミュレーション用の分子構造モデル１（図９に示す）が設定される（ステップＳ１）。

図２（ａ）には、cis-１，４ポリブタジエンのモノマー｛−［ＣＨ_２−ＣＨ＝ＣＨ−ＣＨ_２］−｝の構造図を示す。cis-１，４ポリブタジエンは、メチレン基（−ＣＨ_２−）とメチン基（−ＣＨ−）とから構成される。

また、図２（ｂ）には、上記モノマーをシミュレーション上で取り扱うためにその分子構造に基づいてモデル化（疑似化）したモノマーモデル２ａを、さらに図２（ｃ）には、モノマーモデル２ａを有限個連続して連結することにより形成された１本の分子鎖モデル２の部分平面図がそれぞれ示される。

前記分子鎖モデル２は、原子又はその集合体を表す複数かつ有限個の粒子３と、この粒子３、３間を結合して各粒子の相対位置を定義する結合鎖４とを含む三次元構造を有する。また、粒子３は、メチレン基に対応する粒子３ａと、メチン基に対応する粒子３ｂとからなる（なお、これらの粒子を総称するとき、単に「粒子３」と表示する。）。

本実施形態の分子鎖モデル２は、"United Atom Model"に従ってモデル化されている。即ち、ポリブタジエンの水素原子は陽に扱わず、ＣＨ、ＣＨ_２は、それぞれ水素原子の効果を織り込んだ原子の集合体としての一つの粒子３として取り扱われる。ただし、全ての原子をそれぞれ粒子３として捉えるいわゆる "Full Atom model"により、解析対象の分子構造がモデル化されても良い。このように、分子鎖モデル２の粒子３は、cis-１，４ポリブタジエンの分子構造を表すことができる。

また、前記粒子３は、分子動力学計算に基づいたシミュレーションにおける運動方程式の質点として取り扱われる。従って、粒子３には、質量はもとより、体積、直径、電荷及び／又は初期座標などのパラメータが定義される。これらの各パラメータは、数値情報としてコンピュータ装置に入力される。また、前記結合鎖４は、粒子３、３間の相対位置を特定しかつ拘束するもので、これは例えばベクトル情報としてコンピュータ装置に入力される。

また、分子鎖モデル２は、３次元構造を有し、図３（ａ）〜（ｃ）に示されるように、各粒子３、３間の結合長さ（「ボンドストレッチ」と言う場合がある。）である結合長ｒ、同図（ｂ）に示されるように、結合鎖４を介して連続する３つの粒子３がなす角度（単に「ベンディング」と言う場合がある。）である結合角θ、及び同図（ｃ）に示されるように、結合鎖４を介して連続する４つの粒子３において、隣り合う３つの粒子が作る第１の平面Ｐ１と第２の平面Ｐ２のなす角度（以下、単に「トーション」と言う場合がある）である二面角φがそれぞれ定義される。なお、分子鎖モデル２は、外力又は内力によって、上記結合長、結合角及び二面角が変化する。

さらに、分子動力学計算では、分子鎖モデル２に対して、図３（ａ）の結合長ｒのポテンシャル（略記号：Ｅbs）、同図（ｂ）の結合角θ（略記号：Ｅbe）、同図（ｃ）の二面角φのポテンシャル（略記号Ｅto）及び図３（ｄ）で示されるように、互いに連結されていない粒子間のファンデルワールス力のポテンシャルｒ’（略記号：Ｅvw）の相互作用が定義される。これらのポテンシャル関数は、種々のものが採用できるが、本実施形態では、下記式（１）〜（５）が採用される。

なお、ポリブタジエンのモノマーは−［ＣＨ_２−ＣＨ＝ＣＨ−ＣＨ_２］一であるから、結合長は、ＣＨ＝ＣＨ、ＣＨ_２−ＣＨ_２及び両者を繋ぐＣＨ−ＣＨ_２の３種類になる。いずれも調和振動子として長さの変化に対するポテンシャルが表現される。また、上記式（２）中のばね定数ｋ_ｒならびに平衡長ｒ_０は、本実施形態では、各結合長毎にそれぞれ異なる値が以下のように定められる。
ＣＨ_２−ＣＨ_２：
ｋ_ｒ＝６６３（kJ／mol・nm^２）、ｒ_０＝１．５４（nm）
ＣＨ−ＣＨ_２：
ｋ_ｒ＝７６９（kJ／mol・nm^２）、ｒ_０＝１．５０（nm）
ＣＨ＝ＣＨ：
ｋ_ｒ＝１０３３（kJ／mol・nm^２）、ｒ_０＝１．３３５（nm）

同様に、ポリブタジエンでは、３つの粒子間の結合角は、ＣＨ_２−ＣＨ_２−ＣＨと、ＣＨ＝ＣＨ−ＣＨ_２の２種類になる。いずれも調和振動子として角度の変化に対するポテンシャルが表現される。上記式（３）中のばね定数ｋθ及び平衡角度θ_０は、本実施形態では各結合角毎にそれぞれ異なる値が以下のように設定される。
ＣＨ_２−ＣＨ_２−ＣＨ：
ｋ_θ＝４８２（kJ／mol・rad^２）、θ_０＝１１１．６（deg）
ＣＨ＝ＣＨ−ＣＨ2
ｋ_θ＝３７４（kJ／mol・rad^２）、θ_０＝１２４．０（deg）

さらに、４つの粒子で作られる二面角φは、図３（ｃ）に示したように、一つの結合鎖４を中心軸ＣＬとする分子鎖内の回転を表すポテンシャルであるが、その中心軸ＣＬで区別すると、ＣＨ＝ＣＨ、ＣＨ_２−ＣＨ_２、及びＣＨ−ＣＨ_２の３種類が存在する。これらは、下記のように、上記式（４）中の安定点ならびにエネルギー障壁が異なっている。
ＣＨ_２−ＣＨ_２：
ｋ₁＝３．３５（ｋJ／mol）
ｋ₃＝１３．４（ｋJ／mol）
ｋ₂＝ｋ₄＝０（ｋJ／mol）
ＣＨ−ＣＨ_２：
ｋ₁＝ｋ₂＝ｋ₃＝０（ｋJ／mol）
ｋ₄＝７．９５（ｋJ／mol）
ＣＨ＝ＣＨ：
ｋ₁＝ｋ₃＝ｋ₄＝０（ｋJ／mol）
ｋ₂＝５８．６（ｋJ／mol）

また、図４には、二面角φによる１粒子当たりのトーションポテンシャル（φto）の変化を示す。また、図５には、中心軸ＣＬをなす結合鎖４に沿って見たときの両端の２粒子の位置関係（Newman投影図）を示す。

図４、図５及び図６に示されるように、中心軸ＣＬをなす結合鎖４が二重結合のＣＨ＝ＣＨの場合、二面角φは０゜と１８０゜が安定点になる。そして、０゜及び１８０゜の二面角を持つ４つの粒子がなす単結合の立体構造、即ち立体配座（又はコンフォメーションとも呼ばれる）はそれぞれcis及びtransと呼ばれる。

また、図４、図５及び図７に示されるように、中心軸ＣＬがＣＨ−ＣＨ_２の場合、二面角は０゜及び±１２０゜が安定点になる。前者の立体配座はcis、後者の立体配座はanticlinalと呼ばれる。

さらに、図４、図５及び図８に示されるように、中心軸ＣＬがＣＨ_２−ＣＨ_２の場合、二面角φは、１８０゜の安定点（trans）に加え、transに比べわずかに高いエネルギーを極小値とする準安定な二面角として±６７．５゜が存在する。この立体配座は、gaucheと呼ばれる。

また、図９に示されるに、本実施形態の分子構造モデル１は、予め定められた微小体積を持ったセルＳの中に複数本の分子鎖モデル２を配置することで設定される。該セルＳは、分子鎖モデル２の可動範囲を定める境界をなす（換言すれば、セルＳは、解析対象のポリブタジエンの微小部分に相当する。）。本実施形態のセルＳは、２０nm×２０nm×２０nmの微小な立方体として定義されるが、これに限定されるものではなく、種々の形状及び／又は大きさで定めることができる。

分子構造モデル１を設定する手順は、特に限定されるものではない。例えば、セルＳ内の任意の位置に開始点となる初期の立体配座（例えばcisに固定されたモノマー：ＣＨ_２−ＣＨ＝ＣＨ−ＣＨ_２）を定義し、その一端又は両端に、乱数を用いて二面角φをランダムに設定することにより、ランダムコイル状の分子鎖モデル２をセル中Ｓに成長させることで設定できる。分子鎖モデル２を成長させる途中において、セルＳの壁に遮られ目的とする長さの分子鎖モデル設定できなくなった場合、当該分子鎖モデルを破棄し、開始点を新たに設定して新たに分子鎖モデル２を設定することができる。

これらの処理を繰り返すことにより、セルＳ内に必要な本数の分子鎖モデル２を配置した分子構造モデル１を作ることができる。即ち、解析対象のポリブタジエンの密度になるように、セルＳ内の分子鎖モデル２の本数ないし粒子数が調整される。本実施形態において、セルＳと分子鎖モデル２とからなるポリブタジエンの構造モデルは、次のように設定された。
密度：０．８６ｇ／cm^３
分子鎖モデルの総数：１０００本
総粒子数：３００７９２個
また、本実施形態では、前記セルＳの中に、長さ（ここで言う長さは、１本の分子鎖モデル２に含まれる粒子３の数とし、以下同様である。）が異なる複数種類の分子鎖モデル２が定義されている。具体的には、各分子鎖モデル２は、その平均長さが３００、かつ、標準偏差が５０の正規分布に従って定められた。

次に、本実施形態では、分子動力学計算に基づいて該分子鎖モデル２を安定化させる構造緩和のシミュレーションが行われる（ステップＳ２）。

この構造緩和のシミュレーションでは、初期設定された分子構造モデル１について、全方向自由境界条件の下、予め定めた時間で分子動力学計算が行われる。これにより、各分子鎖モデル２は、セルＳ内を古典力学に従って自由に移動・変形する。これは、分子鎖モデル２の初期構造の影響（人為的影響）を無くし、分子構造モデル１の安定化ないし準安定化を図ることができる。

構造緩和のシミュレーションにおいて、積分の時間ステップ、温度及び／又は緩和時間等は任意に定めることができる。本実施形態の構造緩和のシミュレーションは、次のような条件で行われた。なお、緩和時間には、例えば、ポテンシャルエネルギー及び応力が、実質的に変化しなくなるのに必要かつ十分な時間が設定されることが望ましい。
積分：Verlet法
積分の時間ステップ：０．１fs
温度：３００Ｋ（速度スケーリング法によって制御）
緩和時間：４０ps（＝４００００fs）

図１０には、構造緩和中の分子鎖モデル２のｚ軸方向の応力及びポテンシャルエネルギーの時系列変化を示す。ここで、「分子鎖モデル２のｚ軸方向の応力」とは、セルＳ内の全ての分子鎖モデル２のｚ軸方向の応力を示す。具体的には、図１１に模式的に示されるように、セルＳ内の各粒子３において、ｚ軸方向の正方向に作用する力（切断法による内力）Ｐ_z ⁱ、Ｐ_z ^j、Ｐ_z ^k…を加算することにより該セルＳに働くｚ方向の軸力を求め、それを可動部分の体積で徐すことにより下式で計算されるものとする。
Ｐ＝ΣＰ_z ⁱ／（Ｌ_ｘ×Ｌ_z×Ｌ_ｙ）

図１０から明らかなように、分子鎖モデル２には、ｔ＝０近傍の初期構造時には、極めて大きな圧縮応力が生じている。しかし、構造緩和開始から２０ps以降、応力は、振動してはいるもののその平均が約４０ＭＰａ程度でほとんど変化していないことが分かる。また、ポテンシャルエネルギーも、初期構造時の極めて高いエネルギー状態から急激に低下した後、緩やかに減少して４０psでは実質的に０となって変化が止まっている。つまり、本実施形態の分子鎖モデル２は、４０psの時間で十分な構造緩和が行われたことが分かる。ただし、ポテンシャルエネルギーは、構造緩和後、必ずしも零になるとは限らない。

従来のシミュレーション方法では、上記構造緩和を行った後、変形シミュレーションが行われる。しかし、この変形シミュレーションは、分子構造モデル１に変形を与える必要があるため、非常に多くの計算時間を必要とすることは上述の通りである。

発明者らは、構造緩和後の分子鎖モデル２を各分子鎖部分の応力に分けて分析し、これを変形シミュレーションの結果と突き合わせることを行った。その結果、高分子材料の分子構造の中で主要なエネルギーロスを発生する箇所は、意外にも分子鎖モデル２を構造緩和した後に特定の応力状態を示すことを知見した。具体的には、主要なエネルギーロスを発生する箇所は、構造緩和後に圧縮応力が残存している部分であることを知見した。従って、構造緩和後の分子鎖モデル２の中で、このような圧縮応力が残存する部分を特定できれば、従来のような変形シミュレーションを行うことなしに分子鎖モデル２の中のエネルギーロスの主たる発生箇所を能率良く決定でき、大幅に開発能率を高めることができる。以下、さらに詳細に述べる。

先ず、粒子３に作用する力は、前記式（１）の空間微分から得られる。従って、式（２）〜（５）から、各ポテンシャル成分、即ち結合長（ボンドストレッチ）、結合角（ベンディング）、二面角（トーション）及びファンデルワールスのポテンシャルへの寄与分を計算することができる。図１３には、このようにして得られた上述の応力σ_ｚｚの変化に対する結合長の寄与分を示す。図１３を図１０の応力変化と比較すれば、構造緩和後の分子鎖モデル２の残留応力σ_ｚｚは、実質的に結合長のポテンシャルによるものであることが分かる。なお、ベンディング、トーション及びファンデルワールスについて計算したが、いずれもこれらの寄与分はほぼ零であった。

次に、発明者らは、構造緩和後の残留応力を担っている結合長をさらに種別化し、その応力状態を調べた。具体的には、ポリブタジエンに含まれる３種類の結合長ＣＨ＝ＣＨ、ＣＨ_２−ＣＨ_２及びＣＨ−ＣＨ_２をそれぞれが中心軸ＣＬとなる立体配座によって分類し、それらの応力の分布が調べられた。例えば、図６に示したように、結合長ＣＨ＝ＣＨを中心軸ＣＬとする場合、その結合長の応力が、安定的な立体配座であるcis及びtransに分類して計算された。また、図７に示したように、結合長ＣＨ−ＣＨ_２を中心軸ＣＬとして有する場合、その安定的な立体配座であるcis及びanticlinalについて、個々の結合長の応力が計算された。さらに、図８に示したように、結合長ＣＨ_２−ＣＨ_２を中心軸ＣＬとして有する場合、その安定的な立体配座であるgauche及びtransについて、個々の結合長の応力が計算された。

図１４（ａ）〜（ｃ）には、上記で計算された各応力を構造緩和中の時系列変化で示す。なお、これらの各応力の和が図１３の系全体のボンドストレッチの応力に相当するものである。また、構造緩和中での立体配座の決定は、各安定点の角度±３０゜以内にある二面角をcisやtransといった安定的な立体配座に含めて評価し、それ以外は"other"として分類した。また、図１４では、ＣＨ＝ＣＨ、ＣＨ−ＣＨ_２、ＣＨ_２−ＣＨ_２の結合長の応力が立体配座毎に計算されているが、これらの応力は、結合長の両端の粒子のみを対象として図１１の場合と同様に計算される。

図１４から明らかなように、構造緩和後、大きい引張応力（正の応力）は、中心軸ＣＬがＣＨ＝ＣＨの結合長ではtrans及びother、ＣＨ−ＣＨ_２の結合長ではcis及びotherである。これに対して、構造緩和後に大きい圧縮応力（負の応力）を生じるのは、中心軸ＣＬがＣＨ＝ＣＨの結合長ではcis、及び、同ＣＨ−ＣＨ_２の結合長ではanticlinalである。なお図１４（ｃ）のＣＨ_２−ＣＨ_２の結合長では、同図（ａ）、（ｂ）に比べて、発生している引張及び圧縮の応力は無視しうる程度に十分小さいと言える。

発明者らは、構造緩和後に大きい圧縮応力（負の応力）が生じている結合長−立体配座の組み合わせ（ＣＨ＝ＣＨのcis及びＣＨ−ＣＨ_２のanticlinal）に着目し、これらをさらに種別化して応力状態を調べることを試みた。具体的には、前記結合長−立体配座の組み合わせ（これは４つの粒子３から構成される）に、その両端に結合される２つの粒子３をさらに考慮した合計６つの粒子群を特定し、この粒子群を３つの連続する立体配座の並び別にグループ化しかつ各グループ毎の応力を計算した。なお、構造緩和中での立体配座の決定は、上述の通り、各安定点±３０゜以内にある二面角をcisやtransといった安定的な立体配座として評価し、それ以外は"other"として分類した。その結果を図１５（ａ）、（ｂ）に示す。

図１５（ａ）に示されるように、ＣＨ＝ＣＨの結合長を中心に有する６つの粒子群では、（anticlinal）−（cis）−（anticlinal）の３つの立体配座の並びのグループで大きな圧縮応力を生じていることが分かる。このときの粒子３の並びを図１６（ａ）に示す。

同様に、図１５（ｂ）に示されるように、ＣＨ−ＣＨ_２の結合長を中心に有する６つの粒子群では、cis−anticlinal−gauche及びcis−anticlinal−transの３つの立体配座の並びのグループがそれぞれ圧縮応力を示すことが分かった。このときの粒子群の並びを図１６（ｂ）及び（ｃ）にそれぞれ示す。

図１６（ｂ）及び（ｃ）を精査すると、それらの中心の４粒子の構造（ＣＨ＝ＣＨ−ＣＨ_２−ＣＨ_２）は、図１６（ａ）の粒子群の右側４つの並びと共通していることが分かる。また、図１５（ａ）、（ｂ）の比較においても、図１６（ａ）の粒子の並びのグループの方が大きい圧縮応力を生じる。これらより、図１６（ａ）の並びのグループ（粒子群）に、圧縮応力が生じていることが分かる。そこで、発明者らは、このようなcis１，４のモノマー構造の両端の二面角がanticlinalをなす粒子群（つまり、図１６（ａ）の実線で表された６つの粒子群）を圧縮ノードとして特定した。

さらに、発明者らは、上述の構造緩和後の分子構造モデル１を用いて変形シミュレーションを行い、それから得られる応力−ひずみ曲線のヒステリシスロスにこのような圧縮ノードがどのように影響しているのかを調べた。

変形シミュレーションは、図９で示した前記セルＳの上下面から１nm以内の領域に存在する粒子３をつかみ部とし、それ以外の粒子３は自由に運動できる（但し、境界を越えることはできない）条件として、セルＳの上下方向（ｚ軸方向）に最大ひずみε_ｚｚ＝１．０となるまで引っ張ることにより行われた。この際、積分の時間ステップは１．０fsとし、毎ステップΔε_ｚｚ＝１．０×１０^-5のひずみ増分を与えてひずみを増加させた（ひずみ速度に換算すると１．０×１０¹⁰／ｓ）。なおセルｓの上側のつかみ部から下側のつかみ部まで連続する分子鎖モデル２は存在していない。

また、ひずみはその全粒子３のｚ座標をスケーリング（アフィン変形）することで与えたが、ひずみ増分が小さく系の内部構造緩和が逐次生じるため、変形後期の内部構造は、初期のそれをアフィン変形させたものとは異なるものとなる。そして、最大ひずみε_ｚｚ＝１．０に達した後、保持時間０でひずみを反転させ、同一ひずみ速度で除荷が行われた。この引張一除荷を１サイクルとして、２サイクルの繰り返し変形シミュレーションが行われた。変形シミュレーションにおいて、系の温度は３００Ｋとし、速度スケーリングにより制御した。

図１７には、変形シミュレーションによるｚ軸方向の応力一ひずみの関係を示す。濃い色の線が１サイクル目、薄い色の線が２サイクル目の変形を示す。図では、負荷及び除荷サイクルでの応力経路の相違により、複数本の分子鎖モデル２からなるポリブタジエンの分子構造モデル１には大きなエネルギーロスが再現されている。

発明者らは、図１７で得られた応力−ひずみの関係について、各ポテンシャル毎の寄与分を計算した。その結果を図１８に示す。図１８から明らかなように、変形シミュレーションで得られたヒステリシスロス（エネルギーロス）も、大部分が結合長によって担われており、残りはファンデルワールスであることが分かった。結合角及び二面角は、系の応力変化にあまり寄与しておらず、ヒステリシスロスにも実質的に影響していないことも分かった。

そこで、発明者らは、図１７の１サイクル目の応力変化を、圧縮ノードに属する結合長に生じる応力、及びそれ以外（非圧縮ノード）の結合長に生じる応力にそれぞれ分けて計算を行った。その結果を図１９に示す。図１９のグラフの負側には、「圧縮ノードに属する結合長に生じる応力」が記載されている。つまり、cis１，４のモノマー構造の両端の二面角がanticlinalをなす粒子群の応力が記載されている。また、図１９の正側には、上記圧縮ノードに属さない立体配座をとる粒子、即ち３種類の結合長ＣＨ_２−ＣＨ_２、ＣＨ−ＣＨ_２及びＣＨ＝ＣＨの応力寄与分がそれぞれ記載される。

図１９から明らかなように、非圧縮ノードに属するＣＨ＝ＣＨ及びＣＨ_２−ＣＨ_２の結合長では、初期引張応力のままほとんど変化せず実質的に一定の応力を示し、除荷変形過では、僅かに低い応力側にずれているだけである。つまり、ヒステリシスヘの寄与は殆ど無いと言える。同様に、非圧縮ノードに属するＣＨ−ＣＨ_２の結合長は、引張りを与えると初期引張応力から５ＭＰａ程度の上昇を示すが、すぐにほぼ一定となり変化していないことが分かる。なお、この結合長も、先のＣＨ−ＣＨ、ＣＨ_２−ＣＨ_２と同様、ヒステリシスロスへの寄与は小さい。

他方、圧縮ノードは、大きなヒステリシスロスを示していることが分かる。初期状態で残存していた圧縮ノードの圧縮応力は、引張に対して”バッファ”となり、引張過程ではひずみの増加につれて緩やかに応力が回復し上に凸の曲線となる。一方、除荷初期の応答は、やや下に凸の曲線となる。また、ひずみε_ｚｚが０．５を下回ると、初期応力よりも圧縮応力が大きく変形抵抗が急激に上昇する。このため、圧縮ノードの応力曲線は、上に凸（圧縮ひずみ、応力を正とすると下に凸の２次曲線状）になる。以上より、変形シミュレーションによって得られたヒステリシスロスの大部分は、圧縮ノードによって生じることが判明した。

また、図２０（ａ）には、圧縮ノードの割合が大きい分子鎖モデル２（高圧縮モデル２Ｘ）の一例を、また図２０（ｂ）には、圧縮ノードの割合が小さい分子鎖モデル２（低圧縮モデル２Ｙ）の一例をそれぞれ視覚化して示す。圧縮ノードを薄い灰色で、それ以外の構造をとる分子鎖部分は黒色で着色している。なお、圧縮ノードの割合は、分子鎖モデル２内で圧縮ノードに分類された粒子数を、分子鎖モデル２の全粒子数で除すことにより計算される。

図２０から明らかなように、高圧縮モデル２Ｘは、粒子が非常に密に分布し、その全体形状は糸まり状を呈している。このため、圧縮ノードの結合長に圧縮応力が生じたと考えられる。他方、低圧縮モデル２Ｙは、粒子３が比較的疎に分布しかつ直線状の部分が多く、高圧縮モデル２Ｘに比べるとやや広がった形状を呈していることがわかる。

以上を整理すると、次の通りである。
ａ）ランダムに成長させたポリブタジエンの分子鎖モデル２は、十分な構造緩和を行っても一部の結合長に応力が残留する。
ｂ）残留応力を生じる結合長のうち、cis−１、４の立体配座のモノマーにおいて、両端のＣＨ−ＣＨ_２の二面角がanticlinalとなるノードの結合長は、大きな圧縮応力が生じる圧縮ノードとなっている。また、他の結合長にはすべて引張応力が生じる。
ｃ）圧縮ノード以外の結合長は、引張変形時にほとんど応力上昇せず、ヒステリシスヘの寄与は小さい。
ｄ）圧縮ノードに属する結合長は、引張によって圧縮応力が大きく回復するとともに、除荷時には圧縮側の経路をとって大きなヒステリシスを示す。
以上のように、構造緩和後の分子鎖モデル２に残留圧縮応力が生じている分子鎖部分（分子鎖モデルの局所構造）を突き止めることにより、ヒステリシスロスの主要な発生箇所を特定することができる。

従って、本発明では、このような新規な知見に基づき、図１に示されるように、構造緩和後の分子鎖モデル２を分子鎖部分に分解しかつそれぞれの分子鎖部分に生じる応力を計算し（応力計算ステップＳ３）、該応力計算ステップＳ３の結果から圧縮応力の大きい分子鎖部分をエネルギーロスの主要な発生箇所として決定する（決定ステップＳ４）。

また、応力計算ステップＳ３では、上述のように、分子鎖モデル２を段階的に分子鎖部分に分解することが望ましい。即ち、図１２にまとめたように、先ず、各分子鎖モデル２を結合長毎に分解しそれぞれの結合長を中心軸とする立体配座毎の応力を計算し（第１の計算ステップＳ３１、図１４）、この第１の計算ステップの結果に基づいて、大きな圧縮応力が生じる結合長−立体配座の組み合わせを特定し（ステップＳ３２）、さらにステップＳ３２で特定された結合長−立体配座の組み合わせについて、その両端に結合される粒子をさらに考慮して分子鎖モデルを３つの連続する立体配座の並び別にグループ化しかつ各グループ毎の応力を計算することが望ましい（第２の計算ステップＳ３３、図１５）。そして、決定ステップＳ４では、前記第２の計算ステップＳ３３の結果に基づいて、圧縮応力が大きいグループを、エネルギーロスの主要な発生箇所として決めれば良い。これにより、より精度良くエネルギーロスの主要な発生箇所を特定できる。

なお、前記「特定」は、例えばコンピュータ装置上に記憶されている該当する粒子３にフラグ等の何らかの情報を書き込むことによって行うことができる。

次に、本実施形態では、エネルギーロスの発生箇所として特定された結合鎖４を、他の結合鎖４と差別化して前記分子鎖モデルを視覚化して表示するステップをさらに含む（ステップＳ５）。このステップでは、例えば図２０に示したように、圧縮ノードと非圧縮ノードとを色を変えてモニタ等に表示又は印刷することにより行われる。これにより、分子構造中のどの部分でエネルギーロスが多く生じているのかを、容易に確認かつ全体的に判断できる。そして、そのようにエネルギーロス発生箇所として特定された部分に種々の改良を加えることで高分子材料の設計、試作を能率的に行うことができる。

以上の本発明の実施形態についてポリブタジエンを例に挙げて説明したが、他の高分子材料についても同様の手法で実施しうるのは言うまでもない。

本発明の処理手順を示すフローチャートである。 （ａ）はcis１，４ポリブタジエンのモノマーの構造式、（ｂ）はそれをモデル化した図、（ｃ）は分子鎖モデルの部分平面図である。 （ａ）〜（ｄ）はポテンシャルを説明する分子鎖モデルの部分図である。 トーションのポテンシャル関数である。 中心軸をなす結合鎖に沿って見たときの両端の２粒子の位置関係を示すNewman投影図である。 結合鎖を介して連続する４つの粒子の立体配座を説明する図である。 結合鎖を介して連続する４つの粒子の立体配座を説明する図である。 結合鎖を介して連続する４つの粒子の立体配座を説明する図である。 セル内に分子鎖モデルを配置した高分子材料の構造モデルの斜視図である。 構造緩和中の分子鎖モデルのｚ軸方向の応力及びポテンシャルエネルギーの時系列変化を示すグラフである。 セルの模式図である。 応力計算ステップの具体例を示すフローチャートである。 分子鎖部応力の変化に対する結合長の寄与分を示す。 （ａ）〜（ｃ）には、分子鎖モデルの各結合長について、応力を立体配座毎に分解して示すグラフである。 （ａ）及び（ｂ）は、各結合長について、３つの連続する立体配座の並び別にグループ化しかつ各グループ毎の応力を計算したグラフである。 （ａ）〜（ｃ）は分子鎖モデルの部分平面図である。 変形シミュレーションによるｚ軸方向の応力一ひずみの関係を示す。 図１７で得られた応力−ひずみ応答の各ポテンシャル毎の寄与分を示すグラフである。 図１７の１サイクル目の応力変化を、圧縮ノードに属する結合長に生じる応力、及びそれ以外（非圧縮ノード）の結合長に生じる応力にそれぞれ分けて調べた結果を示すグラフである。 （ａ）は高圧縮モデル、（ｂ）は低圧縮モデルの一例をそれぞれ視覚化して示す斜視図である。

符号の説明

２ 分子鎖モデル
３ 粒子
４ 結合鎖
Ｓ セル

Claims (5)

1. 分子動力学計算を用いて高分子材料のシミュレーションを行う方法であって、
   解析対象の高分子材料の分子構造に基づき、任意の空間を有するセル内に、原子又はその集合体を表す有限個の粒子と、この粒子間を結合して各粒子の相対位置を定義する結合鎖とからなる三次元構造の分子鎖モデルを配置した分子構造モデルを設定するステップ、
   前記分子構造モデルを分子動力学計算に基づいて安定化させる構造緩和のシミュレーションを行うステップ、
   前記構造緩和後の分子鎖モデルを分子鎖部分に分解しかつそれぞれの分子鎖部分に生じる応力を計算する応力計算ステップ、及び
   前記応力計算ステップで計算された各分子鎖部分の応力に基づいてエネルギーロスの主要な発生箇所を決定する決定ステップを含むことを特徴とする高分子材料のシミュレーション方法。
2. 前記決定ステップは、圧縮応力が大きい分子鎖部分をエネルギーロスの主要な発生箇所として決定する請求項１記載の高分子材料のシミュレーション方法。
3. 前記応力計算ステップは、結合鎖を介して連続する４粒子が構成する立体配座に基づいて前記分子鎖モデルを分子鎖部分に分解する請求項１又は２記載の高分子材料のシミュレーション方法。
4. 前記応力計算ステップは、前記分子鎖モデルを２粒子の結合長に基づいて分解しかつそれぞれの結合長を中心軸とする立体配座毎に区別して応力を計算する第１の計算ステップと、
   この第１の計算ステップの結果に基づいて、大きな圧縮応力が生じる結合長−立体配座の組み合わせを特定するステップと、
   前記特定された結合長−立体配座の組み合わせについて、その両端に結合される粒子をさらに考慮して分子鎖モデルを３つの連続する立体配座の並び別にグループ化しかつ各グループ毎の応力を計算する第２の計算ステップとを含むとともに、
   前記決定ステップは、前記第２の計算ステップの結果に基づいて、圧縮応力が大きいグループを、エネルギーロスの主要な発生箇所として特定する請求項１記載の高分子材料のシミュレーション方法。
5. 前記エネルギーロスの主要な発生箇所として決定された部分を他の部分と差別化して表示するステップをさらに含む請求項１乃至４のいずれかに記載の高分子材料のシミュレーション方法。