卒 業 論 文(題 目)

分子動力学法による
非晶・結晶 PP/PE とグラファイト層界面のせん断強度評価

平成 28 年度

岐阜大学工学部
機械工学科

機械コース

氏名 村上 遥亮

目次
1

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

12

解析手法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

2.1

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

2

2.2

原子間ポテンシャル

3

2.3

高速化手法

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

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

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

14

4

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

175

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

32

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

34

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

35

i

1

緒言
近年，単一材料よりも優れた性能を持つ複合材料の需要が高まっている．例えば，

CFRP は金属材料と比べ比重，比強度に優れており，自動車や航空，宇宙など様々な分野に応用されている．このような複合材料の破壊現象を評価する上で重要となるのは，二材料間の界面強度である．
実験では炭素繊維をエポキシ樹脂に埋め込んだ試験片を用いた引張試験により，界面強度が温度やひずみ速度に依存することを示している．(1), (2) Wang らはマイク
ロボンド試験により，単一炭素繊維とエポキシ樹脂の間の界面特性を評価している．
(3)

また，Zu らもマイクロボンド試験により，ガラス繊維-ポリプロピレン (PP) 樹脂

界面のせん断強度が PP の剛性に依存することを予想している．(4) この他にも，プッ
シュイン試験 (5)–(7) や繊維引抜き試験 (8), (9) 等によって繊維-樹脂間の界面特性が評価されているが，界面といったそのものがナノメートル単位の現象を扱う際には，マクロよりむしろミクロの視点からの知見が必要と考える．
ミクロな視点から知見を与えるものとして，分子動力学 (MD) シミュレーションがある．MD シミュレーションを用いて高分子を扱ったものとして，Frankland らはカーボンナノチューブ (CNT) をポリエチレン (PE) 分子鎖中に埋め込んだシミュ
レーションにより，その応力-ひずみ応答を評価している．(10) 界面を評価したものとしては，Chowdhury らが PE 分子鎖中から CNT を引き抜くシミュレーションにより，界面せん断強度がマトリックスの密度に依存することを示している．(11) また，

Coto らもエポキシ樹脂から CNT を引き抜くシミュレーションにより，CNT の半径の増加に従い界面せん断強度が減少することを示してる．(12) 一方，複合材料中の高分子/炭素繊維界面を想定したシミュレーションはなされていない．そこで本研究では，グラファイト層間に非晶・結晶 PP/PE 分子鎖を充填し，圧着した後でせん断す
る MD シミュレーションを行い，界面せん断強度を評価する．

1

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

ri (t + ∆t) = ri (t) + ∆t()
dri (t) (∆t)2 d2 ri (t)
3
+
+
O
(∆t)
dt
2
dt 2

(2.3)

ri (t − ∆t) = ri (t) − ∆t()
dri (t) (∆t)2 d2 ri (t)
3
+
−
O
(∆t)
dt
2
dt 2

(2.4)

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

dr i= v i (t)
dt

(2.5)

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

r i (t + ∆t) = r i (t) + ∆tv i (t) +

2

(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)

が得られる．∆t3 以上の高次項は微小量のため無視できるとすると，時刻 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 での座標を求めるには 2 つの時刻 t と t − ∆t での座標が必要である．第一ステップの計算 (t = 0) では，t = ∆t での座標 r i (∆t) は式 (2.6) と初速度から得ることができる．これより r i (∆t) と r i (0) が既知となり，式 (2.10) を繰
り返し適用することで各粒子の座標を求められる．

原子間ポテンシャル

2.2

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

PP
Etot = ΦBS (r) + ΦBE (θ) + ΦTO (ϕ) + ΦVW (r̄) + ΦINV (Θ)

3(2.12)

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

右辺各項は，分子内の炭素間の結合長 r，結合角 θ，2 面角 ϕ，非共有原子間距離 r̄，

