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

Description

本発明は、高分子材料のシミュレーション方法に関し、詳しくは、フィラーとマトリックスとを含む高分子材料の破壊特性を計算するためのシミュレーション方法に関する。

下記特許文献１は、コンピュータを用いて、フィラーとマトリックスとを含む高分子材料の破壊特性を計算するためのシミュレーション方法を提案している。下記特許文献１のシミュレーション方法では、先ず、高分子材料モデルがコンピュータに設定される。この高分子材料モデルのフィラーとマトリックスとの間には、下記非特許文献１に基づいて、斥力が作用するＬＪポテンシャルが定義されている。そして、下記特許文献１のシミュレーション方法では、高分子材料モデルに歪みを与えて、高分子材料の物理量が計算される。

特開２０１６−８１２９７号公報 Kurt Kremer & Gary S. Grest 著、「Dynamics of entangled linear polymer melts: A molecular-dynamics simulation」、J. Chem Phys. vol.92, No.8, 15 April 1990、p5057-5086

フィラーとマトリックスとの間には、ミクロレベルで引力的な相互作用が存在し、フィラーにマトリックスが物理吸着することが実験的に知られている。しかしながら、上記非特許文献１のＬＪポテンシャルは、斥力的な相互作用であるため、実際の高分子材料の描像とは相反するという問題があった。

本発明は、以上のような実状に鑑み案出されたもので、実際の高分子材料の描像に沿った破壊特性を計算することができる高分子材料モデルのシミュレーション方法を提供することを主たる目的としている。

本発明は、コンピュータを用いて、フィラーとマトリックスとを含む高分子材料の破壊特性を計算するためのシミュレーション方法であって、前記高分子材料に基づいて、数値計算用の高分子材料モデルを前記コンピュータに設定する工程と、前記フィラーと前記マトリックスとの間に斥力系ポテンシャルが作用するとの条件の下、前記高分子材料モデルの応力と歪との関係を、少なくとも歪が零の状態から前記応力が降伏するまでの期間で計算する第１シミュレーション工程と、前記応力と歪との関係に基づいて、前記応力が降伏するまで前記高分子材料モデルに加えられた外部エネルギーを計算する工程と、前記外部エネルギーに基づいて、前記フィラーと前記マトリックスとの間の引力系ポテンシャルを決定する決定工程と、前記フィラーと前記マトリックスとの間に前記引力系ポテンシャルが作用するとの条件の下、前記高分子材料モデルの破壊特性を計算する第２シミュレーション工程とを含むことを特徴とする。

本発明に係る前記高分子材料のシミュレーション方法において、前記決定工程は、前記フィラーと前記マトリックスとの間に前記引力系ポテンシャルが作用するとの条件の下、前記フィラーの表面から前記マトリックスを引き剥がすのに必要なエネルギーを定義する工程と、前記外部エネルギーと、前記引き剥がすのに必要なエネルギーとが等しくなるとの条件の下、前記引力系ポテンシャルを決定する工程とを含んでもよい。

本発明の高分子材料のシミュレーション方法は、前記フィラーと前記マトリックスとの間に前記引力系ポテンシャルが作用するとの条件の下、前記高分子材料モデルの破壊特性を計算する第２シミュレーション工程を含んでいる。従って、本発明の高分子材料のシミュレーション方法は、前記フィラーと前記マトリックスとの間の引力的な相互作用に基づいて、実際の高分子材料の描像に沿った破壊特性を計算することができる。

前記引力系ポテンシャルは、前記フィラーと前記マトリックスとの間に斥力系ポテンシャルが作用するとの条件の下で計算された外部エネルギーに基づいて決定される。従って、本発明の高分子材料のシミュレーション方法は、上記非特許文献１に基づく斥力系ポテンシャルを、引力系ポテンシャルに容易に変換することができる。

高分子材料のシミュレーション方法を実行するためのコンピュータの一例を示す斜視図である。 マトリックスの一例を示す構造式である。 高分子材料のシミュレーション方法の処理手順の一例を示すフローチャートである。 モデル設定工程の処理手順の一例を示すフローチャートである。 高分子材料モデルの一例を示す概念図である。 フィラーモデルの一例を示す概念図である。 マトリックスモデルの一例を示す概念図である。 フィラーモデル及びマトリックスモデルを拡大して示す概念図である。 第１シミュレーション工程の処理手順の一例を示すフローチャートである。 高分子材料モデルに歪を与える工程の一例を説明する概念図である。 高分子材料モデルの応力−歪曲線の一例を示すグラフである。 フィラーモデルと、粗視化粒子モデルとの間の動径分布関数の一例を示すグラフである。 決定工程の処理手順の一例を示すフローチャートである。 ＬＪポテンシャルＵ_LJ（ポテンシャルＵ_FP）と、フィラーモデルの表面を構成するフィラー粒子モデルの重心からの距離ｒとの関係を示すグラフである。 第２シミュレーション工程の処理手順の一例を示すフローチャートである。 （ａ）〜（ｅ）は、フィラーとマトリックスとの間に引力系ポテンシャルが作用するとの条件の下で、高分子材料モデルの破壊特性を計算した実施例の結果を示す概念図である。 （ａ）〜（ｅ）は、フィラーとマトリックスとの間に斥力系ポテンシャルが作用するとの条件の下で、高分子材料モデルの破壊特性を計算した比較例の結果を示す概念図である。

