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

Abstract

【課題】計算落ちを効果的に防ぐ。
【解決手段】フィラーと高分子材料との反応を、コンピュータ１を用いて評価するためのシミュレーション方法である。このシミュレーション方法では、コンピュータ１が、予め定めた仮想空間６に配置されたポリマーモデル２とフィラーモデル１２とを用いて分子動力学計算を行うシミュレーション工程Ｓ７を含む。シミュレーション工程Ｓ７は、伸びきり長Ｒ₀を定義しない第１結合ポテンシャルＰ１を、ポリマー粒子４、４間に設定して、分子動力学計算を行う第１計算工程Ｓ７１と、第１計算工程Ｓ７１の後に、伸びきり長Ｒ₀以上の粒子間距離ｒ_ijを許容しない第２結合ポテンシャルＰ２を、ポリマー粒子４、４間に設定して、分子動力学計算を行う第２計算工程Ｓ７４とを含む。
【選択図】図８

Description

本発明は、計算落ちを効果的に防ぐことができる高分子材料のシミュレーション方法に関する。

ゴム等の高分子材料には、カーボンブラックやシリカといったフィラーが配合されている。これらのフィラーの分散性は、ゴムの強度に大きな影響を与えることが明らかになっている。

近年、ゴム配合の開発のために、フィラーと高分子材料との反応を、コンピュータを用いて評価するためのシミュレーション方法（数値計算）が種々提案されている。図１２（ａ）に示されるように、この種のシミュレーション方法では、先ず、フィラーをモデル化したフィラーモデルａが設定される。次に、高分子材料を、複数のポリマー粒子ｂ１でモデル化したポリマーモデルｂが設定される。さらに、フィラーモデルａとポリマー粒子ｂ１との間に、ポテンシャルｆ１が設定される。そして、予め定めた仮想空間に配置されたフィラーモデルａとポリマーモデルｂとを用いて分子動力学計算が行われる。

また、ポリマーモデルｂは、複数のポリマー粒子ｂ１が、ボンドｃで結合されている。このボンドｃには、伸びきり長が定義された結合ポテンシャルｆ２が設定される。一般に、結合ポテンシャルｆ２は、伸びきり長以上の粒子間距離を許容させていない。これは、高分子材料の分子運動に基づくものである。関連する技術としては、次のものがある。

特開２００６−６４６５８号公報

上記のような結合ポテンシャルｆ２が設定されたボンドｃでは、ポリマー粒子ｂ１、ｂ１の粒子間距離を伸びきり長以上に広げることができない。このため、図１２（ｂ）に示されるように、計算当初は、人為的に配置された隣接するポリマーモデルｂ、ｂが、ポリマー粒子ｂ１、ｂ１間で絡まりやすく、一旦絡まると解けにくい。このため、分子動力学計算において、ポリマーモデルｂ、ｂ間で、非常に大きな力が作用し、解析開始当初に計算落ちが発生しやすいという問題があった。

本発明は、以上のような実状に鑑み案出されたもので、伸びきり長を定義しない第１結合ポテンシャルＰ１を、ポリマー粒子間に設定して、分子動力学計算を行う第１計算工程と、第１計算工程の後に、伸びきり長以上の粒子間距離を許容しない第２結合ポテンシャルＰ２を、ポリマー粒子間に設定して、分子動力学計算を行う第２計算工程とを含むことを基本として、計算落ちを効果的に防ぐことができる高分子材料のシミュレーション方法を提供することを主たる目的としている。

本発明のうち請求項１記載の発明は、フィラーと高分子材料との反応を、コンピュータを用いて評価するためのシミュレーション方法であって、前記コンピュータに、前記高分子材料を、複数のポリマー粒子と、該ポリマー粒子間を結合するボンドとでモデル化したポリマーモデルを設定する工程、前記コンピュータに、前記フィラーをモデル化したフィラーモデルを設定する工程、並びに、前記コンピュータが、予め定めた仮想空間に配置された前記ポリマーモデルと前記フィラーモデルとを用いて分子動力学計算を行うシミュレーション工程を含み、前記ボンドには、結合ポテンシャルＰが設定され、前記結合ポテンシャルＰは、伸びきり長を定義しない第１結合ポテンシャルＰ１と、前記伸びきり長以上の粒子間距離を許容しない第２結合ポテンシャルＰ２とを含み、前記シミュレーション工程は、前記第１結合ポテンシャルＰ１を、前記ポリマー粒子間に設定して、前記分子動力学計算を行う第１計算工程と、前記第１計算工程の後に、前記第２結合ポテンシャルＰ２を、前記ポリマー粒子間に設定して、前記分子動力学計算を行う第２計算工程とを含むことを特徴とする。

