修士論文

分子動力学法によるアモルファスポリエチレンの押し込み
挙動評価
押し込み速度・架橋点の効果

指導教員： 田中

藤麻

克志

成貴

2015 年 2 月
神戸大学大学院 工学研究科

博士課程前期課程

機械工学専攻Master thesis

Molecular Dynamics Studies
on indentation behavior evaluation of theamorphous polyethylene

Shigeki FUJIMA
Feburary 2015
Department of Mechanical Engineering,
Graduate School of Engineering,
Kobe University, Kobe, Japan

要約
本研究では，高分子の分子鎖の変形挙動を原子レベルから検討する為，アモルファス
ポリエチレン (PE) に対するナノインデンテーションを分子動力学法を用いて行った．
まず，周期境界下でバルクとして作成したアモルファス構造の PE について，一方向の周期境界を外して表面とした無限薄膜を横方向応力 0 と横ひずみ 0 の２種類で緩和を行った後，架橋の有無を考慮しながら球体及びピラミッド型圧子の押し込みシミュ
レーションを行った．横応力 0 では厚さ方向の変形抵抗が小さいため急激に圧縮され，全体の密度が高くなり自由表面側に分子鎖が流動し，周期セルが厚さ方向に細長くなる問題点が明らかになった．押し込みを行った結果，架橋を有する系は含まない系よりも反力が上昇し，特にピラミッド型圧子ではこれが顕著であった．また，内部エネ
ルギーが減少する要因は van der Waals であり，押し込みにより高密度となった分子鎖が流動変形した「粘弾性挙動」によるものであることを明らかにした．横ひずみ 0
では横のセル長は固定されているものの，やはり分子鎖は自由表面方向に流動することや，架橋を含む系ではそれが抑制されることが判明した．同様に押し込みを行うと，これまでよりも高速押し込みを行ったことで，押し込み時は全体のエネルギー変動を支配する van der Waals は増加し，押し込み保持直後減少した．この現象も高速押し込みによる粘弾性挙動が顕著に表れた為発生したことを明らかにした．

Summary
Nanoindentation is a powerful tool to evaluate the mecanical properties in the extremely small region. However, polymeric materials have the complicated chain structure in the scale of nanoindentation, and show viscoelastic behavior, so that the deformation mechanism under indentation is still open for discussion. In order to clarify
the deformation behavior of molecular chains under indenter, nanoindentations are
performed on an amorphous block of polyethylene(PE) by molecular dynamics simulations. First, bulk amorphous structure is made in a cubic periodic cell by growing
random PE chains freely. Then the periodicity in the z-direction is removed and relaxed during 100,000fs controlling horizontal stress 0. Here we set two PE ﬁlms, one
is the original, the other adding the chain bridging randomly. However, the lateralstress control shrink the parallelepiped periodic cell endlessly, due to the gap between
very thin thickness and inﬁnite width in the lateral direction. Thus we stopped thestress control and performed indentation by a spherical and pyramid indenter of rigid
diamond on to the amorphous PE ﬁlm surface. The potential energy of the PE chainsdecreases in the later stage of the indentation, while it increases during the tip holding
at the maximum depth. This strange behavior can be explained by the change of thevan der Waals potential (VDW). That is, the VDW is certainly increased and gives
repulsive forces between near-miss two particles under high compression; however, thetotal of VDW is decreased since many far-ﬂung particles come into the VDW interaction. During the holding process, the repulsive motion of near-miss particles spreads
the chains, so that the VDW increases without external work. It should be notedthat we revealed the role and mechanism of the VDW in the viscoelastic motion ofpolymeric materials. We also made inﬁnite thin PE ﬁlm under the condition of thelateral strain ϵ =0 and performed indentation in the same manner. The indentation
response is changed since the chain structure of the surface is diﬀerent, e.g. diﬀerent
density and entanglements; however, the delayed response by the VDW is still key
mechanism.

目次
第1章 緒 論3

第 2 章 解析手法の基礎

62.1

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

62.2

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

7

2.3

架橋点のモデル化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15

2.4

高速化手法

16

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

第 3 章 横応力 0 条件下で作成した薄膜への押し込みシミュレーション

3.1

3.2

3.3

横応力 0 条件下での薄膜の作成 . . . . . . . . . . . . . . . . . . . . . .

18

3.1.1

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

18

3.1.2

初期構造緩和結果及び考察 . . . . . . . . . . . . . . . . . . . . .

22

球体圧子での押し込み . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

3.2.1

押し込み反力の変化 . . . . . . . . . . . . . . . . . . . . . . . .

28

3.2.2

内部エネルギー変化 . . . . . . . . . . . . . . . . . . . . . . . .

29

3.2.3

ポテンシャル成分の変化 . . . . . . . . . . . . . . . . . . . . . .

31

3.2.4

局所密度及び２面角の変化 . . . . . . . . . . . . . . . . . . . . .

34

ピラミッド型圧子での押し込み . . . . . . . . . . . . . . . . . . . . . .

38

3.3.1

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

38

3.3.2

押し込み反力の変化 . . . . . . . . . . . . . . . . . . . . . . . .

40

3.3.3

アモルファスの内部エネルギー変化 . . . . . . . . . . . . . . . .

41

3.3.4

局所密度変化 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

46

第 4 章 横ひずみ 0 条件下で作成した薄膜への押し込みシミュレーション

4.1

18

48

横ひずみ 0 条件下での薄膜の作成 . . . . . . . . . . . . . . . . . . . . .

48

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

48

4.1.1

i

4.1.24.2

目次ii

初期構造緩和結果及び考察 . . . . . . . . . . . . . . . . . . . . .

49

球体及びピラミッド型圧子の押し込み. . . . . . . . . . . . . . . . . .

55

4.2.1

押し込み反力の変化 . . . . . . . . . . . . . . . . . . . . . . . .

55

4.2.2

内部エネルギー変化 . . . . . . . . . . . . . . . . . . . . . . . .

56

4.2.3

ポテンシャル成分の変化 . . . . . . . . . . . . . . . . . . . . . .

58

4.2.4

局所密度及び２面角の変化 . . . . . . . . . . . . . . . . . . . . .

61

第5章 結 論

64

参考文献

66

学術論文・学術講演

69

謝辞

73

第1章
緒論
近年，ナノテクノロジーの進歩に伴い，新しい機能特性を有する材料の研究・開発が盛んにおこなわれている．中でも，半導体デバイスや微小電子機械システム (MEMS)などの設計・開発では，精密機械が微小化・極薄化になるに従い，その表面近傍における微小領域での力学特性評価の必要性が増している [1] ．この微小領域における局所的な力学特性の実験的評価法として，ナノインデンテーション試験が挙げられる．ナ
ノインデンテーション試験は，局所の力学特性を容易に評価できるという利点があり，高分子の分野でも数々の適用例がある．例えば，AFM（原子間力顕微鏡）を用いたナ
ノインデンテーションでは，µN から pN レベルの荷重域での押し込み試験が可能なだけでなく，形成された圧痕形状を直接観察することができる．Jee らは，これを用いて Berkovich 圧子の押し込みによるポリエチレンやポリ塩化ビニルなど様々な高分子の機械特性を測定している[2] ．また，この AFM を用いて人工膝関節置換術で用いられる超高分子量ポリエチレンの硬度の検証もされている [3] ．このほかにも，Loubet らは，同様に Berkovich 圧子の押し込みで，低密度ポリエチレンの表面の線形粘弾性を測定している [5] ．中上らはナノインデンテーション法による金属や高分子の薄膜の硬さ，ヤング率の評価・解析の測定方法と，その適用例について紹介している [4] ．
一方，仮想的な原子レベルの現象観察手法として分子動力学法 (MD) が注目され，高分子材料への適用も古くからみられている．MD を用いたポリエステル層状構造作成[6]や，折り畳み分子鎖によるポリマー結晶のモデル作成 [7] ，ポリエチレン結晶の薄膜モデル形成 [8] など数多くのシミュレーションが行われている．また，数モノマーを 1
ユニットとした粗視化 MD では大規模な系を長時間シミュレーションすることができ，

3

第1章

緒論

4

表面解析やエネルギーの研究 [9][10][11] や，ガラス転移点 [12][13][14] ，Langmuir-Blodgett
膜 [15] や高分子の溶融 [16] の研究のほか，架橋の効果 [17][18] ，reptation[19] についての研究も行われている．
増渕 [20] のレビューによれば，高分子系のミクロシミュレーションには，高分子を構成する粒子間の相互作用を考える「エンタルピー」ベースのシミュレーション (本研究で用いる粗視化分子動力学法もその一つである) と，高分子の構成粒子の種類を問題とせず，ひもとみなした分子鎖の長さ (セグメント長さ) と絡み点のみで高分子の挙動を表現する「エントロピー」ベースのシミュレーションに分けられる．
「エンタルピー」によるシミュレーションは生体高分子の構造解析などの化学的反応が主となる現象には重要であるが，高分子やナノ粒子の動きは非常に複雑かつ緩慢なので，実験観察可能
な時間スケールでのシミュレーションが困難で大きな系への適用も難しい．このため，
エンタルピーベースでの「メカニカル」な変形破壊シミュレーションは現実とのギャッ
プから敬遠されてきたが，素過程を原子レベルで議論する場合，単純かつむしろ極端な条件のシミュレーションの方がメカニズムの解明のヒントとなることが多い [21] ．このような信念のもとで我々の研究グループでは粗視化分子動力学によるポリエチレン・
ポリブタジエンの「メカニカル」なシミュレーションを早くから行ってきた．[22][23][24]
アモルファスポリエチレン (PE) へのナノインデンテーションを行った例では，押し込みを与えているにもかかわらずエネルギーが減少するエントロピー弾性的な挙動を示す事を報告した [25] ．エントロピーベースのシミュレーションは当然エントロピー弾性を再現可能であるが，現象論的に導入されたからみ点は，エンタルピーベースのシ
ミュレーションではその存在が立証できない．従って，エンタルピーベースのシミュ
レーションでエントロピー弾性的挙動を追求する事は，両者を関係付ける架け橋となる可能性を秘めている．
本研究では，前述のナノインデンテーションの更なる考察を目的とし，横方向無限周期境界条件下のポリエチレン薄膜に対して，球体・三角錐・平板の圧子を押し込む
押し込みシミュレーションを行う．架橋を導入した場合や押し込み速度・押し込み深さを変化した場合のシミュレーションも行い，得られた内部エネルギーの変化を第 2章で示す各ポテンシャルエネルギーの変化から考察するとともに，鎖状高分子の 2 面角の変化や局所密度の観点から検討する．

