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

Abstract

【課題】 高分子成分とフィラーとの相互作用を精度良く解析する。【解決手段】 高分子材料に含まれる高分子成分とフィラーとの相互作用を、コンピュータを用いて解析するための方法である。このシミュレーション方法は、コンピュータが、ポリマー粒子モデル７、及び、フィラー粒子モデル１４を対象として分子動力学計算を行う工程とを含んでいる。フィラーモデル１３は、第１フィラーモデル１３Ａと第２フィラーモデル１３Ｂとを含んでいる。第１フィラーモデル１３Ａは、複数のフィラー粒子モデル１４が、見かけ上の球面を構成するように配置されかつ互いに連結されている。第２フィラーモデル１３Ｂは、第１フィラーモデル１３Ａのフィラー粒子モデル１４の少なくとも一つが欠損したものである。第２フィラーモデル１３Ｂは、セル４内に定義された分子鎖モデル６と相互作用するように、セル４の中心側に位置して第１フィラーモデル１３Ａと重ねられている。【選択図】図６

Description

本発明は、高分子材料のシミュレーション方法に関し、詳しくは、高分子材料に含まれる高分子成分とフィラーとの相互作用を、コンピュータを用いて解析するための方法に関する。

近年、フィラーが配合された高分子材料の開発のために、コンピュータを用いたシミュレーションが行われている。下記特許文献１のシミュレーション方法では、フィラーをモデル化したフィラーモデル、及び、高分子成分の分子鎖をモデル化した分子鎖モデルを用いた分子動力学計算が行われ、両者の反応状態等を推定する試みがなされている。

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

上記特許文献１のフィラーモデルは、仮想空間であるセルを構成する一対の面に、均一な平面を構成するように配置された複数のフィラー粒子モデルで定義されている。このようなフィラーモデルを用いたシミュレーション方法では、例えば、フィラーの表面状態等に起因するフィラーと分子鎖との相互作用の変化を考慮することができず、高分子成分とフィラーとの相互作用を精度良く解析できないという問題があった。

本発明は、以上のような実状に鑑み案出されたもので、高分子成分とフィラーとの相互作用を精度良く解析することができる高分子材料のシミュレーション方法を提供することを主たる目的としている。

本発明は、高分子材料に含まれる高分子成分とフィラーとの相互作用を、コンピュータを用いて解析するための方法であって、前記高分子材料の一部に対応する仮想空間であるセルを、前記コンピュータに入力する工程と、前記セル内に、複数のポリマー粒子モデルを用いて前記高分子成分の分子鎖をモデル化した分子鎖モデルを定義する工程と、前記セル内に、複数のフィラー粒子モデルを用いて前記フィラーをモデル化したフィラーモデルを定義する工程と、前記ポリマー粒子モデル間、及び、前記フィラー粒子モデルと前記ポリマー粒子モデルとの間に、引力及び斥力が作用するポテンシャルを定義する工程と、前記コンピュータが、前記ポリマー粒子モデル、及び、前記フィラー粒子モデルを対象として分子動力学計算を行う工程とを含み、前記フィラーモデルは、第１フィラーモデルと第２フィラーモデルとを含んで定義され、前記第１フィラーモデルは、複数の前記フィラー粒子モデルが、見かけ上の球面を構成するように配置されかつ互いに連結されており、前記第２フィラーモデルは、前記第１フィラーモデルの前記フィラー粒子モデルの少なくとも一つが欠損したものであり、前記第２フィラーモデルは、前記セル内に定義された前記分子鎖モデルと相互作用するように、前記セルの中心側に位置して前記第１フィラーモデルと重ねられていることを特徴とする。

本発明に係る前記高分子材料のシミュレーション方法において、前記第１フィラーモデルの前記フィラー粒子モデル、及び、前記第２フィラーモデルの前記フィラー粒子モデルは、それらの半径方向で互いに結合されていてもよい。

本発明に係る前記高分子材料のシミュレーション方法において、前記高分子材料は、前記フィラーに対する前記分子鎖の親和性を高める官能基を含み、前記官能基を少なくとも一つの官能基粒子モデルでモデル化した官能基モデルを、前記フィラーモデルの前記フィラー粒子モデルに結合させる工程と、前記官能基粒子モデルと前記ポリマー粒子モデルとの間に、前記ポテンシャルを定義する工程とをさらに含み、前記分子動力学計算を行う工程は、前記ポリマー粒子モデル、前記フィラー粒子モデル、及び、前記官能基粒子モデルを対象として、前記分子動力学計算を行ってもよい。

本発明に係る前記高分子材料のシミュレーション方法において、前記官能基粒子モデルは、前記第２フィラーモデルの前記欠損した部分に配置されてもよい。

本発明の高分子材料のシミュレーション方法の前記フィラーモデルは、前記第１フィラーモデルと前記第２フィラーモデルとを含んで定義される。前記第１フィラーモデルは、複数の前記フィラー粒子モデルが、見かけ上の球面を構成するように配置されかつ互いに連結されている。前記第２フィラーモデルは、前記第１フィラーモデルの前記フィラー粒子モデルの少なくとも一つが欠損したものである。前記第２フィラーモデルは、前記セル内に定義された前記分子鎖モデルと相互作用するように、前記セルの中心側に位置して前記第１フィラーモデルと重ねられている。

発明者らは、鋭意研究を重ねた結果、フィラーの表面に、フィラーの原子に置き換わる官能基が存在し、フィラーと分子鎖との相互作用に変化が生じることを知見した。本実施形態のシミュレーション方法では、前記第２フィラーモデルの前記フィラー粒子モデルが欠損した部分と、欠損していない部分とにより、前記分子鎖モデルに対する前記ポテンシャルを変化させることができる。従って、本発明のシミュレーション方法は、前記フィラーの表面状態等に起因する前記フィラーと前記分子鎖との相互作用の変化を考慮して、前記高分子成分と前記フィラーとの相互作用を精度良く解析することが可能となる。

