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

Abstract

【課題】 フィラーと高分子材料とが共存する系において、フィラーの周囲に形成される界面層を特定する。
【解決手段】 コンピュータ１を用いて、高分子材料とフィラーとの反応を解析するための方法である。この方法では、粗視化モデル２とフィラーモデル６とを予め定められた空間８内に配置して高分子材料モデル１０を設定する工程Ｓ３、予め定めた条件に基づいて、高分子材料モデル１０の構造緩和を計算する工程Ｓ４、構造緩和の計算後、高分子材料モデル１０を予め定められた歪で変形させて、各ボンド４の配列状態を表す配向パラメータを求める計算工程Ｓ６、配向パラメータに基づいて、フィラーの周囲に形成され高分子材料の界面層を特定する界面特定工程Ｓ８を含む。また、粗視化粒子３、３間、及び粗視化粒子３とフィラー粒子７との間には、粒子間の距離が小さくなるほど無限に大きくなる斥力を含むポテンシャルが定義される。
【選択図】図３

Description

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

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

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

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

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

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

本発明のうち請求項１記載の発明は、コンピュータを用いて、高分子材料とフィラーとの反応を解析するための方法であって、前記コンピュータに、前記高分子材料の高分子鎖を、複数の粗視化粒子と、前記粗視化粒子間を結合するボンドとを用いてモデル化した粗視化モデルを設定する工程、前記コンピュータに、前記フィラーをフィラー粒子でモデル化したフィラーモデルを設定する工程、前記コンピュータに、前記粗視化モデルと前記フィラーモデルとを予め定められた空間内に配置して高分子材料モデルを設定する工程、 前記コンピュータが、予め定めた条件に基づいて、前記高分子材料モデルの構造緩和を計算する工程、前記構造緩和の計算後、前記コンピュータが、前記高分子材料モデルを予め定められた歪で変形させて、前記各ボンドの配列状態を表す配向パラメータを求める計算工程、前記配向パラメータに基づいて、前記コンピュータが、前記フィラーの周囲に形成されかつ前記高分子材料のバルク部分とは異なる力学的特性を示す前記高分子材料の界面層を特定する界面特定工程を含み、前記粗視化粒子間、及び前記粗視化粒子と前記フィラー粒子との間には、前記粒子間の距離が小さくなるほど無限に大きくなる斥力を含むポテンシャルが定義されることを特徴とする。

また、請求項２記載の発明は、前記計算工程は、前記粗視化モデルが配置された前記空間を、前記フィラーモデルの外面に沿った境界面で複数の領域に区分する工程と、前記高分子材料モデルを前記歪で変形させる工程と、前記領域毎に、前記配向パラメータを計算する工程とを含む請求項１に記載の高分子材料のシミュレーション方法である。

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

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

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

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



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

また、請求項６記載の発明は、前記各領域には、前記ボンドが１００個以上配置されている請求項２乃至５のいずれかに記載の高分子材料のシミュレーション方法である。

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

ここで、各定数及び変数は、次のとおりである。
ｒ_ij：粒子間の距離
ｒ_c：カットオフ距離
ε：子間に定義されるＬＪポテンシャルの強度
σ：粒子の直径に相当

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

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

さらに、本発明では、粗視化粒子間、及び粗視化粒子とフィラー粒子との間に、粒子間の距離が小さくなるほど無限に大きくなる斥力を含むポテンシャルが定義される。このようなポテンシャルは、高分子材料の分子運動に基づくものである。従って、本発明では、実際の高分子材料の分子運動に近似させて、高分子材料の界面層を精度よく特定することができる。

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

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

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

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

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

図２及び図４に示されるように、高分子材料の高分子鎖がポリイソプレンである場合には、例えば１．７３個分のモノマーを構造単位（図示省略）として、該構造単位を１個の粗視化粒子３に置換される。これにより、各粗視化モデル２には、複数（例えば、１０〜５０００個）の粗視化粒子３が設定される。

なお、１．７３個分のモノマーを構造単位としたのは、論文（ 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）の記載に基づき求めたものである。また、高分子鎖がポリイソプレン以外の場合でも、例えば、上記論文と論文（L,J.Fetters ,D.J.Lohse and R.H.Colby 著、「Chain Dimension and Entanglement Spacings 」Physical Properties of Polymers Handbook Second Edition P448）から構造単位を設定することができる。

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