第1章

緒論

5

第 2 章では解析手法の基礎として，分子動力学法を簡単に説明し，分子動力学計算で最も重要となるポテンシャルエネルギーについて述べる．また，大規模シミュレー
ションに必要な計算の高速化手法について示す．
第 3 章では，架橋を含む系と含まない系それぞれの薄膜を作成し横応力 0 条件下で構造緩和を施した後その内部構造を解析する．また，球体及びピラミッド型圧子による押し込みシミュレーションを行い，押し込み反力，内部エネルギーの変化，局所密度および 2 面角の変化から分子鎖の構造変化を議論する．
第 4 章では，3 章と同様に作成した薄膜を横ひずみ 0 条件下で構造緩和を施し解析をした後圧子の押し込みを行い，3 章と同様の観点から分子鎖に与える影響を考察する．
最後に, 第 5 章で本研究の総括を述べる．

第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 の方法，予測子–修正子法等がよく用いられる [27] ．
本研究では，以下に示す Verlet の方法を用いた．
時刻 t + ∆t と t − ∆t での粒子 i の位置ベクトル r i (t ± ∆t) を taylor 展開すると，

dr i (t) (∆t)2 d2 r i (t) (∆t)3 d3 r i (t)r i (t + ∆t) = r i (t) + ∆t
+
+
+・
・
・
dt
2
dt2
3!
dt3
dr i (t) (∆t)2 d2 r i (t) (∆t)3 d3 r i (t)
+
−
+・
・
・
r i (t − ∆t) = r i (t) − ∆t
dt
2
dt2
3!
dt3

(2.3)
(2.4)

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

dr i= v i (t)
dt
6

(2.5)

第2章

解析手法の基礎

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

(2.6)
(2.7)

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

F i (t)+・・・mi(∆t)3 d3 r i (t)r i (t + ∆t) − r i (t − ∆t) = 2∆tv i (t) + 2
+・
・
・
3!
dt3
r i (t + ∆t) + r i (t − ∆t) = 2r i (t) + (∆t)2

(2.8)
(2.9)

3

が得られる．これより，(∆t) 以上の高次項を無視すると，時刻t + ∆t での位置ベク
トルと t での速度は

r i (t + ∆t) = 2r i (t) − r i (t − ∆t) + (∆t)2v i (t) =F i (t)mi1{r i (t + ∆t) − r i (t − ∆t)}2∆t(2.10)(2.11)と求められる．これより，時刻 t + ∆t での座標は，時刻 t と t − ∆t での座標，ならびに時刻 t における力 F 0 がわかれば計算できる．最初の計算 (t = 0) では，t = ∆t での座標 r i (∆t) は式 (2.6) と初速度から決定される．これより r i (∆t) と r i (0) が既知となり，式 (2.10) を繰り返し適用することにより各粒子の座標を更新していく．

2.2

原子間ポテンシャル

粒子に作用する力は系のポテンシャルエネルギーにより決定される．従って，系のポテンシャルエネルギーの評価が分子シミュレーションにおいて重要となる．本研究で扱うポリマー材料は強い共有結合と弱い共有結合からなり，また共有結合には分子鎖内部の結合角度や回転等の内部自由度があるため，ポテンシャルエネルギー Etot は
次式のように表される．

Etot = ΦBS (r) + ΦBE (θ) + ΦTO (ϕ) + ΦVW (r̄)(2.12)

右辺各項は，分子内の炭素間の結合長 r，結合角 θ，2 面角 ϕ の共有結合に対するポテ
ンシャルと，異分子鎖間および同一分子鎖内の十分離れた非共有原子間距離 r̄ に対す

第2章

解析手法の基礎8

る van der Waals ポテンシャルを表す．本研究で解析対象とするポリエチレンに対しては，水素原子の運動は陽に扱わず，メチレン炭素 (-CH2 -) を単一粒子として粗視化する united atom model により，以下の関数形及びパラメータが提案されている [26] ．

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

kr (r − r0 )2}(2.13)

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

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

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

各関数のパラメーターの値を表 2.1〜2.4 に，ポテンシャル曲線を図 2.1〜2.4 にそれぞれ示した．なお，式中の総和は対象とする高分子の全ノードについて行うが，ダイヤ
モンド圧子内部の炭素原子は剛体として扱うため MD 計算を行わない．van der Waals
ポテンシャルの計算は異なる分子鎖の粒子間，4 原子団以上離れた粒子間，C 原子と-CH2 -に対して行う．また，今回の押し込みシミュレーションにおける圧子-粒子間は，
van der Waals ポテンシャルの斥力項のみを用いている
式 (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) =∑1{2

kr (r − r0 )2}

Table 2.1 Potential parameter for bond stretch.
r0kr[nm] [kJ/ (mol · nm2 )]0.15331.373 × 105

CH2 – CH2

Φ BS [eV]

100

50

r
0

0.2

r

0.4

[nm]

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

9

第2章

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

ΦBE (θ) =∑1{2kθ (θ − θ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

180

360

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

10

第2章

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

ΦTO (ϕ) =∑(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]100φgauchegauche

-10

trans

0

180

360

φ [deg]Fig.2.3 Relationship between torsion potential ΦTO and torsion angle ϕ.

11

第2章

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

ΦVW (r̄) =

∑

{

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

}

Table 2.4 Potential parameter for van der Waals.

CH2 – CH2

A
[kJ/ (mol · nm12 )]2.972 × 1019

V
[kJ/ (mol · nm6 )]6.907 × 109

ΦVW [eV]

0.01

r0

4

-r [nm]6r-

8

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

12

第2章

解析手法の基礎13

2 体間相互作用r = r αβΦ′ (r) αβrr
F β = −F α

Fα = −(2.17)(2.18)

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

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

4 体間相互作用
r αγ · r γηr γη|r γη |2r βη · r γηAβ = r βη −r γη|r γη |2
Aα · Aβcos ϕ =
|Aα | AβAα = r αγ −Φ′ (ϕ)
1F =sin ϕ |Aα | Aβ
α

Fβ =

′

{

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

}

β





(2.22)Φ (ϕ)
A ·A β
1
α
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.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 bendig angle θ.

j
rjl

φ

k

Ai rik

Aj

rkl
l

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

14

第2章

2.3

解析手法の基礎

15

架橋点のモデル化

先に述べた分子鎖内の共有結合（bond stretch,bending,torsion）および非共有結合（van der Waals）ポテンシャルだけでなく，異なる分子鎖間で近接した 2 粒子を無作為に選びだし，bond stretch と同様の 2 体間ポテンシャル

ΦCL (r′ ) =∑1{2

kr′ (r′ − r0′ )2}(2.26)

で接続することで架橋点をモデル化する．ばね定数 kr′ を弱く設定した場合，引張時に架橋点間距離が開き他の分子鎖がすり抜けてしまうため，変形途中に架橋点間が開かないように bond stretch の約 4 分の 1 とした．安定長 r0′ は van der Waals の LJ ポテ
ンシャルにおける安定位置と同じ値とした．パラメーターの値を表 2.5 に，ポテンシャ
ル曲線を図 2.8 に示す．

Table 2.5 Potential Parameter for crosslink.

CH2 – CH2

r0′
[nm]
0.45

kr ′
[kJ/ (mol · nm2 )]3.474 × 104

ΦCL [eV]

200

100

0

r'

0.4

0.5

r'

0.6

[nm]

Fig.2.8 Relationship between crosslink potential ΦCL and length r′ .

第2章

2.4

解析手法の基礎

16

高速化手法

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

rfc

rc

Fig.2.9 Schematic of bookkeeping method.

第2章

解析手法の基礎

17

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

y
by
x

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

第3章
横応力 0 条件下で作成した薄膜への押し込みシミュレーション
3.1

横応力 0 条件下での薄膜の作成

3.1.1

シミュレーション条件

図 3.1(a) に示すように，14nm × 14nm × 14nm の立方体セル中に，全方向周期境界条件の下で乱数を用いて多数のランダムコイル状分子鎖を成長させ，ポリエチレンア
モルファス構造を作成した．図 3.2 に模式的に示すように，結合長 r および結合角 θ は最安定値 0.1533nm と 113.3◦ に固定した上で，二面角は最安定な ϕ = 180◦ (trans 点) ま
たは準安定な ϕ = ±67.5◦ (gauche 点) のいずれかランダムに設定することで不規則に分子鎖成長させている．分子鎖は周期境界を越えて自由に成長できる条件としており，
図 3.1(a) は周期境界を考慮して分子鎖を１つのセルに折り畳んで表示したもの，(b) は折り畳まずに表示したものである．総粒子数は 100,000 に固定し，分子鎖は 100 本，分子鎖長さは 1000CH2 に統一している．作成時の密度は ρ=0.85g/cm3 となる．
まず，得られた初期構造に対し，全方向周期境界境界条件のもとで 10,000fs の初期緩和計算を行い内部構造緩和を行った．数値積分には Verlet 法を用い，積分の時間ス
テップは 0.1fs とした．温度は 300K とし，速度スケーリングにより制御した．
架橋の効果を議論する為，第 2 章でモデル化した架橋ポテンシャルを用いて，近接
した 2 粒子をランダムに選び出し 20ヶ所架橋を導入したモデルと，作成したそのままのアモルファスポリエチレンの 2 つの系について解析を行った．いずれも，押し込み

18

第3章

横応力 0 条件下で作成した薄膜への押し込みシミュレーション

19を行う上面を自由表面とするために，バルク状態である全方向周期境界から，z 方向上面の周期境界を越えた分子鎖をブロック内部に折り畳んだ後，圧子を薄膜上部に配置する為に z 方向の空間領域を 2 倍にした．その後，横方向周期境界条件で，底面は粒子が越えられない仮想壁の条件下で 100,000fs の緩和計算を行った．
球体圧子を押し込むシミュレーションは，図 3.3 に示すような直径 11nm のダイヤモ
ンド球体圧子で行った．押し込み開始時のスナップショットを図 3.4 に示す．圧子の C
原子数は 9786 であり，先に述べたように圧子は剛体として扱うため，シミュレーショ
ン中は C 原子の運動 (=圧子の変形) は考慮しない．圧子の押し込みは変位制御で行い，
図 3.5 に示すような 2 つの押し込みプロファイルで行った．いずれも，(1) 圧子先端の最大押し込み深さが 4nm となるまで一定速度で押し込み，(2) 圧子を保持，することにより押し込みを行う．押し込み速度は，図 3.5(a) の遅い押し込み (40m/s) と，(b) の速
い押し込み (100m/s) で行った．温度は構造緩和時と同様 300K とし，速度スケーリン
グにより制御した．緩和計算では横方向応力を０に制御したが，押し込みシミュレー
ションではセル寸法を固定する条件とした．

m

14n

14nm

14nm

m
14n

14n
m

z
y

14n

m

x

(a) folded image under periodic cell.

(a) Unfolded image of PE chains.

Fig.3.1 Dimensions of simulation cell.

第3章

横応力 0 条件下で作成した薄膜への押し込みシミュレーション

n=b×agauche
o
67.5

trans
180

o

b
a
gaucheFig.3.2 Conformation of dihedral angles.

11.0nm10.0nmLoading direction

Fig.3.3 Dimensions of diamond sphere indenter.

20

第3章

横応力 0 条件下で作成した薄膜への押し込みシミュレーション

21

Fig.3.4 Snapshots of PE chains in amorphous block and C atoms of diamond
sphere indenter.

Loading Hold

Hold

Indentation depth , d , [nm]Indentation depth , d , [nm]Loading
4

3

2

1

0
0

50000

100000150000

4

3

2

1

0

200000050000

100000150000Time , t , [fs]Time , t , [fs]（a) v = 40m/s（b) v = 100m/s

