JP6294613B2 - 高分子材料のシミュレーション方法 - 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内）に配置して、フィラーモデル３の凝集塊が形成される。これにより、仮想空間８には、解析対象のゴムポリマー内に形成されがちなフィラーの凝集塊を再現でき、現実に近いシミュレーションを行うことができる。

次に、分子動力学（ Molecular Dynamics : ＭＤ）計算が実施される（緩和計算工程Ｓ１４）。図１０は、本実施形態の緩和計算工程Ｓ１４の処理手順の一例を示すフローチャートである。

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

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

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

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

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

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

次に、フィラー粒子４とポリマー粒子６とが連結鎖で結合される（結合工程Ｓ１４６）。図１１（ａ）、（ｂ）は、結合工程Ｓ１４６を説明する概念図である。この結合工程Ｓ１４６では、予め定められた粒子間距離Ｌ２以下に接近したフィラー粒子４とポリマー粒子６とを、連結鎖１０で結合させている。これにより、緩和計算工程Ｓ１４では、フィラーとポリマーとの化学結合を再現することができる。

本実施形態の結合工程Ｓ１４６では、各ポリマー粒子６ａ、６ｂが、フィラーモデル３の表面粒子４ｓのみに結合させている。上述したように、表面粒子４ｓは、中心粒子４ｃに相対位置を保って連結されるため、フィラーモデル３に対するポリマーモデル５の位置を一意に定めることができる。従って、フィラーモデル３は、ポリマーモデル５が滑り運動するのを防ぐことができ、シミュレーション精度を確実に向上しうる。また、粒子間距離Ｌ２は、表面粒子４ｓが関連するポテンシャルの各カットオフ距離ｒ_c（図７に示す）の１０％〜２００％が望ましい。

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

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

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

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

次に、フィラーモデル３の分散状態が良好か否か判断される（評価工程Ｓ１５）。この評価工程Ｓ１５では、先ず、フィラーモデル３の中心粒子４ｃを対象とした動径分布関数ｇ（ｒ）が計算される。動径分布関数ｇ（ｒ）とは、ある中心粒子４ｃ（図３に示す）から距離ｒ離れた位置において、他の中心粒子４ｃが存在する確率密度を表す関数であり、下記の式（２）で定義される。なお、距離ｒは、中心粒子４ｃの中心間の距離として定義される。

上記式（２）において、n(r)は、距離rと、距離r＋Δrとの間に存在する粒子の数を示している。〈n(r)〉は、距離rと、距離r＋Δrとの間に存在する粒子の数n(r)について、粒子平均および時刻平均をとったもの（アンサンブル平均したもの）を示している。ρは計算系全体の数密度を示している。

図１３は、本実施形態の中心粒子４ｃ、４ｃ間の動径分布を示すグラフである。この動径分布では、動径分布関数ｇ（r）がピークとなるrp付近から、距離rが最大値となるrmにかけて、万遍なく分布されていることが確認できる。これは、フィラーモデル３が、仮想空間８の中で、万遍なく分散されていることを示している。このように、本実施形態では、中心粒子４ｃのみを対象に動径分布関数ｇ（r）を求めることより、フィラーモデル３の分散状態を確認することができるため、計算コストの増大を抑制することができる。

なお、動径分布関数の距離rの範囲は、その最小値rsが０に設定されるとともに、最大値rmが仮想空間８の一辺の長さＬ１（図９に示す）の半分に設定されるのが望ましい。これにより、動径分布関数ｇ（r）を求める距離rの範囲が、仮想空間８内に限定されるため、計算コストの増大を抑制しつつ、正確な評価を行うことができる。なお、仮想空間８の一辺の長さＬ１がそれぞれ同一でない場合、前記最大値rmは、全ての辺の最小の長さの半分に設定されるのが望ましい。

また、動径分布関数ｇ（r）の取得間隔rdは、好ましくは、前記最大値rmを５で除した距離以下、さらに好ましくは、前記最大値rmを１０で除した距離以下に設定されるのが望ましい。これにより、動径分布関数ｇ（r）の精度を高めることかでき、正確な評価を行うことができる。なお、取得間隔rdが小さすぎても、計算コストが増大するおそれがある。このような観点より、取得間隔rdは、最大値rmを１０００で除した距離以上に設定されるのが望ましい。なお、図１３に示される動径分布において、取得間隔rdは、最大値rmを１５５で除した距離である。

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

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

図１４は、本実施形態のモデル化工程Ｓ２の処理手順の一例を示すフローチャートである。本実施形態のモデル化工程Ｓ２では、高分子材料モデル１１（図１２に示す）に基づいて、有限個の要素でモデル化した二次元の有限要素モデルが作成される。

