分子動力学法によるポリマー内部
分子鎖構造のモデル化と変形挙動解析

平成 15 年 1 月
神戸大学大学院自然科学研究科
機械工学専攻011T360N
伊藤智啓

Molecular Dynamics Study on
Modeling and Deformation Behaviorof Molecular Chains in Polymer Materials

January 2003
Division of Mechanical Engineering,
Graduate School of Science and Technology,
Kobe University, Kobe, Japan
Tomohiro ITOH

要約
高分子材料は，内部に無秩序に存在する高分子鎖のからみ点数の変化や分子鎖セグメ
ントの回転といった微視的構造変化に依存して複雑な変形挙動を示す．本研究では，
不規則かつ長い高分子鎖の変形挙動について，原子レベルから新たな知見を得ることを目的として，ポリエチレンを対象とした種々の分子動力学シミュレーションを行った．まず，200 万の粒子数からなる非晶性ポリエチレンナノ試験片を用いた引張シミュ
レーションを行い，長く不規則な非晶高分子鎖の複雑な変形挙動を，周期境界等の非物理的な要因を排除して検討した．その結果，引張ごく初期は (1) 結合長の微小変化による変形吸収，(2) 結合角の微小変化による変形吸収，(3)2 面角の微小変化による応力緩和，というメカニズムにより線形弾性応答を示すこと，また，その後は一定応力下で変形が進行する降伏挙動があり，引張後期に応力が急増すること，などが明らかとなった．また，内部分子鎖構造変化の詳細な観察により，降伏⇒硬化過程では分子鎖が複雑にからみあって変形しにくいクラスター状の部分 (からみ点) と，それらを直線状の分子鎖が接続したネットワーク状構造が形成されていること，応力上昇はこのからみ点間の分子鎖が引張方向に配向・延伸したことによるものであること，などを明らかにした．次に，ラメラ構造など，実際のポリマー内部に多く存在する結晶・非晶界面の微視的変形挙動について，結晶相分子鎖方位の違いを考慮して検討した．まず，シミュレーションにより得られた無負荷平衡状態の界面近傍の分子鎖構造を詳細◦に観察し，(1) 結晶相の分子鎖が引張方向に平行なモデル (Model 0 ) では，結晶相の界面近傍に gauche–gauche の構造欠陥による結晶の乱れを多く生じるが，非晶相に大きな変化はなく，結晶・非晶界面は引張軸に対して垂直な平面に保たれること，(2)45◦◦

傾いた界面モデル (Model 45 ) では，界面近傍の非晶分子鎖が結晶化し，引張軸に対◦◦して垂直な平面であった結晶・非晶界面が鋸歯状となること，(3)Model 0 ，45 のい
ずれにおいても，gauche–gauche の構造欠陥が分子鎖方向に伝ぱして結晶の乱れを生じ，それが分子鎖に対して垂直方向に広がりを有すること，などを明らかにした．そ◦の後，得られた界面構造に対して引張シミュレーションを行い，(1)Model 0 では結
晶相はほとんど変形せず，界面近傍の非晶相に変形が集中し，そこに先のナノ試験片◦のシミュレーションと同様の硬化機構を生じて応力が急増すること，(2)Model 45 では van der Waals による弱い凝集力で結合している結晶相分子鎖の横方向に力が作用するため，非晶相だけでなく結晶相もひずみに応じた変形を生じ，引張初期から応力が線形的に増加する挙動を示したこと，などを明らかにした．

Summary
Polymeric materials show complex deformation behavior depending on the internal
microscopic phenomena such as the cross-link or the chain network structure. In thepresent study, several molecular dynamics simulations are conducted on the amorphousand crystal/amorphous polyethylene under uniaxial tension, for fundamental understanding of the mechanical behavior of the irregular, long and massive molecular chainsin polymeric materials. A tensile simulation is implemented with a nanoscopic specimen of amorphous polyethylene involving about 2 million methylene groups. After
showing linear elastic relation at the initial stage, the material “yields” by elongating
without stress increase until strain hardening appears. Careful investigation of changes
in the dihedral angle and morphology of all molecular chains reveals that the gauche⇒ trans transition takes place during yielding, and generating new network-like struc-ture composed of entangled molecular clusters and oriented chains bridging them. The
strain hardening is due to the directional orientation and stretch of molecular chains
between entanglements in the nucleated structure. Two other simulations are conducted on the deformation behavior of the crystal/amorphous interface normal to thetensile direction. One has the crystal of which chain direction is parallel to the tensile◦

axis (Model 0 ). The other has the chain direction rotated 45 degrees against the◦tensile direction (Model 45 ). Detail observation of the chain structure at the inter◦face under unloaded equilibrium condition reveals that the interface in Model 45 is◦serrated by the re-crystallization of amorphous chains while Model 0 has no remarkable change in the morphology of the interface . It is also observed in both Model 0◦◦

and Model 45 that the gauche–gauche defects propagate in the chain direction andmake disorders in the crystal phase. Then the tensile strain is applied on the models,◦leading the following results: (a) Model 0 shows the yielding and (b) strain hardening
as similar as the nanoscopic specimen above mentioned. The chain orientation andelongation take place in the amorphous phase at the crystal/amorphous interface. (b)◦Model 45 does not show the yielding behavior (c) but linear stress increase from the
beginning of the straining. It is because the molecular chains in the crystal phase isalso subjected to the force normal to the chain direction, so that the crystal deforms
in the direction and shows relatively large stress increase.

目次
第1章 緒 論1

第 2 章 解析手法の基礎

32.1

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

3

2.2

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

4

2.3

高速化手法

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

12

2.4

応力の評価

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

14

第 3 章 非晶性ポリマー内部における分子鎖配向およびからみ点形成機構

15

3.1

シミュレーション方法 . . . . . . . . . . . . . . . . . . . . . . . . . . .

15

3.2

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

17

3.2.1

初期降伏 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

3.2.2

降伏⇒硬化挙動 . . . . . . . . . . . . . . . . . . . . . . . . . . .

21

3.2.3

分子鎖配向とからみ点形成 . . . . . . . . . . . . . . . . . . . . .

23

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

27

3.3

第 4 章 結晶・非晶界面における微視的変形機構

28

4.1

シミュレーション方法 . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

4.2

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

30

4.2.1

初期平衡構造 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

30

4.2.2

初期緩和における構造欠陥発生・伝ぱ. . . . . . . . . . . . . .

32

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

42

4.3

第5章 結 論

44

参考文献

47

i

目次

ii

関連発表論文・講演論文

51

謝辞

89

第1章
緒論
高分子 (ポリマー) 材は，その加工性，軽量性，耐久性から金属に代わる新たな材料として日常製品から機械構造部材まで広範に使用されている．しかしながら，その力学的特性は，材料内部に無秩序に存在する高分子鎖のからみ点数の変化や分子鎖セグ
メントの回転といった，複雑な微視的内部挙動に大きく依存し，未だ不明な部分が多い [1] ．ポリマー材の力学的特性を生かした製品の設計および製造のために，材料内部の微視的変形挙動を反映した精密なシミュレーションモデルの構築が急がれている．
ポリマー内部のからみ点の変化を考慮した分子鎖網目論など，微視的内部状態の変形を考慮した構成式モデルが提案され [2]–[6] ，実験結果を良好に再現することが示されている．しかしながら，現象論的に導入されたこれら構成式において，無秩序な分子鎖集合体における個々の分子鎖の変形挙動，および，多数の分子鎖の集団的挙動がどのように反映されているかは未だ不明である．
一方，近年の計算機能力の飛躍的向上を背景に，分子動力学をはじめとする原子シ
ミュレーションによって，材料内部で生じる nm ∼ µm のスケールの複雑な変形挙動を直接追及することが可能になりつつある [7], [8] ．ポリマー材料への分子動力学法の適用も盛んに試みられており，ポリマーの結晶化に関する分子シミュレーションとして，結晶化の一次核形成過程，結晶成長過程および流動場での結晶化等の多彩な研究 [9]–[16]が行われている．他にも，Madkour は統計力学と分子動力学を合わせた新たな手法を提案し，溶媒中のポリマーの混和性を評価することを試みている [17] ．また，青柳ら
は原子結合を考慮した詳細な分子モデルから，ビーズスプリングモデルのような粗視

1

第1章

緒論

2