Fig.3.5 Load waveform.

200000

第3章

3.1.2

横応力 0 条件下で作成した薄膜への押し込みシミュレーション

22

初期構造緩和結果及び考察

架橋を含まない系の 100,000fs の初期緩和中のエネルギー変化を図 3.6(a) に，最初の急激なエネルギー変化を外して 1,000fs 以降のエネルギー変化を拡大したものを図3.6(b) に示す．上面を自由表面にし分子鎖を折り畳んだ無限薄膜にした際，表面形成による高エネルギーが急激に緩和され 1,000fs までは大きく振動した．その後，エネ
ルギーは安定しているように見えるが，1,000fs 以降の拡大図 3.6(b) ではエネルギーは一旦減少した後，緩やかに上昇している事が確認される．ポテンシャルエネルギーの各成分を，1,000fs における値からのエネルギー変化量を縦軸にとり図 3.7 に示す．図3.7 より，100,000fs 経過しても van der Waals 及び torsion が安定値を取っておらず，緩和は終了していないことがわかる．van der Waals は 20,000fs 近傍で底打ちし，以降わずかながら上昇し続けている．torsion は，これもエネルギーの上昇がみられるが，緩和後期につれて上昇は鈍化している．周期セルの底面積の変化を図 3.8(a) に，薄膜の厚さ変化を図 3.8(b) に示す．底面積の変化をみると，緩和直後から横方向の応力制
御により面積の急減がみられるが，20,000fs 近傍から減少が緩やかになっている．図3.8(b) の厚さ変化のグラフでは，z 方向へ増加していることがわかる．図 3.8(c) に系全体の密度変化のグラフを示す．緩和初期に全体の密度が上昇するが，これは図 3.8(a)に示したように横方向応力制御のため急激に圧縮されたことによる．その後緩和するに従い密度が低下しており，図 3.8(b) の高さ変化とあわせて考えると自由表面側に分子鎖が流動していることが示唆される．図 3.9(a) に，各粒子位置においてカットオフ
半径内の粒子数から評価した局所密度 ρ のノード変化を示す．各時間ステップにおけ33る局所密度を，ρ>1.00g/cm の高密度，0.92g/cm <ρの低密度，その間と３種類に分けてカウントしている．なお，1, 000fs の時の値を基準としてそこからの変化量で示している．緩和初期に高密度を示す赤色のグラフが上昇し，密度が疎な青色のグラフは減少している．その後，20, 000fs 近傍で高密度粒子の数が負になり，以降全体的に密度が減少していることが見て取れる．図3.9(b) は全ての 2 面角ノードを調べ，その角度が 175◦ <ϕ<180◦ のものを trans，65◦ <ϕ<70◦ のものを gauche，それ以外を others として数をカウントし，変化を調べたものである．大きな変化は認められないが，最安定点の trans がゆるやかに減少し，準安定点の gauche 及びエネルギー的に不安定な others

第3章

横応力 0 条件下で作成した薄膜への押し込みシミュレーション

23が増加の傾向を示している．

[x103][x103]
Potential energy , E , [eV]

Potential energy , E , [eV]

-0.03

-0.04

-0.05

-0.06

-0.07

-0.08
0

-0.078-0.0785-0.079-0.0795-0.08

20000 40000 60000 80000 100000

020000 40000 60000 80000 100000

Time , t , [fs]Time , t , [fs]（a) Throughout the simulation.

（b) Magnification.

㻱㼠㼛

100

㻱㼎㼑
0

㻱㼎㼟
㻱㼢㼣

-100
0

20000 40000 60000 80000 100000

Time , t , [fs]Change in potential energy , E-E0 , [eV]Change in potential energy , E-E0 , [eV]Fig.3.6 Relationship between potential energy and time (without crosslink).

100

0

-100
0

20000 40000 60000 80000 100000

Time , t , [fs]Fig.3.7 Change in components of potential energy (without crosslink).

横応力 0 条件下で作成した薄膜への押し込みシミュレーション

24

15.2

1.04

20000

19500

Density , [g/cm3]

Thickness , [nm]

15
14.8
14.6
14.4
14.2

1.02

1

0.98

0.96
14

19000

0

20000 40000 60000 80000 100000

20000 40000 60000 80000 100000

020000 40000 60000 80000 100000

Time , t , [fs]Time , t , [fs]Time , t , [fs](a) Cross sectional area.

(b) Thickness at z directon.

(c) Density of film.

Fig.3.8 Change in the crosssection of periodic cell, thickness and density(without crosslink).

10000

0.92<ρ<1.0

5000

0

ρ<0.92

-5000

1.0>ρ
-10000Change in number of node , N-N00

Change in number of node , N-N0

Cross sectional Area , LxxLy , [nm2]

第3章

10000

㟷 㻩gauche㉥ 㻩trans㯮 㻩others

5000

0

-5000

-10000020000 40000 60000 80000 100000

020000 40000 60000 80000 100000

Time , t , [fs]Time , t , [fs](a) Number of density nodes(b) Number of tosional nodes

Fig.3.9 Change in local density and tosional nodes (without crosslink).

第3章

横応力 0 条件下で作成した薄膜への押し込みシミュレーション

25

架橋を導入した系の 100,000fs の初期緩和中でのエネルギー変化を図 3.10(a) に，1,000fs
以降のエネルギー変化を拡大したものを図 3.10(b) に示す．全体的な変化は先の架橋の無い系と同じであり大きな差は見られない．ポテンシャルエネルギーの各成分を，1,000fs における値からのエネルギー変化量を縦軸にとり図 3.11 に示す．架橋なしの系と異なり，van der Waals の緩和後期のエネルギー上昇は抑えられているが，それ以外に顕著な差はない．周期セルの底面積の変化を図 3.12(a) に，厚さ変化を図 3.12(b) に示す．底面積の変化を見ると，架橋未導入とは異なる変動をしている．表面形成直後の横方向の急激な応力制御に対し，一度反発してから底面積が減少している．これは，架橋によってアモルファスの剛性が上昇したためと考えられる．しかしながら，長時間の緩和中に横方向の応力制御のために，自由表面方向に流動変形する傾向は同じである (図 3.12(b))．図 3.12(c) に全体の密度変化，図 3.13 に局所密度と 2 面角のノード変化を示す．架橋による反発で緩和ごく初期は架橋なしのものと異なるが，緩和後期の傾向には大きな差はない．
以上をまとめると，「無限の」厚さを持つ横方向の応力を 0 とするために，圧縮変形が与えられ続けている状態となっている．このまま緩和を継続しても極めて細い棒状のセルとなって押し込みシミュレーションができない可能性があるので，以後の計算では 100,000fs 計算後の状態について押し込みシミュレーションを行う．

第3章

横応力 0 条件下で作成した薄膜への押し込みシミュレーション[x103][x103]-0.0785-0.03

Potential energy , E , [eV]

Potential energy , E , [eV]

26

-0.04

-0.05

-0.06

-0.07

-0.08
0

-0.079-0.0795-0.08-0.0805020000 40000 60000 80000 100000

20000 40000 60000 80000 100000

Time , t , [fs]Time , t , [fs]（a) Throughout the simulation.

（b) Magnification.

200

㻱㼠㼛

100

㻱㼎㼑
0

㻱㼎㼟
-100
0

20000 40000 60000 80000 100000

Time , t , [fs]㻱㼢㼣Change in potential energy , E-E0 , [eV]Change in potential energy , E-E0 , [eV]Fig.3.10 Relationship between potential energy and time (with crosslink).

200

100

0

-100
0

20000 40000 60000 80000 100000

Time , t , [fs]Fig.3.11 Change in component of potential energy (with crosslink).

横応力 0 条件下で作成した薄膜への押し込みシミュレーション

27

15.2

1.04

20100

20000

19900

15

Density , [g/cm3]

Thickness , [nm]

20200

14.8
14.6
14.4
14.2

1.02

1

0.98

0.96
19800

14
20000 40000 60000 80000 100000

020000 40000 60000 80000 100000

Time , t , [fs](a) Cross sectional area.

020000 40000 60000 80000 100000

Time , t , [fs]Time , t , [fs](b) Thickness at z directon.

(c) Density of film.

Fig.3.12 Change in the crosssection of periodic cell, thickness and density(with crosslink).

20000

0.92<ρ<1.010000

ρ<0.920-10000Change in number of node , N-N00

Change in number of node , N-N0

Cross sectional Area , LxxLy , [nm2]

第3章

20000

㟷 㻩gauche㉥ 㻩trans㯮 㻩others

10000

0

-100001.0>ρ
-20000-20000020000 40000 60000 80000 100000