本実施形態のモデル化工程Ｓ２では、先ず、高分子材料モデル１１（図１２に示す）の二次元断面が取得される（工程Ｓ２１）。図１５は、高分子材料モデル１１の二次元断面を示す部分断面図である。二次元断面１２は、コンピュータ１に、三次元の高分子材料モデル１１の断面の位置を指定することによって、任意の断面形状を取得することができる。このような二次元断面１２では、高分子材料モデル１１が所定の位置で切断されたフィラーモデル３及びポリマーモデル５が含まれている。

また、本実施形態の工程Ｓ２１では、フィラーモデル３の輪郭３ｓが、フィラーモデル３の重心点を中心とする単一球Ｂの表面Ｂｓ（図１２に示す）で設定されている。これにより、フィラーモデル３は、中心粒子４ｃと表面粒子４ｓとを含む多面体状（図１２に示す）から、現実のフィラーの輪郭に近似させることができる。このような二次元断面１２は、構造緩和された高分子材料モデル１１に基づいて設定されるため、例えば、現実の高分子材料の電子線透過画像からトモグラフィー法によって得られた二次元断面に近似させることができる。しかも、工程Ｓ２１では、電子線透過画像から作成される場合とは異なり、高分子材料部分と充填剤部分とを識別する画像処理が必要ないため、計算時間を短縮することができる。この二次元断面１２でのフィラーモデル３の座標値、及び、ポリマーモデル５の座標値が、コンピュータ１に記憶される。

次に、コンピュータ１が、二次元断面１２のフィラーモデル３及びポリマーモデル５をモデル化して、第１の有限要素モデル２Ａ（図１８に示す）を作成する（第１工程Ｓ２２）。図１６は、本実施形態の第１工程Ｓ２２の処理手順の一例を示すフローチャートである。図１７（ａ）は、二次元空間を示す平面図である。図１７（ｂ）は、（ａ）を拡大して示す平面図である。

図１７に示されるように、第１工程Ｓ２２では、先ず、複数の要素１４で格子状に区分された二次元空間１５に、図１５に示したフィラーモデル３の輪郭３ｓの座標があてはめられる（工程Ｓ２２１）。

本実施形態の二次元空間１５は、高分子材料モデル１１の二次元断面１２（図１５に示す）と同一形状に設定されている。これにより、二次元断面１２を二次元空間１５に重ね合わせることにより、フィラーモデル３の輪郭３ｓを、二次元空間１５に容易に定義することができる。なお、二次元空間１５の形状は、高分子材料モデル１１の二次元断面１２の形状と異ならせてもよい。このような場合は、隣接するフィラーモデル３、３の間隔を維持したまま、任意のフィラーモデル３の輪郭３ｓが、二次元空間１５にあてはめられる。これにより、フィラーモデル３の分散状態を維持したまま、フィラーモデル３の輪郭３ｓを、二次元空間１５に定義することができる。二次元空間１５は、例えば一辺の長さＬ３（図１８に示す）が５０nm〜１５０nm程度の正方形に設定されるのが望ましい。

本実施形態の要素１４は、図１７（ｂ）に示されるように、４つの節点１７と、各節点１７、１７間を直線でつなぐ４本の辺１８とで形成される。また、要素１４は、その一辺の長さＬ４が、０．１nm〜３．０nm程度の正方形に形成される。なお、要素１４は、本実施形態のような正方形に限定されるわけではなく、例えば長方形ものでもよい。

次に、各要素１４の中心１４ｃが、フィラーモデル３の輪郭３ｓよりもフィラーモデル３の中心３ｃ側に配されるか否かが判断される（工程Ｓ２２２）。なお、フィラーモデル３が非円形である場合、フィラーモデル３の中心３ｃは、図心として定義される。

この工程Ｓ２２２では、各要素１４の中心１４ｃが、フィラーモデル３の輪郭３ｓよりもフィラーモデル３の中心３ｃ側に配されていると判断された場合（工程Ｓ２２２で「Ｙ」）、該要素１４が、フィラーモデル３をなすフィラーモデル部２１の要素１４として定義される（工程Ｓ２２３）。一方、各要素１４の中心１４ｃが、フィラーモデル３の輪郭３ｓよりもフィラーモデル３の中心３ｃ側に配されていないと判断された場合（工程Ｓ２２２で「Ｎ」）は、該要素１４が、ポリマーモデル５をなすポリマーモデル部２２の要素１４として定義される（工程Ｓ２２４）。そして、フィラーモデル部２１及びポリマーモデル部２２の各要素１４の座標が、コンピュータ１に記憶される（工程Ｓ２２５）。