また、本実施形態の第１フィラーモデル及び第２フィラーモデルは、見かけ上の球面を構成するように前記フィラー粒子モデルが配置されているため、現実のフィラーの表面形状に近似させることができる。従って、本発明の高分子材料のシミュレーション方法は、フィラーモデル及び分子鎖モデルを用いた分子動力学計算が実施されることにより、現実のフィラー近傍での分子鎖の分布を精度良く再現することができる。

高分子材料のシミュレーション方法を実行するためのコンピュータの一例を示す斜視図である。 ポリイソプレンの構造式である。 高分子材料のシミュレーション方法の処理手順の一例を示すフローチャートである。 セル、分子鎖モデル及びフィラーモデルの一例を示す概念図である。 分子鎖モデルの一例を示す概念図である。 フィラーモデルの一例を示す概念図である。 第１フィラーモデルの一例を示す部分斜視図である。 第２フィラーモデルの一例を示す部分斜視図である。 ポテンシャルの一例を示す概念図である。 本発明の他の実施形態の高分子材料のシミュレーション方法の処理手順を示すフローチャートである。 本発明の他の実施形態の第２フィラーモデル及び官能基モデルを示す部分斜視図である。

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

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

高分子材料は、少なくとも１種類、本実施形態では１種類の高分子成分を含んでいる。高分子成分の一例は、天然ゴム、イソプレンゴム、ブタジエンゴム、又は、スチレンブタジエンゴム等である。本実施形態の高分子成分としては、イソプレンゴム（cis-１，４ポリイソプレン（以下、単に「ポリイソプレン」ということがある。））である場合が例示される。図２は、ポリイソプレンの構造式である。

ポリイソプレンを構成する分子鎖２は、メチン基等（例えば、−ＣＨ＝、＞Ｃ＝）、メチレン基（−ＣＨ_２−）、及び、メチル基（−ＣＨ_３）によって構成されるイソプレンのモノマー（イソプレン分子）３が、重合度ｎで連結されて構成されている。なお、高分子材料には、ポリイソプレン以外の高分子材料が用いられてもよい。

フィラーの一例としては、グラファイト（カーボンブラック）、又は、シリカ等である。本実施形態のフィラーとしては、グラファイトである場合が例示される。また、高分子鎖には、フィラーに対する分子鎖の親和性を高める官能基が含まれている。官能基の一例としては、アルデヒド、ケトン、カルボン酸、又は、カルボン酸塩である。本実施形態の官能基としては、アルデヒド、及び、ケトンである場合が例示される。

図３は、シミュレーション方法の処理手順の一例を示すフローチャートである。本実施形態のシミュレーション方法では、先ず、高分子材料の一部に対応する仮想空間であるセル４が、コンピュータ１に入力される（工程Ｓ１）。図４は、セル、分子鎖モデル及びフィラーモデルの一例を示す概念図である。

セル４は、少なくとも互いに向き合う一対の面５、５、本実施形態では、互いに向き合う三対の面５、５を有しており、直方体又は立方体（本実施形態では、立方体）として定義されている。各面５、５には、周期境界条件が定義されている。このようなセル４が用いられることにより、後述の分子動力学計算において、例えば、高分子成分の分子鎖２（図２に示す）をモデル化した後述の分子鎖モデル６について、一方側の面５ａから出て行った分子鎖モデル６の一部が、他方側の面５ｂから入ってくるように計算することができる。従って、一方側の面５ａと、他方側の面５ｂとが連続している（繋がっている）ものとして取り扱うことができる。

セル４の一辺の各長さＬ１（Ｌ１ａ、Ｌ１ｂ及びＬ１ｃ）は、適宜設定することができる。本実施形態の長さＬ１は、後述の分子鎖モデル６の拡がりを示す量である慣性半径（図示省略）の３倍以上が望ましい。これにより、セル４は、後述の分子動力学計算において、周期境界条件による自己のイメージとの衝突の発生を防げるため、分子鎖モデル６の空間的拡がりを適切に計算することができる。また、セル４の大きさは、例えば１気圧で安定な体積に設定される。これにより、セル４は、解析対象の高分子材料の少なくとも一部の体積を定義することができる。セル４は、コンピュータ１に記憶される。

次に、本実施形態のシミュレーション方法では、セル４内に、高分子成分の分子鎖をモデル化した分子鎖モデル６が定義される（工程Ｓ２）。本実施形態の分子鎖モデル６は、全原子モデルとして構成されている。

分子鎖モデル６は、複数のポリマー粒子モデル７を用いてモデル化され、セル４内に移動可能に定義されている。移動可能とは、後述の分子動力学計算において、分子鎖モデル６がセル４内を移動することができることを意味している。本実施形態の分子鎖モデル６は、複数のポリマー粒子モデル７と、隣接するポリマー粒子モデル７、７間を結合する結合鎖モデル８とを含んで構成されている。図５は、分子鎖モデルの一例を示す概念図である。

ポリマー粒子モデル７及び結合鎖モデル８は、図２に示した分子鎖２のモノマー３をなす単位構造に基づいて連結される。これにより、モノマーモデル１１が設定される。本実施形態の分子鎖モデル６は、一つのモノマーモデル１１のみで定義されている。これにより、本実施形態のシミュレーション方法は、後述の分子動力学計算において、計算対象のポリマー粒子モデル７の鎖長を大幅に短くすることができるため、構造緩和に要する計算時間を短縮しうる。しかも、分子鎖モデル６は、モノマー３をなす単位構造に基づいて定義されるため、計算精度の低下も抑えることができる。なお、分子鎖モデル６は、計算精度をさらに高めるために、モノマーモデル１１を、分子量（重合度）Ｍnに基づいて連結されることで定義されてもよい。

