卒 業 論 文(題 目)

分子動力学法によるアモルファスポリ
エチレン/グラファイト界面のはく離シミュレーション

平成 27 年度

岐阜大学工学部
機械システム工学科

氏名 堀 広志目次
1

緒言 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

解析手法の基礎 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

2.1

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

2

2.2

原子間ポテンシャル

3

2.3

高速化手法

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

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

3

シミュレーション条件

4

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

15
18

4.1

はく離強度

4.2

応力成分と界面の詳細 . . . . . . . . . . . . . . . . . . . . . . . . . . 29

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

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

39

参考文献 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

41

謝辞 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42

5

i

1

緒言
CFRP(炭素繊維強化プラスチック) は，鋼やアルミニウムと比べ比強度，比弾性率で優れており，航空機や自動車など様々な分野に応用されている．

CFRP の破壊で重要となるのが繊維-樹脂の界面強度である．実験では，高速で生じる界面の破壊現象を高速度カメラによる観察などが行われている (1) ．また，界面の強度を実験的に評価するものとして繊維引き抜き試験 (2) 等による評価がなされているが，様々な樹脂・繊維界面の組み合わせを考えると，高分子-炭素繊維表面のミクロ挙動に関する基礎的な知見が必要と考える．原子シミュレーションで高分子-炭素界面を検討したものとして，Frankland らはカーボンナノチューブ (CNT) をポリエチレン (PE) 分子鎖中に含有させた分子動力学 (MD) シミュレーションによりその応力-ひずみ応答を評価している (3) ．Zhu らは CNT とエポキシ樹脂の混合物の強度を検証している (4) ．Arash らも CNT と高分子の混合物について引っ張り時の各種
パラメータを 検討している (5) ．Li らは高分子と CNT の引き抜きシミュレーションを行い，引き抜き力は CNT の直径に比例することなどを示している (6) ．
上記のように，フラーレン・ポリマー界面の挙動は盛んに検討されているものの，基本は van der Waals(VDW) 力による結合のため, フラーレンのようにカットオフ距
離レベルの構造体では大きな差はないと考える．CFRP 中のポリマー・炭素界面の知見にはむしろメゾスケールの表面凹凸が重要である．本研究では，グラファイト
層間にアモルファス PE 分子鎖を充填し，圧着したあとで剥離する MD シミュレー
ションをグラファイト壁の表面凹凸に着目して行う．

1

2

解析手法の基礎2.1

分子動力学法

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

mi

d2 ri
= Fi
dt2

(2.1)

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

Fi = −

∂Etot
∂ri

(2.2)

式 (2.2) の数値積分には，Verlet の方法，予測子-修正子法等がよく用いられる．本
研究では，以下に示す Verlet の方法を用いた．
時刻 t + ∆t と t − ∆t での粒子 i の位置ベクトル ri (t ± ∆t) を taylor 展開すると，

ri (t + ∆t) = ri (t) + ∆t

dri (t) (∆t)2 d2 ri (t) (∆t)3 d3 ri (t)
+
+
+・
・
・
dt
2
3!
dt3
dt 2

(2.3)

ri (t − ∆t) = ri (t) − ∆tdri (t) (∆t)2 d2 ri (t) (∆t)3 d3 ri (t)
+
+・
・
・
−
dt
2
3!
dt3
dt 2

(2.4)

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

dri= vi (t)dt(2.5)

であり，式 (2.1) と式 (2.5) を式 (2.3) と式 (2.4) に代入すると，

ri (t + ∆t) = ri (t) + ∆tvi (t) +(∆t)2 Fi (t) (∆t)3 d3 ri (t)
+
+・
・
・
2
mi
3!
dt3
2

(2.6)

ri (t − ∆t) = ri (t) − ∆tvi (t) +(∆t)2 Fi (t) (∆t)3 d3 ri (t)
−
+・
・
・
2
mi
3!
dt3

(2.7)

Fi (t)
+・
・
・
mi

(2.8)

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