次に、全ての要素１４がフィラーモデル部２１又はポリマーモデル部２２としてモデル化されたか否かが判断される（工程Ｓ２２６）。この工程Ｓ２２６では、全ての要素１４が、フィラーモデル部２１又はポリマーモデル部２２として定義されたと判断された場合（工程Ｓ２２６で「Ｙ」）、次の第２工程Ｓ２３が実施される。一方、全ての要素１４が、フィラーモデル部２１又はポリマーモデル部２２として定義されていないと判断される場合（工程Ｓ２２６で「Ｎ」）は、工程Ｓ２２２〜工程Ｓ２２６が再度実施される。

図１８は、第１の有限要素モデルを示す平面図である。上記のような工程Ｓ２２１〜工程Ｓ２２６を経ることにより、二次元空間１５には、図１８に示されるように、フィラーモデル部２１とポリマーモデル部２２とを含む第１の有限要素モデル２６Ａが定義される。この例では、理解しやすいように、フィラーモデル部２１が着色されて表示されている。このような第１の有限要素モデル２６Ａは、コンピュータ１に記憶される。

次に、フィラーモデル部２１とポリマーモデル部２２との境界の節点１７ａを含む要素１４を変形させた第２の有限要素モデルが作成される（第２工程Ｓ２３）。図１９は、本実施形態の第２工程Ｓ２３の処理手順の一例を示すフローチャートである。図２０（ａ）は、要素を拡大して示す平面図である。図２０（ｂ）は、変形後の要素を示す平面図である。

本実施形態の第２工程Ｓ２３では、先ず、図２０（ａ）に示されるように、コンピュータ１が、フィラーモデル３の中心３ｃ（図１８に示す）と、境界の節点１７ａとを通る直線２７が定義される（工程Ｓ２３１）。次に、直線２７とフィラーモデル３の輪郭３ｓとの交点２８が求められる（工程Ｓ２３２）。これらの直線２７及び交点２８は、コンピュータ１に記憶される。

次に、コンピュータ１が、境界の節点１７ａを、交点２８に移動させる（工程Ｓ２３３）。本実施形態では、境界の節点１７ａを、直線２７上で移動させている。この境界の節点１７ａの移動に伴い、図２０（ｂ）に示されるように、境界の節点１７ａを含む要素１４のみを変形させることができる。そして、境界の節点１７ａを含む要素１４の各座標が、コンピュータ１に記憶される（工程Ｓ２３４）。

次に、全ての境界の節点１７ａが、直線２７の交点２８に移動されたか否かが判断される（工程Ｓ２３５）。この工程Ｓ２３５では、全ての境界の節点１７ａが、直線２７の交点２８に移動されたと判断された場合（工程Ｓ２２５で「Ｙ」）、次の変形工程Ｓ３が実施される。一方、全ての境界の節点１７ａが、直線２７の交点２８に移動されていないと判断された場合（工程Ｓ２２５で「Ｎ」）は、新たな境界の節点１７ａが選択され、工程Ｓ２３１〜工程Ｓ２３５が再度実施される。これにより、第２工程Ｓ２３では、全ての境界の節点１７ａ（図２０（ｂ）に示す）を、フィラーモデル部２１の表面の凹凸２９（図１８に示す）が小さくなる向きに移動させることができる。

図２１は、第２の有限要素モデルを示す平面図である。第２工程Ｓ２３では、上記のような工程Ｓ２３１〜工程Ｓ２３４を経ることにより、境界の節点１７ａを含む要素１４を変形させた第２の有限要素モデル２６Ｂ（図２１に示す）を作成することができる。

第２の有限要素モデル２６Ｂでは、図２０（ｂ）に示されるように、境界の節点１７ａ間をのびる辺１８が滑らかに連なり、フィラーモデル部２１の表面の凹凸２９（図１８に示す）が小に形成される。一方、境界の節点１７ａを含まない要素１４は、変形されることなく、第１の有限要素モデル２６Ａ（図１８に示す）のときと同じ形状（本実施形態では、正方形）が維持されるため、処理時間を大幅に短縮しうる。

また、第２の有限要素モデル２６Ｂのフィラーモデル部２１及びポリマーモデル部２２の各要素１４には、シミュレーションによる数値解析に必要な情報が定義される。この数値解析とは、例えば有限要素法等の数値解析法を意味する。また、解析に必要な情報としては、各要素１４を構成する節点１７の番号や座標値が少なくとも含まれる。