また、請求項２記載の発明は、前記仮想空間は、互いに向き合う一対の面を少なくとも含み、前記フィラーモデルは、前記仮想空間の前記面としてモデル化されている請求項１に記載の高分子材料のシミュレーション方法である。

また、請求項３記載の発明は、前記第１結合ポテンシャルＰ１は、下記式（１）で定義され、前記第２結合ポテンシャルＰ２は、下記式（２）で定義される斥力ポテンシャルＲと、下記式（３）で定義される引力ポテンシャルＧとの和で定義される請求項１又は２に記載の高分子材料のシミュレーション方法である。



ここで、各定数及び変数は次のとおりである。
ｒ’：平衡長
k、ｋ’：各粒子間のばね定数
ｒ_ij：各粒子間距離
Ｒ₀：伸びきり長
ε：各粒子間に作用する斥力の強度に対応する定数
σ：各粒子間に斥力が作用する距離に関する定数

本発明の高分子材料のシミュレーション方法は、複数のポリマー粒子をボンドで結合させてモデル化したポリマーモデルを設定する工程、各ポリマー粒子間に、結合ポテンシャルＰを設定する工程、及び、予め定めた仮想空間に配置されたフィラーモデルとポリマーモデルとを用いて分子動力学計算を行うシミュレーション工程が含まれる。

ポリマー粒子間の結合ポテンシャルＰは、伸びきり長を定義しない第１結合ポテンシャルＰ１と、伸びきり長以上の粒子間距離を許容しない第２結合ポテンシャルＰ２とを含む。シミュレーション工程は、第１結合ポテンシャルＰ１を、ポリマー粒子間に設定して、分子動力学計算を行う第１計算工程と、第１計算工程の後に、第２結合ポテンシャルＰ２を、ポリマー粒子間に設定して、分子動力学計算を行う第２計算工程とを含む。

第１計算工程では、ポリマー粒子間の距離が伸びきり長に制限されないため、隣接するポリマーモデルがポリマー粒子間で絡まっても解けやすい。これにより、第１計算工程では、ポリマーモデルの絡まりによって作用する大きな力を抑制することができるため、計算落ちの発生を防ぐことができる。したがって、本発明のシミュレーション方法では、円滑に緩和計算することができる。

さらに、第２計算工程では、第１計算工程の後に、伸びきり長以上のポリマー粒子間距離を許容しないで、分子動力学計算が行われる。これにより、第２計算工程では、高分子材料の分子運動に近似させて、緩和計算することができるため、シミュレーション精度を維持することができる。

本実施形態のシミュレーション方法を実行するコンピュータ装置の斜視図である。 本実施形態のシミュレーション方法の処理手順の一例を示すフローチャートである。 ポリマーモデルの概念図である。 図３の部分拡大図である。 仮想空間の概念斜視図である。 フィラーモデルとポリマーモデルとの間の相互ポテンシャルを説明する概念図である。 複数のポリマーモデルが配置された仮想空間を示す概念図である。 シミュレーション工程の処理手順の一例を示すフローチャートである。 （ａ）、（ｂ）は、第１計算工程のポリマーモデルを示す概念図である。 第１計算工程後の仮想空間を示す概念図である。 第２計算工程後の仮想空間の一部を示す概念図である。 （ａ）、（ｂ）は、従来の高分子材料のシミュレーション方法を説明する概念図である。

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

フィラーとしては、例えば、カーボンブラック、シリカ、又は、アルミナ等が含まれる。また、高分子材料としては、例えば、ゴム、樹脂、又は、エラストマー等が含まれる。さらに、本実施形態の高分子材料には、フィラーに対する親和性を高める変性剤が配合される。この変性剤としては、例えば、水酸基、カルボニル基、又は、原子団の官能基等が含まれる。

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

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

ポリマー粒子４は、一つのポリマーモデル２に対して、例えば、１０〜５０００個設けられる。また、図４に示されるように、ボンド５には、結合ポテンシャルＰが設定される。