Time , t , [fs](a) Number of density nodes.

020000 40000 60000 80000 100000

Time , t , [fs](b) Number of tosional nodes.

Fig.3.13 Change in density and tosional nodes (with crosslink).

第3章

横応力 0 条件下で作成した薄膜への押し込みシミュレーション

3.2

球体圧子での押し込み

3.2.1

押し込み反力の変化

28

球体圧子の押し込みシミュレーションにより得られた，押し込み力 R と押し込み時
間 t の関係を，時間スケールを揃えて図 3.14 に示す．押し込み力 R は圧子の C 原子が受ける力の z 方向成分により評価した．圧子の押し込みに伴い反力は増加するが，押し込み速度が遅い (b) は，速度が速い (a) と比べて反力の大きさは小さい．圧子を保持すると，(a)(b) 共に反力は減少するが，両方とも 0.02µN 程度の反力を保持している．
架橋による差は，速い押し込みの (a) ではほとんどないように見えるが，架橋がある方がわずかながら高い値を取っている．遅い押し込みの (b) ではその傾向が明白に認められる．これは，架橋点の効果によりブロック内部の構造が変形しにくいからと考えられる．押し込み保持時でも架橋を含む方の反力の減少は少なく，押し込み速度が速い (a) よりも同様の架橋による効果ははっきりとみられている．

Loading

0.08

Without crosslink

0.06

0.04

Crosslink
0.02

0
0

50000

100000150000Time , t , [fs](a) v = 100m/s

200000

Repulsive force , R , [µN]

Repulsive force , R , [µN]

Loading Hold

Hold

0.08

0.06

0.04

0.02

0
0

50000

100000150000200000Time , t , [fs](b) v = 40m/s

Fig.3.14 Relationship between indentation force and time (spherical indenter).

第3章

3.2.2

横応力 0 条件下で作成した薄膜への押し込みシミュレーション

29

内部エネルギー変化

押し込み時の全ポテンシャルエネルギーの時間変化を図 3.15 及び図 3.16 に示す．それ
ぞれ，架橋未導入のグラフが (a)，架橋導入のグラフが (b) であり，図 3.15 は v=100m/s，
図 3.16 は v=40m/s のものである．いずれも，初期緩和計算 100,000fs 時の値からの変化量 (E -E0 ) をとって示している．ただし，押し込み-保持過程では横方向応力制御は行っていない．100000fs 緩和時に横方向の応力制御を外すとエネルギー上昇が止まる
ことを確認しているため，押し込み時には先に述べた横応力制御によるエネルギー上昇の寄与はない．
図中の実線は押し込み初期ならびに圧子保持時のエネルギー変化を最小二乗フィッ
ティングした勾配である．いずれのグラフも押し込み開始時においてエネルギーの上昇が確認される．特に，架橋を含まない v=100m/s の系では他の系よりも実線の勾配よりが大きい．架橋を含まない系はエネルギーの上昇は含む系よりも小さいものとなっている．一方，v=40m/s の遅い押し込みでは架橋の有無にかかわらずエネルギーの上
昇幅は変わらない．その後，いずれの系も押し込みを継続していくとエネルギーが減少に転じる．圧子保持時にはエネルギーが上昇していること，また押し込み中期にエネルギーが減少していることについては次節で考察する．

横応力 0 条件下で作成した薄膜への押し込みシミュレーション[x103]Loading Hold
0.001

0.00050-0.0005-0.001-0.0015

050000

100000150000200000Time , t , [fs]Change in potential energy , E-E0 , [eV]Change in potential energy , E-E0 , [eV]

第3章

30

[x103]Loading Hold
0.001

0.00050-0.0005-0.001-0.0015

050000

100000150000200000Time , t , [fs](a) Without crosslink(b) With crosslink[x103]

Loading

Hold

0.001

0.00050-0.0005[x103]

Loading

Hold

0.001

0.00050-0.0005-0.001-0.0015

Change in potential energy , E-E0 , [eV]Change in potential energy , E-E0 , [eV]Fig.3.15 Change in total energy under spherical indenter, (v=100m/s).

0

50000

100000150000Time , t , [fs](a) Without crosslink

200000-0.001-0.0015

050000

100000150000Time , t , [fs](b) With crosslink

Fig.3.16 Change in total energy under spherical indenter, (v=40m/s).

200000

第3章

3.2.3

横応力 0 条件下で作成した薄膜への押し込みシミュレーション

31

ポテンシャル成分の変化

各ポテンシャルエネルギー成分の時間変化を図 3.17〜図 3.20 に示す．それぞれ，

v=100m/s で架橋なしを図 3.17，架橋ありを図 3.18，v=40m/s で架橋なしを図 3.19，架橋ありを図 3.20 としている．いずれも，初期緩和計算 100,000fs 時の値からの変化量(E -E0 ) をとって示している．
いずれの系も，全ステップを通して bending と bond stretch のエネルギー変化はほぼ 0 である．押し込み後期のエネルギー減少は van der Waals によるものであることがわかる．式 2.4 および図 2.4 に示した２粒子間の van der Waals エネルギーには平衡距離が存在するものの，粒子が極めて接近した状態にならない限り，系の van der Waals
エネルギーは相互作用距離内に存在する粒子が多ければ多いほど大きな負の値となる．
具体的には，最も近い粒子同士が反力を感じる距離より高密度に圧縮されると，その粒子間には強い反力を生じるが，それより遠方の粒子が van der Waals 半径内に入ってくることにより系の van der Waals エネルギーは減少する．したがって，van der Waals
エネルギーの低下は押し込みによって高密度となったためであり，これが押し込み時の反力を担っているものと考えられる．圧子保持時には反力が低下していることから，
高密度となった分子鎖が流動変形している「粘弾性挙動」によるものと考えられる．図3.17 および図 3.18 の一見奇異にみえるエネルギー変化は，このような高分子特有のメカニズムで説明される．

横応力 0 条件下で作成した薄膜への押し込みシミュレーション

Loading Hold
100

㻱㼠㼛
50

㻱㼎㼑
0

㻱㼎㼟

-50

-100
0

50000

100000150000200000㻱㼢㼣Time , t , [fs]Change in potential energy , E-E0 , [eV]Change in potential energy , E-E0 , [eV]

第3章

32

Loading Hold
100

50

0

-50

-100
0

50000

100000150000200000Time , t , [fs]Loading Hold
100

㻱㼠㼛
50

㻱㼎㼑

0

㻱㼎㼟

-50

-100
0

50000

100000150000Time , t , [fs]200000㻱㼢㼣Change in potential energy , E-E0 , [eV]Change in potential energy , E-E0 , [eV]Fig.3.17 Change in total energy under spherical indenter (without crosslink),(v=100m/s).

Loading Hold
100

50

0

-50

-100
0

50000

100000150000Time , t , [fs]Fig.3.18 Change in total energy under spherical indenter (with crosslink),(v=100m/s).

200000横応力 0 条件下で作成した薄膜への押し込みシミュレーション

Loading

Hold

100

㻱㼠㼛
50

㻱㼎㼑

0

㻱㼎㼟

-50

-100
0

50000

100000150000200000㻱㼢㼣Time , t , [fs]Change in potential energy , E-E0 , [eV]Change in potential energy , E-E0 , [eV]第3章

Loading

33

Hold

100

50

0

-50

-100
0

50000

100000150000200000Time , t , [fs]Loading

Hold

100

㻱㼠㼛
50

㻱㼎㼑

0

㻱㼎㼟

-50

-100
0

50000

100000150000Time , t , [fs]200000㻱㼢㼣Change in potential energy , E-E0 , [eV]Change in potential energy , E-E0 , [eV]Fig.3.19 Change in total energy under spherical indenter (without crosslink),(v=40m/s).

Loading

Hold

100

50

0

-50

-100
0

50000

100000150000Time , t , [fs]Fig.3.20 Change in total energy under spherical indenter (with crosslink),(v=40m/s).

200000

第3章

3.2.4

横応力 0 条件下で作成した薄膜への押し込みシミュレーション

34

局所密度及び２面角の変化

各粒子位置においてカットオフ半径内の粒子数から局所密度 ρを評価した．図3.21及び図3.22は，密度が緩和時より少し高い 0.92<ρ<1.00g/cm3 のノード数を黒色で，密度3が非常に高い ρ>1.00g/cm のノード数を赤色で，それ以外の比較的緩和時と変わらな
3

い ρ<0.92g/cm のノード数を青色で示したものである．なお，2 つのグラフ共に，縦軸は初期緩和計算 100,000fs 時の値からの粒子数の変化量 (N -N０ ) をとって示している．
3

いずれもの系においても，ρ>1.00g/cm の密度が高いノードははじめ減少する．ま3た，ρ<0.92g/cm の低密度ノードも増加している．これは分子鎖が表面流動した効果である．遅い押し込みでは表面流動に十分時間があるため大きく減少する．その後は
3

3

押し込みによって ρ<0.92g/cm の低密度のノードは著しく低下し，ρ>1.00g/cm の高密度が著しく増加する．押し込み保持時は横方向の応力制御が無く上にふたをされた状態で内部緩和が行われるため，低密度への変化が生じている．特に，押し込み速度が遅い図 3.22 では保持時間が長いため，最終的に押し込み前とほぼ同じ状態まで戻っている．
架橋による差を見ると，架橋した (b) は架橋のない (a) に比べて押し込み終了時に低密度 (青色) の減少値が小さく，黒色の 0.92<ρ<1.00g/cm3 の減少が大きい．これは架橋導入により分子鎖同士が束縛されている為，高速での押し込みによる局所的な変形集中が分散されたためと考える．この傾向は，速度が遅い図 3.22 でも同様に見られる．

横応力 0 条件下で作成した薄膜への押し込みシミュレーション

Change in number of node , N-N0

Loading Hold
20000

ρ>1.0
10000

0

0.92<ρ<1.0-100000.92<ρ
0

50000

100000150000200000Change in number of node , N-N0

第3章

35

Loading Hold
20000

ρ>1.0
10000

0

0.92<ρ<1.0
0.92<ρ

-10000050000

100000150000Time , t , [fs]Time , t , [fs](a) Without crosslink(b) Crosslink

200000

Loading

Loading

Hold

20000

ρ>1.0
0.92<ρ<1.0
10000

