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

Abstract

【課題】フィラーの周囲に形成される界面層の厚さを計算する。
【解決手段】コンピュータ１を用いて、高分子材料とフィラーとの反応を解析するための方法である。このシミュレーション方法では、コンピュータ１が、高分子材料の界面層１４ａの厚さＷ１を計算する界面厚計算工程Ｓ８を含む。粗視化モデル２のビーズ３、３間には、該ビーズ３、３間の距離が小さくなるほど無限に大きくなる斥力を含むポテンシャルが定義される。界面厚計算工程Ｓ８は、粗視化モデル２が配置された空間７を、フィラーモデル１１の外面１１ｏに沿った境界面１７で複数の領域１６に区分する工程Ｓ８１と、各領域１６の緩和弾性率を計算する緩和弾性率計算工程Ｓ８２と、各領域１６の緩和弾性率に基づいて、界面層１４ａの厚さＷ１を求める特定工程Ｓ８３とを含む。
【選択図】図３

Description

本発明は、フィラーと高分子材料とが共存する系において、フィラーの周囲に形成される界面層の厚さを計算することができる高分子材料のシミュレーション方法に関する。

近年、ゴム等の高分子材料の設計・開発のために、コンピュータを用いた高分子材料のシミュレーション方法が種々提案されている（例えば、下記特許文献１参照）。この種のシミュレーション方法では、高分子材料の構造や、フィラーの配合率等の諸条件を、計算に織り込むことができる。従って、このシミュレーション方法では、実際に高分子材料を試作することなく、その物性値を計算することができる。

特開２０１２−２３８１６８号公報

ところで、フィラーと高分子材料との共存する系においては、フィラーの周囲に界面層が形成されることが判明している。この界面層は、高分子材料本来の部分（以下、バルク部分という。）とは異なる力学的特性を示す層として知られており、ガラス層などと称されることもある。

従来のシミュレーション方法でも、高分子材料の物性値を計算するに先立ち、界面層の厚さが定義されていた。しかしながら、界面層の厚さは、例えば、高分子材料及びフィラーを用いた実験結果に基づいて定義されていた。このため、従来のシミュレーション方法では、多くの時間や費用が必要になるという問題点があった。従って、界面層の厚さを計算することができる高分子材料のシミュレーション方法が強く求められていた。

本発明は、以上のような実状に鑑み案出されたもので、フィラーと高分子材料とが共存する系の界面層の厚さを計算することができる高分子材料のシミュレーション方法を提供することを主たる目的としている。

本発明のうち請求項１記載の発明は、コンピュータを用いて、高分子材料とフィラーとの反応を解析するための方法であって、前記コンピュータに、前記高分子材料の高分子鎖を、モノマー又はモノマーの一部分をなす構造単位をビーズで置換した粗視化モデルを設定する工程、前記コンピュータに、前記フィラーの外面を少なくとも含むフィラーモデルを設定する工程、前記コンピュータが、予め定められた空間内で前記粗視化モデルと前記フィラーモデルとを用いて分子動力学計算を行うシミュレーション工程、及び前記コンピュータが、前記シミュレーション工程の結果から、前記フィラーモデルの周囲に形成されかつ前記高分子材料のバルク部分とは異なる力学的特性を示す前記高分子材料の界面層の厚さを計算する界面厚計算工程を含み、前記粗視化モデルの前記ビーズ間には、該ビーズ間の距離が小さくなるほど無限に大きくなる斥力を含むポテンシャルが定義され、前記界面厚計算工程は、前記粗視化モデルが配置された前記空間を、前記フィラーモデルの外面に沿った境界面で複数の領域に区分する工程と、前記各領域の緩和弾性率を計算する緩和弾性率計算工程と、前記各領域の緩和弾性率に基づいて、前記界面層の厚さを求める特定工程とを含むことを特徴とする。

また、請求項２記載の発明は、前記特定工程は、隣り合う前記領域間での前記緩和弾性率が予め定めた範囲内になった領域以外を前記界面層として決定し、その厚さを計算する請求項１記載の高分子材料のシミュレーション方法である。

また、請求項３記載の発明は、前記空間は、その外周を囲む一つの平面を含み、前記フィラーモデルは、前記平面でモデル化され、前記境界面は、前記フィラーモデルと平行である請求項１又は２に記載の高分子材料のシミュレーション方法である。

また、請求項４記載の発明は、前記空間は、その外周を囲み、かつ互いに向き合う一対の平面を含み、前記フィラーモデルは、前記一対の平面でモデル化され、前記境界面は、一対の前記フィラーモデルと平行である請求項１又は２に記載の高分子材料のシミュレーション方法である。

また、請求項５記載の発明は、前記粗視化モデルは、一対の前記フィラーモデルの間に配置される請求項４に記載の高分子材料のシミュレーション方法である。