PP の CH 周りの 3 つの二面角の和 Θ に対するポテンシャルを表す．本研究で解析対象とする PP，PE に対しては，水素原子の運動は陽に扱わず，CH，CH2 ，CH3 を単一粒子として粗視化する united atom model により，以下の関数系およびパラメー
タが提案されている．(13)bond stretch potential (2 体間ポテンシャル)ΦBS (r ) =∑kr (r − r0 )2(2.14)

bending potential (3 体間ポテンシャル)ΦBE (θ) =∑kθ (θ − θ0 )2(2.15)

torsion potential (4 体間ポテンシャル)ΦTO (ϕ) =∑(V1 cos ϕ + V2 cos 2ϕ + V3 cos 3ϕ)(2.16)

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

inversion potential (面外角ポテンシャル)
ΦINV (Θ) =∑K2 (Θ − Θ0 )2 + K8 (Θ − Θ0 )8 (Θ = θ1 + θ2 + θ3 )(2.18)

各関数のパラメータの値を表 2.1〜2.5 に，ポテンシャル曲線を図 2.1〜2.5 にそれぞれ示した．
式 (2.2) による力の評価において，2 体間相互作用の bond stretch および van der

Waals は式 (2.19)〜(2.20)，3 体間の bending は式 (2.21)〜(2.23)，4 体間の torsion は式 (2.24)〜(2.27) のように導出される．式中の記号はそれぞれ図 2.6〜2.8 に示した．

4

⃝bond stretch ポテンシャルΦBS (r ) =∑kr (r − r0 )2

Table 2.1

Potential parameter for bond stretch.

CH2 − CHCH2 − CH2
CH3 − CHr0kr[nm] [kJ/(mol · nm2 )]0.15041.306 × 105
0.15331.373 × 105
0.15351.481 × 105

100

BS [eV]

CH2-CH
CH2-CH2
CH3-CH

50



r

0

0.2

0.4

r [nm]

Fig. 2.1

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

5⃝bending ポテンシャル
ΦBE (θ) =∑kθ (θ − θ0 )2

Table 2.2

Potential parameter for bending.

C − CH − C
C − CH2 − C

θ0
[deg]
111.0
113.3

kθ
[kJ/(mol · rad2 )]
334.7
374.7

100

BE [eV]C-CH-CHC-CH2-CH

θ

50



0

180

360

θ [deg]Fig. 2.2

Relationship between bending potential ΦBE and bending angle θ.

6⃝torsion ポテンシャル
ΦTO (ϕ) =∑(V1 cosϕ + V2 cos2ϕ + V3 cos3ϕ)

Table 2.3

Potential parameter for torsion.

V1[kJ/mol]
C − CH2 − CH − C
4.686
C − CH2 − CH2 − C
3.935

C-CH2-CH-C
C-CH2-CH2-C

10

TO [eV]

V2
V3
V6[kJ/mol] [kJ/mol] [kJ/mol]
2.428
4.435
0.0
2.177
7.786
0.0

0



gauche-10gauche
trans

0

180

360

φ [deg]Fig. 2.3

Relationship between torsion potential ΦTO and torsion angle ϕ.

7⃝van der Waals ポテンシャル
ΦVW (r̄ ) =∑(Ar̄ −12 − C r̄ −6 )Table 2.4

Potential parameter for van der Waals.

CH − CHCH2 − CH2
CH3 − CH3A
[kJ/(mol · nm12 )]2.971 × 1019
2.971 × 1019
2.971 × 1019

V[kJ/(mol · nm6 )]3.318 × 109
6.904 × 109
10.042 × 109

0.01

VW [eV]

CH-CH
CH2-CH2
CH3-CH3

0

r



r

-0.01

4

6

8

r [nm]

Fig. 2.4

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

8⃝inversion ポテンシャル
ΦINV (Θ) =∑K2 (Θ − Θ0 )2 + K8 (Θ − Θ0 )8

Table 2.5CH(Θ = θ1 + θ2 + θ3 )Potential parameter for inversion.

