修士論文

繰り返し変形を受ける
ポリブタジエンのヒステリシス発現メカニズム：
分子動力学法による検討

指導教員： 冨田佳宏

平成 20 年 2 月
神戸大学大学院自然科学研究科
機械工学専攻062T375N
古賀 裕

On Hysteresis
of Polybutadiene under Cyclic Deformations :Molecular Dynamics Study

Feburary 2008
Department of Mechanical Engineering,
Graduate School of Science and Technology,
Kobe University, Kobe, Japan
Yutaka KOGA

要約
ゴム材は，内部の三次元網目構造に依存して複雑な変形挙動を示す．本研究では，ゴムのヒステリシスを生じるメカニズムについて分子鎖レベルの知見を得るために，ポリブタジエンアモルファスブロックへの繰り返し変形シミュレーションを行った．既報のポリエチレン (PE) の場合と異なり，PB では十分緩和計算を行っても引張応力が bondstretch(結合長) に残留した．bond stretch の残留応力について，隣接する結合長と，さらにその先の結合長まで含めた 6 つのメチン，メチレン基からなる立体配座で整理した
結果，cis-1,4 の立体配座のモノマーにおいて，両端の=CH-CH2 -の二面角が anticlinalとなるノードの bond stretch に圧縮応力が生じていることを見出した．我々はこの構
造を “compression ノード”と名づけた．繰り返し変形を与えると，bond stretch および van der Waals(非共有結合) にヒステリシスを生じることは PE と同じであったが，bending には引張り-除荷で可逆的に解消する圧縮応力を生じること，この可逆圧縮応力により引張後期では系の応力−ひずみ応答が下に凸になる「ひずみ軟化」を生じることを示した．さらに bond stretch に生じるヒステリシスについて，compression ノードに属する bond stretch が圧縮側で大きくヒステリシスを生じているのに対し，それ以外のノードでは変形中ほぼ一定の引張り応力のままでありヒステリシスには寄与しない，等を明らかにした．また compression ノードを多く有する分子鎖は糸まり状の形態を取り，それらが凝集した部分は「からみ点」となる一方，それらの間の分子鎖は配向しながら大きく延伸してネットワーク構造を形成することを明らかにした．次に平均分子鎖長を倍としたアモルファスブロックで緩和計算ならびに繰り返し変形シミュレーションを行ったが，初期残留応力，繰り返し変形時の応力-ひずみ曲線等はほとんど変わらなかった．さらに，引張りによる局所密度の変化を調べ，PE の場合と異なり引張変形下でも高密度側のピークは不変であり，糸まり状の部分がほとんど解消
されないことを示した．

Summary
Rubber materials show complicated deformation behavior depending on the internal network structure．In this study，cyclic deformation is applied on an amorphous
nano-block of cis-1,4 polybutadiene (PB) by molecular dynamics simulation. Contrary to polyethylene (PE) previously reported, the PB block originates the residual
tensile stress on the bond stretches of chains, despite of the suﬃcient relaxation of
initial structure. We have investigated the residual stress of bond stretch separately,
by considering the possible three-dimensional conformation of successive six methineor methylene groups. It reveals that the bonds belong to the cis-1,4 conformation,of which both ends has the conformation of “anticlinal” dihedral angle of CH2 -CH
center, feel compressive stress while all the other bonds do tensile stress. We have referred to the configuration as “compression node”. As same as the PE case，the bondstretch and van der Waals dominate the hysteresis of the stress-strain curve of the PB
block; however, we have also find that compressive stress emerges on bending angles
during the tension, leading subtle “strain softening” in the stress-strain response. This
negative stress reversibly vanishes during the unloading process. On the hysteresis ofbond stretches, only the compression nodes show hysteresis while the other does almost constant tensile stress during the cyclic deformation. We have then visualized
the change in the chain structure under the cyclic loading, revealing that the compression nodes don’t distribute uniformly but aggregate. The structure show network-likemorphology under large elongation, since the aggregated compression nodes don’t fan
out but the other nodes are straightened and oriented between them. We have then,doubled the average chain length in the amorphous block and performed same simulations．However，there is no remarkable diﬀerence，between long and short chainsin the initial residual stress，stress-strain response，and so on．We have also investigated the change in the local density in the block under the tension．As the PE case，
the distribution of local density branches into two peaks by the formation of network
structure，however, the magnitude of the high density peak doesn’t change in the PB，
while that of PE shifts lower side. This support the observed fact that the aggregatedchains don’t dissolve.

目次
第1章 緒 論1

第 2 章 解析手法の基礎

42.1

分子動力学法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

2.2

原子間ポテンシャル. . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

2.3

高速化手法

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14

第 3 章 初期残留応力とその局所構造

16

3.1

シミュレーション条件 . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

3.2

シミュレーション結果および考察 . . . . . . . . . . . . . . . . . . . . .

19

3.2.1

緩和計算中のエネルギーおよび応力変化 . . . . . . . . . . . . .

19

3.2.2

bond stretch の応力とコンフォメーション . . . . . . . . . . . .

22

3.2.3

初期圧縮応力を生じる局所構造 . . . . . . . . . . . . . . . . . .

24

3.2.4

compression ノードの多い分子鎖構造 . . . . . . . . . . . . . . .

26

結言 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

3.3

第 4 章 繰り返し変形シミュレーション

28

4.1

シミュレーション条件 . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

4.2

シミュレーション結果および考察 . . . . . . . . . . . . . . . . . . . . .

30

4.3

4.2.1

応力-ひずみ関係. . . . . . . . . . . . . . . . . . . . . . . . . .

30

4.2.2

ポテンシャル成分毎のヒステリシスへの寄与 . . . . . . . . . . .

32

4.2.3

bond stretch のヒステリシス発現メカニズム . . . . . . . . . . .

34

4.2.4

compression ノードによるネットワーク構造 . . . . . . . . . . .

35

結言 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

38

第 5 章 分子量の違いによる変化

39

i

目次

ii

5.1

シミュレーション条件 . . . . . . . . . . . . . . . . . . . . . . . . . . .

39

5.2

シミュレーション結果及び考察 . . . . . . . . . . . . . . . . . . . . . .

40

5.2.1

緩和計算中のエネルギーおよび応力変化 . . . . . . . . . . . . .

40

5.2.2

bond stretch の応力とコンフォメーション . . . . . . . . . . . .

42

5.2.3

compression ノードの少ない分子鎖形態 . . . . . . . . . . . . . .

43

5.2.4

繰り返し変形時における応力-ひずみ関係 . . . . . . . . . . . . .

44

5.2.5

局所密度およびネットワーク構造 . . . . . . . . . . . . . . . . .

46

結言 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

48

5.3

第6章 結 論

49

参考文献

52

学術論文・学術講演

56

謝辞

58

第1章
緒論
ゴム材は精密機器に用いられる微細なものから，防舷材や橋梁の免震材といった巨大なものにまで用いられる，現代社会には欠かせない重要な工業材料である．ゴムは三次元の網目構造を持つ高分子であり，履歴現象を伴う Mullins 効果 [1] や Payne 効果[2]といったゴム特有の力学特性は複雑な一次構造や高次構造に起因する [3] ．履歴現象(ヒステリシス) とは，変形負荷時の応力と除荷時の応力経路が異なる現象である．代
表的なゴム製品である自動車用タイヤにとって，ヒステリシスによる損失 (ヒステリシ
スロス) を制御することにより，転がり抵抗の低減とグリップ性能の向上という相反する機能を両立させることが可能となる．しかしながら，タイヤ設計にはこれまで開発者の勘と経験に依存する部分が大きく，希望通りの性能を備えた製品を得ることが難しかった．そこで近年，コンピュータシミュレーションを援用した設計が盛んに行われるようになり，多くの成果を残している [4] ．
ゴムのヒステリシスを表現するために，分子鎖網目理論 [5] に基づき分子鎖の架橋点変化を考慮した非アフィン分子鎖網目モデル [6] や，網目の各短鎖に管模型 [7] に基づく
粘性抵抗 [8] を導入し粘弾性を考慮したモデル [9] が提案され，実験結果を良好に再現できることが知られている．管模型を分子鎖に適用した reptation(爬行) 理論 [10],[11] は，高分子の緩和挙動を取り扱うレオロジーの分野で現在でも活発に拡張されており，鎖長の揺らぎ・鎖の伸縮・管の拘束の緩和など，様々な効果を考慮できるモデル [12]−[15]や，異なる鎖に属する二つの絡み合い点のカップリングを考える sliplink 模型 [16] 等が提案されている．

1

第1章

緒論

2

しかしながら，これらのモデルのキーとなる「からみ点」やその間の「セグメント」[17]

は，現象論的に導入されたものであるため，それよりも「微視的な」メカニズムを追求することはできない．実際の複雑な高分子鎖の凝集構造においてこれらの「からみ点」やその間の「セグメント」が一体どのような局所構造に対応するのか，その物理的描像の解明には個々のモノマーの立体配座（コンフォメーション）を考慮できる
分子スケールのシミュレーションが必要である．
近年の計算機能力の飛躍的向上を背景に，ナノ〜サブミクロンスケールにおける材料の変形挙動を直接追求できる手法として分子動力学法が注目されており，高分子材料への適用も盛んに試みられている．一本の分子鎖の応力-ひずみ挙動 [18]–[21] や，結晶化についての報告 [22] の他，ガラス転移温度についての研究 [23], [24] などが数多くなされている．一方，屋代ら [25], [26] は，ゴムのヒステリシスについて分子レベルからの知見を得るべく，ポリエチレン（PE）／ポリブタジエン（PB）の分子動力学シミュレー
ションを行ってきた．PE の繰り返し変形シミュレーション [25] では，系の応力への寄与を，結合長，結合角，二面角そして非共有結合（van der Waals）に分けて評価し，引張時と除荷時の応力−ひずみ曲線の違い，すなわちヒステリシスは結合長と van der