また、請求項６記載の発明は、前記一対の平面間の距離は、前記粗視化モデルの慣性半径の２倍以上である請求項４又は５に記載の高分子材料のシミュレーション方法である。

また、請求項７記載の発明は、前記境界面間の各距離は、前記粗視化モデルの慣性半径に基づいて設定される請求項４乃至６のいずれかに記載の高分子材料のシミュレーション方法である。

また、請求項８記載の発明は、前記緩和弾性率計算工程は、前記各フィラーモデルから等距離にある一対の前記領域を合わせて、前記緩和弾性率を計算する請求項４乃至７のいずれかに記載の高分子材料のシミュレーション方法である。

また、請求項９記載の発明は、前記界面厚計算工程は、前記緩和弾性率計算工程に先立ち、前記各粗視化モデルの前記ビーズ間の距離を同一に揃える工程を含む請求項１乃至８のいずれかに記載の高分子材料のシミュレーション方法である。

また、請求項１０記載の発明は、前記斥力を含むポテンシャルは、下記式（１）で定義される請求項１乃至９のいずれかに記載の高分子材料のシミュレーション方法である。

ここで、各定数及び変数は次のとおりである。
ｒ_ij：ビーズ間の距離
ε：ビーズ間に定義される斥力ポテンシャルＲの強度に関する定数
σ：ビーズ間に定義される斥力ポテンシャルＲが作用する距離に関する定数（分子動力学の分野では、ＬＪ球の直径と呼ばれる）

本発明の高分子材料のシミュレーション方法は、コンピュータが、予め定められた空間に配置された粗視化モデルとフィラーモデルとを用いて分子動力学計算を行うシミュレーション工程、及び高分子材料の界面層の厚さを計算する界面厚計算工程を含む。

界面厚計算工程は、粗視化モデルが配置された空間を複数の領域に区分する工程と、各領域の緩和弾性率を計算する緩和弾性率計算工程と、緩和弾性率に基づいて、界面層の厚さを求める特定工程とを含む。

一般に、高分子材料のバルク部分では、緩和弾性率は、高分子材料に固有の一定の値を示す。従って、本発明では、粗視化モデルを用いたシミュレーションにおいて、空間が区分された各領域の緩和弾性率を計算することにより、界面層を特定し、さらにはその厚さを確実に計算することができる。

さらに、本発明では、粗視化モデルのビーズ間に、該ビーズ間の距離が小さくなるほど無限に大きくなる斥力を含むポテンシャルが定義される。このようなポテンシャルは、高分子材料の分子運動に基づくものである。従って、本発明では、実際の高分子材料の分子運動に近似させて、粗視化モデルの初期配置を精度よく緩和することができるため、界面層の厚さを精度よく計算することができる。

本実施形態のシミュレーション方法を実行するコンピュータの斜視図である。 ポリブタジエンの構造式である。 本実施形態のシミュレーション方法の処理手順の一例を示すフローチャートである。 粗視化モデルを示す概念図である。 粗視化モデルのポテンシャルを説明する概念図である。 空間及びフィラーモデルの概念図である。 複数の粗視化モデルが配置された空間を示す概念図である。 フィラーモデルと粗視化モデルとの間の相互ポテンシャルを説明する概念図である。 平衡状態の粗視化モデルが配置された空間を示す概念図である。 界面層とバルク部分とを含む高分子材料の線図である。 本実施形態の界面厚計算工程の処理手順の一例を示すフローチャートである。 界面厚計算工程を説明する概念図である。 第一領域〜第五領域の各緩和弾性率Ｇ（ｔ）の自然対数をとった値と、時間幅ｔの自然対数をとった値との関係を示すグラフである。 他の実施形態のフィラーモデルを示す概念図である。

以下、本発明の実施の一形態が図面に基づき説明される。
本実施形態の高分子材料のシミュレーション方法（以下、単に「シミュレーション方法」ということがある）は、コンピュータを用いて、高分子材料とフィラーとの反応を解析するためのものである。

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

フィラーとしては、例えば、カーボンブラック、シリカ又はアルミナ等が含まれる。また、高分子材料としては、例えば、ゴム、樹脂又はエラストマー等が含まれる。本実施形態では、高分子材料として、図２に示されるように、cis-１，４ポリブタジエン（以下、単に「ポリブタジエン」ということがある）が例示される。このポリブタジエンを構成する高分子鎖は、メチレン基（−ＣＨ_２−）とメチン基（−ＣＨ−）とからなるモノマー｛−［ＣＨ_２−ＣＨ＝ＣＨ−ＣＨ_２］−｝が、重合度ｎで連結されて構成されている。なお、高分子材料には、ポリブタジエン以外の高分子材料が用いられてもよい。

