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

Description

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

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

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

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

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

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

本発明は、コンピュータを用いて、高分子材料とフィラーとの反応を解析するための方法であって、前記コンピュータに、前記高分子材料の高分子鎖を、複数の粒子モデルと、前記粒子モデル間を結合するボンドとを含む全原子モデル又はユナイテッドアトムモデルとしてモデル化したポリマーモデルを設定する工程、前記コンピュータに、前記フィラーをフィラー粒子モデルでモデル化したフィラーモデルを設定する工程、前記コンピュータに、前記ポリマーモデルと前記フィラーモデルとを予め定められた空間内に配置して高分子材料モデルを設定する工程、前記コンピュータが、予め定めた条件に基づいて、前記高分子材料モデルの構造緩和を計算する工程、前記構造緩和の計算後、前記コンピュータが、前記高分子材料モデルを予め定められた歪で変形させて、前記各ボンドの配列状態を表す配向パラメータを求める計算工程、前記配向パラメータに基づいて、前記コンピュータが、前記フィラーの周囲に形成されかつ前記高分子材料のバルク部分とは異なる力学的特性を示す前記高分子材料の界面層を特定する界面特定工程を含むことを特徴とする。

本発明に係る前記高分子材料のシミュレーション方法は、前記計算工程は、前記ポリマーモデルが配置された前記空間を、前記フィラーモデルの外面に沿った境界面で複数の領域に区分する工程と、前記高分子材料モデルを前記歪で変形させる工程と、前記領域毎に、前記配向パラメータを計算する工程とを含むのが望ましい。

本発明に係る前記高分子材料のシミュレーション方法は、前記界面特定工程は、前記領域毎に前記配向パラメータの平均を求める工程と、前記配向パラメータの平均を、前記歪に対して同位相で応答する貯蔵成分と、前記歪に対して異なる位相で応答する散逸成分とに分解する工程と、前記貯蔵成分が前記散逸成分よりも大きい領域を前記界面層と決定する工程とを含むのが望ましい。

本発明に係る前記高分子材料のシミュレーション方法は、前記配向パラメータχ_i（ｍ，ｔ）は、前記ボンドに沿ったボンドベクトルｕ_i（ｍ，ｔ）と、前記ボンドの重心と前記フィラーモデルの重心とを結ぶ重心ベクトルｎ_i（ｍ，ｔ）とに基づいて、下記式（１）で定義されるのが望ましい。

ここで、各定数及び変数は次のとおりである。
ｉ：各ボンドのインデックス
ｍ：各領域のインデックス
ｔ：時刻
θ_i：ボンドベクトルｕ_i（ｍ，ｔ）と重心ベクトルｎ_i（ｍ，ｔ）とがなす角度

本発明に係る前記高分子材料のシミュレーション方法は、前記配向パラメータの平均Χ（ｍ，ｔ）は、下記式（２）で定義され、前記歪γは、下記式（３）で定義され、
前記配向パラメータの平均Χ（ｍ，ｔ）は、下記式（４）に基づいて、前記貯蔵成分Χ’と前記散逸成分Χ”とに分解されるのが望ましい。

ここで、各定数及び変数は次のとおりである。
ｉ：各ボンドのインデックス
ｍ：各領域のインデックス
ｔ：時刻
Ｎ：各領域に配置されるボンドの総数
γ₀：歪振幅
ω：角速度

本発明に係る前記高分子材料のシミュレーション方法は、前記各領域には、前記ボンドが１００個以上配置されているのが望ましい。

本発明の高分子材料のシミュレーション方法では、コンピュータが、高分子材料モデルを予め定められた歪で変形させて、各ボンドの配列状態を表す配向パラメータを求める計算工程、及び配向パラメータに基づいて、高分子材料の界面層を特定する界面特定工程を含んでいる。

一般に、高分子材料の界面層では、高分子鎖がフィラーに拘束されるため、高分子鎖がフィラーに拘束されないバルク部分に比べて変形しにくい傾向がある。本発明では、このような傾向に着目し、高分子材料モデルの変形に相関のあるボンドの配向パラメータを用いている。本発明では、配向パラメータに基づいて、高分子材料の界面層が特定される。従って、本発明では、高分子材料を用いた実験等を実施することなく、高分子材料の界面層を特定することができる。

