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

Abstract

【課題】高分子材料の緩和弾性率を短時間で計算する。【解決手段】コンピュータ１を用いて、高分子材料の緩和弾性率を計算するためのシミュレーション方法である。シミュレーション方法は、コンピュータ１が、高分子材料モデル１１を用いて、予め定められた第１温度の緩和弾性率Ｇ（ｔ）を計算する緩和弾性率計算工程Ｓ３を含んでいる。緩和弾性率計算工程Ｓ３は、高分子材料モデル１１に、第１温度よりも高い第２温度を設定する工程Ｓ３１と、第２温度が設定された高分子材料モデル１１の緩和弾性率Ｇ（ｔ）を計算する計算工程Ｓ３３と、高分子材料の時間−温度換算則に基づいて、第２温度の高分子材料モデル１１の緩和弾性率Ｇ（ｔ）から第１温度の高分子材料モデル１１の緩和弾性率Ｇ（ｔ）に変換する工程Ｓ３４とを含む。【選択図】図３

Description

本発明は、緩和弾性率を計算するための高分子材料のシミュレーション方法に関する。

近年、ゴム等の高分子材料の開発のために、高分子材料の性質を、コンピュータを用いて評価するためのシミュレーション方法（数値計算）が種々提案されている（例えば、下記特許文献１参照）。

この種のシミュレーション方法では、先ず、高分子材料の分子鎖をモデル化した分子鎖モデルが設定される。次に、予め定められた空間内に分子鎖モデルが配置され、高分子材料モデルが設定される。そして、分子動力学（ Molecular Dynamics : ＭＤ ）に基づいて、高分子材料モデルの緩和計算が行われる。

特開２０１３−１９５２２０号公報

上記のようなシミュレーション方法を利用して、例えば、緩和弾性率が計算される場合がある。緩和弾性率とは、歪が与えられた粘弾性体の弾性率の変化を示す指標である。緩和弾性率を計算するためには、予め定められた温度（例えば、室温）において、高分子材料モデルに緩和現象を発現させる必要がある。

しかしながら、高分子材料モデルの緩和現象が発現するまでに多くの時間が必要であった。従って、緩和弾性率の計算を短時間で計算することができるシミュレーション方法が強く求められていた。

本発明は、以上のような実状に鑑み案出されたもので、緩和弾性率を短時間で計算することができる高分子材料のシミュレーション方法を提供することを主たる目的としている。

本発明は、コンピュータを用いて、高分子材料の緩和弾性率を計算するためのシミュレーション方法であって、前記コンピュータに、前記高分子材料の分子鎖をモデル化した分子鎖モデルを設定する工程と、前記コンピュータに、前記分子鎖モデルを予め定められた空間内に配置して高分子材料モデルを設定する工程と、前記コンピュータが、前記高分子材料モデルを用いて、予め定められた第１温度の前記緩和弾性率を計算する緩和弾性率計算工程とを含み、前記緩和弾性率計算工程は、前記高分子材料モデルに、前記第１温度よりも高い第２温度を設定する工程と、前記第２温度が設定された前記高分子材料モデルの緩和弾性率を計算する計算工程と、前記高分子材料の時間−温度換算則に基づいて、前記第２温度の前記高分子材料モデルの緩和弾性率から前記第１温度の前記高分子材料モデルの緩和弾性率に変換する工程とを含むことを特徴とする。

本発明に係る前記高分子材料のシミュレーション方法において、前記計算工程は、前記高分子材料モデルに緩和現象が発現する前の前記緩和弾性率を取得する工程と、前記緩和現象が発現した後の前記緩和弾性率を取得する工程とを含むのが望ましい。

本発明に係る前記高分子材料のシミュレーション方法において、前記第２温度は、５００Ｋ以下であるのが望ましい。

本発明に係る前記高分子材料のシミュレーション方法において、前記高分子材料モデルは、少なくとも８本の前記分子鎖モデルを含むのが望ましい。

本発明に係る前記高分子材料のシミュレーション方法において、前記分子鎖モデルは、全原子モデル又はユナイテッドアトムモデルであるのが望ましい。

本発明に係る前記高分子材料のシミュレーション方法において、前記高分子材料モデルの重合度は、５０以下であるのが望ましい。

本発明は、コンピュータが、高分子材料モデルを用いて、予め定められた第１温度での緩和弾性率を計算する緩和弾性率計算工程とを含んでいる。緩和弾性率計算工程では、高分子材料モデルに、第１温度よりも高い第２温度を設定する工程と、第２温度が設定された高分子材料モデルの緩和弾性率を計算する計算工程とが含まれる。高分子材料モデルは、設定される温度が高くなるほど、緩和現象が短時間で発現する傾向がある。このため、計算工程では、第２温度において、緩和現象が発現した高分子材料モデルの緩和弾性率が、短時間で計算されうる。