ポリマー粒子モデル７は、後述する分子動力学計算において、運動方程式の質点として取り扱われる。即ち、ポリマー粒子モデル７は、質量、直径、電荷、又は、初期座標などのパラメータが定義される。本実施形態のポリマー粒子モデル７は、図２に示した分子鎖２の炭素原子をモデル化した炭素粒子モデル７Ｃ、及び、水素原子をモデル化した水素粒子モデル７Ｈを含んでいる。

結合鎖モデル８は、主鎖８ａと側鎖８ｂとを含んでいる。主鎖８ａは、炭素粒子モデル７Ｃ、７Ｃを連結するためのものである。側鎖８ｂは、炭素粒子モデル７Ｃと水素粒子モデル７Ｈとの間を連結するためのものである。また、結合鎖モデル８には、一対の原子間の一重結合を定義した一重結合ボンドモデル８Ａと、一対の原子間の二重結合を定義した二重結合ボンドモデル８Ｂとを含んでいる。

結合鎖モデル８を介して隣り合うポリマー粒子モデル７、７間には、引力及び斥力が作用するポテンシャルＰ１が定義される。一重結合ボンドモデル８Ａ及び二重結合ボンドモデル８Ｂには、異なる値のポテンシャルＰ１がそれぞれ定義される。

ポテンシャルＰ１については、従来と同様に、適宜定義することができる。本実施形態のポテンシャルＰ１は、論文１（Stephen L. Mayo, Barry D. Olafson, and William A. Goddard III 著 「DREIDING: A Generic Force Field for Molecular Simulations」、J. Chem Phys. vol.94, No.26 ,8897-8909, 1990）に基づいて、図２に示した分子鎖２の構造に応じて設定されている。

図４に示したセル４内に定義される分子鎖モデル６の個数については、適宜設定することができる。本実施形態の分子鎖モデル６の個数は、セル４内に配置される全ての分子鎖モデル６の密度が、ポリイソプレンの密度に一致するように調整されるのが望ましい。

本実施形態の分子鎖モデル６は、全原子モデルとして構成されたが、このような態様に限定されない。分子鎖モデル６は、例えば、粗視化モデル（Kremer-Grestモデル等）として定義されてもよい。さらに、分子鎖モデル６は、特開２０１４−２２５２２６号公報の記載に基づいて、緩和された粗視化モデルにモノマーモデル１１を割り当てるリバースマッピングが実施されることによって、全原子モデルが構成されてもよい。分子鎖モデル６は、コンピュータ１に記憶される。

次に、本実施形態のシミュレーション方法では、セル４内に、フィラーをモデル化したフィラーモデルが定義される（工程Ｓ３）。図６は、フィラーモデル１３の一例を示す概念図である。フィラーモデル１３は、複数のフィラー粒子モデル１４を用いてモデル化され、セル４内に移動可能に定義されている。ここで、「移動可能」とは、後述の分子動力学計算において、フィラーモデル１３がセル４内を移動することができることを意味している。なお、フィラーモデル１３は、移動不能に定義されていてもよい。

本実施形態のフィラーモデル１３は、第１フィラーモデル１３Ａと、第２フィラーモデル１３Ｂとを含んで定義されている。図７は、第１フィラーモデル１３Ａの一例を示す部分斜視図である。図７では、第２フィラーモデル１３Ｂを省略して表示している。

図７に示されるように、第１フィラーモデル１３Ａは、複数のフィラー粒子モデル１４が、見かけ上の球面を構成するように配置されており、互いに連結されている。なお、見かけ上の球面としては、隣接するフィラー粒子モデル１４によって形成される凹凸が許容されるものとする。

フィラー粒子モデル１４は、後述の分子動力学計算において、運動方程式の質点として取り扱われる。即ち、フィラー粒子モデル１４には、質量、直径、電荷、又は、初期座標などのパラメータが定義される。フィラー粒子モデル１４の粒子径については、適宜設定することができる。本実施形態のフィラー粒子モデル１４の粒子径は、ポリマー粒子モデル７の粒子径と同一、又は、やや大きく設定されている。

第１フィラーモデル１３Ａは、フィラー粒子モデル１４を連結した複数のユニット１６が、見かけ上の球面を構成するように連結されるのが望ましい。ユニット１６は、例えば、フィラー（本実施形態では、グラファイト）の構造に基づいて、六角形の頂点に配されたフィラー粒子モデル１４が連結されることで定義される。また、隣接するユニット１６、１６のフィラー粒子モデル１４、１４間も同様に、連結されている。これにより、現実のフィラーの構造に基づくフィラーモデル１３を、容易に定義することができる。

隣接するフィラー粒子モデル１４、１４は、結合鎖モデル１５によって連結されている。結合鎖モデル１５を介して隣り合うフィラー粒子モデル１４、１４間には、引力及び斥力が作用するポテンシャルＰ２が定義される。ポテンシャルＰ２については、従来と同様に、適宜定義することができる。本実施形態のポテンシャルＰ２は、上記論文１に基づいて、フィラーの構造に応じて設定されている。

第１フィラーモデル１３Ａを構成するフィラー粒子モデル１４の半径方向内側は、他のフィラー粒子モデル１４が配置されない空間として形成されている。これにより、本実施形態のシミュレーション方法は、例えば、フィラー粒子モデル１４が半径方向内側まで密に配置されたフィラーモデル（図示省略）が用いられるシミュレーション方法に比べて、後述の分子動力学計算での計算対象のフィラー粒子モデル１４の数を大幅に減らすことができるため、計算時間を短縮しうる。