さらに、本発明のポリマーモデルは、高分子材料の高分子鎖を、複数の粒子モデルと、粒子モデル間を結合するボンドとを含む全原子モデル又はユナイテッドアトムモデルとしてモデル化されている。このようなポリマーモデルは、例えば、高分子材料のモノマーを一つの粗視化粒子でモデル化した粗視化モデルに比べて、現実の高分子鎖に近似するため、高分子材料の界面層を精度よく特定することができる。

本実施形態のシミュレーション方法を実行するコンピュータの斜視図である。 ポリイソプレンの構造式である。 本実施形態の計算方法の処理手順の一例を示すフローチャートである。 全原子モデルとしてモデル化されたポリマーモデルの概念図である。 （ａ）〜（ｃ）は、ポテンシャルを説明するポリマーモデルの部分図である。 ポリマーモデルのポテンシャルを説明する概念図である。 ポリマーモデル及びフィラーモデルのポテンシャルを説明する概念図である。 高分子材料モデルの概念図である。 本実施形態の計算工程の処理手順の一例を示すフローチャートである。 空間が区分された領域の概念図である。 図１０の拡大図である。 界面層とバルク部分とを含む高分子材料の線図である。 本実施形態の界面特定工程の処理手順の一例を示すフローチャートである。 領域の貯蔵成分Χ’及び散逸成分Χ”、並びに、フィラーモデルの外面から各領域の境界面までの距離の関係を示すグラフである。 界面層、バルク部分及びそれらの境界の概念図である。 ユナイテッドアトムモデルとしてモデル化されたポリマーモデルの概念図である。

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

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

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

図３には、本実施形態のシミュレーション方法の具体的な処理手順が示されている。本実施形態のシミュレーション方法では、先ず、コンピュータ１に、高分子材料の高分子鎖をモデル化したポリマーモデルが設定される（工程Ｓ１）。図２及び図４に示されるように、ポリマーモデル２は、複数の粒子モデル３と、粒子モデル３、３間を結合するボンド４とを含んでいる。

本実施形態の粒子モデル３は、高分子鎖の炭素原子をモデル化した炭素モデル３Ｃ、及び、高分子鎖の水素原子をモデル化した水素モデル３Ｈを含んでいる。従って、本実施形態のポリマーモデル２は、高分子鎖の全ての炭素原子及び水素原子をモデル化した全原子モデルとして構成されている。

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

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

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

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

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

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

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

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

図６に示されるように、隣接するポリマーモデル２、２の粒子モデル３、３間には、相互作用ポテンシャルＰ１が定義される。相互作用ポテンシャルＰ１は、下記式（８）で定義されるＬＪポテンシャルＵ_LJ（ｒ_ij）である。

ここで、各定数及び変数は、Lennard-Jones ポテンシャルのパラメータであり、次のとおりである。
ｒ_ij：粒子モデル間の距離
ｒ_c：カットオフ距離
ε：粒子モデル間に定義されるＬＪポテンシャルの強度
σ：粒子モデルの直径に相当
なお、距離ｒ_ij及びカットオフ距離ｒ_cは、各粒子モデル３、３の中心３ｃ、３ｃ（図４に示す）間の距離として定義される。

相互作用ポテンシャルＰ１は、粒子間の距離ｒ_ijが2^1/6σよりも小さくなるほど、粒子間に作用する斥力が大きくなる。また、相互作用ポテンシャルＰ１は、粒子間の距離ｒ_ijが2^1/6σになるときに最小となり、粒子間の距離ｒ_ijが2^1/6σより長くなると粒子間に引力が作用する。さらに、相互作用ポテンシャルＰ１は、粒子間の距離ｒ_ijが2^6/7σよりも大になるほど、粒子間に作用する引力が小さくなる。このように、相互作用ポテンシャルＰ１は、粒子間の距離ｒ_ijに応じて、斥力及び引力を定義することができる。

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

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

次に、コンピュータ１に、フィラーをモデル化したフィラーモデルが設定される（工程Ｓ２）。図７に示されるように、フィラーモデル６は、例えば、一つのフィラー粒子モデル７でモデル化される。フィラー粒子モデル７は、粒子モデル３と同様に、分子動力学計算において、運動方程式の質点として取り扱われる。

