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

Description

本発明は、緩和弾性率を計算するための高分子材料のシミュレーション方法に関する。

近年、ゴム、樹脂等のエラストマーを含む高分子材料の開発のために、高分子材料の性質を、コンピュータを用いて評価するためのシミュレーション方法（数値計算）が種々提案されている（例えば、下記特許文献１参照）。

この種のシミュレーション方法では、先ず、高分子材料の分子鎖をモデル化した分子鎖モデルが設定される。次に、予め定められたセル内に分子鎖モデルが配置され、高分子材料モデルが設定される。そして、分子動力学（ Molecular Dynamics ： ＭＤ ）に基づいて、高分子材料モデルの緩和計算が行われる。

特開２０１３−１９５２２０号公報

上記のようなシミュレーション方法を利用して、例えば、緩和弾性率が計算される場合がある。緩和弾性率とは、歪が与えられた粘弾性体の弾性率の変化を示す指標である。一般に、高分子材料の緩和弾性率の計算には多くの時間が必要とされている。

ところで、複数種類の高分子成分が混合されてなる混合系高分子材料においても、セル内に複数種類の分子鎖モデルを配置した高分子材料モデルを設定することにより、混合系高分子材料の全体の緩和弾性率を計算することが可能である。

しかしながら、上述した手法により混合系高分子材料の緩和弾性率を計算する場合、各高分子成分の混合比率すなわち混合系高分子材料の全体に対する各高分子成分の体積分率が変更されるたびに、高分子材料モデルを再設定し、緩和弾性率を再計算する必要があることから、膨大な時間を要していた。

本発明は、以上のような実状に鑑み案出されたもので、混合系高分子材料の緩和弾性率を短時間で計算できる高分子材料のシミュレーション方法を提供することを主たる目的としている。

本発明は、複数種類の高分子成分が混合されてなる混合系高分子材料の全体の緩和弾性率Ｇ_all（ｔ）を予測する高分子材料のシミュレーション方法であって、各高分子成分毎に、単体での緩和弾性率Ｇ_i（ｔ）を計算する第１工程と、前記各高分子成分毎に、単位分子量における緩和時間τ_iを計算する第２工程と、前記第２工程によって計算された前記緩和時間τ_iを用いて、下記式（１）で、前記各高分子成分毎に前記緩和弾性率Ｇ_i（ｔ）の重み付けパラメータＫ_i を計算する第３工程と、
前記各高分子成分毎に、前記混合系高分子材料の全体に対する体積分率φ_iを設定する第４工程と、前記第１工程によって計算された前記緩和弾性率Ｇ_i（ｔ）、前記第３工程によって計算された前記重み付けパラメータＫ_i 及び前記第４工程によって設定された前記体積分率φ_iを用いて、下記式（２）で、前記混合系高分子材料の全体の緩和弾性率Ｇ_all（ｔ）を計算する第５工程とを含むことを特徴とする。
Ｋ_i ＝（１／τ_i）／Σ（１／τ_i ） （１）
ただし、
i ：各高分子成分に割り当てられた整数番号
Ｋ_i ：第i番目の高分子成分の重み付けパラメータ
τ_i ：第i番目の高分子成分の単体での緩和時間
Ｇ_all（ｔ）＝Σ（Ｋ_i×φ_i×Ｇ_i（ｔ）） （２）
ただし、
Ｇ_all（ｔ） ：混合系高分子材料の全体の緩和弾性率
φ_i ：第i番目の高分子成分の混合系高分子材料全体に対する
体積分率
Ｇ_i（ｔ） ：第i番目の各高分子成分の単体での緩和弾性率

本発明に係る前記高分子材料のシミュレーション方法において、前記各高分子成分の単体での緩和弾性率Ｇ_i（ｔ）は、前記各高分子成分がガラス状態にあるガラス領域での緩和弾性率Ｇ_i（ｔ）を含み、前記第１工程において、少なくとも前記ガラス領域での前記緩和弾性率Ｇ_i（ｔ）は、全原子モデル又はユナイテッドアトムモデルでモデル化された分子鎖モデルを用いた分子動力学によって計算されることが望ましい。

本発明に係る前記高分子材料のシミュレーション方法において、前記第１工程は、５００Ｋ以下の温度条件の下で行なわれることが望ましい。

本発明に係る前記高分子材料のシミュレーション方法において、前記各高分子成分の単体での緩和弾性率Ｇ_i（ｔ）は、前記各高分子成分がゴム状態にあるゴム領域での緩和弾性率Ｇ_i（ｔ）を含み、前記第１工程において、少なくとも前記ゴム領域での前記緩和弾性率Ｇ_i（ｔ）は、管模型でモデル化された分子鎖モデルを用いて計算されることが望ましい。

本発明に係る前記高分子材料のシミュレーション方法において、前記各高分子成分の単体での緩和弾性率Ｇ_i（ｔ）は、前記ガラス領域と前記ゴム領域との中間の領域である転移領域での緩和弾性率Ｇ_i（ｔ）を含み、前記第１工程において、少なくとも前記転移領域の前記緩和弾性率Ｇ_i（ｔ）は、粗視化モデルでモデル化された分子鎖モデルを用いた分子動力学によって計算されることが望ましい。

本発明に係る前記高分子材料のシミュレーション方法において、前記第１工程では、１つのセル内に８本以上の前記分子鎖モデルを配置した材料モデルを用いて周期境界条件の下で、分子動力学計算を行なう工程を含むことが望ましい。

本発明に係る前記高分子材料のシミュレーション方法において、前記第２工程では、前記緩和時間τ_iは、全原子モデル又はユナイテッドアトムモデルでモデル化された分子鎖モデルを用いた分子動力学によって計算されることが望ましい。