第１フィラーモデル１３Ａにおいて、隣接するフィラー粒子モデル１４、１４間に形成される隙間は、図４に示したポリマー粒子モデル７の粒子径よりも小さく設定されるのが望ましい。これにより、第１フィラーモデル１３Ａは、第１フィラーモデル１３Ａの半径方向内側に、ポリマー粒子モデル７が浸入するのを阻止できるため、後述の分子動力学計算において、分子鎖モデル６がフィラーモデル１３の内部に浸入するような非現実的な計算を防ぐことができる。

図８は、第２フィラーモデル１３Ｂの一例を示す部分斜視図である。なお、図８では、第１フィラーモデル１３Ａを省略して表示している。第２フィラーモデル１３Ｂは、図７に示した第１フィラーモデル１３Ａのフィラー粒子モデル１４の少なくとも一つが欠損したものである。このため、図８に示されるように、第２フィラーモデル１３Ｂは、第１フィラーモデル１３Ａと同様に、複数の前記フィラー粒子モデル１４が、見かけ上の球面を構成するように配置されかつ互いに連結されている。本実施形態の第２フィラーモデル１３Ｂでは、例えば、フィラーに結合している官能基の数や種類に基づいて、フィラー粒子モデル１４を欠損させている。なお、フィラー粒子モデル１４は、ランダムに欠損させてもよい。官能基の数や種類については、例えば、実験等によって求められる。

図６に示されるように、第２フィラーモデル１３Ｂは、第１フィラーモデル１３Ａに対してセル４の中心側に位置するように、第１フィラーモデル１３Ａと重ねられている。本実施形態の第２フィラーモデル１３Ｂは、第１フィラーモデル１３Ａの半径方向外側において、第１フィラーモデル１３Ａよりも大きな外径を有する見かけ上の球面を構成している。これにより、本実施形態のフィラーモデル１３は、第１フィラーモデル１３Ａと第２フィラーモデル１３Ｂとが重なる２層で構成される。フィラーモデル１３は、セル４の対角線上で互いに向き合って配置されている。

第２フィラーモデル１３Ｂは、第１フィラーモデル１３Ａに対してセル４の中心側に位置しているため、後述の分子動力学計算において、第１フィラーモデル１３Ａよりも分子鎖モデル６と優先的に相互作用するように計算することができる。

このように、フィラーモデル１３は、第１フィラーモデル１３Ａと、第１フィラーモデル１３Ａのフィラー粒子モデル１４を欠損させた第２フィラーモデル１３Ｂとが重ねられることにより、フィラーの表面に形成される凹凸を表現することができる。しかも、本実施形態の第１フィラーモデル１３Ａ及び第２フィラーモデル１３Ｂは、見かけ上の球面を構成するようにフィラー粒子モデル１４が配置されているため、現実のフィラーの表面形状（球に近い形状）に近似させることができる。

セル４内に定義されるフィラーモデル１３の個数としては、適宜設定することができる。本実施形態では、１個のフィラーモデル１３が定義されているが、複数個のフィラーモデル１３が定義されてもよい。なお、複数個のフィラーモデル１３が定義される場合、フィラーモデル１３、１３間の距離は、フィラーの体積分率に基づいて設定されるのが望ましい。

第１フィラーモデル１３Ａ及び第２フィラーモデル１３Ｂの半径方向において、第１フィラーモデル１３Ａと第２フィラーモデル１３Ｂとの間の距離Ｌ２（図９に示す）については、適宜設定することができる。本実施形態の距離Ｌ２は、高分子成分とフィラーとの相互作用を精度良く解析するために、例えば、グラファイトの壁間距離に基づいて設定することができる。なお、距離Ｌ２は、フィラー粒子モデル１４、１４の中心１４ｃ（図８に示す）間の距離として定義される。

第１フィラーモデル１３Ａのフィラー粒子モデル１４、及び、第２フィラーモデル１３Ｂのフィラー粒子モデル１４は、半径方向において互いに結合されてもよい。連結方法については、上述したフィラー粒子モデル１４、１４の連結方法（本実施形態では、図８に示した結合鎖モデル１５及びポテンシャルＰ２）のとおりである。これにより、フィラーモデル１３は、シリカのような三次元の網目構造に形成される。従って、フィラーモデル１３は、半径方向内側に空間が形成されていても、その剛性が大きく定義されるため、後述の分子動力学計算において、解析精度を高めることができる。フィラーモデル１３は、コンピュータ１に記憶される。

次に、本実施形態のシミュレーション方法では、ポリマー粒子モデル７、７間（即ち、結合鎖モデル８を介さずに隣り合うポリマー粒子モデル７、７間）、及び、フィラー粒子モデル１４とポリマー粒子モデル７との間に、引力及び斥力が作用するポテンシャルが定義される（工程Ｓ４）。図９は、ポテンシャルの一例を示す概念図である。

先ず、ポリマー粒子モデル７、７間のポテンシャルＰ３は、下記式（１）のＬＪポテンシャルＵ_LJ（ｒ_ij）によって定義される。

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

上記式（１）において、ポテンシャルＰ３は、ポリマー粒子モデル７、７間の距離ｒ_ijがカットオフ距離ｒ_c未満になる場合のみ、斥力及び引力を生じさせる。なお、ポリマー粒子モデル７、７間の距離ｒ_ijがカットオフ距離ｒ_c以上になった場合は、ポテンシャルＰ３がゼロとなり、斥力及び引力が生じない。