ボンド４には、相互ポテンシャルＰ１が設定される。本実施形態の相互ポテンシャルＰ１は、下記式（５）で定義される斥力を含むポテンシャル（以下、「ＬＪポテンシャルＵ_LJ（ｒ_ij）」ということがある。）と、下記式（６）で定義される結合ポテンシャルＵ_FINEとの和で定義される。

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

上記式（５）において、粗視化粒子３、３間の距離ｒ_ijが小さくなると、斥力が作用するＬＪポテンシャルＵ_LJ（ｒ_ij）が大きくなる。一方、上記式（６）において、粗視化粒子３、３間の距離ｒ_ijが大きくなると、引力が作用する結合ポテンシャルＵ_FINEが大きくなる。従って、相互ポテンシャルＰ１は、距離ｒ_ijを、ＬＪポテンシャルＵ_LJ（ｒ_ij）と結合ポテンシャルＵ_FINEとが互いに釣り合う位置に戻そうとする復元力が定義される。

また、上記式（５）では、粗視化粒子３、３間の距離ｒ_ijが小さくなるほど、ＬＪポテンシャルＵ_LJ（ｒ_ij）が無限に大きくなる。一方、上記式（６）では、距離ｒ_ijが伸びきり長Ｒ₀以上となる場合に、結合ポテンシャルＵ_FINEが∞に設定される。従って、相互ポテンシャルＰ１は、粗視化粒子３、３間の重なりと、伸びきり長Ｒ₀以上の距離ｒ_ijを許容しない。

このような相互ポテンシャルＰ１は、分子動力学計算において、粗視化モデル２を、高分子材料の分子運動に近似させることができる。なお、ＬＪポテンシャルＵ_LJ（ｒ_ij）及び結合ポテンシャルＵ_FINEの各定数及び変数の値としては、適宜設定することができる。本実施形態では、論文（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）に基づいて、次の値が設定される。
ばね定数ｋ：３０
伸びきり長Ｒ₀：１．５
定数ε：１．０
定数σ：１．０
カットオフ距離ｒ_c：２．５

また、相互ポテンシャルＰ１としては、上記のようなＬＪポテンシャルＵ_LJ（ｒ_ij）と、結合ポテンシャルＵ_FINEとの和で定義されるものの他、下記式（７）で定義される結合ポテンシャルＵ_Harmonicでもよい。

ここで、各定数及び変数は、次に示すものを除いて、上記式（５）、（６）と同一である。
ｒ：粗視化粒子間の距離
ｒ₀：平衡長
なお、平衡長ｒ₀は、各粗視化粒子３、３の中心３ｃ、３ｃ間の距離として定義される。

上記式（７）で定義される結合ポテンシャルＵ_Harmonicは、Harmonic型である。Harmonic型とは、いわゆる線形バネが定義され、平衡長ｒ₀からのずれの大きさに比例した復元力が働くポテンシャルである。

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

また、結合ポテンシャルＵ_Harmonicは、ＬＪポテンシャルＵ_LJ（ｒ_ij）と結合ポテンシャルＵ_FINEとの和とは異なり、距離ｒ_ijの上限値が、伸びきり長Ｒ₀に限定されない。従って、結合ポテンシャルＵ_Harmonicは、伸びきり長Ｒ₀以上の距離ｒを許容することができる。このような結合ポテンシャルＵ_Harmonicは、分子動力学計算において、応力集中などによる分子鎖の伸長や、粗視化モデル２、２同士の絡まりによって作用する大きな力を抑制することができ、計算落ちの発生を防ぐことができる。

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

図５に示されるように、隣接する粗視化モデル２、２の粗視化粒子３、３間には、相互ポテンシャルＰ２が定義される。相互ポテンシャルＰ２は、上記式（５）のＬＪポテンシャルＵ_LJ（ｒ_ij）によって定義される。なお、相互ポテンシャルＰ２の強度ε及び定数σも、適宜設定することができるが、上記式（５）と同一範囲が望ましい。

このような粗視化モデル２を用いたシミュレーションでは、例えば全原子モデルを用いたシミュレーションに比べて、計算時間を大幅に短縮でき、大規模かつ長時間の現象を扱うことができる。なお、これらの粗視化モデル２は、コンピュータ１で取り扱い可能な数値データであり、コンピュータ１に記憶される。

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