粒子モデル３とフィラー粒子モデル７との間、及びフィラー粒子モデル７、７間には、相互作用ポテンシャルＰ２が定義される。この相互作用ポテンシャルＰ２は、上記式（８）で定義されるＬＪポテンシャルＵ_LJ（ｒ_ij）である。このような相互作用ポテンシャルＰ２は、粒子間の距離ｒ_ijに応じて、斥力及び引力を定義することができる。

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

次に、予め定められた空間内に、ポリマーモデル２及びフィラーモデル６が配置されて、高分子材料モデル１０が設定される（工程Ｓ３）。図８に示されるように、本実施形態の空間８は、互いに向き合う三対の平面１１、１１を有する立方体として定義されている。各平面１１には、周期境界条件が定義されている。これにより、空間８は、一方の平面１１が、反対側の平面１１と繋がっているものとして取り扱われる。

空間８の一辺の長さＬ１は、適宜設定することができるが、ポリマーモデル２の慣性半径（図示省略）の２倍以上が望ましい。慣性半径は、分子動力学計算において、ポリマーモデル２の拡がりを示すパラメータである。このような空間８では、分子動力学計算において、ポリマーモデル２の運動をスムーズに計算することができる。さらに、空間８の大きさは、例えば１atmで安定な体積に設定される。このような空間８は、コンピュータ１で取り扱い可能な数値データであり、コンピュータ１に記憶される。

空間８には、ポリマーモデル２及びフィラーモデル６が複数個配置される。本実施形態では、ポリマーモデル２が１００〜１０００本程度、フィラーモデル６が４〜３２個程度配置されている。このように、ポリマーモデル２及びフィラーモデル６が配置された空間８は、解析対象の高分子材料の微小構造部分である高分子材料モデル１０として構成される。このような高分子材料モデル１０も数値データであり、コンピュータ１に記憶される。

次に、コンピュータ１が、予め定められた条件に基づいて、高分子材料モデル１０の構造緩和を計算する（工程Ｓ４）。本実施形態の分子動力学計算では、例えば、空間８について所定の時間、ポリマーモデル２及びフィラーモデル６が古典力学に従うものとして、ニュートンの運動方程式が適用される。そして、各時刻での粒子モデル３（図７に示す）及びフィラー粒子モデル７（図７に示す）の動きが、単位時間毎に追跡される。このような構造緩和の計算は、例えばソフトマテリアル総合シミュレーター（ＯＣＴＡ）に含まれるＣＯＧＮＡＣを用いて処理することができる。

構造緩和の計算は、空間８において、圧力及び温度が一定、又は体積及び温度が一定に保たれる。これにより、工程Ｓ４では、実際の高分子材料の分子運動に近似させて、ポリマーモデル２及びフィラー粒子モデル７の初期配置を精度よく緩和することができる。

本実施形態のポリマーモデル２は、全原子モデルとしてモデル化されているため、例えば、モノマーを一つの粗視化粒子でモデル化した粗視化モデルに比べて、現実の高分子鎖に近似する。従って、本実施形態の工程Ｓ４では、ポリマーモデル２の初期配置を精度よく緩和することができる。

次に、コンピュータ１が、ポリマーモデル２及びフィラーモデル６の初期配置を十分に緩和できたか否かを判断する（工程Ｓ５）。この工程Ｓ５では、ポリマーモデル２及びフィラーモデル６の初期配置を十分に緩和できたと判断された場合、次の計算工程Ｓ６が実施される。一方、ポリマーモデル２及びフィラーモデル６の初期配置を十分に緩和できていないと判断された場合は、単位時間を進めて（工程Ｓ７）、工程Ｓ４が再度実施される。これにより、工程Ｓ４では、ポリマーモデル２及びフィラーモデル６の平衡状態（構造が緩和した状態）を確実に計算することができる。なお、構造緩和の計算時間は、上記条件の下で、末端間ベクトルの緩和時間をτ_ｎとした場合、少なくともτ_ｎ以上、より好ましくは１０τ_ｎ以上が望ましい。

次に、構造緩和の計算後、コンピュータ１が、高分子材料モデル１０を予め定められた歪で変形させて、各ボンド４の配向パラメータを求める（計算工程Ｓ６）。図９には、本実施形態の計算工程Ｓ６の処理手順の一例が示されている。