さらに、緩和弾性率計算工程では、高分子材料の時間−温度換算則に基づいて、第２温度の高分子材料モデルの緩和弾性率から第１温度の高分子材料モデルの緩和弾性率に変換する工程が含まれている。一般に、高分子材料は、時間−温度換算則に基づいて、時間と温度とが等価に変換できるという性質を有している。このため、第２温度の高分子材料モデルの緩和弾性率に基づいて、第１温度の高分子材料モデルの緩和弾性率が容易に求められる。従って、本発明は、高分子材料モデルの緩和弾性率を、短時間で計算されうる。

本発明のシミュレーション方法を実行するためのコンピュータの斜視図である。 ポリスチレンの構造式である。 本実施形態のシミュレーション方法の処理手順の一例を示すフローチャートである。 本実施形態の分子鎖モデルの概念図である。 （ａ）〜（ｃ）は、ポテンシャルを説明する分子鎖モデルの部分図である。 本実施形態の高分子材料モデルの概念図である。 本実施形態の緩和弾性率計算工程の処理手順の一例を示すフローチャートである。 本実施形態の計算工程の一例を示すフローチャートである。 第２温度の緩和弾性率の対数と、時間幅の対数との関係を示すグラフである。 図９の一部を拡大して示すグラフである。 緩和された高分子材料モデルの概念図である。 本実施形態の変換工程の処理手順の一例を示すフローチャートである。 第２温度の緩和弾性率から、第１温度の緩和弾性率に変換する工程を説明するグラフである。 比較例の第１温度の緩和弾性率の対数と、時間幅の対数との関係を示すグラフである。

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

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

高分子材料としては、例えば、ゴム、樹脂又はエラストマー等が含まれる。本実施形態の高分子材料としては、ポリスチレンが例示される。図２は、ポリスチレンの構造式である。このポリスチレンを構成する分子鎖Ｍｃは、スチレンから構成されるモノマーが、重合度Ｍｎで連結されている。なお、本発明のシミュレーション方法では、ポリスチレン以外の高分子材料が用いられてもよい。

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

本実施形態の分子鎖モデル２は、複数の粒子モデル３と、粒子モデル３、３間を結合するボンドモデル４とを含む全原子モデルとして構成されている。これらの粒子モデル３及びボンドモデル４は、分子鎖Ｍｃのモノマーをなす単位構造６（図２に示す）に基づいて、互いに連結されることにより、モノマーモデル７が設定される。このモノマーモデル７が、重合度Ｍｎに基づいて連結されることにより、分子鎖モデル２が設定される。

粒子モデル３は、後述する分子動力学計算に基づいたシミュレーションにおいて、運動方程式の質点として取り扱われる。即ち、粒子モデル３は、質量、直径、電荷、又は、初期座標などのパラメータが定義される。本実施形態の粒子モデル３は、分子鎖Ｍｃの炭素原子をモデル化した炭素粒子モデル３Ｃ、及び、分子鎖Ｍｃの水素原子をモデル化した水素粒子モデル３Ｈを含んでいる。

ボンドモデル４は、粒子モデル３、３間を拘束するものである。本実施形態のボンドモデル４は、炭素粒子モデル３Ｃ、３Ｃを連結する主鎖４ａ、及び、炭素粒子モデル３Ｃと水素粒子モデル３Ｈとの間を連結する側鎖４ｂとを含んでいる。

粒子モデル３、３間には、相互作用（斥力及び引力を含む）が生じさせるポテンシャルが定義される。図４に示されるように、ポテンシャルには、ボンドモデル４を介して隣り合う粒子モデル３、３間に定義される第１ポテンシャルＰ１と、ボンドモデル４を介さずに隣り合う粒子モデル３、３間に定義される第２ポテンシャルＰ２とが定義される。なお、複数の分子鎖モデル２が定義される場合は、分子鎖モデル２、２間の粒子モデル３、３間にも、第２ポテンシャルＰ２が定義される。

図５（ａ）〜（ｃ）は、第１ポテンシャルＰ１を説明する分子鎖モデル２の部分図である。図５（ａ）〜（ｂ）に示されるように、分子鎖モデル２には、各粒子モデル３、３間の結合長さである結合長ｒ、及び、ボンドモデル４を介して連続する３つの粒子モデル３がなす角度である結合角θが定義されている。さらに、図５（ｃ）に示されるように、分子鎖モデル２には、ボンドモデル４を介して連続する４つの粒子モデル３において、隣り合う３つの粒子モデル３が作る一方の平面５Ａと他方の平面５Ｂとのなす角度ある二面角φが定義される。結合長ｒ、結合角θ及び二面角φは、分子鎖モデル２に作用する外力又は内力によって変化する。

結合長ｒ、結合角θ及び二面角φは、下記式（１）で定義される結合ポテンシャルＵ_bond（ｒ）、下記式（２）で定義される結合角ポテンシャルＵ_Angle（θ）、及び、下記式（３）で定義される結合二面角ポテンシャルＵ_torsion（φ）によって設定される。