以下、本発明の実施の一形態が図面に基づき説明される。
本実施形態の高分子材料のシミュレーション方法（以下、単に「シミュレーション方法」ということがある。）では、コンピュータを用いて、フィラーとマトリックスとを含む高分子材料の破壊特性が計算される。

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

フィラーとしては、例えば、シリカ、カーボンブラック、又は、アルミナ等が採用される。マトリックスとしては、例えば、天然ゴム、合成ゴム、又は、樹脂等が採用される。図２は、マトリックスの一例を示す構造式である。本実施形態のマトリックスとしては、cis-１，４ポリブタジエン（以下、単に「ポリブタジエン」ということがある）が例示される。このポリブタジエンを構成する分子鎖は、メチレン基（−ＣＨ_２−）とメチン基（−ＣＨ−）とからなるモノマー｛−［ＣＨ_２−ＣＨ＝ＣＨ−ＣＨ_２］−｝が、重合度ｎで連結されて構成されている。なお、マトリックスには、ポリブタジエン以外のものが用いられてもよい。

図３は、本実施形態のシミュレーション方法の処理手順の一例を示すフローチャートである。本実施形態のシミュレーション方法は、先ず、高分子材料に基づいて、数値計算用の高分子材料モデルが、コンピュータ１に設定される（モデル設定工程Ｓ１）。図４は、モデル設定工程Ｓ１の処理手順の一例を示すフローチャートである。

本実施形態のモデル設定工程Ｓ１では、先ず、予め定められた空間（セル）が定義される（工程Ｓ１１）。図５は、高分子材料モデルの一例を示す概念図である。空間２は、後述する分子動力学計算において、計算対象となる領域である。本実施形態の空間２は、例えば、互いに向き合う少なくとも一対、本実施形態では３対の平面３、３を有する立方体として定義される。空間２の内部には、後述するフィラーモデル５、及び、マトリックスモデル１１が複数個配置される。

空間２の各平面３には、上記特許文献１の空間と同様に、周期境界条件が定義されている。また、一対の平面３、３の間隔（即ち、１辺の長さＬ）は、上記特許文献１の空間と同様に設定されるのが望ましい。空間２の内部には、上記特許文献１の空間と同様に、後述するフィラーモデル５及びマトリックスモデル１１の追跡等に用いるための複数の小領域（図示省略）が定義されている。空間２は、コンピュータ１に記憶される。

次に、本実施形態のモデル設定工程Ｓ１では、フィラーをモデル化したフィラーモデル５が設定される（工程Ｓ１２）。図６は、フィラーモデル５の一例を示す概念図である。

フィラーモデル５は、上記特許文献１と同様に、複数のフィラー粒子モデル６と、隣接するフィラー粒子モデル６、６間を結合する結合鎖モデル（図示省略）とが含まれている。フィラー粒子モデル６は、分子動力学計算において、運動方程式の質点として取り扱われる。即ち、フィラー粒子モデル６には、質量、体積、直径、電荷又は初期座標などのパラメータが定義される。

工程Ｓ１２では、上記特許文献１と同様に、走査型透過電子顕微鏡によって取得された高分子材料の三次元画像に基づいて、フィラー粒子モデル６（フィラーモデル５）が空間２（図５に示す）内に配置されるのが望ましい。これにより、工程Ｓ１２では、実際のフィラーの分布に近似させて、フィラー粒子モデル６を空間２内に配置することができる。

結合鎖モデル（図示省略）は、上記特許文献１と同様に、下記式（１）で定義されるポテンシャル（以下、「ＬＪポテンシャルＵ_LJ（ｒ_ij）」ということがある。）と、下記式（２）で定義される結合ポテンシャルＵ_FENEとの和で示されるポテンシャルで定義される。これにより、結合鎖モデルには、フィラー粒子モデル６、６間の距離ｒ_ijを、ＬＪポテンシャルＵ_LJ（ｒ_ij）と結合ポテンシャルＵ_FENEとが互いに釣り合う位置に戻そうとする復元力が定義される。

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

ＬＪポテンシャルＵ_LJ（ｒ_ij）及び結合ポテンシャルＵ_FENEの各パラメータ（ＬＪポテンシャルの強度ε、伸びきり長Ｒ₀、粒子の直径σ及びカットオフ距離ｒ_c）については、上記特許文献１と同様に、上記非特許文献１に基づいて設定される。フィラーモデル５は、コンピュータ１に記憶される。

次に、本実施形態のモデル設定工程Ｓ１では、マトリックス（分子鎖）をモデル化したマトリックスモデル１１が設定される（工程Ｓ１３）。図７は、マトリックスモデル１１の一例を示す概念図である。図８は、フィラーモデル５及びマトリックスモデル１１を拡大して示す概念図である。

本実施形態のマトリックスモデル１１は、粗視化モデルとして構成される。マトリックスモデル１１は、上記特許文献１と同様に、マトリックス（分子鎖）のモノマー又はモノマーの一部分をなす構造単位１２（図２に示す）が、直径を持った球で表現される粒子モデル（以下、「粗視化粒子モデル」ということがある。）１３に置換される。粗視化粒子モデル１３は、分子動力学計算において、運動方程式の質点として取り扱われる。即ち、粗視化粒子モデル１３には、例えば、質量、体積、直径又は電荷などのパラメータが定義される。

