分子動力学による理想ゴムの力学物性解析A Study on Mechanical Properties of Ideal Rubber Using Molecular Dynamics

hp140018 「京」産業利用（実証利用） K Industrial Use

分子動力学による理想ゴムの力学物性解析

A Study on Mechanical Properties of Ideal Rubber Using Molecular Dynamics

日野 理

Osamu Hino 東洋ゴム工業株式会社 TOYO TIRE & RUBBER CO., LTD.

要旨 完全な四分岐構造を持つ理想ゴムの力学物性を，分子動力学シミュレーションを用いて調べた． ゴムを構成する分子鎖としては，Kremer-Grest のビーズスプリングモデル[1]を使用した．架橋点 間ビーズ数を変化させていくつかの理想ゴムモデルを作成し，単軸引張り，純粋せん断および等 二軸引張りの分子動力学シミュレーションを行った．シミュレーション結果より各モデルのモジ ュラスを算出し網目鎖密度と比較し，両者が比例関係にあることがわかった．さらに，引張特性 の温度依存性を計算した結果，理想ゴムモデルにおいては張力が高伸長比までエントロピー弾性 によって支配されていることがわかった．さらに，Green-Kubo 公式を用いて各種ゴムモデルの動 的粘弾性を計算した． キーワード：分子動力学，理想ゴム，引張特性，粘弾性，架橋構造 Abstract

Mechanical properties of the ideal rubber models with perfect tetra networks have been investigated using molecular dynamics simulations. The rubber models are composed of the coarse-grained bead-spring chain models of Kremer and Grest [1]. A number of ideal rubber models with various number of beads between the cross-linking points are constructed. The stress-strain relations of the ideal rubber models are obtained from uniaxial elongation, pure shear, and equibiaxial extension deformation simulations by molecular dynamics. The moduli for the ideal rubber models are estimated and it is found that they are proportional to their nominal cross-link densities. It is also found that the entropy contributions to the stresses are dominant over wide range of deformations. Finally, the viscoelasticities are calculated for the various rubber models using the Green-Kubo formula.

Keywords: Molecular Dynamics, Ideal Rubber, Tensile Property, Viscoelasticity, Crosslinking Network

© 2018 Research Organization for Information Science and Technology All rights reserved. Received: 6 October 2017

Accepted: 20 March 2018 Available online: 12 April 2018

1. 研究の背景と目的 ゴム材料の架橋構造とその力学物性の関係を明らかにすることは，ゴム研究の基本的テーマの 一つである．近年では中性子および放射光による散乱実験や原子間力顕微鏡によるゴム（ゲル） の架橋構造観察が盛んに行われ，力学物性との関連性が議論されている．こうした研究のうち特 に興味深いのは，Tetra-PEG ゲルの架橋構造と力学物性に関する研究[2], [3]である．この研究によ れば，架橋構造の均一性すなわち架橋点の空間分布均一性と部分鎖の単分散性がゴム（ゲル）の 良好な力学物性を実現するための重要な鍵になるとされている． 本研究では，完全な四分岐構造を持つ架橋点によって結合された単一分散部分鎖からなる粗視 化理想ゴムモデル（理想ゴムモデル）を用いて，それらの引張特性と動的粘弾性を分子動力学シ ミュレーションに基づいて計算する．本研究の目標は，上述のモデルが Tetra-PEG ゲルや“教科 書的な”ゴムの力学物性を再現できるかどうかを確認することである．ゴム材料に関する分子動 力学シミュレーションは数多く行われているが，単軸引張りにおける力学特性以外にも焦点を当 てた研究はほとんど存在しない．この意味で，本研究はビーズスプリング鎖によって構成される ネットワークを採用する分子動力学シミュレーションのゴム材料物性研究への適用可能性の是 非を問う重要な試行と言える．また，Tetra-PEG ゲルの動的粘弾性測定によれば，それは通常の ゲルに比べて低損失性において極めて優れていることが報告されている．こうした性質は，転が り抵抗をできるだけ小さく抑えたいタイヤ用ゴム材料設計の観点からも興味深い． そこで本研究では，Green-Kubo 公式を用いて理想ゴムモデルの動的粘弾性を計算し，実際に低 損失性が見出されるかどうか，およびそれがどのような構造的特徴に関係しているかを，分子動 力学シミュレーションを用いて解析した． 2. 計算モデル 本研究においては，Kremer-Grest のビーズスプリングモデルに従って，長さ，質量およびエネ ルギーに関する代表量, m,  を用いて物理量を無次元化して議論する（補足説明）．理想ゴムモデ ルは，以下の手順で作成した．①最初に立方対称性を持つ四分岐架橋点分布を空間群 Fd3m（ダ イヤモンド構造）に従って生成し，②それらの間を特定の個数（10, 20, 30, 40 および 50 個）のビ ーズで結合したものを初期構造とし，③この後，温度 1.0，圧力 4.85 の条件下で平衡化して作成 した．このプロセスの概略を図 1 に示す． また理想ゴムモデルとの比較のため，以下の手順でも架橋ゴムモデルを作成した．①理想ゴム モデルの総部分鎖ビーズ数に（ほぼ）一致するように 200 個のビーズからなる分子鎖で構成され た高分子メルトモデルを作成し，温度 1.0，圧力 4.85 で平衡化する．②理想モデルにおける架橋 点数と同じ個数のビーズ（これを架橋剤ビーズと呼ぶ）を高分子メルトにランダムに添加し，全 ての架橋剤ビーズをその近傍の分子鎖ビーズと結合させて四分岐構造を生成する．③再度，温度 1.0，圧力 4.85 で平衡化する．以上のような手順で作成されたモデルを，実際の加硫工程によっ て得られたゴムをモデル化したものという意味で，加硫ゴムモデルと呼ぶ．

