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

Description

本発明は、フィラーの分散性を評価及び改善に役立つ高分子材料のシミュレーション方法に関する。

ゴム等の高分子材料には、カーボンブラックやシリカといったフィラーが配合されている。これらのフィラーの分散性は、ゴムの強度に大きな影響を与えることが明らかになっている。近年、ゴム配合の開発のために、高分子材料中に配合されたフィラーの分散性を、コンピュータを用いて評価するためのシミュレーション方法（数値計算）が種々提案されている。

この種のシミュレーション方法では、先ず、フィラーを複数のフィラー粒子モデルでモデル化したフィラーモデル、及び、高分子材料を複数のポリマー粒子モデルでモデル化したポリマーモデルがコンピュータに設定される。そして、コンピュータが、フィラー粒子モデルにポリマー粒子モデルを結合して、分子動力学計算が実施される。

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

上記のようなシミュレーション方法では、各フィラー粒子モデルに結合させるポリマー粒子モデルが、コンピュータによってランダムに決定されていた。このため、各フィラー粒子モデルに結合されるポリマー粒子モデルの数に大きな偏りが生じていた。各フィラー粒子モデルに結合されるポリマー粒子モデルの数が大きくなると、各フィラー粒子モデルに作用するポリマー粒子モデルの力の総和が大きくなり、分子動力学計算に計算落ちが発生するという問題があった。

本発明は、以上のような実状に鑑み案出されたもので、各結合対象フィラー粒子モデルに結合させる結合対象ポリマー粒子モデルの数を、予め定められた上限値以下に限定することを基本として、分子動力学計算に計算落ちが発生するのを防ぐことができる高分子材料のシミュレーション方法を提供することを主たる目的としている。

本発明に係る高分子材料のシミュレーション方法は、高分子材料中に配合されたフィラーの分散性を、コンピュータを用いて評価するためのシミュレーション方法であって、前記コンピュータに、前記フィラーを、複数のフィラー粒子モデルでモデル化したフィラーモデルを設定する工程、前記コンピュータに、前記高分子材料を、複数のポリマー粒子モデルでモデル化したポリマーモデルを設定する工程、及び前記コンピュータが、予め定められた仮想空間に配置された前記フィラーモデルと前記ポリマーモデルとを用いて分子動力学計算を行うシミュレーション工程を含み、前記シミュレーション工程は、前記フィラー粒子モデルから予め定められた結合対象フィラー粒子モデルに、前記ポリマー粒子モデルから予め定められた結合対象ポリマー粒子モデルを結合させる結合工程を含み、前記各結合対象フィラー粒子モデルに結合させる前記結合対象ポリマー粒子モデルの数は、予め定められた上限値Ｓｕ以下であることを特徴とする。

本発明に係る前記高分子材料のシミュレーション方法において、前記上限値Ｓｕは、下記の関係を満たすのが望ましい。
Ｎａ／Ｎｂ−１≦Ｓｕ≦Ｎａ／Ｎｂ＋１
ここで、各パラメータは、次のとおりである。
Ｎａ：結合対象ポリマー粒子モデルの総数
Ｎｂ：結合対象フィラー粒子モデルの総数

本発明に係る前記高分子材料のシミュレーション方法において、前記結合工程は、前記上限値Ｓｕを設定する工程、前記各結合対象フィラー粒子モデルに結合させる前記結合対象ポリマー粒子モデルを特定する特定工程、及び特定された前記結合対象ポリマー粒子モデルを、前記結合対象フィラー粒子モデルに結合させる工程を含むのが望ましい。

本発明に係る前記高分子材料のシミュレーション方法において、前記特定工程は、前記各結合対象ポリマー粒子モデルを、夫々粒子間距離が最も小さい前記結合対象フィラー粒子モデルに関連付ける第１工程、前記各結合対象フィラー粒子モデルに結合させる前記結合対象ポリマー粒子モデルを、前記上限値Ｓｕ以下を満たすように、関連付けられた結合対象ポリマー粒子モデル群から粒子間距離が小さい順に特定する第２工程、前記第２工程で特定されなかった前記各結合対象ポリマー粒子モデルを、前記関連付けられた結合対象フィラー粒子モデルの次に粒子間距離が小さい前記結合対象フィラー粒子モデルに関連付ける第３工程、及び前記各結合対象フィラー粒子モデルに結合させる前記結合対象ポリマー粒子モデルを、前記上限値Ｓｕ以下を満たすように、前記第３工程で関連付けられた結合対象ポリマー粒子モデル群から粒子間距離が小さい順に特定する第４工程を含み、全ての前記結合対象ポリマー粒子モデルが前記結合対象フィラー粒子モデルに特定されるまで、前記第３工程及び前記第４工程が繰り返し実施されるのが望ましい。

本発明は、フィラーと、前記フィラーをポリマーに結合させるためのカップリング剤とが配合された高分子材料について、前記フィラーの分散性を、コンピュータを用いて評価するためのシミュレーション方法であって、前記コンピュータに、前記フィラーを、複数のフィラー粒子モデルでモデル化したフィラーモデルを設定する工程、前記コンピュータに、前記高分子材料を、複数のポリマー粒子モデルでモデル化したポリマーモデルを設定する工程、前記コンピュータに、前記カップリング剤を、少なくとも一つのカップリング粒子モデルでモデル化したカップリングモデルを設定する工程、並びに前記コンピュータが、予め定められた仮想空間に配置された前記フィラーモデル、前記ポリマーモデル、及び、前記カップリングモデルを用いて分子動力学計算を行うシミュレーション工程を含み、前記シミュレーション工程は、前記フィラー粒子モデルから予め定められた結合対象フィラー粒子モデルに、前記ポリマー粒子モデルから予め定められた結合対象ポリマー粒子モデル、又は、前記カップリング粒子モデルから予め定められた結合対象カップリング粒子モデルを結合させる結合工程を含み、前記各結合対象フィラー粒子モデルに結合させる前記結合対象ポリマー粒子モデルの数と、前記各結合対象フィラー粒子モデルに結合させる前記結合対象カップリング粒子モデルの数との和は、予め定められた上限値Ｎｓ以下であるのが望ましい。

本発明に係る前記高分子材料のシミュレーション方法において、前記上限値Ｎｓは、下記の関係を満たすのが望ましい。
Ｎｓ＝［（Ｎａ＋Ｎｃ）／Ｎｂ＋１］
ここで、各パラメータは、次のとおりである。
Ｎａ：結合対象ポリマー粒子モデルの総数
Ｎｂ：結合対象フィラー粒子モデルの総数
Ｎｃ：結合対象カップリング粒子モデルの総数

本発明に係る前記高分子材料のシミュレーション方法において、前記結合工程は、前記上限値Ｎｓを設定する工程、前記各結合対象フィラー粒子モデルに結合させる前記結合対象ポリマー粒子モデル、又は、前記各結合対象フィラー粒子モデルに結合させる前記結合対象カップリング粒子モデルを特定する特定工程、並びに特定された前記結合対象ポリマー粒子モデル又は前記結合対象カップリング粒子モデルを、前記結合対象フィラー粒子モデルに結合させる工程を含むのが望ましい。

本発明に係る前記高分子材料のシミュレーション方法において、前記特定工程は、前記各結合対象ポリマー粒子モデル及び前記各結合対象カップリング粒子モデルを、夫々粒子間距離が最も小さい前記結合対象フィラー粒子モデルに関連付ける第１工程、前記各結合対象フィラー粒子モデルに結合させる前記結合対象ポリマー粒子モデル又は前記各結合対象カップリング粒子モデルを、前記上限値Ｎｓ以下を満たすように、関連付けられた結合対象ポリマー粒子モデル又は結合対象カップリング粒子モデルから粒子間距離が小さい順に特定する第２工程、前記第２工程で特定されなかった前記各結合対象ポリマー粒子モデル又は前記各結合対象カップリング粒子モデルを、前記関連付けられた結合対象フィラー粒子モデルの次に粒子間距離が小さい前記結合対象フィラー粒子モデルに関連付ける第３工程、並びに前記各結合対象フィラー粒子モデルに結合させる前記結合対象ポリマー粒子モデル又は前記結合対象カップリング粒子モデルを、前記上限値Ｎｓ以下を満たすように、前記第３工程で関連付けられた結合対象ポリマー粒子モデル又は結合対象カップリング粒子モデルから粒子間距離が小さい順に特定する第４工程を含み、全ての前記結合対象ポリマー粒子モデル又は前記結合対象カップリング粒子モデルが前記結合対象フィラー粒子モデルに特定されるまで、前記第３工程及び前記第４工程が繰り返し実施されるのが望ましい。

本発明に係る前記高分子材料のシミュレーション方法において、前記フィラー粒子モデルは、予め定められた多面体の頂点に少なくとも配置され、前記結合対象フィラー粒子モデルは、前記頂点に配置されたフィラー粒子モデルであるのが望ましい。

請求項１記載の高分子材料のシミュレーション方法は、コンピュータが、予め定められた仮想空間に配置されたフィラーモデルとポリマーモデルとを用いて分子動力学計算を行うシミュレーション工程を含んでいる。シミュレーション工程は、フィラー粒子モデルから予め定められた結合対象フィラー粒子モデルに、ポリマー粒子モデルから予め定められた結合対象ポリマー粒子モデルを結合させる結合工程を含んでいる。各結合対象フィラー粒子モデルに結合させる結合対象ポリマー粒子モデルの数は、予め定められた上限値Ｓｕ以下である。

請求項１記載の高分子材料のシミュレーション方法では、各結合対象フィラー粒子モデルに結合される結合対象ポリマー粒子モデルの数に大きな偏りが生じるのを抑制することができる。これにより、請求項１記載の高分子材料のシミュレーション方法では、各結合対象フィラー粒子モデルに作用する結合対象ポリマー粒子モデルの力の総和が大きくなるのを防ぐことができるため、計算落ちが発生するのを防ぐことができる。従って、請求項１記載の高分子材料のシミュレーション方法では、分子動力学計算を安定して実行することができる。

請求項５記載の高分子材料のシミュレーション方法は、コンピュータが、予め定められた仮想空間に配置されたフィラーモデル、ポリマーモデル、及び、カップリングモデルを用いて分子動力学計算を行うシミュレーション工程を含んでいる。シミュレーション工程は、フィラー粒子モデルから予め定められた結合対象フィラー粒子モデルに、ポリマー粒子モデルから予め定められた結合対象ポリマー粒子モデル、又は、カップリング粒子モデルから予め定められた結合対象カップリング粒子モデルを結合させる結合工程を含んでいる。各結合対象フィラー粒子モデルに結合させる結合対象ポリマー粒子モデルの数と、各結合対象フィラー粒子モデルに結合させる結合対象カップリング粒子モデルの数との和は、予め定められた上限値Ｎｓ以下である。

請求項５記載の高分子材料のシミュレーション方法では、各結合対象フィラー粒子モデルに結合される結合対象ポリマー粒子モデルの数と、各結合対象フィラー粒子モデルに結合させる結合対象カップリング粒子モデルの数との和に、大きな偏りが生じるのを抑制することができる。これにより、請求項５記載の高分子材料のシミュレーション方法では、各結合対象フィラー粒子モデルに作用する結合対象ポリマー粒子モデルの力、又は、結合対象カップリング粒子モデルの力の総和が大きくなるのを防ぐことができるため、計算落ちが発生するのを防ぐことができる。従って、請求項５記載の高分子材料のシミュレーション方法では、分子動力学計算を安定して実行することができる。