0

Change in number of node , N-N0

Change in number of node , N-N0

Fig.3.21 Time change in number of node and atomic density under sphericalindenter, (v=100m/s).

Hold

20000

ρ>1.0
0.92<ρ<1.0
10000

0

-10000-100000.92<ρ
0

50000

100000150000200000Time , t , [fs](a) Without crosslink

0.92<ρ
0

50000

100000150000Time , t , [fs](b) Crosslink

Fig.3.22 Change in the local density under spherical indenter, (v=40m/s).

200000

第3章

横応力 0 条件下で作成した薄膜への押し込みシミュレーション36

２面角ノードの変化を，押し込み速度が速い (v=100m/s) ものを図 3.23 に，押し込み速度が遅い (v=40m/s) ものを図 3.24 に示す．前節同様，175◦ <ϕ<180◦ のものを trans，65◦ <ϕ<70◦ のものを gauche，その他の high energy ノードを others として示している．
なお，両図とも初期緩和計算 100,000fs 時の値からの粒子数の変化量 (N -N0 ) として示している．
いずれも，押し込み→圧子保持の切り換えの際に急激な変化は認められず，圧子下での分子鎖の直線化などは生じていない．特に遅い速度での押し込みは，押し込み保持時に圧子に押しのけられた分子鎖が自由表面に流動する為，先の図 3.10(b) 及び図3.15(b) と同様な変化で見えなくなっているものと推測される．

横応力 0 条件下で作成した薄膜への押し込みシミュレーション

Loading Hold500

others0gauche
-500

trans
0

50000

100000150000200000Time , t , [fs]Change in number of node , N-N0

Change in number of node , N-N0

第3章

37

Loading Hold

500

0

-500

0

50000

100000150000200000Time , t , [fs](a) Without crosslink(b) With crosslink

Hold

Loading

Loading

500

others0gauche
-500

trans
0

50000

100000150000Time , t , [fs](a) Without crosslink

200000Change in number of node , N-N0

Change in number of node , N-N0

Fig.3.23 Change in the conformation of torsion under spherical indenter,(v=100m/s).

Hold

500

0

-500

0

50000

100000150000200000Time , t , [fs](b) With crosslink

Fig.3.24 Time change in number of node and torsion angle under sphericalindenter, (v=40m/s).

第3章

横応力 0 条件下で作成した薄膜への押し込みシミュレーション

3.3

ピラミッド型圧子での押し込み

3.3.1

シミュレーション条件

38

前節と同じ，架橋導入・架橋未導入のアモルファスポリエチレンに対し，横方向周期境界条件のもと構造緩和を施し，ピラミッド型圧子を押し込むシミュレーションを行った．図 3.25 にピラミッド圧子の形状及び寸法を示す．圧子の原子数は 23180 であり，前節同様シミュレーション中は C 原子の運動 (=圧子の変形) は考慮しない．図 3.26に押し込み前のスナップショットを示す．押し込み速度，温度等は前節と同じとした．

第3章

横応力 0 条件下で作成した薄膜への押し込みシミュレーション

10n
m

7.0

nm

m
10n

70.5°

Loading Direction

Fig.3.25 Dimensions of diamond pyramid indenter.

Fig.3.26 Snapshots of PE chains in amorphous block and C atoms of diamondpyramid indenter.

39

第3章

横応力 0 条件下で作成した薄膜への押し込みシミュレーション

40

押し込み反力の変化

3.3.2

ピラミッド型圧子の押し込みより得られた押し込み力 R と押し込み時間 t の関係を図3.27 に示す．球体圧子の時と同様，押し込み速度が速い (a) は速度が遅い (b) と比べて大きな反力を示している．最大押し込み時の反力は球体圧子の時と比べ，どのグラフとも約２割程度低い．圧子を保持すると，(a)(b) 共に反力は減少し，いずれも 0.015µN 近傍の値を示す．この値は球体圧子と比べ 0.005µN 程小さい値である．架橋による差は，球体圧子の時に比べ顕著となっており，局所的に押し込まれたときに周囲の分子鎖へのアンカーとしての効果が表れているものと考える．

Loading

Loading

Hold

Without crosslink

0.03

0.02

Crosslink

0.01

0
50000

100000150000200000

Repulsive force , R , [µN]

Repulsive force , R , [µN]

0.04

0

Hold

0.05

0.05

0.04

0.03

0.02

0.01

0
0

50000

100000150000Time , t , [fs]Time , t , [fs](a) v = 100m/s(b) v = 40m/s

Fig.3.27 Relationship between indentation force and time (pyramid indenter).

200000

第3章

3.3.3

横応力 0 条件下で作成した薄膜への押し込みシミュレーション

41

アモルファスの内部エネルギー変化

押し込み時の全ポテンシャルエネルギーの時間変化を図 3.28 及び図 3.29 に示す．
実線は押し込み時，保持時の変化を最小二乗近似したものである．押し込み開始時のエネルギー上昇をみると，球体圧子と異なり，v=100m/s の架橋を含む系のほうが架橋を含まない系よりも上昇幅が大きいことがわかる．一方，v=40m/s の遅い押し込みでは，架橋を含まない系 (a) が含む系 (b) よりもやや大きい．いずれの系でも押し込みを継続していくとエネルギーは減少に転じるが，押し込み速度が遅い系ではこの減少
幅は小さい．押し込み保持時は全ての系でエネルギー上昇が見らえるが，このエネル
ギーの勾配は押し込み時の上昇よりも小さい．

横応力 0 条件下で作成した薄膜への押し込みシミュレーション[x103]Loading Hold
0.002

0.001

0

0

50000

100000150000200000Time , t , [fs]Change in potential energy , E-E0 , [eV]Change in potential energy , E-E0 , [eV]

第3章

42

[x103]Loading Hold
0.002

0.001

0

0

50000

100000150000200000Time , t , [fs](b) With crosslink.

(a) Without crosslink.

[x103]

Loading

Hold

0.002

0.001

0

0

50000

100000150000Time , t , [fs](a) Without crosslink

200000Change in potential energy , E-E0 , [eV]Change in potential energy , E-E0 , [eV]Fig.3.28 Change in total energy under pyramid indenter , (v=100m/s).

Loading

Hold

0.002

0.001

0

0

50000

100000150000Time , t , [fs](b) With crosslink

Fig.3.29 Change in total energy under pyramid indenter , v=40m/s.

200000

第3章

横応力 0 条件下で作成した薄膜への押し込みシミュレーション

43

各ポテンシャルエネルギーの成分ごとの時間変化を図 3.30〜図 3.33 に示す．それぞれ，v=100m/s で架橋なしを図 3.30，架橋ありを図 3.31，v=40m/s で架橋なしを図

3.32，架橋ありを図 3.33 としている．いずれも，初期緩和計算 100,000fs 時の値からの変化量 (E -E0 ) をとって示している．
いずれの系も，全ステップを通して bending と bond stretch のエネルギー変化はほぼ 0 である．ピラミッド型圧子の場合も球体圧子同様，いずれの系でも押し込み開始
時に van der Waals のエネルギー上昇がみられる．架橋がある PE に速い押し込みを行った場合は押し込み終了時に負の値をとっており，押し込み前より高密度の所が存在する事が示唆される．その他の系では，圧子直下の高密度部分の影響よりも，低密度である自由表面に流動した効果が大きいため，圧子保持時には正値に留まっている．
押し込み保持時に van der Waals は再度上昇するが，架橋を含む系はエネルギー上昇の勾配が架橋を含まない系と比べやや小さい．

横応力 0 条件下で作成した薄膜への押し込みシミュレーション

Loading Hold

100

㻱㼠㼛
50

㻱㼎㼑
0

㻱㼎㼟

-50

-100

0

50000

100000150000200000㻱㼢㼣Time , t , [fs]Change in potential energy , E-E0 , [eV]Change in potential energy , E-E0 , [eV]

第3章

44

Loading Hold
100

50

0

-50

-100

0

50000

100000150000200000Time , t , [fs]Loading Hold

100

㻱㼠㼛
50

㻱㼎㼑
0

㻱㼎㼟

-50

-100

0

50000

100000150000Time , t , [fs]200000㻱㼢㼣Change in potential energy , E-E0 , [eV]Change in potential energy , E-E0 , [eV]Fig.3.30 Change in total energy pyramid indenter (without crosslink), (v=100m/s).

Loading Hold

100

50

0

-50

-100

0

50000

100000150000200000Time , t , [fs]Fig.3.31 Change in total energy pyramid indenter (with crosslink), (v=100m/s).

横応力 0 条件下で作成した薄膜への押し込みシミュレーション

Hold

Loading

100

㻱㼠㼛
50

㻱㼎㼑
0

㻱㼎㼟

-50

-100

0

50000

100000150000200000㻱㼢㼣Time , t , [fs]Change in potential energy , E-E0 , [eV]Change in potential energy , E-E0 , [eV]第3章

Loading

45

Hold

100

50

0

-50

-100

0

50000

100000150000200000Time , t , [fs]Loading

Hold

100

㻱㼠㼛
50

㻱㼎㼑
0

㻱㼎㼟

-50

-100

0

50000

100000150000Time , t , [fs]200000㻱㼢㼣Change in potential energy , E-E0 , [eV]Change in potential energy , E-E0 , [eV]Fig.3.32 Change in total energy pyramid indenter (without crosslink), (v=40m/s).

Loading

Hold

100

50

0

-50

-100

0

50000

100000150000200000Time , t , [fs]Fig.3.33 Change in total energy pyramid indenter (with crosslink), (v=40m/s).

第3章

3.3.4

横応力 0 条件下で作成した薄膜への押し込みシミュレーション

46

局所密度変化

図 3.34 及び図 3.35 に局所密度 ρの変化を示す．高，中，低密度の分け方は前章と同じである．
前節で述べた，押し込み初期の高密度ノードの減少と低密度ノードの増加は顕著になっている．ピラミッド型圧子の先端は球体圧子に比べて先端が鋭く，面に沿って分子鎖が z 方向に流動しているためである．特に，押し込み速度が遅い図3.35ではその傾向が顕著である．先の図3.30および図3.31の vanderWaals 変化と同じく，上昇→下降の変化が明白である．これはピラミッド型圧子の投影面積が周期セルを覆う時に対応するものと考える．

