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

Abstract

【課題】高分子材料のエネルギーロス発生箇所を分子レベルで特定しうるシミュレーション方法を提供する。
【解決手段】分子動力学計算を用いて高分子材料のシミュレーションを行う方法であって、前記高分子材料の分子構造を、原子又はその集合体を表す有限個の粒子と、この粒子間の相対位置を特定する結合鎖とを含んでモデル化したシミュレーション用の分子構造モデルを設定するステップＳ１、前記分子構造モデルに、任意の軸方向の歪を増大させる負荷変形工程と、前記歪を減少させる除荷変形工程とを含む変形シミュレーションを行うステップＳ２、及び特定の歪状態における前記負荷変形時の前記結合鎖の前記軸方向に対する角度αi と、同一の歪状態における前記除荷変形時の前記結合鎖の前記角度α'iとの差に基づいてエネルギーロスの発生箇所を特定するステップ（Ｓ３からＳ６）を含むことを特徴とする。
【選択図】図１

Description

本発明は、高分子材料のエネルギーロス発生箇所を分子レベルで特定するのに役立つシミュレーション方法に関する。

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

本発明は、以上のような実情に鑑み案出なされたもので、コンピュータを用いて分子動力学計算に基づいた高分子材料のシミュレーションを用いて行い、その結果からエネルギーロスの発生箇所を特定することを基本として、能率的な高分子材料の設計ないし開発を可能としうる高分子材料のシミュレーション方法を提供することを目的としている。

本発明のうち請求項１記載の発明は、分子動力学計算を用いて高分子材料のシミュレーションを行う方法であって、前記高分子材料の分子構造を、原子又はその集合体を表す有限個の粒子と、この粒子間の相対位置を特定する結合鎖とを含んでモデル化したシミュレーション用の分子構造モデルを設定するステップ、前記分子構造モデルに、任意の軸方向の歪を増大させる負荷変形工程と、前記歪を減少させる除荷変形工程とを含む変形シミュレーションを行うステップ、及び特定の歪状態かつ前記負荷変形時における前記軸方向に対する前記結合鎖の角度αi と、同一の歪状態かつ前記除荷変形時の前記結合鎖の角度α'iとの差に基づいてエネルギーロスの発生箇所を特定するステップを含むことを特徴とする。

また請求項２記載の発明は、前記エネルギーロスの発生箇所を特定するステップは、前記角度の差（α'i−αi ）が正かつ大きい結合鎖をエネルギーロスの発生箇所として特定する請求項１記載の高分子材料のシミュレーション方法である。

また請求項３記載の発明は、前記エネルギーロスの発生箇所として特定された結合鎖を他の結合鎖と差別化して、前記分子構造モデルを表示するステップをさらに含む請求項１又は２に記載の高分子材料のシミュレーション方法である。

本発明では、分子動力学を用いたシミュレーションに基づいて、高分子材料のエネルギーロスの発生箇所を原子又は分子レベルで特定できる。従って、高分子材料の開発を能率的に行うことができる。

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

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

前記高分子材料は、高分子を含む材料であれば特に限定されないが、好ましくは天然ゴム、合成ゴムゴム及び合成樹脂を少なくとも含む有機高分子材料が好適である。本実施形態では、高分子材料としてアモルファスポリエチレンを用いて以後の説明が行われる。

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

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

図２の上段にはポリエチレンの分子構造が示され、その下段には、本実施形態のポリエチレンの分子構造モデル２が視覚化されて略示される。分子構造モデル２は、原子又はその集合体を表す複数かつ有限個の粒子３と、この粒子３、３間の相対位置を特定する結合鎖４とを含んでモデル化されている。

本実施形態において、ポリエチレンの水素原子は陽に扱わず、メチレン基（−ＣＨ2 −）を原子の集合体として一つの粒子３として扱ういわゆる "United Atom Model"が用いられる。水素原子からの系への影響が比較的小さいため、このような水素原子の効果を炭素原子に取り込んだ "United Atom Model"を用いても計算結果への影響は少ないと考えられる。ただし、全ての原子をそれぞれ粒子３として捉えるいわゆる "Full Atom model"によってモデル化が行われても良いのは言うまでもない。このように、分子構造モデル２の粒子３は、原子自体又はその集合体を表し得る。

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