本実施形態の高分子材料のシミュレーション方法を実行するコンピュータの斜視図である。 本実施形態の高分子材料のシミュレーション方法の処理手順の一例を示すフローチャートである。 本実施形態のフィラーモデルの概念図である。 本実施形態のポリマーモデルの概念図である。 本実施形態の条件設定工程の処理手順の一例を示すフローチャートである。 フィラー粒子モデルとポリマー粒子モデルとのポテンシャルを説明する概念図である。 フィラー粒子モデルのカットオフ距離を説明する概念図である。 シミュレーションモデルの仮想空間の概念斜視図である。 本実施形態のシミュレーション工程の処理手順の一例を示すフローチャートである。 本実施形態の結合工程の処理手順の一例を示すフローチャートである。 本実施形態の特定工程の処理手順の一例を示すフローチャートである。 結合対象ポリマー粒子モデルを関連付ける工程を説明する概念図である。 結合対象ポリマー粒子モデルを特定する工程を説明する概念図である。 特定されなかった結合対象ポリマー粒子モデルの関連付けを説明する概念図である。 第３工程で関連付けられた結合対象ポリマー粒子モデルを特定する工程を説明する概念図である。 結合対象フィラー粒子モデルと結合対象ポリマー粒子モデルとを結合する工程を説明する概念図である。 分子動力学計算後のシミュレーションモデルの概念斜視図である。 本発明の他の実施形態のシミュレーション方法の処理手順の一例を示すフローチャートである。 カップリングモデルの概念図である。 本発明の他の実施形態の条件設定工程の処理手順の一例を示すフローチャートである。 カップリング粒子モデルのポテンシャルを説明する概念図である。 本発明の他の実施形態のシミュレーションモデルの仮想空間９の概念斜視図である。 本発明の他の実施形態のシミュレーション工程の処理手順の一例を示すフローチャートである。 （ａ）は、ポリマー粒子モデルとカップリング粒子モデルとが結合される前の状態を示す概念図、（ｂ）は、リマー粒子モデルとカップリング粒子モデルとが結合された状態を示す概念図である。 本発明の他の実施形態の結合工程の処理手順の一例を示すフローチャートである。 本発明の他の実施形態の特定工程の処理手順の一例を示すフローチャートである。 結合対象ポリマー粒子モデル又は結合対象カップリング粒子モデルを関連付ける工程を説明する概念図である。 結合対象ポリマー粒子モデル又は結合対象カップリング粒子モデルを特定する工程を説明する概念図である。 特定されなかった結合対象ポリマー粒子モデル又は結合対象カップリング粒子モデルの関連付けを説明する概念図である。 第３工程で関連付けられた結合対象ポリマー粒子モデル及び結合対象カップリング粒子モデルを特定する工程を説明する概念図である。 特定された結合対象ポリマー粒子モデル又は結合対象カップリング粒子モデルを、結合対象フィラー粒子モデルに結合する工程を説明する概念図である。 実施例２及び実施例３の動径分布関数を示すグラフである。 実施例２及び実施例４の動径分布関数を示すグラフである。

以下、本発明の実施の一形態が図面に基づき説明される。
本実施形態の高分子材料のシミュレーション方法（以下、単に「シミュレーション方法」ということがある）は、高分子材料中に配合されたフィラーの分散性を、コンピュータを用いて評価するための方法である。フィラーとしては、例えば、カーボンブラック、シリカ又はアルミナ等が含まれる。高分子材料としては、例えば、ゴム、樹脂又はエラストマー等が含まれる。

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

図２は、本実施形態のシミュレーション方法の処理手順の一例を示すフローチャートである。本実施形態のシミュレーション方法では、先ず、フィラーをモデル化したフィラーモデルが設定される（工程Ｓ１）。図３は、本実施形態のフィラーモデルの概念図である。

本実施形態のフィラーモデル３は、複数のフィラー粒子モデル４を用いて、フィラーがモデル化される。フィラー粒子モデル４は、予め定められた多面体１３に基づいて配置されている。本実施形態の多面体１３は、立方体として設定されている。

フィラー粒子モデル４は、多面体１３の重心に配置される１個の中心粒子モデル４ｃと、多面体１３の表面に配置される複数の表面粒子モデル４ｓとを含んでいる。本実施形態の表面粒子モデル４ｓは、多面体１３の頂点１３ａに配置される頂点粒子モデル７と、頂点粒子モデル７、７間に配置される中間粒子モデル８とを含んで構成されている。

中心粒子モデル４ｃ及び表面粒子モデル４ｓは、いずれも径を持った球で表現されている。さらに、中心粒子モデル４ｃと表面粒子モデル４ｓ（頂点粒子モデル７）との間、及び、各表面粒子モデル４ｓ、４ｓ（頂点粒子モデル７及び中間粒子モデル８）の間には、それぞれ平衡長が定義された結合鎖４ｊが設けられている。このようなフィラーモデル３は、コンピュータ１で取り扱い可能な数値データであり、コンピュータ１に入力される。数値データには、例えば、フィラー粒子モデル４の質量、体積、直径及び初期座標などが含まれる。

ここで、「平衡長」とは、多面体１３での各表面粒子モデル４ｓの相対位置を安定して保ちうる中心粒子モデル４ｃと表面粒子モデル４ｓとの間、及び、各表面粒子モデル４ｓ、４ｓの間の結合距離である。これにより、フィラーモデル３は、中心粒子モデル４ｃ及び表面粒子モデル４ｓの配置を、多面体１３の形状に安定して維持することができる。

次に、図２に示されるように、本実施形態のシミュレーション方法では、高分子材料をモデル化したポリマーモデルが設定される（工程Ｓ２）。図４は、本実施形態のポリマーモデル５の概念図である。

ポリマーモデル５は、複数のポリマー粒子モデル６を用いて、高分子材料がモデル化される。本実施形態のポリマー粒子モデル６は、それぞれ異なるポテンシャルが定義される未変性粒子モデル６ａ、及び、変性基粒子モデル６ｂを含んでいる。これらのポリマー粒子モデル６ａ、６ｂは、いずれも径を持った球で表現されている。

未変性粒子モデル６ａ、６ａ間、及び、未変性粒子モデル６ａと変性基粒子モデル６ｂとの間には、平衡長が定義された結合鎖６ｊが定義されている。これにより、ポリマーモデル５は、直鎖状の三次元構造に形成される。このようなポリマーモデル５は、コンピュータ１で取り扱い可能な数値データであり、コンピュータ１に入力される。数値データには、例えば、ポリマー粒子モデル６の質量、体積、直径及び初期座標などが含まれる。

次に、本実施形態のシミュレーション方法では、シミュレーション条件が設定される（条件設定工程Ｓ３）。図５は、本実施形態の条件設定工程Ｓ３の処理手順の一例を示すフローチャートである。

この条件設定工程Ｓ３では、先ず、各フィラー粒子モデル４ｃ、４ｓと各ポリマー粒子モデル６ａ、６ｂとの間に、粒子間の距離の関数であるポテンシャルが定義される（工程Ｓ３１）。図６は、フィラー粒子モデル４とポリマー粒子モデル６とのポテンシャルを説明する概念図である。ポテンシャルは、２つの粒子モデルの間に作用する力を計算する際に用いられるものである。このポテンシャルは、下記の式（１）で定義される。

上記式（１）において、ａ_ijは、粒子モデル間ごとに定義されるポテンシャルの強度を示している。ｒ_ijは、粒子間の距離を示している。ｒ_cは、カットオフ距離を示している。なお、距離ｒ_ij、及び、カットオフ距離ｒ_cは、各粒子モデルの中心間の距離として定義される。

上記式（１）では、２つの粒子間の距離が、予め定められたカットオフ距離ｒ_c以下になったときに、相互作用（本実施形態では、斥力）が生じるポテンシャルが定義される。なお、粒子間の距離ｒ_ijが予め定めたカットオフ距離ｒ_cよりも大きい場合には、ポテンシャルＵがゼロとなり、粒子モデル間に斥力は生じない。

本実施形態では、次の２つの粒子モデルの組合せについて、それぞれポテンシャルＵ１乃至ポテンシャルＵ１０が定義される。
ポテンシャルＵ１：中心粒子モデル４ｃ−未変性粒子モデル６ａ
ポテンシャルＵ２：中心粒子モデル４ｃ−変性基粒子モデル６ｂ
ポテンシャルＵ３：中心粒子モデル４ｃ−表面粒子モデル４ｓ
ポテンシャルＵ４：表面粒子モデル４ｓ−未変性粒子モデル６ａ
ポテンシャルＵ５：表面粒子モデル４ｓ−変性基粒子モデル６ｂ
ポテンシャルＵ６：未変性粒子モデル６ａ−変性基粒子モデル６ｂ
ポテンシャルＵ７：中心粒子モデル４ｃ−中心粒子モデル４ｃ
ポテンシャルＵ８：表面粒子モデル４ｓ−表面粒子モデル４ｓ
ポテンシャルＵ９：未変性粒子モデル６ａ−未変性粒子モデル６ａ
ポテンシャルＵ１０：変性基粒子モデル６ｂ−変性基粒子モデル６ｂ

ポテンシャルの強度ａ_ijとしては、論文（J. Chem Phys. 107(11) 4423-4435 (1997) ）において、同種粒子間の強度ａ_ijを２５とすることが提唱されている。しかし、その後、多くの研究がなされ、同種粒子間ではａ_ij＝５０、異種粒子間ではａ_ij＝７２とするものが出てきた（例えば、論文（Macromolcule vol.39 6744(2006)））。本実施形態では、上記パラメータを参考として、ポテンシャルＵ１乃至ポテンシャルＵ１０の各強度ａ_ijは、次のように設定される。
ポテンシャルＵ１：ａ_ij＝７２
ポテンシャルＵ２：ａ_ij＝２５
ポテンシャルＵ３：ａ_ij＝５０
ポテンシャルＵ４：ａ_ij＝７２
ポテンシャルＵ５：ａ_ij＝２５
ポテンシャルＵ６：ａ_ij＝７２
ポテンシャルＵ７：ａ_ij＝５０
ポテンシャルＵ８：ａ_ij＝５０
ポテンシャルＵ９：ａ_ij＝５０
ポテンシャルＵ１０：ａ_ij＝５０

上記のように、変性基粒子モデル６ｂと中心粒子モデル４ｃとの間のポテンシャルＵ２、及び、変性基粒子モデル６ｂと表面粒子モデル４ｓとの間のポテンシャルＵ５の強度ａ_ij（＝２５）は、未変性粒子モデル６ａと中心粒子モデル４ｃとの間のポテンシャルＵ１、及び、未変性粒子モデル６ａと表面粒子モデル４ｓとの間のポテンシャルＵ４の強度ａ_ij（＝７２）よりも小さく設定される。これにより、変性基粒子モデル６ｂは、未変性粒子モデル６ａに比べて斥力が小さくなる。

このような変性基粒子モデル６ｂは、各フィラー粒子モデル４ｃ、４ｓとの親和性を高く表現できるため、実際の高分子材料に配合される変性剤として再現される。これにより、ポリマーモデル５に変性基粒子モデル６ｂを含ませる（変性ポリマー化する）ことにより、実際の高分子材料と同様に、ポリマーモデル５中のフィラーモデル３の分散状態を変化させることができる。従って、本実施形態のシミュレーション方法では、シミュレーション精度を向上しうる。

さらに、これらのポテンシャルＵ１乃至ポテンシャルＵ１０は、前記式（１）のカットオフ距離ｒ_cが次のように設定される。
ポテンシャルＵ１：ｒ_c＝３
ポテンシャルＵ２：ｒ_c＝３
ポテンシャルＵ３：ｒ_c＝３
ポテンシャルＵ４：ｒ_c＝１
ポテンシャルＵ５：ｒ_c＝１
ポテンシャルＵ６：ｒ_c＝１
ポテンシャルＵ７：ｒ_c＝５
ポテンシャルＵ８：ｒ_c＝１
ポテンシャルＵ９：ｒ_c＝１
ポテンシャルＵ１０：ｒ_c＝１

図７は、フィラー粒子モデルのカットオフ距離を説明する概念図である。フィラーモデル３の中心粒子モデル４ｃが関連するポテンシャル（例えば、ポテンシャルＵ１）の各カットオフ距離ｒ_cは、フィラーモデル３の表面粒子モデル４ｓが関連するポテンシャル（例えば、ポテンシャルＵ４）の各カットオフ距離ｒ_cよりも大に設定されている。さらに、本実施形態の中心粒子モデル４ｃが関連するポテンシャル（例えば、ポテンシャルＵ１）の各カットオフ距離ｒ_cは、表面粒子モデル４ｓが関連するポテンシャル（例えば、ポテンシャルＵ４）の各カットオフ距離ｒ_c、及び、表面粒子モデル４ｓと中心粒子モデル４ｃとの粒子間距離Ｌ_cの和（ｒ_c＋Ｌ_c）よりも大に設定されている。

上記のように、カットオフ距離ｒ_cを定義することにより、フィラーモデル３は、中心粒子モデル４ｃが関連するポテンシャルＵ１〜ポテンシャルＵ３、及び、ポテンシャルＵ７（図５に示す）を、表面粒子モデル４ｓが関連するポテンシャルＵ４、ポテンシャルＵ５、及び、ポテンシャルＵ８よりも優先的に作用させることができる。しかも、中心粒子モデル４ｃは、径を持った球で表現されるため、ポテンシャルＵ１〜ポテンシャルＵ３、及び、ポテンシャルＵ７が放射状に作用する。従って、コンピュータ１は、分子動力学計算において、フィラーモデル３を、実際のフィラーと近似する球として扱えるため、シミュレーション精度を向上しうる。