本発明に係る前記高分子材料のシミュレーション方法において、前記分子鎖モデルの重合度は、５０以下であることが望ましい。

本発明のシミュレーション方法は、各高分子成分毎に、単体での緩和弾性率Ｇ_i（ｔ）を計算する第１工程と、各高分子成分毎に、単体での単位分子量における緩和時間τ_iを計算する第２工程と、第２工程によって計算された緩和時間τ_iを用いて、下記式（１）で、各高分子成分毎に緩和弾性率の重み付けパラメータＫ_i を計算する第３工程とを含む。
Ｋ_i ＝（１／τ_i）／Σ（１／τ_i ） （１）
ただし、
i ：各高分子成分に割り当てられた整数番号
Ｋ_i ：第i番目の高分子成分の重み付けパラメータ
τ_i ：第i番目の高分子成分の単体での緩和時間
そして、本発明のシミュレーション方法は、各高分子成分毎に、混合系高分子材料の全体に対する体積分率φ_iを設定する第４工程と、第１工程によって計算された前記緩和弾性率Ｇ_i（ｔ）、第３工程によって計算された重み付けパラメータＫ_i 及び第４工程によって設定された前記体積分率φ_iを用いて、下記式（２）で、混合系高分子材料の全体の緩和弾性率Ｇ_all（ｔ）を計算する第５工程をさらに含む。
Ｇ_all（ｔ）＝Σ（Ｋ_i×φ_i×Ｇ_i（ｔ）） （２）
ただし、
Ｇ_all（ｔ） ：混合系高分子材料の全体の緩和弾性率
φ_i ：第i番目の高分子成分の混合系高分子材料全体に対する
体積分率
Ｇ_i（ｔ） ：第i番目の各高分子成分の単体での緩和弾性率
すなわち、第５工程では、第１工程によって計算された緩和弾性率Ｇ_i（ｔ）、第３工程によって計算された重み付けパラメータＫ_i 及び第４工程によって設定された体積分率φ_iが下記式（１）に代入され、混合系高分子材料の全体の緩和弾性率Ｇ_all（ｔ）が計算される。

上記式（２）において、各高分子成分の単体での緩和弾性率Ｇ_i（ｔ）は、各高分子成分の緩和時間に基づく重み付けパラメータＫ_i 及び各高分子成分の体積分率φ_iを考慮したうえで処理され、混合系高分子材料の全体の緩和弾性率Ｇ_all（ｔ）が近似的に計算されうることが、発明者によって見出された。

ここで、第１工程によって計算された緩和弾性率Ｇ_i（ｔ）は、各高分子成分毎に固有のものであり、各高分子成分の混合比率に依存しない。従って、各高分子成分の混合比率を変更する場合にあっても、変更前の第１工程で既に算出している緩和弾性率Ｇ_i（ｔ）をそのまま用いることができる。第２工程によって計算された緩和時間τ_i及び第３工程によって計算された重み付けパラメータＫ_i についても同様である。これにより、第２工程及び第３工程と比較して、膨大な計算時間を必要とする第１工程乃至第３工程を再度実行することなく、混合比率を変更した混合系高分子材料の全体の緩和弾性率Ｇ_all（ｔ）が計算可能となる。

本発明のシミュレーション方法を実行するためのコンピュータの斜視図である。 本発明のシミュレーション方法を用いて緩和弾性率が計算される混合系高分子材料モデルの１つのセルを示す斜視図である。 ポリスチレンの構造式である。 分子鎖モデルの概念図である。 混合系高分子材料の全体の緩和弾性率を示すグラフである。 本実施形態のシミュレーション方法の処理手順の一例を示すフローチャートである。 図７の第１工程の処理手順の一例を示すフローチャートである。 混合系高分子材料を構成する第１高分子成分の高分子材料モデルの１つのセルを示す斜視図である。 図７の第１工程で計算された各高分子成分の単体での緩和弾性率を示すグラフである。 末端間ベクトルの長さの時間相関関数と時間との関係を示すグラフである。 本実施形態のシミュレーション方法によって計算された混合系高分子材料の全体の緩和弾性率と、従来のシミュレーション方法によって計算された混合系高分子材料の全体の緩和弾性率とを比較するグラフである。 緩和弾性率のガラス領域、転移領域及びゴム領域を示すグラフである。

以下、本発明の実施の一形態が図面に基づき説明される。
本実施形態の高分子材料のシミュレーション方法（以下、単に「シミュレーション方法」ということがある）は、コンピュータを用いて、高分子材料の緩和弾性率を計算するための方法である。

図１は、本発明のシミュレーション方法を実行するためのコンピュータ１の斜視図である。コンピュータ１は、本体１ａ、キーボード１ｂ、マウス１ｃ及びディスプレイ装置１ｄを含んでいる。この本体１ａには、例えば、演算処理装置（ＣＰＵ）、ＲＯＭ、作業用メモリ、磁気ディスクなどの記憶装置、及びディスクドライブ装置１ａ１、１ａ２が設けられる。また、記憶装置には、本実施形態のシミュレーション方法を実行するための処理手順（プログラム）が予め記憶される。

高分子材料としては、例えば、ゴム、樹脂等のエラストマーが含まれる。本実施形態の高分子材料は、複数種類の高分子成分が混合されてなる混合系高分子材料である。

図２は、本実施形態のシミュレーション方法を用いて緩和弾性率が計算される混合系高分子材料モデル１０の１つのセル（空間）１１を示している。混合系高分子材料モデル１０の１つのセル１１内には、複数種類の高分子成分の分子鎖モデル２が配置されている。複数種類の高分子成分の例としては、ポリスチレンの含有量が異なるＳＢＲ（スチレン・ブタジエンゴム）が挙げられる。より具体的には、ポリスチレンの含有量が２０％のＳＢＲとポリスチレンの含有量が２５％のＳＢＲとが、予め定められた混合比率に従って配置されている。