K2[kJ/(mol · rad−2 )]
144.3

K8
[kJ/(mol · rad−8 )]1.757 × 105

ΦINV [eV]

10

5

0
300

320

340

360

Θ [deg]Fig. 2.5

Relationship between inversion potential ΦINV and inversion angle Θ.

9

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

Φ′ (r) αβrr(2.19)

F β = −F α(2.20)

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

F γ = −F α − F β(2.21)(2.22)(2.23)

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

(2.24)(2.25)

F γ = −F α −(2.26)

Fη

(2.27)

10

i
rij

j

ri
rj
Fig. 2.6

Two molecules i, j and intermolecule vector r ij .

i
j
rik

θ

rjk

k
Fig. 2.7

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

j
rjl

φ

k

Aj

rkl
l

Ai riki
Fig. 2.8

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

2.3

高速化手法

粒子数 N の系において粒子間の全相互作用を評価すると，1step に N × (N − 1)回の計算が必要となるため，N が大きくなると極めて膨大な計算量となる．実際には，一定距離以上離れた粒子は影響を及ぼさないので，作用を及ぼす範囲 (カットオフ半径 rc ) 内の粒子からの寄与を効率よく計算することによって高速化できる. 従来よく用いられてきた高速化手法に粒子登録法がある．これは，図に示したように，

rc より一回り大きい半径 rfc 内の粒子をメモリーに記憶し，その中で rc 内の相互作用を評価する方法であり，N × (rc 内粒子数 ≪ N − 1) に計算負荷が減少される．しかし，粒子登録法では rfc 半径より外の粒子が rc 内に達すると力の評価が適切でなくなるので，一定のステップごとに登録粒子の更新 (N × (N − 1) 回の探査) を行わなければならない．このため，系がある程度の規模以上になると，粒子登録による高速化は登録更新の計算負荷により打ち消される．

rc

rfc

Fig. 2.9

Schematic of bookkeeping method.

12

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

by

y

bx
(a)Domain decomposition
Fig. 2.10x(b)Block position

Schematic of domain decomposition method.

13

3

シミュレーション方法
x，y 方向は周期境界，z 方向は仮想壁面とした立方体セル中に PP/PE 分子鎖をそれぞれランダムに成長させた非晶 PP/PE 無限薄板，ならびに，図 3.1 の単位格子を並べた結晶 PP に，剛体グラファイト壁を配置して圧着・せん断シミュレーションを行う(図 3.2 参照)．表 3.1 に各シミュレーションモデルモデルの分子鎖数と分子鎖 1 本当たりの粒子数を示す．いずれも 293K で構造緩和したときの寸法は 15nm×15nm×15nmに制御した．
まず，1step あたりに 3.0 × 10−5 nm の割合でグラファイト壁を壁の厚さ hwall 分接
近 (上の壁を −hwall /2，下の壁を hwall /2 移動) させ，分子鎖とグラファイト壁を圧着させた．グラファイト壁は (1) 表面凹凸なし，(2)x 方向に正弦波で凹凸をつけたもの (1D パターン)，(3)x，y 方向に正弦波で凹凸をつけたもの (2D パターン) について検討した．凹凸周期は 2π ，4π ，6π とした．図 3.2 に分子鎖の初期配置とグラファ
イト壁の凹凸の例を，図 3.3 および図 3.4 に 1D パターン，2D パターンの例を示す．
圧着後，上のグラファイト壁を圧着過程と同じ速度で x 方向に移動させ，せん断シ
ミュレーションを行った．

Table 3.1

Simulation model’s conditions.

Model
aPP100 Amorphous
aPP32
aPP12
aPE500aPE200cPP70
Crystal

Number ofchainsNumber of UA
in 1 chain

329
1000
2500
200
500
33

100C3 H6
32C3 H6
12C3 H6500CH2200CH2
70C3 H6

PP
PP
PP
PE
PE
PP

14

0.65nm

2.1nm

Fig. 3.1(A) Amorphous PP