さらに、コンピュータ１は、表面粒子モデル４ｓのカットオフ距離ｒ_c内にまで浸入する粒子を除いて、実質的に中心粒子モデル４ｃのみの各ポテンシャルＵ１〜ポテンシャルＵ３、及び、ポテンシャルＵ７で、フィラーモデル３を分子動力学計算することができる。従って、本発明のシミュレーション方法では、計算効率を向上させることができる。このようなポテンシャルＵは、数値データとしてコンピュータ１に入力される。

次に、コンピュータ１内で予め定められた体積を持つ仮想空間の中に、複数のフィラーモデル３及びポリマーモデル５が配置される（工程Ｓ３２）。図８は、シミュレーションモデルの仮想空間９の概念斜視図である。

仮想空間９は、解析対象のゴムポリマーの微小構造部分に相当するものである。本実施形態の仮想空間９は、１辺の長さＬ１が、例えば３０［σ］の立方体として定義されている。この仮想空間９には、例えば、１００個のフィラーモデル３、及び、２０００本のポリマーモデル５が、ランダムに配置される。

次に、フィラーモデル３の凝集塊が形成される（工程Ｓ３３）。この工程Ｓ３３では、少なくとも２つのフィラーモデル３を、相互作用が互いに生じる距離（即ち、カットオフ距離ｒ_c内）に配置して、フィラーモデル３の凝集塊が形成される。これにより、仮想空間９には、解析対象のゴムポリマー内に形成されがちなフィラーの凝集塊を再現でき、現実に近いシミュレーションを行うことができる。

次に、図２に示されるように、本実施形態のシミュレーション方法では、コンピュータ１が、仮想空間９に配置されたフィラーモデル３とポリマーモデル５とを用いて、分子動力学（ Molecular Dynamics : ＭＤ）計算を実施する（シミュレーション工程Ｓ４）。図９は、本実施形態のシミュレーション工程Ｓ４の処理手順の一例を示すフローチャートである。

図８及び図９に示されるように、シミュレーション工程Ｓ４では、先ず、各フィラー粒子モデル４ｃ、４ｓの動きが拘束される（フィラー拘束工程Ｓ４１）。このフィラー拘束工程Ｓ４１では、コンピュータ１によって、仮想空間９でのフィラー粒子モデル４の座標が固定され、フィラー粒子モデル４が移動不能に拘束される。

次に、シミュレーション工程Ｓ４では、コンピュータ１が、ポリマーモデル５のみを対象として、分子動力学計算が実施される（工程Ｓ４２）。工程Ｓ４２での分子動力学計算では、例えば、設定された仮想空間９について所定の時間、拘束されたフィラーモデル３を除いた全てのポリマーモデル５が、古典力学に従うものとしてニュートンの運動方程式が適用される。そして、各時刻における各ポリマー粒子モデル６ａ、６ｂの動きが追跡される。なお、分子動力学計算では、各フィラー粒子モデル４ｃ、４ｓの数、各ポリマー粒子モデル６ａ、６ｂの数、仮想空間９の体積、及び、仮想空間の温度等の諸条件が一定に保たれる。

このように、本実施形態の工程Ｓ４２では、ポリマーモデル５のみが仮想空間９内で分散されるため、フィラーモデル３の凝集塊が保たれた上で、ポリマーモデル５の安定配置を求めることができる。

次に、シミュレーション工程Ｓ４では、コンピュータ１が、ポリマーモデル５が十分に分散したか否かを判断する（工程Ｓ４３）。本実施形態では、ポリマーモデル５が十分に分散したと判断された場合（工程Ｓ４３で「Ｙ」）、次の工程Ｓ４４が実施される。一方、ポリマーモデル５が十分に分散していないと判断された場合（工程Ｓ４３で「Ｎ」）は、単位ステップを進めて（工程Ｓ４５）、工程Ｓ４２が再度実施される。従って、本実施形態のシミュレーション工程Ｓ４では、ポリマーモデル５を効果的に分散させることができる。

工程Ｓ４２においてポリマーモデル５の安定配置を求めるために、１ステップあたりの分子動力学計算の時間幅を０．０５［τ］したときのステップ数は、１００ステップ以上が望ましい。なお、ステップ数が１００ステップ未満であると、計算時間が短くなり、ポリマーモデル５を十分に分散させることができない。このため、上記作用を十分に発揮させることができないおそれがある。このような観点より、ステップ数は、より好ましくは１０００ステップ以上である。一方、ステップ数が多すぎても、計算コストの増加分に対して十分な分散結果を得ることができないおそれがある。このような観点より、ステップ数は、好ましくは１００万ステップ以下である。

次に、シミュレーション工程Ｓ４では、フィラーモデル３の各フィラー粒子モデル４の拘束が解除される（工程Ｓ４４）。この工程Ｓ４４では、コンピュータ１によって、フィラー粒子モデル４の座標の固定が解除される。これにより、シミュレーション工程Ｓ４では、後述する工程Ｓ４７において、フィラーモデル３を、仮想空間９内で自由に移動させることができる。

次に、シミュレーション工程Ｓ４では、フィラー粒子モデル４とポリマー粒子モデル６とが結合される（結合工程Ｓ４６）。この結合工程Ｓ４６では、複数のフィラー粒子モデル４から予め定められた結合対象フィラー粒子モデル１５に、複数のポリマー粒子モデル６から予め定められた結合対象ポリマー粒子モデル１６が、連結鎖で結合される。これにより、本実施形態のシミュレーション方法では、フィラーとポリマーとの化学結合を再現することができる。

図３に示されるように、結合対象フィラー粒子モデル１５としては、フィラーモデル３を構成する複数のフィラー粒子モデル４から、任意のフィラー粒子モデル４を選択することができる。本実施形態の結合対象フィラー粒子モデル１５は、多面体１３の頂点１３ａに配置されたフィラー粒子モデル４（頂点粒子モデル７）の８個から構成される。

図４に示されるように、結合対象ポリマー粒子モデル１６としては、ポリマーモデル５を構成する複数のポリマー粒子モデル６から、任意のポリマー粒子モデル６を選択することができる。本実施形態では、直鎖状に配置される複数のポリマー粒子モデル６のうち、中央に配置されるポリマー粒子モデル６が１個選択される。

図１０は、本実施形態の結合工程Ｓ４６の処理手順の一例を示すフローチャートである。本実施形態の結合工程Ｓ４６では、先ず、各結合対象フィラー粒子モデル１５に結合させる結合対象ポリマー粒子モデル１６の数（以下、単に「結合数」ということがある。）Ｓａの上限値Ｓｕが設定される（工程Ｓ１１）。上限値Ｓｕは、例えば、結合対象フィラー粒子モデル１５の数、及び、結合対象ポリマー粒子モデル１６の数に応じて、適宜設定することができる。本実施形態の上限値Ｓｕは、下記式（２）の関係を満たすように定義される。
Ｎａ／Ｎｂ−１≦Ｓｕ≦Ｎａ／Ｎｂ＋１ …（２）
ここで、各パラメータは、次のとおりである。
Ｎａ：結合対象ポリマー粒子モデルの総数
Ｎｂ：結合対象フィラー粒子モデルの総数

総数Ｎａは、仮想空間９（図８に示す）内に配置された結合対象ポリマー粒子モデル１６の総数である。図４に示したように、本実施形態の各ポリマーモデル５には、結合対象ポリマー粒子モデル１６が１個含まれている。従って、総数Ｎａは、仮想空間９内に配置されたポリマーモデル５の総数（例えば、２０００）と同一である。

総数Ｎｂは、仮想空間９内に配置された結合対象フィラー粒子モデル１５の総数である。図３に示したように、本実施形態の各フィラーモデル３には、結合対象フィラー粒子モデル１５が複数（本実施形態では、８個）含まれている。従って、総数Ｎｂは、仮想空間９内に配置されたフィラーモデル３の総数（例えば１００）と、各フィラーモデル３に含まれる結合対象フィラー粒子モデル１５の数（例えば、８）を乗じた値である。

従って、Ｎａ／Ｎｂは、各結合対象フィラー粒子モデル１５に、結合対象ポリマー粒子モデル１６が均等に結合された場合の個数を示している。例えば、結合対象ポリマー粒子モデルの総数Ｎａが２０００、かつ、結合対象フィラー粒子モデルの総数Ｎｂが８００である場合、Ｎａ／Ｎｂは、２．５となる。Ｎａ／Ｎｂが２．５の場合、上記式（２）の関係に基づくと、上限値Ｓｕは、１．５以上かつ３．５以下に限定される。結合数Ｓａは、正の整数に限定される。このため、上限値Ｓｕは、２又は３に設定される。

次に、結合工程Ｓ４６では、コンピュータ１が、結合対象フィラー粒子モデル１５に結合させる結合対象ポリマー粒子モデル１６を特定する（特定工程Ｓ１２）。図１１は、本実施形態の特定工程Ｓ１２の処理手順の一例を示すフローチャートである。

本実施形態の特定工程Ｓ１２は、先ず、各結合対象ポリマー粒子モデル１６が、夫々粒子間距離Ｌｐが最も小さい結合対象フィラー粒子モデル１５に関連付けられる（第１工程Ｓ２１）。図１２は、結合対象ポリマー粒子モデル１６を関連付ける工程を説明する概念図である。

第１工程Ｓ２１では、先ず、仮想空間９（図８に示す）において、各結合対象ポリマー粒子モデル１６から粒子間距離Ｌｐが最も小さい結合対象フィラー粒子モデル１５が探索される。なお、粒子間距離Ｌｐは、各粒子モデル１５、１６の中心間の距離として定義される。

次に、第１工程Ｓ２１では、各結合対象ポリマー粒子モデル１６と、それぞれ探索された結合対象フィラー粒子モデル１５とが関連付けられる。この関連付けに関する情報は、コンピュータ１に記憶される。なお「関連付けに関する情報」には、例えば、結合対象ポリマー粒子モデル１６及び結合対象フィラー粒子モデル１５の識別番号や、仮想空間９での各座標値が含まれる。

図１２に示される例では、結合対象フィラー粒子モデル１５として、第１結合対象フィラー粒子モデル１５ａ、及び、第２結合対象フィラー粒子モデル１５ｂが配置されている。第１結合対象フィラー粒子モデル１５ａには、３個の結合対象ポリマー粒子モデル１６が関連付けられている。本例では、３個の結合対象ポリマー粒子モデル１６が、結合対象フィラー粒子モデル１５からの粒子間距離Ｌｐが小さい順に、第１結合対象ポリマー粒子モデル１６ａ、第２結合対象ポリマー粒子モデル１６ｂ及び第３結合対象ポリマー粒子モデル１６ｃとして区別される。

一方、第２結合対象フィラー粒子モデル１５ｂには、１個の結合対象ポリマー粒子モデル１６が関連付けられている。本例では、１個の結合対象ポリマー粒子モデル１６が、第４結合対象ポリマー粒子モデル１６ｄとして区別される。図１２において、結合対象フィラー粒子モデル１５と、結合対象ポリマー粒子モデル１６との関連付けが、二点鎖線で示されている。

次に、本実施形態の特定工程Ｓ１２は、各結合対象フィラー粒子モデル１５に結合させる結合対象ポリマー粒子モデル１６が、関連付けられた結合対象ポリマー粒子モデル１６の集団（以下、単に「結合対象ポリマー粒子モデル群」ということがある。）から特定される（第２工程Ｓ２２）。この第２工程Ｓ２では、前記結合数Ｓａが上限値Ｓｕ以下を満たすように、結合対象ポリマー粒子モデル１６を、結合対象ポリマー粒子モデル群から粒子間距離Ｌｐが小さい順に特定している。

図１３は、結合対象ポリマー粒子モデル１６を特定する工程を説明する概念図である。本例での上限値Ｓｕには、「２」が設定されている。

図１２に示したように、第１結合対象フィラー粒子モデル１５ａには、上限値Ｓｕ（本例では、２個）よりも多い３個の結合対象ポリマー粒子モデル１６（結合対象ポリマー粒子モデル群）が関連付けられている。このため、本例では、第１結合対象フィラー粒子モデル１５ａに関連付けられた結合対象ポリマー粒子モデル群のうち、粒子間距離Ｌｐが小さい順に、上限値Ｓｕ（本例では、２個）分の結合対象ポリマー粒子モデル１６ａ、１６ｂが選択される。

図１３に示されるように、選択された第１結合対象ポリマー粒子モデル１６ａ及び第２結合対象ポリマー粒子モデル１６ｂは、第１結合対象フィラー粒子モデル１５ａに結合されるものとして特定される。なお、図１３において、結合対象フィラー粒子モデル１５と結合対象ポリマー粒子モデル１６との特定は、実線で示されている。