また、フィラー粒子７、７間、及び、粗視化粒子３とフィラー粒子７との間には、相互ポテンシャルＰ３が設定される。相互ポテンシャルＰ３は、上記式（５）のＬＪポテンシャルＵ_LJ（ｒ_ij）によって定義される。また、相互ポテンシャルＰ３の各定数及び変数の値としては、適宜設定することができる。フィラーモデル６は、コンピュータ１で取り扱い可能な数値データであり、コンピュータ１に記憶される。

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

空間８の一辺の長さＬ１は、適宜設定することができるが、粗視化モデル２の慣性半径Ｒｇの２倍以上に設定されるのが望ましい。慣性半径Ｒｇは、粗視化モデル２の拡がりを示すパラメータである。このような空間８では、分子動力学計算において、粗視化モデル２のダイナミクス（動き）をスムーズに計算することができる。さらに、空間８の大きさは、例えば、上記論文に基づいて、粒子数密度が０．８５個／σ³程度に設定されるのが望ましい。

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

次に、コンピュータ１が、予め定めた条件に基づいて、高分子材料モデル１０の構造緩和を計算する（工程Ｓ４）。本実施形態の分子動力学計算では、例えば、空間８について所定の時間、粗視化モデル２及びフィラーモデル６が古典力学に従うものとして、ランジュバンの運動方程式が適用される。そして、各時刻での粗視化粒子３及びフィラー粒子７の動きが、単位時間毎に追跡される。

また、本実施形態の構造緩和の計算は、空間８において、圧力及び温度が一定、又は体積及び温度が一定に保たれる。これにより、工程Ｓ４では、実際の高分子材料の分子運動に近似させて、粗視化モデル２及びフィラー粒子７の初期配置を精度よく緩和することができる。このような構造緩和の計算は、例えばソフトマテリアル総合シミュレーター（ＯＣＴＡ）に含まれるＣＯＧＮＡＣを用いて処理することができる。

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

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

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

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

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

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

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

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

一方、フィラーモデル６から大きく離間する領域Ｔ（例えば、第６領域Ｔ６〜第Ｍ領域ＴＭ）では、粗視化粒子３とフィラー粒子７との間に作用する相互ポテンシャルＰ３の引力が小さくなるため、上記変形が定義されると、粗視化粒子３の移動が大きくなり、該領域Ｔでの変形が大きくなる。このように、上記歪で変形させられた高分子材料モデル１０は、フィラーが配合された高分子材料の変形を、精度良く再現することができる。

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

上記式（３）の各定数及び変数については、適宜設定することができるが、例えば、下記の値が設定されるのが望ましい。さらに、高分子材料モデル１０の変形計算は、例えば、からみ合い点の緩和時間をτ_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が乱雑になり（即ち、拘束力が弱いため、角度θ_iがランダムな方向に向く）、配向パラメータχ_i（ｍ，ｔ）が−１〜＋１の範囲で様々な値をとるため、配向パラメータの平均Χ（ｍ，ｔ）は０に近づく。

このように、フィラーモデル６近傍の領域Ｔでは、該領域Ｔにおいてボンド４がフィラーモデル６の同心円状に配向するため、配向パラメータχ_i（ｍ，ｔ）の平均Χ（ｍ，ｔ）が１に近づく。一方、フィラーモデル６から大きく離間する領域Ｔでは、該領域Ｔでの変形、及び、配向パラメータχ_i（ｍ，ｔ）が、フィラーモデル６近傍の領域Ｔに比べて乱雑になるため（即ち、配向パラメータχ_i（ｍ，ｔ）が−１〜＋１の範囲で様々な値となる）、配向パラメータの平均Χ（ｍ，ｔ）が０に近づく。従って、配向パラメータχ_i（ｍ，ｔ）は、高分子材料モデル１０の変形に相関がある。

次に、配向パラメータの平均Χ（ｍ，ｔ）が、歪に対して同位相で応答する貯蔵成分と、歪に対して異なる位相で応答する散逸成分とに分解される（工程Ｓ８２）。この工程Ｓ８２では、先ず、各領域Ｔでの歪の応答が、下記式（８）で計算される。

ここで、各定数及び変数は、次に示すものを除き、上記式（２）及び（３）と同一である。
Χ₀：配向パラメータχの振幅
δ：位相差

