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

Abstract

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

Description

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

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

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

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

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

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

本発明のうち請求項１記載の発明は、コンピュータを用いて、高分子材料とフィラーとの反応を解析するための方法であって、前記コンピュータに、前記高分子材料の高分子鎖を、全原子モデル又はユナイテッドアトムモデルでモデル化したポリマーモデルを設定する工程、前記コンピュータに、前記フィラーの外面を少なくとも含むフィラーモデルを設定する工程、前記コンピュータが、予め定められた空間内で前記ポリマーモデルと前記フィラーモデルとを用いて分子動力学計算を行うシミュレーション工程、及び前記コンピュータが、前記シミュレーション工程の結果から、前記フィラーモデルの周囲に形成されかつ前記高分子材料のバルク部分とは異なる力学的特性を示す前記高分子材料の界面層の厚さを計算する界面厚計算工程を含み、前記界面厚計算工程は、前記ポリマーモデルが配置された前記空間を、前記フィラーモデルの外面に沿った境界面で複数の領域に区分する工程と、前記各領域の緩和弾性率を計算する緩和弾性率計算工程と、前記各領域の緩和弾性率に基づいて、前記界面層の厚さを求める特定工程とを含むことを特徴とする。

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

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

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

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

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

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

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

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

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

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

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

また、本発明のポリマーモデルは、高分子材料の高分子鎖を、全原子モデル又はユナイテッドアトムモデルでモデル化されるため、高分子鎖のモノマーの配置や、シス構造又はトランス構造に起因する界面層の厚さへの影響を解析することができる。

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

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

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

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

図３には、本実施形態のシミュレーション方法の具体的な処理手順が示されている。このシミュレーション方法では、先ず、コンピュータ１に、高分子材料の高分子鎖をモデル化したポリマーモデルが設定される（工程Ｓ１）。

図２及び図４に示されるように、ポリマーモデル２は、高分子材料の高分子鎖を、複数の粒子モデル３でモデル化される。本実施形態の粒子モデル３は、高分子鎖の炭素原子をモデル化した炭素モデル３Ｃ、及び、高分子鎖の水素原子をモデル化した水素モデル３Ｈを含む。従って、本実施形態のポリマーモデル２は、全ての炭素原子及び水素原子をモデル化した全原子モデルとして構成される。

粒子モデル３は、分子動力学計算に基づいたシミュレーションにおいて、運動方程式の質点として取り扱われる。即ち、粒子モデル３は、質量、直径、電荷、又は、初期座標などのパラメータが定義される。

また、粒子モデル３、３間には、結合鎖５が定義される。この結合鎖５は、粒子モデル３、３間を拘束するものである。本実施形態の結合鎖５は、炭素モデル３Ｃ、３Ｃを連結する主鎖５ａ、及び、炭素モデル３Ｃと水素モデル３Ｈとの間を連結する側鎖５ｂとを含んでいる。

図５（ａ）〜（ｂ）に示されるように、ポリマーモデル２には、各粒子モデル３、３間の結合長さである結合長ｒ、及び結合鎖５を介して連続する３つの粒子モデル３がなす角度である結合角θが定義される。さらに、図５（ｃ）に示されるように、ポリマーモデル２には、結合鎖５を介して連続する４つの粒子モデル３において、隣り合う３つの粒子モデル３が作る一方の平面９Ａと他方の平面９Ｂとのなす角度ある二面角φが定義される。なお、ポリマーモデル２は、外力又は内力によって、上記結合長ｒ、結合角θ及び二面角φが変化する。

このような結合長ｒ、結合角θ及び二面角φは、下記式（１）で定義される結合ポテンシャルＵ_bond（ｒ）、下記式（２）で定義される結合角ポテンシャルＵ_Angle（θ）、及び下記式（３）で定義される結合二面角ポテンシャルＵ_torsion（φ）によって設定される。



ここで、各定数及び変数は、次のとおりである。
ｒ：結合長
ｒ₀：平衡長
ｋ₁：ばね定数
θ：結合角
θ₀：平衡角度
k₂：二面角ポテンシャルの強度
Ｎ−1：二面角ポテンシャル多項式の次数
φ：二面角
Ａ_n：二面角定数
なお、結合長ｒ及び平衡長ｒ₀は、各粒子モデル３の中心（図示省略）間の距離として定義される。