一方、結合対象ポリマー粒子モデル群のうち、粒子間距離Ｌｐが最も大きい第３結合対象ポリマー粒子モデル１６ｃは、第１結合対象フィラー粒子モデル１５ａに結合されるものとして特定されない。

第２結合対象フィラー粒子モデル１５ｂには、上限値Ｓｕ（本例では、２個）よりも小さい１個の結合対象ポリマー粒子モデル１６（第４結合対象ポリマー粒子モデル１６ｄ）が関連付けられている。このため、第４結合対象ポリマー粒子モデル１６ｄが、第２結合対象フィラー粒子モデル１５ｂに結合されるものとして特定される。

このように、第２工程Ｓ２２では、結合対象フィラー粒子モデル１５に結合させる結合対象ポリマー粒子モデル１６の結合数Ｓａが上限値Ｓｕ以下となるように、結合対象ポリマー粒子モデル１６が特定される。なお、結合対象ポリマー粒子モデル１６を特定する情報は、例えば、結合対象ポリマー粒子モデル１６及び結合対象フィラー粒子モデル１５の識別番号や、仮想空間９での各座標値が含まれる。結合対象ポリマー粒子モデル１６を特定する情報は、コンピュータ１に記憶される。

次に、本実施形態の特定工程Ｓ１２は、第２工程Ｓ２２で特定されなかった各結合対象ポリマー粒子モデル１６が、新たな結合対象フィラー粒子モデル１５に関連付けられる（第３工程Ｓ２３）。図１４は、特定されなかった結合対象ポリマー粒子モデル１６の関連付けを説明する概念図である。

第３工程Ｓ２３では、先ず、仮想空間９（図９に示す）において、第２工程Ｓ２２で特定されなかった各結合対象ポリマー粒子モデル１６から、第１工程Ｓ２１で関連付けられた結合対象フィラー粒子モデル１５の次に粒子間距離Ｌｐが小さい結合対象フィラー粒子モデル１５が探索される。そして、各結合対象ポリマー粒子モデル１６と、探索された結合対象フィラー粒子モデル１５とが関連付けられる。関連付けに関する情報は、コンピュータ１に記憶される。

図１４に示される例において、第２工程Ｓ２２で特定されなかった結合対象ポリマー粒子モデル１６は、第３結合対象ポリマー粒子モデル１６ｃである。第３工程Ｓ２３では、第３結合対象ポリマー粒子モデル１６ｃとの粒子間距離Ｌｐが、第１結合対象フィラー粒子モデル１５ａの次に小さい結合対象フィラー粒子モデル１５（本例では、第２結合対象フィラー粒子モデル１５ｂ）が探索される。そして、第３結合対象ポリマー粒子モデル１６ｃと、第２結合対象フィラー粒子モデル１５ｂとが関連付けられる。

次に、本実施形態の特定工程Ｓ１２は、各結合対象フィラー粒子モデル１５に結合させる結合対象ポリマー粒子モデル１６が、第３工程Ｓ２３で関連付けられた結合対象ポリマー粒子モデル群から特定される（第４工程Ｓ２４）。第４工程Ｓ２４では、第２工程Ｓ２２と同様に、結合対象ポリマー粒子モデル１６を、前記結合数Ｓａが上限値Ｓｕ以下を満たすように、第３工程Ｓ２３で関連付けられた結合対象ポリマー粒子モデル群から粒子間距離Ｌｐが小さい順に特定している。

図１５は、第３工程Ｓ２３で関連付けられた結合対象ポリマー粒子モデル１６を特定する工程を説明する概念図である。第２結合対象フィラー粒子モデル１５ｂには、第３結合対象ポリマー粒子モデル１６ｃのみが関連付けられている。第２結合対象フィラー粒子モデル１５ｂは、１個の結合対象ポリマー粒子モデル１６（第４結合対象ポリマー粒子モデル１６ｄ）が既に特定されている。このため、第２結合対象フィラー粒子モデル１５ｂでは、結合数Ｓａが上限値Ｓｕ以下を満たすために、関連付けられた１個の結合対象ポリマー粒子モデル１６のみを特定することができる。従って、第３結合対象ポリマー粒子モデル１６ｃは、第２結合対象フィラー粒子モデル１５ｂに結合されるものとして特定される。

なお、第２結合対象フィラー粒子モデル１５ｂに、２個以上の結合対象ポリマー粒子モデル１６が関連付けられている場合は、結合数Ｓａが上限値Ｓｕ以下を満たすように、結合対象ポリマー粒子モデル群から粒子間距離Ｌｐが小さい順に特定される。従って、粒子間距離Ｌｐが２番目以降に小さい結合対象ポリマー粒子モデル１６（図示省略）は、第２結合対象フィラー粒子モデル１５ｂに結合されるものとして特定されない。

このように、第４工程Ｓ２４では、結合対象フィラー粒子モデル１５に結合させる結合対象ポリマー粒子モデル１６を、結合数Ｓａが上限値Ｓｕ以下を満たすように、第３工程Ｓ２３で関連付けられた結合対象ポリマー粒子モデル１６から特定することができる。結合対象ポリマー粒子モデル１６を特定する情報は、コンピュータ１に記憶される。

次に、本実施形態の特定工程Ｓ１２では、全ての結合対象ポリマー粒子モデル１６が、結合対象フィラー粒子モデル１５に特定されたか否かが判断される（工程Ｓ２５）。この工程Ｓ２５では、全ての結合対象ポリマー粒子モデル１６が、結合対象フィラー粒子モデル１５に特定されたと判断された場合、次の工程Ｓ１３が実施される。一方、各結合対象フィラー粒子モデル１５に特定されていない結合対象ポリマー粒子モデル１６が存在すると判断された場合は、特定されていない結合対象ポリマー粒子モデル１６を対象に、第３工程Ｓ２３及び第４工程Ｓ２４が実施される。これにより、特定工程Ｓ１２では、全ての結合対象ポリマー粒子モデル１６を、結合対象フィラー粒子モデル１５に特定することができる。なお、本実施形態では、結合対象ポリマー粒子モデル１６が特定されていない結合対象フィラー粒子モデル１５の存在が許容されている。

次に、結合工程Ｓ４６では、特定工程Ｓ１２で特定された結合対象フィラー粒子モデル１５と結合対象ポリマー粒子モデル１６とが結合される（工程Ｓ１３）。図１６は、結合対象フィラー粒子モデル１５と結合対象ポリマー粒子モデル１６とを結合する工程を説明する概念図である。工程Ｓ１３では、特定された結合対象フィラー粒子モデル１５と結合対象ポリマー粒子モデル１６とを連結鎖１０で結合させている。これにより、フィラーモデル３及びポリマーモデル５は、フィラーとポリマーとの化学結合を再現することができる。

次に、シミュレーション工程Ｓ４では、フィラーモデル３及びポリマーモデル５を対象に分子動力学計算が行われる（工程Ｓ４７）。図１７は、分子動力学計算後のシミュレーションモデルの概念斜視図である。図１７では、ポリマーモデル５が省略されている。

工程Ｓ４７での分子動力学計算は、ポリマー計算工程Ｓ４２と同様に、ニュートンの運動方程式が適用される。そして、各時刻におけるフィラーモデル３の各フィラー粒子モデル４ｃ、４ｓ（図３に示す）、及び、ポリマーモデル５の各ポリマー粒子モデル６ａ、６ｂ（図４に示す）の動きが追跡される。なお、分子動力学計算では、ポリマー計算工程Ｓ４２と同様に、各フィラー粒子モデル４ｃ、４ｓの数、各ポリマー粒子モデル６ａ、６ｂの数、仮想空間９の体積、及び、仮想空間９の温度等の諸条件が一定に保たれる。これにより、図１７に示されるように、工程Ｓ４７では、フィラーモデル３及びポリマーモデル５（図示省略）を、仮想空間９内でそれぞれ分散させることができる。

本実施形態のシミュレーション工程Ｓ４では、結合工程Ｓ４６において、結合対象フィラー粒子モデルに結合させる結合対象ポリマー粒子モデル１６の結合数Ｓａが、予め定められた上限値Ｓｕ以下に設定されている。このため、フィラー粒子モデル４に結合されるポリマー粒子モデル６の数に大きな偏りが生じるのを防ぐことができる。これにより、本実施形態のシミュレーション工程Ｓ４では、分子動力学計算において、結合対象フィラー粒子モデル１５に作用する結合対象ポリマー粒子モデル１６の力の総和が大きくなるのを防ぐことができる。従って、本実施形態のシミュレーション方法は、計算落ちが発生するのを防ぐことができ、分子動力学計算を安定して実行することができる。

また、本実施形態の上限値Ｓｕは、各結合対象フィラー粒子モデル１５に、結合対象ポリマー粒子モデル１６が均等に結合された場合の個数に基づいて、上記式（２）の関係を満たすように設定される。このため、結合対象フィラー粒子モデル１５に、結合対象ポリマー粒子モデル１６を均一に（満遍なく）結合させることができるため、結合対象フィラー粒子モデル１５に作用する結合対象ポリマー粒子モデル１６の力の総和が大きくなるのを確実に防ぐことができる。

さらに、本実施形態の結合対象フィラー粒子モデル１５は、多面体１３の頂点１３ａに配置された頂点粒子モデル７に限定されるため、隣接する結合対象フィラー粒子モデル１５、１５にそれぞれ結合される結合対象ポリマー粒子モデル１６を離間させることができる。これにより、結合対象ポリマー粒子モデル１６、１６間に作用する相互作用が大きくなるのを防ぐことができるため、分子動力学計算を安定して実行することができる。

次に、コンピュータ１が、フィラーモデル３及びポリマーモデル５の十分に分散したか否かを判断する（工程Ｓ４８）。本実施形態の工程Ｓ４８では、フィラーモデル３及びポリマーモデル５が十分に分散したと判断された場合（工程Ｓ４８で「Ｙ」）、次の評価工程Ｓ５（図２に示す）が実施される。一方、フィラーモデル３及びポリマーモデル５が十分に分散していないと判断された場合（工程Ｓ４８で「Ｎ」）は、単位ステップを進めて（工程Ｓ４９）、工程Ｓ４７が再度実施される。従って、本実施形態のシミュレーション工程Ｓ４では、フィラーモデル３及びポリマーモデル５を効果的に分散させることができ、構造緩和された高分子材料モデル１１を設定することができる。

なお、フィラーモデル３及びポリマーモデル５を効果的に分散させるために、工程Ｓ４７において、１ステップあたりの分子動力学計算の時間幅を０．０５［τ］したときのステップ数は、１０００ステップ以上が望ましい。なお、ステップ数が１０００ステップ未満であると、フィラーモデル３及びポリマーモデル５を十分に分散できないおそれがある。このような観点より、ステップ数は、より好ましくは５０００ステップ以上である。一方、ステップ数が多すぎても、計算コストの増加分に対して十分な分散効果を得ることができないおそれがある。このような観点より、ステップ数は、好ましくは１００万ステップ以下である。

次に、図２に示されるように、本実施形態のシミュレーション方法では、フィラーモデル３の分散状態が良好か否か判断される（評価工程Ｓ５）。この評価工程Ｓ５では、上記特許文献１と同様に、フィラーモデル３の中心粒子モデル４ｃを対象とした動径分布関数が計算される。動径分布関数とは、ある中心粒子モデル４ｃ（図３に示す）から距離ｒ離れた位置において、他の中心粒子モデル４ｃが存在する確率密度を表す関数である。距離ｒは、中心粒子モデル４ｃの中心間の距離として定義される。

このように、本実施形態では、中心粒子モデル４ｃのみを対象に動径分布関数を求めることより、フィラーモデル３の分散状態を確認することができるため、計算コストの増大を抑制することができる。

そして、評価工程Ｓ５では、動径分布関数の結果をもとに、コンピュータ１が予め設定された許容範囲内であるか否かが判断される。評価工程Ｓ５では、フィラーモデル３の分散状態が良好であると判断された場合、高分子材料モデル１１に基づいて、高分子材料が製造される（工程Ｓ６）。一方、フィラーモデル３の分散状態が良好でないと判断された場合は、動径分布関数の結果に基づいて、例えば、フィラーモデル３及びポリマーモデル５の設定条件を変更して（工程Ｓ７）、工程Ｓ１〜工程Ｓ５が実施される。これにより、本実施形態のシミュレーション方法では、高分子材料モデル１１（図１７に示す）の構造を、効果的に緩和することができ、実際の高分子材料の構造に、確実に近似させることができる。

また、本実施形態のシミュレーション方法では、フィラーモデル３及びポリマーモデル５の構造やポテンシャルを任意に設定することができるため、現実には存在しない未知の高分子材料をモデル化した高分子材料モデル１１を評価することができる。