Waals が担っていること，除荷時の trans → gauche 遷移により，繰り返し変形後には「たるんだ」分子鎖が増加し，再引張り時には初回よりも低い応力を示すこと，等を明らかにしている．一方，PB 中の二重結合の役割を明確にすべく，PE と PB を比較したシミュレーション [26] では，初期残留応力の違いがあるものの，両者のヒステリシス
ループには大きな違いがなく，「巨視的な」応答は原子間の相互作用に鈍感であると報告している．しかしながら，引張変形後期にわずかに応力が減少し応力−ひずみ曲線
がわずかに下に凸になる「ひずみ軟化」などの明確な違いも見出されている．
本研究では，PB の繰り返し変形シミュレーションを行い，その内部メカニズムについてより詳細に検討する．PE の解析では，二面角が１種類しかなく分子鎖構造変化の考察が容易であったが，PB の場合は 3 種類の二面角が存在し，それぞれの立体配座の組み合わせによって様々な分子鎖構造をとる．結合長，結合角等の状態をそのような分子鎖構造毎に分類して議論することで複雑な PB の内部メカニズムを明らかにする．
第 2 章では解析手法の基礎として，分子動力学法を簡単に説明し，分子動力学計算

第1章

緒論

3

で最も重要となるポテンシャルエネルギーについて述べる．また，大規模計算を行うための高速化手法および応力の評価方法を示す．第 3 章では，一辺 20nm の cis-1,4 ポ
リブタジエンのアモルファスブロックを作成して初期緩和シミュレーションを行い，応力変化と立体配座の観点から初期残留応力を生じる分子鎖構造について検討する．第4 章では，繰り返し引張変形シミュレーションを行い，ヒステリシスを生じるメカニ
ズムを，初期残留応力を生じる局所構造に着目して明らかにする．第 5 章では，分子鎖長を倍としたアモルファスブロックについて同様のシミュレーションを行い，分子鎖長の変化がヒステリシスに与える寄与を議論する．最後に，第 6 章で本研究の総括
を述べる．

第2章
解析手法の基礎2.1

分子動力学法

分子動力学法 (molecular dynamics method，略して MD 法) は，系を構成する各粒子についてニュートンの運動方程式

mid2 r i
= Fi
dt2

(2.1)

を作成し，これを数値積分することにより粒子の軌跡を求める方法である．ここで，mi ，r i はそれぞれ粒子 i の質量および位置ベクトルである．粒子 i に作用する力 F i は，系のポテンシャルエネルギー Etot の各位置における空間勾配として次式により求められる．

Fi = −

∂Etot
∂r i

(2.2)

式 (2.1) の数値積分には，Verlet の方法，予測子–修正子法等がよく用いられる [33] ．
本研究では，以下に示す Verlet の方法を用いた．
時刻 t + ∆t と t − ∆t での粒子 i の位置ベクトル r i (t ± ∆t) を taylor 展開すると，()
dr i (t) (∆t)2 d2 r i (t)3r i (t + ∆t) = r i (t) + ∆t
+
+
O
(∆t)
dt
2
dt2
2 2
(
)
dr i (t) (∆t) d r i (t)3r i (t − ∆t) = r i (t) − ∆t
+
+
O
(∆t)
dt
2
dt2

(2.3)
(2.4)

となる．ここで，v i を時刻 t における粒子 i の速度とすると，

dr i= v i (t)
dt
4

(2.5)

第2章

解析手法の基礎

5であり，式 (2.1) と式 (2.5) を式 (2.3) と式 (2.4) に代入すると，()(∆t)2 F i (t)r i (t + ∆t) = r i (t) + ∆tv i (t) ++ O (∆t)3
2
mi
2
(
)
(∆t) F i (t)r i (t − ∆t) = r i (t) − ∆tv i (t) ++ O (∆t)3
2
mi

(2.6)
(2.7)

となる．両式の和と差をとると，

r i (t + ∆t) + r i (t − ∆t) = 2r i (t) + (∆t)2()
F i (t)+ O (∆t)4mi(

r i (t + ∆t) − r i (t − ∆t) = 2∆tv i (t) + O (∆t)3

)

(2.8)
(2.9)

が得られる．これより，時刻 t + ∆t での位置ベクトルと t での速度は

r i (t + ∆t) = 2r i (t) − r i (t − ∆t) + (∆t)2v i (t) =()
F i (t)+ O (∆t)4
mi()
1{r i (t + ∆t) − r i (t − ∆t)} + O (∆t)22∆t(2.10)(2.11)と求められる．t + ∆t での座標を求めるには 2 つの時刻 t と t − ∆t での座標が必要である．初期の計算 (t = 0) では t = ∆t での座標 r i (∆t) は式 (2.6) と初速度から求められる．r i (∆t) と r i (0) が既知であれば，式 (2.10) を繰り返し適用することにより各粒子の座標を求められる．

2.2

原子間ポテンシャル

粒子に作用する力は系のポテンシャルエネルギーにより決定される．したがって，系のポテンシャルエネルギーの評価が分子シミュレーションにおいて重要となる．Lennard–

Jones ポテンシャル [27] や Morse ポテンシャル [28] などの経験的 2 体ポテンシャルは取り扱いが簡単であるため，従来からよく用いられてきたが，金属材料表面や異相界面など粒子配置が急激に変化する現象を扱う近年の研究では，多体力を考慮したより精度の高いポテンシャル関数 [29], [30] が用いられている．
本研究で解析対象とするポリブタジエンはメチレン炭素 (-CH2 -) とメチン炭素 (-CH-)から構成されるが，水素原子は陽に扱わず，水素原子の効果を繰り込んで CH，CH2

第2章

解析手法の基礎6をそれぞれ一粒子として扱う united atom model[31] を用いて分子鎖の挙動を表現する．
疎視化した粒子に対する分子鎖内結合および分子鎖間 van der Waals 相互作用は，Gee
ら [32] に提案されているものを用いた．

bond stretch ポテンシャル（2 体間ポテンシャル）ΦBS (r) =∑1{

kr (r − r0 )22}(2.12)

bending ポテンシャル（3 体間ポテンシャル）ΦBE (θ) =∑1{2kθ (θ − θ0 )2}(2.13)

torsion ポテンシャル（4 体間ポテンシャル）ΦTO (ϕ) =∑12{k1 (1 + cos ϕ) + k2 (1 − cos 2ϕ) + k3 (1 + cos 3ϕ) + k4 (1 − cos 3ϕ)}(2.14)

van der Waals ポテンシャル（2 体間ポテンシャル）ΦVW (r̄) =

∑

{

σ
σ
4ε ( )12 − ( )6
r̄
r̄

}

(2.15)

各関数のパラメーターの値を表 2.1〜2.5 に，ポテンシャル曲線を図 2.1〜2.4 にそれぞれ示した．なお，式中の総和は対象とする高分子集合体中の全ノードについて行うが，van der Waals ポテンシャルの計算は異なる分子鎖の粒子間，4 原子団以上離れた粒子
間，C 原子と-CH2 -粒子間に対して行う．
式 (2.2) による力の評価において，2 体間相互作用の bond stretch および van der Waals
は式 (2.16)〜(2.17)，3 体間の bending は式 (2.18)〜(2.20)，4 体間の torsion は式 (2.21)〜(2.24) を用いて評価される．式中の記号はそれぞれ図 2.5〜2.7 に示した．

第2章

解析手法の基礎⃝ bond stretch ポテンシャルΦBS (r) =∑1{2

kr (r − r0 )2}

Table 2.1 Potential Parameter for bond stretch.

CH2 – CH2
CH2 – CH
CH = CHΦ ΒS [eV]

150

r0
[nm]
0.154
0.150
0.133

kr
[kJ/ (mol · nm2 )]2.650 × 103
3.075 × 103
4.130 × 103

CH2 -CH2
CH2 -CH
CH=CH

100

r

50

0

0.1

0.2

r

0.3

0.4

[nm]

Fig.2.1 Relationship between bond stretch potential ΦBS and bond length r.

7

第2章

解析手法の基礎⃝ bending ポテンシャル

ΦBE (θ) =∑1{2kθ (θ − θ0 )2}

Table 2.2 Potential Parameter for bending.
θ0[deg.]C – CH2 – C 111.6
C – CH = C 124.0

Φ BE [eV]150kθ[kJ/( mol・rad2 )]482
374-CH2 -CH=

100

θ

50

0

90

180

270

360

θ [deg.]Fig.2.2 Relationship between bending potential ΦBE and bending angle θ.

8

第2章

解析手法の基礎⃝ torsion ポテンシャル

CH2 -CH2
ΦTO (ϕ) =∑12{k1 (1 + cos ϕ) + k3 (1 + cos 3ϕ)}

CH=CH
ΦTO (ϕ) =∑12

k2 (1 − cos 2ϕ)

CH2 -CHΦTO (ϕ) =∑12

k4 (1 − cos 3ϕ)

Table 2.3 Potential Parameter for torsion.
k1
k2
k3
k4
[kJ/mol] [kJ/mol] [kJ/mol] [kJ/mol]
3.35
58.6
13.4
7.95

Table 2.4 Stable and metastable torsion angles of CH=CH, CH-CH2 , CH2 -CH2
center.

CH2 -CH2
CH2 -CH
CH=CH

Stable180◦ (trans)
0◦
(cis)
◦
0
(cis)

Metastable
67.5◦(gauche)◦120(anticlinal)◦
180(trans)9

第2章

CH2 -CH
CH=CH
CH2 -CH2

0.6

Φ TO [eV]

解析手法の基礎

0.4

0.2

φ
0
0

90

180

270

360

Torsion angle, θ , [deg.]

Fig.2.3 Relationship between torsion potential ΦTO and torsion angle ϕ.

10

第2章

解析手法の基礎⃝ van der Waals ポテンシャル

ΦVW (r̄) =

∑

{

σ
σ
4ε ( )12 − ( )6
r̄
r̄

}

Table 2.5 Potential Parameter for van der Waals.