図３には、本実施形態のシミュレーション方法の具体的な処理手順が示されている。このシミュレーション方法では、先ず、コンピュータ１に、高分子材料の高分子鎖を複数のビーズでモデル化した粗視化モデルが設定される（工程Ｓ１）。図２及び図４に示されるように、高分子材料の高分子鎖がポリブタジエンである場合には、例えば１．５５個分のモノマーを構造単位４として、該構造単位４を１個のビーズ３に置換される。これにより、複数（例えば、１０〜５０００個）のビーズ３でモデル化された粗視化モデル２が設定される。

なお、１．５５個分のモノマーを構造単位４としたのは、論文（L,J.Fetters ,D.J.Lohse and R.H.Colby 著、「Chain Dimension and Entanglement Spacings 」Physical Properties of Polymers Handbook Second Edition P448」）と、論文（ Kurt Kremer & Gary S. Grest 著 「Dynamics of entangled linear polymer melts: A molecular-dynamics simulation」、J. Chem Phys. vol.92, No.8, 15 April 1990）の記載に基づき求めたものである。また、高分子鎖がポリブタジエン以外の場合でも、例えば、上記論文に基づいて、構造単位を設定することができる。

ビーズ３は、分子動力学計算において、運動方程式の質点として取り扱われる。即ち、ビーズ３には、質量、体積、直径、電荷又は初期座標などのパラメータが定義される。

また、粗視化モデル２は、ビーズ３、３間を結合するボンド５が設定される。このボンド５には、結合ポテンシャルＰ１が設定される。本実施形態の結合ポテンシャルＰ１は、下記式（１）で定義される斥力を含むポテンシャル（本実施形態では、斥力のみを作用させるポテンシャル（以下、「斥力ポテンシャル」ということがある。））Ｒと、下記式（２）で定義される引力ポテンシャルＧとの和（Ｒ＋Ｇ）で定義される。

ここで、各定数及び変数は、Lennard-Jonesポテンシャルのパラメータであり、次のとおりである。
ｒ_ij：ビーズ間の距離
ｋ：ビーズ間のばね定数
ε：ビーズ間に定義される斥力ポテンシャルＲの強度に関する定数
σ：ビーズ間に定義される斥力ポテンシャルＲが作用する距離に関する定数（分子動力学の分野では、ＬＪ球の直径と呼ばれる）
Ｒ₀：伸びきり長
なお、距離ｒ_ij及び伸びきり長Ｒ₀は、各ビーズ３、３の中心３ｃ、３ｃ間の距離として定義される。

上記式（１）において、ビーズ３、３間の距離ｒ_ijが小さいほど、斥力ポテンシャルＲが大きくなる。さらに、上記式（１）では、ビーズ３、３間の距離ｒ_ijが小さくなるほど、斥力ポテンシャルＲが無限に大きくなる。このような斥力ポテンシャルＲは、高分子材料の分子運動に基づくものである。

一方、上記式（２）において、距離ｒ_ijが大きいほど、引力ポテンシャルＧが大きくなる。従って、結合ポテンシャルＰ１は、距離ｒ_ijを、斥力ポテンシャルＲと引力ポテンシャルＧとが互いに釣り合う位置に戻そうとする復元力が定義される。

また、上記式（２）では、距離ｒ_ijが伸びきり長Ｒ₀以上となる場合には、引力ポテンシャルＧが∞に設定される。従って、結合ポテンシャルＰ１は、伸びきり長Ｒ₀以上の距離ｒ_ijを許容しない。このような結合ポテンシャルＰ１は、粗視化モデル２を、高分子材料の分子運動に近似させることができる。

なお、斥力ポテンシャルＲ及び引力ポテンシャルＧの各定数及び変数の値としては、適宜設定することができる。本実施形態では、論文（Kurt Kremer & Gary S. Grest 著 「Dynamics of entangled linear polymer melts: A molecular-dynamics simulation」、J. Chem Phys. vol.92, No.8, 15 April 1990）に基づいて、次の値が設定される。
ばね定数ｋ：３０
伸びきり長Ｒ₀：１．５
定数ε：１．０
定数σ：１．０

また、ボンド５には、上記のような結合ポテンシャルＰ１のほか、下記式（３）で定義される結合ポテンシャルＰ２が設定されてもよい。

ここで、各定数及び変数は、次に示すものを除いて、上記式（１）、（２）と同一である。
ｒ：ビーズ間の距離
ｒ₀：平衡長

上記式（３）で定義される結合ポテンシャルＰ２は、Harmonic型である。Harmonic型とは、いわゆる線形バネが定義され、平衡長ｒ₀からのずれの大きさに比例した復元力が働くポテンシャルである。なお、平衡長ｒ₀は、各粒子の中心間の距離として定義される。

上記式（３）において、距離ｒと平衡長ｒ₀とが等しい場合には、結合ポテンシャルＰ２が０となる。また、距離ｒが、平衡長ｒ₀よりも大である場合は、該距離ｒが大きくなるほど、ビーズ３、３が近づく方向に、結合ポテンシャルＰ２（引力ポテンシャル）が無限に大きくなる。一方、距離ｒが、平衡長ｒ₀よりも小である場合は、該距離ｒが小さくなるほど、ビーズ３、３が離れる方向に、結合ポテンシャルＰ２（斥力ポテンシャル）が大きくなる。このように、結合ポテンシャルＰ２では、距離ｒを平衡長ｒ₀に戻そうとする復元力が定義される。