さらに、各要素１４には、各々が代表する部分の材料特性（物性値）などが定義される。即ち、フィラーモデル部２１及びポリマーモデル部２２の各要素１４には、高分子材料及びフィラーの物性に応じた材料定数が、それぞれ定義される。そして、これらの情報は、いずれもコンピュータ１に記憶される。

図２２（ａ）は、他の要素を拡大して示す平面図、（ｂ）は変形後の要素を示す平面図である。図２３は、本発明の他の実施形態の第２工程Ｓ２３の処理手順の一例を示すフローチャートである。図２２（ａ）に示されるように、第１の有限要素モデル２６Ａの複数の境界の節点１７ａ、１７ａが、同一の前記直線２７上に位置する場合、各境界の節点１７ａ、１７ａを、直線２７とフィラーモデル３の輪郭３ｓとの交点２８に移動させると、互いに重複する。このような重複は、その境界の節点１７ａを含む要素１４が、例えば四辺形から三角形へと潰れてしまい、それらの剛性を計算できないおそれがある。

このような重複を抑制するために、第２工程Ｓ２３では、図２２（ａ）及び図２３に示されるように、コンピュータ１が、境界の節点１７ａを交点２８に移動させる工程Ｓ２３３の前に、同一の直線２７上で複数の境界の節点１７ａが重複するか否かを確認する工程Ｓ２３６が実施されるのが望ましい。

工程Ｓ２３６では、同一の直線２７上で複数の境界の節点１７ａが重複すると判断された場合（工程Ｓ２３６で「Ｙ」）、重複する境界の節点１７ａ、１７ａと交点２８との距離Ｌ５、Ｌ５をそれぞれ計算して（工程Ｓ２３７）、該距離Ｌ５、Ｌ５に係数を乗じた修正距離で、境界の節点１７ａ、１７ａを直線２７に沿って交点２８側へ移動させる（工程Ｓ２３８）のが望ましい。本実施形態の係数は、０よりも大かつ１よりも小の数値が設定される。これにより、図２２（ｂ）に示されるように、境界の節点１７ａの重複が抑制され、シミュレーション精度が向上する。

また、工程Ｓ２３６では、同一の直線２７上で複数の境界の節点１７ａが重複しないと判断された場合（工程Ｓ２３６で「Ｎ」）、図１９に示した第２工程Ｓ２３のフローチャートと同様に、境界の節点１７ａを交点２８に移動させる工程Ｓ２３３が実施される。このように、第２工程Ｓ２３では、境界の節点１７ａの重複を防ぎつつ、フィラーモデル部２１の表面に形成される凹凸２９（図１８に示す）を小さくすることができる。

図２４には、図２１に示した有限要素モデル２６（第２の有限要素モデル２６Ｂ）を、一方向（例えば、Ｙ軸方向）に０％〜２０％伸長させたときの応力及び伸びの計算結果が実線で示されている（実施例１）。さらに、図２４には、高分子材料モデル１１の構造に基づいて製造された高分子材料（加硫ゴム）を、０％〜２０％伸長させたときの応力と伸びとの測定結果が破線で示されている（実験例）。

また、高分子材料にフィラーが規則的に分散されていると仮定し、二次元空間に任意の輪郭を設定して、実施例１と同様のモデル化工程Ｓ２を経て、有限要素モデルが作成された。そして、図２４には、この有限要素モデルを、０％〜２０％伸長させたときの応力と伸びとの計算結果が、一点鎖線で示されている（比較例）。なお、実施例１、比較例及び実験例の共通仕様は、次のとおりである。また、各パラメータは、明細書に記載しているとおりである。
シミュレーションソフト：有限要素解析アプリケーションソフト（LSTC社製のLS-DYNA）
実施例１及び比較例１の有限要素モデル：
要素の数：１６００個
二次元空間の一辺の長さＬ３：１００nm
要素の一辺の長さＬ４：２．５nm

これらの結果より、実施例１は、比較例１に比べて、実験例に近似していることが確認できる。これは、実施例１が、シミュレーション工程Ｓ１での構造緩和により、高分子材料モデル１１を、実際の高分子材料の分散状態に近似させることができたことによるものと考えられる。従って、本実施形態のシミュレーション方法（実施例１）は、比較例１に比べて、シミュレーション精度を向上しうる。