本実施形態では、炭素粒子モデル７Ｃ、７Ｃ間、水素粒子モデル７Ｈ、７Ｈ間、及び、炭素粒子モデル７Ｃと水素粒子モデル７Ｈとの間に、それぞれ異なるポテンシャルＰ３が設定されている。各ポテンシャルＰ３は、上記式（１）の定数がそれぞれ異なっている。なお、各定数は、例えば、上記論文１に基づいて、適宜設定することができる。ポテンシャルＰ３は、コンピュータ１に記憶される。

次に、フィラー粒子モデル１４とポリマー粒子モデル７との間のポテンシャルＰ４は、ポリマー粒子モデル７、７間のポテンシャルＰ３と同様に、上記式（１）のＬＪポテンシャルＵ_LJ（ｒ_ij）によって定義される。このようなポテンシャルＰ４により、フィラー粒子モデル１４とポリマー粒子モデル７との間に作用する引力及び斥力を定義することができる。ポテンシャルＰ４の各定数及び各変数の値としては、適宜設定することができるが、上記論文１に基づいて設定されるのが望ましい。

ポテンシャルＰ４は、第１フィラーモデル１３Ａの各フィラー粒子モデル１４、及び、第２フィラーモデル１３Ｂの各フィラー粒子モデル１４に、それぞれ定義されている。本実施形態において、第１フィラーモデル１３Ａの各フィラー粒子モデル１４に定義されるポテンシャルＰ４、及び、第２フィラーモデル１３Ｂの各フィラー粒子モデル１４に定義されるポテンシャルＰ４の各定数及び各変数は、同一に設定されている。上述したように、第２フィラーモデル１３Ｂは、第１フィラーモデル１３Ａよりもセル４の中心側（即ち、ポリマー粒子モデル７の近く）に位置している。このため、第２フィラーモデル１３Ｂの各フィラー粒子モデル１４に定義されるポテンシャルＰ４は、第１フィラーモデル１３Ａの各フィラー粒子モデル１４に定義されるポテンシャルＰ４に比べて、ポリマー粒子モデル７に対して優先的に作用させることができる。

さらに、フィラーモデル１３において、第２フィラーモデル１３Ｂのフィラー粒子モデル１４が欠損していない部分（以下、単に「非欠損部」ということがある。）１７では、フィラーモデル１３の半径方向において、第１フィラーモデル１３Ａのフィラー粒子モデル１４のポテンシャルＰ４と、第２フィラーモデル１３Ｂのフィラー粒子モデル１４のポテンシャルＰ４との双方を、ポリマー粒子モデル７に対して作用させることができる。他方、第２フィラーモデル１３Ｂのフィラー粒子モデル１４が欠損している部分（以下、単に「欠損部」ということがある。）１８では、半径方向において、第１フィラーモデル１３Ａのフィラー粒子モデル１４のポテンシャルＰ４のみを、ポリマー粒子モデル７に対して作用させることができる。

このように、本実施形態のシミュレーション方法では、非欠損部１７と、欠損部１８とで、分子鎖モデル６に対するポテンシャルＰ４の強度を変化させることができる。ポテンシャルＰ４は、コンピュータ１に記憶される。

次に、本実施形態のシミュレーション方法では、コンピュータ１が、図４に示したポリマー粒子モデル７、及び、フィラー粒子モデル１４を対象として分子動力学計算を行う（工程Ｓ５）。

本実施形態の分子動力学計算では、例えば、セル４について所定の時間、分子鎖モデル６が古典力学に従うものとして、ニュートンの運動方程式が適用される。そして、各時刻でのポリマー粒子モデル７の動きが、単位時間ステップ毎に追跡される。

分子動力学計算では、セル４において、圧力及び温度が一定、又は体積及び温度が一定に保たれる。これにより、工程Ｓ５では、実際の高分子材料の分子運動に近似させて、分子鎖モデル６の人為的な初期配置を精度よく緩和することができる。これにより、高分子材料モデル１０が設定される。構造緩和の計算は、例えば、従来と同様に、（株）ＪＳＯＬ社製のソフトマテリアル総合シミュレーター（Ｊ−ＯＣＴＡ）に含まれるＣＯＧＮＡＣを用いて処理することができる。

上述したように、本実施形態のシミュレーション方法では、図９に示したフィラーモデル１３の非欠損部１７と欠損部１８とにより、分子鎖モデル６とフィラーモデル１３との間のポテンシャルＰ４の強度を変化させることができる。さらに、本実施形態のフィラーモデル１３は、現実のフィラーの表面形状に近似させることができる。これにより、本実施形態の工程Ｓ５では、フィラーの表面状態（例えば、フィラーの原子に置き換わる官能基の存在）等に起因するフィラーと分子鎖との相互作用の変化を考慮して、分子鎖モデル６を緩和することができる。

また、フィラーモデル１３は、第２フィラーモデル１３Ｂの欠損部１８に分子鎖モデル６が浸入しても、欠損部のない第１フィラーモデル１３Ａによって、分子鎖モデル６の浸入を阻止できる。これにより、分子鎖モデル６がフィラーモデル１３の内部に浸入するような非現実的な計算を防ぐことができる。

次に、図３に示されるように、本実施形態のシミュレーション方法では、コンピュータ１が、分子鎖モデル６の初期配置を十分に緩和できたか否かを判断する（工程Ｓ６）。緩和の判断基準については、従来と同様に、適宜設定することができる。

工程Ｓ６において、分子鎖モデル６の初期配置を十分に緩和できたと判断された場合（工程Ｓ６において、「Ｙ」）、次の工程Ｓ７が実施される。他方、工程Ｓ６において、分子鎖モデル６の初期配置を十分に緩和できていないと判断された場合（工程Ｓ６において、「Ｎ」）、単位ステップを一つ進めて（工程Ｓ８）、工程Ｓ５及び工程Ｓ６が再度実施される。これにより、本実施形態のシミュレーション方法では、分子鎖モデル６の平衡状態（構造が緩和した状態）を確実に計算することができる。単位ステップは、例えば、０．５〜２．０fsに設定される。