また、図３及び図４に示されるように、本実施形態のポリマーモデル２は、変性剤をモデル化した変性剤モデル７が含まれる。変性剤モデル７は、少なくとも一つ、本実施形態では１つの変性基粒子８でモデル化される。この変性基粒子８も、ポリマー粒子４にボンド５を介して結合されている。

このようなポリマーモデル２は、ポリマー粒子４、４間、及び、ポリマー粒子４と変性基粒子８との間の各粒子間距離ｒ_ij（図４に示す）が、ボンド５によって安定して保持される直鎖状の三次元構造をなしている。また、ポリマーモデル２は、高分子材料を分子動力学で取り扱うための数値データ（ポリマー粒子４、及び、変性基粒子８の質量、体積、直径、並びに、初期座標などを含む）である。これらの数値データは、コンピュータ１に入力される。なお、粒子間距離ｒ_ijは、各粒子の中心間の距離として定義される。

次に、結合ポテンシャルＰが設定される（工程Ｓ２）。本実施形態の結合ポテンシャルＰは、伸びきり長Ｒ₀を定義しない第１結合ポテンシャルＰ１と、伸びきり長Ｒ₀以上の粒子間距離ｒ_ijを許容しない第２結合ポテンシャルＰ２とが含まれる。

第１結合ポテンシャルＰ１は、下記式（１）で定義される。

ここで、各定数及び変数は、次のとおりである。
ｒ_ij：各粒子間距離
ｋ：各粒子間のばね定数
ｒ’：平衡長

上記式（１）で定義される第１結合ポテンシャルＰ１は、Harmonic型である。Harmonic型とは、いわゆる線形バネが定義され、平衡長ｒ'からのずれの大きさに比例した復元力
が働くポテンシャルである。なお、平衡長ｒ'は、各粒子の中心間の距離として定義される。

上記式（１）において、粒子間距離ｒ_ijと平衡長ｒ'とが等しい場合には、第１結合ポテンシャルＰ１が０となる。また、粒子間距離ｒ_ijが、平衡長ｒ'よりも大である場合は、該粒子間距離ｒ_ijが大きくなるほど、第１結合ポテンシャルＰ１も大きくなる。一方、粒子間距離ｒ_ijが、平衡長ｒ'よりも小である場合は、該粒子間距離ｒ_ijが小さくなるほど、第１結合ポテンシャルＰ１も大きくなる。このように、第１結合ポテンシャルＰ１では、粒子間距離ｒ_ijを平衡長ｒ'に戻そうとする復元力が定義される。

また、第１結合ポテンシャルＰ１は、後述する第２結合ポテンシャルＰ２とは異なり、粒子間距離ｒ_ijの上限値が、伸びきり長Ｒ₀に限定されない。従って、第１結合ポテンシャルＰ１（結合ポテンシャルＰ）は、伸びきり長Ｒ₀以上の粒子間距離ｒ_ijを許容することができる。

上記式（１）のばね定数ｋ及び平衡長ｒ'の各値は、適宜設定することができる。本実施形態では、第１結合ポテンシャルＰ１と第２結合ポテンシャルＰ２とが近似するように設定される。本実施形態では、例えば、ばね定数ｋが９００、平衡長ｒ'が０．９に設定される。

第２結合ポテンシャルＰ２は、下記式（２）で定義される斥力ポテンシャルＲと、下記式３で定義される引力ポテンシャルＧとの和で定義される。

ここで、各定数及び変数は、次に示すものを除いて、上記式（１）と同一である。また、下記の各定数は、Lennard-Jonesポテンシャルのパラメータである。
k’：各粒子間のばね定数
R₀：伸びきり長
ε：各粒子間に定義される斥力ポテンシャルＲの強度に対応する定数
σ：各粒子間に定義される斥力ポテンシャルＲが作用する距離に関する定数（分子動力学の分野では、ＬＪ球の直径と呼ばれる）

上記式（２）において、粒子間距離ｒ_ijが小さいほど、斥力ポテンシャルＲが大きくなる。一方、上記式（３）において、粒子間距離ｒ_ijが大きいほど、引力ポテンシャルＧが大きくなる。従って、第２結合ポテンシャルＰ２では、粒子間の距離ｒ_ijを、斥力ポテンシャルＲと引力ポテンシャルＧとが互いに釣り合う位置に戻そうとする復元力が定義される。