本実施形態のセル１１は、互いに向き合う三対の平面１２、１２を有する立方体として定義されている。各平面１２には、周期境界条件が定義されている。これにより、セル１１では、例えば、一方の平面１２ａから出て行った分子鎖モデル２の一部が、反対側の平面１２ｂから入ってくるように計算されうる。従って、一方の平面１２ａと、反対側の平面１２ｂとが連続している（繋がっている）ものとして扱われる。

図３は、ポリスチレンの構造式である。このポリスチレンを構成する分子鎖Ｍｃは、スチレンから構成されるモノマーが、重合度（分子量）Ｍｎで連結されている。なお、本発明のシミュレーション方法では、ポリスチレン以外の高分子材料が用いられてもよい。

本実施形態のシミュレーション方法では、コンピュータ１に、高分子材料の分子鎖Ｍｃ（図２に示す）をモデル化した分子鎖モデル２が設定され、その緩和弾性率が計算される。

図４は、本実施形態の分子鎖モデル２の概念図である。分子鎖モデル２は、例えば、複数の粒子モデル３と、粒子モデル３、３間を結合するボンドモデル４とを含む全原子モデルとして構成されている。これらの粒子モデル３及びボンドモデル４は、分子鎖Ｍｃのモノマーをなす単位構造６（図３に示す）に基づいて、互いに連結されることにより、モノマーモデル７が設定される。このモノマーモデル７が、重合度Ｍｎに基づいて連結されることにより、分子鎖モデル２が設定される。

粒子モデル３は、後述する分子動力学計算に基づいたシミュレーションにおいて、運動方程式の質点として取り扱われる。即ち、粒子モデル３は、質量、電荷、又は、初期座標などのパラメータが定義される。本実施形態の粒子モデル３は、分子鎖Ｍｃの炭素原子をモデル化した炭素粒子モデル３Ｃ、及び、分子鎖Ｍｃの水素原子をモデル化した水素粒子モデル３Ｈを含んでいる。

ボンドモデル４は、粒子モデル３、３間を拘束するものである。本実施形態のボンドモデル４は、炭素粒子モデル３Ｃ、３Ｃを連結する主鎖４ａ、及び、炭素粒子モデル３Ｃと水素粒子モデル３Ｈとの間を連結する側鎖４ｂとを含んでいる。

粒子モデル３、３間には、相互作用（斥力及び引力を含む）が生じさせるポテンシャルが定義される。ポテンシャルには、ボンドモデル４を介して隣り合う粒子モデル３、３間に定義される第１ポテンシャルＰ１と、ボンドモデル４を介さずに隣り合う粒子モデル３、３間に定義される第２ポテンシャルＰ２とが定義される。なお、複数の分子鎖モデル２が定義される場合は、分子鎖モデル２、２間の粒子モデル３、３間にも、第２ポテンシャルＰ２が定義される。

図５は、本実施形態のシミュレーション方法を用いて計算される混合系高分子材料の全体の緩和弾性率Ｇ_all（ｔ）である。緩和弾性率Ｇ_all（ｔ）は、時間ｔの関数である。図５では、緩和弾性率Ｇ_all（ｔ）の対数をとった値ｌｏｇＧ_all（ｔ）と、時間ｔの対数をとった値ｌｏｇｔとの関係が示されている（以下、図９、１１、１２においても同様）。混合系高分子材料がタイヤに用いられる場合、緩和弾性率Ｇ_all（ｔ）は、タイヤのグリップ性能や転がり抵抗性能と関わりがあるとして注目されている。

混合系高分子材料の場合、各高分子成分の混合比率によって全体の緩和弾性率Ｇ_all（ｔ）は変化する。すなわち、発明者らは、混合比率の大きい高分子成分は、全体の緩和弾性率Ｇ_all（ｔ）に及ぼす影響が大きいと考えた。

緩和弾性率Ｇ（ｔ）は、所定の時間幅ｔにおいて、歪が与えられた粘弾性体の弾性率の変化を示す指標である。発明者らは、緩和弾性率Ｇ（ｔ）は、外力によって変形された分子鎖が、分子の運動によって元の自然平衡な配置状態に戻る過程に関する物理量であると考えた。一方、分子鎖の形状が元に戻ることを分子レベルで定量的に表すことができる物理量として分子の緩和時間が存在する。そこで、発明者らは、緩和弾性率Ｇ（ｔ）と緩和時間との関係について研究を重ねたところ、混合系高分子材料の場合、各高分子成分の単体での単位分子量における緩和時間τ_iも、全体の緩和弾性率Ｇ_all（ｔ）に影響を及ぼすことが判明した。例えば、緩和時間τ_aが大きい高分子成分ａと緩和時間τ_bが小さい高分子成分ｂとが、混合系高分子材料に含まれる場合、緩和時間τ_bが小さい高分子成分ｂの方が全体の緩和弾性率Ｇ_all（ｔ）に及ぼす影響が大きい。

そして発明者らは、さらなる研究を重ねた結果、ついに下記式（１）によって計算される各高分子成分の緩和時間に基づく重み付けパラメータＫ_i を創出した。
Ｋ_i ＝（１／τ_i）／（１／τ₁ ＋１／τ₂ ＋…＋１／τ_i）
＝（１／τ_i）／Σ（１／τ_i ） （１）
ただし、
i ：各高分子成分に割り当てられた整数番号
Ｋ_i ：第i番目の高分子成分の重み付けパラメータ
τ_i ：第i番目の高分子成分の単体での緩和時間
なお、上記iは、１から始まり、混合系高分子材料を構成する高分子成分の数nで終わるn個の整数番号である。

