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

Abstract

【課題】高分子材料の破壊特性を予測する。
【解決手段】コンピュータ１を用いて、高分子材料の破壊特性を予測するためのシミュレーション方法である。このシミュレーション方法では、高分子材料に基づいて、高分子材料モデル１６をコンピュータ１に設定するモデル設定工程Ｓ１と、コンピュータ１が、高分子材料モデル１６に予め定められた条件に基づいて歪みを与え、高分子材料の物理量を計算するシミュレーション工程Ｓ２とが含まれる。シミュレーション工程Ｓ２は、微小時間あたりの歪みが１０^−５以上の速度Ｖ１で、高分子材料モデル１６をアフィン変形に基づいて伸長させる高速伸長工程Ｓ２２と、高速伸長工程Ｓ２２で伸ばされた高分子材料モデル１６の大きさを拘束し、その歪みを一定に保持した状態で物理量を計算する物理量計算工程Ｓ２３とを含む。
【選択図】図３

Description

本発明は、高分子材料の破壊特性を予測するのに役立つ高分子材料のシミュレーション方法に関する。

近年、ゴム等の高分子材料の開発のために、高分子材料の性質を、コンピュータを用いて評価するためのシミュレーション方法（数値計算）が種々提案されている（例えば、下記非特許文献１参照）。この種のシミュレーション方法では、先ず、高分子材料に基づいて、高分子材料モデルが設定される。次に、一軸引張試験に基づいて、高分子材料モデルの一端及び他端が、互いに反対側へ引張されて、高分子材料モデルの伸長が計算される。そして、伸長された高分子材料モデルに基づいて、物理量が計算される。

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 」、The Journal of Chemical Physics、VOL.131、014904、 2009年

高分子材料モデルの伸長は、例えば、実際の一軸引張試験に基づいて、徐々に引っ張られる。このようなシミュレーション方法では、高分子材料モデルが予め定められた歪みを持つまでに、多くの時間を要するという問題があった。また、引っ張られる一端及び他端に歪みが集中し、高分子材料モデルに、局所的な破壊（例えば、空孔の形成）が生じやすいという問題があった。このような局所的な破壊は、実際の高分子材料の破壊特性に反するものと考えられている。

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

本発明は、コンピュータを用いて、高分子材料の破壊特性を微小時間毎に計算するためのシミュレーション方法であって、前記高分子材料に基づいて、高分子材料モデルを前記コンピュータに設定するモデル設定工程と、前記コンピュータが、前記高分子材料モデルに予め定められた条件に基づいて歪みを与え、前記高分子材料の物理量を計算するシミュレーション工程とを含み、前記シミュレーション工程は、微小時間あたりの歪みが１０^−５以上の速度で、前記高分子材料モデルをアフィン変形に基づいて伸長させる高速伸長工程と、前記高速伸長工程で伸ばされた前記高分子材料モデルの大きさを拘束し、その歪みを一定に保持した状態で前記物理量を計算する物理量計算工程とを含むことを特徴とする。

本発明に係る前記高分子材料のシミュレーション方法において、前記微小時間は、５．０×１０^−１５s〜５．０×１０^−１０sであるのが望ましい。

本発明に係る前記高分子材料のシミュレーション方法において、前記モデル設定工程は、前記高分子材料の分子鎖をモデル化した分子鎖モデルを設定する工程と、前記分子鎖モデルを予め定められた空間内に配置する工程とを含み、前記物理量計算工程は、分子動力学計算に基づいて、前記分子鎖モデルの熱運動を計算するのが望ましい。

本発明に係る前記高分子材料のシミュレーション方法において、前記シミュレーション工程は、前記高速伸長工程に先立ち、前記分子鎖モデルの熱運動を、前記物理量計算工程時の熱運動よりも小さくする熱運動抑制工程を含むのが望ましい。

本発明に係る前記高分子材料のシミュレーション方法において、前記分子鎖モデルは、複数の粒子モデルを含み、前記モデル設定工程は、前記分子鎖モデルに、摩擦係数を含む運動方程式を定義する工程を含み、前記熱運動抑制工程は、前記分子鎖モデルに、前記物理量計算工程時の前記摩擦係数よりも大きい前記摩擦係数を定義するのが望ましい。

本発明に係る前記高分子材料のシミュレーション方法において、前記高速伸長工程時の前記摩擦係数は、前記物理量計算工程時の前記摩擦係数の２倍〜１０００倍であるのが望ましい。

本発明に係る前記高分子材料のシミュレーション方法において、前記熱運動抑制工程は、前記高分子材料モデルに、前記物理量計算工程時の温度よりも小さい温度を定義するのが望ましい。

本発明に係る前記高分子材料のシミュレーション方法において、前記高速伸長工程は、前記高分子材料を用いた一軸引張試験時のポアソン比よりも小さい値のポアソン比で、前記高分子材料を伸長させるのが望ましい。

本発明に係る前記高分子材料のシミュレーション方法において、前記高速伸長工程は、前記高分子材料モデルの応力が、前記高分子材料の引張強さの１％〜１００％になる歪を与えるのが望ましい。

本発明の高分子材料のシミュレーション方法は、コンピュータが、高分子材料モデルに予め定められた条件に基づいて歪みを与え、前記高分子材料の物理量を計算するシミュレーション工程を含んでいる。シミュレーション工程は、微小時間あたりの歪みが１０^−５以上の速度で、前記高分子材料モデルをアフィン変形に基づいて伸長させる高速伸長工程が含まれる。このようなアフィン変形により、高分子材料モデルの局所的な破壊が計算されるのを防ぐことができる。さらに、高速伸長工程では、従来のシミュレーション方法よりも大きな速度で高分子材料モデルが伸長されるため、計算時間が短縮されうる。