また、結合ポテンシャルＰ２は、結合ポテンシャルＰ１とは異なり、距離ｒ_ijの上限値が、平衡長ｒ₀に限定されない。従って、結合ポテンシャルＰ２は、平衡長ｒ₀以上の距離ｒを許容することができる。このような結合ポテンシャルＰ２は、分子動力学計算において、粗視化モデル２、２同士の絡まりによって作用する大きな力を抑制することができるため、計算落ちの発生を防ぐことができる。

上記式（３）のばね定数ｋ及び平衡長ｒ₀の各値は、適宜設定することができる。本実施形態では、結合ポテンシャルＰ１と同様に、上記論文に基づいて、次の値が設定される。
ばね定ｋ：９００
平衡長ｒ₀：０．９

図５に示されるように、隣接する粗視化モデル２、２のビーズ３、３間には、上記式（１）の斥力ポテンシャルＲが定義される。なお、斥力ポテンシャルＲの強度ε及び定数σも、上記式（１）と同一が望ましい。

このような粗視化モデル２は、コンピュータ１で取り扱い可能な数値データであり、コンピュータ１に記憶される。

次に、コンピュータ１に、予め定められた体積をもつ空間が設定される（工程Ｓ２）。図６に示されるように、空間７は、その外周を囲む一つの平面、本実施形態では、互いに向き合う三対の平面８、８を有する立方体として定義される。

また、一対の平面８、８に直交する方向において、該一対の平面８、８間の距離Ｄ１（即ち、１辺の長さＬ１）は、粗視化モデル２（図４）の慣性半径（図示省略）の２倍以上、より好ましくは４倍以上に設定されるのが望ましい。慣性半径は、分子動力学計算において、粗視化モデル２の拡がりを示すパラメータである。このような空間７では、分子動力学計算において、粗視化モデル２の回転運動をスムーズに計算することができる。また、空間７の大きさは、上記論文に基づいて、粒子数密度が０．８５個／σ³程度に設定される。このような空間７は、コンピュータ１で取り扱い可能な数値データであり、コンピュータ１に記憶される。

次に、コンピュータ１に、フィラーの外面を少なくとも含むフィラーモデルが設定される（工程Ｓ３）。フィラーモデル１１は、空間７の少なくとも１枚の平面８、本実施形態では、空間７の上下に配された一対の平面８、８でモデル化される。これにより、一対のフィラーモデル１１、１１は、空間７において、移動不能に固定された外面１１ｏ、１１ｏ（図８に示す）のみが定義される。このようなフィラーモデル１１も、コンピュータ１で取り扱い可能な数値データであり、コンピュータ１に記憶される。なお、フィラーモデル１１、１１の各外面１１ｏ、１１ｏは、粗視化モデル２（図４に示す）が通過不能に定義される。一方、フィラーモデル１１が定義されない平面８には、周期境界条件が定義されている。

次に、図７に示されるように、空間７内に、複数の粗視化モデル２が配置される（工程Ｓ４）。これにより、空間７は、解析対象の高分子材料の微小構造部分として構成される。また、本実施形態の粗視化モデル２は、空間７内において、一対のフィラーモデル１１、１１の間に配置される。なお、粗視化モデル２の本数は、例えば、２００〜４００本程度が望ましい。

次に、コンピュータ１に、粗視化モデル２とフィラーモデル１１との間の相互ポテンシャルＴが設定される（工程Ｓ５）。本実施形態の工程Ｓ５では、図８に示されるように、粗視化モデル２のビーズ３とフィラーモデル１１の外面１１ｏとの間に、引力及び斥力が作用する相互ポテンシャルＴが設定される。相互ポテンシャルＴは、例えば、下記式（４）で定義される。

ここで、各定数及び変数は次のとおりである。
ρ_wall：相互ポテンシャルＴの壁面の密度に関する定数
ε_wall：相互ポテンシャルＴの強度に対応する定数
σ_wall：空間の平面（フィラーモデル）に直交する方向の斥力長さに関する定数
ｒ：フィラーモデルとビーズとの間及びフィラーモデルと変性基粒子との間の距離
ｒ_c：カットオフ距離
なお、距離ｒ及びカットオフ距離ｒ_cは、フィラーモデル１１の平面８と、粗視化モデル２のビーズ３の中心３ｃ（図４に示す）との間の最短距離で定義される。

上記式（４）で定義される相互ポテンシャルＴは、上記式（１）で定義される斥力ポテンシャルＲを、フィラーモデル１１の平面８で積分したものである。従って、相互ポテンシャルＴは、フィラーモデル１１の外面１１ｏ全体に定義される。