本実施形態のシミュレーション方法では、仮想空間９に配置されたフィラーモデル３とポリマーモデル５を用いて、高分子材料中に配合されたフィラーの分散性が評価されるものが例示されたが、これに限定されるわけではない。例えば、本発明のシミュレーション方法は、フィラーと、フィラーをポリマーに結合させるためのカップリング剤とが配合された高分子材料について、フィラーの分散性が評価されるものでもよい。

図１８は、本発明の他の実施形態のシミュレーション方法の処理手順の一例を示すフローチャートである。なお、本明細書において、前実施形態と同一処理が実施される工程、及び、同一のモデル等については、前実施形態と同一の符号を付して、説明を省略することがある。

この実施形態のシミュレーション方法では、先ず、前実施形態と同様に、フィラーをモデル化したフィラーモデル３を設定する工程Ｓ１、及び、高分子材料をモデル化したポリマーモデル５を設定する工程Ｓ２が実施される。

次に、この実施形態のシミュレーション方法では、カップリング剤をモデル化したカップリングモデルが設定される（工程Ｓ８）。図１９は、カップリングモデル１７の概念図である。

カップリングモデル１７は、少なくとも一つ、本実施形態では、複数のカップリング粒子モデル１８を用いて、カップリング剤がモデル化される。カップリング粒子モデル１８は、径を持った球で表現されている。カップリング粒子モデル１８の個数は、ポリマー粒子モデル６の個数よりも小に設定される。これにより、カップリングモデル１７は、高分子材料の分子鎖よりも分子量が小さい実際のカップリング剤が表現されうる。

カップリング粒子モデル１８、１８間には、平衡長が定義された結合鎖２１が定義されている。これにより、本実施形態のカップリングモデル１７は、ポリマーモデル５（図４に示す）と同様に、直鎖状の三次元構造に形成される。なお、一つのカップリング粒子モデル１８で構成されるカップリングモデル１７には、結合鎖２１が定義されない。

本実施形態のカップリングモデル１７は、同一個数のカップリング粒子モデル１８によってモデル化されているが、これに限定されるわけではない。例えば、カップリング剤の鎖長の差異によって生じるフィラーの分散性の影響を評価する場合は、カップリング粒子モデル１８の個数が異なる複数種類のカップリングモデル１７が定義されても良い。

カップリングモデル１７は、コンピュータ１で取り扱い可能な数値データであり、コンピュータ１に入力される。数値データには、例えば、カップリング粒子モデル１８の質量、体積、直径、及び、初期座標などが含まれる。

次に、この実施形態のシミュレーション方法では、シミュレーション条件が設定される（条件設定工程Ｓ３）。図２０は、本発明の他の実施形態の条件設定工程Ｓ３の処理手順の一例を示すフローチャートである。

この実施形態の条件設定工程Ｓ３では、先ず、前実施形態と同様に、各フィラー粒子モデル４ｃ、４ｓと各ポリマー粒子モデル６ａ、６ｂとの間に、ポテンシャルが定義される（工程Ｓ３１）。図６に示したように、フィラー粒子モデル４ｃ、４ｓと、各ポリマー粒子モデル６ａ、６ｂとの間に定義されるポテンシャルＵ１乃至ポテンシャルＵ１０については、前実施形態のポテンシャルと同一である。

次に、この実施形態の条件設定工程Ｓ３では、各カップリング粒子モデル１８と、各フィラー粒子モデル４ｃ、４ｓとの間、及び、各カップリング粒子モデル１８と各ポリマー粒子モデル６ａ、６ｂとの間に、ポテンシャルが定義される（工程Ｓ３４）。ポテンシャルは、上記式（１）で定義される。図２１は、カップリング粒子モデル１８のポテンシャルを説明する概念図である。

本実施形態では、カップリング粒子モデル１８と各フィラー粒子モデル４ｃ、４ｓとの間、及び、カップリング粒子モデル１８と各ポリマー粒子モデル６ａ、６ｂとの間に、下記のポテンシャルＵ１１乃至ポテンシャルＵ１４が定義される。
ポテンシャルＵ１１：カップリング粒子モデル１８−中心粒子モデル４ｃ
ポテンシャルＵ１２：カップリング粒子モデル１８−表面粒子モデル４ｓ
ポテンシャルＵ１３：カップリング粒子モデル１８−未変性粒子モデル６ａ
ポテンシャルＵ１４：カップリング粒子モデル１８−変性基粒子モデル６ｂ

本実施形態のポテンシャルＵ１１乃至ポテンシャルＵ１４の各強度ａ_ijは、上記論文のパラメータを参考として、次のように設定される。
ポテンシャルＵ１１：ａ_ij＝５０
ポテンシャルＵ１２：ａ_ij＝５０
ポテンシャルＵ１３：ａ_ij＝５０
ポテンシャルＵ１４：ａ_ij＝５０

ポテンシャルＵ１１〜ポテンシャルＵ１４では、同一の強度ａ_ijが設定されている。これにより、フィラーモデル３に対するカップリングモデル１７の斥力と、ポリマーモデル５に対するカップリングモデル１７の斥力とが同一に設定されている。従って、ポテンシャルＵ１１〜ポテンシャルＵ１４は、フィラーに対する親和性と、ポリマーに対する親和性が同一であるカップリング剤を再現することができる。

本実施形態のポテンシャルＵ１１〜ポテンシャルＵ１４は、同一の強度ａ_ijが定義されるものが例示されたが、これに限定されるわけではない。例えば、下記のように、異なる強度ａ_ijが設定されても良い。これにより、フィラーに対する親和性と、ポリマーに対する親和性とを異ならせたカップリング剤が再現されうる。
ポテンシャルＵ１１：ａ_ij＝１０
ポテンシャルＵ１２：ａ_ij＝１０
ポテンシャルＵ１３：ａ_ij＝１００
ポテンシャルＵ１４：ａ_ij＝１００

上記のように、フィラーモデル３の各粒子モデル４ｃ、４ｓが関連するポテンシャルＵ１１及びポテンシャルＵ１２の強度ａ_ij（＝１０）は、ポリマーモデル５の各粒子モデル６ａ、６ｂが関連するポテンシャルＵ１３及びポテンシャルＵ１４の強度ａ_ij（＝１００）よりも小に設定される。これにより、フィラーモデル３とカップリングモデル１７との斥力は、ポリマーモデル５とカップリングモデル１７との斥力に比べて小さくなる。このようなカップリングモデル１７は、フィラーモデル３との親和性を、ポリマーモデルとの親和性よりも高く設定することができるため、例えば、カップリング剤の相互作用の差異によって生じるフィラーの分散性の影響を評価することができる。

さらに、ポテンシャルＵ１１乃至ポテンシャルＵ１４は、上記式（１）のカットオフ距離ｒ_cが次のように設定される。
ポテンシャルＵ１１：ｒ_c＝３
ポテンシャルＵ１２：ｒ_c＝１
ポテンシャルＵ１３：ｒ_c＝１
ポテンシャルＵ１４：ｒ_c＝１

フィラーモデル３の中心粒子モデル４ｃが関連するポテンシャルＵ１１のカットオフ距離ｒ_cは、フィラーモデル３の表面粒子モデル４ｓが関連するポテンシャルＵ１２のカットオフ距離ｒ_cよりも大に設定されている。さらに、ポテンシャルＵ１１のカットオフ距離ｒ_cは、ポテンシャルＵ１２のカットオフ距離ｒ_c、及び、図７に示した表面粒子モデル４ｓと中心粒子モデル４ｃとの粒子間距離Ｌ_cの和（ｒ_c＋Ｌ_c）よりも大に設定されている。

このようなカットオフ距離ｒ_cが定義されることにより、フィラーモデル３は、中心粒子モデル４ｃが関連するポテンシャルＵ１１（図２１に示す）を、表面粒子モデル４ｓが関連するポテンシャルＵ１２よりも優先的に作用させることができる。しかも、中心粒子モデル４ｃは、径を持った球で表現されるため、ポテンシャルＵ１１を放射状に作用させることができる。従って、シミュレーション工程Ｓ４では、フィラーモデル３を、実際のフィラーと近似する球として扱えるため、シミュレーション精度を向上しうる。

次に、この実施形態の条件設定工程Ｓ３では、予め定められた仮想空間９の中に、複数のフィラーモデル３、ポリマーモデル５、及び、カップリングモデル１７が配置され（工程Ｓ３２）、フィラーモデル３の凝集塊が形成される（工程Ｓ３３）。図２２は、本発明の他の実施形態のシミュレーションモデルの仮想空間９の概念斜視図である。

仮想空間９は、前実施形態の仮想空間（図８に示す）と同一のものが用いられる。この仮想空間９には、例えば、１００個のフィラーモデル３、１０００本のポリマーモデル５、及び、１０００個のカップリングモデル１７が、ランダムに配置される。また、フィラーモデル３の凝集塊の形成方法については、前実施形態と同一の方法が採用される。

次に、この実施形態のシミュレーション方法では、コンピュータ１が、仮想空間９に配置されたフィラーモデル３とポリマーモデル５とを用いて、分子動力学（ Molecular Dynamics : ＭＤ）計算を実施する（シミュレーション工程Ｓ４）。図２３は、本発明の他の実施形態のシミュレーション工程Ｓ４の処理手順の一例を示すフローチャートである。

この実施形態のシミュレーション工程Ｓ４では、先ず、前実施形態と同様に、各フィラー粒子モデル４ｃ、４ｓ（図２２に示す）の動きが拘束される（フィラー拘束工程Ｓ４１）。次に、この実施形態のシミュレーション工程Ｓ４では、ポリマーモデル５及びカップリングモデル１７を対象として、分子動力学計算が実施される（工程Ｓ５１）。

工程Ｓ５１での分子動力学計算では、例えば、設定された仮想空間９について所定の時間、拘束されたフィラーモデル３を除いた全てのポリマーモデル５、及び、カップリングモデル１７が、古典力学に従うものとしてニュートンの運動方程式が適用される。そして、各時刻における各ポリマー粒子モデル６ａ、６ｂ及びカップリング粒子モデル１８の動きが追跡される。なお、分子動力学計算では、各フィラー粒子モデル４ｃ、４ｓの数、各ポリマー粒子モデル６ａ、６ｂの数、カップリング粒子モデル１８の数、仮想空間９の体積、及び、仮想空間の温度等の諸条件が一定に保たれる。

このように、本実施形態の工程Ｓ５１では、ポリマーモデル５及びカップリングモデル１７のみが、仮想空間９内で分散されるため、フィラーモデル３の凝集塊が保たれた上で、ポリマーモデル５及びカップリングモデル１７の安定配置を求めることができる。

次に、この実施形態のシミュレーション工程Ｓ４では、コンピュータ１が、ポリマーモデル５及びカップリングモデル１７が十分に分散したか否かを判断する（工程Ｓ５２）。本実施形態では、ポリマーモデル５及びカップリングモデル１７が十分に分散したと判断された場合（工程Ｓ５２で「Ｙ」）、次の工程Ｓ４４が実施される。一方、ポリマーモデル５及びカップリングモデル１７が十分に分散していないと判断された場合（工程Ｓ５２で「Ｎ」）は、単位ステップを進めて（工程Ｓ４５）、工程Ｓ５１が再度実施される。従って、本実施形態のシミュレーション工程Ｓ４では、ポリマーモデル５及びカップリングモデル１７を効果的に分散させることができる。なお、ポリマーモデル５及びカップリングモデル１７の安定配置を求めるために、１ステップあたりの分子動力学計算の時間幅を０．０５［τ］したときのステップ数は、前実施形態のポリマー計算工程Ｓ４２と同一範囲が望ましい。

次に、この実施形態のシミュレーション工程Ｓ４では、フィラーモデル３の各フィラー粒子モデル４の拘束が解除される（工程Ｓ４４）。フィラー粒子モデル４の拘束の解除方法については、前実施形態と同一である。

次に、この実施形態のシミュレーション工程Ｓ４は、ポリマー粒子モデル６とカップリング粒子モデル１８とが結合される（工程Ｓ５３）。図２４（ａ）は、ポリマー粒子モデル６とカップリング粒子モデル１８とが結合される前の状態を示す概念図、図２４（ｂ）は、ポリマー粒子モデル６とカップリング粒子モデル１８とが結合された状態を示す概念図である。

本実施形態の工程Ｓ５３では、予め定められた粒子間距離Ｌ３以下に接近したポリマー粒子モデル６及びカップリング粒子モデル１８が、連結鎖２０で結合される。これにより、ポリマーとカップリング剤との化学結合が、シミュレーション工程Ｓ４において、再現することができる。粒子間距離Ｌ３は、カップリング粒子モデル１８が関連するポテンシャルの各カットオフ距離ｒ_c（図７に示す）の０％〜２００％が望ましい。