Unit lattice of 31 helix iPP.

(B) Crystal PP(a) flat(a) flat(b) 2

(c) 4

Fig. 3.2(d) 6

Simulation models.

15

z

z

y

x
Fig. 3.3

Graphite wall pattern (1D).

z

z

y

x
Fig. 3.4

Graphite wall pattern (2D).

16

4

シミュレーション結果および考察
図 4.1〜図 4.3 に 1D パターンでのグラファイト壁に生じた x 方向の力の時間変化を示す．図の上部には圧着過程からせん断過程に移行した直後に生じた反力をせん
断強度として評価した点のスナップショットを示している．図 4.4 にグラファイト壁の凹凸によるせん断強度の違いをまとめて示した．まず，非晶 PP/PE，結晶 PP の違いによるせん断強度の変化を見るため，aPP100，aPE200，cPP70 のせん断強度を比較する．今回検証したどのシミュレーションモデルも，非晶 PP/PE は表面凹凸がない系はせん断に対する抵抗はほぼ 0 で，凹凸の周期が小さくなるに従ってせん
断強度が上昇している．本モデルでは凹凸の高さを同じとしているので，せん断方向の壁の勾配にせん断抵抗が依存するというごく自然な結果となる．また，高分子の種類に着目してせん断強度を比較した場合，グラファイト壁の凹凸周期に関わらず非晶 PE が最も低く，非晶 PP，結晶 PP の順に高くなるという結果が得られた．
図 4.5〜図 4.7 に 2D パターンでのグラファイト壁に生じた x 方向の力の時間変化を，図 4.8 に 1D パターンと比較したせん断強度の違いを示す．2D パターンでも，表面凹凸が小さい方がせん断強度が上昇する傾向は変わらない．1D パターンと比較
すると，どのシミュレーションモデルにおいてもせん断強度は壁面の凹凸周期にかかわらず 2D パターンが 1D パターンを下回っている．これは壁面の y 軸方向にも凹凸をつけたことによりせん断方向の壁面の面積が小さくなったことによるものと考える．
図 4.9〜図 4.11 に非晶モデルで分子鎖長さを変更した場合のグラファイト壁に生じた x 方向の力の時間変化を，図 4.12 に分子鎖長さの違いによるせん断強度の違いを示す．グラファイト壁が flat の場合，非晶 PP/PE はすべてせん断に対する抵抗は
ほぼ 0 である．凹凸をつけた場合，PP，PE 共に分子鎖長さが短い方がせん断強度が高くなる傾向が認められた．このことから，壁面から離れたバルク部分での絡み合いはせん断強度に寄与しないことがわかる．一方，凹凸がある系で分子鎖長さが17短い方が高いせん断強度を示したのは，圧着時の凹部への入り込みやすさが考えられる．これについて考察するため，上のグラファイト壁の最下層 (表面の凸部位置)
より上にある粒子数を壁面凹部に充填された粒子数として，圧着後の構造について
分子鎖長さによる粒子数の違いを評価した．図 4.13 に凹部に充填された粒子数と表面凹凸の関係を示す．PP，PE いずれも分子鎖長さの短い方が凹部に充填された粒
子数が多いことがわかる．壁面凹部に充填される粒子量が多いほど分子間の距離が縮まり，van der Waals 力が増加するため，せん断方向の抵抗力が高くなると考える．しかし，分子鎖の入り込みやすさ以外にも壁面の勾配に沿った形状変化のしやすさなどが考えられるため，より詳細な検討が必要と考える．

18press steppress step

shear step0
Shear force, F, N

Shear force, F, N

0





-0.05

-0.1

-0.05

-0.1

-0.15

-0.15
0

20000

40000 60000 80000
MD step, t, 10-1fs

010000020000(a) flatpress step

40000 60000 80000
MD step, t, 10-1fs

100000(b) 2πpress step

shear step

shear step0
Shear force, F, N

0
Shear force, F, N

shear step