また、シミュレーション工程では、高速伸長工程で伸ばされた高分子材料モデルの大きさを拘束し、その歪みを一定に保持した状態で物理量を計算する物理量計算工程が含まれる。これにより、物理量計算工程では、高分子材料モデルの大きさ（計算領域サイズ）を拘束した状態で、例えば、分子動力学計算に基づく熱運動が計算されることにより、高分子材料モデル全体を対象に、その歪みに基づいた破壊が計算されうる。従って、本発明では、物理量計算工程において、アフィン変形のような人為的な操作が加わるのを防ぎつつ、従来のシミュレーション方法のような局所的な破壊が生じるのを防ぐことができるため、高分子材料の破壊特性が正確に予測されうる。

本発明のシミュレーション方法を実行するためのコンピュータの斜視図である。 高分子材料の構造式である。 本実施形態のシミュレーション方法の処理手順の一例を示すフローチャートである。 本実施形態のモデル設定工程の処理手順の一例を示すフローチャートである。 本実施形態の高分子材料モデルの概念図である。 図５のＡ部拡大図である。 フィラーモデルの概念図である。 フィラー粒子モデル及び結合鎖モデルの拡大図である。 分子鎖モデルの概念図である。 フィラーモデル及び分子鎖モデルを拡大して示す概念図である。 本実施形態のシミュレーション工程の処理手順の一例を示すフローチャートである。 高速伸長工程を説明する高分子材料モデルの概念図である。 （ａ）は、アフィン変形前の小領域の側面図である。（ｂ）は、アフィン変形後の小領域の側面図である。 高分子材料モデルの空孔を示す側面図である。 破壊前の高分子材料モデルの部分斜視図である。 破壊後の高分子材料モデルの部分斜視図である。 本発明の他の実施形態のシミュレーション工程の処理手順の一例を示すフローチャートである。

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

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

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

本実施形態の高分子材料には、充填剤が含有されている。充填剤としては、例えば、カーボンブラック、シリカ、又は、アルミナ等が含まれる。

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

本実施形態のモデル設定工程Ｓ１では、先ず、予め定められた空間が定義される（工程Ｓ１１）。図５は、本実施形態の高分子材料モデルの概念図である。空間２は、後述する分子動力学計算において、計算対象となる領域である。

本実施形態の空間２は、例えば、互いに向き合う少なくとも一対、本実施形態では３対の平面３、３を有する立方体として定義される。空間２の内部には、後述するフィラーモデル５、及び、分子鎖モデル１１が複数個配置される。また、各平面３には、周期境界条件が定義されている。これにより、空間２では、例えば、一方の平面３ａから出て行った分子鎖モデル１１の一部が、反対側の平面３ｂから入ってくるように計算されうる。従って、一方の平面３ａと、反対側の平面３ｂとが連続している（繋がっている）ものとして扱われる。一対の平面３、３の間隔（即ち、１辺の長さＬ１）については、例えば、５０nm〜１０００nm（分子動力学計算の単位では、７６σ〜１５１５σ）に設定されるのが望ましい。

図６は、図５のＡ部拡大図である。空間２の内部には、例えば、立方体状に区分された複数の小領域４が定義されている。各小領域４には、節点４ｔが設定されている。このような小領域４は、後述する分子動力学計算において、各フィラー粒子モデル６（図７に示す）及び粗視化粒子モデル１３（図１０に示す）の追跡等に用いられる。小領域４の１辺の長さＬ２は、例えば、０．５σ〜１０σ程度に設定されるのが望ましい。このような空間２は、コンピュータ１に記憶される。

次に、モデル設定工程Ｓ１では、充填材をモデル化したフィラーモデルが設定される（工程Ｓ１２）。図７は、フィラーモデル５の概念図である。図８は、フィラー粒子モデル６及び結合鎖モデル７の拡大図である。フィラーモデル５は、複数のフィラー粒子モデル６と、隣接するフィラー粒子モデル６、６間を結合する結合鎖モデル７とが含まれている。

フィラー粒子モデル６は、分子動力学計算において、運動方程式の質点として取り扱われる。即ち、フィラー粒子モデル６には、質量、体積、直径、電荷又は初期座標などのパラメータが定義される。

工程Ｓ１２では、フィラー粒子モデル６が配置されるのに先立ち、走査型透過電子顕微鏡によって取得された高分子材料の三次元画像が、空間２（図５に示す）に重ね合わされる。そして、空間２において、三次元画像の充填剤が配置されている領域に、フィラー粒子モデル６が配置される。これにより、フィラー粒子モデル６は、実際の充填剤の分布に近似させて、空間２内に配置されうる。図８に示されるように、フィラー粒子モデル６は、充填剤が配置されている領域に、例えば、結晶格子状に配置されるのが望ましい。これにより、フィラーモデル５は、フィラー粒子モデル６の動きが強固に拘束されうるため、高い剛性を持つ充填剤の物性に近似させることができる。

次に、工程Ｓ１２では、結合鎖モデル７が定義される。本実施形態の結合鎖モデル７は、ボンド関数に基づいて定義される。即ち、結合鎖モデル７は、例えば、下記式（１）で定義されるポテンシャル（以下、「ＬＪポテンシャルＵ_LJ（ｒ_ij）」ということがある。）と、下記式（２）で定義される結合ポテンシャルＵ_FENEとの和で示されるポテンシャルＰ１で定義される。

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

上記式（１）において、フィラー粒子モデル６、６間の距離ｒ_ijが小さくなると、斥力が作用するＬＪポテンシャルＵ_LJ（ｒ_ij）が大きくなる。他方、上記式（２）において、フィラー粒子モデル６、６間の距離ｒ_ij が大きくなると、引力が作用する結合ポテンシャルＵ_FENEが大きくなる。従って、ポテンシャルＰ１は、フィラー粒子モデル６、６間の距離ｒ_ijを、ＬＪポテンシャルＵ_LJ（ｒ_ij）と結合ポテンシャルＵ_FENEとが互いに釣り合う位置に戻そうとする復元力が定義されうる。