横応力 0 条件下で作成した薄膜への押し込みシミュレーション

Loading Hold

Loading Hold
10000

0.92<ρ<1.0
ρ>1.0

5000

0.92<ρ

0-5000-10000050000

100000150000200000Time , t , [fs]Change in number of node , N-N0

Change in number of node , N-N0

第3章

10000

0.92<ρ<1.0
ρ>1.0

5000

0.92<ρ

0-5000-10000050000

100000150000200000Time , t , [fs](a) Without crosslink(b) With crosslink

Loading

Hold

10000

0.92<ρ
5000

0.92<ρ<1.0
0

-5000

ρ>1.0

-10000050000

100000150000200000Change in number of node , N-N0

Change in number of node , N-N0

Fig.3.34 Time change in number of node and atomic density under pyramidindenter , (v=100m/s).

Loading

Hold

10000

0.92<ρ

5000

0.92<ρ<1.0
0

-5000

ρ>1.0
-10000050000

100000150000Time , t , [fs]Time , t , [fs](a) Without crosslink(b) With crosslink

200000

Fig.3.35 Time change in number of node and atomic density under pyramidindenter , (v=40m/s)

47

第4章
横ひずみ 0 条件下で作成した薄膜への押し込みシミュレーション
4.1

横ひずみ 0 条件下での薄膜の作成

4.1.1

シミュレーション条件

先の 14nm × 14nm × 14nm のアモルファス構造について，バルクでの緩和計算後に表面を緩和して作成する際に，横方向応力制御を行わずセル寸法を固定とした．その後，作成した薄膜に対し球体圧子及びピラミッド型圧子の押し込みシミュレーションを行った．寸法，押し込み深さは前節同様であるが，押し込み速度は 500m/s とした．
押し込みプロファイルを図 4.1 に示す．制御温度は構造緩和時と同様 300K である．

Indentation depth , d , [nm]Loading

Hold

4

3

2

1

0
0

4000

8000

12000

Time , t , [fs]Fig.4.1 Load waveform.

48

16000

第4章

4.1.2

横ひずみ 0 条件下で作成した薄膜への押し込みシミュレーション

49

初期構造緩和結果及び考察

前章同様，50,000fs の初期緩和中のエネルギー変化を，架橋なしを図 4.2(a)，架橋ありを図 4.3(a) に示す．また，1,000fs 以降のエネルギー変化を拡大した架橋なしを図

4.2(b)，架橋ありを図 4.3(b) に示す．比較のため，全てのグラフの時間スケールを前章で扱った 100,000fs に揃えている．1,000fs 以降の拡大図 4.2(b) でもやはりエネルギーは減少せず，緩やかに上昇している．一方，架橋ありの図 4.3(b) はエネルギーが大きく振動した後一旦減少し，横応力 0 条件下と同様のエネルギー概形を示した．ポテンシャ
ルエネルギーの各成分を，1,000fs における値からのエネルギー変化量を縦軸にとり架橋なしを図 4.4 に，架橋ありを図 4.5 に示す．架橋ありの図 4.5 では，van der Waals は
架橋なしの時と異なり緩和開始直度に急減し 20,000fs 近傍で底打ちする前章の変化と同様の傾向を示す一方，架橋なしの図 4.4 では 10,000fs 近傍で底打ちし，かつ減少も小さい．薄膜全体の密度変化を，架橋なしを図 4.6(a)，架橋ありを図 4.7(a) に，厚さ変化を図 4.6(b)，架橋ありを図 4.7(b) に示す．なお，横方向のセル長を固定している為，底面積は初期値 (14nm × 14nm) で不変としている．密度変化をみると，いずれも前章の様な緩和直後からの密度増加は見られず終始単調に減少し，薄膜の厚さは z 方向へ増加している事が分かる．つまり，緩和するに従い密度が低下し厚さが増加している事から，横方向応力制御を外しても自由表面側に分子鎖が流動していることが示唆される．各粒子位置においてカットオフ半径内の粒子数から評価した局所密度 ρ のノー
ド変化を，架橋なしを図 4.8(a)，架橋ありを図 4.9(a) に示す．高，中，低密度の分け方
は前章と同じである．いずれも緩和全体を通して赤色のグラフが減少し青色のグラフは増加していることから，全体的に密度が減少していることが見て取れる．図 4.9(a)では架橋なしの系と比べやや大きく減少するが，薄膜の厚さは架橋なしの系よりも増えていない．架橋の効果により粒子が自由表面へ流動しにくくなった為，薄膜内部でより緩和が進んでいる．全ての 2 面角ノード変化を架橋なしを図 4.8(b)，架橋ありを
図 4.9(b) に示す．trans，gauche，others の分け方は前章と同様である．大きな変化は認められないが，前章とは異なり最安定点の trans が増加し，準安定点の gauche はほぼ変わらず，エネルギー的に不安定な others が減少している事から，2 面角は安定の方向へ向かっている．これと図 4.4 の torsion 変化より，自由表面へ流動している粒子

第4章

横ひずみ 0 条件下で作成した薄膜への押し込みシミュレーション

50は全体の個数からは少数であるものの，１つ１つが高エネルギーとなったため全体のtorsion が増加したが，内部は安定となる 2 面角が形成されていることが分かる．架橋ありを図 4.9(b) は架橋による反発で緩和ごく初期は架橋なしの系と異なるが，緩和後期はほぼ同様の傾向を示す．

第4章

横ひずみ 0 条件下で作成した薄膜への押し込みシミュレーション[x103]
Potential energy , E , [eV]

Potential energy , E , [eV][x103]

51

-0.04

-0.05

-0.06

-0.07

-0.08
0

-0.078-0.0785-0.079-0.0795-0.08

20000 40000 60000 80000 100000

0Time , t , [fs]20000 40000 60000 80000 100000

Time , t , [fs]（a) Throughout the simulation.

（b) Magnification.

Fig.4.2 Relationship between potential energy and time (without crosslink).

[x103]
Potential energy , E , [eV]

Potential energy , E , [eV][x103]
-0.04

-0.05

-0.06

-0.07

-0.08
0

20000 40000 60000 80000 100000-0.078-0.0785-0.079-0.0795-0.08

020000 40000 60000 80000 100000

Time , t , [fs]Time , t , [fs]（a) Throughout the simulation.

（b) Magnification.

Fig.4.3 Relationship between potential energy and time (with crosslink).

横ひずみ 0 条件下で作成した薄膜への押し込みシミュレーション

㻱㼠㼛

100

㻱㼎㼑
0

㻱㼎㼟
-100
0

20000 40000 60000 80000 100000

㻱㼢㼣Time , t , [fs]Change in potential energy , E-E0 , [eV]Change in potential energy , E-E0 , [eV]

第4章

52

100

0

-100
0

20000 40000 60000 80000 100000

Time , t , [fs]

㻱㼠㼛

100

㻱㼎㼑
0

㻱㼎㼟
-100
0

20000 40000 60000 80000 100000

Time , t , [fs]㻱㼢㼣Change in potential energy , E-E0 , [eV]Change in potential energy , E-E0 , [eV]Fig.4.4 Change in components of potential energy (without crosslink).

100

0

-100
0

20000 40000 60000 80000 100000

Time , t , [fs]Fig.4.5 Change in component of potential energy (with crosslink).

第4章

横ひずみ 0 条件下で作成した薄膜への押し込みシミュレーション

53

15.2

Thickness , [nm]Density , [g/cm3]

1.04

1.02

1

0.98

15
14.8
14.6
14.4
14.2

0.96

14
0

0

20000 40000 60000 80000 100000

20000 40000 60000 80000 100000

Time , t , [fs]Time , t , [fs](a) Change in the density.

(b) Thickness at z direction.

Fig.4.6 Change in the density and thickness (without crosslink).

15.2

Thickness , [nm]Density , [g/cm3]

1.04

1.02

1

0.98

15
14.8
14.6
14.4
14.2

0.96

14
0

20000 40000 60000 80000 100000

Time , t , [fs](a) Change in the density.

020000 40000 60000 80000 100000

Time , t , [fs](b) Thickness at z direction.

Fig.4.7 Change in the density and thickness (with crosslink).

横ひずみ 0 条件下で作成した薄膜への押し込みシミュレーション

ρ<0.92

10000

0.92<ρ<1.0
0-100001.0>ρ

Change in number of node , N-N0

Change in number of node , N-N0

第4章

54

㟷 㻩gauche㉥ 㻩trans㯮 㻩others

10000

0

-10000-20000020000 40000 60000 80000 100000

020000 40000 60000 80000 100000

Time , t , [fs]Time , t , [fs](a) Number of density nodes.

(b) Number of tosional nodes.

0.92<ρ<1.0
ρ<0.9210000

0-10000Change in number of node , N-N0

Change in number of node , N-N0

Fig.4.8 Change in local density and tosional nodes (without crosslink).

㟷 㻩gauche㉥ 㻩trans㯮 㻩others

10000

0

-10000-200001.0>ρ
020000 40000 60000 80000 100000

020000 40000 60000 80000 100000

Time , t , [fs]Time , t , [fs](a) Number of density nodes.

(b) Number of tosional nodes.

Fig.4.9 Change in density and tosional nodes (with crosslink).

第4章

横ひずみ 0 条件下で作成した薄膜への押し込みシミュレーション

4.2

球体及びピラミッド型圧子の押し込み

4.2.1

押し込み反力の変化

55

押し込みシミュレーションにより得られた押し込み力 R と押し込み時間 t の関係を，
球体圧子を図 4.10(a) に，ピラミッド型圧子を図 4.10(b) に示す．押し込み力 R は圧子の C 原子が受ける力の z 方向成分により評価した．押し込み速度が 500m/s と速いた
め，反力はいずれも前章よりも高い値を取っている．圧子の押し込みに伴い反力は増加するが，球体圧子の方が大きな反力を示している．圧子を保持すると，(a)(b) 共に
反力は減少するが，球体圧子の方は同じ押し込み深さでの体積が大きいためピラミッ
ド型圧子よりも高い反力を示す．
架橋による差は，球体圧子ではほとんどないように見えるが，ピラミッド型圧子では押し込み時の反力の勾配に大きな差があり，最大値も著しく異なる．圧子保持時も架橋を含む系が含まない系よりも大きな反力を示す．やはり，架橋点の効果は局所的に押し込まれた時の周囲の分子鎖へのアンカー効果で顕著となる．