-0.05

-0.1

-0.05

-0.1

-0.15

-0.15
0

20000

40000 60000 80000
MD step, t, 10-1fs

0100000(c) 4πFig. 4.1

20000

40000 60000 80000
MD step, t, 10-1fs

100000(d) 6πChange in the shear force on the GF wall (amorphous PP model，
100C3 H6 x329，flat & 1D patten surface).

19press steppress step

shear step0
Shear force, F, N

Shear force, F, N

0





-0.05

-0.1

-0.05

-0.1

-0.15

-0.15
0

20000

40000 60000 80000
MD step, t, 10-1fs

010000020000(a) flatpress step

40000 60000 80000
MD step, t, 10-1fs

100000(b) 2πpress step

shear step

shear step0
Shear force, F, N

0
Shear force, F, N

shear step





-0.05

-0.1

-0.05

-0.1

-0.15

-0.15
0

20000

40000 60000 80000
MD step, t, 10-1fs

0100000(c) 4πFig. 4.2

20000

40000 60000 80000
MD step, t, 10-1fs

100000(d) 6πChange in the shear force on the GF wall (amorphous PE model，
500CH2 x200，flat & 1D patten surface).

20

press steppress step

shear step0
Shear force, F, N

Shear force, F, N

0





-0.05

-0.1

-0.05

-0.1

-0.15

-0.15
0

20000

40000 60000 80000
MD step, t, 10-1fs

010000020000

press step

40000 60000 80000
MD step, t, 10-1fs

100000(b) 2π(a) flatpress step

shear step

0

shear step0
Shear force, F, N

Shear force, F, N

shear step





-0.05

-0.1

-0.15

-0.05

-0.1

-0.15
0

20000

40000 60000 80000
MD step, t, 10-1fs

1000000(c) 4πFig. 4.3

20000

40000 60000 80000
MD step, t, 10-1fs

100000(d) 6πChange in the shear force on the GF wall (crystal PP model，70C3 H6 x33，
flat & 1D patten surface).

21

aPP100 100C3H6 x329
aPE500 500CH2 x200
cPP70 70C3H6 x33

Maximum shear force, F, N

0.08



0.06

0.04

0.02

0
flat

Fig. 4.42
4
Model Name

6

Change in the maximum shear force.

22

press steppress step

shear step0
Shear force, F, N

0
Shear force, F, N

shear step





-0.05

-0.1

-0.05

-0.1

-0.15

-0.15
0

20000

40000 60000 80000
MD step, t, 10-1fs

010000020000(a) 2π40000 60000 80000
MD step, t, 10-1fs

100000(b) 4πpress step

shear step

Shear force, F, N

0



-0.05

-0.1

-0.15
0

20000

40000 60000 80000
MD step, t, 10-1fs

100000(c) 6πFig. 4.5

Change in the shear force on the GF wall (amorphous PP model，
100C3 H6 x329，2D patten surface).

23

press steppress step

shear step0
Shear force, F, N

Shear force, F, N

0shear step





-0.05

-0.1

-0.15

-0.05

-0.1

-0.15
0

20000

40000 60000 80000
MD step, t, 10-1fs

100000020000(a) 2π40000 60000 80000
MD step, t, 10-1fs

100000(b) 4πpress step

shear step

Shear force, F, N

0



-0.05

-0.1

-0.15
0

20000

40000 60000 80000
MD step, t, 10-1fs

100000(c) 6πFig. 4.6

Change in the shear force on the GF wall (amorphous PE model，
500CH2 x200，2D patten surface).

24press step

shear steppress step0
Shear force, F, N

Shear force, F, N

0shear step





-0.05

-0.1

-0.15

-0.05

-0.1

-0.15
0

20000

40000 60000 80000
MD step, t, 10-1fs

100000020000(a) 2π40000 60000 80000
MD step, t, 10-1fs

100000(b) 4πpress step

shear step

Shear force, F, N

0



-0.05

-0.1