CH2 – CH2
CH = CHε[kJ/mol]
0.69
0.50

σ
[nm]
0.35
0.33

0.3

-CH2-CH=

Φ VW [eV]

0.2

r-

0.1

r-

0

-0.1

0.2

0.4

-r [nm] 0.6

0.8

Fig.2.4 Relationship between van der Waals potential ΦVW and internuclear
separation r̄

11

第2章

解析手法の基礎12

2 体間相互作用r = r αβFα = −

Φ′ (r) αβrr(2.16)

F β = −F α(2.17)

3 体間相互作用
r αγ · r βγ
|r αγ | |r βγ |{}Φ′ (θ)
1r αγ · r βγ
α
βγ
r −
F =
r αγ
sin θ |r αγ | |r βγ ||r αγ |2cos θ ={1Φ′ (θ)r αγ · r βγ βγ
αγ
F =
r
−
r
sin θ |r αγ | |r βγ ||r βγ |2

(2.18)}βF γ = −F α − F β(2.19)(2.20)

4 体間相互作用
r αγ · r γηA =r −r γη
|r γη |2r βη · r γηr γηAβ = r βη −|r γη |2
Aα · Aβcos ϕ =
|Aα | Aβααγ

1Φ′ (ϕ)F =α
sin ϕ |A | Aβ
α

Fβ =

′

{

Aα · Aβ αA
A −|Aα |2

}

β





(2.21)




Φ (ϕ)
1
A ·A β
Aα −
2 A 
α
β

sin ϕ |A | A 

Aβ
α

β

r αγ · r γη α r βη · r γη βF −F|r γη |2|r γη |2r αγ · r γη α r βη · r γη βF η = −F β +F +F|r γη |2|r γη |2

F γ = −F α −(2.22)(2.23)(2.24)

第2章

解析手法の基礎

i
rij

j

ri
rj
Fig.2.5 Two molecules i, j and intermolecule vector r ij

i
j
rik

θ

rjk

k
Fig.2.6 Three molecules i, j, k and bendig angle θ

j
rjl

φ

k

Ai rik

Aj

rkl
l

i
Fig.2.7 Four molecules i, j, k, l and torsion angle ϕ

13

第2章

2.3

解析手法の基礎

14

高速化手法

粒子数 N の系において粒子間の全相互作用を評価すると，1step に N × (N − 1) 回の計算が必要となり，N が大きくなると極めて膨大な計算量となる．実際には，一定距離以上離れた粒子は影響を及ぼさないので，作用を及ぼす範囲 (カットオフ半径 rc )内の粒子からの寄与を効率よく計算することにより高速化できる．従来よく用いられてきた高速化手法に粒子登録法 [33] がある．これは，図 2.8 に示したように，rc よりひとまわり大きい半径 rfc 内の粒子をメモリーに記憶し，その中で rc 内の相互作用を評価する方法であり，N × (rc 内粒子数 ≪ N − 1) に計算負荷が減少される．しかし，粒子登録法では rfc 半径より外の粒子が rc 内に達すると力の評価が適切でなくなるので，一定のステップ毎に登録粒子の更新 (N × (N − 1) 回の探査) を行わなければならない．
このため，系がある程度の規模以上になると，粒子登録による高速化は登録更新の計
算負荷により打ち消される．

rfc

rc

Fig.2.8 Schematic of bookkeeping method.

第2章

解析手法の基礎

15

別の高速化手法としてブロック分割法がある．図 2.9(a) に示すように，シミュレートする系をカットオフ距離程度の格子状に分割し，各ブロックに属する粒子をメモリーに記憶する．着目している粒子に作用する力を評価する際には，その粒子が属するブロックおよび隣接するブロックから相互作用する粒子を探索して行う (図 2.9(b))．粒子が属するブロックは，粒子の位置座標をブロックの辺長 bx，by で除した際の整数により判断できるので，ブロック登録時の計算負荷は粒子数 N のオーダーとなる．したがって，粒子登録法では登録更新の負荷が大きくなるような大規模な系でも高速化が可能である．
ポリマーのポテンシャルでは，共有結合部の bond stretch，bending，torsion ポテン
シャルは相互作用する粒子が同一分子鎖内で予め決まっているため，原子対を探索す
る必要はなく分子鎖単位での並列化による高速化も容易である．一方，van der Waals
ポテンシャルは異分子鎖間，あるいは，同一分子鎖内の 4 粒子以上離れた全粒子に対して相互作用を評価する必要があり，本研究で扱うような大規模な系では，ブロック分割による高速化が必要となる．

y
by
x

0
(a) Domain decompositionbx(b) Block serchingFig.2.9 Schematic of domain decomposition method.

第3章
初期残留応力とその局所構造
3.1

シミュレーション条件

図 3.1 に示すように，20nm × 20nm × 20nm の立方体セル中に乱数を用いて多数のランダムコイル状分子鎖を成長させ，ポリブタジエンアモルファス構造を作成した．分
子鎖成長時には，図 3.2(a) に模式的に示すように，CH2 -CH2 ボンドを中心とする二面角は ϕ = 180◦ (trans 点) または ϕ = ±67.5◦ (gauche 点) のいずれか，CH2 -CH を中心とする二面角は ϕ = 0◦ (cis 点) または ϕ = ±120◦ (anticlinal 点) のいずれかにランダムに設定し，CH=CH を中心とする二面角は ϕ = 0◦ (cis 点) に固定して cis-1,4 ポリブタジエンの立体配座にする．ただし，実際に二面角をこれらの安定点にのみ限定すると，アモルファス構造が目的の密度に達するまで各分子鎖を成長させることが困難となるため，CH2 -CH2 ，CH2 -CH の二面角については trans，gauche，anticlinal から ±15◦ のずれを許容した．ここで，平均分子鎖長さを 300，標準偏差を 50 に設定し，分子鎖長の分布が正規分布となるように分子鎖成長をコントロールした．作成したポリブタジ
エンアモルファス構造の密度は 0.86g/cm3 ，総分子鎖数は 1000，総粒子数は 300,792 である．得られた初期構造に対し，全方向自由表面の境界条件のもとで 40ps(=400000fs)の分子動力学計算を行い内部構造緩和を行った．数値積分には Verlet 法を用い，積分の時間ステップは 0.1fs とした．温度は 300K とし，速度スケーリングにより制御した．

16

第3章

初期残留応力とその局所構造

nm

20nm

20

20nm

Fig.3.1 Simulation cell

17

第3章

trans

gauche
o
67.5

n=b×a

初期残留応力とその局所構造

180

o

CH2

b
a

gaucheCH(a)

Conformation of CH2-CH2 center120o
anticlinal
n=b×a
CH2
CH

b
cis

a

(b)

anticlinal
0o

Conformation of CH-CH2 center

n=b×a

b
CH
a

o

cis

0

CH2

CH
CH2

(c)

Conformation of CH=CH center

Fig.3.2 Conformation of CH2 -CH2 ，CH2 -CH, CH=CH center dihedral angles.

18

第3章

初期残留応力とその局所構造

3.2

シミュレーション結果および考察

3.2.1

緩和計算中のエネルギーおよび応力変化19

図 3.3(a),(b) に初期緩和計算におけるポテンシャルエネルギーとその成分 (bond stretch，

bending，torsion，van der Waals) の時間変化をそれぞれ示す．系全体のエネルギーは初期構造の極めて高いエネルギー状態から急激に低下した後，緩やかに減少して 40ps
ではほぼ 0 となっており，構造緩和が十分行われたことを示す．図 3.4 は対応する系の z 軸方向の応力変化を示したものである．応力は，図 3.5 に模式的に示すように，セル内の各原子位置において，正方向に作用する力 (切断法による内力) を加算することによりセルに働く x，y ，z 方向の軸力を求め，セルの体積で徐すことにより評価した．
初期配置における応力は約-4000MPa と極めて大きな圧縮応力であり，本ポテンシャ
ルに関しては設定した密度が高かったことが分かる．すなわち図 3.3(a) の高いエネル
ギーは初期配置の高密度に起因するものである．このため，自由境界の下で緩和計算を行うと急激に膨張して応力は正の値になり，20ps 以降には振動はしているもののその平均は約 40MPa でほとんど変化しなくなる．原子に作用する力はポテンシャルエネ
ルギーの空間微分から得られるので，式 2.12〜2.15 から各ポテンシャル成分の応力への寄与を別々に評価することができる．図 3.6 は σzz の変化に対し，各ポテンシャル成分の寄与分を示したものである．van der Waals と bending，ならびに torsion については，いずれもその値はほぼ 0 であるため，図 3.4 の残留応力は bond stretch によるものであることが分かる．

初期残留応力とその局所構造

10000

10000

0

0

0

10

20
30
Time, t, ps

040(a) Total

20

bond stretch
bending
torsionvan der Waals-10000-100001020
30Time, t, ps40(b) bond stretch, bending,torsion,and van der Waals

Fig.3.3 Change in the potential energy and its components of bond stretch，
bending，torsion，and van der Waals.

Stress, σzz , MPa

Potential energy, E, eV

第3章

50

0
0

10

20
30
Time, t, ps40

Fig.3.4 Change in the stress during the initial relaxation.

第3章

初期残留応力とその局所構造

P=ΣPzi /Lz

Pzi

j

Pz

Lz

k

Pz

Ly

Lx

Stress, σzz, MPa

Stress, σzz, MPa

Fig.3.5 Evaluation of stress on the simulation cell.

50

0

50

0
0

10

20

30

40

0

10

Time, t, ps

50

0
10

20

30

30

40

(b) Van der Waals

Stress, σzz, MPa

Stress, σzz, MPa(a) Bond stretch

020Time, t, ps

40

50

0
0

10

20

30

Time, t, ps

Time, t, ps(c) Bending(d) Torsion

40

Fig.3.6 Change in the stresses generated on bond stretch, bending, torsion,and van der Waals during the initial relaxation.

21

第3章

3.2.2