上記式（４）において、距離ｒが、予め定められたカットオフ距離ｒ_c以上になる場合は、相互ポテンシャルＴは作用しない。また、距離ｒが、相互ポテンシャルＴの定数σ_wall×２^1/6未満になると、相互ポテンシャルＴは、粗視化モデル２のビーズ３とフィラーモデル１１との間に斥力のみを作用させることができる。一方、距離ｒが相互ポテンシャルの定数σ_wall×２^1/6を超えると、相互ポテンシャルＴは、粗視化モデル２のビーズ３とフィラーモデル１１との間に引力を作用させることができる。このように、上記式（４）では、粗視化モデル２のビーズ３とフィラーモデル１１との間に、引力又は斥力を作用させることができる。

また、相互ポテンシャルＴの定数ρ_wall、定数σ_wall、強度ε_wall、カットオフ距離ｒ_cの各値としては、適宜設定することができる。本実施形態では、上記論文に基づいて、次の値が設定される。
定数ρ_wall＝定数σ_wall＝１．０
強度ε_wall＝１
カットオフ距離ｒ_c＝２^1/6

次に、コンピュータ１が、フィラーモデル１１と、粗視化モデル２とを用いて分子動力学計算による緩和計算を行う（シミュレーション工程Ｓ６）。

本実施形態の分子動力学計算では、例えば、空間７について所定の時間、粗視化モデル２が古典力学に従うものとして、ニュートンの運動方程式が適用される。そして、各時刻でのビーズ３の動きが、単位時間毎に追跡される。また、空間７内のビーズ３の数、体積及び温度は、夫々一定に保たれる。

本実施形態のシミュレーション工程Ｓ６では、フィラーモデル１１が固定された状態で、粗視化モデル２のみを対象に緩和計算が行われる。従って、本実施形態では、フィラーモデル１１及び粗視化モデル２の双方を対象に緩和計算が行なわれていた従来の方法に比べて、計算時間を短縮しうる。

また、フィラーモデル１１は、外面１１ｏの全体において、該外面１１ｏと直角に交わる方向に、各相互ポテンシャルＴ（図８に示す）を作用させることができる。これにより、フィラーモデル１１は、粗視化モデル２に対して、同一方向の相互ポテンシャルＴを平面８全体で作用させることができるため、粗視化モデル２の緩和計算を効率的に行うことができる。

さらに、本実施形態では、粗視化モデル２のビーズ３、３間に、該ビーズ３、３間の距離ｒ_ijが小さくなるほど無限に大きくなる斥力ポテンシャルが定義されるため、実際の高分子材料の分子運動に近似させて、粗視化モデル２の初期配置を精度よく緩和することができる。

次に、コンピュータ１が、粗視化モデル２の初期配置を十分に緩和できたか否かを判断する（工程Ｓ７）。この工程Ｓ７では、粗視化モデル２の初期配置を十分に緩和できたと判断された場合、次の界面厚計算工程Ｓ８が実施される。一方、粗視化モデル２の初期配置を十分に緩和できていないと判断された場合は、単位ステップを進めて（工程Ｓ９）、シミュレーション工程Ｓ６（分子動力学計算）が再度実施される。これにより、シミュレーション工程Ｓ６では、図９に示されるように、粗視化モデル２の平衡状態（構造が緩和した状態）を確実に計算することができる。

次に、コンピュータ１が、シミュレーション工程Ｓ６の結果から、高分子材料の界面層の厚さを計算する（界面厚計算工程Ｓ８）。図１０に示されるように、界面層１４ａは、フィラー１３と高分子材料１４との共存する系において、フィラー１３の周囲に形成されることが知られている。また、界面層１４ａでは、高分子材料１４の本来の部分であるバルク部分１４ｂとは異なる力学的特性が示される。例えば、高分子材料の微小構造部分として構成された空間７（図９に示す）を用いて、高分子材料の物性値を計算する場合には、界面層１４ａの厚さＷ１を予め定義しておくことが重要である。

図１１には、本実施形態の界面厚計算工程Ｓ８の具体的な処理手順が示されている。本実施形態の界面厚計算工程Ｓ８では、先ず、コンピュータ１が、粗視化モデル２が配置された空間７を複数の領域に区分する（工程Ｓ８１）。

図１２に示されるように、本実施形態の領域１６は、一方のフィラーモデル１１Ａと他方のフィラーモデル１１Ｂとの間において、各フィラーモデル１１Ａ、１１Ｂの外面１１ｏに沿った複数（例えば、Ｎ−１個）の境界面１７で区分される。これにより、空間７には、Ｎ個分の領域１６が定義される。本実施形態の領域１６は、一方のフィラーモデル１１Ａ側の第一領域１６Ａ〜他方のフィラーモデル１１Ｂ側の第Ｎ領域１６Ｎを含んでいる。
なお、図１２では、粗視化モデル２を省略して表示している。

