(57)【要約】【課題】フィラーと高分子材料とが共存する系において、フィラーの周囲に形成される界面層を特定する。【解決手段】コンピュータ１を用いて、高分子材料とフィラーとの反応を解析するための方法である。この方法では、ポリマーモデル２とフィラーモデル６とを予め定められた空間８内に配置して高分子材料モデル１０を設定する工程Ｓ３、予め定めた条件に基づいて、高分子材料モデル１０の構造緩和を計算する工程Ｓ４、構造緩和の計算後、高分子材料モデル１０を予め定められた歪で変形させて、各ボンド４の配列状態を表す配向パラメータを求める計算工程Ｓ６、配向パラメータに基づいて、フィラーの周囲に形成され高分子材料の界面層を特定する界面特定工程Ｓ８を含む。ポリマーモデル２は、全原子モデル又はユナイテッドアトムモデルとしてモデル化される。【選択図】図３
2018112525.tif 000002
【特許請求の範囲】
【請求項１】
コンピュータを用いて、高分子材料とフィラーとの反応を解析するための方法であって、
前記コンピュータに、前記高分子材料の高分子鎖を、複数の粒子モデルと、前記粒子モデル間を結合するボンドとを含む全原子モデル又はユナイテッドアトムモデルとしてモデル化したポリマーモデルを設定する工程、
前記コンピュータに、前記フィラーをフィラー粒子モデルでモデル化したフィラーモデルを設定する工程、
前記コンピュータに、前記ポリマーモデルと前記フィラーモデルとを予め定められた空間内に配置して高分子材料モデルを設定する工程、
前記コンピュータが、予め定めた条件に基づいて、前記高分子材料モデルの構造緩和を計算する工程、
前記構造緩和の計算後、前記コンピュータが、前記高分子材料モデルを予め定められた歪で変形させて、前記各ボンドの配列状態を表す配向パラメータを求める計算工程、
前記配向パラメータに基づいて、前記コンピュータが、前記フィラーの周囲に形成されかつ前記高分子材料のバルク部分とは異なる力学的特性を示す前記高分子材料の界面層を特定する界面特定工程を含むことを特徴とする高分子材料のシミュレーション方法。
【請求項２】
前記計算工程は、前記ポリマーモデルが配置された前記空間を、前記フィラーモデルの外面に沿った境界面で複数の領域に区分する工程と、
前記高分子材料モデルを前記歪で変形させる工程と、
前記領域毎に、前記配向パラメータを計算する工程とを含む請求項１に記載の高分子材料のシミュレーション方法。
【請求項３】
前記界面特定工程は、前記領域毎に前記配向パラメータの平均を求める工程と、
前記配向パラメータの平均を、前記歪に対して同位相で応答する貯蔵成分と、前記歪に対して異なる位相で応答する散逸成分とに分解する工程と、
前記貯蔵成分が前記散逸成分よりも大きい領域を前記界面層と決定する工程とを含む請求項２に記載の高分子材料のシミュレーション方法。
【請求項４】
前記配向パラメータχi（ｍ，ｔ）は、前記ボンドに沿ったボンドベクトルｕi（ｍ，ｔ）と、前記ボンドの重心と前記フィラーモデルの重心とを結ぶ重心ベクトルｎi（ｍ，ｔ）とに基づいて、下記式（１）で定義される請求項１乃至３のいずれかに記載の高分子材料のシミュレーション方法。
【数１】
2018112525.tif 000016
ここで、各定数及び変数は次のとおりである。
ｉ：各ボンドのインデックス
ｍ：各領域のインデックス
ｔ：時刻
θi：ボンドベクトルｕi（ｍ，ｔ）と重心ベクトルｎi（ｍ，ｔ）とがなす角度
【請求項５】
前記配向パラメータの平均Χ（ｍ，ｔ）は、下記式（２）で定義され、
前記歪γは、下記式（３）で定義され、
前記配向パラメータの平均Χ（ｍ，ｔ）は、下記式（４）に基づいて、前記貯蔵成分Χ’と前記散逸成分Χ”とに分解される請求項４に記載の高分子材料のシミュレーション方法。
【数２】
2018112525.tif 000017
【数３】
2018112525.tif 000018
【数４】
2018112525.tif 000019
ここで、各定数及び変数は次のとおりである。
ｉ：各ボンドのインデックス
ｍ：各領域のインデックス
ｔ：時刻
Ｎ：各領域に配置されるボンドの総数
γ0：歪振幅
ω：角速度
【請求項６】
前記各領域には、前記ボンドが１００個以上配置されている請求項２乃至５のいずれかに記載の高分子材料のシミュレーション方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、フィラーと高分子材料とが共存する系において、フィラーの周囲に形成される界面層を特定することができる高分子材料のシミュレーション方法に関する。
【背景技術】
【０００２】
近年、ゴム等の高分子材料の設計・開発のために、コンピュータを用いた高分子材料のシミュレーション方法が種々提案されている（例えば、下記特許文献１参照）。この種のシミュレーション方法では、高分子材料の構造や、フィラーの配合率等の諸条件を、計算に織り込むことができる。従って、このシミュレーション方法では、実際に高分子材料を試作することなく、その物性値を計算することができる。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開２０１２−２３８１６８号公報
【発明の概要】
【発明が解決しようとする課題】
【０００４】
ところで、フィラーと高分子材料とが共存する系においては、フィラーの周囲に界面層が形成されることが判明している。この界面層は、高分子材料本来の部分（以下、「バルク部分」という。）とは異なる力学的特性を示すものとして知られており、「ガラス層」などと称されることもある。
【０００５】
従来のシミュレーション方法でも、高分子材料の物性値を計算するに先立ち、界面層が定義されていた。しかしながら、従来のシミュレーション方法では、例えば、高分子材料及びフィラーを用いた実験結果に基づいて、界面層が定義されていたため、多くの時間や費用が必要になるという問題点があった。従って、界面層を定量的に特定することができる高分子材料のシミュレーション方法が強く求められていた。
【０００６】
本発明は、以上のような実状に鑑み案出されたもので、フィラーと高分子材料とが共存する系において、フィラーの周囲に形成される界面層を特定することができる高分子材料のシミュレーション方法を提供することを主たる目的としている。
【課題を解決するための手段】
【０００７】
本発明は、コンピュータを用いて、高分子材料とフィラーとの反応を解析するための方法であって、前記コンピュータに、前記高分子材料の高分子鎖を、複数の粒子モデルと、前記粒子モデル間を結合するボンドとを含む全原子モデル又はユナイテッドアトムモデルとしてモデル化したポリマーモデルを設定する工程、前記コンピュータに、前記フィラーをフィラー粒子モデルでモデル化したフィラーモデルを設定する工程、前記コンピュータに、前記ポリマーモデルと前記フィラーモデルとを予め定められた空間内に配置して高分子材料モデルを設定する工程、前記コンピュータが、予め定めた条件に基づいて、前記高分子材料モデルの構造緩和を計算する工程、前記構造緩和の計算後、前記コンピュータが、前記高分子材料モデルを予め定められた歪で変形させて、前記各ボンドの配列状態を表す配向パラメータを求める計算工程、前記配向パラメータに基づいて、前記コンピュータが、前記フィラーの周囲に形成されかつ前記高分子材料のバルク部分とは異なる力学的特性を示す前記高分子材料の界面層を特定する界面特定工程を含むことを特徴とする。
【０００８】
本発明に係る前記高分子材料のシミュレーション方法は、前記計算工程は、前記ポリマーモデルが配置された前記空間を、前記フィラーモデルの外面に沿った境界面で複数の領域に区分する工程と、前記高分子材料モデルを前記歪で変形させる工程と、前記領域毎に、前記配向パラメータを計算する工程とを含むのが望ましい。
【０００９】
本発明に係る前記高分子材料のシミュレーション方法は、前記界面特定工程は、前記領域毎に前記配向パラメータの平均を求める工程と、前記配向パラメータの平均を、前記歪に対して同位相で応答する貯蔵成分と、前記歪に対して異なる位相で応答する散逸成分とに分解する工程と、前記貯蔵成分が前記散逸成分よりも大きい領域を前記界面層と決定する工程とを含むのが望ましい。
【００１０】
本発明に係る前記高分子材料のシミュレーション方法は、前記配向パラメータχi（ｍ，ｔ）は、前記ボンドに沿ったボンドベクトルｕi（ｍ，ｔ）と、前記ボンドの重心と前記フィラーモデルの重心とを結ぶ重心ベクトルｎi（ｍ，ｔ）とに基づいて、下記式（１）で定義されるのが望ましい。
【数１】
2018112525.tif 000003
ここで、各定数及び変数は次のとおりである。
ｉ：各ボンドのインデックス
ｍ：各領域のインデックス
ｔ：時刻
θi：ボンドベクトルｕi（ｍ，ｔ）と重心ベクトルｎi（ｍ，ｔ）とがなす角度
【００１１】
本発明に係る前記高分子材料のシミュレーション方法は、前記配向パラメータの平均Χ（ｍ，ｔ）は、下記式（２）で定義され、前記歪γは、下記式（３）で定義され、
前記配向パラメータの平均Χ（ｍ，ｔ）は、下記式（４）に基づいて、前記貯蔵成分Χ’と前記散逸成分Χ”とに分解されるのが望ましい。
【数２】
2018112525.tif 000004
【数３】
2018112525.tif 000005
【数４】
2018112525.tif 000006
ここで、各定数及び変数は次のとおりである。
ｉ：各ボンドのインデックス
ｍ：各領域のインデックス
ｔ：時刻
Ｎ：各領域に配置されるボンドの総数
γ0：歪振幅
ω：角速度
【００１２】
本発明に係る前記高分子材料のシミュレーション方法は、前記各領域には、前記ボンドが１００個以上配置されているのが望ましい。
【発明の効果】
【００１３】
本発明の高分子材料のシミュレーション方法では、コンピュータが、高分子材料モデルを予め定められた歪で変形させて、各ボンドの配列状態を表す配向パラメータを求める計算工程、及び配向パラメータに基づいて、高分子材料の界面層を特定する界面特定工程を含んでいる。
【００１４】
一般に、高分子材料の界面層では、高分子鎖がフィラーに拘束されるため、高分子鎖がフィラーに拘束されないバルク部分に比べて変形しにくい傾向がある。本発明では、このような傾向に着目し、高分子材料モデルの変形に相関のあるボンドの配向パラメータを用いている。本発明では、配向パラメータに基づいて、高分子材料の界面層が特定される。従って、本発明では、高分子材料を用いた実験等を実施することなく、高分子材料の界面層を特定することができる。
【００１５】
さらに、本発明のポリマーモデルは、高分子材料の高分子鎖を、複数の粒子モデルと、粒子モデル間を結合するボンドとを含む全原子モデル又はユナイテッドアトムモデルとしてモデル化されている。このようなポリマーモデルは、例えば、高分子材料のモノマーを一つの粗視化粒子でモデル化した粗視化モデルに比べて、現実の高分子鎖に近似するため、高分子材料の界面層を精度よく特定することができる。
【図面の簡単な説明】
【００１６】
【図１】本実施形態のシミュレーション方法を実行するコンピュータの斜視図である。
【図２】ポリイソプレンの構造式である。
【図３】本実施形態の計算方法の処理手順の一例を示すフローチャートである。
【図４】全原子モデルとしてモデル化されたポリマーモデルの概念図である。
【図５】（ａ）〜（ｃ）は、ポテンシャルを説明するポリマーモデルの部分図である。
【図６】ポリマーモデルのポテンシャルを説明する概念図である。
【図７】ポリマーモデル及びフィラーモデルのポテンシャルを説明する概念図である。
【図８】高分子材料モデルの概念図である。
【図９】本実施形態の計算工程の処理手順の一例を示すフローチャートである。
【図１０】空間が区分された領域の概念図である。
【図１１】図１０の拡大図である。
【図１２】界面層とバルク部分とを含む高分子材料の線図である。
【図１３】本実施形態の界面特定工程の処理手順の一例を示すフローチャートである。
【図１４】領域の貯蔵成分Χ’及び散逸成分Χ”、並びに、フィラーモデルの外面から各領域の境界面までの距離の関係を示すグラフである。
【図１５】界面層、バルク部分及びそれらの境界の概念図である。
【図１６】ユナイテッドアトムモデルとしてモデル化されたポリマーモデルの概念図である。
【発明を実施するための形態】
【００１７】
以下、本発明の実施の一形態が図面に基づき説明される。
本実施形態の高分子材料のシミュレーション方法（以下、単に「シミュレーション方法」ということがある）は、コンピュータを用いて、高分子材料とフィラーとの反応を解析するための方法である。
【００１８】
図１に示されるように、コンピュータ１は、本体１ａ、キーボード１ｂ、マウス１ｃ及びディスプレイ装置１ｄを含む。この本体１ａには、例えば、演算処理装置（ＣＰＵ）、ＲＯＭ、作業用メモリ、磁気ディスクなどの記憶装置、及びディスクドライブ装置１ａ１、１ａ２が設けられる。また、記憶装置には、本実施形態のシミュレーション方法を実行するための処理手順（プログラム）が予め記憶される。
【００１９】
高分子材料としては、例えば、ゴム、樹脂又はエラストマー等が含まれる。本実施形態では、図２に示されるように、高分子材料として、cis-１，４ポリイソプレン（以下、単に「ポリイソプレン」ということがある。）が例示される。このポリイソプレンを構成する高分子鎖は、メチン基等（−ＣＨ=、＞Ｃ＝）、メチレン基（−ＣＨ２−）、及び、メチル基（−ＣＨ３）によって構成されるイソプレンのモノマー（イソプレン分子）５が、重合度ｎで連結されて構成されている。なお、高分子材料には、ポリイソプレン以外の高分子材料が用いられてもよい。また、フィラーとしては、例えば、カーボンブラック、シリカ又はアルミナ等が含まれる。
【００２０】
図３には、本実施形態のシミュレーション方法の具体的な処理手順が示されている。本実施形態のシミュレーション方法では、先ず、コンピュータ１に、高分子材料の高分子鎖をモデル化したポリマーモデルが設定される（工程Ｓ１）。図２及び図４に示されるように、ポリマーモデル２は、複数の粒子モデル３と、粒子モデル３、３間を結合するボンド４とを含んでいる。
【００２１】
本実施形態の粒子モデル３は、高分子鎖の炭素原子をモデル化した炭素モデル３Ｃ、及び、高分子鎖の水素原子をモデル化した水素モデル３Ｈを含んでいる。従って、本実施形態のポリマーモデル２は、高分子鎖の全ての炭素原子及び水素原子をモデル化した全原子モデルとして構成されている。
【００２２】
粒子モデル３は、分子動力学計算に基づいたシミュレーションにおいて、運動方程式の質点として取り扱われる。即ち、粒子モデル３は、質量、直径、電荷、又は、初期座標などのパラメータが定義される。
【００２３】
ボンド４は、粒子モデル３、３間を拘束するものである。本実施形態のボンド４は、炭素モデル３Ｃ、３Ｃを連結する主鎖４ａ、及び、炭素モデル３Ｃと水素モデル３Ｈとの間を連結する側鎖４ｂとを含んでいる。
【００２４】
図５（ａ）〜（ｂ）に示されるように、ポリマーモデル２には、各粒子モデル３、３間の結合長さである結合長ｒ、及びボンド４を介して連続する３つの粒子モデル３がなす角度である結合角θが定義されている。さらに、図５（ｃ）に示されるように、ポリマーモデル２には、ボンド４を介して連続する４つの粒子モデル３において、隣り合う３つの粒子モデル３が作る一方の平面５Ａと他方の平面５Ｂとのなす角度ある二面角φが定義される。なお、ポリマーモデル２は、外力又は内力によって、上記結合長ｒ、結合角θ及び二面角φが変化する。
【００２５】
このような結合長ｒ、結合角θ及び二面角φは、下記式（５）で定義される結合ポテンシャルＵbond（ｒ）、下記式（６）で定義される結合角ポテンシャルＵAngle（θ）、及び下記式（７）で定義される結合二面角ポテンシャルＵtorsion（φ）によって設定される。
【数５】
2018112525.tif 000007
【数６】
2018112525.tif 000008
【数７】
2018112525.tif 000009
ここで、各定数及び変数は、次のとおりである。
ｒ：結合長
ｒ0：平衡長
ｋ1：ばね定数
θ：結合角
θ0：平衡角度
k2：二面角ポテンシャルの強度
Ｎ−1：二面角ポテンシャル多項式の次数
φ：二面角
Ａn：二面角定数
なお、結合長ｒ及び平衡長ｒ0は、各粒子モデル３の中心（図示省略）間の距離として定義される。
【００２６】
ポリイソプレンのモノマーは−［ＣＨ２−ＣＨ（ＣＨ3）＝ＣＨ−ＣＨ２］−である。従って、結合ポテンシャルＵbond（ｒ）は、ＣＨ＝ＣＨ、ＣＨ２−ＣＨ２、Ｃ−Ｈ、ＣＨ−ＣＨ3及びＣＨ−ＣＨ２にそれぞれ定義される。各結合ポテンシャルＵbond（ｒ）は、いずれも調和振動子として長さの変化に対するポテンシャルが表現される。また、各結合ポテンシャルＵbond（ｒ）の定数は、適宜設定することができる。本実施形態では、論文（J. Phys. Chem. 94, 8897 (1990））に基づいて、以下のように設定されるのが望ましい。
ＣＨ＝ＣＨ：
ｋ1＝１４００（kcal／mol・Å２）、ｒ0＝１．３３（Å）
ＣＨ２−ＣＨ２：
ｋ1＝７００（kcal／mol・Å２）、ｒ0＝１．５３（Å）
ＣH3−Ｈ：
ｋ1＝７００（kcal／mol・Å２）、ｒ0＝１．０９（Å）
ＣＨ２−Ｈ：
ｋ1＝７００（kcal／mol・Å２）、ｒ0＝０．９９（Å）
ＣＨ−ＣＨ3：
ｋ1＝７００（kcal／mol・Å２）、ｒ0＝１．４３（Å）
ＣＨ−ＣＨ２：
ｋ1＝７００（kcal／mol・Å２）、ｒ0＝１．４３（Å）
【００２７】
また、結合角ポテンシャルＵAngle（θ）は、ＣＨ２−ＣＨ２−ＣＨ、及びＣＨ＝ＣＨ−ＣＨ２等にそれぞれ定義される。いずれの結合角ポテンシャルＵAngle（θ）も、調和振動子として角度の変化に対するポテンシャルが表現される。各結合角ポテンシャルＵAngle（θ）の定数は、適宜設定することができるが、上記論文に基づいて、以下のように設定されるのが望ましい。
ＣＨ２−ＣＨ２−ＣＨ：
ｋ1＝１００（kcal／mol・rad２）、θ0＝１０９．５（deg）
ＣＨ＝ＣＨ−ＣＨ２：
ｋ1＝１００（kcal／mol・rad２）、θ0＝１２０（deg）
【００２８】
さらに、結合二面角ポテンシャルＵtorsion（φ）は、一つのボンド４を中心軸ＣＬ（図５（ｃ）に示す）とする分子鎖内の回転を表すポテンシャルである。この結合二面角ポテンシャルＵtorsion（φ）は、その中心軸ＣＬで区別すると、ＣＨ＝ＣＨ、ＣＨ２−ＣＨ２等にそれぞれ定義される。各結合二面角ポテンシャルＵtorsion（φ）の定数は、適宜設定することができるが、上記論文に基づいて、例えば、以下のように設定されるのが望ましい。
ＣＨ２−ＣＨ２：
k2=０．１１１ (kcal/mol)、Ｎ=４、Ａ0＝１、Ａ1＝３、Ａ2＝０、Ａ3＝−４
【００２９】
図６に示されるように、隣接するポリマーモデル２、２の粒子モデル３、３間には、相互作用ポテンシャルＰ１が定義される。相互作用ポテンシャルＰ１は、下記式（８）で定義されるＬＪポテンシャルＵLJ（ｒij）である。
【００３０】
【数８】
2018112525.tif 000010
ここで、各定数及び変数は、Lennard-Jones ポテンシャルのパラメータであり、次のとおりである。
ｒij：粒子モデル間の距離
ｒc：カットオフ距離
ε：粒子モデル間に定義されるＬＪポテンシャルの強度
σ：粒子モデルの直径に相当
なお、距離ｒij及びカットオフ距離ｒcは、各粒子モデル３、３の中心３ｃ、３ｃ（図４に示す）間の距離として定義される。
【００３１】
相互作用ポテンシャルＰ１は、粒子間の距離ｒijが21/6σよりも小さくなるほど、粒子間に作用する斥力が大きくなる。また、相互作用ポテンシャルＰ１は、粒子間の距離ｒijが21/6σになるときに最小となり、粒子間の距離ｒijが21/6σより長くなると粒子間に引力が作用する。さらに、相互作用ポテンシャルＰ１は、粒子間の距離ｒijが26/7σよりも大になるほど、粒子間に作用する引力が小さくなる。このように、相互作用ポテンシャルＰ１は、粒子間の距離ｒijに応じて、斥力及び引力を定義することができる。
【００３２】
相互作用ポテンシャルＰ１は、炭素モデル３Ｃ、３Ｃ間に設定される第１ポテンシャルＰ１ａ、水素モデル３Ｈ、３Ｈ間に設定される第２ポテンシャルＰ１ｂ、及び、炭素モデル３Ｃと水素モデル３Ｈとの間に設定される第３ポテンシャルＰ１ｃを含んでいる。上記式（４）中の各定数は、適宜設定することができる。本実施形態では、上記論文に基づいて、ポテンシャルＰ１ａ〜Ｐ１ｃ毎に以下のように設定される。
第１ポテンシャルＰ１ａ：
ε＝０．０９５１(kcal/mol)、σ＝３．４７３(Å)、ｒc=８．６８(Å)
第２ポテンシャルＰ１ｂ：
ε＝０．０１５２(kcal/mol)、σ＝２．８４６(Å)、ｒc=７．１２(Å)
第３ポテンシャルＰ１ｃ：
ε＝０．０３８１(kcal/mol)、σ＝３．１６０(Å)、ｒc=７．９０(Å)
【００３３】
このようなポリマーモデル２は、例えば（株）ＪＳＯＬ社製のＪ−ＯＣＴＡというソフトウエアを用いて作成することができる。ポリマーモデル２は、コンピュータ１で取り扱い可能な数値データであり、コンピュータ１に記憶される。
【００３４】
次に、コンピュータ１に、フィラーをモデル化したフィラーモデルが設定される（工程Ｓ２）。図７に示されるように、フィラーモデル６は、例えば、一つのフィラー粒子モデル７でモデル化される。フィラー粒子モデル７は、粒子モデル３と同様に、分子動力学計算において、運動方程式の質点として取り扱われる。
【００３５】
粒子モデル３とフィラー粒子モデル７との間、及びフィラー粒子モデル７、７間には、相互作用ポテンシャルＰ２が定義される。この相互作用ポテンシャルＰ２は、上記式（８）で定義されるＬＪポテンシャルＵLJ（ｒij）である。このような相互作用ポテンシャルＰ２は、粒子間の距離ｒijに応じて、斥力及び引力を定義することができる。
【００３６】
相互作用ポテンシャルＰ２は、フィラー粒子モデル７と炭素モデル３Ｃと間に設定される第１ポテンシャルＰ２ａと、フィラー粒子モデル７と水素モデル３Ｈとの間に設定される第２ポテンシャルＰ２ｂとが含まれる。上記式（８）中の各定数は、適宜設定することができる。本実施形態では、上記論文に基づいて、各ポテンシャルＰ２ａ、Ｐ２ｂ毎に以下のように設定される。
第１ポテンシャルＰ２ａ：
ε＝０．０９５１(kcal/mol)、σ＝３．４７３(Å)、ｒc=８．６８(Å)
第２ポテンシャルＰ２ｂ：
ε＝０．３８１(kcal/mol)、σ＝３．１６０(Å)、ｒc=７．９０(Å)
【００３７】
次に、予め定められた空間内に、ポリマーモデル２及びフィラーモデル６が配置されて、高分子材料モデル１０が設定される（工程Ｓ３）。図８に示されるように、本実施形態の空間８は、互いに向き合う三対の平面１１、１１を有する立方体として定義されている。各平面１１には、周期境界条件が定義されている。これにより、空間８は、一方の平面１１が、反対側の平面１１と繋がっているものとして取り扱われる。
【００３８】
空間８の一辺の長さＬ１は、適宜設定することができるが、ポリマーモデル２の慣性半径（図示省略）の２倍以上が望ましい。慣性半径は、分子動力学計算において、ポリマーモデル２の拡がりを示すパラメータである。このような空間８では、分子動力学計算において、ポリマーモデル２の運動をスムーズに計算することができる。さらに、空間８の大きさは、例えば１atmで安定な体積に設定される。このような空間８は、コンピュータ１で取り扱い可能な数値データであり、コンピュータ１に記憶される。
【００３９】
空間８には、ポリマーモデル２及びフィラーモデル６が複数個配置される。本実施形態では、ポリマーモデル２が１００〜１０００本程度、フィラーモデル６が４〜３２個程度配置されている。このように、ポリマーモデル２及びフィラーモデル６が配置された空間８は、解析対象の高分子材料の微小構造部分である高分子材料モデル１０として構成される。このような高分子材料モデル１０も数値データであり、コンピュータ１に記憶される。
【００４０】
次に、コンピュータ１が、予め定められた条件に基づいて、高分子材料モデル１０の構造緩和を計算する（工程Ｓ４）。本実施形態の分子動力学計算では、例えば、空間８について所定の時間、ポリマーモデル２及びフィラーモデル６が古典力学に従うものとして、ニュートンの運動方程式が適用される。そして、各時刻での粒子モデル３（図７に示す）及びフィラー粒子モデル７（図７に示す）の動きが、単位時間毎に追跡される。このような構造緩和の計算は、例えばソフトマテリアル総合シミュレーター（ＯＣＴＡ）に含まれるＣＯＧＮＡＣを用いて処理することができる。
【００４１】
構造緩和の計算は、空間８において、圧力及び温度が一定、又は体積及び温度が一定に保たれる。これにより、工程Ｓ４では、実際の高分子材料の分子運動に近似させて、ポリマーモデル２及びフィラー粒子モデル７の初期配置を精度よく緩和することができる。
【００４２】
本実施形態のポリマーモデル２は、全原子モデルとしてモデル化されているため、例えば、モノマーを一つの粗視化粒子でモデル化した粗視化モデルに比べて、現実の高分子鎖に近似する。従って、本実施形態の工程Ｓ４では、ポリマーモデル２の初期配置を精度よく緩和することができる。
【００４３】
次に、コンピュータ１が、ポリマーモデル２及びフィラーモデル６の初期配置を十分に緩和できたか否かを判断する（工程Ｓ５）。この工程Ｓ５では、ポリマーモデル２及びフィラーモデル６の初期配置を十分に緩和できたと判断された場合、次の計算工程Ｓ６が実施される。一方、ポリマーモデル２及びフィラーモデル６の初期配置を十分に緩和できていないと判断された場合は、単位時間を進めて（工程Ｓ７）、工程Ｓ４が再度実施される。これにより、工程Ｓ４では、ポリマーモデル２及びフィラーモデル６の平衡状態（構造が緩和した状態）を確実に計算することができる。なお、構造緩和の計算時間は、上記条件の下で、末端間ベクトルの緩和時間をτｎとした場合、少なくともτｎ以上、より好ましくは１０τｎ以上が望ましい。
【００４４】
次に、構造緩和の計算後、コンピュータ１が、高分子材料モデル１０を予め定められた歪で変形させて、各ボンド４の配向パラメータを求める（計算工程Ｓ６）。図９には、本実施形態の計算工程Ｓ６の処理手順の一例が示されている。
【００４５】
本実施形態の計算工程Ｓ６では、先ず、図１０に示されるように、高分子材料モデル１０において、ポリマーモデル２が配置された空間８が、複数の領域Ｔに区分される（工程Ｓ６１）。各領域Ｔは、フィラーモデル６の外面６ｓに沿った複数の境界面Ｔｓで区分されている。また、各境界面Ｔｓは、フィラーモデル６の重心６ｇを中心として、一定の厚さｄｗで同心円状に配置されている。
【００４６】
本実施形態の領域Ｔは、フィラーモデル６に最も隣接する第１番目の第１領域Ｔ１から第Ｎ番目の第Ｎ領域ＴＮまでのＮ個分設定される。これらの領域Ｔは、座標値等の数値データであり、コンピュータ１に記憶される。
【００４７】
次に、高分子材料モデル１０（図８に示す）が予め定められた歪で変形される（工程Ｓ６２）。本実施形態の工程Ｓ６２では、図８に示される初期の高分子材料モデル１０に引張変形を与えた後、該引張変形と逆方向に同一の歪で圧縮変形を与えて初期の状態に戻す工程を１周期とする引張・圧縮変形がシミュレーションされる。本実施形態の歪γは、下記式（３）で定義される。
【００４８】
【数３】
2018112525.tif 000011
ここで、各定数及び変数は次のとおりである。
γ0：歪振幅
ω：角速度
ｔ：時刻
【００４９】
上記歪γは、正弦波形の振動歪である。このような振動歪は、例えば、動的粘弾性測定において、解析対象の高分子材料に加えられる歪に近似させることができる。
【００５０】
空間８において、フィラーモデル６近傍の領域Ｔ（例えば、第１領域Ｔ１〜第５領域Ｔ５）では、粒子モデル３とフィラー粒子モデル７との間に作用する相互作用ポテンシャルＰ２の引力が大きくなる。このため、フィラーモデル６近傍の領域Ｔでは、粒子モデル３が拘束され、上記歪γが定義されても、該領域Ｔの変形が小さくなる。
【００５１】
フィラーモデル６から大きく離間する領域Ｔ（例えば、第６領域Ｔ６〜第Ｎ領域ＴＮ）では、粒子モデル３とフィラー粒子モデル７との間に作用する相互作用ポテンシャルＰ２の引力が小さくなる。このため、フィラーモデル６から大きく離間する領域Ｔでは、粒子モデル３の動きが大きくなり、該領域Ｔでの変形が大きくなる。
【００５２】
このように、工程Ｓ６２では、フィラーが配合された高分子材料の変形を、高分子材料モデル１０を用いて精度良く再現することができる。また、本実施形態のポリマーモデル２は、全原子モデルとしてモデル化されているため、前記粗視化モデルに比べて、高分子材料の変形をより精度良く再現しうる。
【００５３】
工程Ｓ６２では、時刻ｔ毎に、高分子材料モデル１０の粒子モデル３、ボンド４及びフィラー粒子モデル７の座標値等が計算される。これらの座標値等は、数値データであり、コンピュータ１に記憶される。このような高分子材料モデル１０の変形計算も、例えば、ソフトマテリアル総合シミュレーター（ＯＣＴＡ）に含まれるＣＯＧＮＡＣを用いて処理することができる。
【００５４】
上記式（３）の各定数及び変数については、適宜設定することができるが、例えば、下記の値が設定されるのが望ましい。さらに、高分子材料モデル１０の変形計算は、例えば、からみ合い点の緩和時間をτeとした場合、一周期を０．０１〜１０τeの何れかで定義して計算されるのが望ましく、また、一回の変形計算で２０〜１００個の座標データが出力されるのが望ましい。
歪振幅γ0：５％〜２０％
角速度ω：２π／（１０τe）〜２π／（０．０１τe）
【００５５】
次に、領域Ｔ毎に、各ボンド４の配向パラメータが計算される（工程Ｓ６３）。図１１に拡大して示されるように、工程Ｓ６３では、工程Ｓ６２で計算された粒子モデル３、ボンド４及びフィラー粒子モデル７の座標値に基づいて、配向パラメータが計算される。配向パラメータχi（ｍ，ｔ）は、下記式（１）で定義される。
【００５６】
【数１】
2018112525.tif 000012
ここで、各定数及び変数は次のとおりである。
ｉ：各ボンドのインデックス
ｍ：各領域のインデックス
ｔ：時刻
θi：ボンドベクトルｕi（ｍ，ｔ）と重心ベクトルｎi（ｍ，ｔ）とがなす角度（狭角）
【００５７】
配向パラメータχi（ｍ，ｔ）は、例えば、第ｍ番目の領域Ｔｍに配置されるボンド４に沿ったボンドベクトルｕi（ｍ，ｔ）、及び、該ボンド４の重心４ｇとフィラーモデル６の重心６ｇとを結ぶ重心ベクトルｎi（ｍ，ｔ）に基づいて定義される。このような配向パラメータχi（ｍ，ｔ）は、時刻ｔ毎に、フィラーモデル６に対するポリマーモデル２のボンド４の配列状態を定義することができる。
【００５８】
工程Ｓ６３では、空間８に配置される全てのボンド４（主鎖４ａ及び側鎖４ｂを含む）について、配向パラメータχi（ｍ，ｔ）が、時刻ｔ毎に計算される。これらの配向パラメータχi（ｍ，ｔ）は、数値データであり、コンピュータ１に記憶される。なお、ボンド４が隣り合う領域Ｔに跨って配置される場合は、ボンド４の重心４ｇが配置される領域Ｔにおいて、配向パラメータχi（ｍ，ｔ）が計算されるものとする。
【００５９】
次に、配向パラメータχi（ｍ，ｔ）の平均Χ（ｍ，ｔ）に基づいて、コンピュータ１が、高分子材料の界面層を特定する（界面特定工程Ｓ８）。
【００６０】
図１２に示されるように、フィラー１３と高分子材料１４との共存する系において、界面層１４ａは、フィラー１３の周囲に形成されることが知られている。また、界面層１４ａでは、高分子材料１４の本来の部分であるバルク部分１４ｂとは異なる力学的特性が示される。このため、図８に示される高分子材料モデル１０を用いて、高分子材料の物性値を計算する場合には、界面層１４ａの厚さＷ２を予め定義しておくことが重要である。
【００６１】
一般に、高分子材料１４の界面層１４ａでは、高分子鎖がフィラー１３に拘束されるため、高分子鎖がフィラー１３に拘束されないバルク部分１４ｂに比べて変形しにくい傾向がある。
【００６２】
本実施形態では、高分子材料モデル１０の変形に相関のあるボンド４の配向パラメータχi（ｍ，ｔ）の平均Χ（ｍ，ｔ）を用いて、高分子材料１４の界面層１４ａを特定している。図１３には、本実施形態の界面特定工程Ｓ８の具体的な処理手順が示されている。
【００６３】
本実施形態の界面特定工程Ｓ８では、先ず、領域Ｔ毎に、配向パラメータχi（ｍ，ｔ）の平均が求められる（工程Ｓ８１）。配向パラメータの平均Χ（ｍ，ｔ）は、下記式（２）で定義される。
【００６４】
【数２】
2018112525.tif 000013
ここで、各定数及び変数は次のとおりである。
ｉ：各ボンドのインデックス
ｍ：各領域のインデックス
ｔ：時刻
ｄｗ：各領域の厚さ
Ｎ：各領域に配置されるボンドの総数
【００６５】
上記式（２）では、第ｍ番目の領域Ｔの配向パラメータχi（ｍ，ｔ）の平均が、時刻ｔ毎に計算される。本実施形態では、第１領域Ｔ１から第Ｎ領域ＴＮまでの全ての領域Ｔにおいて、配向パラメータの平均Χ（ｍ，ｔ）が計算される。これらの配向パラメータの平均Χ（ｍ，ｔ）は、数値データであり、コンピュータ１に記憶される。
【００６６】
上述したように、フィラーモデル６近傍の領域Ｔでは、粒子モデル３とフィラー粒子モデル７との間に作用する相互作用ポテンシャルＰ２の引力が大きくなる。このため、粒子モデル３、３間のボンド４に定義されるボンドベクトルｕi（ｍ，ｔ）と重心ベクトルｎi（ｍ，ｔ）とがなす角度θiが大となり、配向パラメータχi（ｍ，ｔ）が１に近づく。
【００６７】
一方、フィラーモデル６から大きく離間する領域Ｔでは、粒子モデル３とフィラー粒子モデル７との間に作用する相互作用ポテンシャルＰ２の引力が小さくなる。このため、粒子モデル３、３間のボンド４に定義されるボンドベクトルｕi（ｍ，ｔ）と重心ベクトルｎi（ｍ，ｔ）とがなす角度θiが乱雑となり（即ち、拘束力が弱いため、角度θｉがランダムな方向に向く）、配向パラメータが−１〜＋１の間で様々な値をとるため、配向パラメータの平均Χ（ｍ，ｔ）が０に近づく。
【００６８】
このように、フィラーモデル６近傍の領域Ｔでは、該領域Ｔにおいて、ボンド４がフィラーモデル６の同心円状に配向するため、配向パラメータの平均Χ（ｍ，ｔ）が１に近づく。一方、フィラーモデル６から大きく離間する領域Ｔでは、該領域Ｔでの変形、及び、
配向パラメータχi（ｍ，ｔ）が、フィラーモデル６近傍の領域Ｔに比べて乱雑になるため（即ち、配向パラメータχi（ｍ，ｔ）が−１〜＋１の範囲でランダムな値となる）、配向パラメータの平均Χ（ｍ，ｔ）が０に近づく。従って、配向パラメータの平均Χ（ｍ，ｔ）は、高分子材料モデル１０の変形に相関がある。
【００６９】
次に、配向パラメータの平均Χ（ｍ，ｔ）が、歪に対して同位相で応答する貯蔵成分と、歪に対して異なる位相で応答する散逸成分とに分解される（工程Ｓ８２）。この工程Ｓ８２では、先ず、高分子材料モデル１０の応答が、下記式（９）で計算される。
【００７０】
【数９】
2018112525.tif 000014
ここで、各定数及び変数は、次に示すものを除き、上記式（２）及び（３）と同一である。
Χ0：配向パラメータχの振幅
δ：位相差
【００７１】
上記式（９）は、上記式（２）で定義される各領域Ｔの配向パラメータの平均Χ（ｍ，ｔ）、及び上記式（３）で定義される歪に基づいて、各領域Ｔでの歪の応答が定義される。
【００７２】
次に、上記式（９）で定義される高分子材料モデル１０の応答が、下記式（４）のように分解される。
【００７３】
【数４】
2018112525.tif 000015
ここで、各定数及び変数は、上記式（２）、上記式（３）及び上記式（９）と同一である。
【００７４】
上記式（４）は、上記式（９）を、加法定理に基づいて分解したものである。このような上記式（４）によれば、領域Ｔの歪に対する応答（配向パラメータの平均Χ（ｍ，ｔ））を、歪に対して同位相で応答する貯蔵成分Χ’と、歪に対して異なる位相で応答する散逸成分Χ”とに分解することができる。また、上記式（４）では、貯蔵成分Χ’及び散逸成分Χ”を、時刻ｔから分離することができるため、貯蔵成分Χ’及び散逸成分Χ”を一意に決定することができる。図１４には、各領域Ｔの貯蔵成分Χ’及び散逸成分Χ”、並びにフィラーモデル６の外面６ｓから各領域Ｔの境界面Ｔｓまでの距離の関係を示すグラフが示される。
【００７５】
次に、貯蔵成分Χ’及び散逸成分Χ”に基づいて、界面層が決定される（工程Ｓ８３）。上述したように、図１２に示した高分子材料１４の界面層１４ａでは、バルク部分１４ｂに比べて変形しにくい傾向がある。このため、図１０に示した高分子材料モデル１０の界面層では、歪に対して同位相で応答する貯蔵成分Χ’が、歪に対して異なる位相で応答する散逸成分Χ”よりも大きくなると仮定できる。従って、工程Ｓ８３では、貯蔵成分Χ’が散逸成分Χ”よりも大きい領域Ｔが、界面層１７（図１５に示す）として決定される。
【００７６】
また、高分子材料モデル１０のバルク部分では、貯蔵成分Χ’が散逸成分Χ”よりも小さくなると仮定できる。従って、工程Ｓ８３では、貯蔵成分Χ’が散逸成分Χ”よりも小さい領域Ｔが、バルク部分１８（図１５に示す）として決定される。
【００７７】
さらに、界面層とバルク部分との境界は、貯蔵成分Χ’と散逸成分Χ”とが同一になると仮定できる。従って、工程Ｓ８３では、貯蔵成分Χ’と散逸成分Χ”とが交わる交点１９が、界面層１７とバルク部分１８との境界２０（図１５に示す）の位置として決定される。
【００７８】
なお、境界２０が隣り合う境界面Ｔｓ、Ｔｓ（図１０に示す）の間に配置される場合は、境界２０を基準として、界面層１７及びバルク部分１８が特定されるのが望ましい。また、工程Ｓ８３では、交点１９のフィラーモデル６の外面６ｓからの距離が計算されることにより、界面層１７の厚さＷ３（図１５に示す）が求められる。
【００７９】
このように、本発明のシミュレーション方法では、従来のシミュレーション方法のように、高分子材料及びフィラーを用いた実験を実施することなく、界面層１７を特定し、さらにはその厚さＷ３を確実に計算することができる。従って、本発明のシミュレーション方法では、多くの時間や費用を必要とすることなく、界面層１７を正確に特定することができる。
【００８０】
しかも、本発明のシミュレーション方法では、ポリマーモデル２を用いたシミュレーションにおいて、界面層１７を特定することができるため、現実に存在しない高分子鎖の界面層１７を特定でき、未知の高分子材料の開発に役立つ。
【００８１】
さらに、本実施形態のポリマーモデル２は、全原子モデルとしてモデル化されているため、粗視化モデルに比べて、現実の高分子材料の変形に近似させることができる。従って、本実施形態のシミュレーション方法では、界面層１７を精度よく特定することができる。
【００８２】
界面層１７をより精度良く特定するために、各領域Ｔには、ボンド４が１００個以上配置され、かつ、各領域Ｔの厚さｄｗが単位長さ１σ以下の範囲で決定されるのが望ましい。各領域Ｔに配置されるボンド４が１００個未満であると、各領域Ｔの配向パラメータχi（ｍ，ｔ）を十分に均せないおそれがある。また、厚さｄｗが単位長さ１σを超えると、分解能が低下するおそれがある。
【００８３】
次に、コンピュータ１が、高分子材料モデル１０の物性値を計算する（工程Ｓ９）。この工程Ｓ９では、空間８に、界面層１７の厚さＷ３等の所定のパラメータが設定され、高分子材料の物性値（例えば、複素弾性率）が計算される。
【００８４】
次に、コンピュータ１が、高分子材料モデル１０の物性値が、許容範囲内であるかを判断する（工程Ｓ１０）。この工程Ｓ１０では、物性値が許容範囲内であると判断された場合、ポリマーモデル２を含む空間８の条件等に基づいて、高分子材料が製造される（工程Ｓ１１）。一方、物性値が許容範囲内でないと判断された場合は、ポリマーモデル２を含む空間８の諸条件を変更して（工程Ｓ１２）、工程Ｓ４〜Ｓ１１が再度行われる。このように、本実施形態のシミュレーション方法では、高分子材料モデル１０の物性値が許容範囲内になるまで、ポリマーモデル２を含む空間８の諸条件が変更されるため、所望の性能を有する高分子材料を、効率よく設計することができる。
【００８５】
図４に示されるように、本実施形態のポリマーモデル２は、全原子モデルとして構成されるものが例示されたが、これに限定されるわけではない。図１６に示されるように、例えば、炭素原子に結合した水素原子を、該炭素原子と一体化して一つの粒子モデル３として扱うユナイテッドアトムモデル（ united atom model ）として構成されてもよい。このようなポリマーモデル２は、高分子鎖のモノマーの配置を維持しつつ、水素原子を省略することができる。従って、このようなポリマーモデル２が用いられたシミュレーション方法では、界面層１７の特定精度を維持しつつ、計算時間を短縮することができる。なお、各ポテンシャルのパラメータは、公開されている文献等に基づいて設定されるのが望ましい。
【００８６】
以上、本発明の特に好ましい実施形態について詳述したが、本発明は図示の実施形態に限定されることなく、種々の態様に変形して実施しうる。
【実施例】
【００８７】
図３、図９及び図１３に示される手順に従って、全原子モデルのボンドの配向パラメータが求められ、高分子材料の界面層が特定された（実施例）。そして、実施例において、界面層の厚さを求めるのに必要な時間が測定された。なお、ポテンシャル、及び歪の各パラメータ等は、明細書中の記載通りであり、共通仕様は次のとおりである。
実施例：
ポリマーモデル：
構造：cis-１，４ポリイソプレン
重合度：２０２
慣性半径：１．６nm
フィラーモデルの半径：３．０nm
空間：
１辺の長さＬ１：１０．３nm
ポリマーモデルの本数：３９本
フィラーモデルの個数：１個
各領域：
厚さｄｗ：０．５０nm
個数：１０個
各領域に配置されるボンドの個数：１２１８〜３９３１６個
【００８８】
テストの結果、実施例で計算された界面層の厚さは、４．２nmであった。従って、実施例では、高分子材料を用いた実験等を実施することなく、界面層の厚さを特定することができた。
【符号の説明】
【００８９】
２      ポリマーモデル
４      ボンド
６      フィラーモデル
１０      高分子材料モデル
【図１】
2018112525.tif 000020
【図２】
2018112525.tif 000021
【図３】
2018112525.tif 000022
【図４】
2018112525.tif 000023
【図５】
2018112525.tif 000024
【図６】
2018112525.tif 000025
【図７】
2018112525.tif 000026
【図８】
2018112525.tif 000027
【図９】
2018112525.tif 000028
【図１０】
2018112525.tif 000029
【図１１】
2018112525.tif 000030
【図１２】
2018112525.tif 000031
【図１３】
2018112525.tif 000032
【図１４】
2018112525.tif 000033
【図１５】
2018112525.tif 000034
【図１６】
2018112525.tif 000035