マトリックスモデル１１には、隣接する粗視化粒子モデル１３、１３間を結合する結合鎖モデル１４が含まれている。結合鎖モデル１４は、例えば、上記式（１）のＬＪポテンシャルＵ_LJ（ｒ_ij）と、上記式（２）の結合ポテンシャルＵ_FENEとの和で示されるポテンシャルＰ２（図８に示す）で設定される。ＬＪポテンシャルＵ_LJ（ｒ_ij）及び結合ポテンシャルＵ_FENEの各パラメータについては、上記特許文献１と同様に、上記非特許文献１に基づいて設定される。このような結合鎖モデル１４により、粗視化粒子モデル１３を伸縮自在に拘束した直鎖状のマトリックスモデル１１を設定することができる。

工程Ｓ１３では、上記特許文献１と同様に、高分子材料の三次元画像に基づいて、マトリックスモデル１１が空間２（図５に示す）内に配置されるのが望ましい。これにより、工程Ｓ１３では、実際のマトリックスの分布に近似させて、マトリックスモデル１１を空間２内に配置することができる。マトリックスモデル１１は、コンピュータ１に記憶される。

次に、本実施形態のモデル設定工程Ｓ１では、隣接するマトリックスモデル１１、１１間に、ポテンシャル（引力系ポテンシャル）Ｐ３が定義される（工程Ｓ１４）。本実施形態の工程Ｓ１４では、隣接するマトリックスモデル１１、１１の粗視化粒子モデル１３、１３間に、ポテンシャルＰ３がそれぞれ定義される。ポテンシャルＰ３は、上記式（１）のＬＪポテンシャルＵ_LJ（ｒ_ij）によって定義される。

本実施形態のポテンシャルＰ３のカットオフ距離ｒ_cは、例えば、論文（Ali Makke, Michel Perez, Olivier Lame and Jean-Louis Barrat著、 「Mechanical testing of glassy and rubbery polymers in numerical simulations: Role of boundary conditions in tensile stress experiments 」J. Chem. Phys. vol.131, 014904 (2009)）に基づいて、２．５σに設定される。従って、ポテンシャルＰ３は、マトリックス（分子鎖）間に作用する引力系ポテンシャルとして定義される。これにより、後述の決定工程Ｓ４において、後述の第１シミュレーション工程Ｓ２で実施される高分子材料モデル１６の伸長計算において、隣接するマトリックスモデル１１、１１の離間を防ぐことができるため、フィラーモデル５と、マトリックスモデル１１との間で引き剥がされるときのエネルギー（外部エネルギーＷ１）を優先的に計算することができる。また、ポテンシャルＰ３のＬＪポテンシャルの強度εは、上記論文に基づいて、１．０に設定される。ポテンシャルＰ３は、コンピュータ１に記憶される。

次に、本実施形態のモデル設定工程Ｓ１では、隣接するフィラーモデル５、５間に、ポテンシャル（斥力系ポテンシャル）Ｐ４が定義される（工程Ｓ１５）。本実施形態の工程Ｓ１５では、隣接するフィラーモデル５のフィラー粒子モデル６、６間に、ポテンシャルＰ４がそれぞれ定義される。ポテンシャルＰ４は、上記式（１）のＬＪポテンシャルＵ_LJ（ｒ_ij）によって定義される。

本実施形態のポテンシャルＰ４のカットオフ距離ｒ_cは、上記論文に基づいて、２^1/6σに設定される。従って、ポテンシャルＰ４は、フィラー間に作用する斥力系ポテンシャルとして定義される。また、ポテンシャルＰ４のＬＪポテンシャルの強度εは、上記論文に基づいて、１．０に設定される。ポテンシャルＰ４は、コンピュータ１に記憶される。

次に、本実施形態のモデル設定工程Ｓ１では、フィラーモデル５と、マトリックスモデル１１との間に、ポテンシャル（斥力系ポテンシャル）Ｐ５が定義される（工程Ｓ１６）。本実施形態の工程Ｓ１６では、フィラーモデル５のフィラー粒子モデル６と、マトリックスモデル１１の粗視化粒子モデル１３との間に、ポテンシャルＰ５がそれぞれ定義される。ポテンシャルＰ５は、上記式（１）のＬＪポテンシャルＵ_LJ（ｒ_ij）によって定義される。

本実施形態のポテンシャルＰ５のカットオフ距離ｒ_cは、上記論文に基づいて、２^1/6σに設定される。従って、ポテンシャルＰ５は、フィラーとマトリックスとの間に作用する斥力系ポテンシャルとして定義される。また、ポテンシャルＰ５のＬＪポテンシャルの強度εは、上記論文に基づいて、１．０に設定される。ポテンシャルＰ４は、コンピュータ１に記憶される。

次に、本実施形態のモデル設定工程Ｓ１では、フィラーモデル５及びマトリックスモデル１１に、摩擦係数を含んだ運動方程式が定義される（工程Ｓ１７）。運動方程式は、上記特許文献１に基づいて定義される。運動方程式は、コンピュータ１に記憶される。