また、分子構造モデル２は、３次元構造を有し、図３（ａ）〜（ｃ）に示されるように、各粒子３、３間の結合（「ボンド」とも呼ばれる。）の長さである結合長ｒ、同図（ｂ）に示されるように、隣り合う３つの粒子３がなす角度である結合角θ、同図（ｃ）に示されるように、隣り合う４つの粒子３において，隣り合う３つの粒子が作る第１の平面Ｐ１と、その中の２つの粒子が共通する３つの粒子３が作る第２の平面Ｐ２のなす角度であるトーションφが定義される。

前記分子構造モデル２を設定する手順は、特に限定されるものではないが、本実施形態では次のような手順を経て行われる。先ず、図４に示されるように、分子構造モデル２の開始点となる一つの粒子３ａが、三次元空間内かつ乱数によって決定されたランダムな座標位置に定められる。次に、粒子３、３間の結合長ｒを初期値ｒ0 （＝０．１５３３nm）及び結合角を初期値θ0 （１１３．３deg ）にそれぞれ固定して、順次、粒子３ｂ、３ｃ、３ｄ…をランダムに配置することにより、分子構造モデル２を成長させる。この際、トーションφには、１８０゜（平坦を意味するものとして、以下、"trans" なる用語を用いる場合がある。）、６７．５゜又は２９２．５゜（いずれも、非平坦を意味するものとして、以下、"gauche"なる用語を用いる場合がある。）の３種類のうちの一つがランダムに選択される。

このように分子構造モデル２を順次成長させることにより、図５で視覚的に示されるように、３次元構造を有したランダムコイル状の分子構造モデル２が設定される。

また、本実施形態では、複数本の分子構造モデル２が設定され、それらは図６に示されるような予め定められた体積を持ったセルＳの中に配置される。該セルＳは、解析対象のポリエチレンの微小部分に相当する。本実施形態では、セルＳは、２０nm×２０nm×２０nmの微小な立方体として定義される。

また、各分子構造モデル２の長さ（ここで言う長さは、１本の分子構造モデル２に含まれる粒子３の数とし、以下同様である。）は、後述の変形シミュレーションを行うのに足りる十分な長さであれば特に限定はされない。本実施形態では、前記セルＳの中に、長さが異なる複数種類の分子構造モデル２が定義される。また、各分子構造モデル２の長さは、例えば、その平均長さが３００、かつ、図７に示されるように標準偏差が５０の正規分布となるように定めることができ。

また、もし分子構造モデル２を成長させる途中において、セルＳの壁に遮られ目的とする長さの分子構造モデル設定できなくなった場合、当該分子構造モデルを破棄し、開始点を新たに設定して新たに分子構造モデル２が設定される。これらの処理を繰り返すことにより、セルＳは、複数本の分子構造モデル２で満たされる。そして、セルＳが、目的とするポリエチレンの密度になるように、分子構造モデル２の本数ないし粒子数が微調整される。本実施形態では、セルＳに配置された分子構造モデル２の本数は９２６本、それらの粒子３の総数は２７７３３７個及びセルＳの密度は０．８５ｇ／cm³ とした。

また、セルＳは、その８つの壁面Ｓａを含みかつ該壁面Ｓａから内側に１nmの厚さｔを有する外殻部Ｓｏと、該外殻部Ｓｏで囲まれるコア部Ｓｉとに仮想区分される。本実施形態では、前記外殻部Ｓｏに一部が含まれる粒子３には、移動不能に拘束される条件が定義される。これは、セルＳに例えば１軸方向の引張変形等を与えた場合、セルＳの変形が歪になって現実の変形状態から大きく逸脱するのを防止するのに役立つ。

次に、本実施形態では、設定されたセルＳ（分子構造モデル２）と分子動力学計算とに基づいて変形シミュレーションが行われる（ステップＳ２）。分子動力学計算では、分子構造モデル２に対して、図３（ａ）で示される連結された粒子３、３間のボンドストレッチのポテンシャル（略記号：Ｅbs）、同図（ｂ）の連続する３つ粒子３で構成されるベンディングのポテンシャル（略記号：Ｅbe）、同図（ｃ）のトーションφのポテンシャル（略記号Ｅto）及び図３（ｄ）で示されるように、互いに連結されていない粒子間のファンデルワールス力のポテンシャル（略記号：Ｅvw）の相互作用が考慮される。これらのポテンシャル関数は、種々のものが採用できるが、本実施形態では、下記式（１）〜（５）が採用される。また、それらの各パラメータとしては、表１に示される値が採用される。