また、各境界面１７は、一対のフィラーモデル１１Ａ、１１Ｂと平行に設定される。さらに、各領域１６Ａ〜１６Ｎの隣り合う境界面１７、１７間において、該境界面１７に直交する方向の各距離Ｌ２は、夫々同一に設定されている。これにより、各領域１６Ａ〜１６Ｎの体積は、それぞれ同一に設定される。これらの領域１６Ａ〜１６Ｎは、いずれも数値データであり、コンピュータ１に記憶される。

次に、コンピュータ１が、各領域１６Ａ〜１６Ｎの緩和弾性率を計算する（緩和弾性率計算工程Ｓ８２）。緩和弾性率Ｇ（ｔ）は、所定の時間幅ｔにおいて、歪が与えられた粘弾性体の弾性率の変化を示す指標である。各領域１６Ａ〜１６Ｎの緩和弾性率Ｇ（ｔ）は、下記式（５）で計算される。

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

上記式（５）において、＜σ_xy（ｔ＋τ）×σ_xy（τ）＞は、所定の時間内において、時刻τの応力σ_xyと、時刻（ｔ＋τ）の応力σ_xyとの積を、あらゆる時刻τについて平均（アンサンブル平均）したものである。また、緩和弾性率Ｇ（ｔ）が計算される時間としては、シミュレーションにおける単位時間τ’を用いて、１０００τ’〜１０００万τ’が望ましい。この緩和弾性率Ｇ（ｔ）は、領域１６Ａ〜１６Ｎ毎に、コンピュータ１に記憶される。

次に、コンピュータ１が、各領域１６Ａ〜１６Ｎの緩和弾性率に基づいて、界面層１４ａの厚さＷ１を求める（特定工程Ｓ８３）。界面層１４ａでは、バルク部分１４ｂとは異なる力学的特性（例えば、緩和弾性率Ｇ（ｔ）を含む）が示される。このため、本実施形態の特定工程Ｓ８３は、隣り合う領域１６、１６間での緩和弾性率Ｇ（ｔ）が予め定めた範囲内になった領域１６以外を、界面層１４ａとして決定する。

図１３には、一方のフィラーモデル１１Ａ側において、第一領域１６Ａ〜第五領域１６Ｅの各緩和弾性率Ｇ（ｔ）の自然対数をとった値と、時間幅ｔの自然対数をとった値との関係を示すグラフが示される。このグラフでは、第四領域１６Ｄの緩和弾性率Ｇ（ｔ）の自然対数をとった値と、第五領域１６Ｅの緩和弾性率Ｇ（ｔ）の自然対数をとった値とが略同一であることを確認することができる。従って、図１２に示されるように、第四領域１６Ｄ以降の各領域は、バルク部分１４ｂであると決定される。一方、第一領域１６Ａ〜第三領域１６Ｃは、界面層１４ａであると決定される。

また、他方のフィラーモデル１１Ｂ側においても同様に、図１３のグラフを用いて、界面層１４ａとバルク部分１４ｂとが決定される。なお、緩和弾性率Ｇ（ｔ）の自然対数をとった値が略同一か否かの判断は、隣り合う領域において、緩和弾性率Ｇ（ｔ）の自然対数をとった値の比が、０．７〜１．３の範囲内であるか否かで判断されるのが望ましい。

そして、図１２に示されるように、界面層１４ａとバルク部分１４ｂとの境界面１７（本実施形態では、第三領域１６Ｃと第四領域１６Ｄとの境界面１７）と、フィラーモデル１１との最短距離Ｌ３が計算される。この最短距離Ｌ３により、界面層１４ａの厚さＷ１が求められる。この界面層１４ａの厚さＷ１は、コンピュータ１に記憶される。

このように、本発明のシミュレーション方法では、従来のシミュレーション方法のように、高分子材料及びフィラーを用いた実験を実施すること無く、界面層１４ａを特定し、さらにはその厚さＷ１を確実に計算することができる。従って、本発明のシミュレーション方法では、多くの時間や費用を必要とすることなく、界面層１４ａの厚さＷ１を正確に求めることができる。

しかも、本発明のシミュレーション方法では、粗視化モデル２を用いたシミュレーションにおいて、界面層１４ａの厚さＷ１を求めることができるため、現実に存在しない高分子鎖の界面層１４ａの厚さＷ１を求めることができ、未知の高分子材料の開発に役立つ。

さらに、本実施形態のシミュレーション方法では、粗視化モデル２のビーズ３、３間に、上記のような斥力ポテンシャルが定義されるため、シミュレーション工程Ｓ６において、実際の高分子材料の分子運動に近似させて、粗視化モデルの初期配置を効果的に緩和することができる。従って、本実施形態のシミュレーション方法では、界面層１４ａの厚さＷ１を精度よく計算することができる。