式（１）では、緩和時間τ_iが小さい高分子成分iの重み付けパラメータＫ_i は大きく計算される。従って、この重み付けパラメータＫ_i を用いることにより、混合系高分子材料の全体の緩和弾性率Ｇ_all（ｔ）の計算にあたって、各高分子成分の緩和時間τ_iが緩和弾性率Ｇ_all（ｔ）に及ぼす影響を考慮することが可能となる。

そして、発明者らは、ついに、混合系高分子材料の全体の緩和弾性率Ｇ_all（ｔ）が、重み付けパラメータＫ_i 及び各高分子成分の体積分率φ_iによる重み付けを考慮した上で、下記式（２）によって近似的に計算されうることを見出した。
Ｇ_all（ｔ）＝Ｋ₁×φ₁×Ｇ₁（ｔ）＋Ｋ₂×φ₂×Ｇ₂（ｔ）＋…＋Ｋ_i×φ_i×Ｇ_i（ｔ）
＝Σ（Ｋ_i×φ_i×Ｇ_i（ｔ）） （２）
ただし、
Ｇ_all（ｔ） ：混合系高分子材料の全体の緩和弾性率
φ_i ：第i番目の高分子成分の混合系高分子材料全体に対する
体積分率
Ｇ_i（ｔ） ：第i番目の各高分子成分の単体での緩和弾性率

ここで、各高分子成分の単体での緩和弾性率Ｇ_i（ｔ）は、その計算に膨大な計算時間を必要とするが、各高分子成分毎に固有のものであり、各高分子成分の混合比率に依存しない。従って、最初に各高分子成分の単体での緩和弾性率Ｇ_i（ｔ）を算出しておけば、各高分子成分の混合比率を変更する場合にあっても、既に算出している緩和弾性率Ｇ_i（ｔ）をそのまま用いることができる。緩和時間τ_i及び重み付けパラメータＫ_i についても同様である。すなわち、各高分子成分の体積分率φ_iを変更しつつ、最初に算出した緩和弾性率Ｇ_i（ｔ）、緩和時間τ_i及び重み付けパラメータＫ_i を上記式（２）に代入することにより、混合比率を変更した混合系高分子材料の全体の緩和弾性率Ｇ_all（ｔ）が容易かつ短時間で計算可能となる。

図６は、本実施形態のシミュレーション方法の処理手順の一例を示すフローチャートである。本シミュレーション方法は、第１工程Ｓ１乃至第５工程Ｓ５を含んでいる。第１工程Ｓ１では、各高分子成分毎に、単体での緩和弾性率Ｇ_i（ｔ）が計算される。第２工程Ｓ２では、各高分子成分毎に、単位分子量における緩和時間τ_iが計算される。第３工程Ｓ３では、第２工程Ｓ２によって計算された緩和時間τ_iを用いて、各高分子成分毎に緩和弾性率Ｇ_i（ｔ）の重み付けパラメータＫ_i が計算される。第４工程Ｓ４では、各高分子成分毎に、混合系高分子材料の全体に対する体積分率φ_iが設定される。そして、第５工程Ｓ５では、第１工程Ｓ１によって計算された緩和弾性率Ｇ_i（ｔ）、第３工程によって計算された重み付けパラメータＫ_i 及び第４工程Ｓ４によって設定された前記体積分率φ_iを用いて、混合系高分子材料の全体の緩和弾性率Ｇ_all（ｔ）が計算される。

図７は、第１工程Ｓ１で各高分子成分の緩和弾性率Ｇ_i（ｔ）を計算する処理手順の一例を示している。本実施形態の混合系高分子材料は、第１高分子成分及び第２高分子成分とが混合されてなる。第１高分子成分として、例えば、ポリスチレンの含有量が２０％のＳＢＲが適用されている。また、第２高分子成分として、例えば、ポリスチレンの含有量が２５％のＳＢＲが適用されている。以下、第１高分子成分の緩和弾性率Ｇ_１（ｔ）を計算する処理手順について説明する。

工程Ｓ１１では、コンピュータ１上で、複数の分子鎖モデル２が配置された高分子材料モデル２０が設定され、コンピュータ１に記憶される。

図８は、工程Ｓ１１において設定される第１高分子成分の高分子材料モデル２０の１つのセル２１を示している。高分子材料モデル２０の概要は、図２に示される混合系高分子材料モデル１０と同様である。すなわち、セル２１は、互いに向き合う三対の平面２２、２２を有する立方体として定義されている。各平面２２には、周期境界条件が定義され、一方の平面２２ａと、反対側の平面２２ｂとが連続しているものとして扱われる。高分子材料モデル２０の１つのセル２１内には、第１高分子成分の分子鎖モデル２がランダムに配置されている。従って、高分子材料モデル２０は、緩和計算がなされていない初期の高分子材料モデルである。

図７において、工程Ｓ１２では、分子動力学計算における温度条件が設定され、コンピュータ１に記憶される。さらに工程Ｓ１３では、粒子モデル３（図４参照）のポテンシャルに基づいて、分子鎖モデル２の人為的な初期配置が緩和される。そして工程Ｓ１４では、体積及び温度が一定、又は圧力及び温度が一定の定常状態が計算され、工程Ｓ１５では、分子動力学計算にて高分子材料モデル２０に生ずる応力が計算される。

分子動力学計算では、例えば、所定の時間、分子鎖モデル２が古典力学に従うものとして、ニュートンの運動方程式が適用される。そして、各時刻での粒子モデル３の動きが、単位時間毎に追跡され、高分子材料モデル２０に生ずる応力が計算される。このような構造緩和の計算は、例えば（株）ＪＳＯＬ社製のソフトマテリアル総合シミュレーター（Ｊ−ＯＣＴＡ）に含まれるＣＯＧＮＡＣを用いて処理することができる。