化した分子モデルまでを対象にした，汎用の分子動力学シミュレーションプログラム(COGNAC)[18] を提案し，定温条件，定圧条件，伸張やせん断等の変形，電場印加等の様々な条件の下でシミュレーションを行っている．これらのポリマーシミュレーショ
ンで通常用いられる united atom model(C 原子に結合する H 原子を陽には扱わず，C原子に繰り込んでしまう力場) は，非晶性高分子材料では問題を生じることが指摘されていたが [19]–[22] ，桑島らにより，密度，体膨張率，圧縮率，拡散係数等の実験値を良好に再現するポリエチレンやポリプロピレン等の非晶性高分子のポテンシャルパラメー
ターが提案された [23] ．渋谷らは，桑島らのポテンシャルを用いてポリエチレンの一軸引張シミュレーションを周期セルを用いて行うことにより，分子鎖集合体の弾性・粘性特性を評価し，これをばねとダッシュポットで表される Voigt モデルに反映させる
ミクロ–メゾスケール連結モデルを提案した [24] ．しかしながら，数百程度の粒子からなるセルに周期境界を適用していたため，長く複雑な非晶性高分子鎖の力学挙動を反映しているとは言えず，特に大ひずみ領域ではその影響が顕著に表れるという問題があった．また，実際のポリマー材は 100 % 結晶性・非晶性で存在することはまれであり，通常結晶相と非晶相が入り混じった 2 相構造となっている．このため，ポリマー材を結晶相と非晶相からなる複合材料と考えたモデル化もなされているが，結晶・非晶界面における微視的変形挙動はいまだ検討されていない．
本研究では，非晶相内部で生じる分子鎖の配向やからみ点形成のメカニズム，および，結晶・非晶界面領域の微視的挙動等について新たな知見を得ることを目的として，種々の分子動力学シミュレーションを行った．以下に各章の概要を示す．
第 2 章では解析手法の基礎として，分子動力学法を簡単に説明し，分子動力学計算で最も重要となるポテンシャルエネルギーについて述べた．また，大規模計算を行うための高速化手法および応力の評価方法を示した．第 3 章では，長く不規則な非晶性高分子鎖の変形挙動を，周期境界等による非物理的な要因を排して検討するため，約200 万粒子からなるポリエチレンナノ試験片を用いて引張シミュレーションを行い，引張時の分子鎖配向やからみ点形成等のメカニズムについて検討した．第 4 章では，結晶相と非晶相の界面モデルを作成し，界面領域特有の構造変化について検討した．また，結晶相の chain 方向の影響についても調べた．

第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 の方法，予測子–修正子法等がよく用いられる [29] ．
本研究では，以下に示す 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
3

(2.5)

第2章

解析手法の基礎