各カップリングモデル１７において、粒子間距離Ｌ３以下に接近したポリマー粒子モデル６が複数存在する場合は、粒子間距離Ｌｐが最も小さいポリマー粒子モデル６のみと結合される。これにより、各カップリングモデル１７には、複数のポリマーモデル５が結合されるのを防ぐことができるため、カップリング粒子モデル１８に作用するポリマー粒子モデル６の力の総和が大きくなるのを防ぐことができるため、計算落ちが発生するのを防ぐことができる。

次に、この実施形態のシミュレーション工程Ｓ４は、フィラー粒子モデル４に、ポリマー粒子モデル６又はカップリング粒子モデルが結合される（結合工程Ｓ５４）。この実施形態の結合工程Ｓ５４では、複数のフィラー粒子モデル４から予め定められた結合対象フィラー粒子モデル１５に、複数のポリマー粒子モデル６から予め定められた結合対象ポリマー粒子モデル１６、又は、複数のカップリング粒子モデル１８から予め定められた結合対象カップリング粒子モデル１９が、連結鎖で結合される。これにより、この実施形態のシミュレーション工程Ｓ４では、フィラー、ポリマー及びカップリング剤の化学結合を再現することができる。

図３に示されるように、結合対象フィラー粒子モデル１５としては、任意のフィラー粒子モデル４から選択することができる。この実施形態の結合対象フィラー粒子モデル１５は、前実施形態と同様に、フィラー粒子モデル４（頂点粒子モデル７）の８個から構成される。

図４に示されるように、結合対象ポリマー粒子モデル１６としては、任意のポリマー粒子モデル６を選択することができる。この実施形態の結合対象ポリマー粒子モデル１６は、前実施形態と同様に、直鎖状に配置される複数のポリマー粒子モデル６のうち、中央に配置されるポリマー粒子モデル６が１個選択される。

図１９に示されるように、この実施形態の結合対象カップリング粒子モデル１９は、直鎖状に配置される複数のカップリング粒子モデル１８のうち、一端に配置されるカップリング粒子モデル１８が１個選択される。

図２５は、本発明の他の実施形態の結合工程Ｓ５４の処理手順の一例を示すフローチャートである。結合工程Ｓ５４では、先ず、各結合対象フィラー粒子モデル１５に結合させる結合対象ポリマー粒子モデル１６の数（以下、単に「結合数」ということがある。）Ｓａと、各結合対象フィラー粒子モデル１５に結合させる結合対象カップリング粒子モデル１９の数（以下、単に「結合数」ということがある。）Ｓｂとの和（Ｓａ＋Ｓｂ）の上限値Ｎｓが設定される（工程Ｓ６１）。上限値Ｎｓは、例えば、結合対象フィラー粒子モデル１５の数、結合対象ポリマー粒子モデル１６の数、及び、結合対象カップリング粒子モデル１９の数に応じて、適宜設定することができる。本実施形態の上限値Ｎｓは、下記式（３）の関係を満たすように定義される。
Ｎｓ＝［（Ｎａ＋Ｎｃ）／Ｎｂ＋１］ …（３）
ここで、各パラメータは、次のとおりである。
Ｎａ：結合対象ポリマー粒子モデルの総数
Ｎｂ：結合対象フィラー粒子モデルの総数
Ｎｃ：結合対象カップリング粒子モデルの総数

総数Ｎａは、前実施形態の上限値Ｓｕを示した上記式（２）と同様に、仮想空間９（図２２に示す）内に配置された結合対象ポリマー粒子モデル１６（図４に示す）の総数である。従って、総数Ｎａは、仮想空間９内に配置されたポリマーモデル５の総数（例えば、１０００）と同一である。

総数Ｎｂは、上記式（２）と同様に、仮想空間９（図２２に示す）内に配置された結合対象フィラー粒子モデル１５（図３に示す）の総数である。従って、総数Ｎｂは、フィラーモデル３の総数（例えば１００）と、各フィラーモデル３に含まれる結合対象フィラー粒子モデル１５の数（例えば、８）を乗じた値である。

総数Ｎｃは、仮想空間９（図２２に示す）内に配置された結合対象カップリング粒子モデル１９の総数である。図１９に示したように、本実施形態の各カップリングモデル１７には、結合対象カップリング粒子モデル１９が１個含まれている。従って、総数Ｎｃは、仮想空間９内に配置されたカップリングモデル１７の総数（例えば、１０００）と同一である。

従って、上記式（３）において、（Ｎａ＋Ｎｃ）／Ｎｂは、各結合対象フィラー粒子モデル１５に、結合対象ポリマー粒子モデル１６又は結合対象カップリング粒子モデル１９が均等に結合された場合の個数を示している。また、上記式（３）において、［（Ｎａ＋Ｎｃ）／Ｎｂ＋１］は、（Ｎａ＋Ｎｃ）／Ｎｂを切り上げた整数を示している。例えば、結合対象ポリマー粒子モデルの総数Ｎａが１０００、結合対象フィラー粒子モデルの総数Ｎｂが８００、及び、結合対象カップリング粒子モデルの総数Ｎｃが１０００である場合、Ｎａ＋Ｎｃ）／Ｎｂは、２．５である。また、上記式（３）において、［（Ｎａ＋Ｎｃ）／Ｎｂ＋１］では、２．５が切り上げられ、３となる。従って、上限値Ｎｓは、３が設定される。

次に、この実施形態の結合工程Ｓ５４では、各結合対象フィラー粒子モデル１５に結合させる結合対象ポリマー粒子モデル１６、又は、各結合対象フィラー粒子モデル１５に結合させる結合対象カップリング粒子モデル１９が特定される（特定工程Ｓ６２）。図２６は、本発明の実施形態の特定工程Ｓ６２の処理手順の一例を示すフローチャートである。

本実施形態の特定工程Ｓ６２は、先ず、各結合対象ポリマー粒子モデル１６及び各結合対象カップリング粒子モデル１９が、夫々粒子間距離Ｌｐが最も小さい結合対象フィラー粒子モデル１５に関連付けられる（第１工程Ｓ７１）。図２７は、結合対象ポリマー粒子モデル１６又は結合対象カップリング粒子モデル１９を関連付ける工程を説明する概念図である。図２７及び後述する図２８〜図３１では、カップリングモデル１７に結合されたポリマーモデル５を省略して表示している。

第１工程Ｓ７１では、先ず、仮想空間９（図２２に示す）において、各結合対象ポリマー粒子モデル１６から粒子間距離Ｌｐが最も小さい結合対象フィラー粒子モデル１５が探索される。そして、第１工程Ｓ７１では、結合対象ポリマー粒子モデル１６及び結合対象カップリング粒子モデル１９が、それぞれ探索された結合対象フィラー粒子モデル１５に関連付けられる。

次に、第１工程Ｓ７１では、仮想空間９（図２２に示す）において、各結合対象カップリング粒子モデル１９から粒子間距離Ｌｐが最も小さい結合対象フィラー粒子モデル１５が探索される。そして、第１工程Ｓ７１では、各結合対象カップリング粒子モデル１９と、それぞれ探索された結合対象フィラー粒子モデル１５とが関連付けられる。なお、粒子間距離Ｌｐは、結合対象フィラー粒子モデル１５、結合対象ポリマー粒子モデル１６又は結合対象カップリング粒子モデル１９の中心間の距離として定義される。また、「関連付けに関する情報」には、前実施形態と同一の情報が含まれる。

図２７に示される例では、結合対象フィラー粒子モデル１５として、第１結合対象フィラー粒子モデル１５ａ、及び、第２結合対象フィラー粒子モデル１５ｂが配置されている。

第１結合対象フィラー粒子モデル１５ａには、３個の結合対象ポリマー粒子モデル１６及び２個の結合対象カップリング粒子モデル１９が関連付けられている。３個の結合対象ポリマー粒子モデル１６は、結合対象フィラー粒子モデル１５からの粒子間距離Ｌｐが小さい順に、第１結合対象ポリマー粒子モデル１６ａ、第２結合対象ポリマー粒子モデル１６ｂ、及び、第３結合対象ポリマー粒子モデル１６ｃとして区別される。また、２個の結合対象カップリング粒子モデル１９は、結合対象フィラー粒子モデル１５からの粒子間距離Ｌｐが小さい順に、第１結合対象カップリング粒子モデル１９ａ、及び、第２結合対象カップリング粒子モデル１９ｂとして区分される。

第２結合対象フィラー粒子モデル１５ｂには、１個の結合対象ポリマー粒子モデル１６が関連付けられている。１個の結合対象ポリマー粒子モデル１６は、第４結合対象ポリマー粒子モデル１６ｄとして区別される。図２７において、結合対象フィラー粒子モデル１５と、結合対象ポリマー粒子モデル１６又は結合対象カップリング粒子モデル１９との関連付けが、二点鎖線で示されている。

次に、本実施形態の特定工程Ｓ６２は、各結合対象フィラー粒子モデル１５に結合させる結合対象ポリマー粒子モデル１６又は結合対象カップリング粒子モデル１９が、関連付けられた結合対象ポリマー粒子モデル１６又は結合対象カップリング粒子モデル１９から特定される（第２工程Ｓ７２）。

この第２工程Ｓ７２では、結合対象ポリマー粒子モデル１６の結合数Ｓａと、結合対象カップリング粒子モデル１９の結合数Ｓｂとの和（Ｓａ＋Ｓｂ）が上限値Ｎｓ以下を満たすように、関連付けられた結合対象ポリマー粒子モデル１６又は結合対象カップリング粒子モデル１９から粒子間距離Ｌｐが小さい順に特定される。図２８は、結合対象ポリマー粒子モデル１６又は結合対象カップリング粒子モデル１９を特定する工程を説明する概念図である。本例での上限値Ｎｓには、「３」が設定されている。

図２７に示したように、第１結合対象フィラー粒子モデル１５ａに関連付けられている結合対象ポリマー粒子モデル１６の数と、結合対象カップリング粒子モデル１９の数との和（本例では、５個）は、上限値Ｎｓ（本例では、３個）よりも大きい。このため、第１結合対象フィラー粒子モデル１５ａに関連付けられた結合対象ポリマー粒子モデル１６及び結合対象カップリング粒子モデル１９のうち、粒子間距離Ｌｐが小さい順に、上限値Ｎｓ（本例では、３個）分の結合対象ポリマー粒子モデル１６又は結合対象カップリング粒子モデル１９が選択される。本例では、第１結合対象ポリマー粒子モデル１６ａ、第２結合対象ポリマー粒子モデル１６ｂ、及び、第１結合対象カップリング粒子モデル１９ａが選択される。

図２８に示されるように、選択された第１結合対象ポリマー粒子モデル１６ａ、第２結合対象ポリマー粒子モデル１６ｂ及び第１結合対象カップリング粒子モデル１９ａは、第１結合対象フィラー粒子モデル１５ａに結合されるものとして特定される。なお、図２８において、結合対象フィラー粒子モデル１５に結合されるものとして特定された結合対象ポリマー粒子モデル１６及び結合対象カップリング粒子モデル１９は、実線で示されている。

一方、選択されなかった第３結合対象ポリマー粒子モデル１６ｃ及び第２結合対象カップリング粒子モデル１９ｂは、第１結合対象フィラー粒子モデル１５ａに結合されるものとして特定されない。

図２７に示したように、第２結合対象フィラー粒子モデル１５ｂに関連付けられている第４結合対象ポリマー粒子モデル１６ｄの数（本例では、１個）は、上限値Ｎｓ（本例では、３個）よりも小さい。このため、図２８に示されるように、第４結合対象ポリマー粒子モデル１６ｄは、第２結合対象フィラー粒子モデル１５ｂに結合されるものとして特定される。

このように、第２工程Ｓ７２では、結合対象フィラー粒子モデル１５に結合させる結合対象ポリマー粒子モデル１６の数と、結合対象フィラー粒子モデル１５に結合させる結合対象カップリング粒子モデル１９の数との和が、上限値Ｎｓ以下となるように、結合対象ポリマー粒子モデル１６又は結合対象カップリング粒子モデル１９が特定される。なお、結合対象ポリマー粒子モデル１６及び結合対象カップリング粒子モデル１９を特定する情報は、例えば、結合対象ポリマー粒子モデル１６、結合対象カップリング粒子モデル１９及び結合対象フィラー粒子モデル１５の識別番号や、仮想空間９での各座標値が含まれる。これらの結合対象ポリマー粒子モデル１６及び結合対象カップリング粒子モデル１９を特定する情報は、コンピュータ１に記憶される。