そして、工程Ｓ１６では、工程Ｓ１５にて計算された応力の変化（揺らぎ）に基づいて、緩和弾性率Ｇ_１（ｔ）が計算される。緩和弾性率Ｇ_１（ｔ）は、所定の時間幅ｔにおいて、歪が与えられた粘弾性体の弾性率の変化を示す指標である。緩和弾性率Ｇ_１（ｔ）は、下記式（３）で計算される。

ここで、
Ｖ ：空間の体積
ｋ_Ｂ：ボルツマン定数
Ｔ ：絶対温度
σ_xy ：応力
ｘｙ：任意の直交する２方向
τ ：時刻
ｔ ：時間幅

上記式（３）において、＜σ_xy（ｔ＋τ）×σ_xy（τ）＞は、所定の時間内において、時刻τの応力σ_xyと、時刻（ｔ＋τ）の応力σ_xyとの積を、あらゆる時刻τについて平均（アンサンブル平均）したものである。工程Ｓ１６で計算された緩和弾性率Ｇ_１（ｔ）は、時間幅ｔごとにコンピュータ１に記憶される。

その後、工程Ｓ１７では、予め定められた終了時刻に到達したか否か、が判断される。終了時刻は、緩和現象の発現及び緩和弾性率Ｇ_１（ｔ）の計算時間を考慮して、適宜定められる。

工程Ｓ１７において終了時刻に到達していない場合は（工程Ｓ１７においてＮ）、工程Ｓ１８で時間幅ｔを１つ増加させ、工程Ｓ１５に戻って、高分子材料モデル２０に生ずる応力が計算される。一方、工程Ｓ１７において終了時刻に到達している場合は（工程Ｓ１７においてＹ）、第１高分子成分の緩和弾性率Ｇ_１（ｔ）の計算を終了する。これにより、工程Ｓ１４で高分子材料モデル２０に変形が与えられてから、上記終了時刻までの緩和弾性率Ｇ_１（ｔ）が得られる。第２高分子成分以降の緩和弾性率Ｇ_２（ｔ）、…についても、同様に計算される。計算された緩和弾性率Ｇ_１（ｔ）は、コンピュータ１に記憶され、第５工程Ｓ５で緩和弾性率Ｇ_all（ｔ）の計算に用いられる。

図９は、第１工程Ｓ１（図６参照）においてコンピュータ１によって計算された第１高分子成分の単体での緩和弾性率Ｇ_１（ｔ）及び第２高分子成分の単体での緩和弾性率Ｇ_２（ｔ）が示されている。図９において、第１高分子成分の緩和弾性率Ｇ_１（ｔ）は実線で、第２高分子成分の緩和弾性率Ｇ_２（ｔ）は一点鎖線で示されている。なお、混合される高分子成分が３以上である場合も、同様に示される。

図６に示される第２工程Ｓ２では、各高分子成分毎に、単体での単位分子量における緩和時間τ_iが計算される。

図１０は、全原子モデルの末端間ベクトルの長さの時間相関関数Ｃ_iと時間ｔとの関係を示すグラフである。緩和時間τ_iは、例えば、特開２０１４−２０６９１３号公報に開示されている時間のパラメータと同様に、末端間ベクトルの相関がなくなる時間、すなわち時間相関関数Ｃ_iが１／ｅに減少するのに要する時間（時定数）で定義される。これは、原子論の崩壊に基づくものである。時間相関関数Ｃ_iは、計算系内の個々の分子について末端間ベクトルの時間相関を算出し、その平均値である計算系全体での時間相関を、下記式（４）に示されるexp関数で近似することにより求められる。
Ｃ_i＝Ａ×exp（−ｔ／τ_i）…（４）
ここで、
Ｃ_i：第i番目の高分子成分の全原子モデルの末端間ベクトルの長さの
時間相関関数
Ａ：フィッティングパラメータ
ｔ：時間
τ_i：第i番目の高分子成分の粗視化モデルの時間のパラメータ
緩和時間τ_iは、exp関数で近似された時間相関関数Ｃ_iから容易に計算されうる。こうして求められた緩和時間τ_iは、コンピュータ１に記憶され、第３工程Ｓ３で重み付けパラメータＫ_i の計算に用いられる。

図６に示される第３工程Ｓ３では、第２工程Ｓ２によって計算された緩和時間τ_iを用いて、各高分子成分毎に緩和弾性率Ｇ_i（ｔ）の重み付けパラメータＫ_i が計算される。重み付けパラメータＫ_i は、上記式（１）を用いて計算され、コンピュータ１に記憶され、第５工程Ｓ５で緩和弾性率Ｇ_all（ｔ）の計算に用いられる。

図６に示される第４工程Ｓ４では、各高分子成分毎に、混合系高分子材料の全体に対する体積分率φ_iが設定される。体積分率φ_iの総和Σφ_i＝１となるように、体積分率φ_iが設定される。例えば、第１高分子成分及び第２高分子成分が同一の混合比で混合される場合、φ_１＝φ_２＝０．５となる。設定された体積分率φ_iは、コンピュータ１に入力され、記憶される。

そして、図６に示される第５工程Ｓ５では、上記式（２）に重み付けパラメータＫ_１、Ｋ_２、…Ｋ_i 、体積分率φ_１、φ_２ 、…φ_i 、及び緩和弾性率Ｇ_１（ｔ）、Ｇ_２（ｔ）、…Ｇ_i （ｔ）、が代入され、混合系高分子材料の全体の緩和弾性率Ｇ_all（ｔ）が算出される。第４工程Ｓ４における体積分率φ_iの設定及び第５工程Ｓ５における緩和弾性率Ｇ_all（ｔ）の計算は、簡素な処理によって短時間で完了する。これらに要する時間は、第１工程Ｓ１における緩和弾性率Ｇ_i（ｔ）の計算に要する時間と比較すると無視できるレベルである。