また、上記式（１）では、フィラー粒子モデル６、６間の距離ｒ_ijが小さくなるほど、ＬＪポテンシャルＵ_LJ（ｒ_ij）が無限に大きくなる。他方、上記式（２）では、距離ｒ_ijが伸びきり長Ｒ₀以上となる場合に、結合ポテンシャルＵ_FENEが∞に設定されている。従って、ポテンシャルＰ１では、伸びきり長Ｒ₀以上の距離ｒ_ijが許容されない。

なお、ＬＪポテンシャルＵ_LJ（ｒ_ij）及び結合ポテンシャルＵ_FENEの強度ε、伸びきり長Ｒ₀、粒子の直径σ及びカットオフ距離ｒ_cについては、適宜設定されうる。これらの定数は、例えば、論文１（ 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）に基づいて、下記のように設定されるのが望ましい。また、ばね定数ｋは、フィラーモデル５の剛性を決定するパラメータである。ばね定数ｋは、例えば２０〜３０００程度に設定されるのが望ましい。このような結合鎖モデル７が定義されることにより、フィラーモデル５の剛性を高めることができる。フィラーモデル５は、コンピュータ１に記憶される。
強度ε：１．０
伸びきり長Ｒ₀：１．５
距離σ：１．０
カットオフ距離ｒ_c：2^1/6σ

次に、本実施形態のモデル設定工程Ｓ１では、高分子材料の分子鎖をモデル化した分子鎖モデルが設定される（工程Ｓ１３）。図９は、分子鎖モデル１１の概念図である。本実施形態の分子鎖モデル１１は、粗視化モデルとして構成される。即ち、分子鎖モデル１１は、図２及び図９に示されるように、分子鎖のモノマー又はモノマーの一部分をなす構造単位１２が、直径を持った球で表現される粒子モデル（以下、「粗視化粒子モデル」ということがある。）１３に置換される。例えば、分子鎖がポリブタジエンである場合には、１．５５個分のモノマーを構造単位１２として、１個の粗視化粒子モデル１３に置換される。これにより、分子鎖モデル１１には、複数個（例えば、１０個〜５０００個）の粗視化粒子モデル１３が設定される。

なお、１．５５個分のモノマーを構造単位１２としたのは、上記論文１と、論文２（L,J.Fetters ,D.J.Lohse and R.H.Colby 著、「Chain Dimension and Entanglement Spacings 」Physical Properties of Polymers Handbook Second Edition P448」）との記載に基づいて求められうる。また、分子鎖がポリブタジエン以外の場合でも、例えば、論文１及び論文２に基づいて、構造単位１２を設定することができる。

粗視化粒子モデル１３は、分子動力学計算において、運動方程式の質点として取り扱われる。即ち、粗視化粒子モデル１３には、例えば、質量、体積、直径又は電荷などのパラメータが定義される。

分子鎖モデル１１には、隣接する粗視化粒子モデル１３、１３間を結合する結合鎖モデル１４が含まれている。図１０は、フィラーモデル５及び分子鎖モデル１１を拡大して示す概念図である。

結合鎖モデル１４は、ボンド関数に基づいて定義される。即ち、結合鎖モデル１４は、例えば、上記式（１）のＬＪポテンシャルＵ_LJ（ｒ_ij）と、上記式（２）の結合ポテンシャルＵ_FENEとの和で示されるポテンシャルＰ２で設定される。本実施形態では、上記論文１に基づいて、次の値が設定される。
強度ε：１．０
伸びきり長Ｒ₀：１．５
距離σ：１．０
カットオフ距離ｒ_c：2^1/6σ
ばね定数ｋ：３０

このような結合鎖モデル１４により、粗視化粒子モデル１３を伸縮自在に拘束された直鎖状の分子鎖モデル１１を設定することができる。分子鎖モデル１１は、コンピュータ１に記憶される。

図５に示されるように、分子鎖モデル１１は、空間２内に配置される。本実施形態では、分子鎖モデル１１が配置されるのに先立ち、高分子材料の三次元画像が、空間２に重ね合わされる。そして、空間２において、三次元画像の分子鎖が配置されている領域に、分子鎖モデル１１が配置される。これにより、実際の分子鎖の分布に近似させて、分子鎖モデル１１が、空間２内に配置されうる。なお、分子鎖モデル１１（粗視化粒子モデル１３）の個数については、空間２の大きさや、小領域４（図６に示す）の個数等に基づいて、適宜設定されうる。

次に、本実施形態のモデル設定工程Ｓ１では、図１０に示されるように、隣接する分子鎖モデル１１、１１の粗視化粒子モデル１３、１３間に、ポテンシャルＰ３が定義される（工程Ｓ１４）。ポテンシャルＰ３は、上記式（１）のＬＪポテンシャルＵ_LJ（ｒ_ij）によって定義される。なお、ポテンシャルＰ３の各定数及び各変数としては、例えば、上記論文１に基づいて、適宜設定されうる。

次に、本実施形態のモデル設定工程Ｓ１では、隣接するフィラーモデル５のフィラー粒子モデル６、６間、及び、フィラー粒子モデル６と粗視化粒子モデル１３との間に、ポテンシャルＰ４が定義される（工程Ｓ１５）。ポテンシャルＰ４は、上記式（１）のＬＪポテンシャルＵ_LJ（ｒ_ij）によって定義される。また、ポテンシャルＰ４の各定数及び各変数の値としては、例えば、上記論文１に基づいて、適宜設定されうる。

次に、本実施形態のモデル設定工程Ｓ１は、分子鎖モデル１１に、摩擦係数を含んだ運動方程式が定義される（工程Ｓ１６）。本実施形態の運動方程式としては、上記論文１に基づいて、下記式（３）が採用されている。