次に、図３に示されるように、本実施形態のシミュレーション方法では、コンピュータ１が、図９に示したフィラーモデル１３と分子鎖モデル６との相互作用の強さを評価して、高分子成分とフィラーとの相互作用を解析する（工程Ｓ７）。

本実施形態の工程Ｓ７では、先ず、ポリマー粒子モデル７及びフィラー粒子モデル１４を対象に、予め定められた時間（例えば、２〜１０ns）分の分子動力学計算が、単位ステップ毎に計算される。各単位ステップにおいて、各分子鎖モデル６のポリマー粒子モデル７の座標が、コンピュータ１に記憶される。

上述したように、本実施形態のシミュレーション方法では、フィラーモデル１３の非欠損部１７と欠損部１８とにより、分子鎖モデル６とフィラーモデル１３との間のポテンシャルＰ４の強度を変化させることができる。これにより、本実施形態の工程Ｓ７での分子動力学計算では、フィラーの表面状態（例えば、フィラーの原子に置き換わる官能基の存在）等に起因するフィラーと分子鎖との相互作用の変化を考慮して、分子鎖モデル６とフィラーモデル１３との相互作用を、精度良く計算することができる。さらに、本実施形態のフィラーモデル１３は、現実のフィラーの表面形状に近似させることができるため、現実のフィラー近傍での分子鎖の分布を精度良く再現することができる。

また、フィラーモデル１３は、第２フィラーモデル１３Ｂの欠損部１８に分子鎖モデル６が浸入しても、欠損部のない第１フィラーモデル１３Ａによって、分子鎖モデル６の浸入を防ぐことができる。これにより、分子鎖モデル６がフィラーモデル１３の内部に浸入するような非現実的な計算を防ぐことができる。

次に、本実施形態の工程Ｓ７では、コンピュータ１に記憶されているポリマー粒子モデル７の座標値に基づいて、フィラーモデル１３と分子鎖モデル６との相互作用の強さが評価される。相互作用の強さの評価については、適宜行うことができる。本実施形態では、分子鎖モデル６のポリマー粒子モデル７の平均二乗変位（ Mean-square Displacement：ＭＳＤ）に基づいて、フィラーモデル１３と分子鎖モデル６との相互作用の強さが評価される。本実施形態では、ポリマー粒子モデル７のうち、炭素粒子モデル７Ｃのみを対象に、平均二乗変位が計算される。平均二乗変位は、下記式（２）で定義される。

ＭＳＤ＝＜｜ｒ（τ）−ｒ（０）｜^２＞ …（２）
なお、符号は次のとおりである。
ｒ（０）：任意の時刻におけるポリマー粒子モデルの座標値
ｒ（τ）：任意の時刻から時間τ後におけるポリマー粒子モデルの座標値
＜＞：アンサンブル平均

上記式（２）において、上記｜ｒ（τ）−ｒ（０）｜は、ポリマー粒子モデル７が時間幅τで移動した距離である。この距離｜ｒ（τ）−ｒ（０）｜を二乗した値について、評価対象となるポリマー粒子モデル７に対してアンサンブル平均をとることにより、時間幅τにおける平均二乗変位が求められる。なお、各座標値ｒ（τ）、ｒ（０）は、各ポリマー粒子モデル７の中心７ｃで特定されるものとする。なお、「τ」は、シミュレーションに用いられる時間のパラメータである。

平均二乗変位は、各ポリマー粒子モデル７が時間幅τで移動した距離の二乗を、評価対象のポリマー粒子モデル７でアンサンブル平均した値である。このような平均二乗変位の値を調べることにより、ポリマー粒子モデル７の運動量を把握することができる。

平均二乗変位が小さい場合、ポリマー粒子モデル７の運動量が小さいことを示している。このため、フィラーモデル１３と分子鎖モデル６との相互作用が強いと評価することができる。これにより、工程Ｓ７では、高分子成分とフィラーとの相互作用が大きいと推定することができる。

他方、平均二乗変位が大きい場合、ポリマー粒子モデル７の運動量が大きいことを示している。このため、フィラーモデル１３と分子鎖モデル６との相互作用が弱いと評価することができる。これにより、工程Ｓ７では、高分子成分とフィラーとの相互作用が小さいと推定することができる。

このように、工程Ｓ７では、平均二乗変位に基づいて、フィラーモデル１３と分子鎖モデル６との相互作用の強さ、及び、高分子成分とフィラーとの相互作用を、容易かつ正確に評価することができる。

上述したように、フィラーモデル１３と分子鎖モデル６との相互作用の強さの評価は、フィラーの表面状態（例えば、フィラーの原子に置き換わる官能基の存在）等に起因するフィラーと分子鎖との相互作用の変化を考慮した分子動力学計算の結果に基づいて行われている。このため、工程Ｓ７では、フィラーモデル１３と分子鎖モデル６との相互作用の強さに基づいて、高分子成分とフィラーとの相互作用を、精度良く解析することができる。高分子成分とフィラーとの相互作用の解析結果は、高分子材料の開発に役立たせることができる。

これまでの実施形態では、図９に示したフィラーモデル１３の非欠損部１７と欠損部１８とにより、フィラーと分子鎖との相互作用の変化が考慮されたが、このような態様に限定されない。例えば、フィラーモデル１３の非欠損部１７と欠損部１８だけでなく、さらに官能基をモデル化した官能基モデルを用いて、フィラーと分子鎖との相互作用の変化が考慮されてもよい。