図 1 理想ゴムモデル作成過程（左がダイヤモンド構造の格子点間に直線状にビーズを配置した初期 状態．右は，それを温度 1.0，圧力 4.85 で平衡化した後の状態を示している．） 3. 並列計算の方法と効果（性能） 分子動力学計算は LAMMPS[4]を使用して行った．LAMMPS では OpenMP によるハイブリッド 並列計算も可能だが，主に使用した 48 ノード，384 コアレベルではフラット MPI による計算が 最も効率的であった．ベンチマークとして，使用した最大のモデル（100 万粒子系）において 384 ノードまでを使用した並列計算性能を測定した結果，アムダール則による評価法において 99.88% の並列化率を得た． 4. 研究成果 4.1. 引張特性計算 2 章に示した手順で作成した理想ゴムモデルを用いて，単軸引張り，純粋せん断および等二軸 引張りシミュレーションを行い，そこからひずみと応力に関するデータを抽出し，引張特性を計 算した．これらの引張特性データを 3 項からなる Ogden 式[5]：





3 1 3 n n n n x y z n n W       



  

(1)

をひずみエネルギー密度関数とする超弾性体モデルによって，3 種類の変形モードに関して，













3 3 1 0.5 1 1 1 1 1 3 1 2 1 1 , n n n n n n SE n SE SE PS n PS PS n n EB n EB EB n f f f                                 







(2)

という応力ひずみ関係を用いてフィッティングした．式(1)において，W はひずみエネルギー密 度， , ,  _{x y z}

はそれぞれ，

x y z, ,

軸方向の伸長比，

 _n, _n



n1, 2,3



は Ogden 式に用いられる数

値パラメータである．式(2)

において，fSE,fPS, fEB

はそれぞれ，単軸引張り，純粋せん断，等二

軸引張り変形における公称応力，

  SE, PS, EB

はそれぞれ，単軸引張り，純粋せん断，等二軸

引張り変形における伸長比を表す．

図 2 は，20 個のビーズからなる部分鎖で構成された理想ゴ ムモデルにおける 3 種類の引張特性と式

(2)の

Ogden 式によるフィッティング結果を合わせて示 したものである．図 2 から，一組のパラメータセットで異なる変形モードにおける引張特性をよ く再現できることが分かる．これは，他の部分鎖長を持つ理想ゴムモデルの引張特性においても 共通である．各モデルの部分鎖の個数はその作成時に決まっているので，それを平衡化して得ら れた体積で割ることにより，厳密な網目鎖密度が計算される．一方，各理想ゴムモデルの引張特 性は，低ひずみ領域においていずれも理想鎖のそれによってよく近似できる． 図 2 ビーズ 20 から構成された部分鎖を持つ理想ゴムモデルの引張特性とその Ogden 式によるフィッ ティング結果（いずれも縦軸が公称応力，横軸が伸長比である．上左：単軸引張り，上右：純粋せん 断，下：等二軸引張りにおける引張特性を示している．青□：計算結果，赤○：フィッティング結果） 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 1 2 3 4 5 6 公 称応⼒ ／εσ -3 伸⻑⽐ Ogden(単軸引張り) MD（単軸引張り） 0 0.004 0.008 0.012 0.016 0.02 1 1.1 1.2 1.3 1.4 1.5 公称応⼒ ／ε σ -3 伸⻑⽐ Ogden(等⼆軸引張り) MD（等⼆軸引張り） 0 0.01 0.02 0.03 0.04 0.05 0.06 1 1.5 2 2.5 3 3.5 4 公称 応⼒／εσ -3 伸⻑⽐ Ogden(純粋せん断) MD（純粋せん断）

そこで，各モデルについて得られた Ogden 式より得られる次式： 2 3 0 1 1 2 1 1 lim 2 2n n n W G    _       



(3)