しかも、本実施形態の第２工程Ｓ２３では、フィラーモデル部２１の表面に形成される凹凸２９（図１８に示す）を小さくして、滑らかな表面形状に形成できる。従って、本実施形態では、フィラーモデル部２１に接触するポリマーモデル部２２の剛性が大に見積られることを抑制できる。これにより、実施例１は、実験例に効果的に近似させることができ、シミュレーション精度を向上しうる。

次に、変形工程Ｓ３では、コンピュータ１が、有限要素モデル２６（第２の有限要素モデル２６Ｂ）の変形計算を行い、有限要素モデル２６の物理量を計算する。本実施形態の変形工程Ｓ３では、上述のように、有限要素モデル２６を、一方向に（例えば、Ｙ軸方向に０％〜２０％）伸長させて、有限要素モデル２６の応力及び伸びが計算される。なお、有限要素モデル２６を変形させる方法については、上記のような方法に限定されるわけではない。例えば、有限要素モデル２６を１０％初期伸張させた後に、周期的な歪を±１％与えて変形させる方法や、有限要素モデル２６を圧縮又はせん断変形させる方法でもよい。

次に、有限要素モデル２６（第２の有限要素モデル２６Ｂ）の物理量が、予め設定された許容範囲内であるか否かが判断される（工程Ｓ４）。この工程Ｓ４では、第２の有限要素モデル２６Ｂの物理量が、許容範囲内であると判断された場合（工程Ｓ４で「Ｙ」）、高分子材料モデル１１に基づいて、高分子材料が製造される（工程Ｓ５）。一方、第２の有限要素モデル２６Ｂの物理量が、許容範囲内でないと判断された場合（工程Ｓ４で「Ｎ」）は、フィラーモデル３及びポリマーモデル５の設定等を変更して（工程Ｓ６）、工程Ｓ１〜工程Ｓ４が再度実施される。これにより、本実施形態のシミュレーション方法では、許容範囲の物理量を有する高分子材料を製造することができる。

このように、本実施形態のシミュレーション方法では、高分子材料モデル１１（図１２に示す）からモデル化された有限要素モデル２６（図２１に示す）を用いて、変形計算及びその評価をすることができるため、高分子材料を実際に作成する必要がない。従って、本実施形態のシミュレーション方法では、高分子材料を効率よく開発することができる。

しかも、本実施形態のシミュレーション方法では、高分子材料モデル１１として、現実には存在しない未知の高分子材料をモデル化することができる。従って、本発明の高分子材料のシミュレーション方法は、未知の高分子材料の開発に役立つ。さらに、本実施形態の高分子材料モデルは、分子動力学計算に基づいて、効果的に構造緩和されるため、図２４のグラフに示したように、シミュレーション精度が低下することもない。

本実施形態のモデル化工程Ｓ２では、高分子材料モデル１１（図１２に示す）から、二次元の有限要素モデル２６が作成されるものが例示されたが、これに限定されるわけではない。モデル化工程Ｓ２では、例えば、高分子材料モデル１１に基づいて、有限個の要素でモデル化した三次元の有限要素モデルが作成されてもよい。図２５は、この実施形態のモデル化工程Ｓ２の処理手順の一例を示すフローチャートである。図２６は、図１２に示した高分子材料モデル１１の部分斜視図である。なお、この実施形態では、モデル化工程Ｓ２を除く他の工程Ｓ１〜工程Ｓ６が、前実施形態と同一である。

この実施形態のモデル化工程Ｓ２では、先ず、高分子材料モデル１１のフィラーモデル３及びポリマーモデル５をモデル化して、第１の有限要素モデル３１Ａ（図３０に示す）が作成される（第１工程Ｓ７１）。図２７は、この実施形態の第１工程Ｓ７１の処理手順の一例を示すフローチャートである。

図２８は、三次元空間を示す断面斜視図である。図２６及び図２８に示されるように、第１工程Ｓ７１では、先ず、複数の六面体要素３３で区分された三次元空間３４に、高分子材料モデル１１のフィラーモデル３の輪郭３ｓの座標があてはめられる（工程Ｓ７１１）。

図２６に示されるように、フィラーモデル３の輪郭３ｓは、前実施形態と同様に、高分子材料モデル１１のフィラーモデル３の重心点を中心とする単一球Ｂの表面Ｂｓ（図１２に示す）で設定されている。これにより、フィラーモデル３は、中心粒子４ｃと表面粒子４ｓとを含む多面体状（図１２に示す）から、現実のフィラーの輪郭に近似させることができる。