初期残留応力とその局所構造

22bond stretch の応力とコンフォメーション

残留応力がどの bond stretch で発生しているのかをさらに詳細に調べるため，CH=CH，

CH2 -CH2 ，CH-CH2 の 3 つの bond stretch を，さらにそれぞれが中心となるコンフォ
メーションにより分類して調べた．例えば CH=CH の bond stretch ならば，図 3.7 に模式的に示すように cis，trans のコンフォメーション毎に別々に応力を評価する．図

3.8 は bond stretch の応力を，それぞれのコンフォメーション毎に分けて評価し，緩和
計算中の変化を示したものである．シミュレーション中でのコンフォメーションの決定は，各安定点から± 30◦ 以内にある二面角を cis や trans として評価し，それ以外はother として分類した．大きな正の応力を生じるのは，図 3.8(a) の CH=CH における
trans ならびに other，図 3.8(b)CH-CH2 の other ならびに cis である．一方，緩和計
算中に強い圧縮応力を生じているのは CH=CH の cis ならびに CH-CH2 の anticlinalである．図 3.8(c) の CH2 -CH2 は図 3.8(a)，(b) に比べ生じている応力は小さいが全体的に正の値をとる．これらの和が図 3.6(a) の系全体の bond stretch で平均した変化で
ある．

r

180

bond stretch

o

trans

CH

0o

cis

CH2

Fig.3.7 Categorizing the bond stretch by its conformation.

第3章

初期残留応力とその局所構造

trans

other

20

20

, MPa

20

cis

0

0

0

Stress,

zz

other

cis

-20
0

10

23

anticlinal

-20

0

40

10

20

30

gauche
trans
other

-20

0

40

10

20

30

20

30

Time,

ps

Time, , ps

Time, , ps(a) CH=CH(b) CH-CH2(c) CH2-CH2

40

Fig.3.8 Change in the stress on each bond stretch categorized with its conformation (initial relaxation).

(anticlinal)-(cis)-(cis)(anticlinal)-(cis)-(other)(other)-(cis)-(other)(trans)-(anticlinal)-(gauche)(trans)-(anticlinal)-(other)

0

0

-20

-20

(cis)-(anticlinal)-(gauche)(cis)-(anticlinal)-(trans)(anticlinal)-(cis)-(anticlinal)

0

10

20

30

Time, t, ps(a) CH=CH

40

0

10

20

30

Time, t, ps(b) CH-CH2

40(anticlinal)-(trans)-(other)(anticlinal)-(gauche)-(other)

20
Stress, σzz , MPa

20
Stress, σzz , MPa

Stress, σzz , MPa

20

0

-20

(anticlinal)-(trans)-(anticlinal)(anticlinal)-(gauche)-(anticlinal)

0

10

20

30

Time, t, ps(c) CH2-CH2

Fig.3.9 Change in the stress on bond stretch separately evaluated by its combination of successive three conformations.

40

第3章

3.2.3

初期残留応力とその局所構造

24

初期圧縮応力を生じる局所構造

1 種類の二面角しか存在しないポリエチレンと異なり，ポリブタジエンの場合は隣接した結合の二面角も考慮しなければ分子鎖の局所構造を明らかにすることができない．図 3.8(a) ならびに図 3.8(b) において圧縮応力を生じた cis の CH=CH，ならびに anticlinal の CH-CH2 の bond stretch に生じる応力を，さらにその両端の結合が取るコンフォメーションを考慮し分類して評価したものが図 3.9(a)，(b) である．図3.9(a) の cis のコンフォメーションをとる CH=CH では，青線で示した (anticlinal)(cis)-(anticlinal) の組み合わせのみが圧縮応力を示している．このときの CH=CH のcis を中心とする局所構造を模式的に図 3.10- に示す．なお，図中で CH 粒子を黒色，
CH2 粒子を青色で，応力を評価する bond stretch を赤色で示した．

図 3.9(b) の anticlinal の CH-CH2 については，黒ならびに青線で示した (cis)-(anticlinal)-(trans) と (cis)-(anticlinal)-(gauche) の 2 つの組み合わせのときに圧縮応力を示すが，このときの局所構造は図 3.10-

′

ならびに

となり，

の構造と中心の 4 粒子 (CH=CH-CH2 -CH2 ) の構造が共通している．なお，図 3.8(c) の CH2 -CH2 の二面角のみ考慮した分類では圧縮応力を示さなかったが，CH2 -CH2 の trans，gaucheのコンフォメーションも両側の 2 つを考慮して評価すると，図 3.9(c) の黒ならびに青線で示したように，anticlinal に挟まれた trans と gauche が圧縮となった．CH2 -CH2 がgauche で両端が anticlinal の構造は図 3.10- に相当する．CH2 -CH2 が trans で両端が anticlinal の構造を図 3.10- に示す．これは図 3.10-′の構造と CH-CH2 -CH2 -CHの部分が同じである．
図 3.10- の右側で再び anticlinal となり，CH=CH と結合していることを考えると，

3.10- に関連した bond stretch が圧縮応力を生じていることが分かる．cis-1,4 のモノマー構造の両端の二面角が anticlinal であるこの構造を，以下では compression ノードと称することにする．なお，上記以外の全ての組み合わせについても同様に調べたが，他に圧縮応力を示すノードは存在しなかった．

第3章

l CH

na
cli

CH2

ti
an

CH

初期残留応力とその局所構造

an
ti

cli
na

cis

CH2

l

CH2

CH 2

CH

an
ti

CH

CH

cli
na

l

cis

an
ti

CH2 trans

cli
n

CH

CH

gauche

CH 2

CH2

al

cis

CH2

CH2

CH

CH2

a

C

trans an
CH2 tic
cli
lin
na
al
l

CH nti

CH2

CH

C

Fig.3.10 Schematic of local configurations that generate negative residual
stress on its bond.

25

第3章

初期残留応力とその局所構造

26

compression ノードの多い分子鎖構造

3.2.4

緩和計算終了時の系において，分子鎖内の粒子間の結合数をカウントし，その総数
で compression ノードに分類された結合数を除すことにより compression ノードの 1
分子鎖内に占める割合を調べた．その割合が高かった分子鎖と低かった分子鎖の例を
図 3.11(a)，(b) にそれぞれ示す．両者を把握し易くするため，図中左上に示すように

compression ノードは緑色と黄色で，それ以外の構造をとる部分は紫色と青色でそれぞれ CH 粒子と CH2 粒子を着色している．compression ノードが分子鎖内に多く含まれる場合には分子鎖は糸まり形状を取っており，そのため compression ノードの bondstretch に圧縮応力を生じたものと考えられる．一方 compression ノードが少ない場合には直線状の部分が増え，やや広がった形態となっている．

CH2
CH 2
CH

CH

CH2

CH2
CH

CH2
CH2

CH

CH2
CH2

compression-node(anticlinal)-(cis)-(anticlinal)(a) compression-node richother-node(b) less compression-node

Fig.3.11 Example of molecular chain structure.

第3章

3.3

初期残留応力とその局所構造

27

結言

cis-1,4 ポリブタジエンのアモルファス構造を作成して初期構造緩和シミュレーショ
ンを行い，得られた平衡状態における残留応力およびそれを生じる局所構造について
詳細に検討した．得られた結果を要約して以下に示す．

(1) 系のエネルギーがほぼ収束したにもかかわらず，約 40MPa の引張応力が残留した．ポテンシャル成分毎の応力評価を行った結果，引張り残留応力は bond stretchに生じていることが示された．

(2) CH=CH，CH2 -CH2 ，CH-CH2 が中心となるコンフォメーション毎に分類し，それぞれの bond stretch に生じる応力を調べた結果，圧縮を生じているノードがあることを明らかにした．

(3) 圧縮残留応力を生じる bond stretch を，隣接する bond stretch と，さらにその先の bond stretch まで含めた 6 つのメチン，メチレン基からなる立体配座で整理した結果，cis-1,4 の立体配座のモノマーにおいて，両端の=CH-CH2 -の二面角がanticlinal となるノードの bond stretch に圧縮応力が生じていることが明らかにした．

(4) 上記の cis-1,4 のモノマー構造の両端の二面角が anticlinal であるこの構造を，compression ノードと称する．この compression ノードを多く有する分子鎖は糸まり状となり，反対に少ない場合には分子鎖は広がった形態をとることを明らか
にした．

第4章
繰り返し変形シミュレーション
4.1

シミュレーション条件

前章で得た初期平衡状態のアモルファスブロックに対し，ひずみ制御で z 方向に繰り返し変形を与えるシミュレーションを行った．
図 4.1 に示すように立方体セルの上下端から 1nm 以内の粒子をつかみ部とし，その領域内にある粒子は自由に運動できるが境界を越えることはできない条件で，セルの上下方向 (z 軸方向) に最大ひずみ εzz =1.0 となるまで引っ張るシミュレーションを行った．積分の時間ステップは 1.0fs とし，毎ステップ ∆εzz = 1.0 × 10−5 のひずみ増分を与えてひずみを増加させた．ひずみ速度に換算すると ε̇zz = 1.0 × 1010 /s となる．なおセ
ルの上端つかみ部から下端まで連結した分子鎖は存在しない．横方向については自由境界としたが，壁面から 1nm 以内の粒子を登録し，その重心からアモルファスブロッ
クの壁の平均位置を算出することでセルの体積を評価した．ひずみは全粒子の z 座標
をスケーリング (アフィン変形) することで与えているが，ひずみ増分が小さく系の内部構造緩和が逐次生じるため，変形後期の内部構造は初期のそれをアフィン変形させたものとは異なる．ここで，ひずみは初期緩和計算後のセル長さを基準とする工学ひずみである．最大ひずみ εzz =1.0 に達した後，保持時間 0 でひずみを反転させ，同一ひずみ速度で除荷した．
この引張り−除荷を 1 サイクルとして，2 サイクルの繰り返し変形シミュレーションを行った．系の温度は 300K とし，速度スケーリングにより制御した．応力の評価につ