本実施形態の計算工程Ｓ６では、先ず、図１０に示されるように、高分子材料モデル１０において、ポリマーモデル２が配置された空間８が、複数の領域Ｔに区分される（工程Ｓ６１）。各領域Ｔは、フィラーモデル６の外面６ｓに沿った複数の境界面Ｔｓで区分されている。また、各境界面Ｔｓは、フィラーモデル６の重心６ｇを中心として、一定の厚さｄｗで同心円状に配置されている。

本実施形態の領域Ｔは、フィラーモデル６に最も隣接する第１番目の第１領域Ｔ１から第Ｎ番目の第Ｎ領域ＴＮまでのＮ個分設定される。これらの領域Ｔは、座標値等の数値データであり、コンピュータ１に記憶される。

次に、高分子材料モデル１０（図８に示す）が予め定められた歪で変形される（工程Ｓ６２）。本実施形態の工程Ｓ６２では、図８に示される初期の高分子材料モデル１０に引張変形を与えた後、該引張変形と逆方向に同一の歪で圧縮変形を与えて初期の状態に戻す工程を１周期とする引張・圧縮変形がシミュレーションされる。本実施形態の歪γは、下記式（３）で定義される。

ここで、各定数及び変数は次のとおりである。
γ₀：歪振幅
ω：角速度
ｔ：時刻

上記歪γは、正弦波形の振動歪である。このような振動歪は、例えば、動的粘弾性測定において、解析対象の高分子材料に加えられる歪に近似させることができる。

空間８において、フィラーモデル６近傍の領域Ｔ（例えば、第１領域Ｔ１〜第５領域Ｔ５）では、粒子モデル３とフィラー粒子モデル７との間に作用する相互作用ポテンシャルＰ２の引力が大きくなる。このため、フィラーモデル６近傍の領域Ｔでは、粒子モデル３が拘束され、上記歪γが定義されても、該領域Ｔの変形が小さくなる。

フィラーモデル６から大きく離間する領域Ｔ（例えば、第６領域Ｔ６〜第Ｎ領域ＴＮ）では、粒子モデル３とフィラー粒子モデル７との間に作用する相互作用ポテンシャルＰ２の引力が小さくなる。このため、フィラーモデル６から大きく離間する領域Ｔでは、粒子モデル３の動きが大きくなり、該領域Ｔでの変形が大きくなる。

このように、工程Ｓ６２では、フィラーが配合された高分子材料の変形を、高分子材料モデル１０を用いて精度良く再現することができる。また、本実施形態のポリマーモデル２は、全原子モデルとしてモデル化されているため、前記粗視化モデルに比べて、高分子材料の変形をより精度良く再現しうる。

工程Ｓ６２では、時刻ｔ毎に、高分子材料モデル１０の粒子モデル３、ボンド４及びフィラー粒子モデル７の座標値等が計算される。これらの座標値等は、数値データであり、コンピュータ１に記憶される。このような高分子材料モデル１０の変形計算も、例えば、ソフトマテリアル総合シミュレーター（ＯＣＴＡ）に含まれるＣＯＧＮＡＣを用いて処理することができる。

上記式（３）の各定数及び変数については、適宜設定することができるが、例えば、下記の値が設定されるのが望ましい。さらに、高分子材料モデル１０の変形計算は、例えば、からみ合い点の緩和時間をτ_eとした場合、一周期を０．０１〜１０τ_eの何れかで定義して計算されるのが望ましく、また、一回の変形計算で２０〜１００個の座標データが出力されるのが望ましい。
歪振幅γ₀：５％〜２０％
角速度ω：２π／（１０τ_e）〜２π／（０．０１τ_e）

次に、領域Ｔ毎に、各ボンド４の配向パラメータが計算される（工程Ｓ６３）。図１１に拡大して示されるように、工程Ｓ６３では、工程Ｓ６２で計算された粒子モデル３、ボンド４及びフィラー粒子モデル７の座標値に基づいて、配向パラメータが計算される。配向パラメータχ_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（ｍ，ｔ）を十分に均せないおそれがある。また、厚さｄｗが単位長さ１σを超えると、分解能が低下するおそれがある。

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

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

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

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