上記式（８）は、上記式（２）で定義される各領域Ｔの配向パラメータの平均Χ（ｍ，ｔ）、及び上記式（３）で定義される歪に基づいて、各領域Ｔでの歪の応答が定義される。

次に、上記式（８）で定義される高分子材料モデル１０の応答が、下記式（４）のように分解される。

ここで、各定数及び変数は、上記式（２）、上記式（３）及び上記式（８）と同一である。

上記式（４）は、上記式（８）を、加法定理に基づいて分解したものである。このような上記式（４）によれば、高分子材料モデル１０の応答（配向パラメータの平均Χ（ｍ，ｔ））を、歪に対して同位相で応答する貯蔵成分Χ’と、歪に対して異なる位相で応答する散逸成分Χ”とに分解することができる。また、上記式（４）では、貯蔵成分Χ’及び散逸成分Χ”を、時刻ｔから分離することができるため、貯蔵成分Χ’及び散逸成分Χ”を一意に決定することができる。図１２には、各領域Ｔの貯蔵成分Χ’及び散逸成分Χ”、並びにフィラーモデル６の外面６ｓから各領域Ｔの境界面Ｔｓまでの距離の関係を示すグラフが示される。

次に、貯蔵成分Χ’及び散逸成分Χ”に基づいて、界面層が決定される（工程Ｓ８３）。上述したように、図１０に示した高分子材料１４の界面層１４ａでは、バルク部分１４ｂに比べて変形しにくい傾向がある。このため、図８に示した高分子材料モデル１０の界面層では、歪に対して同位相で応答する貯蔵成分Χ’が、歪に対して異なる位相で応答する散逸成分Χ”よりも大きくなると仮定することができる。このような観点より、工程Ｓ８３では、貯蔵成分Χ’が散逸成分Χ”よりも大きい領域Ｔが、界面層１７（図１３に示す）として決定される。

また、高分子材料モデル１０のバルク部分では、貯蔵成分Χ’が散逸成分Χ”よりも小さくなると仮定することができる。従って、工程Ｓ８３では、貯蔵成分Χ’が散逸成分Χ”よりも小さい領域Ｔが、バルク部分１８（図１３に示す）として決定される。

さらに、界面層とバルク部分との境界は、貯蔵成分Χ’と散逸成分Χ”とが同一になると仮定することができる。従って、貯蔵成分Χ’と散逸成分Χ”とが交わる交点１９が、界面層１７とバルク部分１８との境界２０（図１３に示す）の位置として決定される。

なお、境界２０が隣り合う境界面Ｔｓ、Ｔｓ（図８に示す）の間に配置される場合は、境界２０を基準として、界面層１７及びバルク部分１８が区分されるのが望ましい。また、交点１９のフィラーモデルの外面からの距離を計算することにより、界面層１７の厚さＷ３（図１３に示す）が求められる。

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

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

さらに、本発明のシミュレーション方法では、粗視化粒子３、３間、及び粗視化粒子３とフィラー粒子７との間に、上記のような斥力を含むポテンシャル（ＬＪポテンシャルＵ_LJ（ｒ_ij））が定義されるため、実際の高分子材料の分子運動に近似させて、シミュレーションすることができる。従って、本実施形態のシミュレーション方法では、界面層１７を精度よく特定することができる。

界面層１７を精度良く特定するために、ボンド４が１００以上配置され、かつ、各領域Ｔの厚さｄｗが単位長さσ以下の範囲で決定されるのが望ましい。各領域Ｔに配置されるボンド４が１００個未満であると、各領域Ｔの配向パラメータχ_i（ｍ，ｔ）を十分に均せないおそれがある。また、厚さｄｗが単位長さ１σを超えると、分解能が低下するおそれがある。

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

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

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

図３、図７及び図１１に示される手順に従って、高分子材料モデルのボンドの配向パラメータが求められ、高分子材料の界面層が特定された（実施例）。

そして、実施例において、界面層の厚さを求めるのに要した時間が測定された。なお、ポテンシャル、及び歪の各パラメータ等は、明細書中の記載通りであり、共通仕様は次のとおりである。
粗視化モデル：
粗視化粒子の個数：１０００個／本
慣性半径：１６．９σ
フィラーモデルの半径：１１．２σ
各領域：
厚さｄｗ：１．０σ
個数：１０個
各領域に配置されるボンドの個数：７３０〜２０７３個
空間：
１辺の長さＬ１：６１．１σ
粗視化モデルの本数：１８９本
フィラーモデルの個数：１個

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