なお、界面厚計算工程Ｓ８では、緩和弾性率計算工程Ｓ８２に先立ち、コンピュータ１が、各粗視化モデル２のビーズ３、３間の距離Ｌ４（図４に示す）を同一に揃えられるのが望ましい（工程Ｓ８４）。これにより、粗視化モデル２の各ボンド５において、各ポテンシャルＰ１、Ｐ２のバラツキを抑制でき、緩和弾性率Ｇ（ｔ）を安定的に計算することができる。

さらに、各粗視化モデル２のビーズ３、３間の距離Ｌ４は、結合ポテンシャルＰ１又は結合ポテンシャルＰ２の平衡長と同一に設定されるのが望ましい。これにより、ビーズ３、３間には、各ポテンシャルＰ１、Ｐ２によってかかる力をゼロにでき、緩和弾性率Ｇ（ｔ）を安定的に計算することができる。

また、境界面１７、１７間の各距離Ｌ２は、粗視化モデル２（図４に示す）の慣性半径に基づいて設定されるのが望ましい。上述のとおり、慣性半径は、分子動力学計算において、粗視化モデル２の拡がりを示すパラメータである。界面層１４ａは、粗視化モデル２の鎖長によって変化すると考えられるため、各距離Ｌ２が慣性半径に基づいて設定されるのが望ましい。本実施形態の距離Ｌ２は、粗視化モデル２の慣性半径と同一に設定されている。

本実施形態の緩和弾性率計算工程Ｓ８２では、領域１６毎に、緩和弾性率Ｇ（ｔ）が計算されるものが例示されたが、これに限定されるわけではない。例えば、各フィラーモデル１１Ａ、１１Ｂから等距離にある一対の領域１６、１６を合わせて、緩和弾性率Ｇ（ｔ）が計算されてもよい。

例えば、第一領域１６Ａは、第Ｎ領域１６Ｎと合わせて緩和弾性率Ｇ（ｔ）が求められるのが望ましい。これにより、緩和弾性率Ｇ（ｔ）を計算する領域１６の数は、全ての領域１６Ａ〜１６Ｎの数の１／２に減少するため、計算時間を大幅に短縮することができる。さらに、緩和弾性率Ｇ（ｔ）は、約２倍の粗視化モデル２を対象にアンサンブル平均されるため、該緩和弾性率Ｇ（ｔ）を精度よく求めることができる。

なお、第一領域１６Ａ及び第Ｎ領域１６Ｎでは、Ｚ軸方向において、応力σ_xyは互いに逆向きとなる。このため、第一領域１６Ａ又は第Ｎ領域１６Ｎのいずれか一方の領域において、応力σ_xyのＺ軸方向の成分には「−１」を乗じて、緩和弾性率Ｇ（ｔ）が求められるのが望ましい。

次に、コンピュータ１が、高分子材料の物性値を計算する（工程Ｓ１０）。この工程Ｓ１０では、空間７に、界面層１４ａの厚さＷ１（図１２に示す）等の所定のパラメータが設定され、高分子材料の物性値（例えば、複素弾性率）が計算される。

次に、コンピュータ１が、高分子材料の物性値が、許容範囲内であるかを判断する（工程Ｓ１１）。この工程Ｓ１１では、物性値が許容範囲内であると判断された場合、粗視化モデル２を含む空間７の条件等に基づいて、高分子材料が製造される（工程Ｓ１２）。一方、物性値が許容範囲内でないと判断された場合は、粗視化モデル２を含む空間７の諸条件を変更して（工程Ｓ１３）、工程Ｓ６〜Ｓ１１が再度行われる。このように、本実施形態のシミュレーション方法では、高分子材料の物性値が許容範囲内になるまで、粗視化モデル２を含む空間７の諸条件が変更されるため、所望の性能を有する高分子材料を、効率よく設計することができる。

また、図６に示したように、本実施形態のフィラーモデル１１は、空間７の平面８でモデル化されるものが例示されたが、これに限定されるわけではない。例えば、図１４に示されるように、空間７の平面８に沿って複数の粒子１８が配置されるものでもよい。この場合、フィラーモデル１１の粒子１８と、粗視化モデル２のビーズ３（図７に示す）との間の相互ポテンシャルとしては、例えば、斥力ポテンシャルＲが定義されるのが望ましい。このようなフィラーモデル１１は、凹凸などの表面形状の影響を評価するのに役立つ。

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

図３に示される手順に従って、空間が区分された領域毎に緩和弾性率が計算され、界面層の厚さが求められた（実施例１）。また、各領域のうち、各フィラーモデルから等距離にある一対の領域を合わせて緩和弾性率が計算され、界面層の厚さが求められた（実施例２）。