上記式（３）では、ポテンシャルの和Ｕ_ijから、摩擦力Γｒ_i'を減じ、さらに、ガウスホワイトノイズＷ_i(t)が加えられることにより、各分子鎖モデル１１の加速度ｒ''が定義される。摩擦係数Γは、各粗視化粒子モデル１３について、その周囲の環境から及ぼされる平均的な摩擦力を表すために設定される。また、摩擦係数Γが大きくなるほど、分子鎖モデル１１の加速度ｒ''が小さくなる。なお、各定数及び変数については、高分子材料の種類等に応じて、適宜設定されうる。例えば、本実施形態の摩擦係数Γには、例えば、０．５〜５０００程度が設定されうる。

次に、本実施形態のモデル設定工程Ｓ１は、フィラーモデル５に、摩擦係数を含んだ運動方程式が定義される（工程Ｓ１７）。本実施形態の運動方程式としては、上記式（３）が採用されている。これにより、フィラーモデル５の加速度ｒ''が定義される。

このように、本実施形態のモデル設定工程Ｓ１では、図５に示されるように、フィラーモデル５及び分子鎖モデル１１が空間２に配置され、さらに、各ポテンシャルＰ１〜Ｐ４、及び、運動方程式が定義されることにより、高分子材料モデル１６が設定されうる。このような高分子材料モデル１６は、コンピュータ１に記憶される。

次に、本実施形態のシミュレーション方法では、コンピュータ１が、高分子材料モデル１６に、予め定められた条件に基づいて歪みを与え、高分子材料の物理量を計算する（シミュレーション工程Ｓ２）。図１１は、本実施形態のシミュレーション工程Ｓ２の処理手順の一例を示すフローチャートである。

本実施形態のシミュレーション工程Ｓ２では、先ず、図５に示されるように、分子動力学計算に基づいて、高分子材料モデル１６の構造緩和が計算される（工程Ｓ２１）。本実施形態の工程Ｓ２１では、先ず、フィラーモデル５の位置（座標）を固定して、分子鎖モデル１１のみを対象に、構造緩和が計算される。次に、フィラーモデル５の固定を解除して、フィラーモデル５及び分子鎖モデル１１を対象に、構造緩和が計算される。これにより、高分子材料の充填剤の領域から、フィラーモデル５が大きく位置ずれするのを防ぎつつ、高分子材料モデル１６の構造緩和が計算されうる。

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

工程Ｓ２１では、分子鎖モデル１１の初期配置が、十分に緩和されるまで、シミュレーションの単位時間毎に、高分子材料モデル１６の構造緩和が計算される。これにより、本実施形態では、分子鎖モデル１１の平衡状態（構造が緩和した状態）が、確実に計算されうる。

次に、本実施形態のシミュレーション工程Ｓ２は、構造緩和された高分子材料モデル１６が伸長される（高速伸長工程Ｓ２２）。図１２は、高速伸長工程Ｓ２２を説明する高分子材料モデルの概念図である。高速伸長工程Ｓ２２では、高分子材料モデル１６の一端１６ａ及び他端１６ｂが互いに離間するように、高分子材料モデル１６の伸長が計算される。本例では、ｚ軸方向において、高分子材料モデル１６の一端１６ａ及び他端１６ｂを離間させて、高分子材料モデル１６の伸長が計算される。

本実施形態の高速伸長工程Ｓ２２では、微小時間あたりの歪みが１０^−５以上の速度Ｖ１で、高分子材料モデル１６をアフィン変形に基づいて伸長させている。高分子材料モデル１６の伸長は、上記式（３）の運動方程式が適用された分子動力学計算に基づいて計算されている。本実施形態において、微小時間とは、分子動力学計算の単位ステップ（ＭＤステップ）で計算される時間である。

従来の分子動力学計算では、本実施形態とは異なり、純粋にランジュバン方程式の熱運動に基づいた変形（以下、単に「非アフィン変形」ということがある。）が計算されている。このような従来の変形計算では、伸長される方向（ｚ軸方向）において、高分子材料モデル１６が均一に伸長されない。このため、後述する物理量計算工程Ｓ２３において、一端１６ａ及び他端１６ｂに歪みが集中し、高分子材料モデル１６に局所的な破壊（ボイドの形成）が生じやすい。このような局所的な破壊は、実際の高分子材料の破壊特性に反するものと考えられている。

他方、本実施形態の高速伸長工程Ｓ２２では、上述したように、アフィン変形に基づいて、高分子材料モデル１６が伸長されている。図１３（ａ）は、アフィン変形前の小領域の側面図である。（ｂ）は、アフィン変形後の小領域の側面図である。このようなアフィン変形では、フィラーモデル５のフィラー粒子モデル６（図１０に示す）の座標、及び、分子鎖モデル１１の粗視化粒子モデル１３の座標が強制的に線形変換され、高分子材料モデル１６が相似的に変形されている。

本実施形態では、空間２がｚ軸方向に伸長され、かつ、フィラーモデル５のフィラー粒子モデル６の座標、及び、分子鎖モデル１１の粗視化粒子モデル１３の座標が、ｚ軸方向に強制的に線形変換されている。これにより、高分子材料モデル１６の伸長が計算される。

このように、アフィン変形では、伸長される方向（ｚ軸方向）において、高分子材料モデル１６が均一に伸長されうる。従って、後述する物理量計算工程Ｓ２３において、一端１６ａ及び他端１６ｂに歪みが集中するのを抑制できるため、高分子材料モデル１６に局所的な破壊（ボイドの形成）が生じるのを防ぎうる。