図１１には、図６のシミュレーション方法によって上記式（２）を用いて計算された混合系高分子材料の全体の緩和弾性率Ｇ_all（ｔ）が、実線で示されている。また、図２に示される混合系高分子材料モデル１０から直接的に計算された全体の緩和弾性率Ｇ’_all（ｔ）が破線で示されている。図１１から明らかなように、本実施形態のシミュレーション方法で計算された緩和弾性率Ｇ_all（ｔ）は、直接的に計算された緩和弾性率Ｇ’_all（ｔ）と非常によく一致しており、混合系高分子材料の全体の緩和弾性率Ｇ_all（ｔ）が下記式（１）によって近似的に計算されうる。

さらに図６に示されるように、本シミュレーション方法では、各高分子成分の混合比率を変更する場合にあっても、第１工程Ｓ１で既に算出している緩和弾性率Ｇ_i（ｔ）及び第３工程Ｓ３で既に算出している重み付けパラメータＫ_i をそのまま用いることができる。すなわち、各高分子成分毎の混合系高分子材料の全体に対する体積分率φ_iを変更する場合は（第６工程Ｓ６においてＹ）、第４工程Ｓ４に戻って体積分率φ_iが再設定され、第５工程Ｓ５にて全体の緩和弾性率Ｇ_all（ｔ）が再計算される。

既に述べたように、第４工程Ｓ４における体積分率φ_iの設定及び第５工程Ｓ５における緩和弾性率Ｇ_all（ｔ）の計算は、極めて短時間で完了するため、混合比率を変更した混合系高分子材料の全体の緩和弾性率Ｇ_all（ｔ）は、極めて短時間で効率よく計算することができる。

図１２は、図７における第１工程Ｓ１で計算される第１高分子成分の単体での緩和弾性率Ｇ_１（ｔ）である。第２高分子成分の単体での緩和弾性率Ｇ_２（ｔ）、・・・についても同様であるので、以下、全ての高分子成分の単体での緩和弾性率Ｇ_i（ｔ）を代表して、緩和弾性率Ｇ_１（ｔ）について説明する。

緩和弾性率Ｇ_１（ｔ）は、ガラス領域Ａ１、転移領域Ａ２及びゴム領域Ａ３を含んでいる。ガラス領域Ａ１とは、第１高分子成分がガラス状態にある緩和弾性率Ｇ_１（ｔ）の領域である。転移領域Ａ２とは、第１高分子成分がガラス状態からゴム状態に転移する転移状態にある緩和弾性率Ｇ_１（ｔ）の領域である。ゴム領域Ａ３とは、第１高分子成分がゴム状態にある緩和弾性率Ｇ_１（ｔ）の領域である。図１２に示されるように、ガラス領域Ａ１、転移領域Ａ２、ゴム領域Ａ３の順に時間スケールが大きくなる。

第１高分子成分の状態は、時間の経過に伴い、ガラス状態から転移状態を経てゴム状態に遷移する。図１２中、ガラス状態とは、第１高分子成分に変形を与えた直後に、ｌｏｇＧ_１（ｔ）が大きい値をとる状態である。ガラス状態では、ｌｏｇｔの増加に対して、ｌｏｇＧ_１（ｔ）が略一定の値となる状態を含む。転移状態とは、ガラス状態とゴム状態との中間の状態であって、ｌｏｇｔの増加に対して、ｌｏｇＧ_１（ｔ）が減少する状態を含む。ゴム状態とは、転移状態の後であって、ｌｏｇｔの増加に対して、ｌｏｇＧ_１（ｔ）が略一定の値となる状態を含む。

例えば、第１高分子成分がタイヤに用いられる場合、ガラス領域Ａ１から転移領域Ａ２にわたる領域での第１高分子成分の緩和弾性率Ｇ_１（ｔ）は、タイヤのグリップ性能と関わりがある。一方、転移領域Ａ２からゴム領域Ａ３にわたる領域での第１高分子成分の緩和弾性率Ｇ_１（ｔ）は、タイヤの転がり抵抗と関わりがある。従って、ガラス領域Ａ１からゴム領域Ａ３にわたって緩和弾性率Ｇ_１（ｔ）を計算することにより、グリップ性能及び転がり抵抗性能に優れたタイヤを開発できる。

本実施形態のシミュレーション方法において、第１高分子成分がガラス状態にあるガラス領域Ａ１での緩和弾性率Ｇ_１（ｔ）は、例えば、全原子モデルでモデル化された分子鎖モデル２を用いた分子動力学によって計算されるのが望ましい。

全原子モデルとは、図４に示されるように、全ての炭素原子及び水素原子がモデル化された分子鎖モデル２である。全原子モデルの分子鎖モデル２を用いて取得された緩和弾性率Ｇ（ｔ）は、小さな時間スケールで計算される反面、実際の高分子材料の緩和弾性率に、高い精度で近似させることができる。従って、全原子モデルの分子鎖モデル２を用いてガラス領域Ａ１での緩和弾性率Ｇ_１（ｔ）を計算することにより、正確な緩和弾性率Ｇ_１（ｔ）を取得できる。

本実施形態の分子鎖モデル２は、全原子モデルとして構成されたものが例示されたが、これに限定されるものではない。分子鎖モデル２は、例えば、炭素原子と、炭素原子に結合した水素原子とを一体化して、一つの粒子モデル（図示省略）として扱うユナイテッドアトムモデル（ united atom model ）として構成されてもよい。ユナイテッドアトムモデルでは、全原子モデルと比較して、高分子鎖のモノマーの配置や、シス構造又はトランス構造を維持しつつ、水素原子を省略することができるため、緩和弾性率Ｇ_１（ｔ）の計算時間を短縮することができる。