図１０は、本発明の他の実施形態のシミュレーション方法の処理手順を示すフローチャートである。図１１は、本発明の他の実施形態の第２フィラーモデル１３Ｂ及び官能基モデル２１を示す部分斜視図である。なお、図１１では、第１フィラーモデル１３Ａを省略して表示している。この実施形態において、前実施形態と同一の構成については、同一の符号を付し、説明を省略することがある。

この実施形態のシミュレーション方法では、分子動力学計算を行う工程Ｓ５に先立って、官能基モデル２１を、フィラーモデル１３に結合させる（工程Ｓ９）。官能基モデル２１は、官能基を少なくとも一つの官能基粒子モデル２２でモデル化したものである。官能基モデル２１としては、例えば、アルデヒド（−ＣＨＯ）をモデル化したアルデヒドモデル２１Ａと、ケトン（−ＣＯ−）をモデル化したケトンモデル２１Ｂとを含んでいる。

上述したように、官能基は、フィラーの表面に、フィラーの原子に置き換わって存在しており、フィラーと分子鎖との相互作用に変化を生じさせている。このため、本実施形態の官能基モデル２１は、第２フィラーモデル１３Ｂの欠損部１８に配置される。欠損部１８に配置された官能基モデル２１は、隣接する第２フィラーモデル１３Ｂのフィラー粒子モデル１４に連結されている。

官能基粒子モデル２２は、分子動力学計算において、運動方程式の質点として取り扱われる。即ち、ポリマー粒子モデル７は、質量、直径、電荷、又は、初期座標などのパラメータが定義される。本実施形態の官能基粒子モデル２２は、図１１に示されるように、炭素粒子モデル２２ｃ、水素粒子モデル２２ｈ、及び、酸素粒子モデル２２ｏを含んでいる。

官能基粒子モデル２２、２２間の連結、及び、官能基粒子モデル２２とフィラー粒子モデル１４との連結は、適宜定義することができる。本実施形態の官能基粒子モデル２２、２２間の連結、及び、官能基粒子モデル２２とフィラー粒子モデル１４との連結は、結合鎖モデル２３によって連結されている。結合鎖モデル２３を介して隣り合う官能基粒子モデル２２、２２間、官能基粒子モデル２２とフィラー粒子モデル１４との間には、引力及び斥力が作用するポテンシャルＰ５が定義される。ポテンシャルＰ５については、従来と同様に、適宜定義することができる。本実施形態のポテンシャルＰ５は、上記論文１に基づいて、官能基の構造に応じて設定されている。官能基モデル２１は、コンピュータ１に記憶される。

次に、この実施形態のシミュレーション方法では、官能基粒子モデル２２とポリマー粒子モデル７との間、及び、官能基粒子モデル２２とフィラー粒子モデル１４との間に、引力及び斥力が作用するポテンシャルＰ６が定義される（工程Ｓ１０）。ポテンシャルＰ６は、上記式（１）のＬＪポテンシャルＵ_LJ（ｒ_ij）によって定義される。ポテンシャルＰ６の各定数及び各変数の値としては、適宜設定することができるが、上記論文１に基づいて設定されるのが望ましい。

この実施形態の分子動力学計算を行う工程Ｓ５では、コンピュータ１が、ポリマー粒子モデル７、フィラー粒子モデル１４、及び、官能基粒子モデル２２を対象として、分子動力学計算を行う（工程Ｓ５）。分子動力学計算の詳細については、前実施形態で説明したとおりである。

この実施形態において、フィラーモデル１３と分子鎖モデル６との相互作用の強さを評価する工程Ｓ７では、先ず、前実施形態と同様に、ポリマー粒子モデル７、フィラー粒子モデル１４、及び、官能基粒子モデル２２を対象に、分子動力学計算が行われる。そして、ポリマー粒子モデル７、官能基粒子モデル２２、及び、フィラー粒子モデル１４の距離に基づいて、フィラーモデル１３と分子鎖モデル６との相互作用の強さが評価される。相互作用の強さの評価方法については、前実施形態で説明したとおりである。

この実施形態では、フィラーモデル１３の非欠損部１７と欠損部１８とによって変化するポテンシャルＰ４（図９に示す）だけでなく、官能基粒子モデル２２とポリマー粒子モデル７との間のポテンシャルＰ６（図１１に示す）も、分子鎖モデル６に対して作用している。これにより、本実施形態のシミュレーション方法では、フィラーの原子に置き換わる官能基の存在に起因するフィラーと分子鎖との相互作用の変化を考慮することができるため、分子鎖モデル６とフィラーモデル１３との相互作用（高分子成分とフィラーとの相互作用）を、より精度良く計算することができる。また、官能基粒子モデル２２は、分子動力学計算において、実際の官能基のように運動できるため、フィラーモデル１３と分子鎖モデル６との相互作用の強さを評価する工程Ｓ７において、高分子成分とフィラーとの相互作用をより精度良く評価することができる。

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

図３に示した処理手順に従って、仮想空間であるセル内に、分子鎖モデル、及び、第１フィラーモデルと第２フィラーモデルとを含む２層のフィラーモデルが定義された（実施例１）。第１フィラーモデル及び第２フィラーモデルのフィラー粒子モデルは、見かけ上の球面を構成するように配置された。さらに、図１０に示した処理手順に従い、実施例１の分子鎖モデル及び２層のフィラーモデルに加えて、官能基モデルがさらに定義された（実施例２及び実施例３）。

比較のために、仮想空間であるセル内に、分子鎖モデル、及び、第１フィラーモデルと第２フィラーモデルとを含む２層のフィラーモデルが定義された（比較例１）。比較例１の第１フィラーモデル及び第２フィラーモデルのフィラー粒子モデルは、見かけ上の平面を構成するように配置された。さらに、比較例１の分子鎖モデル及び２層のフィラーモデルに加えて、官能基モデルがさらに定義された（比較例２及び比較例３）。