Hold

Loading

0.2

Without crosslink

0.1

Crosslink
0
0

5000

10000

Time , t , [fs]（a) Spherical indenter.

15000

Repulsive force , R , [µN]

Repulsive force , R , [µN]

Loading

Hold

0.2

0.1

0
0

5000

10000

Time , t , [fs]（b) Pyramid indenter.

Fig.4.10 Relationship between indentation force and time.

15000

第4章

4.2.2

横ひずみ 0 条件下で作成した薄膜への押し込みシミュレーション

56

内部エネルギー変化

押し込み時の全ポテンシャルエネルギーの時間変化を図 4.11 および図 4.12 に示す．
これまで同様，架橋未導入のグラフが (a)，架橋導入のグラフが (b) としている．
いずれも圧子が表面に触れた時にエネルギーが大きく振動し，低いエネルギー状態となる．その後エネルギーは上昇するが，架橋ありの系に球体圧子で押し込んだ図 4.11(b)以外では一定値で頭打ちとなり横ばいとなっている．押し込み保持時には球体圧子は初期にエネルギーの減少がみられたが，ピラミッド型圧子では架橋を含む (b) のみわ
ずかに認められるもののほとんど低下していない．

[x103]横ひずみ 0 条件下で作成した薄膜への押し込みシミュレーション
Hold

Loading0.001

0-0.001

0

5000

10000

15000

Time , t , [fs]Change in potential energy , E-E0 , [eV]Change in potential energy , E-E0 , [eV]

第4章

[x103]

57

Hold

Loading0.001

0-0.001

0

4000

8000

12000

16000

Time , t , [fs]（a) Without crosslink.

（b) With crosslink.

[x103]

Loading

Hold

0.001

0

Change in potential energy , E-E0 , [eV]Change in potential energy , E-E0 , [eV]Fig.4.11 Change in total energy under spherical indenter.

[x103]

Loading

Hold

0.001

0

-0.001-0.001

0

4000

8000

12000

Time , t , [fs]（a) Without crosslink.

16000

0

4000

8000

12000

Time , t , [fs]（b) With crosslink.

Fig.4.12 Change in total energy under pyramid indenter.

16000

第4章

4.2.3

横ひずみ 0 条件下で作成した薄膜への押し込みシミュレーション

58

ポテンシャル成分の変化

各ポテンシャルエネルギー成分の時間変化を図 4.13〜図 4.16 に示す．
いずれの系も，表面に接触した際に bond stretch が急上昇し，torsion が急減している．また，van der Waals も減少している．これは，表面の分子鎖が圧子表面に張り付いた形状となり直線化したためと考える．先の横応力 0 に制御したアモルファスでこの応答がみられなかったのは，分子鎖流動によって自由に動ける分子鎖が表面近傍ではほぼ直線状態になっており，絡んでいない状態になっていた可能性がある．
球体圧子では，押し込み後期に van der Waals の増加がみられ，これが圧子を停止させた時に減少する事でエネルギーの低下をもたらしている．これは先に述べた「粘弾性挙動」に対応するものと考える．ピラミッド型圧子でも同様の傾向がみられるが，同一の押し込み深さでは押し込み体積が小さいため変化が大きく現れなかったものと考える．

横ひずみ 0 条件下で作成した薄膜への押し込みシミュレーション

Loading

Hold

400

㻱㼠㼛

200

㻱㼎㼑

0

㻱㼎㼟

-200

0

4000

8000

12000

16000

㻱㼢㼣

Time , t , [fs]Change in potential energy , E-E0 , [eV]Change in potential energy , E-E0 , [eV]

第4章

59

Loading

Hold

400

200

0

-200

0

4000

8000

12000

16000

Time , t , [fs]Loading

Hold

400

㻱㼠㼛

200

㻱㼎㼑

0

㻱㼎㼟

-200

0

4000

8000

12000

Time , t , [fs]

16000

㻱㼢㼣

Change in potential energy , E-E0 , [eV]Change in potential energy , E-E0 , [eV]Fig.4.13 Change in total energy under sphere indenter , (without crosslink).

Hold

Loading
400

200

0

-200

0

4000

8000

12000

Time , t , [fs]Fig.4.14 Change in total energy under sphere indenter , (with crosslink).

16000

横ひずみ 0 条件下で作成した薄膜への押し込みシミュレーション

Loading

Hold

400

㻱㼠㼛

200

0

㻱㼎㼑

-200

㻱㼎㼟
0

4000

8000

12000

16000

㻱㼢㼣

Time , t , [fs]Change in potential energy , E-E0 , [eV]Change in potential energy , E-E0 , [eV]第4章

Loading

60

Hold

400

200

0

-200

0

4000

8000

12000

16000

Time , t , [fs]Loading

Hold

400

200

㻱㼠㼛

0

㻱㼎㼑
㻱㼎㼟

-200

0

4000

8000

12000

Time , t , [fs]

16000

㻱㼢㼣

Change in potential energy , E-E0 , [eV]Change in potential energy , E-E0 , [eV]Fig.4.15 Change in total energy under pyramid indenter , (without crosslink).

Loading

400

Hold

200

0

-200

0

4000

8000

12000

Time , t , [fs]Fig.4.16 Change in total energy under pyramid indenter , (with crosslink).

16000

第4章

4.2.4

横ひずみ 0 条件下で作成した薄膜への押し込みシミュレーション

61

局所密度及び２面角の変化

図 4.17 および図 4.18 に局所密度 ρの変化を示す．
前章の横応力 0 の薄膜に比べ，押し込み初期の密度変化は小さい．架橋した系では高密度ノードがやや減少しているのが認められる．この初期応答の違いは，やはり前章のアモルファス表面は絡みが解消した形状の抵抗がない構造になっているため，圧子下とそれ以外の密度変化が生じやすかったものと考える．押し込み後期には圧子下での局所密度の増加が認められる．やはり押し込み体積が多い球体圧子の方がこの変化が著しい．架橋による差は押し込み体積が小さいピラミッド型圧子ではほとんどな3

いが，球体圧子では差が生じている．
ρ>1.00g/cm の圧子直下の高密度部分に差はない
3が，架橋があると ρ<0.92g/cm の「疎」になった部分が抑えられている．架橋がある
場合はアンカー効果によって全体的に圧縮されるので黒の 0.92<ρ<1.00g/cm3 も大きく減少している．保持時の流動変形で圧子下の高密度部分が緩和されるが，これは主に圧子下の高密度ノード・圧子下でない表面の低密度ノード以外となる為，黒線が架橋の有無にかかわらず同じように増加する．

横ひずみ 0 条件下で作成した薄膜への押し込みシミュレーション
Hold

Loading

10000

Loading

ρ>1.00.92<ρ<1.0

0

0.92<ρ

-10000

0

4000

8000

12000

16000

Time , t , [fs]Change in number of node , N-N0

Change in number of node , N-N0

第4章

62

Hold

ρ>1.0

10000

0.92<ρ<1.0
0-100000.92<ρ
0

4000

8000

12000

16000

Time , t , [fs](a) Without crosslink(b) Crosslink

Fig.4.17 Change in the atomic density (spherical indenter).

Hold

ρ>1.0
0.92<ρ<1.0

10000

0

-100000.92<ρ
0

4000

8000

12000

Time , t , [fs](a) Without crosslink

Hold

Loading

16000

Change in number of node , N-N0

Change in number of node , N-N0

Loading

ρ>1.00.92<ρ<1.0

10000

0

-100000.92<ρ
0

4000

8000

12000

Time , t , [fs](b) Crosslink

Fig.4.18 Change in the atomic density (pyramid indenter).

16000

第4章

横ひずみ 0 条件下で作成した薄膜への押し込みシミュレーション

63

２面角ノードの変化を，図 4.19 および図 4.20 に示す．
先述の圧子表面への分子鎖はりつきが，trans ノードの著しい増加と others の減少から確認される．一方，押し込み後期や保持への切り換え時に特に変化は見られず，これまでの応答は分子鎖の形態変化よりも van der Waals による時間遅れが支配的であると結論付けられる．
Hold

Loading

1000

others0gauche
-1000

trans
-2000
0

4000

8000

12000

16000

Time , t , [fs]Change in number of node , N-N0

Change in number of node , N-N0

Loading

Hold

1000

0

-1000

-2000
0

4000

8000

12000

16000

Time , t , [fs](a) Without crosslink(b) With crosslink

Loading

Hold

Loading

1000

others0gauche
-1000

trans
-2000
0

4000

8000

12000

Time , t , [fs](a) Without crosslink

16000

Change in number of node , N-N0

Change in number of node , N-N0

Fig.4.19 Change in the conformation of torsional angle (spherical indenter).

Hold

1000

0

-1000

-2000
0

4000

8000

12000

16000

Time , t , [fs](b) With crosslink

Fig.4.20 Change in the conformation of torsional angle (pyramid indenter).

第5章
結論
本研究では，高分子の分子鎖の変形メカニズムを原子レベルから検討する為，アモ
ルファスポリエチレン (PE) に対する押し込みシミュレーションを分子動力学法を用いて行った．周期境界下でバルクとして作成したアモルファス構造に，一方向の周期境界を外して表面とした無限薄膜を横方向応力制御 0 としたものと横ひずみ 0 で緩和
し，それぞれに架橋ポテンシャルを導入して架橋の有無を考慮しながら球体およびピ
ラミッド型圧子の押し込みシミュレーションを行った．以下に，得られた結果を総括する．
第 2 章では，本研究で用いた解析手法の基礎について述べた．まず，分子動力学法の概要ならびに基礎方程式を示し，本研究で用いた数値積分法について説明した．次に，粒子間相互作用の評価に用いられるポテンシャルエネルギーについて述べ，CH2をひとつの粒子として扱う United Atom model によるポリエチレンのポテンシャル関
数と，化学的架橋をモデル化した架橋ポテンシャルについて詳述した．さらに，大規模シミュレーションに必要な計算の高速化手法について述べた．
第 3 章では，横応力 0 の条件下で無限薄膜 PE を作成した後，球体およびピラミッド型圧子を押し込むシミュレーションを行った．まず横方向 0 下での構造緩和について，
内部エネルギー，各ポテンシャルエネルギーの成分，局所密度，2 面角の変化から，架橋の有無にかかわらず (1) 厚さ方向は数 nm で横方向無限のアモルファス PE に，横方
向応力 0 の制御を行うと厚さ方向の変形抵抗が小さいため急激に圧縮される，(2) そのため全体の密度が大きく上昇し自由表面側に分子鎖が流動する，(3) その後も横方向の圧縮が続けられ周期セルが (厚さ方向) に細長い形状になる，などの問題点が明らかに

