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