図２８に示されるように、本実施形態の三次元空間３４は、高分子材料モデル１１（図２６に示す）と同一形状に設定されている。これにより、高分子材料モデル１１を三次元空間３４に重ね合わせることにより、フィラーモデル３の輪郭を、三次元空間３４に容易に定義することができる。なお、三次元空間３４の形状は、高分子材料モデル１１の形状と異ならせてもよい。このような場合は、隣接するフィラーモデル３、３の間隔を維持したまま、任意のフィラーモデル３の輪郭３ｓが、三次元空間３４にあてはめられる。これにより、フィラーモデル３の分散状態を維持したまま、フィラーモデル３の輪郭３ｓを、三次元空間３４に定義することができる。三次元空間３４は、例えば一辺の長さＬ６が５０nm〜１５０nm程度の正方形に形成されるのが望ましい。

図２９は、六面体要素を示す斜視図である。本実施形態の六面体要素３３は、８つの節点３５と、各節点３５、３５間を直線でつなぐ１２本の辺３６とで形成される。また、六面体要素３３は、その一辺の長さＬ７が０．１nm〜３．０nm程度の立方体に形成される。なお、六面体要素３３は、立方体に限定されるわけではなく、例えば、直方体でもよい。

次に、各六面体要素３３の中心３３ｃが、フィラーモデル３の輪郭３ｓよりもフィラーモデル３の中心３ｃ（図２８に示す）側に配されるか否かが判断される（工程Ｓ７１２）。なお、フィラーモデル３が非球形である場合、フィラーモデル３の中心３ｃは、図心として定義される。

この工程Ｓ７１２では、各六面体要素３３の中心３３ｃが、フィラーモデル３の輪郭３ｓよりもフィラーモデル３の中心３ｃ（図２８に示す）側に配されていると判断された場合（工程Ｓ７１２で「Ｙ」）、該六面体要素３３が、フィラーモデル３（図２６に示す）をなすフィラーモデル部３７の六面体要素３３として定義される（工程Ｓ７１３）。一方、各六面体要素３３の中心３３ｃが、フィラーモデル３の輪郭３ｓよりもフィラーモデル３の中心３ｃ側に配されていないと判断された場合（工程Ｓ７１２で「Ｎ」）は、該六面体要素３３が、ポリマーモデル５（図２６に示す）をなすポリマーモデル部３８の六面体要素３３として定義される（工程Ｓ７１４）。そして、フィラーモデル部３７及びポリマーモデル部３８の各六面体要素３３の座標が、コンピュータ１に記憶される（工程Ｓ７１５）。

次に、全ての六面体要素３３がフィラーモデル部３７又はポリマーモデル部３８として定義されたか否かが判断される（工程Ｓ７１６）。この工程Ｓ７１６では、全ての六面体要素３３が、フィラーモデル部３７又はポリマーモデル部３８として定義されたと判断された場合（工程Ｓ７１６で「Ｙ」）、次の第２工程Ｓ７２が実施される。一方、全ての六面体要素３３が、フィラーモデル部３７又はポリマーモデル部３８として定義されていないと判断される場合（工程Ｓ７１６で「Ｎ」）は、工程Ｓ７１２〜工程Ｓ７１６が再度実施される。

図３０は、第１の有限要素モデルを示す断面斜視図である。上記のような工程Ｓ７１２〜工程７１６を経ることにより、三次元空間３４には、フィラーモデル部３７とポリマーモデル部３８とを含む第１の有限要素モデル３１Ａが定義される。この例では、理解しやすいように、フィラーモデル部３７が着色されて表示される。このような第１の有限要素モデル３１Ａは、コンピュータ１に記憶される。

次に、フィラーモデル部３７とポリマーモデル部３８との境界の節点３５ａ（図２９に示す）を含む六面体要素３３を変形させた第２の有限要素モデルが作成される（第２工程Ｓ７２）。図３１は、この実施形態の第２工程Ｓ７２の処理手順の一例を示すフローチャートである。図３２（ａ）は、要素を拡大して示す平面図である。図３２（ｂ）は、変形後の要素を示す平面図である。

本実施形態の第２工程Ｓ７２では、先ず、図３２（ａ）に示されるように、コンピュータ１が、フィラーモデル３の中心３ｃ（図３０に示す）と、境界の節点３５ａとを通る直線４０が定義される（工程Ｓ７２１）。次に、直線４０とフィラーモデル３の輪郭３ｓとの交点４１が求められる（工程Ｓ７２２）。これらの直線４０及び交点４１は、コンピュータ１に記憶される。