また、上記式（３）では、粒子間距離ｒ_ijが伸びきり長Ｒ₀以上となる場合には、引力ポテンシャルＧが∞に設定される。従って、第２結合ポテンシャルＰ２は、第１結合ポテンシャルＰ１とは異なり、伸びきり長Ｒ₀以上の粒子間距離ｒ_ijを許容しない。このような第２結合ポテンシャルＰ２は、第１結合ポテンシャルＰ１に比べて、高分子材料の分子運動に近似させることができる。

なお、引力ポテンシャルＧのばね定数ｋ’、伸びきり長Ｒ₀、及び、斥力ポテンシャルＲの強度εの値としては、適宜設定することができる。本実施形態では、論文（J. Chem Phys. vol.92, No.8, 5057（1990））に基づいて、ばね定数ｋ’が３０、伸びきり長Ｒ₀が１．５、強度εが１．０、及び、定数σが１．０に設定されている。

次に、図３に示されるように、ポリマーモデル２、２間において、ポリマー粒子４、４間、変性基粒子８、８間及びポリマー粒子４と変性基粒子８との間に、斥力ポテンシャルＲが設定される（工程Ｓ３）。本実施形態の斥力ポテンシャルＲは、上記式（２）で定義される。また、斥力ポテンシャルＲの強度ε及び定数σも、上記式（２）と同一が望ましい。

次に、図５に示されるように、予め定められた体積をもつ仮想空間６が設定される（工程Ｓ３）。仮想空間６は、互いに向き合う少なくとも一対の面１１、１１、本実施形態では、３対の面１１、１１を有する立方体として定義される。これらの各面１１は、ポリマーモデル２（図３に示す）が通過不能に定義される。

また、一対の面１１が直交する方向において、該一対の面１１の間隔Ｄ１（即ち、１辺の長さＬ１）が、ポリマーモデル２の慣性半径（図示省略）の２倍以上、より好ましくは４倍以上に設定されるのが望ましい。これにより、後述する分子動力学計算において、ポリマーモデル２の仮想空間６内での回転運動をスムーズに計算することができる。また、仮想空間６の大きさは、上記論文に基づいて、粒子数密度が０．８５個／σ³程度に設定される。

次に、フィラーをモデル化したフィラーモデル１２が設定される（工程Ｓ４）。このフィラーモデル１２は、本実施形態では、仮想空間６の上下に配された一対の面１１、１１でモデル化される。これにより、一対のフィラーモデル１２、１２は、仮想空間６において、移動不能に定義される。

次に、フィラーモデル１２とポリマーモデル２との間に、相互ポテンシャルＴが設定される（工程Ｓ５）。本実施形態の工程Ｓ５では、図６に示されるように、フィラーモデル１２とポリマー粒子４との間、及び、フィラーモデル１２と変性基粒子８との間に、引力、及び、斥力が作用する相互ポテンシャルＴが設定される。相互ポテンシャルＴは、例えば、下記式（４）で定義される。

ここで、各定数及び変数は次のとおりである。
ρ_wall：相互ポテンシャルＴの壁面の密度に関する定数
ε_wall：相互ポテンシャルＴの強度に対応する定数
σ_wall：仮想空間の面（フィラーモデル）に直交する方向の斥力長さに関する定数
ｒ：フィラーモデルとポリマー粒子との間及びフィラーモデルと変性基粒子との間の距離
ｒ_c：カットオフ距離
なお、フィラーモデル１２と各粒子４、８との間の距離ｒ及びカットオフ距離ｒ_cは、フィラーモデル１２の面と、各粒子４、８の中心との間の最短距離で定義される。

上記式（４）で定義される相互ポテンシャルＴは、上記式（２）で定義される斥力ポテンシャルＲを、フィラーモデル１２の面１１で積分したものである。

また、上記式（４）において、前記距離ｒが、予め定められたカットオフ距離ｒ_c以上となる場合は、相互ポテンシャルＴは作用しない。さらに、距離ｒが相互ポテンシャルの定数σ_wall×２^1/6未満になると、相互ポテンシャルＴは、フィラーモデル１２と各粒子４、８との間に斥力のみを生じさせる。一方、距離ｒが相互ポテンシャルの定数σ_wall×２^1/6を超えると、相互ポテンシャルＴは、フィラーモデル１２と各粒子４、８との間に引力を生じさせることができる。このように、上記式（４）では、フィラーモデル１２と各粒子４、８との間に、引力及び斥力を定義することができる。