さらに、フィラー粒子モデルの欠損を有しない第２フィラーモデルを含む２層のフィラーモデルが定義された（比較例４）。また、仮想空間であるセル内に、分子鎖モデル、及び、第２フィラーモデルのみを含む１層のフィラーモデルが定義された（比較例５、比較例６）。

そして、実施例１〜３、及び、比較例１〜６の分子鎖モデル、フィラーモデル、及び、官能基モデルを対象に、分子動力学計算（５ns）がコンピュータを用いて実施され、フィラーと分子鎖との相互作用の強さが評価された。なお、相互作用の強さは、実施例１の平均二乗変位の値を１００とする指数で表示されている。数値が大きいほど、相互作用が強いことを示している。共通仕様は、次のとおりである。

セル：一辺の各長さＬ１：３．６８〜３．８３nm
高分子成分：ポリイソプレン（密度：０．９〜１．０ｇ／cm³）
分子鎖モデル：
個数：セル内に配置される全ての分子鎖モデルの密度が、
ポリイソプレンの上記の密度に一致するように調整
フィラー：グラファイト（カーボンブラック）
フィラーモデルの個数：１個
第２フィラーモデル：
フィラー粒子モデルの欠損の割合：２１％
官能基：アルデヒド、ケトン
テストの結果を表に示す。

一般に、グラファイトの表面が活性化する（即ち、官能基モデルの個数、又は、欠損の割合が大きくなる）ほど、相互作用（界面接着力）が強くなることが知られている（例えば、「カーボンブラック便覧（第３版）」、カーボンブラック協会、１９９５年、Ｐ２１６）。実施例１〜３において、官能基モデルの個数が大きくなるほど、フィラーと分子鎖との相互作用が強くなる傾向を示した。従って、実施例１〜３の結果より、本発明のシミュレーション方法は、上記文献「カーボンブラック便覧（第３版）」に記載のように、高分子成分とフィラーとの相互作用を精度良く解析できた。

他方、比較例１〜３は、官能基モデルの個数が大きくなるほど、フィラーと分子鎖との相互作用を強くすることができたが、上記文献「カーボンブラック便覧（第３版）」に記載の傾向と一致しなかった。このため、比較例１〜３は、実施例に比べて、高分子成分とフィラーとの相互作用を、精度良く解析することができなかった。

比較例４は、第２フィラーモデルに欠損部及び官能基モデルが定義されなかったため、高分子成分とフィラーとの相互作用を、精度良く解析することができなかった。第２フィラーモデルのみを含む１層のフィラーモデルが定義された比較例５及び比較例６は、欠損部分に分子鎖モデルが浸入し、分子鎖モデルが欠損部分からセルの面にはみ出した。このため、比較例５及び比較例６は、高分子成分とフィラーとの相互作用を解析することができなかった。

４ セル
６ 分子鎖モデル
１３ フィラーモデル
１３Ａ 第１フィラーモデル
１３Ｂ 第２フィラーモデル
１４ フィラー粒子モデル

Claims (4)

1. 高分子材料に含まれる高分子成分とフィラーとの相互作用を、コンピュータを用いて解析するための方法であって、
   前記高分子材料の一部に対応する仮想空間であるセルを、前記コンピュータに入力する工程と、
   前記セル内に、複数のポリマー粒子モデルを用いて前記高分子成分の分子鎖をモデル化した分子鎖モデルを定義する工程と、
   前記セル内に、複数のフィラー粒子モデルを用いて前記フィラーをモデル化したフィラーモデルを定義する工程と、
   前記ポリマー粒子モデル間、及び、前記フィラー粒子モデルと前記ポリマー粒子モデルとの間に、引力及び斥力が作用するポテンシャルを定義する工程と、
   前記コンピュータが、前記ポリマー粒子モデル、及び、前記フィラー粒子モデルを対象として分子動力学計算を行う工程とを含み、
   前記フィラーモデルは、第１フィラーモデルと第２フィラーモデルとを含んで定義され、
   前記第１フィラーモデルは、複数の前記フィラー粒子モデルが、見かけ上の球面を構成するように配置されかつ互いに連結されており、
   前記第２フィラーモデルは、前記第１フィラーモデルの前記フィラー粒子モデルの少なくとも一つが欠損したものであり、
   前記第２フィラーモデルは、前記セル内に定義された前記分子鎖モデルと相互作用するように、前記セルの中心側に位置して前記第１フィラーモデルと重ねられている高分子材料のシミュレーション方法。
2. 前記第１フィラーモデルの前記フィラー粒子モデル、及び、前記第２フィラーモデルの前記フィラー粒子モデルは、それらの半径方向で互いに結合されている請求項１記載の高分子材料のシミュレーション方法。
3. 前記高分子材料は、前記フィラーに対する前記分子鎖の親和性を高める官能基を含み、
   前記官能基を少なくとも一つの官能基粒子モデルでモデル化した官能基モデルを、前記フィラーモデルの前記フィラー粒子モデルに結合させる工程と、
   前記官能基粒子モデルと前記ポリマー粒子モデルとの間に、前記ポテンシャルを定義する工程とをさらに含み、
   前記分子動力学計算を行う工程は、前記ポリマー粒子モデル、前記フィラー粒子モデル、及び、前記官能基粒子モデルを対象として、前記分子動力学計算を行う請求項１又は２に記載の高分子材料のシミュレーション方法。
4. 前記官能基粒子モデルは、前記第２フィラーモデルの前記欠損した部分に配置される請求項３記載の高分子材料のシミュレーション方法。