次に、本実施形態のモデル設定工程Ｓ１では、分子動力学計算に基づいて、高分子材料モデル１６の構造緩和が計算される（工程Ｓ１８）。本実施形態の工程Ｓ１８では、先ず、フィラーモデル５（フィラー粒子モデル６）の位置（座標）を固定して、マトリックスモデル１１（粗視化粒子モデル１３）のみを対象に、構造緩和が計算される。次に、工程Ｓ１８では、フィラーモデル５の固定を解除して、フィラーモデル５及びマトリックスモデル１１を対象に、構造緩和が計算される。これにより、工程Ｓ１８では、高分子材料のフィラーの領域から、フィラーモデル５が大きく位置ずれするのを防ぎつつ、高分子材料モデル１６の構造緩和が計算されうる。

本実施形態の分子動力学計算では、例えば、空間２について所定の時間、フィラーモデル５又はマトリックスモデル１１が古典力学に従うものとして、工程Ｓ１８で定義された運動方程式が適用される。そして、各時刻でのマトリックスモデル１１の動きが、単位時間毎に追跡される。本実施形態では、空間２において、圧力及び温度が一定、又は、体積及び温度が一定に保たれている。これにより、工程Ｓ１８では、実際の高分子材料の分子運動に近似させて、マトリックスモデル１１の初期配置が、精度よく緩和される。このような構造緩和の計算には、例えば、（株）ＪＳＯＬ社製のソフトマテリアル総合シミュレーター（Ｊ−ＯＣＴＡ）に含まれている分子動力学エンジンＣＯＧＮＡＣが用いられるのが望ましい。

工程Ｓ１８では、マトリックスモデル１１の初期配置が十分に緩和されるまで、シミュレーションの単位時間毎に、高分子材料モデル１６（図５に示す）の構造緩和が計算される。これにより、本実施形態では、マトリックスモデル１１の平衡状態（構造が緩和した状態）が計算され、数値計算用の高分子材料モデル１６が設定されうる。高分子材料モデル１６は、コンピュータ１に記憶される。

次に、本実施形態のシミュレーション方法は、コンピュータ１が、フィラーとマトリックスとの間に斥力系ポテンシャルが作用するとの条件の下、高分子材料モデル１６の応力と歪との関係を計算する（第１シミュレーション工程Ｓ２）。本実施形態の第１シミュレーション工程Ｓ２では、少なくとも歪が零の状態から、応力が降伏するまでの期間で、高分子材料モデル１６の応力と歪との関係が計算される。図９は、第１シミュレーション工程Ｓ２の処理手順の一例を示すフローチャートである。

本実施形態の第１シミュレーション工程Ｓ２では、先ず、歪が零の状態の高分子材料モデル１６の応力が計算される（工程Ｓ２１）。高分子材料モデル１６の応力は、上記した分子動力学エンジンＣＯＧＮＡＣによって計算することができる。歪が零の状態の高分子材料モデル１６の応力は、コンピュータ１に記憶される。

次に、本実施形態の第１シミュレーション工程Ｓ２では、高分子材料モデル１６に歪が与えられる（工程Ｓ２２）。図１０は、高分子材料モデルに歪を与える工程の一例を説明する概念図である。

本実施形態の工程Ｓ２２では、ｚ軸方向において、高分子材料モデル１６の一端１６ａ及び他端１６ｂを離間させて、高分子材料モデル１６の伸長が計算される。本実施形態では、上記特許文献１と同様に、分子動力学計算の単位時間（ＭＤステップ）あたりの歪が１０−５以上の速度Ｖ１で、高分子材料モデル１６をアフィン変形に基づいて伸長させている。速度Ｖ１の好ましい範囲については、上記特許文献１と同様に設定されうる。また、工程Ｓ２２では、フィラーモデル５及びマトリックスモデル１１の熱運動（分子運動）が計算される。

次に、本実施形態の第１シミュレーション工程Ｓ２では、高分子材料モデル１６の応力が計算される（工程Ｓ２３）。工程Ｓ２３では、工程Ｓ２２で与えられた歪に基づいて、高分子材料モデル１６の応力が単位時間毎に計算される。さらに、工程Ｓ２３では、高分子材料モデル１６の伸長によって計算されるフィラー粒子モデル６、及び、粗視化粒子モデル１３の動きが、単位時間毎に追跡される。

工程Ｓ２３では、高分子材料モデル１６に与えられた歪、各ポテンシャルＰ１〜Ｐ５（図８に示す）及び運動方程式に基づいて、フィラーモデル５及びマトリックスモデル１１の熱運動が計算される。この熱運動が計算されることにより、空間２には、粗視化粒子モデル１３及びフィラー粒子モデル６が存在しない小領域（図示省略）が形成される。このような小領域は、高分子材料モデル１６に形成された空孔（ボイド）として定義される。高分子材料モデル１６の応力、及び、フィラー粒子モデル６、及び、粗視化粒子モデル１３の各座標値は、コンピュータ１に記憶される。

次に、本実施形態の第１シミュレーション工程Ｓ２では、高分子材料モデル１６に与えられた歪が、予め定められた上限値以上か否かが判断される（工程Ｓ２４）。歪の上限値については、少なくとも応力が降伏する際の歪よりも大きければ、適宜設定することができる。