次に、コンピュータ１が、境界の節点３５ａを、交点４１に移動させる（工程Ｓ７２３）。本実施形態では、直線４０に沿って、境界の節点３５ａを交点４１に移動させている。この境界の節点３５ａの移動に伴い、図３２（ｂ）に示されるように、境界の節点３５ａを含む六面体要素３３のみを変形させることができる。そして、境界の節点３５ａを含む各六面体要素３３の座標が、コンピュータ１に記憶される（工程Ｓ７２４）。

次に、全ての境界の節点３５ａが、直線４０の交点４１に移動されたか否かが判断される（工程Ｓ７２５）。この工程Ｓ７２５では、全ての境界の節点３５ａが、直線４０の交点４１に移動されたと判断された場合（工程Ｓ７２５で「Ｙ」）、次の変形工程Ｓ３が実施される。一方、全ての境界の節点３５ａが、直線４０の交点４１に移動されていないと判断された場合（工程Ｓ７２５で「Ｎ」）は、新たな境界の節点３５ａが選択され、工程Ｓ７２１〜工程Ｓ７２４が再度実施される。これにより。第２工程Ｓ７２では、全ての境界の節点３５ａを、フィラーモデル部３７の表面の凹凸４３（図３０に示す）が小さくなる向きに移動させることができる。

図３３は、第２の有限要素モデルを示す斜視図である。第２工程Ｓ７２では、上記のような工程Ｓ７２１〜工程Ｓ７２４を経ることにより、境界の節点３５ａを含む六面体要素３３のみを変形させた第２の有限要素モデル３１Ｂ（図３３に示す）が形成される。一方、境界の節点３５ａを含まない六面体要素３３は、変形されることなく、第１の有限要素モデル３１Ａ（図３０に示す）のときと同じ形状が維持されるため、処理時間を大幅に短縮しうる。

各六面体要素３３には、該六面体要素３３を構成する節点３５の番号や、節点３５の座標値等、シミュレーションによる数値解析に必要な情報や、高分子材料及びフィラーの物性に応じた材料定数が定義され、コンピュータ１に記憶される。

なお、第２工程Ｓ７２において、複数の境界の節点３５ａが、同一の前記直線４０上に位置する場合、前実施形態で説明したように、境界の節点３５ａ、３５ａが重複するおそれがある。

このような重複を抑制するために、第２工程Ｓ７２では、前実施形態と同様に、コンピュータ１が、境界の節点３５ａ（図３２（ａ）に示す）を交点４１（図３２（ａ）に示す）に移動させる工程Ｓ７２３の前に、同一の直線４０（図３２（ａ）に示す）上で複数の境界の節点３５ａが重複するか否かを確認する工程と、上記工程で重複すると判断された場合の境界の節点３５ａと交点４１との距離（図示省略）を計算して、該距離に０よりも大かつ１よりも小の係数を乗じた修正距離で直線４０に沿って移動させる工程が含まれるのが望ましい。これにより、第２工程Ｓ７２では、境界の節点３５ａの重複を防ぎつつ、フィラーモデル部３７の表面の凹凸４３（図３０に示す）を小さくすることができる。

図３４には、前実施形態と同様に、図３３に示した有限要素モデル３１（第２の有限要素モデル３１Ｂ）、一方向（例えば、Ｙ軸方向）に０％〜２０％伸長させて、有限要素モデル３１の応力及び伸びの計算結果を示すグラフ（実線）が示されている（実施例２）。

また、図３４には、高分子材料モデル１１の構造に基づいて製造された高分子材料（加硫ゴム）を、０％〜２０％伸長させたときの応力と伸びとの測定結果が破線で示されている（実験例）。さらに、高分子材料にフィラーが規則的に分散されていると仮定し、三次元空間に任意の輪郭を設定して、実施例１と同様のモデル化工程Ｓ２を経て、有限要素モデルが作成された。そして、図３４には、この有限要素モデルを、０％〜２０％伸長させたときの応力及び伸びの計算結果が、一点鎖線で示されている（比較例２）。なお、実施例２、比較例２及び実験例の共通仕様は、以下に記載される仕様を除いて、前実施形態と同一である。
実施例２及び比較例２の有限要素モデル：
要素の数：６４０００個
二次元空間の一辺の長さＬ６：１００nm
要素の一辺の長さＬ７：２．５nm

これらの結果より、実施例２は、比較例２に比べて、実験例に近似しているため、シミュレーション精度を向上しうることが確認できる。これは、実施例２が、シミュレーション工程Ｓ１での構造緩和により、高分子材料モデル１１を、実際の高分子材料の分散状態に近似させることができたことによるものと考えられる。従って、実施例２は、比較例２に比べて、シミュレーション精度を向上しうる。