また、本実施形態の相互ポテンシャルＴは、下記の相互ポテンシャルＴ１、Ｔ２に区分される。
フィラーモデル１２−ポリマー粒子４：相互ポテンシャルＴ１
フィラーモデル１２−変性基粒子８ ：相互ポテンシャルＴ２

各相互ポテンシャルＴ１、Ｔ２の定数ρ_wall、定数σ_wall、強度ε_wall、カットオフ距離ｒ_cの各値としては、適宜設定することができる。本実施形態では、上記論文に基づいて、次の値が設定される。
相互ポテンシャルＴ１：
ρ_wall＝σ_wall＝１．０
ε_wall＝１
ｒ_c＝１．１２
相互ポテンシャルＴ２：
ρ_wall＝σ_wall＝１．０
ε_wall＝５
ｒ_c＝２．５

上記のように、本実施形態では、相互ポテンシャルＴ２の強度ε_wall（＝５）が、相互ポテンシャルＴ１の強度ε_wall（＝１）よりも大に設定される。これにより、変性基粒子８のフィラーモデル１２に対する力の強度が、ポリマー粒子４のフィラーモデル１２に対する力の強度よりも大に設定される。

また、本実施形態では、相互ポテンシャルＴ２のカットオフ距離ｒ_cが２．５に設定されるため、フィラーモデル１２及び変性基粒子８間において、距離ｒが２^1/6（≒１．１２）から２．５未満の広範囲に亘り、引力を作用させることができる。一方、相互ポテンシャルＴ１は、カットオフ距離ｒ_cが１．１２に設定されるため、フィラーモデル１２及びポリマー粒子４において、斥力のみを作用させることができる。従って、変性基粒子８は、ポリマー粒子４に比べて、引力が優先的に作用する。

このように、各相互ポテンシャルＴ１、Ｔ２は、強度ε_wall及びカットオフ距離ｒ_cが上記の範囲に設定されることにより、変性基粒子８のフィラーモデル１２に対する親和性を、ポリマー粒子４のフィラーモデル１２に対する親和性よりも高くすることができる。

次に、図７に示されるように、仮想空間６の中に、複数のポリマーモデル２が配置される（工程Ｓ６）。本実施形態では、仮想空間６に、２００〜４００本程度のポリマーモデル２が、フィラーモデル１２を構成する一対の各面１１、１１に凝集して配置されている。

次に、フィラーモデル１２と、ポリマーモデル２とを用いて分子動力学計算による緩和計算が行なわれる（シミュレーション工程Ｓ７）。

本実施形態の分子動力学計算では、例えば、仮想空間６について所定の時間、ポリマーモデル２が古典力学に従うものとして、ニュートンの運動方程式が適用される。そして、各時刻でのポリマー粒子４、及び、変性基粒子８の動きが、単位時間毎に追跡される。また、仮想空間６内の粒子、体積及び温度は、夫々一定に保たれる。

本実施形態のシミュレーション工程Ｓ７では、仮想空間６において、フィラーモデル１２を一対の面１１に固定させた状態で、ポリマーモデル２のみを対象に緩和計算が行われる。従って、本実施形態では、フィラーモデル、及び、ポリマーモデルの双方を対象に緩和計算が行なわれていた従来の方法に比べて、計算時間を短縮しうる。

また、フィラーモデル１２は、面１１全体において、該面１１と直角に交わる方向に、各相互ポテンシャルＴ１、Ｔ２（図６に示す）を作用させることができる。これにより、フィラーモデル１２は、ポリマーモデル２に対して、同一方向の相互ポテンシャルＴ１、Ｔ２を面１１全体で作用させることができるため、ポリマーモデル２の緩和計算を効率的に行うことができる。

図８には、本実施形態のシミュレーション工程Ｓ７の具体的な処理手順が示されている。
本実施形態のシミュレーション工程Ｓ７は、先ず、ポリマー粒子４、４間、及び、ポリマー粒子４と変性基粒子８との間に、第１結合ポテンシャルＰ１を設定して、分子動力学計算が行われる（第１計算工程Ｓ７１）。

図７に示したように、ポリマーモデル２がフィラーモデル１２側に凝集して配置される計算当初は、図９（ａ）に示されるように、隣接するポリマーモデル２、２がポリマー粒子４、４間で絡まりやすい。第１計算工程Ｓ７１では、伸びきり長Ｒ₀以上の粒子間距離ｒ_ijを許容して分子動力学計算されるため、図９（ｂ）に示されるように、粒子４、４間の距離を広げることができる。このため、隣接するポリマーモデル２が粒子間で絡まっても解けやすい。