工程Ｓ２４において、高分子材料モデル１６に与えられた歪が上限値以上であると判断された場合（工程Ｓ２４において、「Ｙ」）、次の工程Ｓ２５が実施される。他方、工程Ｓ２４において、高分子材料モデル１６に与えられた歪が上限値未満であると判断された場合（工程Ｓ２４において、「Ｎ」）、単位時間を一つ進めて（工程Ｓ２６）、工程Ｓ２２〜工程Ｓ２４が実施される。従って、第１シミュレーション工程Ｓ２では、少なくとも歪が零の状態から、応力が降伏するまでの期間で、高分子材料モデル１６の応力を、単位時間毎に計算することができる。

次に、本実施形態の第１シミュレーション工程Ｓ２では、高分子材料モデル１６の応力と歪との関係が求められる（工程Ｓ２５）。本実施形態の工程Ｓ２５では、高分子材料モデル１６に与えられた歪ごとに、高分子材料モデル１６の応力がプロットされることによって、高分子材料モデル１６の応力−歪曲線が求められる。図１１は、高分子材料モデル１６の応力−歪曲線の一例を示すグラフである。

この高分子材料モデル１６の応力−歪曲線において、高分子材料モデル１６の応力が降伏する降伏点Ｑ１で、高分子材料モデル１６に空孔（ボイド）が形成されている。高分子材料モデル１６の応力と歪との関係（応力−歪曲線）は、コンピュータ１に記憶される。

次に、本実施形態の第１シミュレーション工程Ｓ２では、高分子材料モデル１６の応力が降伏した時点（降伏点Ｑ１）において、フィラーモデル５の表面を構成するフィラー粒子モデル６と粗視化粒子モデル１３との間の動径分布関数が求められる（工程Ｓ２７）。本実施形態の動径分布関数は、応力が降伏した時点（降伏点Ｑ１）で記憶されたフィラー粒子モデル６及び粗視化粒子モデル１３の座標値に基づいて計算される。図１２は、フィラーモデル５と、粗視化粒子モデル１３との間の動径分布関数の一例を示すグラフである。図１２において、距離ｒ＝０は、フィラーモデル５の表面を構成するフィラー粒子モデル６の重心を示しており、距離ｒ₁は、２^1/6σである。

本実施形態の動径分布関数は、フィラーモデル５の表面を構成するフィラー粒子モデル６からの距離ｒ（σ）において、粗視化粒子モデル１３が存在する確率密度を示す関数である。動径分布関数の計算方法については、例えば、特許文献（特開２０１５−５６００２号公報）に記載されているとおりである。

図１２に示されるように、動径分布関数ｇ（ｒ）の増分が大きくなる変曲点Ｑ２の距離ｒよりも大きい領域に、粗視化粒子モデル１３が主として存在している。このため、変曲点Ｑ２の距離ｒ₂以下の領域に、高分子材料モデル１６に空孔（ボイド）が形成されている。従って、変曲点Ｑ２での距離ｒ₂を、フィラーモデル５の表面の距離（即ち、フィラーモデル５の表面を構成するフィラー粒子モデル６の表面までの距離）ｒ₁で減じた剥離距離Ｒは、フィラーモデル５のフィラー粒子モデル６からマトリックスモデル１１の粗視化粒子モデル１３を引き剥がすのに必要な距離（以下、単に「粗視化粒子モデルの剥離距離」ということがある。）として定義することができる。動径分布関数は、コンピュータ１に記憶される。

次に、本実施形態のシミュレーション方法は、コンピュータ１が、応力と歪との関係に基づいて、応力が降伏するまで高分子材料モデル１６に加えられた外部エネルギーを計算する（工程Ｓ３）。外部エネルギーＷ１は、図１１に示した高分子材料モデル１６の応力−歪曲線において、高分子材料モデル１６の応力σ_strを、歪ε_strが零の状態から応力が降伏するまで（降伏点Ｑ１まで）の期間で積分した値（図において、外部エネルギーＷ１に対応する領域に色付けしている）に等しい。外部エネルギーＷ１は、下記式（３）によって定義することができる。

ここで、各パラメータは、次のとおりである。
Ａ：フィラー粒子モデルの重心を通る断面積の総和
ｒ：フィラーモデルの表面を構成するフィラー粒子モデルの重心からの距離
σ_str：高分子材料モデルの応力
Ｎ_fil：フィラー粒子モデルの合計個数
ｒ₁：フィラーモデルの表面の距離（即ち、フィラーモデルの表面を構成するフィラー粒子モデルの表面までの距離）
ｒ₂：フィラーモデルの表面からポリマーが引き剥がされた距離（図１２に示した変曲点Ｑ２での距離）

上記式（３）において、σ_str（ｒ）は、フィラーモデル５の表面を構成するフィラー粒子モデル６の重心を基準とする距離ｒでの高分子材料モデルの応力を示している。この応力σ_str（ｒ）を、距離ｒ₁から距離ｒ₂までの範囲（即ち、剥離距離Ｒ（図１２に示す））で積分し、さらに、フィラー粒子モデル６の重心を通る断面積の総和Ａ及びフィラー粒子モデル６の合計個数Ｎ_filを乗じることで、外部エネルギーＷ１を求めることができる。上記式（３）の各パラメータは、いずれも既知の値である。このため、外部エネルギーＷ１は、一意に定められる。外部エネルギーＷ１は、コンピュータ１に記憶される。

次に、本実施形態のシミュレーション方法は、コンピュータ１が、外部エネルギーＷ１に基づいて、フィラーとマトリックスとの間の引力系ポテンシャルを決定する（決定工程Ｓ４）。図１３は、決定工程Ｓ４の処理手順の一例を示すフローチャートである。