また、高速伸長工程Ｓ２２では、分子鎖モデル１１及びフィラーモデル５の熱運動（分子運動）が計算される。このような熱運動は、本来、高分子材料モデル１６に、微小なボイド（空孔）が形成される要因となる。しかしながら、本実施形態の速度Ｖ１（微小時間あたりの歪みが１０^−５以上）は、分子鎖モデル１１及びフィラーモデル５の熱運動（分子運動）の速度よりも大きい。これにより、高分子材料モデル１６に、微小なボイド（空孔）が形成されるのを効果的に防ぎうる。従って、高分子材料モデル１６の運動に人為的な操作（即ち、アフィン変形）が適用される高速伸長工程Ｓ２２において、高分子材料モデル１６の破壊が生じるのを抑制できる。これにより、人為的な操作が適用されない後述する物理量計算工程Ｓ２３においてのみ、高分子材料モデル１６に破壊を生じさせることができる。

上述のとおり、高速伸長工程Ｓ２２では、アフィン変形による運動よりも弱いものの、熱運動が計算される。図８に示した各フィラー粒子モデル６に定義された結合鎖モデル７の個数は、図９に示した各粗視化粒子モデル１３に定義された結合鎖モデル１４の個数よりも大きい。結合鎖モデルの個数が大きいほど、結合力が高いため、熱運動が強くなる。従って、フィラーモデル５は、分子鎖モデル１１に比べて、熱運動が強くなり、アフィン変形によって大きくなったフィラー粒子モデル６、６間の距離を、迅速に小さくすることができる。これにより、後述する物理量計算工程Ｓ２３において、フィラー粒子モデル６、６間に大きな引力が作用して、計算が破綻（計算落ち）するのを防ぐことができる。

また、本実施形態では、分子鎖モデル１１及びフィラーモデル５の熱運動（分子運動）よりも大きな速度Ｖ１に基づいて、高分子材料モデル１６の伸長が計算されるため、例えば、徐々に引っ張られていた従来のシミュレーション方法に比べて、計算時間が短縮されうる。

高分子材料モデル１６を伸長させる速度Ｖ１については、微小時間あたりの歪みが１０^−５以上であれば、適宜設定されうる。速度Ｖ１が大きいほど、各分子鎖モデル１１の熱運動が抑制されつつ、計算時間を短縮することができる。従って、速度Ｖ１は、微小時間あたりの歪みが、好ましくは１．０４×１０^−４（音速の１倍）以上、さらに好ましくは１．０４×１０^−２（音速の１００倍）以上である。逆に、速度Ｖ１が大きすぎる場合は、単位ステップあたりの粒子の動きが大きくなり、計算が破綻するおそれがある。従って、速度Ｖ１は、微小時間あたりの歪みが、好ましくは１．０４（音速の１００００倍）以下、さらに好ましくは０．１０４（音速の１０００倍）以下である。

また、微小時間については、適宜設定されうる。なお、微小時間が大きいと、分子の熱振動を再現することができなくなり、計算が破綻するおそれがある。逆に、微小時間が小さいと、計算時間が増大するおそれがある。このため、微小時間は、好ましくは０．００３τ以上であり、また、好ましくは０．０１２τ以下である。なお、１τのＳＩ単位系への変換については、想定しているポリマーあるいは着目する物理量によっても異なる。一般的には、５．０×１０^−１５s〜５．０×１０^−１０sである。

高分子材料モデル１６に与えられる歪みについては、適宜定義されうる。なお、歪みが小さいと、後述する物理量計算工程Ｓ２３において、高分子材料モデル１６に十分な破壊が生じないおそれがある。逆に、歪みが大きいと、各分子鎖モデル１１に大きな力が作用し、計算落ちが生じるおそれがある。このため、本実施形態では、高分子材料モデル１６の応力が、高分子材料の引張強さの１％〜１００％になる歪みが与えられるのが望ましい。なお、高分子材料の引張強さは、実際の伸長試験に基づいて予め測定された値が設定されるのが望ましい。また、高分子材料の引張強さは、例えば、高速伸長工程Ｓ２２において取得された応力−歪み曲線に基づいて決定されてもよい。

高速伸長工程Ｓ２２において、伸長後の高分子材料モデル１６のポアソン比は、高分子材料を用いた一軸引張試験時のポアソン比よりも小さいのが望ましい。これにより、小さな歪みでも、高分子材料モデル１６の体積が増大されうる。この高分子材料モデル１６の体積の増大により、ボイドを短時間で生じさせることができ、計算時間が短縮されうる。なお、本実施形態では、ｚ軸方向の長さＬｚのみが拡大するように、高分子材料モデル１６の変形が計算されるため、ポアソン比がゼロである。

次に、本実施形態のシミュレーション工程Ｓ２では、高速伸長工程Ｓ２２で伸ばされた高分子材料モデル１６の大きさが拘束され、その歪みを一定に保持した状態で物理量が計算される（物理量計算工程Ｓ２３）。本実施形態において、「高分子材料モデル１６の大きさが拘束される」とは、ｚ軸方向に伸長された空間２の大きさが拘束されることを示している。図１２に示されるように、高分子材料モデル１６の大きさが拘束されると、高分子材料モデル１６に与えられた歪みが、全ての分子鎖モデル１１及び全てのフィラーモデル５に均一に伝達される。これらの歪み、各ポテンシャルＰ１〜Ｐ４（図８及び図１０に示す）及び、運動方程式に基づいて、分子鎖モデル１１及びフィラーモデル５の熱運動が計算される。

図１４は、高分子材料モデル１６の空孔を示す側面図である。上記のような熱運動により、粗視化粒子モデル１３とフィラー粒子モデル６との間、及び、隣接する分子鎖モデル１１の粗視化粒子モデル１３、１３間が接近又は離間する。これにより、空間２には、粗視化粒子モデル１３及びフィラー粒子モデル６が存在しない小領域４が形成される。このような小領域４は、高分子材料モデル１６に形成された空孔（ボイド）２１として定義される。このような空孔２１が計算されることにより、高分子材料モデル１６の破壊が計算される。

