以下、本発明の実施の一形態が図面に基づき説明される。
本実施形態の高分子材料のシミュレーション方法(以下、単に「シミュレーション方法」ということがある)は、コンピュータを用いて、高分子材料の破壊特性を予測するための方法である。
図1は、本発明のシミュレーション方法を実行するためのコンピュータ1の斜視図である。コンピュータ1は、本体1a、キーボード1b、マウス1c及びディスプレイ装置1dを含んでいる。この本体1aには、例えば、演算処理装置(CPU)、ROM、作業用メモリ、磁気ディスクなどの記憶装置、及びディスクドライブ装置1a1、1a2が設けられる。また、記憶装置には、本実施形態のシミュレーション方法を実行するための処理手順(プログラム)が予め記憶される。
高分子材料としては、例えば、天然ゴム、合成ゴム、又は、樹脂等が含まれる。図2は、高分子材料の構造式である。本実施形態の高分子材料としては、cis-1,4ポリブタジエン(以下、単に「ポリブタジエン」ということがある)が例示される。このポリブタジエンを構成する高分子鎖は、メチレン基(−CH2−)とメチン基(−CH−)とからなるモノマー{−[CH2−CH=CH−CH2]−}が、重合度nで連結されて構成されている。なお、高分子材料には、ポリブタジエン以外の高分子材料が用いられてもよい。
本実施形態の高分子材料には、充填剤が含有されている。充填剤としては、例えば、カーボンブラック、シリカ、又は、アルミナ等が含まれる。
図3は、本実施形態のシミュレーション方法の処理手順の一例を示すフローチャートである。本実施形態のシミュレーション方法では、先ず、高分子材料に基づいて、高分子材料モデルがコンピュータ1に設定される(モデル設定工程S1)。図4は、本実施形態のモデル設定工程S1の処理手順の一例を示すフローチャートである。
本実施形態のモデル設定工程S1では、先ず、予め定められた空間が定義される(工程S11)。図5は、本実施形態の高分子材料モデルの概念図である。空間2は、後述する分子動力学計算において、計算対象となる領域である。
本実施形態の空間2は、例えば、互いに向き合う少なくとも一対、本実施形態では3対の平面3、3を有する立方体として定義される。空間2の内部には、後述するフィラーモデル5、及び、分子鎖モデル11が複数個配置される。また、各平面3には、周期境界条件が定義されている。これにより、空間2では、例えば、一方の平面3aから出て行った分子鎖モデル11の一部が、反対側の平面3bから入ってくるように計算されうる。従って、一方の平面3aと、反対側の平面3bとが連続している(繋がっている)ものとして扱われる。一対の平面3、3の間隔(即ち、1辺の長さL1)については、例えば、50nm〜1000nm(分子動力学計算の単位では、76σ〜1515σ)に設定されるのが望ましい。
図6は、図5のA部拡大図である。空間2の内部には、例えば、立方体状に区分された複数の小領域4が定義されている。各小領域4には、節点4tが設定されている。このような小領域4は、後述する分子動力学計算において、各フィラー粒子モデル6(図7に示す)及び粗視化粒子モデル13(図10に示す)の追跡等に用いられる。小領域4の1辺の長さL2は、例えば、0.5σ〜10σ程度に設定されるのが望ましい。このような空間2は、コンピュータ1に記憶される。
次に、モデル設定工程S1では、充填材をモデル化したフィラーモデルが設定される(工程S12)。図7は、フィラーモデル5の概念図である。図8は、フィラー粒子モデル6及び結合鎖モデル7の拡大図である。フィラーモデル5は、複数のフィラー粒子モデル6と、隣接するフィラー粒子モデル6、6間を結合する結合鎖モデル7とが含まれている。
フィラー粒子モデル6は、分子動力学計算において、運動方程式の質点として取り扱われる。即ち、フィラー粒子モデル6には、質量、体積、直径、電荷又は初期座標などのパラメータが定義される。
工程S12では、フィラー粒子モデル6が配置されるのに先立ち、走査型透過電子顕微鏡によって取得された高分子材料の三次元画像が、空間2(図5に示す)に重ね合わされる。そして、空間2において、三次元画像の充填剤が配置されている領域に、フィラー粒子モデル6が配置される。これにより、フィラー粒子モデル6は、実際の充填剤の分布に近似させて、空間2内に配置されうる。図8に示されるように、フィラー粒子モデル6は、充填剤が配置されている領域に、例えば、結晶格子状に配置されるのが望ましい。これにより、フィラーモデル5は、フィラー粒子モデル6の動きが強固に拘束されうるため、高い剛性を持つ充填剤の物性に近似させることができる。
次に、工程S12では、結合鎖モデル7が定義される。本実施形態の結合鎖モデル7は、ボンド関数に基づいて定義される。即ち、結合鎖モデル7は、例えば、下記式(1)で定義されるポテンシャル(以下、「LJポテンシャルULJ(rij)」ということがある。)と、下記式(2)で定義される結合ポテンシャルUFENEとの和で示されるポテンシャルP1で定義される。
ここで、各定数及び変数は、Lennard-Jones及びFENEの各ポテンシャルのパラメータであり、次のとおりである。
rij:粒子間の距離
rc:カットオフ距離
k:粒子間のばね定数
ε:粒子間に定義されるLJポテンシャルの強度
σ:粒子の直径に相当
R0:伸びきり長
なお、距離rij、カットオフ距離rc、及び、伸びきり長R0は、各フィラー粒子モデル6の中心6c、6c間の距離として定義される。
上記式(1)において、フィラー粒子モデル6、6間の距離rijが小さくなると、斥力が作用するLJポテンシャルULJ(rij)が大きくなる。他方、上記式(2)において、フィラー粒子モデル6、6間の距離rij が大きくなると、引力が作用する結合ポテンシャルUFENEが大きくなる。従って、ポテンシャルP1は、フィラー粒子モデル6、6間の距離rijを、LJポテンシャルULJ(rij)と結合ポテンシャルUFENEとが互いに釣り合う位置に戻そうとする復元力が定義されうる。
また、上記式(1)では、フィラー粒子モデル6、6間の距離rijが小さくなるほど、LJポテンシャルULJ(rij)が無限に大きくなる。他方、上記式(2)では、距離rijが伸びきり長R0以上となる場合に、結合ポテンシャルUFENEが∞に設定されている。従って、ポテンシャルP1では、伸びきり長R0以上の距離rijが許容されない。
なお、LJポテンシャルULJ(rij)及び結合ポテンシャルUFENEの強度ε、伸びきり長R0、粒子の直径σ及びカットオフ距離rcについては、適宜設定されうる。これらの定数は、例えば、論文1( 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)に基づいて、下記のように設定されるのが望ましい。また、ばね定数kは、フィラーモデル5の剛性を決定するパラメータである。ばね定数kは、例えば20〜3000程度に設定されるのが望ましい。このような結合鎖モデル7が定義されることにより、フィラーモデル5の剛性を高めることができる。フィラーモデル5は、コンピュータ1に記憶される。
強度ε:1.0
伸びきり長R0:1.5
距離σ:1.0
カットオフ距離rc:21/6σ
次に、本実施形態のモデル設定工程S1では、高分子材料の分子鎖をモデル化した分子鎖モデルが設定される(工程S13)。図9は、分子鎖モデル11の概念図である。本実施形態の分子鎖モデル11は、粗視化モデルとして構成される。即ち、分子鎖モデル11は、図2及び図9に示されるように、分子鎖のモノマー又はモノマーの一部分をなす構造単位12が、直径を持った球で表現される粒子モデル(以下、「粗視化粒子モデル」ということがある。)13に置換される。例えば、分子鎖がポリブタジエンである場合には、1.55個分のモノマーを構造単位12として、1個の粗視化粒子モデル13に置換される。これにより、分子鎖モデル11には、複数個(例えば、10個〜5000個)の粗視化粒子モデル13が設定される。
なお、1.55個分のモノマーを構造単位12としたのは、上記論文1と、論文2(L,J.Fetters ,D.J.Lohse and R.H.Colby 著、「Chain Dimension and Entanglement Spacings 」Physical Properties of Polymers Handbook Second Edition P448」)との記載に基づいて求められうる。また、分子鎖がポリブタジエン以外の場合でも、例えば、論文1及び論文2に基づいて、構造単位12を設定することができる。
粗視化粒子モデル13は、分子動力学計算において、運動方程式の質点として取り扱われる。即ち、粗視化粒子モデル13には、例えば、質量、体積、直径又は電荷などのパラメータが定義される。
分子鎖モデル11には、隣接する粗視化粒子モデル13、13間を結合する結合鎖モデル14が含まれている。図10は、フィラーモデル5及び分子鎖モデル11を拡大して示す概念図である。
結合鎖モデル14は、ボンド関数に基づいて定義される。即ち、結合鎖モデル14は、例えば、上記式(1)のLJポテンシャルULJ(rij)と、上記式(2)の結合ポテンシャルUFENEとの和で示されるポテンシャルP2で設定される。本実施形態では、上記論文1に基づいて、次の値が設定される。
強度ε:1.0
伸びきり長R0:1.5
距離σ:1.0
カットオフ距離rc:21/6σ
ばね定数k:30
このような結合鎖モデル14により、粗視化粒子モデル13を伸縮自在に拘束された直鎖状の分子鎖モデル11を設定することができる。分子鎖モデル11は、コンピュータ1に記憶される。
図5に示されるように、分子鎖モデル11は、空間2内に配置される。本実施形態では、分子鎖モデル11が配置されるのに先立ち、高分子材料の三次元画像が、空間2に重ね合わされる。そして、空間2において、三次元画像の分子鎖が配置されている領域に、分子鎖モデル11が配置される。これにより、実際の分子鎖の分布に近似させて、分子鎖モデル11が、空間2内に配置されうる。なお、分子鎖モデル11(粗視化粒子モデル13)の個数については、空間2の大きさや、小領域4(図6に示す)の個数等に基づいて、適宜設定されうる。
次に、本実施形態のモデル設定工程S1では、図10に示されるように、隣接する分子鎖モデル11、11の粗視化粒子モデル13、13間に、ポテンシャルP3が定義される(工程S14)。ポテンシャルP3は、上記式(1)のLJポテンシャルULJ(rij)によって定義される。なお、ポテンシャルP3の各定数及び各変数としては、例えば、上記論文1に基づいて、適宜設定されうる。
次に、本実施形態のモデル設定工程S1では、隣接するフィラーモデル5のフィラー粒子モデル6、6間、及び、フィラー粒子モデル6と粗視化粒子モデル13との間に、ポテンシャルP4が定義される(工程S15)。ポテンシャルP4は、上記式(1)のLJポテンシャルULJ(rij)によって定義される。また、ポテンシャルP4の各定数及び各変数の値としては、例えば、上記論文1に基づいて、適宜設定されうる。
次に、本実施形態のモデル設定工程S1は、分子鎖モデル11に、摩擦係数を含んだ運動方程式が定義される(工程S16)。本実施形態の運動方程式としては、上記論文1に基づいて、下記式(3)が採用されている。
上記式(3)では、ポテンシャルの和Uijから、摩擦力Γri'を減じ、さらに、ガウスホワイトノイズWi(t)が加えられることにより、各分子鎖モデル11の加速度r''が定義される。摩擦係数Γは、各粗視化粒子モデル13について、その周囲の環境から及ぼされる平均的な摩擦力を表すために設定される。また、摩擦係数Γが大きくなるほど、分子鎖モデル11の加速度r''が小さくなる。なお、各定数及び変数については、高分子材料の種類等に応じて、適宜設定されうる。例えば、本実施形態の摩擦係数Γには、例えば、0.5〜5000程度が設定されうる。
次に、本実施形態のモデル設定工程S1は、フィラーモデル5に、摩擦係数を含んだ運動方程式が定義される(工程S17)。本実施形態の運動方程式としては、上記式(3)が採用されている。これにより、フィラーモデル5の加速度r''が定義される。
このように、本実施形態のモデル設定工程S1では、図5に示されるように、フィラーモデル5及び分子鎖モデル11が空間2に配置され、さらに、各ポテンシャルP1〜P4、及び、運動方程式が定義されることにより、高分子材料モデル16が設定されうる。このような高分子材料モデル16は、コンピュータ1に記憶される。
次に、本実施形態のシミュレーション方法では、コンピュータ1が、高分子材料モデル16に、予め定められた条件に基づいて歪みを与え、高分子材料の物理量を計算する(シミュレーション工程S2)。図11は、本実施形態のシミュレーション工程S2の処理手順の一例を示すフローチャートである。
本実施形態のシミュレーション工程S2では、先ず、図5に示されるように、分子動力学計算に基づいて、高分子材料モデル16の構造緩和が計算される(工程S21)。本実施形態の工程S21では、先ず、フィラーモデル5の位置(座標)を固定して、分子鎖モデル11のみを対象に、構造緩和が計算される。次に、フィラーモデル5の固定を解除して、フィラーモデル5及び分子鎖モデル11を対象に、構造緩和が計算される。これにより、高分子材料の充填剤の領域から、フィラーモデル5が大きく位置ずれするのを防ぎつつ、高分子材料モデル16の構造緩和が計算されうる。
本実施形態の分子動力学計算では、例えば、空間2について所定の時間、分子鎖モデル11又はフィラーモデル5が古典力学に従うものとして、上記式(3)の運動方程式が適用される。そして、各時刻での分子鎖モデル11の動きが、単位時間毎に追跡される。本実施形態では、空間2において、圧力及び温度が一定、又は、体積及び温度が一定に保たれている。これにより、工程S21では、実際の高分子材料の分子運動に近似させて、分子鎖モデル11の初期配置が、精度よく緩和される。このような構造緩和の計算には、例えば、(株)JSOL社製のソフトマテリアル総合シミュレーター(J−OCTA)に含まれているCOGNACが用いられるのが望ましい。
工程S21では、分子鎖モデル11の初期配置が、十分に緩和されるまで、シミュレーションの単位時間毎に、高分子材料モデル16の構造緩和が計算される。これにより、本実施形態では、分子鎖モデル11の平衡状態(構造が緩和した状態)が、確実に計算されうる。
次に、本実施形態のシミュレーション工程S2は、構造緩和された高分子材料モデル16が伸長される(高速伸長工程S22)。図12は、高速伸長工程S22を説明する高分子材料モデルの概念図である。高速伸長工程S22では、高分子材料モデル16の一端16a及び他端16bが互いに離間するように、高分子材料モデル16の伸長が計算される。本例では、z軸方向において、高分子材料モデル16の一端16a及び他端16bを離間させて、高分子材料モデル16の伸長が計算される。
本実施形態の高速伸長工程S22では、微小時間あたりの歪みが10−5以上の速度V1で、高分子材料モデル16をアフィン変形に基づいて伸長させている。高分子材料モデル16の伸長は、上記式(3)の運動方程式が適用された分子動力学計算に基づいて計算されている。本実施形態において、微小時間とは、分子動力学計算の単位ステップ(MDステップ)で計算される時間である。
従来の分子動力学計算では、本実施形態とは異なり、純粋にランジュバン方程式の熱運動に基づいた変形(以下、単に「非アフィン変形」ということがある。)が計算されている。このような従来の変形計算では、伸長される方向(z軸方向)において、高分子材料モデル16が均一に伸長されない。このため、後述する物理量計算工程S23において、一端16a及び他端16bに歪みが集中し、高分子材料モデル16に局所的な破壊(ボイドの形成)が生じやすい。このような局所的な破壊は、実際の高分子材料の破壊特性に反するものと考えられている。
他方、本実施形態の高速伸長工程S22では、上述したように、アフィン変形に基づいて、高分子材料モデル16が伸長されている。図13(a)は、アフィン変形前の小領域の側面図である。(b)は、アフィン変形後の小領域の側面図である。このようなアフィン変形では、フィラーモデル5のフィラー粒子モデル6(図10に示す)の座標、及び、分子鎖モデル11の粗視化粒子モデル13の座標が強制的に線形変換され、高分子材料モデル16が相似的に変形されている。
本実施形態では、空間2がz軸方向に伸長され、かつ、フィラーモデル5のフィラー粒子モデル6の座標、及び、分子鎖モデル11の粗視化粒子モデル13の座標が、z軸方向に強制的に線形変換されている。これにより、高分子材料モデル16の伸長が計算される。
このように、アフィン変形では、伸長される方向(z軸方向)において、高分子材料モデル16が均一に伸長されうる。従って、後述する物理量計算工程S23において、一端16a及び他端16bに歪みが集中するのを抑制できるため、高分子材料モデル16に局所的な破壊(ボイドの形成)が生じるのを防ぎうる。
また、高速伸長工程S22では、分子鎖モデル11及びフィラーモデル5の熱運動(分子運動)が計算される。このような熱運動は、本来、高分子材料モデル16に、微小なボイド(空孔)が形成される要因となる。しかしながら、本実施形態の速度V1(微小時間あたりの歪みが10−5以上)は、分子鎖モデル11及びフィラーモデル5の熱運動(分子運動)の速度よりも大きい。これにより、高分子材料モデル16に、微小なボイド(空孔)が形成されるのを効果的に防ぎうる。従って、高分子材料モデル16の運動に人為的な操作(即ち、アフィン変形)が適用される高速伸長工程S22において、高分子材料モデル16の破壊が生じるのを抑制できる。これにより、人為的な操作が適用されない後述する物理量計算工程S23においてのみ、高分子材料モデル16に破壊を生じさせることができる。
上述のとおり、高速伸長工程S22では、アフィン変形による運動よりも弱いものの、熱運動が計算される。図8に示した各フィラー粒子モデル6に定義された結合鎖モデル7の個数は、図9に示した各粗視化粒子モデル13に定義された結合鎖モデル14の個数よりも大きい。結合鎖モデルの個数が大きいほど、結合力が高いため、熱運動が強くなる。従って、フィラーモデル5は、分子鎖モデル11に比べて、熱運動が強くなり、アフィン変形によって大きくなったフィラー粒子モデル6、6間の距離を、迅速に小さくすることができる。これにより、後述する物理量計算工程S23において、フィラー粒子モデル6、6間に大きな引力が作用して、計算が破綻(計算落ち)するのを防ぐことができる。
また、本実施形態では、分子鎖モデル11及びフィラーモデル5の熱運動(分子運動)よりも大きな速度V1に基づいて、高分子材料モデル16の伸長が計算されるため、例えば、徐々に引っ張られていた従来のシミュレーション方法に比べて、計算時間が短縮されうる。
高分子材料モデル16を伸長させる速度V1については、微小時間あたりの歪みが10−5以上であれば、適宜設定されうる。速度V1が大きいほど、各分子鎖モデル11の熱運動が抑制されつつ、計算時間を短縮することができる。従って、速度V1は、微小時間あたりの歪みが、好ましくは1.04×10−4(音速の1倍)以上、さらに好ましくは1.04×10−2(音速の100倍)以上である。逆に、速度V1が大きすぎる場合は、単位ステップあたりの粒子の動きが大きくなり、計算が破綻するおそれがある。従って、速度V1は、微小時間あたりの歪みが、好ましくは1.04(音速の10000倍)以下、さらに好ましくは0.104(音速の1000倍)以下である。
また、微小時間については、適宜設定されうる。なお、微小時間が大きいと、分子の熱振動を再現することができなくなり、計算が破綻するおそれがある。逆に、微小時間が小さいと、計算時間が増大するおそれがある。このため、微小時間は、好ましくは0.003τ以上であり、また、好ましくは0.012τ以下である。なお、1τのSI単位系への変換については、想定しているポリマーあるいは着目する物理量によっても異なる。一般的には、5.0×10−15s〜5.0×10−10sである。
高分子材料モデル16に与えられる歪みについては、適宜定義されうる。なお、歪みが小さいと、後述する物理量計算工程S23において、高分子材料モデル16に十分な破壊が生じないおそれがある。逆に、歪みが大きいと、各分子鎖モデル11に大きな力が作用し、計算落ちが生じるおそれがある。このため、本実施形態では、高分子材料モデル16の応力が、高分子材料の引張強さの1%〜100%になる歪みが与えられるのが望ましい。なお、高分子材料の引張強さは、実際の伸長試験に基づいて予め測定された値が設定されるのが望ましい。また、高分子材料の引張強さは、例えば、高速伸長工程S22において取得された応力−歪み曲線に基づいて決定されてもよい。
高速伸長工程S22において、伸長後の高分子材料モデル16のポアソン比は、高分子材料を用いた一軸引張試験時のポアソン比よりも小さいのが望ましい。これにより、小さな歪みでも、高分子材料モデル16の体積が増大されうる。この高分子材料モデル16の体積の増大により、ボイドを短時間で生じさせることができ、計算時間が短縮されうる。なお、本実施形態では、z軸方向の長さLzのみが拡大するように、高分子材料モデル16の変形が計算されるため、ポアソン比がゼロである。
次に、本実施形態のシミュレーション工程S2では、高速伸長工程S22で伸ばされた高分子材料モデル16の大きさが拘束され、その歪みを一定に保持した状態で物理量が計算される(物理量計算工程S23)。本実施形態において、「高分子材料モデル16の大きさが拘束される」とは、z軸方向に伸長された空間2の大きさが拘束されることを示している。図12に示されるように、高分子材料モデル16の大きさが拘束されると、高分子材料モデル16に与えられた歪みが、全ての分子鎖モデル11及び全てのフィラーモデル5に均一に伝達される。これらの歪み、各ポテンシャルP1〜P4(図8及び図10に示す)及び、運動方程式に基づいて、分子鎖モデル11及びフィラーモデル5の熱運動が計算される。
図14は、高分子材料モデル16の空孔を示す側面図である。上記のような熱運動により、粗視化粒子モデル13とフィラー粒子モデル6との間、及び、隣接する分子鎖モデル11の粗視化粒子モデル13、13間が接近又は離間する。これにより、空間2には、粗視化粒子モデル13及びフィラー粒子モデル6が存在しない小領域4が形成される。このような小領域4は、高分子材料モデル16に形成された空孔(ボイド)21として定義される。このような空孔21が計算されることにより、高分子材料モデル16の破壊が計算される。
図15は、破壊前の高分子材料モデル16の部分斜視図である。図16は、破壊後の高分子材料モデル16の部分斜視図である。図15及び図16は、図12に示した高分子材料モデル16の中央の一部分を拡大したものである。なお、図15及び図16では、分子鎖モデル11が省略して表示されている。図16に示されるように、本実施形態の高分子材料モデル16は、フィラーモデル5の周囲を覆っている空孔(ボイド)21が、z軸方向の広範囲に形成されており、局所的な破壊が防がれうる。
このように、本実施形態の物理量計算工程S23では、高分子材料モデル16の全体に歪みを均一に作用させて、高分子材料モデル16の破壊が計算されうる。従って、本実施形態の物理量計算工程S23では、従来のような局所的な破壊が生じるのを防ぐことができるため、高分子材料の破壊特性が正確に予測されうる。
しかも、本実施形態の高速伸長工程S22では、アフィン変形に基づいて、高分子材料モデル16の伸長が計算されているため、高分子材料モデル16の局所的な破壊が効果的に防がれうる。さらに、物理量計算工程S23では、高分子材料モデル16の大きさが拘束されているため、アフィン変形のような人為的操作が加わることがなく、高分子材料モデル16の破壊が計算されうる。従って、本実施形態では、計算精度が低下するのを防ぎうる。
なお、粗視化粒子モデル13、13間のポテンシャルP3(図10に示す)において、引力を含むように、カットオフ距離rcが設定されるのが望ましい。これにより、高分子材料モデル16において、空孔(ボイド)21を確実に生じさせることができる。なお、ポテンシャルP3のカットオフ距離rcは、例えば2.5σ程度に設定されるのが望ましい。
物理量計算工程S23では、高分子材料モデル16の物理量として、例えば、高分子材料モデル16の応力や、空孔(ボイド)21の偏りを示す指標が計算されるのが望ましい。本実施形態では、空孔(ボイド)21の偏りを示す指標としては、例えば、下記式(4)に示すシャノンエントロピーSが用いられている。
ここで、各定数及び変数は、次のとおりである。
Ndiv:小領域4の総数
i:小領域の添字(i=1〜Ndiv)
Lz:引張後の空間2のZ軸方向の長さ
Ci:(i−1)Lz/Ndiv≦z≦i・Lz/Ndivの範囲において、空孔21が存在する場合に、「1」を返し、空孔21が存在しない場合に、「0」を返す関数である。
上記式(4)で示されるシャノンエントロピーSは、各小領域4に空孔21が万遍なく存在する場合、最大値log(1/Ndiv)を取る。他方、空孔21の偏りが大きいほど、シャノンエントロピーSが小さくなる。従って、このようなシャノンエントロピーSが用いられることにより、空孔21の偏りが容易に判断されうる。
次に、本実施形態のシミュレーション方法では、高分子材料モデル16の物理量(破壊特性)が、良好か否かが判断される(工程S3)。本実施形態では、例えば、高分子材料モデル16の応力、空孔(ボイド)21の総量、又は、空孔(ボイド)21の偏り等に基づいて、高分子材料モデル16の物理量(破壊特性)が、良好か否かが判断される。工程S3では、高分子材料モデル16の物理量が良好であると判断された場合(工程S3で、「Y」)、高分子材料モデル16に基づいて、高分子材料が製造される(工程S4)。他方、工程S3では、高分子材料モデル16の物理量が良好でないと判断された場合は(工程S3で、「N」)、高分子材料モデル16の諸条件が変更されて(工程S5)、工程S1〜工程S3が再度実施される。
このように、高分子材料モデル16の物理量が良好になるまで、高分子材料モデル16の諸条件が変更されるため、良好な破壊特性を有する高分子材料が、効率よく設計されうる。
本実施形態のシミュレーション方法では、高速伸長工程S22において、微小時間あたりの歪みが10−5以上の速度V1で高分子材料モデル16が伸長されることにより、分子鎖モデル11の熱運動の計算が抑制される態様が例示されたが、これに限定されるわけではない。例えば、高速伸長工程S22に先立ち、分子鎖モデル11の熱運動を、物理量計算工程S23時の熱運動よりも小さくする熱運動抑制工程S24が含まれるのが望ましい。図17は、本発明の他の実施形態のシミュレーション工程S2の処理手順の一例を示すフローチャートである。なお、この実施形態において、前実施形態と同一の構成、及び、同一の処理が実施される工程については、同一の符号を付して、説明が省略される場合がある。
この実施形態の熱運動抑制工程S24は、分子鎖モデル11(図10に示す)に、物理量計算工程S23時の摩擦係数よりも大きい摩擦係数が定義される。この実施形態では、上記式(3)の運動方程式において、物理量計算工程S23時よりも大きな摩擦係数Γが定義されうる。
このような大きな摩擦係数Γが定義されることにより、高速伸長工程S22において、分子鎖モデル11の粗視化粒子モデル13の動きが強固に拘束される。これにより、この実施形態の高速伸長工程S22では、分子鎖モデル11の熱運動が、物理量計算工程S23時の熱運動よりも小に設定されうる。これにより、高速伸長工程S22では、高分子材料モデル16の一端16a及び他端16bにおいて計算されがちな局所的な破壊が、より効果的に防がれうる。
また、この実施形態のシミュレーション工程S2では、物理量計算工程S23に先立ち、分子鎖モデル11の熱運動が元に戻される(工程S25)。これにより、物理量計算工程S23では、本来の摩擦係数Γに基づいて、分子鎖モデル11の熱運動が計算されうる。従って、高分子材料モデル16の物理量が、精度良く計算されうる。
高速伸長工程S22時の摩擦係数Γについては、適宜設定されうる。高速伸長工程S22時の摩擦係数Γが小さいと、分子鎖モデル11の熱運動を十分に小さくできないおそれがある。このため、高速伸長工程S22の摩擦係数Γは、好ましくは、物理量計算工程S23時の摩擦係数Γの2倍以上であり、より好ましくは10倍以上である。他方、高速伸長工程S22時の摩擦係数Γが大きいと、アフィン変形によるフィラーモデル5の自発的な熱運動が阻害されるため、物理量計算工程S23において、フィラー粒子モデル6、6間の引力が大きくなり、計算落ちが発生するおそれがある。このため、高速伸長工程S22の摩擦係数Γは、好ましくは、物理量計算工程S23時の摩擦係数Γの1000倍以下であり、より好ましくは200倍以下である。
前実施形態の熱運動抑制工程S24では、分子鎖モデル11に、物理量計算工程S23時の摩擦係数Γよりも大きい摩擦係数Γが定義される場合が例示されたが、これに限定されるわけではない。例えば、高分子材料モデル16に、物理量計算工程S23時の温度よりも小さい温度が定義されてもよい。
このような小さな温度が定義されることにより、高速伸長工程S22において、分子鎖モデル11の熱運動が、物理量計算工程S23時の熱運動よりも小さくなる。これにより、高速伸長工程S22では、高分子材料モデル16の一端16a及び他端16bにおいて計算されがちな局所的な破壊が、より効果的に防がれうる。
高速伸長工程S22時の温度については、適宜設定されうる。高速伸長工程S22時の温度が大きいと、分子鎖モデル11の熱運動を十分に小さくできないおそれがある。このため、高速伸長工程S22の温度は、好ましくは、好ましくは、物理量計算工程S23時の温度の0.5倍以下であり、より好ましくは0.1倍以下である。
これまでの熱運動抑制工程S24では、物理量計算工程S23時の摩擦係数Γよりも大きい摩擦係数Γや、物理量計算工程S23時の温度よりも小さい温度が、高分子材料モデル16にそれぞれ定義される場合が例示されたが、これに限定されるわけではない。例えば、物理量計算工程S23時の摩擦係数Γよりも大きい摩擦係数Γと、物理量計算工程S23時の温度よりも小さい温度との双方が、高分子材料モデル16に定義されてもよい。これにより、高速伸長工程S22において、分子鎖モデル11の熱運動を、物理量計算工程S23時の熱運動よりも効果的に小さくすることができる。
図5に示されるように、本実施形態の分子鎖モデル11は、粗視化粒子モデル13と結合鎖モデル14とを含む粗視化モデルとして構成される場合が例示されたが、これに限定されるものではない。分子鎖モデル11は、図2に示した分子鎖のモノマーに基づいて、複数の粒子モデルとボンドモデルとで定義された全原子モデルや、ユナイテッドアトムモデルとして構成されてもよい。このような全原子モデル又はユナイテッドアトムモデルは、粗視化モデルに比べて、分子鎖の挙動が正確に表現されうるため、高分子材料の破壊特性がより正確に予測されうる。また、高分子材料モデル16は、空間2の内部に、フィラーモデル5及び分子鎖モデル11が配置されたモデルに限定されるわけではなく、例えば、連続体モデル(一例として、有限要素モデル等)として構成されてもよい。
以上、本発明の特に好ましい実施形態について詳述したが、本発明は図示の実施形態に限定されることなく、種々の態様に変形して実施しうる。
図4に示した処理手順に従って、予め定められた空間に、フィラーモデル及び分子鎖モデルが配置され、高分子材料モデルが設定された(実施例、比較例1、比較例2)。
実施例では、図11に示した処理手順に従って、高分子材料モデルが、微小時間あたりの歪みが10−5以上の速度で伸長される高速伸長工程と、高分子材料モデルの大きさを拘束し、その歪みを一定に保持した状態で物理量を計算する物理量計算工程とが実施された。また、実施例の高速伸長工程では、アフィン変形に基づいて、高分子材料の伸長が計算された。さらに、実施例では、高速伸長工程に先立ち、物理量計算工程時の摩擦係数よりも大きい摩擦係数が定義された(熱運動抑制工程)。
比較例1及び比較例2では、微小時間あたりの歪みが10−5未満の速度に基づいて、高分子材料モデルが徐々に引っ張られ、高分子材料モデルの物理量が逐次計算された。また、比較例1では、非アフィン変形に基づいて、高分子材料モデルの伸長が計算された。他方、比較例2では、アフィン変形に基づいて、高分子材料モデルの伸長が計算された。
実施例、比較例1及び比較例2では、上記式(4)に基づいて、空孔(ボイド)の偏りを示す指標であるシャノンエントロピーSが計算された。シャノンエントロピーSが大きいほど、空孔の偏りが小さく良好である。なお、ポテンシャル等の各数値は、明細書中の記載の通りであり、その他の共通仕様は次のとおりである。結果を表1に示す。
高分子材料モデル:
分子鎖モデル:
粗視化粒子モデルの個数:1000000個
粗視化粒子モデルの数密度:0.9σ−3
分子鎖モデル1個あたりの粗視化粒子モデルの個数:200個
フィラーモデル:
フィラー粒子モデルの体積分率:0.18
フィラーモデルの半径:10σ
フィラーモデル1個あたりのフィラー粒子モデルの個数:5000個
シミュレーション:
コンピュータのCPU:intel Core i7-4770(4コア)
ポアソン比:0
高分子材料モデルのひずみ:0.3(高分子材料の引張強さの100%)
実施例:
高速伸長工程:
速度V1:微小時間あたりの歪み3.12×10−4
(音速の3倍[ひずみ速度:0.053(1/τ)])
微小時間:0.006τ(6×10−13s)
高速伸長工程のシミュレーション時間:5.66τ
物理量計算工程のシミュレーション時間:1181.1τ
高速伸長工程の摩擦係数:物理量計算工程時の摩擦係数の10倍
比較例1、比較例2:
高速伸長工程:
速度V1:微小時間あたりの歪み1.47×10−6
(音速の0.015倍[ひずみ速度:0.00025(1/τ)])
微小時間:0.006τ(6×10−13s)
引張時間:1200τ
テストの結果、実施例は、比較例1に比べて、シャノンエントロピーが大であり、空孔の偏りが小さかった。従って、実施例は、比較例1に比べて、局所的な破壊が生じるのを防ぐことができ、高分子材料の破壊特性を正確に予測しうることが確認できた。
実施例と、比較例2とは、シャノンエントロピーが同一であった。しかしながら、比較例2では、物理量を計算する工程において、分子運動に対して、アフィン変形による人為的な操作が常に加わった。他方、実施例では、物理量計算工程において、高分子材料モデルの大きさが拘束されているため、分子運動に対して、アフィン変形による人為的な操作が加わらなかった。さらに、実施例では、高速伸長工程において、アフィン変形に基づいて、高分子材料モデルを伸長させているため、高分子材料モデルの一端及び他端において計算されがちな局所的な破壊が効果的に防がれた。従って、実施例は、比較例2に比べて、高分子材料の破壊特性を正確に予測できた。
実施例では、高分子材料モデルの大きさを拘束して、その歪みを一定に保持した状態で物理量が計算されている。ここで、歪みは、物理量計算工程において、破壊が比較的急速に進行する程度に大きく設定されている。このため、実施例では、微小時間あたりの歪みが10−5未満の速度で徐々に変形させながら物理量が計算される比較例1及び比較例2に比べて、物理量計算工程のシミュレーション時間の一部を省略することができる。従って、実施例は、比較例1及び比較例2に比べて、計算時間を短縮しうる。