２ 粗視化モデル
４ ボンド
６ フィラーモデル
１０ 高分子材料モデル

Claims (7)

1. コンピュータを用いて、高分子材料とフィラーとの反応を解析するための方法であって、
   前記コンピュータに、前記高分子材料の高分子鎖を、複数の粗視化粒子と、前記粗視化粒子間を結合するボンドとを用いてモデル化した粗視化モデルを設定する工程、
   前記コンピュータに、前記フィラーをフィラー粒子でモデル化したフィラーモデルを設定する工程、
   前記コンピュータに、前記粗視化モデルと前記フィラーモデルとを予め定められた空間内に配置して高分子材料モデルを設定する工程、
   前記コンピュータが、予め定めた条件に基づいて、前記高分子材料モデルの構造緩和を計算する工程、
   前記構造緩和の計算後、前記コンピュータが、前記高分子材料モデルを予め定められた歪で変形させて、前記各ボンドの配列状態を表す配向パラメータを求める計算工程、
   前記配向パラメータに基づいて、前記コンピュータが、前記フィラーの周囲に形成されかつ前記高分子材料のバルク部分とは異なる力学的特性を示す前記高分子材料の界面層を特定する界面特定工程を含み、
   前記粗視化粒子間、及び前記粗視化粒子と前記フィラー粒子との間には、前記粒子間の距離が小さくなるほど無限に大きくなる斥力を含むポテンシャルが定義されることを特徴とする高分子材料のシミュレーション方法。
2. 前記計算工程は、前記粗視化モデルが配置された前記空間を、前記フィラーモデルの外面に沿った境界面で複数の領域に区分する工程と、
   前記高分子材料モデルを前記歪で変形させる工程と、
   前記領域毎に、前記配向パラメータを計算する工程とを含む請求項１に記載の高分子材料のシミュレーション方法。
3. 前記界面特定工程は、前記領域毎に前記配向パラメータの平均を求める工程と、
   前記配向パラメータの平均を、前記歪に対して同位相で応答する貯蔵成分と、前記歪に対して異なる位相で応答する散逸成分とに分解する工程と、
   前記貯蔵成分が前記散逸成分よりも大きい領域を前記界面層と決定する工程とを含む請求項２に記載の高分子材料のシミュレーション方法。
4. 前記配向パラメータχ_i（ｍ，ｔ）は、前記ボンドに沿ったボンドベクトルｕ_i（ｍ，ｔ）と、前記ボンドの重心と前記フィラーモデルの重心とを結ぶ重心ベクトルｎ_i（ｍ，ｔ）とに基づいて、下記式（１）で定義される請求項１乃至３のいずれかに記載の高分子材料のシミュレーション方法。
   
   ここで、各定数及び変数は次のとおりである。
   ｉ：各ボンドのインデックス
   ｍ：各領域のインデックス
   ｔ：時刻
   θ_i：ボンドベクトルｕ_i（ｍ，ｔ）と重心ベクトルｎ_i（ｍ，ｔ）とがなす角度
5. 前記配向パラメータの平均Χ（ｍ，ｔ）は、下記式（２）で定義され、
   前記歪γは、下記式（３）で定義され、
   前記配向パラメータの平均Χ（ｍ，ｔ）は、下記式（４）に基づいて、前記貯蔵成分Χ’と前記散逸成分Χ”とに分解される請求項４に記載の高分子材料のシミュレーション方法。
   
   
   
   ここで、各定数及び変数は次のとおりである。
   ｉ：各ボンドのインデックス
   ｍ：各領域のインデックス
   ｔ：時刻
   Ｎ：各領域に配置されるボンドの総数
   γ₀：歪振幅
   ω：角速度
6. 前記各領域には、前記ボンドが１００個以上配置されている請求項２乃至５のいずれかに記載の高分子材料のシミュレーション方法。
7. 前記斥力を含むポテンシャルは、下記式（５）で定義される請求項１乃至６のいずれかに記載の高分子材料のシミュレーション方法。
   
   ここで、各定数及び変数は、次のとおりである。
   ｒ_ij：粒子間の距離
   ｒ_c：カットオフ距離
   ε：粒子間に定義されるＬＪポテンシャルの強度
   σ：粒子の直径に相当