第１工程Ｓ１において、全原子モデル又はユナイテッドアトムモデルでモデル化された分子鎖モデルを用いた分子動力学によって緩和弾性率Ｇ_１（ｔ）を計算する際には、５００Ｋ以下の温度条件の下で実行されるのが望ましい。温度が５００Ｋを超える場合、分子鎖モデル２の構造が変化したり（例えば、ブタジエンのｃｉｓ−がｔｒａｎｓ−に変化する）、計算落ちを招くおそれがある。

第１工程Ｓ１において、全原子モデル又はユナイテッドアトムモデルでモデル化された分子鎖モデル２を用いた分子動力学によって緩和弾性率Ｇ_１（ｔ）を計算する際には、分子鎖モデル２の重合度Ｍｎは、５０以下であるのが望ましい。緩和現象が発現するまでの時間は、分子鎖モデル２の重合度Ｍｎに依存する。重合度Ｍｎが５０を超える場合、緩和現象が発現するまでの時間が長くなり、緩和弾性率Ｇ_１（ｔ）まで計算できなくなるおそれがある。

本実施形態のシミュレーション方法において、第１高分子成分が転移状態にある転移領域Ａ２での緩和弾性率Ｇ_１（ｔ）は、例えば、粗視化モデルでモデル化された分子鎖モデル２を用いた分子動力学によって計算されるのが望ましい。

粗視化モデルとは、図３において、モノマーをなす単位構造６が置換された粗視化粒子モデル（図示省略）と、粗視化粒子モデル間を連結する結合鎖モデル（図示省略）とを含む粗視化モデル（図示省略）として構成された分子鎖モデル２である。既に述べたように、転移領域Ａ２での緩和弾性率Ｇ_１（ｔ）は、ガラス領域Ａ１での緩和弾性率Ｇ_１（ｔ）に比べて、大きな時間スケールに分布している。粗視化モデルは、全原子モデルやユナイテッドアトムモデルに比べて簡素化されているので、緩和弾性率Ｇ_１（ｔ）の計算時間が短縮される。従って、粗視化モデルは、ガラス領域Ａ１よりも大きな時間スケールで緩和弾性率Ｇ_１（ｔ）を計算することが求められる転移領域Ａ２での計算に適する。

本実施形態のシミュレーション方法において、第１高分子成分がゴム状態にあるゴム領域Ａ３での緩和弾性率Ｇ_１（ｔ）は、例えば、管模型でモデル化された分子鎖モデル２を用いて計算されるのが望ましい。

管模型とは、からみ合った分子鎖が運動できる領域を管に見立てて、管の運動を解くモデルである。管模型では、分子鎖のからみ合い点は、管の折れ曲がり点として扱われる。からみ合い点は、緩和弾性率Ｇ_１（ｔ）の計算中、他の分子鎖とのからみ合いによる力を受けて、セル２１内を運動する。

管模型では、分子鎖Ｍｃの詳細な原子構造には立ち入らず、からみ合い点間を繋ぐ伸縮可能な管によって一つの分子鎖Ｍｃが表現される。このため、管模型は、上述した粗視化モデルと比較しても、より一層簡素化された分子鎖モデル２であり、緩和弾性率Ｇ_１（ｔ）の計算時間が大幅に短縮される。従って、管模型は、転移領域Ａ２よりもさらに大きな時間スケールで緩和弾性率Ｇ_１（ｔ）を計算することが求められるゴム領域Ａ３での計算に適する。

全原子モデルを用いて計算された緩和弾性率Ｇ（ｔ）と、粗視化モデルを用いて計算された緩和弾性率Ｇ（ｔ）と、管模型とを用いて計算された緩和弾性率Ｇ（ｔ）と組み合わせる方法としては、例えば、「第６２回高分子討論会予稿集1J13」に基づいて実施されうる。

ガラス領域Ａ１と転移領域Ａ２との境界及び転移領域Ａ２とゴム領域Ａ３との境界は、例えば、コンピュータ１の計算能力に応じて適宜設定されうる。

図８において、セル２１に配置される分子鎖モデル２の本数については、適宜設定されうる。分子鎖モデル２の本数が少ないと、分子動力学の計算において、一方の平面２２ａから出て行き、かつ、反対側の平面２２ｂから入ってきた分子鎖モデル２の一端側と、この分子鎖モデル２の他端側とがからまって、計算落ちを招くおそれがある。逆に、分子鎖モデル２の本数が多くても、運動方程式の質点として取り扱われる粒子モデルが増大し、多くの計算時間を要するおそれがある。このような観点より、セル２１に配置される分子鎖モデル２の本数は、好ましくは８本以上であり、また、好ましくは、１００本以下である。

緩和時間τ_iを計算する第２工程Ｓ２では、例えば、全原子モデル又はユナイテッドアトムモデルでモデル化された分子鎖モデル２を用いた分子動力学によって計算されるのが望ましい。既に述べたように、全原子モデルの分子鎖モデル２を用いて緩和時間τ_iを計算することにより、正確な緩和時間τ_iを緩和弾性率Ｇ_１（ｔ）を取得できる。ユナイテッドアトムモデルでは、全原子モデルと比較して、緩和時間τ_iの計算時間を短縮することができる。

そして、この第２工程Ｓ２で全原子モデル又はユナイテッドアトムモデルでモデル化された分子鎖モデル２を用いた分子動力学によって緩和時間τ_iを計算する際には、分子鎖モデル２の重合度Ｍｎは、５０以下であるのが望ましい。重合度Ｍｎが５０を超える場合、緩和現象が発現するまでの時間が長くなり、緩和弾性率Ｇ_１（ｔ）まで計算できなくなるおそれがある。