このように、第１計算工程Ｓ７１では、ポリマーモデル２、２同士の絡まりによって作用する大きな力を抑制することができるため、計算落ちの発生を防ぐことができる。従って、本実施形態のシミュレーション方法では、図１０に示されるように、ポリマーモデル２を円滑に緩和することができる。

また、第１計算工程Ｓ７１では、分子動力学計算の単位ステップ当たりの時間刻みを、上記論文に基づく時間刻みよりも小さくするのが望ましい。これにより、単位ステップ毎の各粒子４、８の移動量を小さくすることができるため、ポリマーモデル２、２同士の絡まりによって作用する力が、急激に増大するのを抑制することができる。従って、計算落ちの発生を、効果的に防ぐことができる。

このような作用を効果的に発揮させるために、単位ステップ当たりの時間刻みを、上記論文の時間刻みの１／５〜１／２０に設定されるのが望ましい。なお、単位ステップ当たりの時間刻みが、上記論文の時間刻みの１／５を超えると、上記のような作用を十分に発揮できないおそれがある。逆に、単位ステップ当たりの時間刻みが、上記論文に基づく時間刻みの１／２０未満になると、計算コストが増大するおそれがある。なお、論文に基づく時間刻みが０．００６τであるのに対して、本実施形態の時間刻みが０．０００６τ（論文の１／１０）に設定されている。

次に、図８に示されるように、本実施形態では、第１結合ポテンシャルＰ１の経時変化が単調増加、又は、単調減少しなくなったかが判断される（工程Ｓ７２）。緩和した状態（定常状態）では、第１結合ポテンシャルＰ１が一定の範囲に収束する。従って、第１結合ポテンシャルＰ１の経時変化が単調増加、又は、単調減少しなくなった場合、緩和した状態（定常状態）になったと判断することができる。本実施形態では、例えば、１０ステップ毎に結合ポテンシャルＰ１の平均値を計算し、この平均値に基づいて、第１結合ポテンシャルＰ１の単調増加、又は、単調減少が判断される。

この工程Ｓ７２では、第１結合ポテンシャルＰ１の経時変化が単調増加、又は、単調減少していると判断された場合、単位ステップを進めて（工程Ｓ７３）、第１計算工程Ｓ７１が再度実行される。一方、第１結合ポテンシャルＰ１の経時変化が単調増加、又は、単調減少しなくなったと判断された場合には、次の第２計算工程Ｓ７４が行われる。

次に、第２結合ポテンシャルＰ２を、ポリマー粒子４間に設定して、分子動力学計算が行われる（第２計算工程Ｓ７４）。このような第２計算工程Ｓ７４では、高分子材料の分子運動に近似させて、緩和計算することができるため、シミュレーション精度を向上することができる。

しかも、第２計算工程Ｓ７４に先立ち、第１計算工程Ｓ７１において、ポリマーモデル２、２同士の絡みが解かれて、ポリマーモデル２が効率よく緩和される。これにより、ポリマー粒子４、及び、変性基粒子８が、各相互ポテンシャルＴ１、Ｔ２を作用させて効果的に緩和計算されるため、図１１に拡大して示されるように、該変性基粒子８のフィラーモデル１２に対する近づきやすさ、即ち親和性を正確に評価することができる。

次に、第２計算工程Ｓ７４の処理が、所定のステップ数に達したか否かが判断される（工程Ｓ７５）。この工程Ｓ７５では、コンピュータ１が、所定のステップ数に達していないと判断した場合には、単位ステップを進めて（工程Ｓ７６）、第２計算工程Ｓ７４が再度実行される。一方、所定のステップ数に達したと判断した場合には、次の評価工程Ｓ８が行われる。

次に、シミュレーション工程Ｓ７の結果から、ポリマーモデル２の緩和状態が評価される（評価工程Ｓ８）。本実施形態の評価工程Ｓ８では、変性基粒子８のフィラーモデル１２に対する親和性が評価される。この評価工程Ｓ８において、緩和状態が良好であると判断された場合、本実施形態のシミュレーション方法による一連の処理が終了する。一方、緩和状態が良好でないと判断された場合、例えば、ポリマーモデル２やフィラーモデル１２の設定条件が変更され（工程Ｓ９）、工程Ｓ１〜Ｓ７が行われる。これにより、ポリマーモデル２を確実に緩和させる条件等を、的確に解析することができる。

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