-0.15
0

20000

40000 60000 80000
MD step, t, 10-1fs

100000(c) 6πFig. 4.7

Change in the shear force on the GF wall (clystal PP model，70C3 H6 x33，
2D patten surface).

25

aPP100 100C3H6 x329 (1D)
aPP100 100C3H6 x329 (2D)

aPE500 500CH2 x200 (1D)aPE500 500CH2 x200 (2D)

0.08
Maximum shear force, F, N

Maximum shear force, F, N

0.08



0.06

0.04

0.02

0.06

0.04

0.02

0

0
flat

2
4
Model Name

6(a) Amorphous PP model



0.06

0.04

0.02

0
flat

2
4
Model Name6(c) Crystal PP modelFig. 4.8

6(b) Amorphous PE modelcPP70 70C3H6 x33 (1D)
cPP70 70C3H6 x33 (2D)

0.08
Maximum shear force, F, N

2
4
Model Nameflat

Change in the maximum shear force (2D pattern).

26press steppress step

shear step0
Shear force, F, N

Shear force, F, N

0





-0.05

-0.1

-0.15

-0.05

-0.1

-0.15
0

20000

40000 60000 80000
MD step, t, 10-1fs

100000020000(a) flatpress step

40000 60000 80000
MD step, t, 10-1fs

100000(b) 2πpress step

shear step

0

shear step0
Shear force, F, N

Shear force, F, N

shear step





-0.05

-0.1

-0.15

-0.05

-0.1

-0.15
0

20000

40000 60000 80000
MD step, t, 10-1fs

1000000(c) 4πFig. 4.920000

40000 60000 80000
MD step, t, 10-1fs

100000(d) 6πChange in the shear force on the GF wall (amorphous PP model，
32C3 H6 x1000，flat & 1D patten surface).

27

press step

shear steppress step0
Shear force, F, N

Shear force, F, N

0





-0.05

-0.1

-0.15

-0.05

-0.1

-0.15
0

20000

40000 60000 80000
MD step, t, 10-1fs

100000020000(a) flatpress step

40000 60000 80000
MD step, t, 10-1fs

100000(b) 2πpress step

shear step

0

shear step0
Shear force, F, N

Shear force, F, N

shear step





-0.05

-0.1

-0.15

-0.05

-0.1

-0.15
0

20000

40000 60000 80000
MD step, t, 10-1fs

1000000(c) 4πFig. 4.10

20000

40000 60000 80000
MD step, t, 10-1fs

100000(d) 6πChange in the shear force on the GF wall (amorphous PP model，
12C3 H6 x2500，flat & 1D patten surface).

28

press steppress step

shear step0
Shear force, F, N

Shear force, F, N

0




-0.05

-0.1

-0.05

-0.1

-0.15

-0.15
0

20000

40000 60000 80000
MD step, t, 10-1fs

010000020000

press step

40000 60000 80000
MD step, t, 10-1fs

100000(b) 2π(a) flatpress step

shear step

shear step0
Shear force, F, N

0
Shear force, F, N

shear step





-0.05

-0.1

-0.15

-0.05

-0.1

-0.15
0

20000

40000 60000 80000
MD step, t, 10-1fs

1000000(c) 4πFig. 4.11

20000

40000 60000 80000
MD step, t, 10-1fs

100000(d) 6πChange in the shear force on the GF wall (amorphous PE model，200CH2 x500，flat & 1D patten surface).

29

aPP100 100C3H6 x329
aPP32 32C3H6 x1000
aPP12 12C3H6 x2500
aPE500 500CH2 x200
aPE200 200CH2 x500

Maximum shear force, F, N

0.08



0.06

0.04

0.02

0
flat

Fig. 4.124
2
Model Name

6

Change in the maximum shear force.

30

aPP100 100C3H6 x329
aPP32 32C3H6 x1000
aPP12 12C3H6 x2500
aPE500 500CH2 x200
aPE200 200CH2 x500