によって各々のモジュラスを評価する．古典的ゴム弾性理論によれば，網目鎖密度とモジュラス は比例関係にある．ここで作成した理想ゴムモデルにおいても，図 3 に示すように，これらはほ ぼ比例していると言える．その比例定数は約 0.25 と算出された． 図 3 理想ゴムモデルの網目鎖密度に対する引張計算から得られたモジュラス さらに，温度 0.9 および 1.1 における単軸引張りシミュレーションを各理想ゴムモデル対して実 施し，引張特性を計算した．これらの結果を次式：



1.1,





0.9,



0.2 L f T L f T L f T      _{ } _   

(4)

に代入して温度 1.0 における張力の温度微分を求めた．L は，引張り方向のモデルの長さである． 熱力学関係式： T T L A U f f T L L T          _{ } _{ }  _{ }         

(5)

に基づいて温度 1.0 における引張特性をエネルギー弾性成分（式(5)右辺第 1 項）とエントロピー 弾性成分（式(5)右辺第 2 項）に分離した．A とU は，それぞれモデルのヘルムホルツ自由エネル ギーと内部エネルギーを表す．図 4 にビーズ数 10, 20 および 50 の部分鎖の理想ゴムモデルにお ける結果をプロットしたものを示す． 図 4 で極めて特徴的なことは，理想ゴムモデルの単軸引張り変形に対する応答応力が，伸長比 6.0 という大きなひずみにいたるまで，部分鎖長によらずほぼ完全にエントロピー弾性によって 与えられることである．この性質と網目鎖密度とモジュラスの比例関係を併せ，理想ゴムモデル は非常に“ゴムらしい”力学物性を持つことを，分子動力学計算によって示すことができた． y = 0.25097 x R² = 0.99845 0 0.005 0.01 0.015 0.02 0.025 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 モジュ ラ ス 網⽬鎖密度

図 4 エネルギー弾性およびエントロピー弾性成分に分離した理想ゴムモデルの単軸引張りにおける 引張特性（上左：10 個，上右：20 個，下：50 個のビーズからなる部分鎖を持つモデルによる結果で ある．赤○：全応力，青○：エントロピー弾性成分，緑○：エネルギー弾性成分） 4.2. 動的粘弾性計算 理想ゴムモデルと 2 章で示した手順によって作成された加硫ゴムモデルを用いた分子動力学シ ミュレーションから，Green-Kubo 公式：

 

1

  



3 i j ij ij V G t P P t T    



 (6) に従って応力緩和を計算した．ここで，

 

_{, , ,}

 

_{, , , ,} ij G t V T P t i j x y zは，それぞれ応力緩和，体積， 温度および系の時間によって変動するせん断応力成分を表す．また，L は



による時間平均を 表している．動的粘弾性は，G t の片側フーリエ変換を用いて，

 

 

 





* 0 0 exp G  G i G t i t dt   

_

 (7) -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 1 2 3 4 5 6 公称 応⼒／ε σ -3 伸⻑⽐ ビーズ数10 エントロピー成分 エネルギー成分 -0.002 0.002 0.006 0.01 0.014 0.018 0.022 0.026 1 2 3 4 5 6 公称応⼒／εσ -3 伸⻑⽐ ビーズ数50 エントロピー成分 エネルギー成分 -0.01 0.01 0.03 0.05 0.07 0.09 0.11 1 2 3 4 5 6 公称応⼒／ εσ -3 伸⻑⽐ ビーズ数20 エントロピー成分 エネルギー成分

として与えられる． * 0, , G G  は，それぞれ平衡弾性率，複素弾性率および周波数である．力学損 失を表す指標としてよく用いられる，正接損失は，

 

**

 

_{ }

_{ }

 

Im '' tan Re ' G G G G         (8) で定義される．図 5 は，20 個のビーズからなる部分鎖で構成された理想ゴムモデルと，それを同 じ架橋点数，総部分鎖ビーズ数を持つ加硫ゴムモデルを用いた分子動力学シミュレーション結果 を用いて計算された動的粘弾性を示したものである． 図 5 理想ゴムモデルと加硫ゴムモデルによる分子動力学シミュレーションより得られた動的粘弾性 （左：理想ゴムモデル，右：加硫ゴムモデルの動的粘弾性を表す．横軸が周波数，左縦軸が弾性率， 右縦軸が正接損失で両対数表示している．青○：貯蔵弾性率，赤○：損失弾性率，緑○：正接損失） 部分鎖数を一致させるようにしたので，平衡弾性率は二つのモデルにおいて同等になったが， 異なる動的粘弾性を得た．まず，周波数 10-5_～10-7_{の低周波領域で，理想ゴムモデルの正接損失} は加硫ゴムモデルのそれに比べて非常に小さくなった．これは，理想ゴムモデルで速やかに応力 緩和することを意味する．加硫ゴムモデルにおいて，応力緩和の遅れはダングリング鎖と部分鎖 同士の絡み合いに由来することを，シミュレーション結果を解析することにより示すことができ た．これは，理想ネットワーク構造を持つ Tetra-PEG ゲルが非常に小さな正接損失を持つことに 対応する．一方，周波数 10-2_～10-4_{の中間的な周波領域では逆に理想ゴムモデルの方が大きな正} 接損失を持つと計算された（対応する実験報告はない）．これは，理想ゴムモデルで部分鎖長が 単分散であること，ならびに，それとは対照的に加硫ゴムモデルでは，部分鎖長分布はブロード であり，それが正接損失の周波数分散に反映されていることがシミュレーション結果の解析から わかった．これらの知見は，タイヤの低転がり抵抗と制動性能の両立を目指すタイヤメーカーと しては非常に興味深い． 0.0 0.3 0.6 0.9 1.2 1.5 1.E‐06 1.E‐05 1.E‐04 1.E‐03 1.E‐02 1.E‐01 1.E+00