ri (t + ∆t) + ri (t − ∆t) = 2ri (t) + (∆t)2

ri (t + ∆t) − ri (t − ∆t) = 2∆tvi (t) +・
・
・

(2.9)

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

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

2.2

原子間ポテンシャル

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

Etot は次式のように表される．
Etot = ΦBS (r) + ΦBE (θ) + ΦTO (ϕ) + ΦVW (r̄)
3

(2.12)

右辺各項は，分子内の炭素間の結合長 r，結合角 θ，2 面角 ϕ の共有結合に対するポテンシャルと，異分子鎖間および同一分子鎖内の十分離れた非共有原子間距離 r̄ に対する van der Waals ポテンシャルを表す．本研究で解析対象とするポリエチレン
に対しては，水素原子の運動は陽に扱わず，メチレン炭素 (−CH2 −) を単一粒子として粗視化する united atom model により，以下の関数系およびパラメータが提案されている (7) ．

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

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

torson ポテンシャル (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 にそれぞれ示した．van der Waals ポテンシャルの計算は異なる分子鎖の粒子間，4 原子団以上離れた粒子間，そしてグラファイト壁の C 原子と -CH2 - に対して行う．
式 (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 に示した．

4

⃝ 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.

5⃝ bending ポテンシャル

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

Table 2.2 : Potential parameter for bending.

C – CH2 – C

θ0[deg.]
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 θ.

6⃝ torsion ポテンシャル

ΦTO (ϕ) =∑(V1 cos ϕ + V2 cos 2ϕ + V3 cos 3ϕ + V6 cos 6ϕ)
Table 2.3 : Potential parameter for torsion.

V1[kJ/mol]C – CH2 – CH2 – C
3.935

V2
V3
V6
[kJ/mol] [kJ/mol] [kJ/mol]
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 ϕ.

7⃝ 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̄.

8

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

Φ′ (r) αβrr(2.17)

F β = −F α(2.18)

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

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

4 体間相互作用
α

A

=

Aβ =
cos ϕ =Fα =Fβ =rαγ · rγη γη
rr −|rγη |2rβη · rγη γη
r βγ −
r|rγη |2
Aα · Aβ|Aα ||AβAα · Aβ αΦ′ (ϕ)1β{A−A }sinϕ |Aα ||Aβ ||Aα |2αγ

Φ′ (ϕ)
Aα · Aβ β
1
α
A }
{A
−
sinϕ |Aα ||r β ||Aβ |2r αγ · r γη α r βη · r γη βF −F|r γη |2|r γη |2r αγ · r γη α r βη · r γη β
= −F β +F +F|r γη |2|r γη |2

(2.22)(2.23)

F γ = −F α −(2.24)

Fη

(2.25)

9

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

Aj

rkl
l

Ai riki

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

10

2.3

高速化手法

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

11

rc

rfc

Fig. 2.8 : Schematic of bookkeeping method.

12

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

13

y
by
x

0
bx

(a) Domain decomposition(b) Block position

Fig. 2.9 : Schematic of domain decomposition method.

14

3

シミュレーション条件
x ,y 方向は周期境界，z 方向は壁面の条件下で 15[nm] × 15[nm] × 15[nm] の立方体セル中に 1000CH2 の PE 分子鎖を 100 本ランダムに成長させて PE 無限平板を作成する．図 3.1 に模式的に示すように，作成した PE 板の上下にグラファイト層を配置し，600[K] で 5000[fs] の緩和計算を行ったあと，25000[fs] でグラファイト壁を壁の厚さ hwall 分接近させて (上を −hwall /2，下を hwall /2 移動)，PE 分子鎖とグラファ
イト壁を圧着させる．その後，10000[fs] かけて 600[K] から 300[K] まで冷却して初期配置とした．
グラファイト壁は，模式図 3.1(a) のようにグラファイト層をセル底面に平行に配置したものについて，(1) 表面凹凸なし，(2)x 方向に正弦波で凹凸をつけたもの (1次元パターン)，(3)x ，y 方向に正弦波の凹凸をつけたもの (2 次元パターン)，について検証した．図 3.2〜3.4 にそれぞれ例を示す．正弦波の周期は 2π, 6π, 10π とした．
次に模式図 3.1(b) のようにグラファイトを垂直に配置したものについて同様に (1)，(2)，(3) のケースについて検討した．
圧着後，150 万ステップで 300% の引張となる割合で一定速度でグラファイト壁を上下に移動させ，剥離シミュレーションを行った．準静的な条件とするためグラ
ファイト壁だけでなく，PE 分子鎖にも位置に応じたアフィン変形を与えている．

15

PE

PE

( a ) Ho rizo n t a l( b ) Ve rt ica lFig. 3.1 : Schematic of PE block and graphite walls.

Fig. 3.2 : PE chains and graphite walls(flat surface).

16

Fig. 3.3 : PE chains and graphite walls(1D pattern surface).

Fig. 3.4 : PE chains and graphite walls(2D pattern surface).

17

4

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

4.1

はく離強度

(a) セル底面に平行なグラファイト壁モデル
図 4.1 にグラファイト壁に生じる引力の変化を，また図 4.2 に最大引力を示した時のスナップショットを示す．グラファイト壁に作用する力は上下の C 原子に作用する z
方向の力を外向き正としてそれぞれ算出し，平均して評価している．表面凹凸がな
い系 (flat) と 2π の正弦波凹凸では，わずかではあるが前者のほうが引力が大きい．

10π の表面凹凸では最大引力が著しく低下している．6π の周期でも最大引力は 2πより低下するが 10π よりは大きく，10π, 6π, 2π の順番で最大引力が大きくなっている．つまりこれらから表面凹凸が小さいほうが界面強度が低下するという法則性が導かれる．ただし，表面凹凸を小さくした極限である表面凹凸がない系は最も高い界面強度を示していることから，引張方向に垂直な有効接触面の大きさに起因するものと結論付けられる．図 4.3 に界面強度の変化を，横軸にモデル名をとってまとめて示した．
図 4.4 および図 4.5 に 2D パターンの結果を示す．また，最大引力を flat，1D，2D
全てで比較したものを図 4.6 に示す．x の 1 方向に凹凸をつけた系と x, y の 2 方向に凹凸をつけた系を比べると，6π の周期の時に著しく強度が低下しているが，表面凹凸が小さいほうが強度が低下する傾向は変わらない．

(b) 底面に垂直なグラファイト壁モデル
図 4.7 にグラファイト壁に生じる引力の変化を，図 4.8 に最大引力の時のスナップ
ショットを示す．グラファイトを平行に配置した系と垂直に配置した系を比べると，

18

明らかに垂直に配置した系の方がすべて最大引力は低下している．図 4.9 に最大引力の変化を，先の平行配置の結果と比較して示す．flat，2π ，6π までは平行配置と同様に表面凹凸が小さくなると強度が低下しているが，10π では flat と同程度の強度に戻っている．これは，垂直配置したグラファイト壁では，表面凹凸の間隔が狭くなると flat と差がなくなることを意味している．
図 4.10 に x, y 方向に凹凸をつけた時のグラファイト壁に生じる反力の変化を，図4.11 に最大反力の時のスナップショットを示す．表面凹凸による顕著な傾向を見出すことは難しい．図 4.12 にはく離強度をすべて比較して示す．橙色の印が垂直配向
した 2D パターンのものである．2D パターンでは 6π の界面強度が高く，10π が低い
結果となった．1D パターンと異なり 2D の垂直配向の場合は凸部のグラファイトの構造が周期によって著しく変化するためと考えられるが，この点については後で考
察する．

19

200

Force, F , µN

100

Force,

F,

µN

200

0

100

0

-100

-100
0

Lift, Å0500(a) flat200

F , µN

100

Force,

F,

µN

500

(b) 2π200

Force,

Lift, Å

0

100

0

-100

-100
0

Lift, Å0500(c) 6πLift, Å500(d) 10πFig. 4.1 : Change in the attractive force on rigid wall (horizontal graphite wall, flat
& 1D pattern surface).

20(a) flat(b) 2π(d) 10π(c) 6πFig. 4.2 : Snapshots at the maximum attractive force (horizontal graphite wall,
flat & 1D pattern surface).

21

Maximum attractive force,

F,

µN

250

200

150

100

6π
2π
Model name

flat

10 π

Fig. 4.3 : Maximum attractive force(horizontal graphite wall, flat & 1D patternsurface).

100

0

0

F,

µN

µN

0

100

Force, F ,100

Force,

Force,

200

200

F,

µN

200

-100

-100

-100
0

500
Lift, Å(a) 2π0Lift, Å(b) 6π500

0500
Lift, Å(c) 10πFig. 4.4 : Change in the attractive force on rigid wall (horizontal graphite wall, 2D
pattern surface).

22

(a) 2π(c) 10π(b) 6πFig. 4.5 : Snapshots at the maximum attractive force (horizontal graphite wall, 2D
pattern surface).

250

Maximum attractive force,

F,

µN

1D
2D

200

150

100

flat

6π
2π
Model name

10 πFig. 4.6 : Change in the maximum attractive force(horizontal graphite wall)

23

µN

0

Force,

100

F,

100

Force,

200

F,

µN

200

-100

0

-100
0

Lift, Å5000(a) flat
200

µN

0

Force,

100

F,

100

F,

µN

500

(b) 2π200

Force,

Lift, Å

0

-100

-100
0

Lift, Å0500(c) 6πLift, Å500(d) 10πFig. 4.7 : Change in the attractive force (vertical graphite wall, 1D pattern surface).

24

(a) flat(b) 2π(c) 6π(d) 10πFig. 4.8 : Snapshots at the maximum attractive force (vertical graphite wall, flat
& 1D pattern surface).

25

250

Maximum attractive force,

F,

µN

1D
ver-1D

200

150

100

flat

6π
2π
Model name

10 πFig. 4.9 : Change in the maximum attractive force

26

200

µN

Force,

Force,

100

0

F,
-100

-100
0

Lift, Å(a) 2π

500

100

F,

µN

0

200

Force,

100

F,

µN

200

0

-100

0

Lift, Å(b) 6π500

0Lift, Å500(c) 10πFig. 4.10 : Change in the attractive force on rigid wall (vertical graphite wall, 2D
pattern surface).

27

(a) 2π(b) 6π(c) 10πFig. 4.11 : Snapshots at the maximum attractive force (vertical graphite wall, 2D
pattern surface).

250

Maximum attractive force,

F,

µN

1D
2D
ver-1Dver-2D
200

150

100

flat

6π
2π
Model name

10 πFig. 4.12 : Change in the maximum attractive force

28

4.2

応力成分と界面の詳細(a) セル底面に平行なグラファイト壁モデル
PE 分子鎖に生じる応力を評価し，bond stretch(bs), bending(be), torsion(to), vander Waals(vw) の寄与に分けてその変化を図 4.13 および図 4.14 に示した．応力は z方向成分 (σzz ) についてのみ示している．いずれのグラフにおいても，vw と to がほぼ 0 の値を示していることから，準静的な引張り条件であることが確認できる．すべての成分の合計 (緑色の線，all) は基本的には前節のグラファイト壁に生じる引力と対応するが，PE 分子鎖はランダムに配向しているため，応力の z 方向成分が完全に壁面引力に一致するわけではない．例えば図 4.13 の flat と 2π を比較すると PE 分
子鎖の応力は 2π の方が僅かに大きい．その違いは主に bs がもたらしているが，圧着 → はく離過程で表面凹凸の底部における分子鎖の延伸が flat モデルより大きいため (2π 周期の底部の PE 層厚さ >flat の PE 層厚さ) と考えられる．10π の周期では

1D, 2D ともに bs, be が低下することで応力が低下する．一方，壁面引力が低下した
2D パターンは，いずれも PE 分子鎖の応力は flat, 1D より大きい．
図 4.15 および図 4.16 に最大引力より少し前の時点における界面近傍のスナップ
ショットを示す．flat, 1D パターンでは界面に沿うように PE 分子鎖が van der Waals
力により密着しているが，引張方向に力を効率的に支持できるのは flat や 2π のような垂直な平面が多い場合と考えられる．2D パターンでは界面に密着している分子鎖が少なく，界面近傍では分子鎖が大きく伸ばされている．このため，壁面に作用する引力は小さいが，PE 分子鎖の応力は大きいという結果になる．

29

(b) セル底面に垂直なグラファイト壁モデル
図 4.17〜4.18 に応力の変化を，図 4.19〜4.20 に界面近傍のスナップショットをまとめて示す．flat, 1D より 2D パターンの方が高い応力を示しているのは平行配置の時と同じである．ただし，グラファイト層を平行に配置した系と違い，垂直に配置した系では表面凹凸の寸法に対して応力の変化が小さい．これはグラファイト層が垂直に配置されているため，面ではなく線または点で PE 分子鎖と van der Waals 接合するためである．一方，図 4.19 の 1D-10π をみると他に比べて界面に PE 分子鎖が多く見られる．これは凹凸の側面にグラファイト平面が現れたためで，このため強度が上昇したものと考えられる．2D パターンでは 6π と 10π に大きな差がないので図?? の強度関係についてはさらなる検討が必要と考える．

30

200

200

bs

bs

100

all
100

P , MPabevw

Stress,Stress,P , MPa

all

0

be
vw

0

to

-100

0

to

500

1000
Lift, Å

-100

1500

0

(a) flat500

1000
Lift, Å1500(b) 2π
200

200
bs

bs

all
be

all

P , MPa

100

vw

Stress,Stress,P , MPa

100

0

vw
0
to

to

-100

0

be

500

1000
Lift, Å1500(c) 6π

-100

0

500

1000
Lift, Å1500(d) 10πFig. 4.13 : Change in the stress components of PE chains (horizontal graphite wall,
flat & 1D pattern surface).

31

200

200
bs

bs

100

all
100

be

0
to

-100

0

500

be

P , MPavw

1000
Lift, Åvw

0Stress,Stress,P , MPa

all

to

-100

1500

(a) 2π

0

500

1000
Lift, Å1500(b) 6π200bs
allStress,P , MPa

100
be

vw

0
to

-100

0

500

1000
Lift, Å1500(c) 10πFig. 4.14 : Change in the stress components of PE chains (horizontal graphite wall,
2D pattern surface).

32

(a) flat(b) 1D-2π(c) 1D-6π(d) 1D-10π

Fig. 4.15 : Magnified view at the interface (horizontal graphite wall, flat & 1D
pattern surface).

33

(a) 2D-2π(b) 2D-6π(c) 2D-10π

Fig. 4.16 : Magnified view at the interface (horizontal graphite wall, 2D patternsurface).

34

200

200

bs

bs
be
vw

Stress,Stress,P , MPa100

P , MPa

all

0

all

100

be
vw
0
to

to

-100

0

500

1000
Lift, Å

-100

1500

0

200

1500

200

bs

bs

P , MPa

P , MPa

all

100

be
vw

Stress,Stress,1000
Lift, Å(b) 2π(a) flat

0

all

100

be
vw
0
to

to

-100

500

0

500

1000
Lift, Å

-100

1500

0

500

1000
Lift, Å1500(d) 10π(c) 6πFig. 4.17 : Change in the stress components of PE chains (vertical graphite wall,
flat & 1D pattern surface).

35

200

200

bs
all
vw

0
to

0

500

be

100

be

-100

all

P , MPa

100

Stress,Stress,P , MPa

bs

1000
Lift, Å

vw

0
to

-100

1500

(a) 2π

0

500

1000
Lift, Å1500(b) 6π200

Stress,P , MPa

bs
all

100

be
vw
0
to

-100

0

500

1000
Lift, Å1500(c) 10πFig. 4.18 : Change in the stress components of PE chains (vertical graphite wall,
2D pattern surface).

36(a) flat(b) 1D-2π(c) 1D-6π(d) 1D-10π

Fig. 4.19 : Magnified view at the interface (vertical graphite wall, flat & 1D patternsurface).
37

(a) 2D-2π(b) 2D-6π(c) 2D-10π

Fig. 4.20 : Magnified view at the interface (vertical graphite wall, 2D patternsurface).

38

5

結言
CFRP の繊維ー樹脂界面強度について原子レベルから知見を得るために，グラファ

イト壁間にアモルファス PE 分子鎖を充填し，圧着したあとではく離する分子動力学シミュレーションを行った．グラファイト壁は，はく離方向に平行または垂直に配置し，(a)flat, (b) 一方向に 2π, 6π, 10π の正弦波で表面凹凸をつけたもの (1D パ
ターン)，(c) 同様に二方向に表面凹凸をつけたもの (2D パターン)，として表面凹凸の影響について検討した．得られた結果を以下に示す．

1. グラファイト壁に生じる最大引力 (はく離強度) を調べたところ，いずれも flat が最も大きいが，垂直配置した界面は平行のそれに比べ著しく低下した．

2.1D パターンでは表面凹凸が小さいほうが界面強度が低下するという傾向が見られた．ただし，グラファイト層を垂直に配置した系の 10π の周期では flat と同程度の値を示した．

3.2D パターンも表面凹凸が小さいほうが強度が低下する傾向が見られた．平行配置
では 1D>2D の傾向があったが，垂直では表面凹凸が小さくなると (6π, 10π) 逆転しており統一的な傾向は見られない．

4. 引張シミュレーションにおける PE 分子鎖に生じる応力と，界面近傍の分子鎖の変形挙動から 1〜3 の強度について考察した．PE 分子鎖は界面に沿って van der Waals
力により密着するが，引張方向に力を効率的に支持できるのは垂直な平面が多い場合と考えられる．表面凹凸を 2 方向につけた系は界面に密着している分子鎖が少な
く，界面近傍で分子鎖が大きく伸ばされているため，壁面に作用する引力は小さい
が，PE 分子鎖の応力は大きくなる．

39

5. グラファイト層を垂直に配置した系では線または点で PE 分子鎖と接合するため，
表面凹凸の寸法に対する変化が平行配置した系より顕著でなくなる．

40

参考文献
(1) 平野・草野・青木，日本複合材料学会誌，40，4，(2014)，153-159(2) 小柳・河合・萩原・渡辺，実験力学，10, 4, (2010), 407-412(3) S.J.V. Frankland and V. M. Harik, Com. Sci. and. Tec, 63, (2003), 1655-1661(4) R.Zhu, et al, Mat. Sci. and. Eng, 447, (2007), 51-57(5) B.Arash, H.S.Park, T.Rabczuk, Com. Str, 134, (2015), 981-988(6) Li,Y., et al, Com. Mat. Sci, 50, (2011), 1854-1860(7) 桑島・野間・逢坂，第 4 回計算化学シンポジウム，(1994)，53．

41

謝辞
本研究を遂行するにあたり，屋代如月教授には浅学非才な著者に対し懇切丁寧に指導していただきました．ここに心よりお礼申し上げます．本論文を完成させるにあたり，広い視野から研究全般に対して多くのご助言をいただきました内藤圭史助
教にも心より感謝いたします．また， 本研究において計算機は科学研究費助成事業
により援助していただきました．ここに感謝いたします．同輩としてともに切磋琢磨し合った堤貴文氏，寺田稜氏，西川涼一郎氏，山田泰成氏，吉原大志氏に御礼申し上げます．最後に，4 年間の学生生活を暖かく見守り精神的にも経済的にも支えて頂いた家族に感謝いたします．ありがとうございました．

42