本実施形態の決定工程Ｓ４では、先ず、フィラーの表面からマトリックスを引き剥がす（即ち、高分子材料モデル１６に空孔（ボイド）を形成する）のに必要なエネルギー（以下、単に「剥離エネルギー」ということがある。）Ｗ２が定義される（工程Ｓ４１）。本実施形態の剥離エネルギーＷ２は、フィラーとマトリックスとの間に引力系ポテンシャルが作用するとの条件の下、下記式（４）で定義することができる。

ここで、各パラメータは、次のとおりである。
Ｕ_FP：フィラー粒子モデルと粗視化粒子モデルとの間のポテンシャル（引力系ポテンシャル）
Ｎ_fil：フィラー粒子モデルの合計個数
Ｎ_{pol_particle}：フィラーモデル近傍の粗視化粒子モデルの合計個数
ｒ：フィラーモデルの表面を構成するフィラー粒子モデルの重心からの距離
ｒ₁：フィラーモデルの表面の距離（即ち、フィラーモデルの表面を構成するフィラー粒子モデルの表面までの距離）
ｒ₂：フィラーモデルの表面からポリマーが引き剥がされた距離（図１２に示した変曲点Ｑ２での距離）

上記式（４）において、ポテンシャルＵ_FPは、フィラー粒子モデル６と粗視化粒子モデル１３との間に定義されたポテンシャルＰ５（図８に示す）に対応するものである。ポテンシャルＵ_FPは、上記式（１）のＬＪポテンシャルＵ_LJ（ｒ_ij）で定義される。上記式（１）のカットオフ距離ｒ_cについては、モデル設定工程Ｓ１で設定された斥力系ポテンシャルのカットオフ距離ｒ_c（本例では、２^1/6σ）よりも大きい値（本実施形態では、１．５σ又は２．５σ）に設定される。これにより、ポテンシャルＵ_FPは、引力系ポテンシャルとして定義される。

上記式（４）において、フィラー粒子モデルの合計個数Ｎ_{pol_particle}は、フィラーモデル５の表面を構成するフィラー粒子モデル６の外面から、粗視化粒子モデルの半径の２〜５倍（本実施形態では、３倍）の範囲に存在する粗視化粒子モデルの合計個数として定義される。この合計個数Ｎ_{pol_particle}は、モデル設定工程Ｓ１の工程Ｓ１８（図４）で構造緩和された高分子材料モデル１６に基づいて特定される。

上記式（４）において、ポテンシャルＵ_FP（ＬＪポテンシャルＵ_LJ（ｒ_ij））のカットオフ距離ｒ_cには、上記論文のカットオフ距離ｒ_cとは異なる値（本実施形態では、１．５σ又は２．５σ）に設定されたため、ＬＪポテンシャルの強度εを、上記論文に基づいて設定することができない。従って、上記式（４）において、ＬＪポテンシャルの強度εは、変数（未知のパラメータ）として定義される。上記式（４）は、下記式（５）を経て、下記式（６）のように、変形することができる。

図１４は、ＬＪポテンシャルＵ_LJ（ポテンシャルＵ_FP）と、フィラーモデルの表面を構成するフィラー粒子モデルの重心からの距離ｒとの関係を示すグラフである。図１４から明らかなように、Ｕ_FP(r₁)=−ε、Ｕ_FP(r₂)＝０となる。これにより、一つの粗視化粒子モデル１３をフィラーモデル５から引き剥がすのに必要なエネルギーは、εとなる。そして、このエネルギーεに、フィラー粒子モデルの合計個数Ｎ_fil、及び、フィラーモデル近傍の粗視化粒子モデルの合計個数Ｎ_{pol_particle}を乗じることで、フィラーの表面からマトリックスを引き剥がすのに必要なエネルギー（剥離エネルギー）Ｗ２を定義することができる。

剥離エネルギーＷ２及び外部エネルギーＷ１は、いずれも、高分子材料モデル１６に空孔（ボイド）を形成するのに必要なエネルギーを示している。従って、剥離エネルギーＷ２は、外部エネルギーＷ１と等しくなると仮定することができる。剥離エネルギーＷ２は、コンピュータ１に記憶される。

次に、本実施形態の決定工程Ｓ４は、外部エネルギーＷ１と、フィラーの表面からマトリックスを引き剥がすのに必要なエネルギー（剥離エネルギー）Ｗ２とが等しくなるとの条件の下、引力系ポテンシャルが決定される（工程Ｓ４２）。

上述したように、外部エネルギーＷ１は、上記式（３）の各パラメータがいずれも既知であるため、一意に定められている。他方、剥離エネルギーＷ２の未知のパラメータは、ポテンシャルＵ_FP（ＬＪポテンシャルＵ_LJ（ｒ_ij））の強度εのみである。従って、剥離エネルギーＷ２は、外部エネルギーＷ１と等しくなるとの条件の下で（Ｗ１＝Ｗ２）、引力系ポテンシャルのＬＪポテンシャルの強度εを特定することができる。

そして、ポテンシャルＵ_FP（上記式（１）のＬＪポテンシャルＵ_LJ（ｒ_ij））に、特定された強度ε、及び、カットオフ距離ｒ_c（本実施形態では、１．５σ又は２．５σ）が代入されることにより、引力系ポテンシャルが決定される。