図２及び図８に示す手順に従って、第１結合ポテンシャルＰ１をポリマー粒子間に設定して分子動力学計算を行う第１計算工程と、第１計算工程の後に、第２結合ポテンシャルＰ２をポリマー粒子間に設定して分子動力学計算を行う第２計算工程とを含むシミュレーション工程が実施された（実施例１、２）。なお、実施例２は、実施例１よりも、単位ステップ当たりの時間刻みを小さくしている。

また、比較のために、第２結合ポテンシャルＰ２のみをポリマー粒子間に設定して、分子動力学計算を行うシミュレーション工程が実施された（比較例）。

そして、実施例１、２、及び、比較例において、計算落ちすることなく緩和計算が終了するかが確認された。さらに、緩和計算の開始から終了までの合計計算時間が測定され、実施例１を１００とする指数で表示した。数値が小さい程、合計計算時間が少なく良好である。

なお、実施例１、２、及び、比較例において、第１結合ポテンシャルＰ１、第２結合ポテンシャルＰ２、斥力ポテンシャルＲ、及び、相互ポテンシャルＴの各パラメータは、明細書中の記載通りであり、共通仕様は次のとおりである。
ポリマーモデル：
合計本数：３０本
ポリマーモデル１本当たりのポリマー粒子の個数：１０００個
ポリマーモデル１本当たりの変性基粒子の個数：１個
テストの結果を表１に示す。

テストの結果、実施例１、２では、計算落ちが発生することなく緩和計算が終了することが確認できた。一方、比較例では、計算落ちが頻発した。このため、比較例では、計算落ちが発生する度に、オペレータの操作が必要となった。従って、実施例１、２の合計計算時間は、比較例の合計計算時間よりも大幅に減少することが確認できた。

２ ポリマーモデル
４ ポリマー粒子
１２ フィラーモデル
Ｐ１ 第１結合ポテンシャル
Ｐ２ 第２結合ポテンシャル

Claims (3)

1. フィラーと高分子材料との反応を、コンピュータを用いて評価するためのシミュレーション方法であって、
   前記コンピュータに、前記高分子材料を、複数のポリマー粒子と、該ポリマー粒子間を結合するボンドとでモデル化したポリマーモデルを設定する工程、
   前記コンピュータに、前記フィラーをモデル化したフィラーモデルを設定する工程、並びに、
   前記コンピュータが、予め定めた仮想空間に配置された前記ポリマーモデルと前記フィラーモデルとを用いて分子動力学計算を行うシミュレーション工程を含み、
   前記ボンドには、結合ポテンシャルＰが設定され、
   前記結合ポテンシャルＰは、伸びきり長を定義しない第１結合ポテンシャルＰ１と、前記伸びきり長以上の粒子間距離を許容しない第２結合ポテンシャルＰ２とを含み、
   前記シミュレーション工程は、前記第１結合ポテンシャルＰ１を、前記ポリマー粒子間に設定して、前記分子動力学計算を行う第１計算工程と、
   前記第１計算工程の後に、前記第２結合ポテンシャルＰ２を、前記ポリマー粒子間に設定して、前記分子動力学計算を行う第２計算工程とを含むことを特徴とする高分子材料のシミュレーション方法。
2. 前記仮想空間は、互いに向き合う一対の面を少なくとも含み、
   前記フィラーモデルは、前記仮想空間の前記面としてモデル化されている請求項１に記載の高分子材料のシミュレーション方法。
3. 前記第１結合ポテンシャルＰ１は、下記式（１）で定義され、
   前記第２結合ポテンシャルＰ２は、下記式（２）で定義される斥力ポテンシャルＲと、下記式（３）で定義される引力ポテンシャルＧとの和で定義される請求項１又は２に記載の高分子材料のシミュレーション方法。
   
   
   
   ここで、各定数及び変数は次のとおりである。
   ｒ'：平衡長
   k、ｋ'：各粒子間のばね定数
   ｒ_ij：各粒子間距離
   Ｒ₀：伸びきり長
   ε：各粒子間に作用する斥力の強度に対応する定数
   σ：各粒子間に斥力が作用する距離に関する定数