ここで、各定数及び変数は、次のとおりである。
ｒ：結合長
ｒ₀：平衡長
ｋ₁、ｋ₂：ばね定数
θ：結合角
θ₀：平衡角度
k₃：二面角ポテンシャルの強度
Ｎ−1：二面角ポテンシャル多項式の次数
φ：二面角
Ａ_n：二面角定数
なお、結合長ｒ及び平衡長ｒ₀は、各粒子モデル３の中心（図示省略）間の距離として定義される。

第１ポテンシャルＰ１は、炭素粒子モデル３Ｃ、３Ｃ間、及び、炭素粒子モデル３Ｃと水素粒子モデル３Ｈとの間で、それぞれ大きさが異なっている。これは、各第１ポテンシャルＰ１において、各結合ポテンシャルＵ_bond（ｒ）、結合角ポテンシャルＵ_Angle（θ）、及び、結合二面角ポテンシャルＵ_torsion（φ）の各定数に、それぞれ異なる値が設定されているためである。なお、各定数は、論文（J. Phys. Chem. 94, 8897 (1990））に基づいて、分子鎖Ｍｃの構造に応じて設定されるのが望ましい。

第２ポテンシャルＰ２（図４に示す）は、下記式（４）で定義される ＬＪポテンシャルＵ_LJ（ｒ_ij）である。

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

第２ポテンシャルＰ２は、粒子モデル間の距離ｒ_ijがσよりも小さくなるほど、粒子モデル３、３間に作用する斥力が大きくなる。また、第２ポテンシャルＰ２は、粒子モデル間の距離ｒ_ijがσになるときに最小となり、粒子モデル３、３間に斥力や引力は働かない。さらに、第２ポテンシャルＰ２は、粒子モデル間の距離ｒ_ijがσよりも大になると、粒子モデル３、３間に作用する引力が働く。このように、第２ポテンシャルＰ２は、粒子モデル間の距離ｒ_ijに応じて、斥力及び引力を定義することができる。

第２ポテンシャルＰ２は、炭素粒子モデル３Ｃ、３Ｃ間、炭素粒子モデル３Ｃと水素粒子モデル３Ｈとの間、及び、水素粒子モデル３Ｈ、３Ｈ間で、それぞれ大きさが異なっている。これは、各第２ポテンシャルＰ２において、上記式（４）の定数に、それぞれ異なる値が設定されているためである。なお、各定数は、例えば、上記論文に基づいて、適宜設定されうる。

このような分子鎖モデル２は、例えば（株）ＪＳＯＬ社製のＪ−ＯＣＴＡというソフトウェアを用いて作成することができる。分子鎖モデル２は、コンピュータ１で取り扱い可能な数値データであり、コンピュータ１に記憶される。

次に、本実施形態のシミュレーション方法では、分子鎖モデル２が予め定められた空間９内に配置されることにより、コンピュータ１に高分子材料モデルが設定される（工程Ｓ２）。図６は、本実施形態の高分子材料モデルの概念図である。

本実施形態の空間９は、互いに向き合う三対の平面１０、１０を有する立方体として定義されている。各平面１０には、周期境界条件が定義されている。これにより、空間９では、例えば、一方の平面１０ａから出て行った分子鎖モデル２の一部が、反対側の平面１０ｂから入ってくるように計算されうる。従って、一方の平面１０ａと、反対側の平面１０ｂとが連続している（繋がっている）ものとして扱われる。

空間９の一辺の長さＬ１は、適宜設定されうる。本実施形態の長さＬ１は、分子鎖モデル２の慣性半径（図示省略）の２倍以上が望ましい。慣性半径は、後述する分子動力学計算において、分子鎖モデル２の拡がりを示すパラメータである。このような空間９では、分子鎖モデル２の回転運動がスムーズに計算されうる。さらに、空間９の大きさは、例えば１atmで安定な体積に設定される。このような空間９は、高分子材料の少なくとも一部の体積が定義されうる。

工程Ｓ２では、複数の分子鎖モデル２が、空間９内に配置される。これにより、高分子材料モデル１１が設定される。本実施形態では、例えば、オペレータ等によって、複数の分子鎖モデル２が、空間９内にランダムに配置される。従って、工程Ｓ２において、高分子材料モデル１１は、緩和計算されていない初期の高分子材料モデルである。高分子材料モデル１１は、コンピュータ１に記憶される。

空間９に配置される分子鎖モデル２の本数については、適宜設定されうる。分子鎖モデル２の本数が少ないと、後述する分子動力学計算において、一方の平面１０ａから出て行き、かつ、反対側の平面１０ｂから入ってきた分子鎖モデル２の一端側と、この分子鎖モデル２の他端側とが絡まって、不自然な分子運動を起こすおそれがある。逆に、分子鎖モデル２の本数が多くても、運動方程式の質点として取り扱われる粒子モデル３が増大し、後述する緩和弾性率計算工程Ｓ３において、多くの計算時間を要するおそれがある。このような観点より、空間９に配置される分子鎖モデル２の本数は、好ましくは８本以上であり、また、好ましくは、１００本以下である。

次に、本実施形態のシミュレーション方法では、コンピュータ１が、高分子材料モデル１１を用いて、予め定められた第１温度での緩和弾性率Ｇ（ｔ）を計算する（緩和弾性率計算工程Ｓ３）。

緩和弾性率Ｇ（ｔ）は、所定の時間幅ｔにおいて、歪が与えられた粘弾性体の弾性率の変化を示す指標である。緩和弾性率Ｇ（ｔ）は、下記式（５）で計算される。

ここで、
Ｖ：空間の体積
ｋ_Ｂ：ボルツマン定数
Ｔ：絶対温度
σ_xy：応力
ｘｙ：任意の直交する２方向
τ：時刻
ｔ：時間幅

上記式（５）において、＜σ_xy（ｔ＋τ）×σ_xy（τ）＞は、所定の時間内において、時刻τの応力σ_xyと、時刻（ｔ＋τ）の応力σ_xyとの積を、あらゆる時刻τについて平均（アンサンブル平均）したものである。

緩和弾性率Ｇ（ｔ）の計算では、高分子材料モデル１１に緩和現象を発現させる必要がある。なお、高分子材料は、設定される温度が高くなるほど、緩和現象が短時間で発現する傾向がある。

このような観点より、本実施形態の緩和弾性率計算工程Ｓ３では、第１温度よりも高い第２温度が高分子材料モデル１１に設定されることにより、緩和現象を短時間で発現させている。そして、第２温度の高分子材料モデル１１の緩和弾性率（以下、単に「第２温度の緩和弾性率」ということがある。）Ｇ（ｔ）が、第１温度の高分子材料モデル１１の緩和弾性率（以下、単に「第１温度の緩和弾性率」ということがある。）Ｇ（ｔ）に変換されることにより、第１温度の緩和弾性率Ｇ（ｔ）が計算される。図７は、本実施形態の緩和弾性率計算工程Ｓ３の処理手順の一例を示すフローチャートである。

本実施形態の緩和弾性率計算工程Ｓ３では、先ず、高分子材料モデル１１（図６に示す）に、第１温度よりも高い第２温度が設定される（工程Ｓ３１）。本実施形態の第１温度については、高分子材料が使用される実際の温度に基づいて、適宜設定されうる。本実施形態の第１温度は、例えば、室温に基づいて、２９０Ｋ〜３０５Ｋ程度に設定される。

第２温度は、第１温度よりも高い温度であれば、適宜設定されうる。なお、第２温度が高いと、後述する分子動力学計算に基づいたシミュレーションにおいて、分子鎖モデル２の構造の変化や、計算落ちを招くおそれがある。逆に、第２温度が低いと、高分子材料モデル１１の緩和現象が、短時間で発現しないおそれがある。このような観点より、第２温度は、好ましくは５００Ｋ以下であり、また、好ましくは３５０Ｋ以上である。これらの第１温度及び第２温度は、コンピュータ１に記憶される。

次に、本実施形態の緩和弾性率計算工程Ｓ３では、緩和弾性率Ｇ（ｔ）を計算するのに先立って、分子鎖モデル２の初期配置が緩和される（工程Ｓ３２）。構造緩和の計算では、高分子材料モデル１１を用いた分子動力学計算が行われる。

分子動力学計算では、例えば、空間９について所定の時間、分子鎖モデル２が古典力学に従うものとして、ニュートンの運動方程式が適用される。そして、各時刻での粒子モデル３の動きが、単位時間毎に追跡される。このような構造緩和の計算は、例えば（株）ＪＳＯＬ社製のソフトマテリアル総合シミュレーター（Ｊ−ＯＣＴＡ）に含まれるＣＯＧＮＡＣを用いて処理することができる。

本実施形態の分子動力学計算では、空間９において、圧力（例えば、１atm）及び第２温度が一定（ＮＰＴ一定）に保たれる。これにより、工程Ｓ３２では、実際の高分子材料の分子運動に近似させて、分子鎖モデル２の初期配置が精度よく緩和されうる。さらに、本実施形態では、高分子材料モデル１１に、第１温度よりも高い第２温度が設定されているため、例えば、第１温度が設定された高分子材料モデル１１に比べて、短時間で緩和される。

本実施形態の工程Ｓ３２では、分子鎖モデル２の人為的な初期配置が排除されたとみなすことができるまで、単位時間毎に分子動力学計算が行われる。人為的な初期配置が排除されたか否かの判断は、適宜設定されうる。本実施形態では、例えば従来と同様に、分子鎖モデル２の末端間ベクトルの自己相関関数が１／ｅ以下になったか否かにより、分子鎖モデル２の初期配置が緩和できたか否かが判断されるのが望ましい。これにより、工程Ｓ３２では、分子鎖モデル２の人為的な初期配置が、確実に排除される。

次に、本実施形態の緩和弾性率計算工程Ｓ３では、第２温度が設定された高分子材料モデル１１の緩和弾性率Ｇ（ｔ）が、時間幅ｔ毎に計算される（計算工程Ｓ３３）。図８は、本実施形態の計算工程Ｓ３３の一例を示すフローチャートである。

本実施形態の計算工程Ｓ３３は、先ず、人為的な初期配置が緩和された高分子材料モデル１１を用いて、高分子材料モデル１１の緩和弾性率Ｇ（ｔ）が、時間幅ｔ（単位時間）毎に計算される（工程Ｓ３３１）。工程Ｓ３３１では、先ず、工程Ｓ３２と同様に、空間９において、圧力（例えば、１atm）及び第２温度が一定（ＮＰＴ一定）に保たれる。そして、工程Ｓ３３１では、高分子材料モデル１１を用いた分子動力学計算が行われ、時間幅ｔ毎に、高分子材料モデル１１の緩和弾性率Ｇ（ｔ）が計算される。緩和弾性率Ｇ（ｔ）は、時間幅ｔ毎に、コンピュータ１に記憶される。

図９は、第２温度の緩和弾性率の対数ｌｏｇ（Ｇ（ｔ））と、時間幅の対数ｌｏｇ（ｔ）との関係を示すグラフである。このグラフでは、ｌｏｇ（ｔ）の増加とともに、ｌｏｇ（Ｇ（ｔ））が減少している。これは、高分子材料モデル１１の構造が、時間幅ｔの増加とともに、徐々に緩和されていることを示している。このようなグラフは、時間幅ｔ毎に計算された緩和弾性率Ｇ（ｔ）に基づいて作成されうる。なお、時間幅ｔの間隔は、１ps〜１００ps程度が望ましい。

次に、本実施形態の計算工程Ｓ３３では、高分子材料モデル１１に緩和現象が発現したか否かが判断される（工程Ｓ３３２）。図１０は、図９の一部を拡大して示すグラフである。緩和現象が発現した否かの判断は、図１０のグラフにおいて、緩和弾性率の対数ｌｏｇ（Ｇ（ｔ））の傾きαが、−０．６〜−０．４であるか否かで判断される。これは、高分子鎖のＲｏｕｓｅモデルに基づくものである。Ｒｏｕｓｅモデルによれば、ｌｏｇＧ（ｔ）とｌｏｇｔとの関係において、傾き−１／２のときに、高分子鎖の緩和がみられることが知られている。Ｒｏｕｓｅモデルの詳細については、例えば、論文（Alexei E. Likhtman, Sathish K. Sukumaran, and Jorge Ramirez 著、「Linear Viscoelasticity from Molecular Dynamics Simulation of Entangled Polymers」、Macromolecules, 40, 6748 (2007)）に示されている。

図１０に示されるように、傾きαは、例えば、現在の時間幅の対数ｌｏｇ（ｔａ）、及び、現在よりも前の時間幅の対数ｌｏｇ（ｔｂ）において、各緩和弾性率の対数ｌｏｇ（Ｇ（ｔａ））、ｌｏｇ（Ｇ（ｔｂ））間を結ぶ直線Ｌｓに基づいて定義される。なお、上記数値範囲（−０．６〜−０．４）の傾きαは、緩和現象が発現している領域Ｒａ（図９に示す）内において、各時間幅の対数ｌｏｇ（ｔａ）、ｌｏｇ（ｔｂ）の双方が選択された場合に設定される。

各時間幅の対数ｌｏｇ（ｔａ）及びｌｏｇ（ｔｂ）の差Ａｄについては、適宜設定されうる。なお、差Ａｄが大きいと、各時間幅の対数ｌｏｇ（ｔａ）及びｌｏｇ（ｔｂ）の双方が、領域Ｒａ内で選択されなくなり、緩和現象の発現を正確に判断できなくなるおそれがある。逆に、差Ａｄが小さいと、緩和弾性率Ｇ（ｔ）の誤差等により、緩和現象が発現する前に−０．６〜−０．４の傾きαが設定され、緩和現象の発現を正確に判断できないおそれがある。このような観点より、差Ａｄは、好ましくは１０００ns以下であり、また、好ましくは１ns以上である。

工程Ｓ３３２では、高分子材料モデル１１に緩和現象が発現したと判断された場合（工程Ｓ３３２で、「Ｙ」）、工程Ｓ３３１で計算された緩和弾性率Ｇ（ｔ）が、緩和現象が発現した後の緩和弾性率Ｇ（ｔ）として取得される（工程Ｓ３３３）。一方、高分子材料モデル１１に緩和現象が発現していないと判断された場合（工程Ｓ３３２で、「Ｎ」）、工程Ｓ３３１で計算された緩和弾性率Ｇ（ｔ）が、緩和現象が発現する前の緩和弾性率Ｇ（ｔ）として取得される（工程Ｓ３３４）。そして、工程Ｓ３３３又は工程Ｓ３３４が実施された後は、次の工程Ｓ３３５が実施される。

次に、計算工程Ｓ３３では、予め定められた終了時間Ｅｔが経過したか否かが判断される（工程Ｓ３３５）。工程Ｓ３３５では、予め定められた終了時間Ｅｔが経過したと判断された場合（工程Ｓ３３５で、「Ｙ」）、計算工程Ｓ３３の一連の処理が終了され、次の変換工程Ｓ３４が実施される。他方、終了時間Ｅｔが経過していないと判断された場合は（工程Ｓ３３５で、「Ｎ」）、時間幅ｔが一つ増加されて（工程Ｓ３３６）、工程Ｓ３３１〜工程Ｓ３３５が再度実施される。

これにより、計算工程Ｓ３３では、計算開始から終了時間Ｅｔが経過するまでの間、緩和現象が発現する前の緩和弾性率Ｇ（ｔ）、及び、緩和現象が発現した後の緩和弾性率Ｇ（ｔ）の双方が、時間幅ｔ毎に計算される。そして、図９に示したグラフが作成されうる。本実施形態の計算工程Ｓ３３では、第１温度よりも高い第２温度に基づいて、緩和弾性率Ｇ（ｔ）が短時間で計算されるため、第１温度に基づいて緩和弾性率Ｇ（ｔ）が計算されていた従来の方法に比べて、緩和現象が発現した緩和弾性率Ｇ（ｔ）を、短時間で計算しうる。緩和弾性率Ｇ（ｔ）は、コンピュータ１に記憶される。図１１は、緩和後の高分子材料モデル１１を示す図である。

計算工程Ｓ３３の終了時間Ｅｔについては、適宜選択されうる。なお、終了時間Ｅｔが小さいと、高分子材料モデル１１の緩和現象が発現する前に、計算工程Ｓ３３が終了するおそれがある。逆に、終了時間Ｅｔが大きくても、高分子材料モデル１１の緩和現象が終了した後も、緩和弾性率Ｇ（ｔ）が継続して計算されてしまい、計算コストが増大するおそれがある。このような観点より、終了時間Ｅｔは、好ましくは１０ns以上であり、また、好ましくは１０００ns以下である。

緩和現象が発現するまでの時間は、分子鎖モデル２の重合度Ｍｎに依存している。重合度Ｍｎが大きいと、計算工程Ｓ３３の終了時間Ｅｔまでに、緩和現象が発現しないおそれがある。このような観点より、重合度Ｍｎは、好ましくは５０以下であり、さらに好ましくは３０以下である。他方、重合度Ｍｎが小さいと、分子鎖モデル２の分子量が小さくなり、シミュレーション精度が低下するおそれがある。このような観点より、重合度Ｍｎは、好ましくは３以上であり、さらに好ましくは５以上である。

次に、本実施形態の緩和弾性率計算工程Ｓ３では、時間−温度換算則に基づいて、第２温度の緩和弾性率Ｇ（ｔ）から、第１温度の緩和弾性率Ｇ（ｔ）に変換される（変換工程Ｓ３４）。図１２は、本実施形態の変換工程Ｓ３４の処理手順の一例を示すフローチャートである。図１３は、第２温度の緩和弾性率Ｇ（ｔ）から、第１温度の緩和弾性率Ｇ（ｔ）に変換する工程を説明するグラフである。

本実施形態の変換工程Ｓ３４では、先ず、第２温度の緩和弾性率Ｇ（ｔ）から、第１温度の緩和弾性率Ｇ（ｔ）に変換するためのシフトファクター（移動因子）Ａｔが求められる（工程Ｓ３４１）。

一般に、高分子材料は、時間−温度換算則に基づいて、時間と温度とが等価に変換できるという性質を有している。このため、図１３に示されるように、第２温度の緩和弾性率の対数ｌｏｇ（Ｇ（ｔ））が、時間幅の対数ｌｏｇ（ｔ）の軸方向に平行に移動されることにより、異なる温度の緩和弾性率の対数ｌｏｇ（Ｇ（ｔ））が求められる。

シフトファクターＡｔは、第２温度の時間幅の対数ｌｏｇ（ｔ）を除することにより、第２温度の各緩和弾性率の対数ｌｏｇ（Ｇ（ｔ））を、第１温度の各緩和弾性率の対数ｌｏｇ（Ｇ（ｔ））に移動させるための比率である。このようなシフトファクターＡｔは、従来の方法に基づいて、高分子材料を用いた実験、又は、高分子材料モデル１１を用いたシミュレーションによって求められうる。本実施形態では、ＷＬＦ（Williams-Landel-Ferry）式が用いられることにより、シフトファクターＡｔが求められた。シフトファクターＡｔは、コンピュータ１に記憶される。

次に、本実施形態の変換工程Ｓ３４では、シフトファクターＡｔを用いて、第２温度の緩和弾性率Ｇ（ｔ）から、第１温度の緩和弾性率Ｇ（ｔ）に変換される（工程Ｓ３４２）。

図１３に示されるように、第２温度の各時間幅の対数ｌｏｇ（ｔ）は、シフトファクターＡｔが乗じられることにより、時間幅の対数ｌｏｇ（ｔ）の軸方向に平行に移動される。本実施形態では、第２温度の各緩和弾性率の対数ｌｏｇ（Ｇ（ｔ））が、時間幅の対数ｌｏｇ（ｔ）が大きくなる方向に移動されている。これにより、第１温度の各緩和弾性率の対数ｌｏｇ（Ｇ（ｔ））が求められる。

このように、本実施形態の緩和弾性率計算工程Ｓ３では、シフトファクターＡｔに基づいて、第２温度の各緩和弾性率の対数ｌｏｇ（Ｇ（ｔ））が、第１温度の各緩和弾性率の対数ｌｏｇ（Ｇ（ｔ））に、容易に変換されうる。そして、第１温度の各緩和弾性率の対数ｌｏｇ（Ｇ（ｔ））が、自然数に戻されることにより、第１温度の緩和弾性率Ｇ（ｔ）が求められる。第１温度の各緩和弾性率Ｇ（ｔ）は、コンピュータ１に記憶される。

上述したように、本実施形態のシミュレーション方法では、緩和現象が発現した第２温度の緩和弾性率Ｇ（ｔ）が、短時間で計算されうる。第２温度の緩和弾性率Ｇ（ｔ）は、シフトファクターＡｔに基づいて、第１温度の緩和弾性率Ｇ（ｔ）に容易に変換される。従って、本実施形態のシミュレーション方法では、緩和現象が発現した第１温度の緩和弾性率Ｇ（ｔ）が、短時間で計算されうる。

次に、シミュレーション方法では、第１温度の緩和弾性率Ｇ（ｔ）が、許容範囲内であるか否かが判断される（工程Ｓ４）。この工程Ｓ４では、第１温度の緩和弾性率Ｇ（ｔ）が許容範囲内であると判断された場合（工程Ｓ４において、「Ｙ」）、図６に示した高分子材料モデル１１に基づいて、高分子材料が製造される（工程Ｓ５）。一方、第１温度の緩和弾性率Ｇ（ｔ）が許容範囲外と判断された場合は（工程Ｓ４において、「Ｎ」）、分子鎖モデル２の諸条件が変更されて（工程６）、工程Ｓ１〜工程Ｓ４が再度実施される。このように、本実施形態では、第１温度の緩和弾性率Ｇ（ｔ）が許容範囲内になるまで、分子鎖モデル２の諸条件が変更されるため、所望の性能を有する高分子材料が、効率よく設計されうる。しかも、本実施形態では、第１温度の緩和弾性率Ｇ（ｔ）が短時間で計算されるため、所望の性能を有する高分子材料が、短時間で計算されうる。

本実施形態の分子鎖モデル２は、全原子モデルとして構成されたものが例示されたが、これに限定されるものではない。分子鎖モデル２は、例えば、炭素原子と、炭素原子に結合した水素原子とを一体化して、一つの粒子モデル（図示省略）として扱うユナイテッドアトムモデル（ united atom model ）として構成されてもよい。このような分子鎖モデル２は、高分子鎖のモノマーの配置や、シス構造又はトランス構造を維持しつつ、水素原子を省略することができるため、計算時間を短縮することができる。

また、分子鎖モデル２は、例えば、モノマーをなす単位構造６（図２に示す）が置換された粗視化粒子モデル（図示省略）と、粗視化粒子モデル間を連結する結合鎖モデル（図示省略）とを含む粗視化モデル（図示省略）として構成されてもよい。このような粗視化モデルは、全原子モデルやユナイテッドアトムモデルに比べて、大きな時間スケールで、第１温度の緩和弾性率Ｇａ（ｔ）を計算することができる。

上述したように、粗視化モデルの分子鎖モデル２を用いて取得された緩和弾性率Ｇ（ｔ）は、全原子モデルの分子鎖モデル２を用いて取得された緩和弾性率Ｇ（ｔ）に比べて、大きな時間スケールで計算される。他方、全原子モデルの分子鎖モデル２を用いて取得された緩和弾性率Ｇ（ｔ）は、小さな時間スケールで計算される反面、実際の高分子材料の緩和弾性率に、高い精度で近似させることができる。従って、全原子モデルの分子鎖モデル２を用いて取得された緩和弾性率Ｇ（ｔ）と、粗視化モデルの分子鎖モデル２を用いて取得された緩和弾性率Ｇ（ｔ）とが組み合わされることにより、正確な緩和弾性率Ｇａ（ｔ）が、大きな時間スケールで計算されうる。

全原子モデルの緩和弾性率Ｇ（ｔ）と、粗視化モデルの緩和弾性率Ｇ（ｔ）とを組み合わせる方法としては、例えば、文献（第６２回高分子討論会予稿集1J13）に基づいて実施されうる。また、この実施形態では、全原子モデルの緩和弾性率Ｇ（ｔ）と、粗視化モデルの緩和弾性率Ｇ（ｔ）とが組み合わされるのに先立ち、全原子モデル及び粗視化モデルのそれぞれにおいて、高分子材料モデル１１に緩和現象が発現する前の緩和弾性率Ｇ（ｔ）と、緩和現象が発現した後の緩和弾性率Ｇ（ｔ）との双方が取得される。

なお、この実施形態では、全原子モデルの緩和弾性率Ｇ（ｔ）と、粗視化モデルの緩和弾性率Ｇ（ｔ）とが組み合わされる態様が例示されたが、これに限定されるわけではない。例えば、ユナイテッドアトムモデルの緩和弾性率Ｇ（ｔ）と、粗視化モデルの緩和弾性率Ｇ（ｔ）とが組み合わされることにより、大きな時間スケールの緩和弾性率Ｇ（ｔ）を、比較的短時間で計算されうる。

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

図３に示した処理手順に従って、ポリスチレンの分子鎖モデルが設定された。複数の分子鎖モデルが空間内に配置されることにより、高分子材料モデルが設定された（実施例、比較例）。

実施例では、図７に示した処理手順に従って、第１温度の緩和弾性率を計算する緩和弾性率計算工程が実施された。緩和弾性率計算工程では、図８及び図１２に示した処理手順に従って、第１温度よりも高い第２温度が設定された高分子材料モデルの緩和弾性率を計算する計算工程と、高分子材料の時間−温度換算則に基づいて、第２温度の高分子材料モデルの緩和弾性率から第１温度の高分子材料モデルの緩和弾性率に変換する工程とが実施された。そして、図９及び図１３に示されるように、第１温度の緩和弾性率の対数と、時間幅の対数との関係が求められた。

比較例では、従来のシミュレーション方法と同様に、第１温度が設定された高分子材料モデルに基づいて、緩和弾性率が計算された。そして、図１４に示されるように、第１温度の緩和弾性率の対数と、時間幅の対数との関係が求められた。なお、実施例及び比較例のポテンシャル等については、明細書中の記載のとおりであり、その他の共通仕様は次のとおりである。
高分子材料モデル：
第１温度：３００Ｋ（２７℃）
第２温度：５００Ｋ（２２７℃）
分子鎖モデルの本数：９２本
重合度：１０
緩和弾性率を求めるための分子動力学計算の終了時間Ｅｔ：１０ns
差Ａｄ：５ns
時間幅ｔの間隔：１０ps
空間の長さＬ１：７nm
シフトファクターＡｔ：０．００１

テストの結果、実施例のシミュレーション方法では、第２温度が設定された高分子材料モデルについて、終了時間Ｅｔまでに、緩和現象を発現させることができた。そして、実施例では、第２温度の高分子材料モデルの緩和弾性率に基づいて、第１温度の高分子材料モデルの緩和弾性率が求められた。

他方、比較例のシミュレーション方法では、第１温度が設定された分子材料モデルについて、終了時間Ｅｔまでに、緩和現象を発現させることができなかった。緩和現象を発現させるには、少なくとも実施例の１０倍〜１００倍の時間が必要である。従って、実施例のシミュレーション方法では、比較例のシミュレーション方法に比べて、高分子材料の緩和弾性率を、短時間で計算することができた。

１１ 高分子材料モデル
Ｇ（ｔ） 緩和弾性率

Claims (6)

1. コンピュータを用いて、高分子材料の緩和弾性率を計算するためのシミュレーション方法であって、
   前記コンピュータに、前記高分子材料の分子鎖をモデル化した分子鎖モデルを設定する工程と、
   前記コンピュータに、前記分子鎖モデルを予め定められた空間内に配置して高分子材料モデルを設定する工程と、
   前記コンピュータが、前記高分子材料モデルを用いて、予め定められた第１温度の前記緩和弾性率を計算する緩和弾性率計算工程とを含み、
   前記緩和弾性率計算工程は、前記高分子材料モデルに、前記第１温度よりも高い第２温度を設定する工程と、
   前記第２温度が設定された前記高分子材料モデルの緩和弾性率を計算する計算工程と、
   前記高分子材料の時間−温度換算則に基づいて、前記第２温度の前記高分子材料モデルの緩和弾性率から前記第１温度の前記高分子材料モデルの緩和弾性率に変換する工程とを含むことを特徴とする高分子材料のシミュレーション方法。
2. 前記計算工程は、前記高分子材料モデルに緩和現象が発現する前の前記緩和弾性率を取得する工程と、
   前記緩和現象が発現した後の前記緩和弾性率を取得する工程とを含む請求項１記載の高分子材料のシミュレーション方法。
3. 前記第２温度は、５００Ｋ以下である請求項１又は２記載の高分子材料のシミュレーション方法。
4. 前記高分子材料モデルは、少なくとも８本の前記分子鎖モデルを含む請求項１乃至３のいずれかに記載の高分子材料のシミュレーション方法。
5. 前記分子鎖モデルは、全原子モデル又はユナイテッドアトムモデルである請求項１乃至４のいずれかに記載の高分子材料のシミュレーション方法。
6. 前記高分子材料モデルの重合度は、５０以下である請求項５記載の高分子材料のシミュレーション方法。