このように、本実施形態の決定工程Ｓ４では、フィラーとマトリックスとの間に斥力系ポテンシャルが作用するとの条件の下で計算された外部エネルギーＷ１に基づいて、引力系ポテンシャルを決定することができる。これにより、本実施形態のシミュレーション方法は、従来の上記論文に基づく斥力系ポテンシャルを、引力系ポテンシャルに容易に変換することができる。従って、本実施形態のシミュレーション方法は、実験等に基づいて引力系ポテンシャルを求める必要がないため、計算コストを低減することができる。引力系ポテンシャルは、コンピュータ１に記憶される。

次に、本実施形態のシミュレーション方法は、コンピュータ１が、フィラーとマトリックスとの間に引力系ポテンシャルが作用するとの条件の下、高分子材料モデルの破壊特性を計算する（第２シミュレーション工程Ｓ５）。本実施形態では、モデル設定工程Ｓ１で設定された高分子材料モデル１６において、フィラーモデル５とマトリックスモデル１１との間に引力系ポテンシャルを定義して、高分子材料モデル１６の破壊特性が計算される。図１５は、第２シミュレーション工程Ｓ５の処理手順の一例を示すフローチャートである。

本実施形態の第２シミュレーション工程Ｓ５では、フィラーモデル５のフィラー粒子モデル６と、マトリックスモデル１１の粗視化粒子モデル１３との間に、引力系ポテンシャルが定義される（工程Ｓ５１）。工程Ｓ５１では、フィラーモデル５のフィラー粒子モデル６と、マトリックスモデル１１の粗視化粒子モデル１３との間に定義された斥力系ポテンシャルＰ５（図８に示す）に代えて、決定工程Ｓ４で特定された強度ε、及び、カットオフ距離ｒ_c（本実施形態では、１．５σ又は２．５σ）がそれぞれ代入されたポテンシャル（引力系ポテンシャル）Ｕ_FPが定義される。再定義されたポテンシャルＵ_FPは、コンピュータ１に記憶される。

次に、本実施形態の第２シミュレーション工程Ｓ５では、上記特許文献１のシミュレーション方法と同様の手順で、高分子材料モデル１６を伸長させて（工程Ｓ５２）、高分子材料モデル１６の物理量が計算される（工程Ｓ５３）。

工程Ｓ５２では、高分子材料モデル１６の一端１６ａ及び他端１６ｂが互いに離間するように、高分子材料モデル１６の伸長が計算される。なお、高分子材料モデル１６を伸長させる速度Ｖ１、及び、高分子材料モデル１６に与えられる歪については、上記特許文献１と同一範囲内に定義することができる。

工程Ｓ５３では、伸ばされた高分子材料モデル１６の大きさが拘束され、その歪みを一定に保持した状態で物理量が計算される。工程Ｓ５３では、第１シミュレーション工程Ｓ２と同様に、高分子材料モデル１６に与えられた歪み、各ポテンシャルＰ１〜Ｐ４（図８に示す）、引力系ポテンシャル（フィラーとマトリックスとの間の引力系ポテンシャル）Ｕ_FP、及び、運動方程式に基づいて、フィラーモデル５及びマトリックスモデル１１の熱運動が計算される。この熱運動が計算されることによって、空間２には、粗視化粒子モデル１３及びフィラー粒子モデル６が存在しない小領域（図示省略）が形成される。このような小領域は、高分子材料モデル１６に形成された空孔（ボイド）２１（図１６に示す）として定義され、高分子材料モデル１６の破壊が計算される。

第２シミュレーション工程Ｓ５では、フィラーとマトリックスとの間に引力系ポテンシャル（ポテンシャルＵ_FP）が作用するとの条件の下、高分子材料モデル１６の破壊が計算されている。従って、本実施形態の第２シミュレーション方法は、フィラーとマトリックスとの間の引力的な相互作用に基づいて、実際の高分子材料の描像に沿った破壊特性を計算することができる。

次に、本実施形態のシミュレーション方法では、高分子材料モデル１６の物理量（破壊特性）が、良好か否かが判断される（工程Ｓ６）。工程Ｓ６では、例えば、高分子材料モデル１６の応力、空孔（ボイド）２１（図１６に示す）の総量、又は、空孔（ボイド）２１の偏り等に基づいて、高分子材料モデル１６の物理量（破壊特性）が、良好か否かが判断される。

工程Ｓ６において、高分子材料モデル１６の物理量が良好であると判断された場合（工程Ｓ６で、「Ｙ」）、高分子材料モデル１６に基づいて、高分子材料が製造される（工程Ｓ７）。他方、工程Ｓ６において、高分子材料モデル１６の物理量が良好でないと判断された場合は（工程Ｓ６で、「Ｎ」）、高分子材料の諸条件が変更されて（工程Ｓ８）、工程Ｓ１〜工程Ｓ６が再度実施される。

このように、高分子材料モデル１６の物理量が良好になるまで、高分子材料モデル１６の諸条件が変更されるため、良好な破壊特性を有する高分子材料が、効率よく設計されうる。

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

図３及び図４に示した処理手順に従って、フィラーとマトリックスとを含む高分子材料の数値計算用の高分子材料モデルが設定され、高分子材料の破壊特性が計算された（実施例及び比較例）。