さらに、比較のために、高分子材料及びフィラーを用いた実験結果に基づいて、界面層の厚さが求められた（比較例）。そして、実施例１、実施例２及び比較例において、界面層の厚さを求めるのに要した時間が測定された。なお、各ポテンシャルのパラメータ等は、明細書中の記載通りであり、共通仕様は次のとおりである。
粗視化モデル：
ビーズの個数：５０個
慣性半径：２．８９σ
空間：
１辺の長さＬ１（一対の平面間の距離Ｄ１）：２８．９σ
粗視化モデルの個数：３００

テストの結果、実施例１、実施例２及び比較例で計算された界面層の厚さは、次のとおりであった。実施例１、実施例２及び比較例の界面層の厚さは略同一であり、いずれも界面層の厚さを正確に求めうることが確認できた。
実施例１：８．６７σ
実施例２：８．６７σ
比較例：８．６７σ

また、実施例１の計算時間は、８０時間であった。一方、実施例２の計算時間は４０時間であった。従って、一対の領域を合わせて緩和弾性率が計算された実施例２は、領域毎に緩和弾性率が計算された実施例１に比べて、計算時間を短縮しうることが確認できた。

また、比較例の計算時間は、８００時間であった。従って、実施例１及び実施例２は、比較例に比べて、計算時間を大幅に短縮しうることが確認できた。

１ コンピュータ
２ 粗視化モデル
１１ フィラーモデル
１６ 領域
１７ 境界面

Claims (10)

1. コンピュータを用いて、高分子材料とフィラーとの反応を解析するための方法であって、
   前記コンピュータに、前記高分子材料の高分子鎖を、モノマー又はモノマーの一部分をなす構造単位をビーズで置換した粗視化モデルを設定する工程、
   前記コンピュータに、前記フィラーの外面を少なくとも含むフィラーモデルを設定する工程、
   前記コンピュータが、予め定められた空間内で前記粗視化モデルと前記フィラーモデルとを用いて分子動力学計算を行うシミュレーション工程、及び
   前記コンピュータが、前記シミュレーション工程の結果から、前記フィラーモデルの周囲に形成されかつ前記高分子材料のバルク部分とは異なる力学的特性を示す前記高分子材料の界面層の厚さを計算する界面厚計算工程を含み、
   前記粗視化モデルの前記ビーズ間には、該ビーズ間の距離が小さくなるほど無限に大きくなる斥力を含むポテンシャルが定義され、
   前記界面厚計算工程は、前記粗視化モデルが配置された前記空間を、前記フィラーモデルの外面に沿った境界面で複数の領域に区分する工程と、
   前記各領域の緩和弾性率を計算する緩和弾性率計算工程と、
   前記各領域の緩和弾性率に基づいて、前記界面層の厚さを求める特定工程とを含むことを特徴とする高分子材料のシミュレーション方法。
2. 前記特定工程は、隣り合う前記領域間での前記緩和弾性率が予め定めた範囲内になった領域以外を前記界面層として決定し、その厚さを計算する請求項１記載の高分子材料のシミュレーション方法。
3. 前記空間は、その外周を囲む一つの平面を含み、
   前記フィラーモデルは、前記平面でモデル化され、
   前記境界面は、前記フィラーモデルと平行である請求項１又は２に記載の高分子材料のシミュレーション方法。
4. 前記空間は、その外周を囲み、かつ互いに向き合う一対の平面を含み、
   前記フィラーモデルは、前記一対の平面でモデル化され、
   前記境界面は、一対の前記フィラーモデルと平行である請求項１又は２に記載の高分子材料のシミュレーション方法。
5. 前記粗視化モデルは、一対の前記フィラーモデルの間に配置される請求項４に記載の高分子材料のシミュレーション方法。
6. 前記一対の平面間の距離は、前記粗視化モデルの慣性半径の２倍以上である請求項４又は５に記載の高分子材料のシミュレーション方法。
7. 前記境界面間の各距離は、前記粗視化モデルの慣性半径に基づいて設定される請求項４乃至６のいずれかに記載の高分子材料のシミュレーション方法。
8. 前記緩和弾性率計算工程は、前記各フィラーモデルから等距離にある一対の前記領域を合わせて、前記緩和弾性率を計算する請求項４乃至７のいずれかに記載の高分子材料のシミュレーション方法。
9. 前記界面厚計算工程は、前記緩和弾性率計算工程に先立ち、前記各粗視化モデルの前記ビーズ間の距離を同一に揃える工程を含む請求項１乃至８のいずれかに記載の高分子材料のシミュレーション方法。
10. 前記斥力を含むポテンシャルは、下記式（１）で定義される請求項１乃至９のいずれかに記載の高分子材料のシミュレーション方法。
    
    ここで、各定数及び変数は次のとおりである。
    ｒ_ij：ビーズ間の距離
    ε：ビーズ間に定義される斥力ポテンシャルＲの強度に関する定数
    σ：ビーズ間に定義される斥力ポテンシャルＲが作用する距離に関する定数（分子動力学の分野では、ＬＪ球の直径と呼ばれる）