28

第4章

繰り返し変形シミュレーション

29

いては 3 章と同様に切断法を用い，セル上下端のつかみ部をのぞく可動部分の粒子に作用する内力を，可動部分の体積で除すことにより評価した．

chuck

nm

20nm

loading direction

20

z

y
x

20nm

Fig.4.1 Simulation cell and loading direction

第4章

繰り返し変形シミュレーション

4.2

シミュレーション結果および考察

4.2.1

応力-ひずみ関係

30

図 4.2 に系の z 軸方向の応力-ひずみ曲線を示す．黒線が 1 サイクル目，青色の線が 2
サイクル目の変形を表している．図ではひずみ εzz =0.01 間隔のゆらぎを生じているが，実際には毎ステップ振動している．ただし，毎ステップの応力振動を表示すると図が煩雑となり全体的な傾向が見えにくくなるので，1000 ステップごとの応力変化を示し
た．1 サイクル目の引張り開始時は，3.2.1 節で議論した bond stretch による約 40MPaの残留応力が存在する．引張りを開始するとひずみ約 0.1 までに応力は約 75MPa 程度まで上昇する．その後はほぼ一定応力となるが，引張りを続けているにもかかわらず

εzz = 0.4 近傍から応力はわずかに減少する．εzz = 0.6 近傍から再びわずかに上昇する
ため，結果として下に凸の曲線を描いている．εzz = 1.0 でひずみを反転すると，ごく短い時間 (〜2000fs) で約 65MPa から約 45MPa まで急減する．これはポリエチレンで報告した弾性変形に相当する [25] が，ポリエチレンに比べて減少量は小さく，また垂直に応力低下せずやや傾きを示している．その後，応力が 0 となる εzz = 0.3 まではほぼ線形に減少するが，圧縮側に入ると非線形に抵抗が増加している．ひずみの基準は無負荷平衡状態でのセル辺長としているので，εzz = 0.0 では圧縮応力約 55MPa となった．ポリエチレンの時と同様，引張り・除荷時の応力経路の違いにより大きなヒステ
リシスを生じている．2 サイクル目の応力-ひずみ曲線は 1 サイクル終了時の圧縮状態からの引張りとなるため，εzz = 0.1 までは 1 サイクル目よりもわずかに低い値を取る
が，その後は 1 サイクル目の応力-ひずみ曲線とほぼ一致する．除荷時は応力が 0 となるひずみが 1 サイクル目よりも大きいため，εzz < 0.3 以下のひずみでは 1 サイクル目
より圧縮応力が大きくなり除荷曲線は下側を通るが，それ以外に顕著な違いは見られ
なかった．

第4章

繰り返し変形シミュレーション

σzz , MPa

50

Stress,1000Loading

Unloading-50

1st cycle
2nd cycle
0

0.5

1

Strain , ε zz
Fig.4.2 Stress-strain curves under cyclic deformation.

31

第4章

4.2.2

繰り返し変形シミュレーション

32

ポテンシャル成分毎のヒステリシスへの寄与

図 4.2 の応力-ひずみ応答について，各ポテンシャル成分の寄与毎に評価し，それぞれ
図 4.3(a)〜(d) に示した．ポリエチレンと同様 [25] ，図 4.3(a)，(b) にそれぞれ示す bondstretch ならびに van der Waals がヒステリシスを主に担っていることが分かる．ただし，引張り過程での応力変化を詳細に見ると，bond stretch は εzz = 0.2 近傍でほぼ一定応力となった後，εzz = 0.6 近傍より再びわずかに応力は上昇している．van der

Waals はほとんど変化していない．すなわち 4.2.1 節で示した，ポリエチレンとの違いである変形後期のわずかな応力減少は bond stretch，van der Waals のいずれでもない．
注目すべきは，図 4.3(c) に示した bending である．引張りを与えているにもかかわらずその応力は減少し，引張り過程において常に負の値をとっている．この bending の
圧縮応力が 4.2.1 節で示した「ひずみ軟化」をもたらしている．またこの bending の圧
縮応力は引張り時と同じ経路をとって回復し，εzz = 0.0 では初期とほぼ同じ状態となる可逆的な挙動を示した．torsion については図 4.3(d) に示すように引張り・除荷時ともにその値はほぼ 0 であり，ヒステリシスへの寄与はない．2 サイクル目の変形では，

bond stretch ならびに van der Waals はいずれも 1 サイクル終了時には初期より低い応力となるため，再引張り初期は 1 サイクル目の応力-ひずみ曲線よりも低い値を取る
が，その後は 1 サイクル目とほぼ同じ経路を描く．除荷過程では，van der Waals は反転させた直後から 1 サイクル目よりもわずかに大きな圧縮応力を示し，除荷にともな
いその差が顕著となる．bond stretch は系の応力が 0 となるひずみ付近から 1 サイク
ル目よりも大きな圧縮応力を示し，応力-ひずみ曲線は下側を通っている．bending とtorsion には 1 サイクル目と 2 サイクル目に大きな違いはない．

Stress, σzz, MPa

第4章

100

100

50

50

0

0

-50

1st cycle
2nd cycle
00.5
Strain , ε zz

繰り返し変形シミュレーション

-50

1

0

(a) Bond stretch

1(b) Van der Waals

100

100
Stress, σzz, MPa

0.5Strain , ε zz

50

50

0

0

-50

-50

0

0.5
Strain , ε zz(c) Bending

1

0

0.5
Strain , ε zz

1(d) Torsion

Fig.4.3 Change in the stresses generated on bond stretch, bending, torsion,and van der Waals under cyclic deformation.

33

第4章

4.2.3

繰り返し変形シミュレーション

34

bond stretch のヒステリシス発現メカニズム

前章で示した compression ノードに着目し，それに属する bond stretch に生じる応力，およびそれ以外の bond stretch に生じる応力を分けて評価したものが図 4.4 である．図には 1 サイクル目の変化のみ示している．図 3.10- に示した compression ノー
ドでないコンフォメーションをとる CH=CH と CH2 -CH2 の bond stretch は初期引張応力のままほとんど変化せず一定であり，除荷過程で圧縮となるときに僅かに低い応力側にずれるだけでヒステリシスへの寄与は小さい．compression ノードに属さない

CH-CH2 の bond stretch(図中青線) は，引張りを与えると初期引張応力から 5MPa 程度の上昇を示す．しかしすぐにほぼ一定で変化しなくなる．除荷過程ではこの応力増加分だけ応力-ひずみ経路が下側を通るため，先の CH=CH，CH2 -CH2 よりヒステリシス
が大きくなる．最も顕著なヒステリシスを示すのは図中紫色の線で示した compression
ノードである．初期圧縮応力は引張りに対して “バッファ”となり，compression ノー
ドは引張過程では緩やかではあるが与えたひずみに応じて応力が回復し，上に凸の曲線となる．一方，除荷初期の応答は，引張りのそれを反転させた形で，やや下に凸の曲線となる．εzz < 0.5 で初期応力より圧縮側になると変形抵抗が急激に上昇するため曲線は上に凸 (圧縮ひずみ，応力を正とすると下に凸の 2 次曲線状) になる．このように bond stretch のヒステリシスは主として compression ノードの圧縮応力において発現する．

σzz , MPa

50

Stress,

100

0

CH=CH
CH-CH2
CH2-CH2
Compression-node

-50

0

0.5

1

Strain , ε zz

Fig.4.4 Change in the stress of bond stretches in the compression and theother nodes.

第4章

4.2.4

繰り返し変形シミュレーション

35

compression ノードによるネットワーク構造

1 サイクル目の変形における分子鎖構造の変化を図 4.5 に示す．図 (a) 上に模式的に示したように，中心付近の厚み 2.0nm の領域の分子鎖のみ抽出して示している．3.2.4
節と同様に compression ノードを緑色に，それ以外のノードは青色で着色している．図4.5(a) の引張り前において，compression ノードは均一に分布しているのではなく，所々に凝集している様子が認められる．系が引張りを受け伸張すると，compression ノードが凝集した緑色の部分は分子鎖が密となったままで，青色の他のコンフォメーションの部分が大きく延伸してネットワーク構造を形成している (図 4.5(b)〜(f))．前章の図

3.11(a) で示したように，分子鎖内に compression ノードが多くなると個々の分子鎖は糸まり状の形状をとるが，図 4.5 中に黒の実線の楕円で囲った部分ではこれら糸まり状の分子鎖が複数凝集して互いに絡み合い，延伸しにくくなったからみ点となっているものと考えられる．除荷過程では，図 4.5(f)〜(k) に示したように，からみ点間の分子鎖が疎になった部分が消滅するが，緑色の compression ノードが凝集した部分の間にはやはり青色のその他のコンフォメーションの分子鎖が再び存在している．図 4.6 に 2
サイクル目の分子鎖構造変化を示す．図中の黒線の楕円は 2 サイクル目の compression
ノードの位置から決めたものではなく，そのまま重ねて 1 サイクル目の変形における絡み点の位置 (図 4.5) を示したものである．2 サイクル目の変形においても絡み点は消失することなく，そのため延伸時のネットワーク構造も 1 サイクル目のそれと同様となっている．
前節で議論した様に，凝集してからみ点を形成している compression ノードは初期圧縮応力を示すことから，その bond stretch は平衡長さよりも短くなっており，引張り変形に対する “あそび ”を有している．そのため引張り変形時に，コンフォメーショ
ンを崩すことなく bond stretch の伸張 (回復) により変形を吸収することができる．一方，compression ノード以外の分子鎖構造をとるノードでは，初期状態で bond stretch
が正の引張り応力を示していることから，bond stretch はもともと延伸されており，それ以上伸張させることができない．そのため引張り変形下ではこれらのノードはからみ点間で二面角の回転により直線化し，ネットワーク状の構造を形成したと考えられ

第4章

繰り返し変形シミュレーション

36