64

第5章

結論

65

なった．100,000fs 緩和時の状態で応力制御を外し，押し込みシミュレーションを行った結果，(4) 架橋を有する系は反力が上昇し，特にピラミッド型圧子の押し込みで顕著となる，(5) 押し込みを行うと内部エネルギーが上昇→減少する，(6)(4) の変化は vander Waals によるものであり，圧子の下で押し込みにより高密度となったため，である．(7) 圧子保持時は van der Waals が上昇することによってポテンシャルエネルギーが上昇する，(8) 一見奇異に見える (7) の挙動は，高密度になると近接粒子同士の vander Waals エネルギーが上がってもそれより遠方の粒子が引力に寄与するためエネル
ギーが下がるという (5) の挙動で説明できる．これが高分子の「粘弾性挙動」に対応する，などを明らかにした．
第 4 章では，横ひずみ 0 条件下で無限薄膜 PE を作成し押し込みシミュレーションを行った．まず構造緩和時の変化から，(1) 横ひずみ 0 条件下でも分子鎖は自由表面に流動し密度が低下する，(2) 架橋を有すると (1) の効果は抑制される，ことがわかった．
球体圧子・ピラミッド圧子による押し込みを行った結果，(3) 押し込み深さ-反力は架橋した系が反力が大きくなり，その効果はやはりピラミッド型が顕著である，(4) 横方
向応力 0 のときと異なり圧子が接触した時にエネルギーが大きく減少する，(5) 接触時
は bond stretch が上昇，torsion が減少しており，圧子表面に付着した分子鎖が直線化する為 (4) を生じる，(6) 球体圧子では押し込み後期に van der Waals が上昇する，(7)球体圧子のみ，押し込み保持時にエネルギーが減少したがこれは (6) の高密度部分が緩和されたためである，などが明らかとなった．
アモルファス構造では最安定という状態は存在せず，特に今回は有限セルでの表面構造作成時の状態に依存した結果となった．すなわち，横方向応力 0 とした系では絡み点が解消した状態の側面となり，そのため今回の検証だけでは普遍的なメカニズムを見出す事は難しいが，高密度となっても van der Waals が下がり，それが時間遅れ応答を担うという点では今回の検証で得た重要な知見である．

参考文献[1] 三宅正二郎, ”Evaluation of Nanomechanical Properties of Surfaces.”, 精密工学
会誌, 73, 8, (2007), 859-863.

[2] Jee, Ah-Young, and Minyung Lee, ”Comparative analysis on the nanoindentationof polymers using atomic force microscopy.”, Polymer Testing, 29, 1, (2010), 9599.
[3] Ho, S. P., et al, ”Nanoindentation properties of compression-moulded ultra-highmolecular weight polyethylene.”, Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine, 217, 5, (2003), 357-366.
[4] 中上明光, and 川上信之, ”ナノインデンテーション法による薄膜の機械的特性評価”, RD 神戸製鋼技報, 52, 2, (2002), 74-77.

[5] Loubet, J. L., W. C. Oliver, and B. N. Lucas, ”Measurement of the loss tangent of low-density polyethylene with a nanoindentation technique.”, Journal of
Materials Research, 15, 05, (2000), 1195-1198.
[6] 福田 光完, 玉井 良則, 桑島 聖, ”分子動力学シミュレーションによるアルキル基側鎖を有する液晶ポリエステル層状構造の作成” ,Journal of Computer Chemistry,

Japan, Vol.3, No.1, (2004), pp.13-20.
[7] Yamamoto, Takashi, ”Molecular dynamics simulation of polymer crystallization
through chain folding.”, The Journal of chemical physics, 107, 7, (1997), 26532663.
[8] Ko, Min Jae, et al., ”Characterization of polyethylene crystallization from an oriented melt by molecular dynamics simulation.”, The Journal of chemical physics,
121, 6, (2004), 2823-2832.

66

参考文献

67

[9] Mansﬁeld, Kevin F., and Doros N. Theodorou., ”Molecular dynamics simulationof a glassy polymer surface.”, Macromolecules, 24, 23, (1991), 6283-6294.
[10] Thompson, Peter A., and Sandra M. Troian., ”A general boundary condition for
liquid ﬂow at solid surfaces.”, Nature , 389, 6649, (1997), 360-362.
[11] Fan, Cun Feng, and Tahir Cag̈in., ”Wetting of crystalline polymer surfaces:
A molecular dynamics simulation.”, The Journal of chemical physics, 103, 20,(1995), 9053-9061.
[12] Bulacu, Monica Iulia., ”Molecular dynamics studies of entangled polymer chains”,
University Library of Groningen, Host, 2008.
[13] Takeuchi, Hisao, and Ryong‐Joon Roe., ”Molecular dynamics simulation of local
chain motion in bulk amorphous polymers. I. Dynamics above the glass transition.”, The Journal of chemical physics, 94, 11, (1991), 7446-7457.
[14] Takeuchi, Hisao, and Ryong‐Joon Roe., ”Molecular dynamics simulation of local
chain motion in bulk amorphous polymers. II. Dynamics at glass transition.”, The
Journal of chemical physics, 94, 11, (1991), 7458-7465.
[15] Moller, M. A., et al., ”Molecular dynamics simulation of a Langmuir?Blodgett
ﬁlm.”, The Journal of chemical physics, 94, 12, (1991), 8390-8401.
[16] Harmandaris, Vagelis A., Kostas Ch Daoulas, and Vlasis G. Mavrantzas., ”Molecular dynamics simulation of a polymer melt/solid interface : local dynamics andchain mobility in a thin ﬁlm of polyethylene melt adsorbed on graphite.”, Macromolecules, 38, 13, (2005), 5796-5809.
[17] Bandyopadhyay, Ananyo, et al., ”Molecular modeling of crosslinked epoxy polymers: The eﬀect of crosslink density on thermomechanical properties.”, Polymer,
52, 11, (2011), 2445-2452.

参考文献

68

[18] Li, Ying, Martin Kroger, and Wing Kam Liu., ”Primitive chain network study
on uncrosslinked and crosslinked cis-polyisoprene polymers.”, Polymer, 52, 25,(2011), 5867-5878.
[19] Öttinger, Hans Christian., ”Coarse-graining of wormlike polymer chains for
substantiating reptation.”, Journal of non-newtonian fluid mechanics, 120, 1,(2004), 207-213.
[20] 増渕 雄一, ”高分子系のナノシミュレーション”, CICSJ-Bulletin, 24, (2006), 152.
[21] 屋代 如月, 冨田 佳宏, ”高分子材料の力学特性評価と高機能化 (20) -粗視化分子動力学法による高分子ミクロシミュレーション その 2-”, 機械の研究 Vol.63, No.11,(2011), pp.955-961.
[22] 谷口 誠典, and 屋代如月. ”820 架橋条件の違いによる鎖状高分子のヒステリシス
変化: 分子動力学法による検討 (OS7. 材料の組織・強度に関するマルチスケール
アナリシス/OS8. 電子・原子・マルチシミュレーションに基づく材料特性評価 合
同ポスターセッション, OS・一般セッション講演).” 計算力学講演会講演論文集

2011.24 (2011): 269-270.
[23] 屋代 如月, et al. ”A molecular dynamics study on hysteresis of amorphous polymers.” 日本機械学会論文集 A 編 72.715 (2006):[24] 谷口 誠典, 神戸大学大学院工学研究科修士論文, (2012).
[25] 屋代, 古田, 冨田, ”分子動力学法による結晶/非晶ポリエチレンへのナノインデン
テーション”, 計算力学講演会講演論文集, (2005), pp.673-674.

[26] Rigby, David, and Ryong‐Joon Roe. ”Molecular dynamics simulation of polymer liquid and glass. I. Glass transition.” The Journal of chemical physics 87.12(1987): 7285-7292.
[27] 上田 顯, コンピュータシミュレーション, (1990), 朝倉書店.

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

 粗視化分子動力学によるポリエチレンのエントロピー挙動の検討
藤麻

成貴，屋代

如月

日本材料学会 第 18 回分子動力学シンポジウム, 東京工業大学, (2013.5)69学術論文・学術講演
論文貼り付け用ページ70学術論文・学術講演
論文貼り付け用ページ

71学術論文・学術講演
論文貼り付け用ページ72

謝辞
本研究を遂行するにあたり，日ごろの研究生活に対して幅広い視点から御助言・御指導を賜りました屋代如月教授 (現岐阜大学) に，深甚な感謝の意を表します．研究を進めるにあたり，研究報告会だけでなく日頃から多くの御助言を賜りました田中克志
教授に心より感謝いたします．また，日々の研究活動が円滑に行えるよう便宜を図っていただきました片山雷太技術職員に心から御礼申し上げます．
日々の研究内容や大学院一年次の就職活動において，若輩だった私に多くの有益なアドバイスを与えてくださった表面・界面工学研究室の諸先輩方に厚く御礼申し上げます．特に浅識な著者にプログラム，材料等に関して幅広く教えてくださった藤原正大様 (現 株式会社神戸製鋼所)，山根尭様 (現 コマツ)，寺田駿様 (現 マツダ株式会社) に深く感謝致します．
さらに，常日頃から共に研究に取り組み，御助言，御協力，御討論下さいました同研
究室諸氏に深く感謝致します．同じ研究室に配属されて以来共に研究に取り組み，互いに切磋琢磨し，時には励ましあった伊藤貴彦君，小出卓矢君，菱谷卓也君，森太一
君に御礼申し上げます．特に，同じシミュレーション班で日々議論しあい多くの時間を共有した津川悠太君，毛利友宙君，本馬剛徳君には研究活動においても多大なご協力を頂きました．御礼申し上げます．
最後に，六年間の学生生活を暖かく見守り，経済的にも精神的にも支えて頂いた父
母に最大限の敬意と感謝の意を表します．
平成 27 年 2 月
藤麻73

成貴