以上、本発明の特に好ましい実施形態について詳述したが、本発明は図示の実施形態に限定されることなく、種々の態様に変形して実施しうる。

第１高分子成分として、ポリスチレンの含有量が２０％のＳＢＲの分子鎖モデル２と、第２高分子成分として、ポリスチレンの含有量が２５％のＳＢＲの分子鎖モデル２とが設定された。表１の仕様に基づき、第１高分子成分及び第２高分子成分が３種類の混合比率で混合された混合系高分子材料が設定された。各混合系高分子材料の全体の緩和弾性率Ｇ_all（ｔ）が、図６に示される処理手順に従って計算され、その計算時間が測定された。比較例として、図２に示される混合系高分子材料モデル１０から直接的に全体の緩和弾性率Ｇ’_all（ｔ）が計算され、その計算時間が測定された。

表１から明らかなように、実施例のシミュレーション方法は、比較例と比べて、３種類以上の混合比率で混合された混合系高分子材料が設定される場合に、各混合系高分子材料の全体の緩和弾性率Ｇ_all（ｔ）が、短時間で計算できることが確認された。また、図１１に示されるように、実施例で計算された緩和弾性率Ｇ_all（ｔ）は、比較例で計算された緩和弾性率Ｇ’_all（ｔ）とよく一致していることが確認された。

２ 分子鎖モデル
１１ 高分子材料モデル
２１ セル
Ａ１ ガラス領域
Ａ２ 転移領域
Ａ３ ゴム領域
Ｇ（ｔ） 緩和弾性率

Claims (8)

1. 複数種類の高分子成分が混合されてなる混合系高分子材料の全体の緩和弾性率Ｇ_all（ｔ）を予測する高分子材料のシミュレーション方法であって、
   各高分子成分毎に、単体での緩和弾性率Ｇ_i（ｔ）を計算する第１工程と、
   前記各高分子成分毎に、単体での単位分子量における緩和時間τ_iを計算する第２工程と、
   前記第２工程によって計算された前記緩和時間τ_iを用いて、下記式（１）で、前記各高分子成分毎に前記緩和弾性率Ｇ_i（ｔ）の重み付けパラメータＫ_i を計算する第３工程と、
   前記各高分子成分毎に、前記混合系高分子材料の全体に対する体積分率φ_iを設定する第４工程と、
   前記第１工程によって計算された前記緩和弾性率Ｇ_i（ｔ）、前記第３工程によって計算された前記重み付けパラメータＫ_i 及び前記第４工程によって設定された前記体積分率φ_iを用いて、下記式（２）で、前記混合系高分子材料の全体の緩和弾性率Ｇ_all（ｔ）を計算する第５工程とを含むことを特徴とする高分子材料のシミュレーション方法。
   Ｋ_i ＝（１／τ_i）／Σ（１／τ_i ） （１）
   ただし、
   i ：各高分子成分に割り当てられた整数番号
   Ｋ_i ：第i番目の高分子成分の重み付けパラメータ
   τ_i ：第i番目の高分子成分の単体での緩和時間
   Ｇ_all（ｔ）＝Σ（Ｋ_i×φ_i×Ｇ_i（ｔ）） （２）
   ただし、
   Ｇ_all（ｔ） ：混合系高分子材料の全体の緩和弾性率
   φ_i ：第i番目の高分子成分の混合系高分子材料全体に対する
   体積分率
   Ｇ_i（ｔ） ：第i番目の各高分子成分の単体での緩和弾性率
2. 前記各高分子成分の単体での前記緩和弾性率Ｇ_i（ｔ）は、前記各高分子成分がガラス状態にあるガラス領域での緩和弾性率Ｇ_i（ｔ）を含み、
   前記第１工程において、少なくとも前記ガラス領域での前記緩和弾性率Ｇ_i（ｔ）は、全原子モデル又はユナイテッドアトムモデルでモデル化された分子鎖モデルを用いた分子動力学によって計算される請求項１記載の高分子材料のシミュレーション方法。
3. 前記第１工程は、５００Ｋ以下の温度条件の下で行なわれる請求項２記載の高分子材料のシミュレーション方法。
4. 前記各高分子成分の単体での緩和弾性率Ｇ_i（ｔ）は、前記各高分子成分がゴム状態にあるゴム領域での緩和弾性率Ｇ_i（ｔ）を含み、
   前記第１工程において、少なくとも前記ゴム領域での前記緩和弾性率Ｇ_i（ｔ）は、管模型でモデル化された分子鎖モデルを用いて計算される請求項２記載の高分子材料のシミュレーション方法。
5. 前記各高分子成分の単体での緩和弾性率Ｇ_i（ｔ）は、前記ガラス領域と前記ゴム領域との中間の領域である転移領域での緩和弾性率Ｇ_i（ｔ）を含み、
   前記第１工程において、少なくとも前記転移領域の前記緩和弾性率Ｇ_i（ｔ）は、粗視化モデルでモデル化された分子鎖モデルを用いた分子動力学によって計算される請求項４記載の高分子材料のシミュレーション方法。
6. 前記第１工程では、１つのセル内に８本以上の前記分子鎖モデルを配置した材料モデルを用いて周期境界条件の下で、分子動力学計算を行なう工程を含む請求項２乃至５のいずれかに記載の高分子材料のシミュレーション方法。
7. 前記第２工程では、前記緩和時間τ_iは、全原子モデル又はユナイテッドアトムモデルでモデル化された分子鎖モデルを用いた分子動力学によって計算される請求項１乃至６のいずれかに記載の高分子材料のシミュレーション方法。
8. 前記分子鎖モデルの重合度は、５０以下である請求項２又は７記載の高分子材料のシミュレーション方法。