る．また，4.2.2 節で示した bending の可逆圧縮応力は，丸まった形状を変えないまま，
からみ点の中で結合長が回復した結果と解釈できる．

Loadingzxy(a) ε zz = 0.0(b) ε zz = 0.2(c) ε zz = 0.4(d) ε zz = 0.6(e) ε zz = 0.8(f) ε zz = 1.0

Unoading(f) ε zz = 1.0(g) ε zz = 0.8(h) ε zz = 0.6(i) ε zz = 0.4(j) ε zz = 0.2(k) ε zz = 0.0

Fig.4.5 Snapshots of molecular chains during the 1st cyclic deformation.
Green molecules indicate the “compression node”, or the cis-1,4 conformation with anticlinal dihedral angles of both ends =CH-CH2 -，while blue ones the other conformation.

第4章

繰り返し変形シミュレーション

37

Loading(a) ε zz = 0.0(b) ε zz = 0.2(c) ε zz = 0.4(d) ε zz = 0.6(e) ε zz = 0.8(f) ε zz = 1.0

Unoading(f) ε zz = 1.0(g) ε zz = 0.8(h) ε zz = 0.6(i) ε zz = 0.4(j) ε zz = 0.2(k) ε zz = 0.0

Fig.4.6 Snapshots of molecular chains during the 2nd cyclic deformation.

第4章

4.3

繰り返し変形シミュレーション

38

結言

ヒステリシスをもたらすゴムの分子鎖変形挙動について，原子レベルから知見を得ることを目的としてポリブタジエンの繰り返し引張変形シミュレーションを行った．得
られた結果を要約して以下に示す．

(1) 繰り返し変形を与えると，ポリエチレンの場合と同様に bond stretch および vander Waals にヒステリシスを生じたが，bending に引張過程で負の応力が生じ，除荷時に戻るという可逆的な現象が認められた．この圧縮応力は，引張後期の系の応力−ひずみ応答が下に凸になる「ひずみ軟化」をもたらす．

(2) compression ノードを構成する bond stretch 以外は，引張変形時にほとんど応力
上昇せず，ヒステリシスへの寄与は小さい．一方，compression ノードに属する

bond stretch は，引張によって圧縮応力が大きく回復するとともに，除荷時には圧縮側の経路をとって大きなヒステリシスを示した．

(3) 繰り返し変形時の分子鎖構造変化を観察した結果，compression ノードは均一に分布せず凝集しており，引張時は compression ノード間の他の構造の分子鎖が大きく延伸され，compression ノードの凝集した部分をつなぐようなネットワーク
状の様相を呈した．