変形シミュレーションでは、先ず、設定されたセルＳについて、構造緩和が行われる。この構造緩和は、初期設定されたセルＳについて、所定の時間（本実施形態では、１０００fs）の分子動力学計算を行うことにより、セルＳの初期構造の影響を無くして、その構造を安定化ないし準安定化させる。この構造緩和の条件として、積分の時間ステップや温度等は任意に定めることができる。本実施形態では、積分の時間ステップは０．１fs、温度は速度スケーリング法によって３００Ｋに制御される。

次に構造緩和が行われたセルＳに対し、任意の軸方向の歪を増大させる負荷変形工程と、前記歪を減少させる除荷変形工程とが行われる。本実施形態では、セルＳに、Ｚ軸方向の引張ひずみが与えられる。前記負荷変形工程では、Ｚ軸方向のセルＳのひずみが０．５になるまで、一定のひずみ速度（例えば１．０×１０^-11 ／s ）でセルＳを引き延ばすことにより行われる。また、除荷変形工程では、セルＳのＺ軸方向のひずみが、前記負荷変形工程と同一のひずみ速度で減少するように行われる。なお前記ひずみは、構造緩和後のセルＳのＺ軸方向の長さを基準とする工学ひずみである。また、各粒子３のＺ座標は、下記式（６）のように、毎ステップスケールすることでひずみを与えることができる。ここで、Δε_zzは、例えば１．０×１０^-5とした。また、この分子動力学計算によるシミュレーションにおいても、速度スケーリングによってセルＳの温度が３００Ｋに制御された。

図８には、このような変形シミュレーションにより得られたセルＳの応力−ひずみ線図が示される。この応力は、前記コア部Ｓｉに座標値が含まれる移動可能な粒子３に作用する応力ベクトルを、各粒子３に作用する力の中で、正のものだけを加算することにより求め、前記コア部Ｓｉの体積で除した真応力である。

図８を参照すると、初期設定された密度０．８５ｇ／cm³ のセルＳでは、基準状態において約１５０ＭＰａの静水圧圧縮状態にあることが分かる。また、ひずみの増加とともに引張応力は非線形に増加し、ひずみが０．５になった時点で除荷を開始すると、極めて短時間で応力が約１００ＭＰａから４５ＭＰａまで大きく減少することが分かる。また、ひずみを０まで戻した際、応力は約−２０６ＭＰａとなり、初期状態よりも大きな値を示している。この応力は、ひずみを０に保持した状態のものではなく、動的な除荷過程におけるひずみの０の瞬間の値のため、このような値が得られたと推察される。そして、これらの負荷変形及び除荷変形を経ることにより、応力とひずみとの間にヒステリシスロスが生成される。

また、原子（粒子３）に作用する力は、前記式（１）の空間微分から得られる。従って、応力ベクトルを評価する際、式（２）〜（５）から結合長（ボンドストレッチ）、結合角（ベンディング）、二面角（トーション）及びファンデルワールスのポテンシャルへの寄与を明らかにすることができる。図９には、図８の応力をポテンシャルの各成分に分けて示す。

図９から明らかなように、図８のヒステリシスロス（エネルギーロス）は、大部分がボンドストレッチによって担われており、残りはファンデルワールスであることが判明した。ただし、ファンデルワールスの寄与は引張によって零に漸近しており、ヒステリシスには圧縮としてのみ寄与することがわかる。ベンディングは、負荷変形時および除荷変形時に、それぞれ僅かに一定の正およぴ負の応力を担い、ボンドストレッチと同様、負荷反転時に急速な応力減少を示している。さらに、トーションは、変形中、常にほぼ零であり、系の応力変化に寄与しないことが明らかである。これらは、いずれもヒステリシスロスに大きな影響を与えていない。