実施例では、図９、図１３及び図１５に示した処理手順に従って、フィラーとマトリックスとの間に斥力系ポテンシャルが作用するとの条件の下、高分子材料モデルの応力と歪との関係を計算する第１シミュレーション工程と、応力と歪との関係に基づいて、応力が降伏するまで高分子材料モデルに加えられた外部エネルギーを計算する工程とが実施された。外部エネルギーＷ１は、上記式（３）に基づいて定義された。

さらに、実施例では、外部エネルギーＷ１に基づいて、フィラーとマトリックスとの間の引力系ポテンシャルが決定され、フィラーモデルとマトリックスモデルとの間に定義された。この引力系ポテンシャルの決定工程では、上記式（３）の外部エネルギーＷ１と、引力系ポテンシャルに基づいて定義された上記式（４）の剥離エネルギーＷ２とが等しくなるとの条件の下、ＬＪポテンシャルの強度εが決定された。そして、フィラーとマトリックスとの間に引力系ポテンシャルが作用するとの条件の下、高分子材料モデルの破壊特性が計算された。

比較例では、上記特許文献１と同様に、フィラーとマトリックスとの間の斥力系ポテンシャルが設定された。その他のポテンシャルについては、実施例と同様に定義された。そして、比較例の高分子材料モデルの破壊特性が計算された。共通仕様は次のとおり、その他のパラメータは、上記特許文献１のとおりである。
マトリックスモデル（分子鎖）：
鎖長：１０００
本数：１８９
数密度：０．９
カップリング剤モデルの個数：１０
フィラー粒子モデルとマトリックスモデルとの間を連結する架橋剤モデル：１５６４個
フィラーの半径：１０．２σ
空間（セル）のサイズ：
ｘ軸方向の長さＬ：６０σ
ｙ軸方向の長さＬ：６０σ
ｚ軸方向の長さＬ：６０．２σ
伸長速度：歪１．０／６３００τ
ポアソン比：０．０
引力系ポテンシャル：
カットオフ距離ｒ_c：２．５σ
ＬＪポテンシャルの強度ε：上記決定工程で求められた値
斥力系ポテンシャル
カットオフ距離ｒ_c：２^1/6σ
ＬＪポテンシャルの強度ε：１．０（上記論文に基づく）

図１６（ａ）〜（ｅ）は、フィラーとマトリックスとの間に引力系ポテンシャルが作用するとの条件の下で、高分子材料モデルの破壊特性を計算した実施例の結果を示す概念図である。図１７（ａ）〜（ｅ）は、フィラーとマトリックスとの間に斥力系ポテンシャルが作用するとの条件の下で、高分子材料モデルの破壊特性を計算した比較例の結果を示す概念図である。

図１６に示した実施例の高分子材料モデルは、フィラーモデルとマトリックスモデルとの界面部分において、マトリックスモデルがフィラーモデルに吸着する高分子材料の描像を再現することができた。他方、図１７に示した比較例の高分子材料モデルは、フィラーモデルとマトリックスモデルとの界面部分において、マトリックスモデルがフィラーモデルから剥離しており、実際の高分子材料の描像を再現することができなかった。従って、実施例は、比較例に比べて、実際の高分子材料の描像に沿った破壊特性を計算することができた。

さらに、実施例では、上記論文に基づく斥力系ポテンシャルを用いて、引力系ポテンシャルに容易に変換することができた。従って、実施例では、実験等に基づいて引力系ポ テンシャルを求める必要がないため、計算コストを低減することができた。

Ｓ１ 高分子材料モデルを設定する工程
Ｓ２ 第１シミュレーション工程
Ｓ３ 外部エネルギーを計算する工程
Ｓ４ 引力系ポテンシャルを決定する工程
Ｓ５ 第２シミュレーション工程

Claims (2)

1. コンピュータを用いて、フィラーとマトリックスとを含む高分子材料の破壊特性を計算するためのシミュレーション方法であって、
   前記高分子材料に基づいて、数値計算用の高分子材料モデルを前記コンピュータに設定する工程と、
   前記フィラーと前記マトリックスとの間に斥力系ポテンシャルが作用するとの条件の下、前記高分子材料モデルの応力と歪との関係を、少なくとも歪が零の状態から前記応力が降伏するまでの期間で計算する第１シミュレーション工程と、
   前記応力と歪との関係に基づいて、前記応力が降伏するまで前記高分子材料モデルに加えられた外部エネルギーを計算する工程と、
   前記外部エネルギーに基づいて、前記フィラーと前記マトリックスとの間の引力系ポテンシャルを決定する決定工程と、
   前記フィラーと前記マトリックスとの間に前記引力系ポテンシャルが作用するとの条件の下、前記高分子材料モデルの破壊特性を計算する第２シミュレーション工程と、
   を含む高分子材料のシミュレーション方法。
2. 前記決定工程は、前記フィラーと前記マトリックスとの間に前記引力系ポテンシャルが作用するとの条件の下、前記フィラーの表面から前記マトリックスを引き剥がすのに必要なエネルギーを定義する工程と、
   前記外部エネルギーと、前記引き剥がすのに必要なエネルギーとが等しくなるとの条件の下、前記引力系ポテンシャルを決定する工程とを含む請求項１記載の高分子材料のシミュレーション方法。

Images