次に、本実施形態の特定工程Ｓ６２は、第２工程Ｓ７２で特定されなかった各結合対象ポリマー粒子モデル１６又は各結合対象カップリング粒子モデル１９を、新たな結合対象フィラー粒子モデル１５に関連付けられる（第３工程Ｓ７３）。図２９は、特定されなかった結合対象ポリマー粒子モデル１６又は結合対象カップリング粒子モデル１９の関連付けを説明する概念図である。

第３工程Ｓ７３では、先ず、仮想空間９（図２２に示す）において、第２工程Ｓ７２で特定されなかった各結合対象ポリマー粒子モデル１６及び結合対象カップリング粒子モデル１９から、第１工程Ｓ７１で関連付けられた結合対象フィラー粒子モデル１５の次に粒子間距離Ｌｐが小さい結合対象フィラー粒子モデル１５が探索される。そして、各結合対象ポリマー粒子モデル１６及び結合対象カップリング粒子モデル１９が、探索された結合対象フィラー粒子モデル１５に関連付けられる。関連付けに関する情報は、コンピュータ１に記憶される。

図２８に示した例において、第２工程Ｓ７２で特定されなかった結合対象ポリマー粒子モデル１６及び結合対象カップリング粒子モデル１９は、第３結合対象ポリマー粒子モデル１６ｃ及び第２結合対象カップリング粒子モデル１９ｂである。第３工程Ｓ７３では、図２９に示されるように、第３結合対象ポリマー粒子モデル１６ｃ及び第２結合対象カップリング粒子モデル１９ｂから、第１結合対象フィラー粒子モデル１５ａの次に粒子間距離Ｌｐが小さい結合対象フィラー粒子モデル１５（本例では、第２結合対象フィラー粒子モデル１５ｂ）が探索される。そして、第３結合対象ポリマー粒子モデル１６ｃ及び第２結合対象カップリング粒子モデル１９が、第２結合対象フィラー粒子モデル１５ｂに関連付けられる。

次に、本実施形態の特定工程Ｓ６２では、各結合対象フィラー粒子モデル１５に結合させる結合対象ポリマー粒子モデル１６又は結合対象カップリング粒子モデル１９が、第３工程Ｓ７３で関連付けられた結合対象ポリマー粒子モデル１６又は結合対象カップリング粒子モデル１９から特定される（第４工程Ｓ７４）。第４工程Ｓ７４では、第２工程Ｓ７２と同様に、結合対象ポリマー粒子モデル１６の結合数Ｓａと、結合対象カップリング粒子モデル１９の結合数Ｓｂとの和（Ｓａ＋Ｓｂ）が上限値Ｎｓ以下を満たすように、関連付けられた結合対象ポリマー粒子モデル１６又は結合対象カップリング粒子モデル１９から粒子間距離Ｌｐが小さい順に特定される。

図３０は、第３工程Ｓ７３で関連付けられた結合対象ポリマー粒子モデル１６及び結合対象カップリング粒子モデル１９を特定する工程を説明する概念図である。図２９に示されるように、第２結合対象フィラー粒子モデル１５ｂには、第３結合対象ポリマー粒子モデル１６ｃ及び第２結合対象カップリング粒子モデル１９ｂが関連付けられている。第２結合対象フィラー粒子モデル１５ｂは、１個の結合対象ポリマー粒子モデル１６（第４結合対象ポリマー粒子モデル１６ｄ）が既に特定されている。このため、第２結合対象フィラー粒子モデル１５ｂでは、第３結合対象ポリマー粒子モデル１６ｃ及び第２結合対象カップリング粒子モデル１９ｂを含めた結合数の和（Ｓａ＋Ｓｂ）が、上限値Ｎｓ（本例では、３個）以下を満たす。このため、図３０に示されるように、第３結合対象ポリマー粒子モデル１６ｃ及び第２結合対象カップリング粒子モデル１９ｂは、第２結合対象フィラー粒子モデル１５ｂに結合されるものとして特定される。

なお、第２結合対象フィラー粒子モデル１５ｂに、新たに３個以上の結合対象ポリマー粒子モデル１６又は結合対象カップリング粒子モデル１９が関連付けられる場合は、結合数の和（Ｓａ＋Ｓｂ）が上限値Ｎｓ（本例では、３個）以下を満たすように、関連付けられた結合対象ポリマー粒子モデル１６及び結合対象カップリング粒子モデル１９から粒子間距離Ｌｐが小さい順に特定される。従って、粒子間距離Ｌｐが３番目以降に小さい結合対象ポリマー粒子モデル１６又は結合対象カップリング粒子モデル１９は、第２結合対象フィラー粒子モデル１５ｂに結合されるものとして特定されない。

このように、第４工程Ｓ７４では、結合対象フィラー粒子モデル１５に結合させる結合対象ポリマー粒子モデル１６、及び、結合対象カップリング粒子モデル１９を、結合数の和（Ｓａ＋Ｓｂ）が上限値Ｎｓ以下を満たすように、第３工程Ｓ７３で関連付けられた結合対象ポリマー粒子モデル１６及び結合対象カップリング粒子モデル１９から特定することができる。結合対象ポリマー粒子モデル１６及び結合対象カップリング粒子モデル１９を特定する情報は、コンピュータ１に記憶される。

次に、本実施形態の特定工程Ｓ６２では、全ての結合対象ポリマー粒子モデル１６又は結合対象カップリング粒子モデル１９が、結合対象フィラー粒子モデル１５に特定されたか否かが判断される（工程Ｓ７５）。工程Ｓ７５では、全ての結合対象ポリマー粒子モデル１６又は結合対象カップリング粒子モデル１９が、結合対象フィラー粒子モデル１５に特定されたと判断された場合、次の工程Ｓ６３が実施される。一方、各結合対象フィラー粒子モデル１５に特定されていない結合対象ポリマー粒子モデル１６又は結合対象カップリング粒子モデル１９が存在すると判断された場合は、特定されていない結合対象ポリマー粒子モデル１６又は結合対象カップリング粒子モデル１９を対象に、第３工程Ｓ７３及び第４工程Ｓ７４が実施される。これにより、特定工程Ｓ６２では、全ての結合対象ポリマー粒子モデル１６又は結合対象カップリング粒子モデル１９を、結合対象フィラー粒子モデル１５に特定することができる。なお、本実施形態では、結合対象ポリマー粒子モデル１６又は結合対象カップリング粒子モデル１９が特定されていない結合対象フィラー粒子モデル１５の存在が許容されている。

次に、結合工程Ｓ４６では、特定された結合対象ポリマー粒子モデル１６又は結合対象カップリング粒子モデル１９が、結合対象フィラー粒子モデル１５に結合される（工程Ｓ６３）。図３１は、特定された結合対象ポリマー粒子モデル１６又は結合対象カップリング粒子モデル１９を、結合対象フィラー粒子モデル１５に結合する工程を説明する概念図である。工程Ｓ６３では、特定された結合対象ポリマー粒子モデル１６と、結合対象フィラー粒子モデル１５とを連結鎖２２で結合させている。これにより、フィラーモデル３及びポリマーモデル５は、フィラーとポリマーとの化学結合を再現することができる。さらに、工程Ｓ６３では、特定された結合対象カップリング粒子モデル１９と、結合対象フィラー粒子モデル１５とを、連結鎖２２で結合させている。これにより、フィラーモデル３及びカップリングモデル１７は、フィラーとカップリング剤との化学結合を再現することができる。

次に、シミュレーション工程Ｓ４では、フィラーモデル３、ポリマーモデル５及びカップリングモデル１７を対象に分子動力学計算が実施される（工程Ｓ５５）。工程Ｓ５５は、前実施形態の工程Ｓ４７と同様に、ニュートンの運動方程式が適用される。そして、各時刻におけるフィラーモデル３の各フィラー粒子モデル４ｃ、４ｓ（図３に示す）、ポリマーモデル５の各ポリマー粒子モデル６ａ、６ｂ（図４に示す）、及び、カップリングモデル１７のカップリング粒子モデル１８の動きが追跡される。これにより、図１７に示されるように、工程Ｓ５５では、フィラーモデル３、ポリマーモデル５、及び、カップリングモデル１７を、仮想空間９内でそれぞれ分散させることができる。

本実施形態のシミュレーション工程Ｓ４では、結合工程Ｓ４６において、各結合対象フィラー粒子モデル１５に結合させる結合対象ポリマー粒子モデル１６の数Ｓａと、各結合対象フィラー粒子モデル１５に結合させる結合対象カップリング粒子モデル１９の数Ｓｂとの和（Ｓａ＋Ｓｂ）は、予め定められた上限値Ｎｓ以下に設定されている。このため、フィラー粒子モデル４に結合されるポリマー粒子モデル６及びカップリング粒子モデル１８の数に大きな偏りが生じるのを防ぐことができる。これにより、本実施形態のシミュレーション工程Ｓ４では、分子動力学計算において、結合対象フィラー粒子モデル１５に作用する結合対象ポリマー粒子モデル１６及び結合対象カップリング粒子モデル１９の力の総和が大きくなるのを防ぐことができる。従って、本実施形態のシミュレーション方法は、計算落ちが発生するのを防ぐことができ、分子動力学計算を安定して実行することができる。

また、本実施形態の上限値Ｎｓは、各結合対象フィラー粒子モデル１５に、結合対象ポリマー粒子モデル１６又は結合対象カップリング粒子モデル１９が均等に結合された場合の個数に基づいて、上記式（３）の関係を満たすように設定される。このため、結合対象フィラー粒子モデル１５に、結合対象ポリマー粒子モデル１６又は結合対象カップリング粒子モデル１９を均一に（満遍なく）結合させることができるため、結合対象フィラー粒子モデル１５に作用する結合対象ポリマー粒子モデル１６の力の総和が大きくなるのを確実に防ぐことができる。

さらに、本実施形態の結合対象フィラー粒子モデル１５は、多面体１３の頂点１３ａ（図３に示す）に配置された頂点粒子モデル７に限定されるため、隣接する結合対象フィラー粒子モデル１５、１５にそれぞれ結合される結合対象ポリマー粒子モデル１６又は結合対象カップリング粒子モデル１９を離間させることができる。これにより、結合対象ポリマー粒子モデル１６又は結合対象カップリング粒子モデル１９に作用する相互作用が大きくなるのを防ぐことができるため、分子動力学計算を安定して実行することができる。

次に、シミュレーション工程Ｓ４では、コンピュータ１が、フィラーモデル３、ポリマーモデル５及びカップリングモデル１７の十分に分散したか否かを判断する（工程Ｓ５６）。本実施形態の工程Ｓ５６では、フィラーモデル３、ポリマーモデル５及びカップリングモデル１７が十分に分散したと判断された場合（工程Ｓ５６で「Ｙ」）、次の評価工程Ｓ５（図２に示す）が実施される。一方、フィラーモデル３、ポリマーモデル５及びカップリングモデル１７が十分に分散していないと判断された場合（工程Ｓ５６で「Ｎ」）は、単位ステップを進めて（工程Ｓ４９）、工程Ｓ５５が再度実施される。従って、本実施形態のシミュレーション工程Ｓ４では、フィラーモデル３及びポリマーモデル５を効果的に分散させることができ、構造緩和された高分子材料モデル１１を設定することができる。なお、フィラーモデル３、ポリマーモデル５及びカップリングモデル１７を効果的に分散させるために、工程Ｓ５５において、１ステップあたりの分子動力学計算の時間幅を０．０５［τ］したときのステップ数は、前実施形態の工程Ｓ４７と同一範囲が望ましい。

次に、本実施形態のシミュレーション方法では、図２に示されるように、フィラーモデル３の分散状態が良好か否か判断される（評価工程Ｓ５）。評価工程Ｓ５では、前実施形態の評価工程Ｓ５と同様に、フィラーモデル３の中心粒子モデル４ｃを対象とした動径分布関数が計算される。そして、評価工程Ｓ５では、動径分布関数の結果をもとに、コンピュータ１が予め設定された許容範囲内であるか否かが判断される。

評価工程Ｓ５では、フィラーモデル３の分散状態が良好であると判断された場合、高分子材料モデル１１に基づいて、高分子材料が製造される（工程Ｓ６）。一方、フィラーモデル３の分散状態が良好でないと判断された場合は、動径分布関数の結果に基づいて、例えば、フィラーモデル３及びポリマーモデル５の設定条件を変更して（工程Ｓ７）、工程Ｓ１〜工程Ｓ５が実施される。これにより、本実施形態のシミュレーション方法では、高分子材料モデル１１（図１７に示す）の構造を、効果的に緩和することができ、実際の高分子材料の構造に、確実に近似させることができる。