4であり，式 (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 ポテンシャル [25] や Morse ポテンシャル [26] などの経験的 2 体ポテンシャルは取り扱いが簡単であるため，従来からよく用いられてきたが，金属材料表面や異相界面など粒子配置が急激に変化する現象を扱う近年の研究では，多体力を考慮したより精度の高いポテンシャル関数 [27], [28] が用いられている．
本研究で扱うポリマー材料では，強い共有結合と弱い非共有結合からなる．さらに，共有結合には分子鎖内部の結合角度や回転等の内部自由度があるため，ポテンシャル

第2章

解析手法の基礎

5

エネルギー Etot は以下のような形で表される．

Etot = ΦBS (r) + ΦBE (θ) + ΦT O (ϕ) + ΦV W (r̄)(2.12)

右辺各項はそれぞれ分子内の炭素間の結合長 r (図 2.1)，結合角 θ (図 2.2)，2 面角 ϕ (図2.3)，異分子鎖間および同一分子鎖内の十分離れた非結合原子間距離 r̄ (図 2.4) に対す
るポテンシャルを表す．
本研究で解析対象とするポリエチレンに対しては，メチレン炭素 (CH2 ) を単一粒子として扱う united atom model を用いてパラメーターがフィッティングされた以下の関数形が提案されている [23] ．

bond stretch ポテンシャル（2 体間ポテンシャル）ΦBS (r) =∑{kr (r − r0 )2}(2.13)

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

torsion ポテンシャル（4 体間ポテンシャル）
ΦT O (ϕ) =∑(V1 cos ϕ + V2 cos 2ϕ + V3 cos 3ϕ + V6 cos 6ϕ)(2.15)

van der Waals ポテンシャル（2 体間ポテンシャル）
ΦV W (r̄) =∑{A (r̄)−12 − C (r̄)−6}(2.16)

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

第2章

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

Table 2.1 Potential parameter for bond stretch.

CH2 – CH2r0kr[nm] [kJ/ (mol · nm2 )]0.15331.373 × 105

ΦBS [eV]

100

50

r
0
0

2
r [A]

4

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

6

第2章

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

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

Table 2.2 Potential parameter for bending.
θ0[deg.]C – CH2 – C 113.3kθ
[kJ/( mol・rad2 )]374.7

ΦBE [eV]

100

θ
50

0
0

180
θ [deg.]360

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

7

第2章

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

ΦT O (ϕ) =∑(V1 cos ϕ + V2 cos 2ϕ + V3 cos 3ϕ + V6 cos 6ϕ)

Table 2.3 Potential parameter for torsion.

C – CH2 – CH2 – C

V1
V2
V3
V6
[kJ/mol] [kJ/mol] [kJ/mol] [kJ/mol]
3.935
2.177
7.786
0.0

ΦTO [eV]

0.1

0

–0.1
0

φ

gauchegauchetrans180
φ [deg.]360

Fig.2.3 Relationship between torsion potential ΦT O and torsion angle ϕ.

8

第2章

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

ΦV W (r̄) =∑{A (r̄)−12 − C (r̄)−6}Table 2.4 Potential parameter for van der Waals.

CH2 – CH2

A
[kJ/mol・nm12 ]2.972 × 1019

C
[kJ/mol・nm6 ]6.907 × 109

ΦVW [eV]

0.01

r
r

0

4

6
r [A]

8

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

9

第2章

解析手法の基礎102 体間相互作用
r = |r ij |Fi = −

Φ′ (r)
r ij
r(2.17)

F j = −F i

(2.18)

3 体間相互作用r ik · r jk
|r ik | |r jk |{}Φ′ (θ)
1
r ik · r jk
Fi =r jk −r ik
sin θ |r ik | |r jk ||r ik |2cos θ ={

r ik · r jk
Φ′ (θ)
1
r ik −r jk
Fj =sin θ |r ik | |r jk ||r jk |2

(2.19)}

F k = −F i − F j(2.20)(2.21)

4 体間相互作用r ik · r kl
r kl
|r kl |2
r jl · r kl
r kl
Aj = r jl −|r kl |2
Ai · Ajcos ϕ =
|Ai | |Aj |{}Φ′ (ϕ)
1
Ai · Aj
Fi =
Aj −
Ai
sin ϕ |Ai | |Aj ||Ai |2
Ai = r ik −{1
Φ′ (ϕ)
Ai · Aj
Fj =
Ai −
Aj
sin ϕ |Ai | |Aj ||Aj |2r jl · r kl
r ik · r kl
Fi −F k = −F i −Fj|r kl |2
|r kl |2
r jl · r kl
r ik · r kl
Fj
F l = −F j +2 Fi +|r kl |
|r kl |2

(2.22)}(2.23)(2.24)(2.25)

第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 bending angle θ.

j
rjl

φ

k

Ai rik

Aj

rkl
l

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

11

第2章

2.3

解析手法の基礎

12

高速化手法

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

rfc

rc

Fig.2.8 Schematic of bookkeeping method.

第2章

解析手法の基礎

13

別の高速化手法としてブロック分割法がある．図 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.

第2章

2.4

解析手法の基礎14

応力の評価

材料の変形・破壊において，その力学状態を記述するパラメーターとして応力 σ が用いられる．負荷を受けた材料の内部には内力が発生し，外力と釣り合う．材料内部の任意の面に働く内力をその面積で除したものが応力 σij であり，2 階の対称テンソルとなっている．
構造や負荷が単純ではない場合や，材料が均一でない場合には，応力は均一にならず材料内部に分布が現れる．連続体力学では，これらを局所領域に対する量として扱うことで変形や破壊を論じてきた．これまでに，原子モデルを用いる場合についても様々なものが提案されている [29], [30] ．
本研究では，連続体力学の応力の定義を拡張して以下のような応力を定義する．まず，ある粒子位置で座標軸と垂直な断面を考える．平衡状態にある系では，その面の片側から粒子が受ける力は内力 (図 2.10) であり，それを断面積で除したものをその粒子位置での局所応力と定義する．以下のシミュレーションではこの定義に基づき応力を評価した．

internal forceyx
Fig.2.10 Internalforce acting on plane in one molecule coordinate.

第3章
非晶性ポリマー内部における分子鎖配向およびからみ点形成機構
本章では，長く不規則な非晶性高分子鎖の複雑な力学挙動について，周期境界等の非物理的な要因を排除して検討するために，分子鎖数 3,542，約 200 万の粒子からなる
ポリエチレンナノ試験片を用いて引張分子動力学シミュレーションを行い，内部に生じる個々の分子鎖挙動を詳細に観察する．

3.1

シミュレーション方法

非晶性高分子における不規則かつ長い分子鎖の挙動をシミュレートする場合，周期境界の適用は非物理的な拘束を生じる可能性がある．そこで本研究では，図 3.1 に模式的に示すような形状のポリエチレンナノスケール試験片を解析対象とすることとし，以下の手順で試験片内部の初期分子鎖構造を作成した．

(1) 試験片となる空間上に，分子鎖を成長させる核となる粒子をランダムに配置する．
◦(2) 結合長 r および結合角 θ は最安定値 0.1533 nm と 113.3 ，2 面角は最安定な trans◦または準安定な gauche±5 の値を保つようにしながらランダムに成長させる (図3.2)．
(3) 成長させた分子鎖の端点が試験片領域外に達した場合，または他に成長させた分子鎖の存在により成長できなくなった場合，そこを終点とし別の分子鎖の作成を

15

第3章

非晶性ポリマー内部における分子鎖配向およびからみ点形成機構

16行う．
上記手順を試験片の密度が低密度ポリエチレンの値 (0.87 g/cm3 ) になるまで繰り返した．図 3.1 中に示したように，試験片平行部の長さは 50 nm，平行部断面は 26nm×26nmの正方形であり，上下に 5 nm の厚さのつかみ部を有する．作成した初期構造における分子鎖数は 3,542，1 分子鎖あたりの粒子数は 500 ∼ 1,500，総粒子数は 1,984,434 である．試験片上端から下端まで連結した分子鎖はなく，上下端つかみ部にある分子鎖の他端は試験片内部または表面にある．
まず，上下端のつかみ部にある粒子の z 方向の運動を拘束し，3000 fs の初期緩和計算を行った．数値積分には Verlet 法を用い，積分の時間ステップは 0.1 fs とした．温度は 300 K とし速度スケーリング法により制御した．
その後，得られた初期平衡状態をひずみの基準とし，z 軸方向に ε̇zz = 5.0 × 1011 s−1(5.0 × 10−5 /step) でひずみを増加させる引張シミュレーションを行った．ひずみは z

50nm

5nm

方向の粒子間隔を均等に広げることにより制御した．

5nm

z
y
x

nm
36

5nm 26nm 5nm

Fig.3.1 Dimensions of nanoscopic specimen of amorphous polyethylene

第3章

非晶性ポリマー内部における分子鎖配向およびからみ点形成機構

17φr = 0.1533 nm
r 1 = |r cos θ |r 2 = |r sin θ |θ = 113.3
φ = gauche, trans

θ

X
r2

z

X1
r1

r

y
x

X0

Fig.3.2 Modeling of random coil molecular chain

3.2

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

3.2.1

初期降伏

初期平衡状態における各粒子間の結合長 r，結合角 θ，2 面角 ϕ および各粒子位置での局所密度 ρ (van der Waals エネルギーのカットオフ距離 0.8 nm を半径とする球内粒◦子密度) の分布を図 3.3 に示す．結合長 r，結合角 θ はそれぞれ 0.153 nm，113 に鋭いピークを有し，ほとんどのノード (結合角，結合長) が安定な値をとっている．2 面角 ϕ◦◦は最も安定である trans 点 (180 ) と準安定な gauche 点 (67.5 ) にピークがあるが，その周囲の角度にも多く存在しており，複雑に入り組んだ構造においてやや高いエネル
ギー状態にあるノードが存在していることが示唆される．密度は設定した 0.87 g/cm3にピークがあるが，表面の存在により低密度側の粒子も多く認められる．

第3章

非晶性ポリマー内部における分子鎖配向およびからみ点形成機構

Number of node5
[×10 ]
3.0

6.0

4.0

2.0

0
0.1

Number of node4
[×10 ]
8.0

0.15

trans

6.0

4.0

gauche
2.0

0

60

120

180

Torsion angle, φ , deg.

2.0

1.0

0

0.2

Bond stretch, r, nm

60

120

180

Bending angle, θ , deg.
5
[×10 ]
1.0

Number of node

Number of node4
[×10 ]
8.0

18

0.8
0.6
0.4
0.2

0

0.5

Density,

1.0

1.5

ρ, g/cm3

Fig.3.3 Distribution in bond stretch, bending angle, torsion angle and atomic
density.

第3章

非晶性ポリマー内部における分子鎖配向およびからみ点形成機構

19

引張シミュレーションにより得られた，εzz = 0.0 ∼ 0.08 の引張ごく初期における応力–ひずみ関係を図 3.4 に示す．応力は，上下端つかみ部粒子に作用している力の z 方向成分から試験片に作用する軸力を算出し，これを平行部断面積で除すことにより評価している．ひずみ 0.012 までほぼ線形関係を示した後，εzz = 0.012 (図中矢印 ⃝A)において最初の応力ピークを示し応力はわずかに低下する．その後，再び応力は上昇し，εzz = 0.027 (図中矢印 ⃝B ) においてピーク応力 0.92GPa を示して応力は大きく減少した．以降は，ひずみの増加にも関わらず，応力は振動しながらほぼ一定の値を示し変形が進行している．

εzz = 0.0 ∼ 0.08 における系のポテンシャルエネルギー，および，その成分 (bondstretch, bending, torsion, van der Waals) の変化を図 3.5 に示す．なお，ポテンシャル
エネルギーの各成分の値はオーダーが異なるため，初期平衡状態における値 (E0 ) か
らの変化量を示している．また，図 3.4 で示した二つの応力ピークに対応するひずみを細線 ⃝A および ⃝
B で示している．⃝
A, ⃝B いずれも，ほぼ bending ポテンシャルが極大，

Tensile stress, σzz, GPa

1

B

0.5

A
0

-0.5

0

0.02

0.04

0.06

0.08

Applied strain, εzz
Fig.3.4 Change in tensile stress at the initial stage of the tensile simulation(εzz ≤ 0.08).

非晶性ポリマー内部における分子鎖配向およびからみ点形成機構

Change in potential energy, ∆E = E - E0, eV

第3章[×10 -3]
0.6

∆Etot

20

∆Eto

0.4
∆Ebe

0.2
∆Ebs
∆Evw

0

Peak A Peak B
-0.2

0

0.02

0.04

0.06

0.08

Applied strain, εzz
Fig.3.5 Change in potential energy at the initial stage of the tensile simulation (εzz ≤ 0.08).

かつ bond stretch ポテンシャルが極小点を経て増加しはじめる点となっている．また，系の応力変化としては現れていないが，⃝A, ⃝B の前には上昇した bond stretch ポテン
シャルが減少し，bending ポテンシャルが増加するエネルギーシフトが生じている．以
上を統合して考えると，応力ピーク前後においては次のような内部構造変化を生じていることが示唆される; (1) 変形初期は結合長の微小変化により変形を吸収, (2) 結合長による変形吸収が限界に達すると，結合角が微小変化し bond stretch ⇒ bending のエネルギーシフトを生じる，(3) 結合角による変形吸収限界に達したとき応力低下を生じる．ここで，torsion ポテンシャルが ⃝A において急減，⃝
B において極小点となっていることから，応力低下には 2 面角変化が重要な役割を担っていることがわかる．すべての 2 面角変化に関する詳細な観察 [31] より，(1) ⃝A における変化は，初期内部構造においてやや高いエネルギー状態にあった 2 面角が，引張を駆動力として緩和されたことによるもの，(2) ⃝B における変化は，⃝
A の緩和により最も低いエネルギー状態となった 2 面角が，結合角の変形を吸収するため変形したことによるもの，(3) ⃝
A で生じ

第3章

非晶性ポリマー内部における分子鎖配向およびからみ点形成機構

21◦

た 2 面角変化は ±30 以下であり，安定角間 (gauche ⇔ trans) を遷移するような回転は◦生じていないのに対し，⃝
B では ±60 以上の回転が生じていたこと，などが明らかと
なった．

3.2.2

降伏⇒硬化挙動

ひずみ 3.0 まで行った，引張シミュレーション全体を通しての引張応力−ひずみ関係を図 3.6 に示す．前述の初期降伏以降，引張ひずみ約 1.5 までの領域では応力は振動しながら約 0.75GPa の値を保ち，大きな応力増加を伴なわずに変形が進行した．一方，

εzz > 1.5 の変形後期には著しく応力が上昇した．実際のポリマー材料の引張応力–ひずみ関係と直接対応するものではないが，図 3.6 の応力–ひずみ関係とマクロのそれとの類似性から，εzz ≤ 1.5 を降伏領域，εzz > 1.5 を硬化領域と称することにする．図 3.6の応力変化に対応するポテンシャルエネルギー変化を図 3.7 に示す．εzz ≤ 1.5 の領域
において，bond stretch, bending が単調増加しているのに対し，torsion ポテンシャル

Tensile stress,

σzz , GPa

2.0

1.5

1.0

0.5

0

1.0

Applied strain,

2.0

3.0

ε zz

Fig.3.6 Change in tensile stress in the tensile simulation up to the strain of
εzz = 3.0.

第3章

非晶性ポリマー内部における分子鎖配向およびからみ点形成機構

22

Change in potential energy, ∆ E = E - E0, eV[×10 -3]
4.0

∆Etot
3.0

∆Ebe

∆Ebs

2.0

1.0

∆Evw

0.0

∆Eto
-1.0

0

1

Applied strain,

2

3

ε zz

Fig.3.7 Change in potential energy in the tensile simulation up to the strainof εzz = 3.0.

のみが εzz = 0.5 近傍で増加から減少へと変化している．このことから，降伏領域の変形挙動にも 2 面角変化が重要な役割を果たしていることがわかる．また，硬化が開始
する εzz = 1.5 近傍では torsion ポテンシャルは初期平衡状態より低くなっており，降伏領域において gauche ⇔ trans 遷移の内部構造変化を生じていることが示唆される．
降伏⇒硬化領域において生じた内部構造変化を明らかにするため，各ひずみにおけるすべての 2 面角 ϕ の変化を調べた．図 3.8 は，2 面角が ∆εzz = 0.05 前の状態から ±60◦以上変化して gauche ⇒ trans または trans ⇒ gauche の遷移を生じたノードの数を，横軸にひずみをとって示したものである．εzz = 0.027 の初期降伏直後から εzz = 0.5 の領域において，系のエネルギーを減じる gauche ⇒ trans 遷移は単調に減少しているのに対して，系のエネルギーを増加させる trans ⇒ gauche 遷移は 4000∼4500 程度のほぼ一定の値をとっている．これは，gauche ⇒ trans 遷移による trans 点の増加と，trans ⇒gauche の発生がバランスしたためと考えられる．一方，ひずみ 0.5 以降はいずれの遷移も単調に減少し，かつ gauche ⇒ trans 遷移数は trans ⇒ gauche 遷移数の倍程度生じ

第3章

非晶性ポリマー内部における分子鎖配向およびからみ点形成機構

23

Number of dihedral nodes, Ntransition[×10 4]1.0
gauche => transtrans => gauche

0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

0

1

2

Applied strain,

3

ε zz

Fig.3.8 Change in the number of dihedral nodes showing gauche ⇔ trans transition.

ていることから，内部では gauche 点が減少し trans 点が増加する．すなわち，分子鎖が直線状になる構造変化が降伏⇒硬化領域において生じている．

3.2.3

分子鎖配向とからみ点形成

引張変形時の分子鎖構造変化の例を図 3.9 に示す．図では試験片平行部中央から10 nm × 10 nm × 3 nm の領域を抽出して示している．黒く着色した粒子は，εzz = 0の初期平衡状態において分子鎖折れ曲がり点にあると判断した粒子であるが，これについては後述する．εzz = 0.5 において，全体的にはほぼ均一に変形しているが，アモ
ルファス金属の場合 [32] と異なり部分的に密度の低い領域が現れている．εzz = 1.0 の
図では密度の不均一さが顕著となり，クラスター状に密度の高い部分が存在している．
また，密度の低い部分には直線状に伸びた分子鎖が多く認められる．硬化が開始する

εzz = 1.5 では，円 A および B で囲ったような，分子鎖が複雑にからまりあったクラス
ター状の部分を，長く直線状に伸びた分子鎖が連結している構造になっている．また，

εzz = 1.0 の図と比べてクラスター間の分子鎖は引張方向に配向している．εzz = 1.5 以

第3章

非晶性ポリマー内部における分子鎖配向およびからみ点形成機構

ε zz = 0.0

ε zz = 0.5

ε zz = 1.0
B

ε zz = 1.5
AB

ε zz = 2.0

AB

ε zz = 3.0

A
Fig.3.9 Snapshots of molecular chains.

24

第3章

非晶性ポリマー内部における分子鎖配向およびからみ点形成機構

25

降は，クラスターの解消も認められるが，基本的には εzz = 1.5 で形成されたクラス
ターネットワーク構造のまま，クラスター間の分子鎖が引張方向に大きく延伸して変形している (εzz = 2.0 および εzz = 3.0 の図)．試験片上端から下端まで連結した分子鎖のない本シミュレーションモデルにおいて生じた応力の急増は，このような「からみ点」による分子鎖ネットワーク形成およびからみ点間分子鎖の配向・延伸によるものと結論づけられる．
降伏⇒硬化過程における局所密度変化を定量的に示したものが図 3.10 である．引張により密度分布のピークは低密度側にシフトするが，図 3.9 で観察したような密度があまり変化しないノードの存在により高密度側に分布の裾野が残る (εzz = 0.5, 1.0)．

εzz = 1.5 において ρ = 0.53 g/cm3 近傍の高密度側にもピークが現れており，硬化を開始するひずみ近傍において新しい内部構造が形成されたことが定量的に示される．
図 3.9 で観察した「からみ点」について，初期平衡状態における分子鎖構造との対応を明らかにするため，分子鎖の折れ曲がり点に着目した検討を行った．図 3.11 に模式的に示すように，あるメチレン粒子 A から前後に 10 粒子離れた粒子 B, C への位置ベ

Rate of node number, N/Ntotal

0.15

ε zz =1.5
0.10

ε zz =1.0
ε zz =0.5
ε zz =0.0

0.05

0

0.2

0.4

0.6

0.8

1.0

Atomic density, ρ , g/cm 3
Fig.3.10 Change in distribution of atomic density.

第3章

非晶性ポリマー内部における分子鎖配向およびからみ点形成機構

26

A

10 monomers

θ'
10 monomers

Cmolecular
chain
B

Fig.3.11 Schematic of evaluation of flexion node.

◦

クトルを決定し，両ベクトルの内角 θ が 90 以下のものを折れ曲がり点として評価す
る．先に触れたように，図 3.9 中の黒い粒子は εzz = 0 の初期平衡状態において折れ曲がり点と判断された粒子であり，各ひずみ下で折れ曲がり点を判断して着色したものではない．図 3.9 の εzz = 3.0 の図から，からみ点となる部分には初期配置において折れ曲がり点であったノードが密集していることがわかる．しかしながら，引張方向に長く配向した部分にも黒い粒子が多数認められ，変形の進行とともに解消される折れ曲がり点も多いことを示している．すなわち，初期平衡状態において折れ曲がり点が密集した領域と，降伏⇒硬化過程において現れるからみ点は必ずしも対応しない．モデルのスケールアップの観点から興味深い，からみ点の位置やからみ点間距離などは，初期平衡状態における分子鎖構造で単純に特徴付けられるものではなく，局所の複雑な変形履歴に依存して変化するものと推測される．

第3章

3.3

非晶性ポリマー内部における分子鎖配向およびからみ点形成機構

27

結言

不規則かつ長い非晶性高分子鎖の変形挙動について，周期境界による非物理的な要因を排除して検討することを目的として，内部に分子鎖の終点や複雑に入り組んだ構造を有するポリエチレンナノ試験片の引張分子動力学シミュレーションを行った．得
られた結果を要約して以下に示す．

(1) 本シミュレーションにおける引張条件 (ε̇zz = 5.0 × 1011 s−1 ) では，引張ごく初期に線形弾性応答を示した後，ほぼ一定応力で変形が進行する降伏挙動が見られた．また，εzz = 1.5 以降に再び応力が上昇するひずみ硬化挙動が認められた．

(2) 系のエネルギー変化の詳細な観察から，線形弾性応答から初期降伏にかけては，⃝
1 結合長の微小変化による変形吸収，⃝2 結合角の微小変化による変形吸収，⃝
3 2面角の微小変化による応力緩和，というメカニズムがあることが示された．

(3) 2 面角変化の詳細な観察から，一定応力下で変形が進行した降伏過程においてはgauche ⇒ trans 遷移が優位に生じており，分子鎖が直線状になる内部構造変化が生じていることが示された．

(4) 内部分子鎖構造変化の直接観察ならびに各粒子近傍の局所密度分布の変化から，
降伏⇒硬化過程では分子鎖が複雑にからみあって変形しにくいクラスター状の部
分 (からみ点) と，それらを直線状の分子鎖が接続したネットワーク状構造が形成されていることがわかった．

(5) 試験片上端から下端まで連続した分子鎖が存在しない本シミュレーションモデル
において生じた εzz = 1.5 以降の応力上昇は，(4) の分子鎖ネットワーク構造の形
成およびからみ点間分子鎖の配向・延伸によるものと結論づけられた．

第4章
結晶・非晶界面における微視的変形機構
前章において，完全非晶性のポリエチレン内部の分子鎖変形挙動について検討した．
一方，実際のポリマー材料では，非晶性と称されるものでも内部には一部結晶化したような領域が存在する．また，結晶性ポリマーについても，微視的にみると内部に非晶相を含む 2 相複合構造であることが知られており，これら結晶・非晶界面を考慮したシミュレーションが必要となっている．本章では，ポリマー内部に存在する結晶相と非晶相の界面領域における微視的変形挙動を原子レベルから明らかにすることを目的として，２種類の界面モデルを用いた引張シミュレーションを行った．

4.1

シミュレーション方法

図 4.1 に模式的に示すような 2 つのシミュレーションセルを用いて解析を行った．1 つ◦は引張方向に配向した結晶相を非晶相と接合したモデル (Model 0 ) であり，もう 1 つ◦◦は結晶相の分子鎖方向が引張軸に対して 45 の角度をなすモデル (Model 45 ) である．
いずれのモデルも，x，y 方向には周期境界条件，z 方向はセルの上下 0.5 nm の厚さの領域を変位制御部とし，粒子の z 方向運動を拘束する境界条件を適用している．結晶相は◦◦結合長 r，結合角 θ，2 面角 ϕ がそれぞれ最安定な 0.1533 nm，113.3 ，180 (trans 点) となる平面ジグザグ構造の分子鎖が，格子間隔 a = 0.4 nm で等間隔に林立した構造としている．実際のポリエチレンの結晶構造は斜方晶 (orthorhombic) と単斜晶 (monoclinic)

28

第4章

結晶・非晶界面における微視的変形機構

29

z
x

10nm
z

amorphous

10nm

amorphous

crystalycrystal(a) Model 0

m

x

10nm

10n

◦

z
x

20nm
z

amorphous

20nm

amorphous

crystal

x

20nm

y

crystal(b) Model 45◦

Fig.4.1 Schematic of two simulation models with crystal/amorphous interface: (a) molecular chains in crystalline phase are aligned in the◦◦tensile direction (Model 0 ), (b) the chain direction is rotated 45◦
against the tensile axis (Model 45 ).

m
20n

第4章

結晶・非晶界面における微視的変形機構

30をとることが X 線回折により明らかにされている [33], [34] が，ここではモデル結晶として単純化を行っている．非晶相における分子鎖構造は，前章と同様の手順に加えて，結晶相分子鎖終端からの成長も乱数を用いて確率的に行っている．ただし，結晶相 → 非
晶相 → 結晶相と連続するような分子鎖は考慮していない．いずれのモデルも，セル全体の平均密度が高密度ポリエチレンの値 (0.94 g/cm3 ) となるまで非晶領域の分子鎖成長を行っている．
◦

Model 0 のセル寸法は 10 nm × 10 nm × 10 nm，作成した初期構造における分子鎖数は 958，1 分子鎖あたりの粒子数は 70 ∼ 140，総粒子数は 81,134 である．Model 45◦のセル寸法は 20 nm × 20 nm × 20 nm，初期構造における分子鎖数は 2,538，1 分子鎖あたりの粒子数は 200 ∼ 400，総粒子数は 654,295 である．作成した初期分子鎖構造に対し，垂直応力が零に近づくようにセル寸法を微調整しながら 10000 fs の初期緩和計算を行った．数値積分には Verlet 法を用い，積分の時間ステップは 0.1 fs とした．温度は 300 K とし，速度スケーリング法により制御した．
その後，得られた初期平衡状態におけるセルをひずみの基準とし，z 軸方向にひずみ増分 ∆εzz = 5.0 × 10−4 を与え，10 fs(=100 MD step) の緩和計算を行うステップを繰り返してひずみを増加させた．ひずみは z 方向の粒子間隔を均等に広げることにより制御した．なお，ひずみ速度に換算すると ε̇zz = 5.0 × 1010 s−1 となり，前章より 1
オーダー小さい条件となっている．

4.2

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

4.2.1

初期平衡構造初期緩和シミュレーションにより得られた，無負荷平衡状態における各粒子間の結合長 r，結合角 θ，2 面角 ϕ および各粒子位置での局所密度 ρ(前章と同じ定義とする)の分布を図 4.2 に示す．どちらのモデルにおいても，結合長 r，結合角 θ はそれぞれ◦

0.153 nm，113 に鋭いピークを有し，ほとんどのノード (結合角，結合長) が安定な値をとっている．2 面角 ϕ については，結晶相の直鎖状分子鎖構造により trans 点に多く

第4章

結晶・非晶界面における微視的変形機構

3[×10 ]
2.0

4[×10 ]

Number of node

Number of node

1.0
1.5

1.0

0.5

0.8
0.6
0.4
0.2

0
0.1

0.15

0

0.2

60

Bond stretch, r, nm

5.0

trans

4.0
3.0
2.0
1.0
0

4.0
3.0
2.0
1.0

gauche
60

120

180

0

0.5

Torsion angle, φ , deg.

Density,(a) Model 0

1.0

0.5

0
0.1

0.15

6.0

4.0

2.0

0

0.2

4
[×10 ]
3.0

60

120

180

Bending angle, θ , deg.
4
[×10 ]
1.5

trans

2.5

Number of node

Number of node

1.5◦

Bond stretch, r, nm

2.0
1.5
1.0

0

1.0

ρ, g/cm3

4[×10 ]
8.0

Number of node

Number of node4
[×10 ]
1.5

0.5

180

3
[×10 ]
5.0

Number of node

Number of node3
[×10 ]
6.0

120

Bending angle, θ , deg.

gauche

60

120

180

1.0

0.5

0

0.5

Torsion angle, φ , deg.

(b) Model 45

Density,

1.0

1.5

ρ, g/cm3◦

Fig.4.2 Distribution in bond stretch, bending angle, torsion angle and atomic
density at initial equilibrium state.

31

第4章

結晶・非晶界面における微視的変形機構

32のノードが集中しているが，非晶相内における折れ曲がり構造により gauche 点にもゆ◦るやかなピークが認められる．なお，Model 45 の方が分子鎖数が多いため，gauche
点の数が多く，かつ，trans および gauche をピークとする分布は広がったものとなっ◦◦ている．同様の理由から，局所密度分布も Model 0 に比べ Model 45 の方がなだらかな分布となっている．いずれも結晶相による高密度側のピークと，非晶相による低密度側のピークの 2 つに分岐しており，2 相構造を保っている．なお，上下端表面 (変位制御部) も含めて評価しているため，表面粒子は低い局所密度を示し，各ピークの低◦密度側にも小さなピークを生じている．ただし，Model 45 の結晶相表面部のピークは非晶相と結晶相の広がった分布に隠れて見えない．

4.2.2

初期緩和における構造欠陥発生・伝ぱ◦Model 0 の初期構造緩和時に観察された結晶・非晶界面近傍の構造変化を図 4.3 に示す．図では 2 面角が trans 点近傍にある粒子を白，それ以外のものを黒で示している．結晶から非晶へと成長する際に 2 面角が変化しているので，t = 0 の初期配置において，界面より 1 モノマー (メチレン基 3 つ) 程度結晶相側の領域が黒く着色されている．t = 0 fs と t = 200 fs の図を比較してわかるように，緩和によってまず最初にこの界面近傍の影響領域に乱れを生じる．その後，楕円で囲った分子鎖において黒い部分
が下方に移動していく様子からわかるように，gauche–gauche のような構造欠陥が分子鎖内を伝ぱしていくメカニズムにより，乱れが界面から結晶内部に伝ぱするのが観察された．
◦

Model 45 の初期構造緩和時の結晶・非晶界面近傍の構造変化を図 4.4 に示す．図◦4.3 と同様に trans 点以外の 2 面角を黒く着色して表示している．Model 0 と比べて，
t = 0 fs の初期構造における界面近傍の結晶相に黒い部分が少なく，結晶⇒非晶と連続する分子鎖ノードの 2 面角変化が少ないことがわかる．これは，斜めに配向した結晶相の隣接分子鎖の影響を受けて，分子鎖が直線状に成長したためである．緩和を開始◦

すると，Model 0 の場合と異なり，界面近傍の領域は比較的その構造を保ったまま，

gauche–gauche の構造欠陥が多数の分子鎖で発生し分子鎖方向に伝ぱした．
◦◦

10000 fs の初期緩和シミュレーション後の Model 0 および Model 45 の分子鎖構造

第4章

結晶・非晶界面における微視的変形機構

t = 0 fst = 200 fs

t = 220 fs

t = 240 fs

t = 260 fs

t = 280 fs

Fig.4.3 Propagation of conformational defects from the crystal/amorphous◦interface during initial relaxation (Model 0 ).

33

第4章

結晶・非晶界面における微視的変形機構

t = 0 fs

t = 60 fs

t = 80 fs

t = 100 fs

t = 120 fs

t = 140 fs

Fig.4.4 Propagation of conformational defects from the crystal/amorphous◦interface during initial relaxation (Model 45 ).

34

第4章

結晶・非晶界面における微視的変形機構

35を図 4.5 に示す．図ではセル全体の分子鎖構造を xz 平面に投影して示しており，また，◦◦◦Model 45 は Model 0 の 1/2 の縮尺で表示している．図より，Model 0 では界面より結晶相側 2.5 nm 程度の領域に乱れが多く認められる．また，図 (a) の結晶相中の右側の乱れは，分子鎖の方向に沿って結晶内部に深く伝ぱしているのがわかる．Model 0◦の非晶相には大きな変化はなく，初期配置で作成した結晶・非晶界面をほぼ保ってい◦る．一方，Model 45 では，界面に接する非晶相が結晶化し，結晶・非晶界面は初期配置で作成したそれと大きく異なっている．結晶化による局所的な界面は，分子鎖方向に対して垂直となる傾向にあり，界面は結晶化した部分の長さがそろった領域が階段状に並んだ構造となっている．また，結晶相内においても，界面ごく近傍には乱れた領域が少なく，界面より数 nm 程度距離をおいて分子鎖に垂直方向に広がった乱れを生じている．結晶の深部においても同様の乱れが多数認められる．ここで，乱れていない結晶部分が，ほぼ同程度の分子鎖長さを有するクラスターとして存在している

Crystallization

Initialinterfaceposition
Defects
Defects(a) Model 0o(b) Model 45o

Fig.4.5 Snapshots of molecular chains after the initial relaxation. Model 45◦
is shown in one-half reduced scale against Model 0 .

◦

第4章

結晶・非晶界面における微視的変形機構

36◦のは興味深い．図 4.4 で述べた，構造欠陥の集団的発生・伝ぱが Model 45 における界面の再結晶化および内部欠陥構造形成をもたらしたものと考えられる．しかしながら，初期緩和過程におけるこれらの構造欠陥伝ぱおよび結晶化のメカニズムは，初期配置に強く依存すると考えられるためこれ以上触れないこととする．次節では，これらの初期平衡構造を有する結晶・非晶界面に対して引張シミュレーションを行った結果について検討する．

4.2.2.1

引張応答および内部構造変化◦

引張シミュレーションにより得られた応力–ひずみ関係を図 4.6 に示す．Model 0 は◦

Model 45 に比べ小さい系であるため応力のゆらぎが大きくなっている．全体的な応◦力変化を見ると，Model 0 は前章の完全非晶体の引張と同じように，ひずみの増加にも関わらずほぼ一定応力で変形が進行する領域 (εzz < 0.05) と，その後応力が急激に◦増加する硬化領域が認められる．一方，Model 45 では全く異なった挙動となり，引張初期から応力がほぼ一定の割合で増加した．
◦

引張時における Model 0 の内部分子鎖構造の変化を図 4.7 に示す．εzz = 0.0 ∼ 0.25までの変化を通じて，結晶・非晶界面の位置がほとんど変化していない．このことか◦ら，Model 0 では結晶相はほとんど変形せず，非晶相で変形を吸収していることがわ
かる．非晶相の変化を詳細に観察すると，εzz = 0.05 までは非晶相はほぼ均等に変化しているのに対し，図中実線長方形で囲ったように，εzz = 0.10 以降は界面近傍の非晶分子鎖が引張方向に配向し大きく延伸している．このとき，非晶相の他の部分には顕◦著な変化はない．すなわち，Model 0 の εzz > 0.05 における応力上昇は，界面近傍の応力集中部において生じた非晶分子鎖の配向および延伸によるものであり，基本的には前章の完全非晶体の引張における硬化機構と同一のものであると結論づけられる．
◦

また，図 4.7 中に実線楕円で囲ったように，Model 0 の結晶相では引張ひずみの増加とともに初期緩和過程で生じた界面近傍の乱れが解消しているのが観察された．図4.7 中に小円で示した領域の変化を拡大して図 4.8 に示す．εzz = 0.0 の図を見ると，結晶の乱れた領域では平行に並んだ直線分子鎖間を斜めに交差する分子鎖が多数認められる．これは，gauche–gauche の構造欠陥により，平面ジグザグ構造からずれて隣接格

第4章

結晶・非晶界面における微視的変形機構

0.05

0.10

Tensile stress,

σzz , GPa

4.0

3.0

2.0

1.0

0

-1.0

0

0.15

Applied strain,(a) Model 0

0.20

0.25

ε zz

◦

Tensile stress,

σzz , GPa

4.0

3.0

2.0

1.0

0

-1.0

0

0.05

0.10

0.15

Applied strain,(b) Model 45

0.20

ε zz

◦

Fig.4.6 Stress–strain curves.

0.25

37

第4章

結晶・非晶界面における微視的変形機構

ε zz = 0.00

ε zz = 0.05

ε zz = 0.10

ε zz = 0.15

ε zz = 0.20

ε zz = 0.25

38◦

Fig.4.7 Snapshots of internal molecular chains in the simulation of Model 0 .

第4章

結晶・非晶界面における微視的変形機構

ε zz = 0.00

ε zz = 0.05

ε zz = 0.10

ε zz = 0.15

ε zz = 0.20

ε zz = 0.25

Fig.4.8 Re-crystallization of the disordered molecular chains under tension◦(Model 0 ).

39

第4章

結晶・非晶界面における微視的変形機構

40

子点にシフトしている分子鎖である．引張ひずみの増加とともに，これらの分子鎖は元の格子点に戻り，結晶の乱れが解消されている．このような，引張ひずみを駆動力として，結晶相内部の構造欠陥が消滅するメカニズムが本シミュレーションにより確認された．
◦

図 4.9 に，Model 45 の引張シミュレーションにおいて観察された内部分子鎖構造変化を示す．結晶・非晶界面の位置関係を明確にするため，図では初期配置における界面◦の位置 (セルの z 方向中央) を実線で示している．Model 0 の場合と異なり，結晶・非晶界面の位置は常にセルの中央付近にあり，ひずみに比例して結晶相も変形していること◦がわかる．これは，引張軸に対して分子鎖方向が 45 傾いているため，van der Waals による弱い結合力で凝集している分子鎖垂直方向に力が作用したためである．Model 45◦の応力ひずみ関係が引張初期から線形的な応力上昇を示したのは，主として結晶相の分子鎖間隔が拡張したことによるものと結論づけられる．
図 4.9 の結晶相内部の変化を詳細に観察すると，εzz = 0.0 において乱れていた部分◦が解消していることがわかる．これは Model 0 で示した構造欠陥消滅と同じメカニズムによるものである．また，εzz = 0.15 近傍から，結晶相下部において分子鎖が引張方向に回転しているが，これは結晶相下端の変位制御部の影響によるものであり，そのメカニズムはここでは議論の対象としない．しかしながら，変形に局所的な強い拘束を受けたときには，分子鎖回転を生じる可能性があることを示す興味深い現象である．さらに，初期平衡状態では階段状であった結晶・非晶界面が，変形後期には引張軸に対してほぼ垂直になだらかな面を形成しているなどの興味深い点があるが，そのメカニズムの解明については今後の検討課題としたい．

第4章

結晶・非晶界面における微視的変形機構

ε zz = 0.00

ε zz = 0.05

ε zz = 0.10

ε zz = 0.15

ε zz = 0.20

ε zz = 0.25◦

Fig.4.9 Snapshots of internal molecular chains in the simulation of Model 45 .
Solid line indicates the initial position of the crystal/amorphous interface.

41

第4章

4.3

結晶・非晶界面における微視的変形機構

42

結言

ラメラ構造など，実際のポリマー材料内部に多く存在する結晶・非晶界面について，
その微視的変形挙動を原子レベルから明らかにするため，結晶相と非晶相の 2 相からなる界面モデルを用いた引張分子動力学シミュレーションを行った．ここで，結晶相◦◦◦の分子鎖が引張方向に平行なもの (Model 0 ) と，45 傾いたもの (Model 45 ) の 2 つのモデルを考慮し，引張軸に対する結晶相方位の影響について検討した．まず，無負荷平衡状態のシミュレーションを 10000 fs 行い，界面の分子鎖構造を詳細に観察した．
その結果，◦(1) Model 0 では，結晶相の界面近傍に gauche–gauche の構造欠陥による結晶の乱れが多く認められた．非晶相に大きな変化はなく，結晶・非晶界面は引張軸に対して垂直な平面を保っていた．
◦(2) Model 45 では，非晶分子鎖の結晶化が界面近傍で認められた．結晶化による局所的な界面は，分子鎖方向に対して垂直となる傾向があり，結晶・非晶界面は鋸歯状となった．
◦◦(3) Model 0 ，45 のいずれにおいても，gauche–gauche の構造欠陥が分子鎖方向に伝ぱし，結晶の乱れを生じることが観察された．また，この乱れは分子鎖に対して垂直方向に広がりを有することがわかった．
その後，得られた界面構造に対して引張シミュレーションを行い，引張負荷時の応答
および内部分子鎖構造変化を詳細に観察した．その結果，◦(4) Model 0 の応力–ひずみ関係は，完全非晶体における引張と同様，ほぼ一定の応力で変形が進行した後，応力が急増するひずみ硬化挙動を示した．一方，Model 45◦のそれは引張初期から線形的な応力増加を示した．
◦(5) Model 0 では，結晶相はほとんど変形せず非晶相が変形を吸収している．非晶相は引張初期にはほぼ均一に変形し，その後界面近傍において非晶分子鎖が引張

第4章

結晶・非晶界面における微視的変形機構

43◦

方向に配向・延伸していた．Model 0 の応力上昇は，このような界面近傍の応力集中部に生じた分子鎖配向・延伸によりもたらされていた．
◦(6) Model 45 では，非晶相だけでなく結晶相もひずみに応じた変形を生じていた．
◦

これは，引張方向に対して分子鎖が 45 傾いているために，van der Waals による弱い凝集力で結合している分子鎖横方向に力が作用したためである．引張初期からの線形的な応力増加は，このような結晶相の分子鎖横方向の変形によるもの
であることが示された．

第5章
結論
本研究では，ポリマー材の変形・破壊時に内部で生じる分子鎖の微視的変形挙動について新たな知見を得ることを目的とし，ポリエチレンを対象として種々の分子動力学シミュレーションを行った．特に，長く不規則な非晶分子鎖の配向や分子鎖間のからみ，結晶・非晶界面領域の構造変化等について詳細に検討した．以下に，得られた結果を総括する．
第 2 章では，本研究で用いた解析手法の基礎について述べた．まず，分子動力学法の概要ならびに基礎方程式を示し，本研究で用いた数値積分法について説明した．次に，粒子間相互作用の評価に用いられるポテンシャルエネルギーについて述べ，ポリエチレンのポテンシャル関数を具体的に説明した．さらに，大規模シミュレーションに必要な計算の高速化手法，および，原子系における応力の評価法について述べた．
第 3 章では，200 万の粒子数からなる大規模ポリエチレンモデルによる引張シミュ
レーションを行い，長く不規則な非晶高分子鎖の複雑な変形挙動を周期境界等の非物理的な要因を排除して検討した．ひずみ速度 ε̇zz = 5.0 × 1011 s−1 の本シミュレーション
条件では，引張ごく初期に線形弾性応答を示した後，ほぼ一定応力で変形が進行する
降伏挙動が見られた．また，εzz = 1.5 以降に再び応力が上昇するひずみ硬化挙動が認められた．系のエネルギー変化の詳細な観察により，線形弾性応答から初期降伏にかけては，(1) 結合長の微小変化による変形吸収，(2) 結合角の微小変化による変形吸収，(3)2 面角の微小変化による応力緩和，というメカニズムがあることが示された．また，
2 面角の gauche ⇔ trans 遷移の詳細な観察から，一定応力下で変形が進行した降伏過

44

第5章

結論

45

程においては gauche ⇒ trans 遷移が優位に生じており，分子鎖が直線状になる内部構
造変化が生じていることが示された．さらに，内部分子鎖構造変化の直接観察ならびに各粒子近傍の局所密度分布の変化により，降伏⇒硬化過程では分子鎖が複雑にから
みあって変形しにくいクラスター状の部分 (からみ点) と，それらを直線状の分子鎖が接続したネットワーク状構造が形成されていることがわかった．これらの事実に対する考察より，試験片上端から下端まで連続した分子鎖が存在しない本シミュレーショ
ンモデルにおいて生じた応力上昇は，分子鎖ネットワーク構造の形成，および，からみ点間分子鎖の配向・延伸によりもたらされたということが結論づけられた．
第 4 章では，ラメラ構造など，実際のポリマー材料内部で多く観察される結晶・非晶界面における微視的変形挙動，および，結晶相分子鎖方位の違いによる影響などを
原子レベルから明らかにするため，結晶相の分子鎖が引張方向に平行な界面モデル◦◦◦(Model 0 ) と 45 傾いた界面モデル (Model 45 ) の 2 つのポリエチレン界面モデルについて引張分子動力学シミュレーションを行った．まず，無負荷平衡状態のシミュ◦

レーションを行い，界面の分子鎖構造を詳細に検討した．その結果，(1)Model 0 では，結晶相の界面近傍に gauche–gauche の構造欠陥による結晶の乱れが多く認められたが，非晶相に大きな変化はなく，結晶・非晶界面は引張軸に対して垂直な平面を保っ◦ている，(2)Model 45 では，界面近傍の非晶分子鎖が結晶化し，引張軸に対して垂直◦◦な平面であった結晶・非晶界面は鋸歯状となる，(3)Model 0 ，45 のいずれにおいても，gauche–gauche の構造欠陥が分子鎖方向に伝ぱし，結晶の乱れを生じる，(4) 結晶の乱れは分子鎖に対して垂直方向に広がりを有する，などが明らかとなった．その後，得られた界面構造に対して引張シミュレーションを行い，引張負荷時の応答，および，◦

内部分子鎖構造変化を詳細に観察した．その結果，(1)Model 0 の応力–ひずみ関係は，完全非晶体における引張と同様に，ほぼ一定応力の変形の後，応力が急増するひずみ◦硬化挙動を示す，(2)Model 45 の応力–ひずみ関係は引張初期から線形的に応力が増◦加する，(3)Model 0 では結晶相はほとんど変形せず，非晶相で変形を吸収している，◦(4)Model 0 の非晶相は，引張初期のほぼ均一な変形の後，界面近傍の非晶分子鎖が引張方向に配向・延伸しており，その硬化機構は第 3 章の完全非晶体のそれと基本的◦◦には同じである，(5)Model 45 では，引張方向に対して分子鎖が 45 傾いているため

第5章

結論

46

に，van der Waals による弱い凝集力で結合している分子鎖横方向に力が作用し，非晶◦相だけでなく結晶相もひずみに応じた変形を生じている，(6)Model 45 における引張初期からの線形的な応力増加は，このような分子鎖横方向変形によるものである，などが明らかとなった．
実際の高分子材料における分子鎖折り畳み等の構造変化は極めて緩慢に生じる．したがって，本研究で観察した内部分子鎖構造変化は，通常のポリマー材料の変形試験において内部で生じている現象とは特に時間的スケールに大きな隔たりがあるのは事実である．分子動力学によるほとんどすべての研究が対峙するこの問題の一つの解決策として，遷移状態理論に基づく分子動力学法が提案されている [35] が，分子鎖変形挙動における時間依存性については今後の検討課題としたい．

参考文献[1] 成沢郁夫，高分子材料強度学，(1982)，オーム社．
[2] Boyce, M. C.，Parks, D. M. and Argon, A. S.，Inelastic Deformation of Glassy
Polymers. Part I：Rate Dependent Constitutive Model，Mech. Mater.，7 (1988)，
15-33．
[3] Wu, P. D. and Van der Giessen, E.，On Improved Network Models for Rubber
Elasticity and Their Applications to Orientation Hardening in Glassy Polymers，
J. Mech. Phys. Solids.，41-3 (1993)，427-456．
[4] Arruda, E. M. and Boyce, M. C.，Eﬀects of Strain Rate，Temperature and Thermomechanical Coupling on the Finite Strain Deformation of Glassy Polymers，
J. Mech. Phys. Solid.，41-2 (1993)，389-41．2[5] 冨田佳宏，田中繁之，からみ点変化を考慮した分子鎖網目理論によるガラス状ポ
リマーの構成式と変形挙動のシミュレーション，日本機械学会論文集，61-581A(1995)，120-126．
[6] Tomita, Y.，Adachi, T.，and Tanaka, S.，Modelling and Application of Constitutive Equation for Glassy Polymer Based on Nonaﬃne Network Theory，Eur.
J. Mech.，A/Solids，16-5 (1997)，745-755．
[7] Abraham, F. F.，Atomic Simulation of Crack Propagation，Materials Science
for the 21st Century; The Society of Materials Science，A (2001)，195-202．
[8] Vashishta, P.，Bachilechner, M. E.，Nakano, A.，Campbell, T. J.，Kalia, R.
K.，Kodiyalam, S.，Ogata, S. and Shimojo, F.，Multimillion Atom Simulationof Mechanical Behavior of Nanostructured Materials，Materials Science for the
21st Century; The Society of Materials Science，A (2001)，213-219．
47

参考文献

48

[9] 小山暁，山本隆，高分子結晶化のモレキュラー・ダイナミクス，日本ゴム協会誌，
75-10 (2002)，430-434．
[10] Shimizu, T. and Yamamoto, T.，Melting and Crystallization in Thin Film ofn-Alkanes: A Molecular Dynamics Simulation，J. Chem. Phys.，113 (2000)，
3351-3359．
[11] Kavassalis, T. A. and Sandararajan, P. R.，Macromolecules，26 (1993)，41444150．
[12] Fujiwara, S. and Sato, T.，Molecular Dynamics Simulations of Structural Formation of a Single Polymer Chain: Bond-Orientational Order and Conformational
Defects，J. Chem. Phys.，107 (1997)，613-622．
[13] Liu, C. and Muthukumar, M.，Langevin Dynamics Simulations of Early-Stage
Polymer Nucleation and Crystallization，J. Chem. Phys.，109 (1998)，25362542．
[14] Yamamoto, T.，Molecular Dynamics Simulation of Polymer Ordering I. Crystallization from Vapor Phase，J. Chem. Phys.，109 (1998)，4638-4645．
[15] Chen, C. M. and Higgs, P. G.，Monte-Carlo Simulation of Polymer Crystallisationin Dilute Solution，J. Chem. Phys.，108 (1998)，4305-4314．
[16] Koyama, A.，Yamamoto, T.，Fukao, K. and Miyamoto, Y.，Molecular Dynamics
Simulation of Polymer Crystallization from an Oriented Amorphous State，Phys.
Rev.，E，65 (5) (2002)，050801(R)-1 - 4．
[17] Madkour, T. M.，A Combined Statistical Mechanics and Molecular Dynamics
Approach for the Evaluation of the Miscibility of Polymers in Good, Poor andNon-solvents，Chem. Phys.，274 (2001)，187-198．

参考文献

49

[18] Aoyagi, T.，Sawa, F.，Shoji, T.，Fukunaga, H.，Takimoto, J. and Doi, M.，A
General-Purpose Coarse-Grained Molecular Dynamics Program，Comp. Phys.
Comm.，145-2 (2002)，267-279．
[19] Pant, P. V. K. and Boyd, R. H.，Simulation of Diﬀusion of Small-Molecule Penetrants in Polymers，Macromolecules，25 (1992)，494-495．
[20] Muller-Plathe, F.，Rogers, S. C. and van Gunsteren, W. F.，Diﬀusion Coeﬃcientsof Penetrant Gases in Polyisobutylene can be Calculated Correctly by Molecular
Dynamics Simulations，Macromolecules，25 (1992)，6722-6724．
[21] Pant, P. V. K. and Boyd, R. H.，Molecular Dynamics Simulation of Diﬀusion of
Small Penetrants in Polymers，Macromolecules，26 (1993)，679-686．
[22] Yoon, D. Y.，Smith, G. D. and Matsuda, T.，A Comparison of the United Atomand Explicit Atom Models in Simulations of Polymethylene，J. Chem. Phys.，98(1993)，10037-10043．
[23] 桑島聖，野間裕之，逢坂俊郎，非晶高分子の MD 計算用力場パラメータについて
II:ポリオレフィン，第 4 回計算化学シンポジウム，東京 (1994)，53-56．
[24] 渋谷陽二，中村浩二郎，冨田佳宏，非晶性ポリマーの単軸引張り変形における自由エネルギの変化，日本機械学会論文集，65-629A (1999)，87-92．

[25] Lennard-Jones, J. E.，Proc. Roy. Soc.，London，ser，A，106 (1924)，463．
[26] Morse, P. M.，Diatomic Molecules According to the Wave Mechanics. II. Vibrational Levels，Phys. Rev.，34 (1929)，57-64．
[27] Daw, M. S. and Baskes, M. I.，Semiempirical, Quantum Mechanical Calculationof Hydrogen Embrittlement in Metals，Phys. Rev. Lett.，50-17 (1983)，12851288．

参考文献

50

[28] Daw, M. S. and Baskes, M. I.，Embedded-Atom Method: Derivation and Application to Impurities, Surfaces, and Other Defects in Metals，Phys. Rev. B，29-12 (1984)，6443-6453．
[29] 上田顯，コンピューターシミュレーション，朝倉書店，(1990)．
[30] Martin, J. W.，Many-Body Forces in Metals and the Brugger Elastic Constants，J. Phys.，C，8 (1975)，2837-2857．
[31] 伊藤智啓，屋代如月，冨田佳宏，非晶性ポリマー内部の分子鎖集団挙動に関する分子動力学解析，日本材料学会第 7 回分子動力学シンポジウム講演論文集，(2002)，

50-55．
[32] 松本龍介，北川浩，中谷敬子，中谷彰宏，アモルファス金属中を進展するき裂の分子動力学シミュレーション，日本機械学会論文集，67-653A (2001)，23-29．

[33] Bunn, C. W.，The Crystal Structure of Long-Chain Normal Paraﬃn Hydrocarbons; the ‘Shape’ of the = CH2 Group.，Trans. Faraday Soc.，35 (1939)，
482-491．
[34] Seto, T.，Hara, T. and Tanaka, K.，Phase Transformation and Deformation
Process in Oriented Polyethylene，Jpn. J. Appl. Phys.，7 (1968)，31-38．
[35] Voter, A. F.，A Method for Accelerating the Molecular Dynamics Simulation of
Infrequent Events，J. Chem.Phys.，106 (1997)，4665-4677．

関連発表論文・講演論文
関連発表論文 屋代 如月・伊藤 智啓・冨田 佳宏，単軸引張下の非晶性ポリマー内部における分子鎖配向とからみ点形成の分子動力学解析，材料，(2003)，(投稿中)．

 Kisaragi YASHIRO・Tomohiro ITOH・Yoshihiro TOMITA，Molecular Dynamics
Simulation on deformation behavior in amorphous polymer: Nucleation of chain
entanglements and network structure under uniaxial tension，International Journal of Mechanical Sciences，(2003)，(submitted)．
 屋代 如月・伊藤 智啓・冨田 佳宏，ポリマーの結晶/非晶界面における分子鎖変形挙動の分子動力学解析日本機械学会論文集，A 編，(2003)，(投稿中)．
講演論文

 伊藤 智啓・屋代如月・冨田佳宏，非晶性ポリマー内部の分子鎖集団挙動に関する分子動力学解析，日本材料学会第 7 回分子動力学シンポジウム講演論文集，(2002)，

50-55．
 伊藤 智啓・屋代 如月・冨田 佳宏，ポリマーの結晶・非晶界面における微視的
変形挙動の分子動力学解析，日本機械学会第 15 回計算力学講演会講演論文集，(2002)，161-162．

51

関連発表論文・講演論文
論文貼り付け用ページ52

関連発表論文・講演論文
論文貼り付け用ページ53

関連発表論文・講演論文
論文貼り付け用ページ

54

関連発表論文・講演論文
論文貼り付け用ページ55

関連発表論文・講演論文
論文貼り付け用ページ

56

関連発表論文・講演論文
論文貼り付け用ページ

57

関連発表論文・講演論文
論文貼り付け用ページ

58

関連発表論文・講演論文
論文貼り付け用ページ

59

関連発表論文・講演論文
論文貼り付け用ページ60

関連発表論文・講演論文
論文貼り付け用ページ61

関連発表論文・講演論文
論文貼り付け用ページ62

関連発表論文・講演論文
論文貼り付け用ページ63

関連発表論文・講演論文
論文貼り付け用ページ64

関連発表論文・講演論文
論文貼り付け用ページ

65

関連発表論文・講演論文
論文貼り付け用ページ

66

関連発表論文・講演論文
論文貼り付け用ページ

67

関連発表論文・講演論文
論文貼り付け用ページ68

関連発表論文・講演論文
論文貼り付け用ページ

69

関連発表論文・講演論文
論文貼り付け用ページ

70

関連発表論文・講演論文
論文貼り付け用ページ

71

関連発表論文・講演論文
論文貼り付け用ページ72

関連発表論文・講演論文
論文貼り付け用ページ73

関連発表論文・講演論文
論文貼り付け用ページ

74

関連発表論文・講演論文
論文貼り付け用ページ75

関連発表論文・講演論文
論文貼り付け用ページ76

関連発表論文・講演論文
論文貼り付け用ページ

77

関連発表論文・講演論文
論文貼り付け用ページ

78

関連発表論文・講演論文
論文貼り付け用ページ

79

関連発表論文・講演論文
論文貼り付け用ページ80

関連発表論文・講演論文
論文貼り付け用ページ

81

関連発表論文・講演論文
論文貼り付け用ページ

82

関連発表論文・講演論文
論文貼り付け用ページ

83

関連発表論文・講演論文
論文貼り付け用ページ

84

関連発表論文・講演論文
論文貼り付け用ページ

85

関連発表論文・講演論文
論文貼り付け用ページ

86

関連発表論文・講演論文
論文貼り付け用ページ

87

関連発表論文・講演論文
論文貼り付け用ページ88

謝辞
本研究を遂行するにあたり，終始懇切丁寧な御指導賜りました冨田佳宏教授に心より感謝の意を表します．研究を進めるにあたり，研究報告会などにおいて様々な御助言を頂いた安達泰治助教授に深く感謝致します．研究全般に渡り直接御指導を頂いたのみならず，学会講演や論文作成において多くの御助言・御協力頂きました屋代如月
助手に深く感謝致します．また，研究環境づくりに御尽力下さった古宇田由夫技官に心から御礼申し上げます．
さらに，常日頃から共に研究に取り組み，御助言，御協力，御討論下さいました固体力学研究室諸氏に深く感謝致します．特に，研究室に配属されて以来，同期の市原彰君，内田真君，梶田文昭君，國本宏通君，西海征志君には研究外においても御助言，御協力頂きました．感謝致します．また，同じ研究グループの金井君，山上君，黒瀬
君，田端君，橋本君には研究活動において多大なご協力を頂きました．御礼申し上げます．
最後に，これまでの 24 年間，家族には精神的にも経済的にも多大な支援を頂きました．ここに記して，最大限の感謝の意を表します．
平成 15 年 1 月
伊藤89

智啓