1.E-07 1.E-06 1.E-05 1.E-04 1.E-03 1.E-02 1.E-01

ta nδ モジ ュ ラ ス ／ εσ -3 周波数/τ-1 G' G'' tanδ 10 10-1 10-2 10-3 10-4 10-5 10-6 10-7 10-6 10-5 10-4 10-3 10-2 10-1 1 モジ ュ ラ ス ／ εσ -3 0.0 0.3 0.6 0.9 1.2 1.5 1.E‐06 1.E‐05 1.E‐04 1.E‐03 1.E‐02 1.E‐01 1.E+00

1.E-07 1.E-06 1.E-05 1.E-04 1.E-03 1.E-02 1.E-01

ta nδ モジ ュ ラ ス ／ εσ -3 周波数/τ-1 G' G'' tanδ 10 10-1 10-2 10-3 10-4 10-5 10-6 10-7 ₁₀-6 ₁₀-5 ₁₀-4 ₁₀-3 ₁₀-2 ₁₀-1 0.0 0.3 0.6 0.9 1.2 1.5 1.E‐05 1.E‐04 1.E‐03 1.E‐02 1.E‐01 1.E+00

1.E-07 1.E-06 1.E-05 1.E-04 1.E-03 1.E-02 1.E-01

ta n δ G' G'' tanδ 10 10-1 10-2 10-3 10-4 10-5 10-6 10-7 ₁₀-6 ₁₀-5 ₁₀-4 ₁₀-3 ₁₀-2 ₁₀-1 0.0

4.3. 「京」で計算することの効用 本研究目的においては，大規模モデルによるシミュレーションではなく，小～中規模の計算モ デルを用いた多数のシミュレーションが必要だった．この点，京コンピュータの巨大な資源を細 分化して使用することで，効率的に計算データを得ることができた． 5. まとめと今後の課題 完全な四分岐構造を持つ理想ゴムモデルの引張特性および動的粘弾性を，分子動力学シミュレ ーションを用いて計算した．引張特性計算により，作成された理想ゴムモデルは古典的なゴム理 論から導かれる性質をよく再現していることがわかった．また，動的粘弾性計算により，理想ゴ ムが低損失性を持つばかりでなく，高い制動性能を持つゴムの候補となり得るという結果を得た． 今後は，本研究で使用した定性的ゴムモデルではなく，実在のゴムをモデル化する方法を目指 した研究開発を行う予定である． 補足説明 分子動力学シミュレーションでは，このような無次元化が標準的に行われる．本研究において は長さ，質量およびエネルギーに関する代表量, m,  を用いて物理量を無次元化して議論する． これによれば，時間の代表量は， 1 2 1 2 m



_



 で与えられる．温度はエネルギーと同じ次元を持つとする（従ってボルツマン定数は 1 となる）． また，応力，モジュラス，応力緩和およびひずみエネルギー密度関数の代表量は，いずれも， m2-2-3 = -3 となる．Ogden 式に現れるパラメータについては，_n



n1, 2,3



を無次元，_n



n1, 2,3



は応力 と同じ次元を持つとする．Kremer-Grest のビーズスプリングモデルでは，m1, 1,  と設定1 されている． 参考文献

[1] K. Kremer, G. S. Grest, J. Chem. Phys. 92 (1990) 5057-5086.

[2] T. Sakai, Y. Akagi, T. Matsunaga, M. Kurakazu, U. I. Chung, M. Shibayama, Macromolecular rapid communications, 31 (2010) 1954-1959.

[3] Y. Akagi, T. Katashima, Y. Katsumoto, K. Fujii, T. Matsunaga, U. I. Chung, T. Sakai, Macromolecules, 44 (2011) 5817-5821.

[4] S. J. Plimpton, J Comp Phys, 117, 1-19 (1995).