発明者らは、これらの結果を踏まえ、高分子材料のエネルギーロス発生箇所を特定するためにさらなる研究を行った。先ず、図１０には、ポリエチレンのメチレン基のトーションφのポテンシャル関数が示される。トーションφのポテンシャルは、１８０゜（trans ）で最安定状態となり、６７．５又は２９２．５゜（gauche）で準安定状態になる。発明者らは、前記変形シミュレーションの結果から、負荷変形工程及び除荷変形工程からなる１サイクルについて、全ての分子構造モデル２のノード（トーションを構成する連続する３つの粒子を１組とするもの）のトーションの変化を調べた。具体的には、トーションφが０゜≦φ≦１８０゜の範囲のノードを対象として評価し、前記トーションが１６０゜≦φ≦１８０゜のものを transノード、同４７．５゜≦φ≦８７．５゜のものを gauche ノード、その他のものを high energyノードとして分類し、初期構造におけるノード数Ｎ0 からの変化量（Ｎ−Ｎ₀ ）をそれぞれ調べた。

図１１にはその結果が示され、横軸には時間とひずみが、また縦軸にはノードの変化量（Ｎ−Ｎ₀ ）が示される。また、時間が０〜５０００fsの間は負荷変形工程、５０００〜１００００fsの間は除荷変形工程であることを示している。

図１１から明らかなように、引張初期には、high energy ノードが急減し、それに対応してgaucheノードが増加している。これは、引張を駆動力として、セルＳの内部が構造緩和されたことを意味する。一方、応カ−ひずみの勾配が緩やかになるひずみが０．１５（ｔ＝１５００fs）近傍で、 gauche ノードの増加及び high energyノードの減少が止まり、以降は gauche ノードが減少している。trans ノードは、負荷変形工程では、最大ひずみまで単調に増加している。これは、１軸の引張変形における負荷変形工程では、その軸方向に沿った直線状のノードが増えることを示唆している。

また、ｔ＝５０００fsでひずみ増分Δε_zzが反転された場合、trans ノードが急減し、それに対応して gauche ノードが急増することが分かる。即ち、除荷変形工程では、trans ノードから gauche ノードへの遷移が生じ、分子構造モデル２の中に多数の折れ曲がり部分が導入される。さらに、ｔ＝１００００fs、すなわちひずみが０に戻った時点では、初期構造よりも gauche ノードが増えており、trans ノードが少ない状態にある。以上より、前記１サイクルの変形を行った後のセルＳの構造は、分子構造モデルに多くの“たるみ”（あるいは折れ曲がり）が含まれていることが示唆される。

また、図１２には、変形シミュレーションの結果から得られた、引張の軸方向とのなす角度が２０゜以下である結合鎖（以下、このような結合鎖を「軸方向に配向された結合鎖」と呼ぶ。）の数の変化量を初期構造からの変化量（Ｎ−Ｎ₀ ）を示す。先の解析では、負荷変形工程においてtrans ノードが増加することが判明したが、図１２からも明らかなように、それは軸方向に配向された結合鎖の増加によってもたらされている。このため、ボンドストレッチが担う応力が増加する。しかし、除荷変形工程後は、先で述べた通り、分子構造モデルが弛むことにより、初期構造よりも軸方向に配向された結合鎖が減少し、分子鎖モデル２の中に多数の折れ曲がりないし弛みが導入される。

従って、高分子材料のヒステリシスは、図１３に模式的に示されるように、分子構造モデルを、負荷変形時及び除荷変形時の同一の歪εa の状態で比較した場合、除荷変形時に引張軸方向に対して大きく弛んでいる部分が原因であると特定できる。そこで、本実施形態では、特定の歪状態かつ負荷変形時の引張の軸方向に対する結合鎖４の角度αi と、同一の歪状態かつ除荷変形時の結合鎖４の角度α'iとの差（α'i−αi ）に基づいてエネルギーロスの発生箇所を特定する処理を行う。より具体的には、全ての結合鎖４について、前記角度の差（α'i−αi ）を計算し、その値が正かつ大きい値を示す結合鎖４をエネルギーロスの発生箇所として特定することができる（ステップＳ３〜Ｓ６）。前記各角度αi 及びαi'は、変形シミュレーションの計算結果を用い、例えば引張の軸方向に対する結合鎖４のベクトルのなす角度を求めることで容易に得ることができる。