図１５は、破壊前の高分子材料モデル１６の部分斜視図である。図１６は、破壊後の高分子材料モデル１６の部分斜視図である。図１５及び図１６は、図１２に示した高分子材料モデル１６の中央の一部分を拡大したものである。なお、図１５及び図１６では、分子鎖モデル１１が省略して表示されている。図１６に示されるように、本実施形態の高分子材料モデル１６は、フィラーモデル５の周囲を覆っている空孔（ボイド）２１が、ｚ軸方向の広範囲に形成されており、局所的な破壊が防がれうる。

このように、本実施形態の物理量計算工程Ｓ２３では、高分子材料モデル１６の全体に歪みを均一に作用させて、高分子材料モデル１６の破壊が計算されうる。従って、本実施形態の物理量計算工程Ｓ２３では、従来のような局所的な破壊が生じるのを防ぐことができるため、高分子材料の破壊特性が正確に予測されうる。

しかも、本実施形態の高速伸長工程Ｓ２２では、アフィン変形に基づいて、高分子材料モデル１６の伸長が計算されているため、高分子材料モデル１６の局所的な破壊が効果的に防がれうる。さらに、物理量計算工程Ｓ２３では、高分子材料モデル１６の大きさが拘束されているため、アフィン変形のような人為的操作が加わることがなく、高分子材料モデル１６の破壊が計算されうる。従って、本実施形態では、計算精度が低下するのを防ぎうる。

なお、粗視化粒子モデル１３、１３間のポテンシャルＰ３（図１０に示す）において、引力を含むように、カットオフ距離ｒ_cが設定されるのが望ましい。これにより、高分子材料モデル１６において、空孔（ボイド）２１を確実に生じさせることができる。なお、ポテンシャルＰ３のカットオフ距離ｒ_cは、例えば２．５σ程度に設定されるのが望ましい。

物理量計算工程Ｓ２３では、高分子材料モデル１６の物理量として、例えば、高分子材料モデル１６の応力や、空孔（ボイド）２１の偏りを示す指標が計算されるのが望ましい。本実施形態では、空孔（ボイド）２１の偏りを示す指標としては、例えば、下記式（４）に示すシャノンエントロピーＳが用いられている。

ここで、各定数及び変数は、次のとおりである。
Ｎ_div：小領域４の総数
ｉ：小領域の添字（ｉ＝１〜Ｎ_div）
Ｌｚ：引張後の空間２のＺ軸方向の長さ
Ｃ_ｉ：（ｉ−１）Ｌｚ／Ｎ_div≦ｚ≦ｉ・Ｌｚ／Ｎ_divの範囲において、空孔２１が存在する場合に、「１」を返し、空孔２１が存在しない場合に、「０」を返す関数である。

上記式（４）で示されるシャノンエントロピーＳは、各小領域４に空孔２１が万遍なく存在する場合、最大値ｌｏｇ（１／Ｎ_div）を取る。他方、空孔２１の偏りが大きいほど、シャノンエントロピーＳが小さくなる。従って、このようなシャノンエントロピーＳが用いられることにより、空孔２１の偏りが容易に判断されうる。

次に、本実施形態のシミュレーション方法では、高分子材料モデル１６の物理量（破壊特性）が、良好か否かが判断される（工程Ｓ３）。本実施形態では、例えば、高分子材料モデル１６の応力、空孔（ボイド）２１の総量、又は、空孔（ボイド）２１の偏り等に基づいて、高分子材料モデル１６の物理量（破壊特性）が、良好か否かが判断される。工程Ｓ３では、高分子材料モデル１６の物理量が良好であると判断された場合（工程Ｓ３で、「Ｙ」）、高分子材料モデル１６に基づいて、高分子材料が製造される（工程Ｓ４）。他方、工程Ｓ３では、高分子材料モデル１６の物理量が良好でないと判断された場合は（工程Ｓ３で、「Ｎ」）、高分子材料モデル１６の諸条件が変更されて（工程Ｓ５）、工程Ｓ１〜工程Ｓ３が再度実施される。

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

本実施形態のシミュレーション方法では、高速伸長工程Ｓ２２において、微小時間あたりの歪みが１０^−５以上の速度Ｖ１で高分子材料モデル１６が伸長されることにより、分子鎖モデル１１の熱運動の計算が抑制される態様が例示されたが、これに限定されるわけではない。例えば、高速伸長工程Ｓ２２に先立ち、分子鎖モデル１１の熱運動を、物理量計算工程Ｓ２３時の熱運動よりも小さくする熱運動抑制工程Ｓ２４が含まれるのが望ましい。図１７は、本発明の他の実施形態のシミュレーション工程Ｓ２の処理手順の一例を示すフローチャートである。なお、この実施形態において、前実施形態と同一の構成、及び、同一の処理が実施される工程については、同一の符号を付して、説明が省略される場合がある。

この実施形態の熱運動抑制工程Ｓ２４は、分子鎖モデル１１（図１０に示す）に、物理量計算工程Ｓ２３時の摩擦係数よりも大きい摩擦係数が定義される。この実施形態では、上記式（３）の運動方程式において、物理量計算工程Ｓ２３時よりも大きな摩擦係数Γが定義されうる。

このような大きな摩擦係数Γが定義されることにより、高速伸長工程Ｓ２２において、分子鎖モデル１１の粗視化粒子モデル１３の動きが強固に拘束される。これにより、この実施形態の高速伸長工程Ｓ２２では、分子鎖モデル１１の熱運動が、物理量計算工程Ｓ２３時の熱運動よりも小に設定されうる。これにより、高速伸長工程Ｓ２２では、高分子材料モデル１６の一端１６ａ及び他端１６ｂにおいて計算されがちな局所的な破壊が、より効果的に防がれうる。