また、本実施形態のシミュレーション方法では、フィラーモデル３、ポリマーモデル５及びカップリングモデル１７の構造やポテンシャルを任意に設定することができるため、現実には存在しない未知の高分子材料をモデル化した高分子材料モデル１１を評価することができる。

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

［実施例Ａ］
図２、図９、図１０及び図１１に示す手順に従って、予め定められた上限値Ｓｕ以下で、各結合対象フィラー粒子モデルに結合対象ポリマー粒子モデルを結合させて、分子動力学計算が実施された（実施例１）。

比較のために、上限値Ｓｕを設定することなく、各結合対象フィラー粒子モデルに、結合対象ポリマー粒子モデルをランダムに結合させて、分子動力学計算が実施された（比較例１）。

実施例１及び比較例１の各シミュレーション方法がそれぞれ３０回実施され、分子動力学計算に計算落ちが発生するか否かが確認された。各パラメータは、明細書中の記載通りであり、共通仕様は次のとおりである。
フィラーモデル：
合計個数：１００個
各フィラーモデルの結合対象フィラー粒子モデルの個数：８個
ポリマーモデル：
合計個数：２０００個
各ポリマーモデルの結合対象ポリマー粒子モデルの個数：１個
結合対象ポリマー粒子モデルの総数Ｎａ：２０００
結合対象フィラー粒子モデルの総数Ｎｂ：８００
Ｎａ／Ｎｂ：２．５
上限値Ｓｕ：２

テストの結果、実施例１では、全ての計算において計算落ちが発生しなかった。一方、比較例１では、分子動力学計算に計算落ちが１５回発生した。従って、実施例のシミュレーション方法では、分子動力学計算を安定して実行できることが確認できた。

［実施例Ｂ］
図１８、図２３、図２５及び図２６に示す手順に従って、予め定められた上限値Ｎｓ以下で、各結合対象フィラー粒子モデルに、結合対象ポリマー粒子モデル又は結合対象カップリング粒子モデルを結合させて、分子動力学計算が実施された（実施例２〜実施例４）。

実施例２〜実施例４のカップリング粒子モデルの個数、及び、ポテンシャルＵ１１〜Ｕ１４の強度ａ_ijは、次のとおりである。実施例３は、実施例２に比べて、カップリング粒子モデルの個数が小に設定されている。実施例４は、カップリング粒子モデルのフィラー粒子モデルへの親和性が、ポリマー粒子モデルへの親和性よりも高く設定されている。
実施例２：
カップリング粒子モデルの個数：５個
ポテンシャルＵ１１〜Ｕ１４の強度ａ_ij：５０
実施例３：
カップリング粒子モデルの個数：１個
ポテンシャルＵ１１〜Ｕ１４の強度ａ_ij：５０
実施例４：
カップリング粒子モデルの個数：５個
ポテンシャルＵ１１の強度ａ_ij：１０
ポテンシャルＵ１２の強度ａ_ij：１０
ポテンシャルＵ１３の強度ａ_ij：１００
ポテンシャルＵ１４の強度ａ_ij：１００

比較のために、上限値Ｎｓを設定することなく、各結合対象フィラー粒子モデルに、結合対象ポリマー粒子モデル又は結合対象カップリング粒子モデルをランダムに結合させて、分子動力学計算が実施された（比較例２）。

実施例２〜実施例４及び比較例２の各シミュレーション方法がそれぞれ３０回実施され、分子動力学計算に計算落ちが発生するか否かが確認された。各パラメータは、上記したカップリング粒子モデルの仕様を除き、明細書中の記載通りであり、共通仕様は次のとおりである。図３２には、実施例２及び実施例３の動径分布関数が示される。また、図３３には、実施例２及び実施例４の動径分布関数が示される。
結合対象ポリマー粒子モデルの総数Ｎａ：１０００
結合対象フィラー粒子モデルの総数Ｎｂ：８００
結合対象カップリング粒子モデルの総数Ｎｃ：２５００
（Ｎａ＋Ｎｃ）／Ｎｂ＋１：５．３７５
上限値Ｎｓ：６

テストの結果、実施例２〜４では、全ての計算において計算落ちが発生しなかった。一方、比較例２では、分子動力学計算に計算落ちが５０回発生した。従って、実施例のシミュレーション方法では、分子動力学計算を安定して実行できることが確認できた。

図３２に示されるように、実施例２の動径分布関数は、実施例３の動径分布関数よりも小さくなっている。これは、実施例２のフィラーモデルが、実施例３のフィラーモデルよりも良好に分散していることを示している。一般に、カップリング剤の鎖長が大きいほど、フィラーの分散性が高くなることが知られている。従って、本発明のシミュレーション方法では、カップリング剤の鎖長の差異によって生じるフィラーの分散性の影響を、正確に評価できることが確認できた。

図３３に示されるように、実施例２の動径分布関数は、実施例４の動径分布関数よりも小さくなっている。これは、実施例２のフィラーモデルが、実施例４のフィラーモデルよりも良好に分散していることを示している。一般に、フィラーに対する親和性とポリマーに対する親和性が同一のカップリング剤は、フィラーの分散性が高くなることが知られている。従って、本発明のシミュレーション方法では、カップリング剤の相互作用の差異によって生じるフィラーの分散性の影響を評価することができることが確認できた。

１ コンピュータ
３ フィラーモデル
４ フィラー粒子モデル
５ ポリマーモデル
６ ポリマー粒子モデル
９ 仮想空間
１５ 結合対象フィラー粒子モデル
１６ 結合対象ポリマー粒子モデル

Claims (9)

1. 高分子材料中に配合されたフィラーの分散性を、コンピュータを用いて評価するためのシミュレーション方法であって、
   前記コンピュータに、前記フィラーを、複数のフィラー粒子モデルでモデル化したフィラーモデルを設定する工程、
   前記コンピュータに、前記高分子材料を、複数のポリマー粒子モデルでモデル化したポリマーモデルを設定する工程、及び
   前記コンピュータが、予め定められた仮想空間に配置された前記フィラーモデルと前記ポリマーモデルとを用いて分子動力学計算を行うシミュレーション工程を含み、
   前記シミュレーション工程は、前記フィラー粒子モデルから予め定められた結合対象フィラー粒子モデルに、前記ポリマー粒子モデルから予め定められた結合対象ポリマー粒子モデルを結合させる結合工程を含み、
   前記各結合対象フィラー粒子モデルに結合させる前記結合対象ポリマー粒子モデルの数は、予め定められた上限値Ｓｕ以下であることを特徴とする高分子材料のシミュレーション方法。
2. 前記上限値Ｓｕは、下記の関係を満たす請求項１記載の高分子材料のシミュレーション方法。
   Ｎａ／Ｎｂ−１≦Ｓｕ≦Ｎａ／Ｎｂ＋１
   ここで、各パラメータは、次のとおりである。
   Ｎａ：結合対象ポリマー粒子モデルの総数
   Ｎｂ：結合対象フィラー粒子モデルの総数
3. 前記結合工程は、前記上限値Ｓｕを設定する工程、
   前記各結合対象フィラー粒子モデルに結合させる前記結合対象ポリマー粒子モデルを特定する特定工程、及び
   特定された前記結合対象ポリマー粒子モデルを、前記結合対象フィラー粒子モデルに結合させる工程を含む請求項１又は２に記載の高分子材料のシミュレーション方法。
4. 前記特定工程は、前記各結合対象ポリマー粒子モデルを、夫々粒子間距離が最も小さい前記結合対象フィラー粒子モデルに関連付ける第１工程、
   前記各結合対象フィラー粒子モデルに結合させる前記結合対象ポリマー粒子モデルを、前記上限値Ｓｕ以下を満たすように、関連付けられた結合対象ポリマー粒子モデル群から粒子間距離が小さい順に特定する第２工程、
   前記第２工程で特定されなかった前記各結合対象ポリマー粒子モデルを、前記関連付けられた結合対象フィラー粒子モデルの次に粒子間距離が小さい前記結合対象フィラー粒子モデルに関連付ける第３工程、及び
   前記各結合対象フィラー粒子モデルに結合させる前記結合対象ポリマー粒子モデルを、前記上限値Ｓｕ以下を満たすように、前記第３工程で関連付けられた結合対象ポリマー粒子モデル群から粒子間距離が小さい順に特定する第４工程を含み、
   全ての前記結合対象ポリマー粒子モデルが前記結合対象フィラー粒子モデルに特定されるまで、前記第３工程及び前記第４工程が繰り返し実施される請求項３記載の高分子材料のシミュレーション方法。
5. フィラーと、前記フィラーをポリマーに結合させるためのカップリング剤とが配合された高分子材料について、前記フィラーの分散性を、コンピュータを用いて評価するためのシミュレーション方法であって、
   前記コンピュータに、前記フィラーを、複数のフィラー粒子モデルでモデル化したフィラーモデルを設定する工程、
   前記コンピュータに、前記高分子材料を、複数のポリマー粒子モデルでモデル化したポリマーモデルを設定する工程、
   前記コンピュータに、前記カップリング剤を、少なくとも一つのカップリング粒子モデルでモデル化したカップリングモデルを設定する工程、並びに
   前記コンピュータが、予め定められた仮想空間に配置された前記フィラーモデル、前記ポリマーモデル、及び、前記カップリングモデルを用いて分子動力学計算を行うシミュレーション工程を含み、
   前記シミュレーション工程は、前記フィラー粒子モデルから予め定められた結合対象フィラー粒子モデルに、前記ポリマー粒子モデルから予め定められた結合対象ポリマー粒子モデル、又は、前記カップリング粒子モデルから予め定められた結合対象カップリング粒子モデルを結合させる結合工程を含み、
   前記各結合対象フィラー粒子モデルに結合させる前記結合対象ポリマー粒子モデルの数と、前記各結合対象フィラー粒子モデルに結合させる前記結合対象カップリング粒子モデルの数との和は、予め定められた上限値Ｎｓ以下であることを特徴とする高分子材料のシミュレーション方法。
6. 前記上限値Ｎｓは、下記の関係を満たす請求項５記載の高分子材料のシミュレーション方法。
   Ｎｓ＝［（Ｎａ＋Ｎｃ）／Ｎｂ＋１］
   ここで、各パラメータは、次のとおりである。
   Ｎａ：結合対象ポリマー粒子モデルの総数
   Ｎｂ：結合対象フィラー粒子モデルの総数
   Ｎｃ：結合対象カップリング粒子モデルの総数
7. 前記結合工程は、前記上限値Ｎｓを設定する工程、
   前記各結合対象フィラー粒子モデルに結合させる前記結合対象ポリマー粒子モデル、又は、前記各結合対象フィラー粒子モデルに結合させる前記結合対象カップリング粒子モデルを特定する特定工程、並びに
   特定された前記結合対象ポリマー粒子モデル又は前記結合対象カップリング粒子モデルを、前記結合対象フィラー粒子モデルに結合させる工程を含む請求項５又は６に記載の高分子材料のシミュレーション方法。
8. 前記特定工程は、前記各結合対象ポリマー粒子モデル及び前記各結合対象カップリング粒子モデルを、夫々粒子間距離が最も小さい前記結合対象フィラー粒子モデルに関連付ける第１工程、
   前記各結合対象フィラー粒子モデルに結合させる前記結合対象ポリマー粒子モデル又は前記各結合対象カップリング粒子モデルを、前記上限値Ｎｓ以下を満たすように、関連付けられた結合対象ポリマー粒子モデル又は結合対象カップリング粒子モデルから粒子間距離が小さい順に特定する第２工程、
   前記第２工程で特定されなかった前記各結合対象ポリマー粒子モデル又は前記各結合対象カップリング粒子モデルを、前記関連付けられた結合対象フィラー粒子モデルの次に粒子間距離が小さい前記結合対象フィラー粒子モデルに関連付ける第３工程、並びに
   前記各結合対象フィラー粒子モデルに結合させる前記結合対象ポリマー粒子モデル又は前記結合対象カップリング粒子モデルを、前記上限値Ｎｓ以下を満たすように、前記第３工程で関連付けられた結合対象ポリマー粒子モデル又は結合対象カップリング粒子モデルから粒子間距離が小さい順に特定する第４工程を含み、
   全ての前記結合対象ポリマー粒子モデル又は前記結合対象カップリング粒子モデルが前記結合対象フィラー粒子モデルに特定されるまで、前記第３工程及び前記第４工程が繰り返し実施される請求項７記載の高分子材料のシミュレーション方法。
9. 前記フィラー粒子モデルは、予め定められた多面体の頂点に少なくとも配置され、
   前記結合対象フィラー粒子モデルは、前記頂点に配置されたフィラー粒子モデルである請求項１乃至８のいずれかに記載の高分子材料のシミュレーション方法。