図３、図９及び図１３に示される手順に従って、全原子モデルのボンドの配向パラメータが求められ、高分子材料の界面層が特定された（実施例）。そして、実施例において、界面層の厚さを求めるのに必要な時間が測定された。なお、ポテンシャル、及び歪の各パラメータ等は、明細書中の記載通りであり、共通仕様は次のとおりである。
実施例：
ポリマーモデル：
構造：cis-１，４ポリイソプレン
重合度：２０２
慣性半径：１．６nm
フィラーモデルの半径：３．０nm
空間：
１辺の長さＬ１：１０．３nm
ポリマーモデルの本数：３９本
フィラーモデルの個数：１個
各領域：
厚さｄｗ：０．５０nm
個数：１０個
各領域に配置されるボンドの個数：１２１８〜３９３１６個

テストの結果、実施例で計算された界面層の厚さは、４．２nmであった。従って、実施例では、高分子材料を用いた実験等を実施することなく、界面層の厚さを特定することができた。

２ ポリマーモデル
４ ボンド
６ フィラーモデル
１０ 高分子材料モデル

Claims (6)

1. コンピュータを用いて、高分子材料とフィラーとの反応を解析するための方法であって、
   前記コンピュータに、前記高分子材料の高分子鎖を、複数の粒子モデルと、前記粒子モデル間を結合するボンドとを含む全原子モデル又はユナイテッドアトムモデルとしてモデル化したポリマーモデルを設定する工程、
   前記コンピュータに、前記フィラーをフィラー粒子モデルでモデル化したフィラーモデルを設定する工程、
   前記コンピュータに、前記ポリマーモデルと前記フィラーモデルとを予め定められた空間内に配置して高分子材料モデルを設定する工程、
   前記コンピュータが、予め定めた条件に基づいて、前記高分子材料モデルの構造緩和を計算する工程、
   前記構造緩和の計算後、前記コンピュータが、前記高分子材料モデルを予め定められた歪で変形させて、前記各ボンドの配列状態を表す配向パラメータを求める計算工程、
   前記配向パラメータに基づいて、前記コンピュータが、前記フィラーの周囲に形成されかつ前記高分子材料のバルク部分とは異なる力学的特性を示す前記高分子材料の界面層を特定する界面特定工程を含むことを特徴とする高分子材料のシミュレーション方法。
2. 前記計算工程は、前記ポリマーモデルが配置された前記空間を、前記フィラーモデルの外面に沿った境界面で複数の領域に区分する工程と、
   前記高分子材料モデルを前記歪で変形させる工程と、
   前記領域毎に、前記配向パラメータを計算する工程とを含む請求項１に記載の高分子材料のシミュレーション方法。
3. 前記界面特定工程は、前記領域毎に前記配向パラメータの平均を求める工程と、
   前記配向パラメータの平均を、前記歪に対して同位相で応答する貯蔵成分と、前記歪に対して異なる位相で応答する散逸成分とに分解する工程と、
   前記貯蔵成分が前記散逸成分よりも大きい領域を前記界面層と決定する工程とを含む請求項２に記載の高分子材料のシミュレーション方法。
4. 前記配向パラメータχ_i（ｍ，ｔ）は、前記ボンドに沿ったボンドベクトルｕ_i（ｍ，ｔ）と、前記ボンドの重心と前記フィラーモデルの重心とを結ぶ重心ベクトルｎ_i（ｍ，ｔ）とに基づいて、下記式（１）で定義される請求項１乃至３のいずれかに記載の高分子材料のシミュレーション方法。
   

   

   
   ここで、各定数及び変数は次のとおりである。
   ｉ：各ボンドのインデックス
   ｍ：各領域のインデックス
   ｔ：時刻
   θ_i：ボンドベクトルｕ_i（ｍ，ｔ）と重心ベクトルｎ_i（ｍ，ｔ）とがなす角度
5. 前記配向パラメータの平均Χ（ｍ，ｔ）は、下記式（２）で定義され、
   前記歪γは、下記式（３）で定義され、
   前記配向パラメータの平均Χ（ｍ，ｔ）は、下記式（４）に基づいて、前記貯蔵成分Χ’と前記散逸成分Χ”とに分解される請求項３に記載の高分子材料のシミュレーション方法。
   

   

   
   

   

   
   

   

   
   ここで、各定数及び変数は次のとおりである。
   ｉ：各ボンドのインデックス
   ｍ：各領域のインデックス
   ｔ：時刻
   Ｎ：各領域に配置されるボンドの総数
   γ₀：歪振幅
   ω：角速度
6. 前記各領域には、前記ボンドが１００個以上配置されている請求項２乃至５のいずれかに記載の高分子材料のシミュレーション方法。

Images