ポリブタジエンのモノマーは−［ＣＨ_２−ＣＨ＝ＣＨ−ＣＨ_２］一である。従って、結合ポテンシャルＵ_bond（ｒ）は、ＣＨ＝ＣＨ、ＣＨ_２−ＣＨ_２、Ｃ−Ｈ、及びＣＨ−ＣＨ_２にそれぞれ定義される。各結合ポテンシャルＵ_bond（ｒ）は、いずれも調和振動子として長さの変化に対するポテンシャルが表現される。また、各結合ポテンシャルＵ_bond（ｒ）の定数は、適宜設定することができる。本実施形態では、論文（J. Phys. Chem. 94, 8897 (1990））に基づいて、以下のように設定されるのが望ましい。
ＣＨ＝ＣＨ：
ｋ₁＝１４００（kcal／mol・Å^２）、ｒ₀＝１．３３（Å）
ＣＨ_２−ＣＨ_２：
ｋ₁＝７００（kcal／mol・Å^２）、ｒ₀＝１．５３（Å）
ＣH₃−Ｈ：
ｋ₁＝７００（kcal／mol・Å^２）、ｒ₀＝１．０９（Å）
ＣH₂−Ｈ：
ｋ₁＝７００（kcal／mol・Å^２）、ｒ₀＝０．９９（Å）
ＣＨ−ＣＨ_２：
ｋ₁＝７００（kcal／mol・Å^２）、ｒ₀＝１．４３（Å）

また、結合角ポテンシャルＵ_Angle（θ）は、ＣＨ_２−ＣＨ_２−ＣＨ、及びＣＨ＝ＣＨ−ＣＨ_２等にそれぞれ定義される。いずれの結合角ポテンシャルＵ_Angle（θ）も、調和振動子として角度の変化に対するポテンシャルが表現される。各結合角ポテンシャルＵ_Angle（θ）の定数は、適宜設定することができるが、上記論文に基づいて、以下のように設定されるのが望ましい。
ＣＨ_２−ＣＨ_２−ＣＨ：
ｋ₁＝１００（kcal／mol・rad^２）、θ₀＝１０９．５（deg）
ＣＨ＝ＣＨ−ＣＨ_２：
ｋ₁＝１００（kcal／mol・rad^２）、θ₀＝１２０（deg）

さらに、結合二面角ポテンシャルＵ_torsion（φ）は、一つの結合鎖５を中心軸ＣＬ（図５（ｃ）に示す）とする分子鎖内の回転を表すポテンシャルである。この結合二面角ポテンシャルＵ_torsion（φ）は、その中心軸ＣＬで区別すると、ＣＨ＝ＣＨ、ＣＨ_２−ＣＨ_２等にそれぞれ定義される。各結合二面角ポテンシャルＵ_torsion（φ）の定数は、適宜設定することができるが、上記論文に基づいて、例えば、以下のように設定されるのが望ましい。
ＣＨ_２−ＣＨ_２：
ｋ₂=０．１１１ (kcal/mol)、N=４、A₀＝１、A₁＝３、A₂＝０、A₃＝−４

図６に示されるように、隣接するポリマーモデル２、２の粒子モデル３、３間には、下記式（４）で定義されるポテンシャルＲが定義される。

ここで、各定数及び変数は、Lennard-Jonesポテンシャルのパラメータであり、次のとおりである。
ε：粒子モデル間に定義されるポテンシャルＲの強度に関する定数
σ：粒子モデル間に定義されるポテンシャルＲが作用する距離に関する定数（分子動力学の分野では、ＬＪ球の直径と呼ばれる）
ｒij：粒子モデル間の距離
ｒ_c：カットオフ距離
なお、距離ｒij及びカットオフ距離ｒ_cは、各粒子モデル３の中心（図示省略）間の距離として定義される。

粒子モデル３間のポテンシャルＲは、炭素モデル３Ｃ、３Ｃ間に設定される第１ポテンシャルＲ１、水素モデル３Ｈ、３Ｈ間に設定される第２ポテンシャルＲ２、及び、炭素モデル３Ｃと水素モデル３Ｈとの間に設定される第３ポテンシャルＲ３を含む。上記式（４）中の各定数は、適宜設定することができる。本実施形態では、上記論文に基づいて、ポテンシャルＲ１〜Ｒ３毎に以下のように設定される。
第１ポテンシャルＲ１：
ε＝０．０９５１(kcal/mol)、σ＝３．４７３(Å)、ｒ_c=８．６８(Å)
第２ポテンシャルＲ２：
ε＝０．０１５２(kcal/mol)、σ＝２．８４６(Å)、ｒ_c=７．１２(Å)
第３ポテンシャルＲ３：
ε＝０．０３８１(kcal/mol)、σ＝３．１６０(Å)、ｒ_c=７．９０(Å)

このようなポリマーモデル２は、例えば（株）ＪＳＯＬ社製のＪ−ＯＣＴＡというソフトウエアを用いて作成することができる。ポリマーモデル２は、コンピュータ１で取り扱い可能な数値データであり、コンピュータ１に記憶される。

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

また、一対の平面８、８に直交する方向において、該一対の平面８、８間の距離Ｄ１（即ち、１辺の長さＬ１）は、ポリマーモデル２（図４に示す）の慣性半径（図示省略）の２倍以上、より好ましくは４倍以上に設定されるのが望ましい。慣性半径は、分子動力学計算において、ポリマーモデル２の拡がりを示すパラメータである。このような空間７では、分子動力学計算において、ポリマーモデル２の回転運動をスムーズに計算することができる。また、空間７の大きさは、例えば１(atm)で安定な体積に設定される。このような空間７は、コンピュータ１で取り扱い可能な数値データであり、コンピュータ１に記憶される。

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

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

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

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

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

また、相互ポテンシャルＴは、炭素モデル３Ｃとフィラーモデル１１との間の第１相互ポテンシャルＴ１、及び水素モデル３Ｈとフィラーモデル１１との間の第２相互ポテンシャルを含む。これらの相互ポテンシャルＴ１、Ｔ２の各定数は、適宜設定することができる。本実施形態では、以下のように設定される。
第１相互ポテンシャルＴ１：
ρ_wall＝１．０、σ_wall＝３．４７３(Å)、ε_wall＝０．０９５１(kcal/mol)、ｒ_c＝８．６８(Å)
第２相互ポテンシャルＴ２：
ρ_wall＝１．０、σ_wall＝３．１６０(Å)、ε_wall＝０．０３８１(kcal/mol)、ｒ_c＝７．９０(Å)

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

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

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

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

本実施形態のシミュレーション方法では、ポリマーモデルが全原子モデルとしてモデル化されるため、実際の高分子材料の分子運動に近似させることができ、ポリマーモデル２の初期配置を精度よく緩和することができる。シミュレーション工程Ｓ６の計算結果は、単位時間毎に、コンピュータ１に記憶される。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

さらに、本実施形態のポリマーモデル２は、全原子モデルでモデル化されるため、高分子鎖のモノマーの配置や、シス構造又はトランス構造に起因する界面層の厚さへの影響を解析することができる。

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

さらに、各ポリマーモデル２の粒子モデル３、３間の距離Ｌ４は、結合ポテンシャルＵ_bond（ｒ）の平衡長ｒ₀と同一に設定されるのが望ましい。これにより、結合鎖５は、結合ポテンシャルＵ_bond（ｒ）がゼロに設定され、緩和弾性率Ｇ（ｔ）を安定的に計算することができる。

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

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

この場合、第一領域１６Ａは、第Ｎ領域１６Ｎと合わせて緩和弾性率Ｇ（ｔ）が求められるのが望ましい。これにより、緩和弾性率Ｇ（ｔ）の計算回数が、全ての領域１６Ａ〜１６Ｎを計算する場合に比べて、１／２に減少するため、計算時間を大幅に短縮することができる。さらに、緩和弾性率Ｇ（ｔ）は、約２倍のポリマーモデル２を対象にアンサンブル平均されるため、該緩和弾性率Ｇ（ｔ）のバラツキを抑制することができる。

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

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

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

また、図４に示されるように、本実施形態のポリマーモデル２は、全原子モデルとして構成されるものが例示されたが、これに限定されるわけではない。図１５に示されるように、例えば、炭素原子に結合した水素原子を、該炭素原子と一体化して一つの粒子モデル３として扱うユナイテッドアトムモデル（united atom model ）として構成されてもよい。このようなポリマーモデル２は、高分子鎖のモノマーの配置や、シス構造又はトランス構造を維持しつつ、水素原子を省略することができるため、計算時間を短縮することができる。なお、各ポテンシャルのパラメータは、公開されている文献等に基づいて設定されるのが望ましい。

さらに、図７に示されるように、本実施形態のフィラーモデル１１は、空間７の平面８でモデル化されるものが例示されたが、これに限定されるわけではない。例えば、空間７内に配置される少なくとも一つの粒子（図示省略）でモデル化されるものでもよい。また、図１６に示されるように、空間７の平面８に沿って複数の粒子１８が配置されるものでもよい。この場合、フィラーモデル１１の各粒子１８と、ポリマーモデル２の各粒子モデル３（図８に示す）との間に、相互ポテンシャルＴを定義することができる。このようなフィラーモデル１１は、凹凸などの表面形状の影響を評価するのに役立つ。

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

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

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

テストの結果、実施例１、実施例２及び比較例で計算された界面層の厚さは、次のとおりであった。
実施例１：７．５Å
実施例２：７．５Å
比較例：７．５Å

実施例１、実施例２及び比較例の界面層の厚さは略同一であり、いずれも界面層の厚さを正確に求めうることが確認できた。

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

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

１ コンピュータ
２ ポリマーモデル
１１ フィラーモデル
１６ 領域
１７ 境界面

Claims (9)

1. コンピュータを用いて、高分子材料とフィラーとの反応を解析するための方法であって、
   前記コンピュータに、前記高分子材料の高分子鎖を、全原子モデル又はユナイテッドアトムモデルでモデル化したポリマーモデルを設定する工程、
   前記コンピュータに、前記フィラーの外面を少なくとも含むフィラーモデルを設定する工程、
   前記コンピュータが、予め定められた空間内で前記ポリマーモデルと前記フィラーモデルとを用いて分子動力学計算を行うシミュレーション工程、及び
   前記コンピュータが、前記シミュレーション工程の結果から、前記フィラーモデルの周囲に形成されかつ前記高分子材料のバルク部分とは異なる力学的特性を示す前記高分子材料の界面層の厚さを計算する界面厚計算工程を含み、
   前記界面厚計算工程は、前記ポリマーモデルが配置された前記空間を、前記フィラーモデルの外面に沿った境界面で複数の領域に区分する工程と、
   前記各領域の緩和弾性率を計算する緩和弾性率計算工程と、
   前記各領域の緩和弾性率に基づいて、前記界面層の厚さを求める特定工程とを含むことを特徴とする高分子材料のシミュレーション方法。
2. 前記特定工程は、隣り合う前記領域間での前記緩和弾性率が予め定めた範囲内になった領域以外を前記界面層として決定し、その厚さを計算する請求項１記載の高分子材料のシミュレーション方法。
3. 前記空間は、その外周を囲む一つの平面を含み、
   前記フィラーモデルは、前記平面でモデル化され、
   前記境界面は、前記フィラーモデルと平行である請求項１又は２に記載の高分子材料のシミュレーション方法。
4. 前記空間は、その外周を囲み、かつ互いに向き合う一対の平面を含み、
   前記フィラーモデルは、前記一対の平面でモデル化され、
   前記境界面は、一対の前記フィラーモデルと平行である請求項１又は２に記載の高分子材料のシミュレーション方法。
5. 前記ポリマーモデルは、一対の前記フィラーモデルの間に配置される請求項４に記載の高分子材料のシミュレーション方法。
6. 前記一対の平面間の距離は、前記ポリマーモデルの慣性半径の２倍以上である請求項４又は５に記載の高分子材料のシミュレーション方法。
7. 前記境界面間の各距離は、前記ポリマーモデルの慣性半径に基づいて設定される請求項４乃至６のいずれかに記載の高分子材料のシミュレーション方法。
8. 前記緩和弾性率計算工程は、前記各フィラーモデルから等距離にある一対の前記領域を合わせて、前記緩和弾性率を計算する請求項４乃至７のいずれかに記載の高分子材料のシミュレーション方法。
9. 前記ポリマーモデルは、複数の粒子モデルでモデル化され、
   前記界面厚計算工程は、前記緩和弾性率計算工程に先立ち、前記各ポリマーモデルの前記粒子モデル間の距離を同一に揃える工程を含む請求項１乃至８のいずれかに記載の高分子材料のシミュレーション方法。