さらに、この実施形態の有限要素モデル３１は、高分子材料モデル１１が三次元で構成されるため、二次元の有限要素モデル２６（図２１に示す）に比べてシミュレーション精度を向上しうる。

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

３ フィラーモデル
５ ポリマーモデル
８ 仮想空間
１１ 高分子材料モデル
１４ 要素
２６ 有限要素モデル

Claims (9)

1. コンピュータを用いて、フィラーが配合された高分子材料の変形を計算するシミュレーション方法であって、
   前記コンピュータが、前記フィラーをモデル化したフィラーモデルと、前記高分子材料をモデル化したポリマーモデルとが配された仮想空間を用いて分子動力学計算を行うことにより、構造緩和された高分子材料モデルを設定するシミュレーション工程と、
   前記高分子材料モデルに基づいて、有限個の要素でモデル化した有限要素モデルを設定するモデル化工程と、
   前記有限要素モデルの変形計算を行う工程とを含み、
   前記フィラーモデルは、前記フィラーを複数のフィラー粒子でモデル化され、
   前記シミュレーション工程は、前記フィラーモデルの前記フィラー粒子の動きを拘束するフィラー拘束工程、
   前記ポリマーモデルのみを対象として前記分子動力学計算を行うポリマー計算工程、
   前記フィラーモデルの前記フィラー粒子の拘束を解除する工程、及び、
   前記フィラーモデル及び前記ポリマーモデルを対象に前記分子動力学計算を行うフィラー・ポリマー計算工程をさらに含むことを特徴とする高分子材料のシミュレーション方法。
2. 前記ポリマーモデルは、前記高分子材料を、少なくとも一つのポリマー粒子でモデル化され、
   前記フィラーモデルの前記フィラー粒子は、一つの中心粒子と、該中心粒子を中心とする単一球の表面上に中心を有する少なくとも４つの表面粒子とを含む請求項１に記載の高分子材料のシミュレーション方法。
3. 前記シミュレーション工程は、予め定めた粒子間距離以下に接近した前記フィラー粒子と前記ポリマー粒子とを結合させる結合工程を含み、
   前記ポリマー粒子は、前記表面粒子のみに結合する請求項２に記載の高分子材料のシミュレーション方法。
4. 前記フィラー粒子と前記ポリマー粒子との間には、粒子間の距離が、予め定められたカットオフ距離以下になったとき相互作用が生じるポテンシャルが定義され、
   前記中心粒子が関連するカットオフ距離は、前記表面粒子が関連するカットオフ距離と前記単一球の半径との和よりも大きい請求項２又は３に記載の高分子材料のシミュレーション方法。
5. 前記有限要素モデルは、二次元モデルであり、
   前記モデル化工程は、前記高分子材料モデルの二次元断面を取得する工程と、
   前記二次元断面の前記フィラーモデル及び前記ポリマーモデルを、それぞれ正方形又は長方形の要素でモデル化して第１の有限要素モデルを作成する第１工程と、
   前記第１の有限要素モデルにおいて、前記フィラーモデルをなすフィラーモデル部と、前記ポリマーモデルをなすポリマーモデル部との境界の節点を、前記フィラーモデル部の表面の凹凸が小さくなる向きに移動させ、該境界の節点を含む前記要素を変形させた第２の有限要素モデルを作成する第２工程とを含む請求項１乃至４のいずれかに記載の高分子材料のシミュレーション方法。
6. 前記第２工程は、前記各境界の節点を、該境界の節点と前記フィラーモデル部の中心とを通る直線上で移動させる請求項５記載の高分子材料のシミュレーション方法。
7. 前記第２工程は、前記各境界の節点を、前記直線と前記フィラーモデルの輪郭との交点へ移動させる請求項６記載の高分子材料のシミュレーション方法。
8. 前記輪郭は、前記フィラーモデルの重心点を中心とする単一球の表面で設定される請求項７記載の高分子材料のシミュレーション方法。
9. 前記有限要素モデルは、三次元モデルであり、
   前記モデル化工程は、前記高分子材料モデルの前記フィラーモデル及び前記ポリマーモデルを、それぞれ六面体要素でモデル化して第１の有限要素モデルを作成する第１工程と、
   前記第１の有限要素モデルにおいて、前記フィラーモデルをなすフィラーモデル部と、前記ポリマーモデルをなすポリマーモデル部との境界の節点を、前記フィラーモデル部の表面の凹凸が小さくなる向きに移動させ、該境界の節点を含む六面体要素を変形させた第２の有限要素モデルを作成する第２工程とを含む請求項１乃至４のいずれかに記載の高分子材料のシミュレーション方法。