(4) compression ノードを多く有する糸まり状の分子鎖が凝集した部分は互いに運動しにくい「からみ点」となる．このからみ点は，bond stretch に負の初期応力を有しているため，引張変形に対するバッファとなり，繰り返し変形時のヒステリシスを担う．一方，糸まり状の構造を保ったまま bond stretch の圧縮が回復する
ために bending に (1）の圧縮応力を生じる，という分子鎖の変形メカニズムを明らかにした．

第5章
分子量の違いによる変化
5.1

シミュレーション条件

ポリブタジエン分子鎖の平均分子鎖長を変化させ，第 3 章，第 4 章と同様の緩和計
算ならびに繰り返し変形シミュレーションを行った．分子鎖の平均分子鎖長 (分子鎖の粒子数) が前章までの系の倍の 600，標準偏差が 50 の正規分布となるように制御した．
得られたポリブタジエンアモルファス構造の密度は 0.84g/cm3 ，総分子鎖数は 496，総粒子数は 298144 である．作成した初期構造に対し，前章までと同様の条件で緩和，繰
り返し変形シミュレーションを行った．なお，以降では便宜上平均分子鎖長が 600 のポリブタジエンアモルファスブロックを model(L)，第 3，4 章で述べた平均分子鎖長が 300 のそれを model(S) と称する．

39

第5章

分子量の違いによる変化

5.2

シミュレーション結果及び考察

5.2.1

緩和計算中のエネルギーおよび応力変化40

図 5.1(a),(b) に初期緩和計算におけるポテンシャルエネルギーと z 軸方向の応力の時間変化をそれぞれ示す．図には比較のため第 3 章の model(S) の解析結果を青色で表示している．ポテンシャルエネルギー，応力の変化ともに model(L)，model(S) の間に大きな差は見られず，約 40MPa の引張り応力が残留した．図 5.2 に各ポテンシャル成分の応力への寄与分を示す．model(S) と同様，bending，torsion，van der Waals には

3000

0

Model(L)Model(S)

-3000
0

10

20
30
Time, t, ps(a) Potential energy

40

Stress, σzz , MPa

Potential energy, E, eV

応力が残留しておらず bond stretch にのみ応力が生じている．

50

Model(L)Model(S)

0
0

10

20

30

40

Time, t, ps(b) StressFig.5.1 Change in the potential energy and stress during the initial relaxation(long vs. short chains).

50

Model(L)Model(S)Stress, σzz, MPa

Stress, σzz, MPa

第5章

0分子量の違いによる変化

50

0

0

10

20

30

40

0

10

Time, t, ps

50

0
0

10

20

30

40

(b) Van der Waals

Stress, σzz, MPa

Stress, σzz, MPa(a) Bond stretch

20
30
Time, t, ps

40

50

0
0

10

20

30

Time, t, ps

Time, t, ps(c) Bending(d) Torsion

40

Fig.5.2 Change in the stresses generated on bond stretch, bending, torsion,and van der Waals during the initial relaxation (long vs. shortchains).

41第5章

分子量の違いによる変化

42

bond stretch の応力とコンフォメーション

5.2.2

図 5.2(a) の bond stretch に生じる応力について，第 3 章と同様に，CH2 -CH2 -[CH=CH]-

CH2 -CH2 と CH-CH2 -[CH2 -CH]=CH-CH2 ，ならびに CH=CH-[CH2 -CH2 ]-CH=CH の連続する 6 粒子がとるコンフォメーションの組み合わせによる分類を行った．全ての組み合わせについて調べたが，図が煩雑になるため系への寄与が大きいもののみを図5.3 に示す．図には合わせて model(S) のそれを破線および 〜 の数字で表示している．図に示したいずれの組み合わせにおいても，model(S) と有意な差はみられず，局所構造とその bond stretch に生じる応力については分子鎖長の違いによる変化はない．

(anticlinal)-(cis)-(anticlinal)(anticlinal)-(trans)-(other)(anticlinal)-(other)-(other)Stress, σzz , MPa

20(cis)-(anticlinal)-(gauche)(cis)-(anticlinal)-(trans)(gauche)-(cis)-(other)20(anticlinal)-(gauche)-(anticlinal)(anticlinal)-(gauche)-(other)(anticlinal)-(trans)-(anticlinal)

20

3
3

0

2

0

2

0

3

1

1
2

-20

-20

1 (anticlinal)-(cis)-(anticlinal)
2 (anticlinal)-(trans)-(other)3 (anticlinal)-(other)-(other)

0

10

20

30

40

1
2
3

0

1

-20

(cis)-(anticlinal)-(gauche)(cis)-(anticlinal)-(trans)(gauche)-(cis)-(other)

10

20

30

40

1 (anticlinal)-(gauche)-(anticlinal)
2 (anticlinal)-(gauche)-(other)3 (anticlinal)-(trans)-(anticlinal)

0

10

20

30

Time, t, ps

Time, t, ps

Time, t, ps(a) CH=CH(b) CH-CH2(c) CH2-CH2

Fig.5.3 Change in the stress on bond stretch separately evaluated by its combination of successive three conformations. Solid lines indicate thestress on model(L), while the dashed ones that on model(S).

40

第5章

5.2.3

分子量の違いによる変化

43

compression ノードの少ない分子鎖形態

compression ノードが 1 分子鎖内に占める割合が低い分子鎖の例を図 5.4 に示す．3.2.4
節と同様に，compression ノードを緑色に，それ以外の構造をとる部分は青色で着色している．3 章の model(S) と比べ，分子鎖内に compression ノードノードによる糸まり状の部分を有するようになり，model(S) のときのように糸まり状の部分がほとんどない分子鎖は存在しなかった．

Fig.5.4 Example of molecular chain involving less compression-node (long chain).

第5章

分子量の違いによる変化

44

繰り返し変形時における応力-ひずみ関係

5.2.4

図 5.5(a)，(b) に 1 サイクル目と 2 サイクル目の変形における系の z 軸方向の応力-ひずみ曲線をそれぞれ示す．図には model(S) の結果を青線で合わせて表示している．1
サイクル目の引張り時には，εzz = 0.4 までの範囲において model(S) より平均して約5MPa 高い応力を生じた．また除荷時の εzz < 0.3 の圧縮変形領域でも model(S) より
圧縮応力が大きくなっているが，その他では大きな違いは見られない．2 サイクル目の応力-ひずみ曲線は εzz < 0.25 の圧縮変形領域を除いては model(S) とほぼ一致した．
図 5.5(a) の 1 サイクル目の応力-ひずみ応答について，第 4 章と同様に各ポテンシャ
ル成分の寄与毎に評価しそれぞれ図 5.6(a)〜(d) に示す．model(S) と比較して，引張り過程の εzz = 0.1〜0.5 において図 5.6(a) の bond stretch は引張応力が，除荷過程のεzz < 0.3 の範囲においての図 5.6(b) の van der Waals は圧縮応力がそれぞれ増大している．図 5.6(c)(d) の bending ならびに torsion は model(S) からの変化は見られない．

100

σzz , MPa50

Loading

50

0Stress,0

Stress,σzz , MPa

100

Unloading

-50

-50

Model(L)Model(S)

-100
0

0.5

Strain , ε zz(a) 1st cycle

1Model(L)Model(S)

-100
0

0.5

1

Strain , ε zz(b) 2nd cycle

Fig.5.5 Stress-strain curves under cyclic deformation (long vs. short chains).

第5章

100

Stress, σzz, MPa

Stress, σzz, MPa100

分子量の違いによる変化

50
0
-50

50
0
-50

Model(L)Model(S)

0

0.5

1

0

Strain , ε zz

1

Strain , ε zz(a) Bond stretch(b) Van der Waals

100

Stress, σzz, MPa

100

Stress, σzz, MPa

0.5

50
0

50
0
-50

-50

0

0.5

1

0

0.5

Strain , ε zz

Strain , ε zz(c) Bending(d) Torsion

1

Fig.5.6 Stress-strain curves under cyclic deformation (long vs. short chains).

45

第5章

5.2.5

分子量の違いによる変化

46

局所密度およびネットワーク構造

1 サイクル目の引張り過程において，引張り前 (εzz =0.0)，途中 (εzz = 0.5)，最大ひずみ点 (εzz =1.0) でのアモルファスブロックの局所密度の分布を図 5.7(a)〜(c) に示す．
局所密度は，各粒子を中心に van der Waals エネルギーのカットオフ距離 0.8nm を半径とする球内にある近接粒子数を調べ，粒子質量の合計を球体積で除すことにより評価している．図 5.7(a) の引張り前の状態において，model(L) では model(S) よりも高密度な粒子が多いことが分かる．引張り変形下の分布図 5.7(b)，(c) では差がほとんどなくなっているものの，注意深く見るとピークは僅かではあるが高密度側にあり，またピーク左側の面積が減少している．また，εzz = 1.0 の図 (c) では，ポリエチレンで報告した [25] ように，ネットワーク構造の出現による密度ピークの分岐を生じ始めており，ρα < 0.5 の「疎な」粒子の分布には model(L) と model(S) の間に差はない．また，高密度側のピークの位置がほとんど変化していない (ρα =1.0g/cm3 近傍) 点は，ポリエチレンの場合 [25] と異なり注目すべき点である．図 5.8(a) に 1 サイクル目の変形時における分子鎖構造の変化を示す．図には引張り過程における構造を，前章と同様に

compression ノードを緑色に，それ以外のノードは青色で着色して示している．また，比較のため model(S) のものを図 (b) に示した．model(S) に比べ，compression ノードが凝集した部分とほかの部分がより分離しているように見える (εzz =0.0 の図)．しかしながら，変形下の凝集した部分と疎な部分のネットワーク構造には大きな違いはない．
これは，先の局所密度分布で示したように model(L) と model(S) の差がわずかである
ことからも理解できる．

第5章

Number of node, N

8000

分子量の違いによる変化

8000

8000

6000

6000

6000

4000

4000

4000

2000

2000

2000

Model(L)Model(S)

0

0
0

0.5

1

Local density, ρ α, g/cm3

1.5

47

0
0

0.5

1

1.5

Local density, ρ α, g/cm3(a) ε zz = 0.0

0

0.5

1

1.5

Local density, ρ α, g/cm3(c) ε zz = 1.0(b) ε zz = 0.5

Fig.5.7 Change in the distribution of the local density under tension.

zxy

ε zz = 0.0

ε zz = 0.2

ε zz = 0.4

ε zz = 0.6

ε zz = 0.8

ε zz = 1.0

ε zz = 0.6

ε zz = 0.8

ε zz = 1.0(a) Model (L)ε zz = 0.0

ε zz = 0.2

ε zz = 0.4(b) Model (S)Fig.5.8 Snapshots of molecular chains during the 1st tension (long vs. short chain).

第5章

5.3

分子量の違いによる変化

48

結言

分子量の違いによって生じる力学応答および分子鎖構造を調べるために，平均分子鎖長を前章までの倍としたアモルファスブロック (model(L)) の初期緩和計算ならびに
繰り返し変形シミュレーションを行った．得られた結果を以下に要約して示す．

(1) 分子鎖長を倍とした系 (model(L)) においても，初期緩和時のポテンシャルエネル
ギー，応力の変化に違いはなく，前章までの系 (model(S)) とほぼ同じ約 40MPaの引張応力が残留した．

(2) 分子鎖が長くなると，分子鎖内に compression ノードによる糸まり状の部分を有するようになり，model(S) のときのように糸まり状の部分がほとんどない分子鎖が存在しなくなる．

(3) 繰り返し変形下における応力-ひずみ曲線は，model(L) と model(S) でほとんど同じだったが，引張り時には bond stretch の引張り応力が，除荷時には van der

Waals の圧縮応力がそれぞれ model(S) よりもわずかに増大した．
(4) 局所密度の変化を調べた結果，引張り前の状態で model(L) は model(S) よりも高密度な粒子が多く存在した．また，ポリエチレンのときと同様，引張下での分子鎖ネットワーク構造の出現による局所密度の分岐を生じ始めていたが，高密度側のピークの位置はほとんど変わらない．このことは，糸まり状になった部分が変形によっても解消していないことを示している．

第6章
結論
本研究では，ヒステリシスをもたらすゴム内部の微視的分子鎖変形メカニズムを分子レベルで解明するために，分子動力学法によりポリブタジエンアモルファスブロッ
クへの繰り返し変形シミュレーションを行った．また，分子鎖長の異なる 2 つのアモ
ルファスブロックを比較し，分子鎖長の変化についても言及した．以下に，得られた結果を総括する．
第 2 章では，本研究で用いた解析手法の基礎について述べた．まず，分子動力学法の概要ならびに基礎方程式を示し，本研究で用いた数値積分法について説明した．次に，粒子間相互作用の評価に用いられるポテンシャルエネルギーについて述べ，ポリブタジエンのポテンシャル関数を具体的に説明した．さらに，大規模シミュレーショ
ンに必要な計算の高速化手法について述べた．
第 3 章では，cis-1,4 ポリブタジエンのアモルファス構造を作成して初期構造緩和シ
ミュレーションを行い，得られた平衡状態における残留応力およびそれを生じる局所構造について詳細に検討した．系のエネルギーがほぼ収束した初期平衡状態において
約 40MPa の引張応力が残留し，それが bond stretch によるものであることを明らかにした．さらに，CH2 -CH2 -[CH=CH]-CH2 -CH2 と CH-CH2 -[CH2 -CH]=CH-CH2 ，ならびに CH=CH-[CH2 -CH2 ]-CH=CH の 6 粒子がとるコンフォメーションで bond stretch に生じる残留応力を整理すると，cis-1,4 の立体配座をとるモノマーで，かつ両端の=CH-

CH2 -の二面角が anticlinal となるノードの結合長に圧縮応力が，それ以外の結合長には引張り応力が生じていることを明らかにした．この構造を”compression ノード”と名づ

49

第6章

結論

50

け，アモルファスブロック中の分子鎖の空間的な広がりを調べたところ，compression
ノードを多く有する分子鎖は糸まり状となり，反対に少ない場合には分子鎖は広がった形態をとることを明らかにした．
第 4 章では，ヒステリシスをもたらすゴムの分子鎖変形メカニズムを明らかにする
ため，第 3 章で作成したポリブタジエンアモルファスブロックに繰り返し引張変形を与えるシミュレーションを行った．応力-ひずみ曲線をポテンシャル成分毎に分離して評
価し，(1)bond stretch および van der Waals による応力にヒステリシスを生じる，(2)ポリエチレンと異なり，bending は除荷すると可逆的に解消する圧縮応力を引張時に生じる，(3)(2) の可逆圧縮応力により引張後期では応力−ひずみ応答が下に凸になる「ひずみ軟化」を生じる，等を明らかにした．さらに bond stretch に生じるヒステリシスについて，compression ノードに属するものとそれ以外に分けて評価し，(4)compression
ノードは負の範囲で大きくヒステリシスを描く，(5) それ以外のノードは変形中一定の引張り応力を示しヒステリシスには寄与しない，ことを示した．また繰り返し変形時の分子鎖構造変化を直接観察し，compression ノードは均一に分布せず凝集し，引張り時も密な構造を伴ったままであること，凝集した compression ノード間の分子鎖が配向しながら大きく延伸してネットワーク構造を形成すること，を明らかにした．以上を総合すると，(i)compression ノードを多く有する糸まり状の分子鎖が凝集した部分
は互いに運動しにくい「からみ点」となる，(ii) からみ点は，bond stretch に負の初期応力を有しているため，引張変形に対するバッファとなり，繰り返し変形時のヒステリシスを担う，(iii) 糸まり状の構造を保ったまま bond stretch の圧縮が回復するために bending に圧縮応力を生じる，という分子鎖の変形メカニズムが明らかになった．
第 5 章では，分子量の違いによって生じる分子鎖変形挙動への影響を調べるために，
平均分子鎖長を前章までの倍としたアモルファスブロックの初期緩和計算ならびに繰
り返し変形シミュレーションを行った．緩和計算時，ならびに繰り返し変形における応力の変化は前章までの短い分子鎖のそれとほとんど変わらなかった．ただし，引張り時には bond stretch の引張り応力が，除荷時には van der Waals の圧縮応力がそれ
ぞれ model(S) よりもわずかに増大していることを示した．また，分子鎖が長くなると糸まり状の部分をほとんど有しない線状の分子鎖が存在しなくなり，引張り初期には

第6章

結論

51

compression ノードが凝集した部分と他の部分がよりはっきりと分離した構造となった．さらに，引張りによる局所密度の変化を調べ，分子鎖が長い系の方が初期平衡状態において高密度な粒子が多く存在すること，引張変形下で凝集した部分とそれをつなぐ疎な部分のネットワーク構造を生じても，ポリエチレンの場合と異なり高密度側のピークは不変であり，糸まり状の部分がほとんど解消されないこと，等を明らかにした．
エントロピー弾性が本質であるゴム材料の変形において，本研究で示した現象およびメカニズムは，時間スケールの大きな隔たりがあることは否めないが，分子鎖変形
挙動およびからみ点形成メカニズムの物理的描像を考える際の一助になるものと考える．また，実際のゴムは物理的架橋であるからみ合いに加え，化学的架橋あるいはカーボンブラックやシリカといったフィラーを加えることによってその力学的特性が大きく変化することが知られている．本研究で得た知見をもとにスチレン-ブタジエン等の，より実際に近いゴム分子鎖での検討や，フィラーとゴムの作用などさらなる検討が期
待される．

参考文献
[1] Mullins, L., “ Thixotropic Behavior of Carbon Black in Rubber ”，Rubber Chemistry and Technology，12, (1987), 79-91.
[2] Payne, A. R., Whittaker, R. E., “ Low strain dynamic properties of filled rubbers”
Rubber Chemistry and Technology, 44, (1971), 440-470.
[3] Ferry, J. D., Viscoelastic Properties of Polymers, (1980), Wiley[4] 白石正貴，林公博，“ 自動車走行シミュレーションのための FEM タイヤモデル”，
日本ゴム協会誌，72-4, (2001), 223-251.

[5] 冨田佳宏，“ ガラス状ポリマーの分子鎖網目理論による構成式と変形挙動のシミュ
レーション”，塑性と加工，37-424, (2001), 485-491.

[6] Tomita, Y., Adachi, T. and Tanaka, S., “ Modelling and application of constitutive equation for glassy polymer based on nonaﬃne network theory”, European
journal of mechanics. A. Solids, 16, (1997), 745-755.
[7] Gennes, P. G., “ Reptation of a Polymer Chain in the Presence of Fixed Obstacles”, Journal of Chemical Physics, 55, (1971), 572-579.
[8] Bergstrom, J. S. and Boyce, M. C., “ Constitutive modeling of the large strain
time-dependent behavior of elastomers ”, Journal of the Mechanics and Physics
of Solids, 46-5, (1998), 931-954.
[9] 東圭佑，内藤正登，冨田佳宏, “ ゴムのひずみ速度依存性構成式とカーボンブラッ
ク充填ゴムの変形挙動シミュレーション”, 日本機械学会論文集, A73-729, (2007),

558-566.
[10] 土井正男, 小貴明, 高分子物理・相転移ダイナミクス, (1992), 岩波
52

参考文献

53

[11] Doi, M. and Edwards, S. F., The theory of polymer dynamics, (1986), Oxford
University Press[12] Cloizeaux, J. des, “ Double Reptation vs. Simple Reptation in Polymer Melts”,
Europhysics Letter, 5, (1988), 437-442.
[13] Marruchi, G., “ Dynamics of entanglements: A nonlinear model consistent with
the Cox-Merz rule”, Journal of Non-Newtonian Fluid Mechanics, 62, (1996),
279-289.
[14] Larson, R. G., “ Combinatorial Rheology of Branched Polymer Melts”, Macromolecules, 34, (2001), 4556-4571.
[15] McLeish, T. C. B. and Larson, R. G., “ Molecular constitutive equations for aclass of branched polymers: The pom-pom polymer”, Journal of Rheology, 42,(1998), 81-110.
[16] Hua, C. C. and Schieber, J., “ Segment connectivity, chain-length breathing, segmental stretch, and constraint release in reptation models. I. Theory and singlestep strain predictions”, Journal of Chemical Physics, 109, (1998), 10018-10027.
[17] 成沢郁夫，高分子材料強度学, (1982), オーム社
[18] Gao, J. and Weiner, J. H., “ Computer Simulation of Viscoelasticity in Polymer
Melts”, Macromolecules, 25, (1992), 1348-1356.
[19] Gao, J. and Weiner, J. H., “ Contribution of covalent bond force to pressure inpolymer melts”, Journal of Chemical Physics, 91, (1989), 3168-3173.
[20] Mckechnie, J. I., Haward, R. N., Brown, D. and Clarke, J. H. R., “ Eﬀects ofChain Configurational Properties on the Stress-Strain Behavior of Glassy Linear
Polymers”, Macromolecules, 26, (1993), 198-202.

参考文献

54

[21] Clarke, J. H. R. and Brown, D., “ Molecular Dynamics Modelling of Polymer
Materials”, Molecular Simulation, 3, (1989), 27-47.
[22] Yamamoto, T., “ Interphases and Mesophases in Polymer”, Crystallization III,(Allegra, G., ed.), Advances in Polymer Science, 191, (2005), 37-85.
[23] Bennemann, C., Paul, W., Binder, K., Dunweg, B., “ Molecular-dynamics simulations of the thermal glass transition in polymer melts: α-relaxation behavior”, Physical Revew, E57, (1998), 843-851.
[24] Takeuchi, H. and Roe, R. J., “ Molecular dynamics simulation of local chain
motion in bulk amorphous polymers. I. Dynamics above the glass transition ”,
Journal of Chemical Physics, 94, (1991), 7446-7457.
[25] 屋代如月, 内藤正登, 皆川康久, 冨田佳宏, “ 非晶性高分子材料のヒステリシスに関する分子動力学的研究”, 日本機械学会論文集, A72, (2006), 277-284.

[26] Yashiro, K., Naito, M., Minagawa, Y. and Tomita, Y., “On the Hysteresis of
Polyethylene and Polybutadiene: A Molecular Dynamics Study”, Proceeding of
ICCES2005, Advances in Computational and Experimental Engineering and Science, CDROM[27] Lennard-Jones, J. E., Proceedings of the Royal Society, London, ser, A106,(1924), 463.
[28] Morse, P. M., “ Diatomic molecules according to the wave mechanics. II. Vibrasional laevels”, Physical Review, 34, (1929), 57-64.
[29] Daw, M. S. and Baskes, M. I., “ Semiemprical, Quantum Mechanical Calculationof Hydrogen Embrittlement in Metals”, Physical Review Letters, 50-17, (1983),
1285-1288.

参考文献

55

[30] Daw, M. S., Baskes, M. I., “ Embedded-atom method: Derivation and Application to Impurities, Surfase, and Other Defect in Metals”, Physical Review B,
29-12, (1984), 6443-6453.
[31] Rigby, D., Ree, R. J., “Molecular dynamics simulation of polymer liquid andglass.I.Glass transition”, Journal of Chemical Physics, 87, (1987), 7285-7292.
[32] Richard H. Gee, Richard H. Boyd, “Conformational dynamics in bulk polyethylene: A molecular dynamics simulation study”, Journal of Chemical Physics,
101-1, (1994), 788-797.
[33] 上田顯, コンピューターシミュレーション, (1990), 朝倉書店.

学術論文・学術講演
学術論文

 ヒステリシスをもたらすポリブタジエンの分子鎖凝集構造：
分子動力学による検討
古賀裕・内藤正登・皆川康久・屋代如月・冨田佳宏
材料, 投稿中学術講演

 cis-1,4 ポリブタジエンのヒステリシスに関する分子動力学的研究
古賀裕，内藤正登，屋代如月，冨田佳宏
日本機械学会第 19 回計算力学講演会，名古屋大学，(2006.11)

 繰り返し変形を受けるＣＢ充填ゴムの分子動力学シミュレーション
山本雄太，古賀裕，屋代如月，冨田佳宏
日本機械学会関西支部平成 18 年度卒業研究発表講演会，大阪産業大学, (2007.3)

 繰り返し変形下におけるポリブタジエンの可逆圧縮応力：
分子動力学法による検討
古賀裕，内藤正登，屋代如月，冨田佳宏
日本機械学会第 20 回計算力学講演会，同志社大学，(2007.11)

 ヒステリシスをもたらすポリブタジエンの分子鎖構造：分子動力学法による検討
古賀裕・内藤正登・皆川康久・屋代如月・冨田佳宏
日本材料学会関西支部第 2 回若手シンポジウム，有馬簡易保養保険センター,(2007.12)
 分子動力学法によるポリブタジエンとフィラー界面に関する力学挙動評価：界面強度による影響56

学術論文・学術講演

57

山本雄太，古賀裕，屋代如月，冨田佳宏
日本機械学会関西支部

第 83 期定時総会講演会，大阪大学, (2008.3)

 Molecular Dynamics Simulation on cis-1,4 Polybutadiene under Cyclic Deformation: Chain Entanglements and Hysteresis
K. Yashiro, Y. Koga, M. Naito, Y. Minagawa and Y. Tomita
8th World Congress on Computational Mechanics, WCCM8,(2008.6)

謝辞
本研究を遂行するにあたり，御指導を賜りました固体力学研究室の冨田佳宏教授および屋代如月准教授に，心より感謝の意を表します．固体力学研究室に所属した 3 年間では，学術論文，学術講演会といった学外での研究発表の機会を数多く与えて頂きましたが，両先生の時に厳しく，時に優しく，仔細に渡る懇切丁寧な御指導がなければこのような貴重な経験を得ることは到底適いませんでした．ともすれば方途を見失いないがちな私にとって，このような活発な環境の中で研究生活が行えたことは大変幸せなことでした．また，研究報告会などにおいて様々な御助言を賜りました長谷部忠司准教授に厚く御礼申し上げます．高木知弘准教授 (京都工業繊維大学) には，研究生活における種々の相談を受けて頂き，多くの御助言，御協力を賜りました．深く感謝致します．そして，日々の研究が順調に進んだのは，古宇田由夫技官の委細に渡る御配慮があったからこそです．心より御礼を申し上げます．社会人博士として，研究を進める上で貴重な御助言と，有益な議論を頂いた内藤正登氏 (現 SRI 研究開発株式
会社) と，皆川康久氏 (SRI 研究開発株式会社) に厚く御礼申し上げます．また，後輩として活発な議論を交わした山本雄太氏を始め，研究室諸氏に深く感謝致します．特に研究室に配属されて以来共に研究に取り組み，同輩として互いに切磋琢磨し，時には励ましあった同期の朝西誠君，久馬雅彦君，田中克昌君，樋口昌宏君には厚く御礼申し上げます．最後に，六年間の学生生活において，多大な援助と理解を頂いた祖母，
両親に敬意と感謝の意を表します．
平成 20 年 2 月
古賀58

裕