図１４には、前記変形シミュレーションから１本の分子構造モデル２の一部が取り出されかつ視覚化して示される。左側は負荷変形工程時の特定の歪状態、右側はそれと同一の歪の除荷変形工程時の状態である。各々の角度は、引張の軸方向に対する角度であり、カッコ内には、各々の結合鎖４について前記角度の差（α'i−αi ）が表示される。この例では、符号Ｒで示される結合鎖を含む部分の前記角度差が＋２６゜で最も大きい。従って、少なくともこの結合鎖は、エネルギーロスが大きく発生している箇所として特定される。また、エネルギーロス発生箇所として、その両端の粒子３を含めても良い。ここで、エネルギーロス発生箇所として特定するための前記角度の差（α'i−αi ）の閾値は、材料などに応じてまた、全体のエネルギーロスなどに応じて任意に定めることができる。また、それらは相対的なものとして捉えても良い。

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

上記シミュレーションの結果より、エネルギーロスの発生箇所を特定できるので、その位置に例えばベンゼン環などのような官能基を付けることエネルギーロスを低減させることができる。また、エネルギーロスの大小と、分子鎖の長さ、分子鎖の大きさ及び／又は分子鎖の量との関係などを得ることができ、製品開発に大いに役立つ。特に複雑な分子構造を有するポリマーの場合、該ポリマーを構成するモノマーのどの部分でエネルギーロスが生じているのかを調べることができる。従って、その分子構造を変えることでエネルギーロスがより小さくなるような分子構造設計が可能になる。さらに、本実施形態のシミュレーション方法は、例えばゴム等の配合設計では、ゴムのエネルギーロスの検証ツールとしても利用できる。このようにして得られた新規な配合のゴムで例えばタイヤ等の試作、実解析を行うことで効率の良い製品設計ができる。

以上の実施形態では、ポリエチレンというメチレン基が同じ結合で並んだ材料を用いたが、例えばゴムのように部分的に二重結合を含む場合や、官能基が結合しているような高分子材料にも本発明は適用できる。

本発明の処理手順を示すフローチャートである。 ポリエチレンの分子構造とその分子構造モデルを説明する線図である。 （ａ）〜（ｄ）はポテンシャルを説明する分子構造モデルの部分図である。 分子構造モデルの設定手順を説明する概念図である。 分子構造モデルを視覚化して示す斜視図である。 セルの斜視図である。 セルに配置される分子構造モデルの分布図である。 変形シミュレーションから得られた応力−ひずみ曲線である。 図８の応力をポテンシャル毎に分けた応力−ひずみ曲線である。 トーションのポテンシャル関数である。 ノードの変化量と時間との関係を示すグラフである。 ノードの変化量と時間との関係を示すグラフである。 ヒステリシスのメカニズムを模式的に示す線図である。 分子構造モデルの一部を視覚化して示す図である。

符号の説明

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

Claims (3)

1. 分子動力学計算を用いて高分子材料のシミュレーションを行う方法であって、
   前記高分子材料の分子構造を、原子又はその集合体を表す有限個の粒子と、この粒子間の相対位置を特定する結合鎖とを含んでモデル化したシミュレーション用の分子構造モデルを設定するステップ、
   前記分子構造モデルに、任意の軸方向の歪を増大させる負荷変形工程と、前記歪を減少させる除荷変形工程とを含む変形シミュレーションを行うステップ、及び
   特定の歪状態かつ前記負荷変形時における前記軸方向に対する前記結合鎖の角度αi と、同一の歪状態かつ前記除荷変形時の前記結合鎖の角度α'iとの差に基づいてエネルギーロスの発生箇所を特定するステップを含むことを特徴とする高分子材料のシミュレーション方法。
2. 前記エネルギーロスの発生箇所を特定するステップは、前記角度の差（α'i−αi ）が正かつ大きい結合鎖をエネルギーロスの発生箇所として特定する請求項１記載の高分子材料のシミュレーション方法。
3. 前記エネルギーロスの発生箇所として特定された結合鎖を他の結合鎖と差別化して、前記分子構造モデルを表示するステップをさらに含む請求項１又は２に記載の高分子材料のシミュレーション方法。

Images