Number of particles, N

3000

2500

2000

2

4

6

Model Name
Fig. 4.13

Number of particles filled in the surface pit.

31

5

結言
繊維-樹脂界面のせん断強度についてミクロな視点からの知見を得るために，グラファイト壁間に非晶 PP/PE，結晶 PP 分子鎖を充填し，圧着した後で x 方向に移動させる MD シミュレーションを行った．グラファイト壁は (1)flat，(2)x 方向に 2π ，4π ，6π の正弦波で表面凹凸をつけたもの (1D パターン)，(3) 同様の凹凸周期で x，y 方向に表面凹凸をつけたもの (2D パターン) について検討した．得られた結果を以下に示す．

1. 1D パターンでグラファイト壁に生じるせん断強度を調べたところ，どのシミュ
レーションモデルでも flat が最も小さく，凹凸周期が小さくなるにしたがって上昇する傾向が見られた．

2. 分子鎖の種類別にせん断強度を比較した場合，凹凸周期にかかわらず非晶 PE
モデルが最も小さく，非晶 PP，結晶 PP の順に高くなる傾向が見られた．

3. 2D パターンも凹凸周期が小さい方がせん断強度が上昇する傾向が見られた．
また，1D パターンとせん断強度を比較した場合，凹凸周期にかかわらず 2Dパターンの方が低い値を示した．

4. 分子鎖長さの違いによるせん断強度の変化を比較した結果，PP，PE モデル共に分子鎖長さの短い方がせん断強度は上昇する傾向が見られた．

32

5. 1 および 3 より，せん断強度はせん断方向での壁の勾配，ならびに面積により決定付けられるという自然な結論が導かれた．

6. 4 で分子鎖長さの短い方がせん断強度が上昇した要因として，凹部への入り込
みやすさが考えられるため，圧着時に凹部に充填された粒子数を調べた．その結果，分子鎖長さの短い方が多く充填される傾向が認められた．このため，

van der Waals 力による反力が増加し，短い分子鎖ほどせん断抵抗が増加した
と考える．

33

参考文献
(1) J. Koyanagi, et al, Com. Str., 92, (2010), 150-154(2) 三輪・大沢・足立, 繊維学会誌, 41, (1985), T223-T230(3) C. Wang, et al, Mat. and. Des., 85, (2015), 800-807(4) M. Zu, et al, Carbon, 50, (2012), 1271-1279(5) M. Rodriguez, et al, Com. Sci. and. Tec., 72, (2012), 1924-1932(6) A. S. Childers, et al, Mat. Let., 174, (2016), 106-109(7) D. Esqu-de los Ojos, et al, Com. Mat. Sci., 117, (2016), 330-337(8) 小柳・河合・萩原・渡辺, 実験力学, 10, 4, (2010), 407-412(9) J. Koyanagi, et al, Com. PartA, 43, (2012), 1819-1827(10) S. J. V. Frankland, et al, Com. Sci. and. Tec., 63, (2003), 1655-1661(11) S. C. Chowdhury and T. Okabe, Com. PartA, 38, (2007), 747-754(12) B. Coto, et al, Com. Mat. Sci., 79, (2013), 99-104(13) 桑島・野間・逢坂, 第四回計算化学シンポジウム, (1994), 53

34

謝辞
本論文を作成するにあたり，懇切丁寧に指導していただいた屋代如月教授と内藤圭史助教に心より感謝の意を表します．また， 本研究において計算機は科学研究費助成事業により援助していただきました．ここに感謝いたします．そして，同輩としてともに切磋琢磨し合った足立善英氏，片岡裕司氏，坪井伶以氏，二村晟平氏，山内雅也氏，山下真司氏，山脇一輝氏，先輩として数多くの助言をしてくださった寺田
稜氏，西川涼一郎氏に感謝したします．最後に，4 年間の学生生活を暖かく見守り，
経済面や精神面で支えて頂いた家族に感謝いたします．ありがとうございました．

35