また、この実施形態のシミュレーション工程Ｓ２では、物理量計算工程Ｓ２３に先立ち、分子鎖モデル１１の熱運動が元に戻される（工程Ｓ２５）。これにより、物理量計算工程Ｓ２３では、本来の摩擦係数Γに基づいて、分子鎖モデル１１の熱運動が計算されうる。従って、高分子材料モデル１６の物理量が、精度良く計算されうる。

高速伸長工程Ｓ２２時の摩擦係数Γについては、適宜設定されうる。高速伸長工程Ｓ２２時の摩擦係数Γが小さいと、分子鎖モデル１１の熱運動を十分に小さくできないおそれがある。このため、高速伸長工程Ｓ２２の摩擦係数Γは、好ましくは、物理量計算工程Ｓ２３時の摩擦係数Γの２倍以上であり、より好ましくは１０倍以上である。他方、高速伸長工程Ｓ２２時の摩擦係数Γが大きいと、アフィン変形によるフィラーモデル５の自発的な熱運動が阻害されるため、物理量計算工程Ｓ２３において、フィラー粒子モデル６、６間の引力が大きくなり、計算落ちが発生するおそれがある。このため、高速伸長工程Ｓ２２の摩擦係数Γは、好ましくは、物理量計算工程Ｓ２３時の摩擦係数Γの１０００倍以下であり、より好ましくは２００倍以下である。

前実施形態の熱運動抑制工程Ｓ２４では、分子鎖モデル１１に、物理量計算工程Ｓ２３時の摩擦係数Γよりも大きい摩擦係数Γが定義される場合が例示されたが、これに限定されるわけではない。例えば、高分子材料モデル１６に、物理量計算工程Ｓ２３時の温度よりも小さい温度が定義されてもよい。

このような小さな温度が定義されることにより、高速伸長工程Ｓ２２において、分子鎖モデル１１の熱運動が、物理量計算工程Ｓ２３時の熱運動よりも小さくなる。これにより、高速伸長工程Ｓ２２では、高分子材料モデル１６の一端１６ａ及び他端１６ｂにおいて計算されがちな局所的な破壊が、より効果的に防がれうる。

高速伸長工程Ｓ２２時の温度については、適宜設定されうる。高速伸長工程Ｓ２２時の温度が大きいと、分子鎖モデル１１の熱運動を十分に小さくできないおそれがある。このため、高速伸長工程Ｓ２２の温度は、好ましくは、好ましくは、物理量計算工程Ｓ２３時の温度の０．５倍以下であり、より好ましくは０．１倍以下である。

これまでの熱運動抑制工程Ｓ２４では、物理量計算工程Ｓ２３時の摩擦係数Γよりも大きい摩擦係数Γや、物理量計算工程Ｓ２３時の温度よりも小さい温度が、高分子材料モデル１６にそれぞれ定義される場合が例示されたが、これに限定されるわけではない。例えば、物理量計算工程Ｓ２３時の摩擦係数Γよりも大きい摩擦係数Γと、物理量計算工程Ｓ２３時の温度よりも小さい温度との双方が、高分子材料モデル１６に定義されてもよい。これにより、高速伸長工程Ｓ２２において、分子鎖モデル１１の熱運動を、物理量計算工程Ｓ２３時の熱運動よりも効果的に小さくすることができる。

図５に示されるように、本実施形態の分子鎖モデル１１は、粗視化粒子モデル１３と結合鎖モデル１４とを含む粗視化モデルとして構成される場合が例示されたが、これに限定されるものではない。分子鎖モデル１１は、図２に示した分子鎖のモノマーに基づいて、複数の粒子モデルとボンドモデルとで定義された全原子モデルや、ユナイテッドアトムモデルとして構成されてもよい。このような全原子モデル又はユナイテッドアトムモデルは、粗視化モデルに比べて、分子鎖の挙動が正確に表現されうるため、高分子材料の破壊特性がより正確に予測されうる。また、高分子材料モデル１６は、空間２の内部に、フィラーモデル５及び分子鎖モデル１１が配置されたモデルに限定されるわけではなく、例えば、連続体モデル（一例として、有限要素モデル等）として構成されてもよい。

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

図４に示した処理手順に従って、予め定められた空間に、フィラーモデル及び分子鎖モデルが配置され、高分子材料モデルが設定された（実施例、比較例１、比較例２）。

実施例では、図１１に示した処理手順に従って、高分子材料モデルが、微小時間あたりの歪みが１０^−５以上の速度で伸長される高速伸長工程と、高分子材料モデルの大きさを拘束し、その歪みを一定に保持した状態で物理量を計算する物理量計算工程とが実施された。また、実施例の高速伸長工程では、アフィン変形に基づいて、高分子材料の伸長が計算された。さらに、実施例では、高速伸長工程に先立ち、物理量計算工程時の摩擦係数よりも大きい摩擦係数が定義された（熱運動抑制工程）。

比較例１及び比較例２では、微小時間あたりの歪みが１０^−５未満の速度に基づいて、高分子材料モデルが徐々に引っ張られ、高分子材料モデルの物理量が逐次計算された。また、比較例１では、非アフィン変形に基づいて、高分子材料モデルの伸長が計算された。他方、比較例２では、アフィン変形に基づいて、高分子材料モデルの伸長が計算された。

実施例、比較例１及び比較例２では、上記式（４）に基づいて、空孔（ボイド）の偏りを示す指標であるシャノンエントロピーＳが計算された。シャノンエントロピーＳが大きいほど、空孔の偏りが小さく良好である。なお、ポテンシャル等の各数値は、明細書中の記載の通りであり、その他の共通仕様は次のとおりである。結果を表１に示す。
高分子材料モデル：
分子鎖モデル：
粗視化粒子モデルの個数：１００００００個
粗視化粒子モデルの数密度：０．９σ^−３
分子鎖モデル１個あたりの粗視化粒子モデルの個数：２００個
フィラーモデル：
フィラー粒子モデルの体積分率：０．１８
フィラーモデルの半径：１０σ
フィラーモデル１個あたりのフィラー粒子モデルの個数：５０００個
シミュレーション：
コンピュータのＣＰＵ：intel Core i7-4770（４コア）
ポアソン比：０
高分子材料モデルのひずみ：０．３（高分子材料の引張強さの１００％）
実施例：
高速伸長工程：
速度Ｖ１：微小時間あたりの歪み３．１２×１０^−４
（音速の３倍［ひずみ速度：０．０５３（１／τ）］）
微小時間：０．００６τ（６×１０^−１３s）
高速伸長工程のシミュレーション時間：５．６６τ
物理量計算工程のシミュレーション時間：１１８１．１τ
高速伸長工程の摩擦係数：物理量計算工程時の摩擦係数の１０倍
比較例１、比較例２：
高速伸長工程：
速度Ｖ１：微小時間あたりの歪み１．４７×１０^−６
（音速の０．０１５倍［ひずみ速度：０．０００２５（１／τ）］）
微小時間：０．００６τ（６×１０^−１３s）
引張時間：１２００τ

テストの結果、実施例は、比較例１に比べて、シャノンエントロピーが大であり、空孔の偏りが小さかった。従って、実施例は、比較例１に比べて、局所的な破壊が生じるのを防ぐことができ、高分子材料の破壊特性を正確に予測しうることが確認できた。

実施例と、比較例２とは、シャノンエントロピーが同一であった。しかしながら、比較例２では、物理量を計算する工程において、分子運動に対して、アフィン変形による人為的な操作が常に加わった。他方、実施例では、物理量計算工程において、高分子材料モデルの大きさが拘束されているため、分子運動に対して、アフィン変形による人為的な操作が加わらなかった。さらに、実施例では、高速伸長工程において、アフィン変形に基づいて、高分子材料モデルを伸長させているため、高分子材料モデルの一端及び他端において計算されがちな局所的な破壊が効果的に防がれた。従って、実施例は、比較例２に比べて、高分子材料の破壊特性を正確に予測できた。

実施例では、高分子材料モデルの大きさを拘束して、その歪みを一定に保持した状態で物理量が計算されている。ここで、歪みは、物理量計算工程において、破壊が比較的急速に進行する程度に大きく設定されている。このため、実施例では、微小時間あたりの歪みが１０^−５未満の速度で徐々に変形させながら物理量が計算される比較例１及び比較例２に比べて、物理量計算工程のシミュレーション時間の一部を省略することができる。従って、実施例は、比較例１及び比較例２に比べて、計算時間を短縮しうる。

１ コンピュータ
１６ 高分子材料モデル

Claims (9)

1. コンピュータを用いて、高分子材料の破壊特性を微小時間毎に計算するためのシミュレーション方法であって、
   前記高分子材料に基づいて、高分子材料モデルを前記コンピュータに設定するモデル設定工程と、
   前記コンピュータが、前記高分子材料モデルに予め定められた条件に基づいて歪みを与え、前記高分子材料の物理量を計算するシミュレーション工程とを含み、
   前記シミュレーション工程は、微小時間あたりの歪みが１０^−５以上の速度で、前記高分子材料モデルをアフィン変形に基づいて伸長させる高速伸長工程と、
   前記高速伸長工程で伸ばされた前記高分子材料モデルの大きさを拘束し、その歪みを一定に保持した状態で前記物理量を計算する物理量計算工程とを含むことを特徴とする高分子材料のシミュレーション方法。
2. 前記微小時間は、５．０×１０^−１５s〜５．０×１０^−１０sである請求項１記載の高分子材料のシミュレーション方法。
3. 前記モデル設定工程は、前記高分子材料の分子鎖をモデル化した分子鎖モデルを設定する工程と、
   前記分子鎖モデルを予め定められた空間内に配置する工程とを含み、
   前記物理量計算工程は、分子動力学計算に基づいて、前記分子鎖モデルの熱運動を計算する請求項１又は２記載の高分子材料のシミュレーション方法。
4. 前記シミュレーション工程は、前記高速伸長工程に先立ち、前記分子鎖モデルの熱運動を、前記物理量計算工程時の熱運動よりも小さくする熱運動抑制工程を含む請求項３記載の高分子材料のシミュレーション方法。
5. 前記分子鎖モデルは、複数の粒子モデルを含み、
   前記モデル設定工程は、前記分子鎖モデルに、摩擦係数を含む運動方程式を定義する工程を含み、
   前記熱運動抑制工程は、前記分子鎖モデルに、前記物理量計算工程時の前記摩擦係数よりも大きい前記摩擦係数を定義する請求項４記載の高分子材料のシミュレーション方法。
6. 前記高速伸長工程時の前記摩擦係数は、前記物理量計算工程時の前記摩擦係数の２倍〜１０００倍である請求項５記載の高分子材料のシミュレーション方法。
7. 前記熱運動抑制工程は、前記高分子材料モデルに、前記物理量計算工程時の温度よりも小さい温度を定義する請求項４乃至６のいずれかに記載の高分子材料のシミュレーション方法。
8. 前記高速伸長工程は、前記高分子材料を用いた一軸引張試験時のポアソン比よりも小さい値のポアソン比で、前記高分子材料を伸長させる請求項１乃至７のいずれかに記載の高分子材料のシミュレーション方法。
9. 前記高速伸長工程は、前記高分子材料モデルの応力が、前記高分子材料の引張強さの１％〜１００％になる歪を与える請求項１乃至８のいずれかに記載の高分子材料のシミュレーション方法。