修 士 論 文(題 目)

ポリプロピレン–グラファイト界面
接着強度に関する分子動力学研究

平成 30 年度

岐阜大学大学院
自然科学技術研究科 博士前期課程
物質・ものづくり工学専攻

氏名 坪井伶以

目次

第1章

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

1

第2章

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

4

2.1

分子動力学法の概要. . . . . . . . . . . . . . . . . . . . . . . . . . .

4

2.2

United Atom モデルによる粗視化ポテンシャル . . . . . . . . . . . .

5

2.3

力の導出 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

2.3.1

2 体間項（bond stretch，van der Waals） . . . . . . . . . . .

8

2.3.2

3 体間項（bending，inversion） . . . . . . . . . . . . . . . . .

9

2.3.3

4 体間項（torsion） . . . . . . . . . . . . . . . . . . . . . . . . 10

2.4

速度スケーリング法

2.5

高速化手法

. . . . . . . . . . . . . . . . . . . . . . . . . . . 13

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.5.1

粒子登録法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.5.2

ブロック分割法 . . . . . . . . . . . . . . . . . . . . . . . . . . 15

第3章

グラファイトの表面凹凸および積層方向の影響 . . . . . . . . . .

3.1

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

3.2

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

17

3.2.1

MLG に生じた力の変化 . . . . . . . . . . . . . . . . . . . . . 20

3.2.2

MLG の表面凹凸および積層方向の影響 . . . . . . . . . . . . . 21

3.2.3

界面の局所密度による評価 . . . . . . . . . . . . . . . . . . . . 22

i

第4章

PP の結晶方向および枝分かれの影響 . . . . . . . . . . . . . . .

4.1

結晶 PP の配向方向の影響 . . . . . . . . . . . . . . . . . . . . . . . . 26

4.2

第5章

26

4.1.1

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

4.1.2

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

非晶 PP の長鎖分岐の影響 . . . . . . . . . . . . . . . . . . . . . . . . 32

4.2.1

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

4.2.2

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

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

36

参考文献 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38ii

第 1 章 緒言
炭素繊維強化プラスチック（CFRP）は鋼やアルミニウムに比べ比強度，比弾性率で優れており，航空機や自動車などの様々な分野に応用されている．CFRP の破壊で重要となるのが，繊維と樹脂の界面強度である．
中谷らはショートビーム法により CF–樹脂界面の層間せん断強度を測定し，低温焼成の CF においては酸素官能基による化学的結合が支配的であるのに対して高温焼成のものは van der Waals（VDW）力のような物理的結合が支配的であると説明している (1) ．小柳らはカーボンナノチューブ（CNT）の引き抜き試験を行い，樹脂
メニスカスが従来の界面はく離の起点と想定されていた CNT 埋め込み部の付け根付近の応力集中を緩和させ，埋め込みの最も深い部分がはく離の起点となっている場合があることを確認している (2) ．播磨らは，マレイン酸変性 PP エマルジョンとシランカップリング剤処理を施した CFRP の繊維–樹脂間の界面せん断強度が未処理の CFRP に比べて 12%向上したことを確認している (3) ．Albertsen らは CF の表
面処理が CFRP の層間破壊靭性に及ぼす影響を調査することを目的として CFRP の引張及びせん断試験を行い，表面処理レベルの増加と共に破壊エネルギーも増加することを明らかにした (4) ．Rodriguez らは樹脂中の繊維を押し込む三次元有限要素
シミュレーションを行い，界面はく離開始時の臨界荷重は界面強度と硬化応力に依存するが摩擦係数には依存しないことを示した (5) ．Esqué-de los Ojos らも繊維の直径や繊維–樹脂間の摩擦係数の値を変化させて繊維を押し込む有限要素シミュレー
ションを行い，各パラメータがある値に達すると界面強度に影響しないことを示した (6) ．
このように多くの実験やシミュレーションがなされているが，様々な繊維と樹脂

1

の組合せを考えると，界面強度の研究にはミクロ挙動に関する基礎的な知見も必要と考える．原子シミュレーションで炭素界面及び高分子の界面を検討した例として，

Frankland らは CNT を含有させたポリエチレン（PE）を用いて分子動力学（MD）
シミュレーションを行い，CNT のアスペクト比と応力–ひずみ曲線との関係を検討している (7) ．Han らは MD 法を用いて CNT とポリメタクリル酸メチル樹脂（PMMA）の複合材料の弾性率を調べ，巨視的な混合則による計算との間に大きな差があることを示した (8) ．Yang らは CNT–ポリプロピレン（PP）複合材料中の CNT の欠陥
が材料の強度に与える影響を MD シミュレーションによって検討し，欠陥が存在しているとしても CNT–PP 間の界面接着力により複合材料の横弾性係数と縦方向の剛性は増加することを明らかにした (9) ．
電荷移動を考えない古典 MD の枠組みでは，グラファイトと高分子の相互作用は

VDW のみによって決まる．VDW 相互作用は比較的遠方まで及ぶため，CNT のカ
イラリティ等，グラファイト構造の局所的な違いは平均化されて界面強度にほとん
ど影響しないと予想される．そのためグラファイトの表面凹凸や高分子の分子鎖の界面への接着形態による界面強度への影響を調べる必要がある．Li らはエポキシ樹
脂と多層グラファイト（Multi-Layer Graphite，MLG）の積層周期セルで引張シミュ
レーションを行い，MLG の積層方向（平行または垂直）が界面強度に及ぼす影響を検討した (10) ．Rissanou らは PMMA 分子鎖を MLG で挟んだモデルを用いて MD シ
ミュレーションを行い，MLG から 2∼3[nm] 離れた分子鎖はグラファイト近傍にある分子鎖よりも密度が低いにもかかわらず運動量が少ないことを示した (11) ．
我々の研究グループでは，PE および PP を MLG で挟みこみ圧着・はく離させる
シミュレーションを行い，MLG 表面の凹凸や高分子の長さ，結晶構造の効果を検討してきた (12) ．本研究では，グラファイトと高分子の物理的相互作用に着目し，グラファイト–高分子界面の接着形態が接合強度に及ぼす影響について議論する．3 章ではグラファイトの表面構造に着目し，表面に正弦波状の凹凸をつけたモデルやグラ
ファイトの積層方向を変化させたモデルを用いてシミュレーションを行った．4 章
では，結晶 PP の配向方向を変化させたモデルや長く枝分かれさせた非晶 PP 分子2鎖モデルを用いてシミュレーションを行った．5 章では得られた結果を総括した．

3

第 2 章 分子動力学法
2.1

分子動力学法の概要

分子動力学法 (molecular dynamics method，略して MD 法) は，系を構成する各粒子についてニュートンの運動方程式mαd2 r α= Fαdt2(2.1.1)をたて，これを数値積分することにより粒子の軌跡を求める方法である (13) ．ここで，mα ，r α はそれぞれ粒子 α の質量および位置ベクトルである．粒子 α に作用す
る力 F α は，系のポテンシャルエネルギー Etot の各位置における空間勾配として次式により求められる．

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

r α (t + ∆t) = r α (t) + ∆t

(2.1.3)(2.1.4)となる．ここで，v α を時刻 t における粒子 α の速度とすると，

dr α= v α (t)dt

(2.1.5)であり，式 (2.1.1) と式 (2.1.5) を式 (2.1.3) と式 (2.1.4) に代入すると，)((∆t)2 F α (t)
3
+
O
(∆t)
2
mα
2
(
)
(∆t) F α (t)3r α (t − ∆t) = r α (t) − ∆tv α (t) +
+
O
(∆t)
2
mα

r α (t + ∆t) = r α (t) + ∆tv α (t) +4(2.1.6)(2.1.7)

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

r α (t + ∆t) + r α (t − ∆t) = 2r α (t) + (∆t)2()
F α (t)
4
+
O
(∆t)
mα

(

r α (t + ∆t) − r α (t − ∆t) = 2∆tv α (t) + O (∆t)3(2.1.8))(2.1.9)が得られる．∆t3 以上の高次項は無視できるとすると，時刻 t + ∆t での位置ベクト
ルはr α (t + ∆t) = 2r α (t) − r α (t − ∆t) + (∆t)2

F α (t)mα(2.1.10)と求められる．時刻 t と 1 ステップ前 t − ∆t の座標，および粒子に働く力が分かれば，式 (2.1.10) による座標更新を繰り返すことで原子の運動が追跡できる．

2.2

United Atom モデルによる粗視化ポテンシャル

式 (2.1.2) で示したように，粒子に作用する力は系のポテンシャルエネルギーにより決定される．従って，系のポテンシャルエネルギーをいかに精度よく評価するかが分子シミュレーションにおいて重要となる．本研究では，水素原子の挙動を陽に考えず，メチン基（CH），メチレン基（CH2 ），メチル基（CH3 ）をひとつの粒子として扱う United Atom（UA）モデルを用いて PP の粒子間相互作用を表現する．
具体的には，結合角に対して 3 体間ポテンシャル，二面角に対して 4 体間ポテンシャ
ルを考慮した次式のポテンシャル関数を用いて系のエネルギー Etot を評価する (15) ．

Etot = Φbs (r) + Φbe (θ) + Φto (ϕ) + Φvw (r̄) + Φinv (Θ)(2.2.1)

右辺各項は，分子内の炭素間の結合長 r に対する bond stretch ポテンシャル（2 体間），結合角 θ に対する bending ポテンシャル（3 体間），二面角 ϕ に対する torsion
ポテンシャル（4 体間），非共有原子間距離 r̄ に対する VDW ポテンシャル（2 体間），

CH 周りの 3 つの二面角の和 Θ に対する inversion ポテンシャル（3 体間）を表す．
福田ら (16)(17) によるポテンシャル関数を以下に示し，パラメータを Table 2.2.1 に示す．また各ポテンシャルのグラフを Fig. 2.2.1 に示した．Φvw は 12-6 の Lennard-Jones

5

（LJ）型の VDW ポテンシャルである．MLG の C 原子の VDW は Tsai らの文献 (18) を参照した．異なる粒子間の VDW 相互作用は，結合則 Aαβ =

{(

1/13

A1/13
αα + Aββ

)

}13

/2

，

Cαβ = (Cαα + Cββ )1/2 によって評価した (16) ．ここで A, C は LJ ポテンシャルのパラメータで，α, β は原子種を表す．Fig. 2.2.1(f) に剛体壁の C 原子との相互作用の強さを示した．

Φbs (r) =∑1

kr (r − r0 )2(2.2.2)2∑1
Φbe (θ) =kθ (θ − θ0 )22
Φto (ϕ) =Φvw (r̄) =Φinv (Θ) =∑(V1 cos ϕ + V2 cos 2ϕ + V3 cos 3ϕ)∑(∑(2.2.3)

Ar̄−12 − C r̄−6

(2.2.4))(2.2.5){K2 (Θ − Θ0 )2 + K8 (Θ − Θ0 )8 }

Θ = θ1 + θ2 + θ3 ， Θ0 = 3θ0(2.2.6)

Bond stretch
CH–CH2
CH–CH3

Table 2.2.1
r0 [nm]0.15400.1535

Potential parameter.
kr [kJ mol−1 nm−2 ]261100296200

Bending
C–CH–CC–CH2 –C

θ0 [deg]
111.0
113.3

kθ [kJ mol−1 rad−2 ]
669.4
748.9

V1 [kJ mol−1 ]4.686

V2 [kJ mol−1 ]2.428

Torsion
C–CH–CH2 –C

V3 [kJ mol−1 ]4.435

van der Waals
CH–CH
CH2 –CH2
CH3 –CH3
C–C（graphite）

A [kJ mol−1 nm12 ]2.971 × 10−5
2.971 × 10−5
2.971 × 10−50.222 × 10−5

C [kJ mol−1 nm6 ]3.138 × 10−3
6.904 × 10−3
10.042 × 10−3
1.437 × 10−3

InversionCHK2 [kJ mol−1 rad−2 ]144.3

K8 [kJ mol−1 rad−8 ]1757006

CH - CH2
CH - CH3

Bending potential, Φ be (θ), kJ/mol

Bond stretch potential, Φ bs (r), kJ/mol

400

300

200

100

C - CH - C
C - CH2 - C

6000

4000

2000

0

0
0.1

0.15
Bond length, r, nm

0.2

0

60

(a) Bond stretch

120
180
240
Bending angle, θ, deg

300

360

320
330
340
350
Inversion angle, Θ, deg

360(b) Bending
Inversion potential, Φ inv(Θ), kJ/mol

Torsion potential, Φ to (φ), kJ/mol

500

10

0

400

300

200

100

0

-10
0

60

120
180
240
Dihedral angle, φ, deg

300

360

310

(d) Invrsion

0

VDW potential, Φ vw(r), kJ/mol

VDW potential, Φ vw(r), kJ/mol(c) Torsion

-0.5
CH-CH
CH2-CH2
CH3-CH3
C-C (graphite)
-1

0.4

0.6
0.8
Non-bonded distance, r, nm

-0.5
C-CH
C-CH2
C-CH3
-1

1

(e) VDW for same UAs
Fig. 2.2.1

0

0.4

0.6
0.8
Non-bonded distance, r, nm(f) VDW for interface
Potential curves.

7

1

力の導出

2.3

式 (2.1.2) より，粒子 α の位置における力の i 方向成分は，

Fiα = −

∂Etot
∂Φbs ∂Φbe ∂Φto ∂Φvw ∂Φinv=− α −
−
−
−
α
∂ri
∂ri
∂riα
∂riα
∂riα
∂riα

(2.3.1)ベクトルの定義は MD のそれに従うものとする，すなわち

r αβ ≡ r α − r β(2.3.2)2 体間項（bond stretch，van der Waals）

2.3.1

粒子 α, β 間のポテンシャルが Φ(r) のように 2 体間の距離の関数で表されている
とき，原子 α に働く力はαβ∂rαβ ∂Φ(rαβ )∂Φ′ αβ ri= −Φ (r ) αβ
− α =− α∂ri∂ri ∂rαβr(2.3.3)原子 β に働く力はαβ∂Φ
∂rαβ ∂Φ(rαβ )′ αβ ri− β =− β= Φ (r ) αβr∂ri∂ri ∂rαβ(2.3.4)

また，

∂rαβ
riαβ
= αβ ,
∂riα
r
(

は rαβ

)2

∂rαβ
riαβ
= − αβr∂riβ(2.3.5)= r αβ · r αβ の両辺を riα ，riβ で偏微分することにより得られる．作用・反作用の法則から，多数の原子寄与がある場合でも粒子 α < β の組について

rαβ
Fiα = Fiα − Φ ′ (rαβ ) iαβrrαβ
Fiβ = Fiβ + Φ ′ (rαβ ) iαβr(2.3.6)(2.3.7)のようにプログラム中で加算すれば，カットオフ半径内のすべての原子寄与を考慮できる．riαβ /rαβ はベクトル r αβ の x 軸の方向余弦であり，2 原子間に働く力とその成分を図示すると Fig. 2.3.1 のようになる．

8

z

Fα

α

θ

rαβ
F

β

y

β

θ
rxαβ
cosθ = αβrFig. 2.3.1

2.3.2

Schematic of force vector between two particles α and β.

3 体間項（bending，inversion）µ − α − β の順で原子が連結されているとすると，µ を中心とする結合角 θµ ，α を中心とする結合角 θα ，β を中心とする結合角 θβ に r α は現れるので−∂Φbe∂ cos θ ∂θ ∂Φbe
= −
α
∂ri
∂riα ∂ cos θ ∂θ′′′
∂ cos θµ Φbe(θµ ) ∂ cos θα Φbe(θα ) ∂ cos θβ Φbe
(θβ )
=
+
+
α
α
α
∂ri
sin θµ∂ri
sin θα∂ri
sin θβ

(2.3.8)

ここで，µ < α < β の組について′′∂ cos θµ Φbe(θµ ) ∂ cos θα Φbe(θα )はを評価して Fiβ に加える
βα∂ri
sin θµsin θα∂ri′′∂ cos θβ Φbe (θβ ) ∂ cos θα Φbe(θα )はを評価して Fiµ に加えるµ
α∂risin θβ∂ri
sin θαことにより評価される．

cos θα =rβα rµαrαβ rαµ=rβα rµαrαβ rαµ(2.3.9)であるから，riβ ，riµ ，riα による微分はそれぞれ次のようになる．
(∂ cos θα1
rαβ rαµ αβ
αµ
=
−r
+
ri
i
rαβ rαµ(rαβ )2
∂riβ
(

)

=
)

1
rαβ

(
(

)

riαµ
riαβ
− αµ + cos θα αβ (2.3.10)rr)

1rαβ rαµ αµ
riαβ
riαµ
∂ cos θα1αβ
= αβ αµ −ri + αµ 2 ri
= αµ − αβ + cos θα αµ (2.3.11)
∂riµ
r r
r
r
r
(r )
9

(

∂ cos θα ∂ cos θα∂ cos θα
= −
+
α
∂ri
∂riµ
∂riβ

)

(2.3.12)

したがって，µ < α < β の組から θα を決定し，(fiβ

Φ ′ (θα ) ∂ cos θαΦ ′ (θα )
riαµ
riαβ
=
=
−
+
cos
θ
α
sin θα ∂riβ
rαβ sin θα
rαµ
rαβ

fiµ

Φ ′ (θα ) ∂ cos θαriαβΦ ′ (θα )
riαµ
=
−
+
cos
θ
=
α
sin θα ∂riµ
rαµ sin θα
rαβ
rαµ

(

)

(2.3.13))(2.3.14)Fiβ = Fiβ + fiβ(2.3.15)Fiµ = Fiµ + fiµ(2.3.16)(Fiα = Fiα − fiβ + fiµ)(2.3.17)を計算すればよい．力の成分の方向は Fig. 2.3.2 のようになる．

μ-eαβ
cosθeαμcosθeαβ

αβ

e

α

-eαμ

θ
β

eαμ
Fig. 2.3.2

2.3.3

Schematic of directions between 3-body term.

4 体間項（torsion）

3 体間のとき同様に，Fig. 2.3.3 のように µ < α < β < ε で構成される二面角で αが ε の位置にある場合 (α − 2)，α が β の位置にある場合 (α − 1)，α が α の位置にある場合 (α)，α が µ の位置にある場合 (α + 1) のように考えると

−

′
∂Φto
∂ cos ϕ ∂ϕ ∂Φto∂ cos ϕ Φto
=
−
=
∂riα
∂riα ∂ cos ϕ ∂ϕ∂riα sin ϕ

10

=

′
′
(ϕα−2 ) ∂ cos ϕα−1 Φto(ϕα−1 )∂ cos ϕα−2 Φto
+
α
α
∂ri
sin ϕα−2∂ri
sin ϕα−1′′∂ cos ϕα Φto (ϕα ) ∂ cos ϕα+1 Φto (ϕα+1 )++∂riα
sin ϕα∂riα
sin ϕα+1(2.3.18)′′∂ cos ϕα Φto∂ cos ϕα−2 Φto(ϕα−2 )(ϕα )の評価はを Fiε に加算，
αε∂risin ϕα−2∂ri
sin ϕα′′∂ cos ϕα−1 Φto∂ cos ϕα Φto(ϕα )(ϕα−1 )の評価はを Fiβ に加算，
βα∂risin ϕα−1sin ϕα∂ri′′∂ cos ϕα+2 Φto (ϕα+2 )∂ cos ϕα Φto
(ϕα )
の評価は
を Fiµ に加算，µ
α∂ri
sin ϕα+2∂ri
sin ϕαすることにより評価される．Fig. 2.3.3 のように aµ と aε を定義すると，
(

a

µ

=

aε =
cos ϕα =∂ cos ϕα=∂aεi∂ cos ϕα=∂aµi)

rαβ rµα −rαβrαβ rαµ αβαµ−r + r=−r+r
r αβ r µα r αβ
(r αβ )2
(
)
αβ βε
r
r
rαβ
rαβ rβε αββε−rβε + r βε αβ βε
=
−r
+
r
r r
r αβ
(r αβ )2aε aµaε a{µ}1 aµi
aεi− cos ϕα ε ≡ Aεi
aε aµ
a
{
}
1 aεi
aµi
− cos ϕα µ ≡ Aµi
µ
a
aε
a
αµ

αµ

(2.3.19)(2.3.20)(2.3.21)(2.3.22)(2.3.23)ベクトル Aµ ，Aε は Fig. 2.3.4 の向きとなる．これより∂aεj ∂ cos ϕα ∂aµj ∂ cos ϕα∂ cos ϕα
= α
+ α
∂riα
∂ri ∂aεj∂ri ∂aµj
右辺第一項の計算は以下のようになる．


(2.3.24)


αβ αβ
 r βε r αβ∂aεj ∂ cos ϕαrαβ rβε  ri rj
i
j

=
+
−2
+
δ
Aε
ij
 (r αβ )2 j∂riα ∂aεj(rαβ )2(rαβ )2rβε εrαβ rβε ε=
A = cos θβ αβ Ai(rαβ )2 ir(2.3.25)(2.3.26)

ベクトル Aε と rαβ は直交するので，上式の左大括弧のうちクロネッカー δij の項の
み非零となる（rjαβ Aεj =0）．同様に，式 (2.3.24) の右項は








∂aµj ∂ cos ϕαriαµ rjαβ + riαµ rjαβ rαβ rαµriαβ rjαβ


Aµ
=
−δ
+
+
−2
+
δ
ij
 j
 ij
∂riα ∂aµj(rαβ )2(rαβ )2(rαβ )2(2.3.27)= −Aµi +rαβ rαµ µ
rαµ µ
µ
A
=
−A
+
cos
θ
α αβ Aii(rαβ )2 ir
11(2.3.28)

したがって，cos ϕα の riα ，riβ ，riµ ，riε での微分はそれぞれ次のようになる．

∂ cos ϕα∂riα
∂ cos ϕα∂riβ
∂ cos ϕα∂riµ∂ cos ϕα∂riεrαβ rβε ε rαβ rαµ µ
A + αβ 2 Ai(rαβ )2 i
(r )
αβ βε
r r
rαβ rαµ µεε= −Ai + αβ 2 Ai − αβ 2 Ai
(r )
(r )
= −Aµi −(2.3.29)(2.3.30)= Aµi(2.3.31)= Aεi(2.3.32)

以上を総合すると，µ < α < β < ε の組み合せについて，以下を計算すればよい．
′
Φto
Aεi
sin ϕα′ΦtoAµ≡sin ϕα ifε ≡(2.3.33)fµ(2.3.34)rαβ rβε ε rαβ rαµ µ
f + αβ 2 fi(rαβ )2 i(r )rαβ rβεrαβ rαµ= Fiβ + fiβ = Fiβ − fiε + αβ 2 fiε − αβ 2 fiµ
(r )
(r )

Fiα = Fiα + fiα = Fiα − fiµ −(2.3.35)Fiβ(2.3.36)

Fiµ = Fiµ + fiµ(2.3.37)Fiε = Fiε + fiε(2.3.38)二面角のときと同じく，式 (2.3.22)，(2.3.23) のベクトルは αβ 軸を紙面垂直に投影
すれば Fig. 2.3.4 ような単位ベクトルに関係するものであることが分かる．αβ の軸上に投影した座標において，α 点周りのモーメントのつり合い式は時計回りを正で
考えると
)

(

fxµ (−rαµ cos θα ) − rαβ fxβ − rαβ − rβε cos θβ fxε = 0(2.3.39)

これより

rβε
rαµ
fxβ = −fxε + cos θβ αµ fxε − cos θα αβ fxµrrとなり，式 (2.3.37) と一致することが確認できる．

12(2.3.40)ε-rαμcosθα
φ

aε

α

θβ

μ

αβ

a

r

θα

μ

Fig. 2.3.3

βrβεcosθβSchematic of dihedral angle φ and projection vector aµ and aε .

Aμ-cosφaμ
μ

-aε
φ
-aμ

aμ
αβ
Fig. 2.3.4

2.4

aε

Aε
φ
ε

-cosφaεSchematic of directions between force vectors of 4-body term.

速度スケーリング法

分子動力学における温度制御には時間をスケーリングする自由度を導入し連成させて解く Nose-Hoover 法 (19) と，ここで説明する速度スケーリング法が用いられる．
速度スケーリング法は，統計熱力学より導かれる次の原子の運動エネルギーと温度の関係を用いて，以下のように制御する．

1 α α2 3
m vi = kB T
22

(2.4.1)

ここで mα は粒子 α の質量，viα は温度 T での粒子 α の速度，kB は Boltzmann 定数αとおで kB = 1.38 × 10−23 [J/K] である．目標の温度 T0 における原子 α の速度を vi0

13

くと
√
α
vi0
=

3kB T0mα(2.4.2)と表される．同様に，温度 T の時の原子 α の速度は
√

viα =

3kB T
mα

(2.4.3)と表される．式 (2.4.2)，(2.4.3) より以下の式が得られる．
α
vi0
=
viα

√

T0
T

(2.4.4)

つまり，系の温度を T から T0 にするには，式 (2.4.4) の右辺を現在の速度に掛けてやればよい．ただ，これだけでは数値積分に反映されないので，Verlet 法における座標更新の式 (2.1.10) を以下のように置き換える必要がある．
√
ααr (t + ∆t) = r (t) +

(

α
T0
α
α
2 F (t)r (t) − r (t − ∆t) + (∆t)T
mα)(2.4.5)平衡状態では，Nose-Hoover 法 (19) によって得られるカノニカルアンサンブルに一致することが示されている．

2.5

高速化手法

2.5.1

粒子登録法

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

rc より一回り大きい半径 rfc 内の粒子をメモリーに記憶し，その中で rc 内の相互作用を評価する方法であり，N × (rc 内粒子数 ≪ N − 1) に計算負荷が減少される．しかし，粒子登録法では rfc 半径より外の粒子が rc 内に達すると力の評価が適切でな

14

くなるので，一定のステップごとに登録粒子の更新（N × (N − 1) 回の探査）を行
わなければならない．このため，系がある程度の規模以上になると，粒子登録による高速化は登録更新の計算負荷により打ち消される．

rc

rfc

Fig. 2.5.1

2.5.2

Schematic of bookkeeping method.

ブロック分割法

別の高速化手法としてブロック分割法がある．Fig. 2.5.2(a) に示すように，シミュレートする系をカットオフ距離程度の格子状に分割し，各ブロックに属する粒子をメモリーに記憶する．着目している粒子に作用する力を評価する際には，その粒子が属するブロックおよび隣接するブロックから相互作用する粒子を探索して行う（Fig. 2.5.2(b)）．粒子が属するブロックは，粒子の位置座標をブロックの辺長 bx，

by で除した際の整数により判断できるので，ブロック登録時の計算負荷は粒子数
N のオーダーとなる．従って，粒子登録法では登録更新の負荷が大きくなるような大規模な系でも高速化が可能である．ポリマーのポテンシャルでは，共有結合部のbond stretch，bending，torsion ポテンシャルは相互作用する粒子が同一分子鎖内部
であらかじめ決まっているため，原子対を探索する必要はなく分子鎖単位での並列化による高速化も容易である．一方，VDW ポテンシャルは異分子鎖間，あるいは，

15

同一分子鎖内の 4 粒子以上離れた全粒子に対して相互作用を評価する必要があり，本研究で扱うような大規模な系では，ブロック分割による高速化が必要となる．

y
by
x

0
bx

(a) Domain decomposition

Fig. 2.5.2(b) Block position

Schematic of domain decomposition method.

16

第 3 章 グラファイトの表面凹凸および積層方向の影響
3.1

シミュレーション条件

x，y 方向は周期境界条件下，z 方向は仮想壁面条件下で，一辺 40[nm] の立方体
セル中にランダムウォーク法により PP 分子鎖を配置し，20000[step]（=2000[fs]，1[step]=0.1[fs]）緩和計算を行った．その後，2000[fs] かけてセルサイズを一辺 30[nm]まで圧縮し，さらに 2000[fs] 緩和計算を行いアモルファス PP 構造を作成した．分子鎖は 1 分子鎖あたり 502[united atoms（UA）] のものを 1993 本，総粒子数 1,000,486[UA]としている．Fig. 3.1.1 のように，作成した PP モデルの z 方向上下に MLG を配置しシミュレーションモデルとした．図の PP は分子鎖ごとに色分けしており，右図はランダムに作成した分子鎖の最初の 10 本のみ表示させたものである．MLG は拡大図で示したような六印環を繋げたグラフェンシートを x，y 方向に交互にずらして重ねて作成している．使用した MLG モデルの一覧を Fig. 3.1.2 にまとめて示す．MLG の積層方向は，x 軸に対して（a）水平に配置したもの，（b）54[deg] 傾けたもの，（c）垂直に配置したものの 3 パターンを考慮した．
（b）に関しては，傾けた際に界面に水平に炭素原子が並ぶように調節したために，54[deg] としている（Fig. 3.1.3）．また各積
層方向の MLG にたいして，凹凸を付けていないもの（flat），x 方向に 4π ，8π 周期の正弦波（波長 λ）で凹凸を付けたものも検討している．片方の MLG の総原子数は

148,606∼221,368[UA] である． MLG は剛体壁として扱い，毎ステップ 5.0×10−5 [nm]の速度で MLG の厚さ hwall 分圧着（上の壁を −hwall /2，下の壁を hwall /2 移動）し，その後移動方向を反転させ同じ速度ではく離させる MD 計算を行った．温度は 293[K]としている．

17

2.1[nm]

30[nm]
MLG

30[nm]

PP

top of view

z

z
x

x

y
Fig. 3.1.1

λ = ∞（flat）

Simulation model.

λ = 15[nm]（4π ）λ = 7.5[nm]（8π ）（a）horizontal（0[deg]）（b）diagonal（54[deg]）zx（c）vertical（90[deg]）

Fig. 3.1.2

Stacking direction and surface pattern of MLG.

18

distance between graphenegraphene sheets
[n
.35

m]

0

54[deg]interfaceFig. 3.1.3

Schematic of diagonal direction angle for MLG.

19

3.2

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

3.2.1

MLG に生じた力の変化

Fig. 3.2.1 に積層方向 0[deg]・凹凸 flat の MLG で圧着・はく離を行った時の MLGに作用する z 方向の力の変化の例を示す．横軸がシミュレーションステップ数で縦軸が MLG に働いた力を表している．赤の実線が下の MLG，青の一点鎖線が上の MLGに生じた力であり，圧着時にはそれぞれ反対方向に大きな反力を受けている．はく離過程では，圧着時の力に比べると極めて小さいものの，反対方向にピークを生じ
はく離に対する抵抗力が発生していることが分かる．22500[step] 近傍ではく離抵抗
力のピークを示した後減少し，400000[step] 近傍ではほぼ零となり，分子鎖が MLG
からはく離したことが分かる．Fig. 3.2.2 は（a）22500[step] と（b）700000[step] の時のスナップショットである．MLG の移動に伴って PP が引き伸ばされ，界面近傍の分子鎖が配向している様子が観察できる．
0.5
Maximum force at debonding

Force on upper MLG
Force on MLG, Fz, µN

Force on MLG, Fz, µN

2

Maximum force at debonding
0

Force on lower MLG

0.25

0

Force on lower MLG
-0.25

-2
Start debonding

Start debonding
0
21000

200000
Step, t, 0.1fs（a）throughout

Fig. 3.2.1-0.5400000021000

200000
Step, t, 0.1fs

400000（b）zoom up of the force on lower MLG

Change in the force on MLG（horizontal, flat）.

20

（a）22500[step]xz（b）700000[step]

Fig. 3.2.2

3.2.2

Snapshots at (a) max. force and (b) after debonding.

MLG の表面凹凸および積層方向の影響

Fig. 3.2.1 における上下の MLG に生じたはく離抵抗力の最大値（Max. force at debonding）の平均を MLG 界面の投影面積（30[nm]×30[nm]）で除した値を界面接着強度とし，表面凹凸による変化をまとめたものを Fig. 3.2.3 に示す．横軸が MLG の表面パターン，縦軸が MLG に働いた z 方向の応力を表す．前報 (12) で示したように，flat
が最も大きな値を示し，凹凸が細かくなるほど強度が低下する傾向が見られる．特に，水平に積層したものは他のモデルに比べていずれの表面でも強度が高く，また凹凸による変化が顕著である．一方，MLG の「端部」で PP と接する diagonal とvertical は，凹凸による変化に同様の傾向が見られるもののその変化は小さい．flat，
λ=15[nm] では diagonal が vertical よりわずかに強度が高いが，λ=7.5[nm] では逆転
している．

21

200

100

3.2.3

8

Maximum stress, σzz, MPa

300

0Fig. 3.2.3

horizontal
diagonal
vertical

15
Surface length, λ , nm

7.5

Change in the debonding stress with the surface pattern.

界面の局所密度による評価

Fig. 3.2.4 は，最大応力点近傍である 22500[step] において MLG の炭素原子に接触している粒子のカットオフ半径内に含まれる原子を局所密度として算出し，それを界面全体で平均したものである．horizontal は表面凹凸があると密度が著しく低下
する．一方 diagonal と vertical を見ると，界面強度が逆転した λ=7.5[nm] において，局所密度はいずれも λ=15[nm] よりも大きくなっている．Fig. 3.2.5 は，PP 分子鎖のVDW 半径内にある MLG の炭素原子数を示したものである．horizontal は flat が当
然ながら多いが，λ=15[nm] と 7.5[nm] でまちまちになり，diagonal と vertical は flat
から小さな表面凹凸になるほど接着 C 原子の数が多くなっている．Fig. 3.2.6 に示す
ように，horizontal は表面凹凸によって連続したグラフェンシートの平面部分が削られ MLG 層間の空洞が露出するのに対して，diagonal，vertical では凹凸によってグラフェンシートの平面部分が少し増加する．このため diagonal，vertical では flatの MLG は PP と接着している原子数が少なく，表面凹凸が小さくなるにつれ増加す
る．一方，Fig. 3.2.4 の局所密度の傾向については，PP 分子鎖の入りにくさのため，
表面凹凸をつけると低下する．diagonal，vertical の λ=7.5[nm] で局所密度が増加し22

たのは，露出したグラフェンシート側面の VDW 半径内に入る PP 分子鎖が増えた
ためで，これらは引張に対して抵抗力を持たず，このため局所密度と界面強度の傾向が逆転したものと考える．
1.1

Local density, ρ, g/cm3

horizontal
diagonal
vertical
1.05

1

0.9

Fig. 3.2.4

8

0.95

15
Surface length, λ , nm

7.5

Change in the local density at the PP/MLG interface near max. force.

23

110000

100000

90000

80000

Fig. 3.2.5

8

Number of united atoms, NUA, UA

horizontal
diagonal
vertical

15
Surface length, λ , nm

7.5

Change in the number of MLG atoms at the PP/MLG interface nearmax. force.

24

flat

λ=7.5[nm]（a）horizontal

flat

λ=7.5[nm]（b）diagonalflat

λ=7.5[nm]（c）verticalFig. 3.2.6

Snapshots of the interface of MLG.

25

第 4 章 PP の結晶方向および枝分かれの影響
4.1

結晶 PP の配向方向の影響

4.1.1

シミュレーション条件

Fig. 4.1.1 に示すように 31 らせん構造をとるアイソタクチック PP 分子鎖 4 本を 1
ユニットとし，x，y 方向は周期境界，z 方向は壁面条件下で，20[nm]×20[nm]×20[nm] のセル中にユニットを平行に並べて結晶 PP を作成した．分子鎖の配向方
向は，壁面垂直方向から 0[deg]，30[deg]，60[deg]，90[deg] 傾けた計 4 パターンとした（Fig. 4.1.2）．図では分子鎖毎に着色しており，グラデーションの角度が分子鎖の方向に対応する．また PP0 のように各モデルを称している．作成した結晶 PP の0.613[nm]2.096[nm]配向方向，分子鎖数，分子鎖長さ，総粒子数を Table 4.1.1 にまとめて示す． 作成しyx

0.656[nm]Fig. 4.1.1

31 screw chain

Crystalline PP unit.

26

z
x

（a）PP0（b）PP30（c）PP60（d）PP90

Fig. 4.1.2

Orientation direction pattern of crystalline PP.

Table 4.1.1

Chain length，number of total chains and UAs.
Orientation
Length ChainsUAs
direction[deg]02951080 318,600
30
340
936
318,24060
589
540
318,06090
295
1080 318,600

Model
PP0
PP30
PP60
PP90

た PP の上下に MLG を配置し 60000[step] 緩和した後，3 章と同じ条件で圧着・はく離シミュレーションを行った．MLG には Fig. 4.1.3 に PP0 の例で示すような波長λ = 10[nm]（周期 4π ）の正弦波で凹凸を付けている．

4.1.2

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

はく離過程中に MLG に働いた最大応力を Fig. 4.1.4 にまとめて示す．各モデルでの最大抵抗力は PP0，PP30 では差がなく，PP60，PP90 では大きくなっていった．

Fig. 4.1.5 に示したように，配向方向ごとの局所密度は 30[deg] 以下と 60[deg] 以上で大きく差が生じている．これは，分子鎖が傾くにつれて MLG との接着形態が点での接触から線での接触に変化していき，界面接着性が向上したことを示している．

Fig. 4.1.6 は，圧着終了時の PP 分子鎖を 5 本のみ表示させ界面近傍の接着形態を示したスナップショットである．分子鎖末端付近のみが接着している点での接触（PP0，PP30）は，PP60 や PP90 の線での接触に比べて強度が低下する．Fig. 4.1.7 は，MLG
への VDW 力が零になった瞬間（すなわち完全に離れた時）のスナップショットであ27

20[nm]λ=10[nm]2.1[nm]20[nm]

Fig. 4.1.3

Simulation model.

る．PP90 は他のモデルに比べて大きく引き伸ばされ，Fig 4.1.7(d) に円で示したよ
うに上下の界面ともに PP 分子鎖内部で分子鎖同士が乖離しているのが分かる．分子鎖同士の結合は弱い VDW 力によるものなので，分子鎖が層状に積み重なっている PP90 は他のモデルに比べてはく離方向の力に弱く，大きく変形する．なお，ス
ナップショットは周期境界条件下で表示させているため，分子鎖末端の切れ目が隙間のように見えている（Fig 4.1.7(d)）．

28

Maximum stress, σzz, MPa

250

200

150

100

50

0

Fig. 4.1.4

PP0

PP30
PP60
Orientation direction pattern

PP90

Change in the debonding stress with the orientation direction pattern.

1.3

Local density, ρ, g/cm3

1.25
1.2
1.15
1.1
1.05
1

Fig. 4.1.5

PP0

PP30
PP60
Orientation direction pattern

PP90

Change in the local density at the interface with the orientation direction pattern.

29

Fig. 4.1.6（a）PP0（b）PP30（c）PP60（d）PP90

Snapshots of the PP/MLG interface（5 PP chains）.

30

（a）PP0（b）PP30（c）PP60（d）PP90zxFig. 4.1.7

Snapshots of the crystaline PP after debonding.

31

4.2

非晶 PP の長鎖分岐の影響

4.2.1

シミュレーション条件

3.1 節と同様の条件で PP モデルを作成するが，ここでは Fig. 4.2.1 に示すように
PP 分子鎖を枝分かれさせたモデルも用いる．MLG の積層方向は界面に水平とし，表面凹凸は flat と波長 λ=15，7.5[nm] としている．分子鎖長さや分岐の有無の影響を比較するため，Table 4.2.1 に示すような複数の条件で作成した．分岐させていな
い PP500，PP1000，PP2000 は，1 分子鎖あたりの粒子数を約 500，1000，2000 とし，密度が一定になるように分子鎖数を調節した．分岐させたモデル（モデル名の末尾に「b」を付けたもの）は，Fig. 4.2.1(b) に示すように分子鎖末端から n/4 離れたところから n/4 の長さだけ分岐鎖を設けて作成した．このため 1 分子鎖あたりの粒子数は分岐無しのものと比べて約 1.5 倍になっている．

Table 4.2.1

4.2.2

Chain length，number of total chains and UAs.
Model Length Chains
UAs
PP500
502
1993 1,000,486
PP10001003
998
1,000,994
PP20002002
500
1,001,000
PP500b754
1330 1,002,820
PP1000b1498
667
999,166
PP2000b2998
334
1,001,332

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

Fig. 4.2.2 に分子鎖長さと分岐による最大応力の変化を示す．分岐無しのモデルは分子鎖が長くなるほど最大応力が減少した．分岐有りのモデルではいずれの MLG
凹凸パターンにおいても，分岐無しのモデルより高い値を示している．また分岐有りの場合，長さによる差は flat ではほとんどなく，λ=15，7.5[nm] でも PP1000b 以
外はほぼ同じである．小さい表面凹凸ほど強度が低下する傾向はあるが，PP1000bのみ λ=15[nm] と 7.5[nm] であまり低下していない．

32

branched point
n/4
n/4

n

n

（a）non-branched PP chain

Fig. 4.2.1（b）branched PP chain

Schematic of non-branched/branched PP chain models and snapshotsof PP chain.

Fig. 4.2.3 に各分子鎖ごとの MLG 凹凸パターンによる局所密度の変化を示した．
分岐無しのモデルは，分子鎖長さによる局所密度の差が著しく小さく，それぞれのMLG 凹凸パターンごとにほぼ同じ値をとっている．分岐有りの PP は，枝分かれのため界面近傍の分子鎖が他の分子鎖の接着を阻害する効果がより大きくなり局所密度が分岐無しより低くなっている．ただし，PP1000b だけは flat 表面のときの局所密度が分岐無しの系より高く，λ=15[nm] では分岐無しと同じ値になり λ=7.5[nm] では分岐無しと有りの中間となっている．界面に触れている分子鎖数の傾向は他と変わらなかったため，偶然なのかそれとも分岐数 250[UA] の長さが MLG の表面凹凸の特性長さと関係したものなのか，現時点では不明である．Fig. 4.2.4 は，MLG の

VDW 半径内に 300 粒子以上存在している分子鎖のうち，10 本のみ表示させたスナッ
プショットである．分子鎖が長くなると，一本の分子鎖が MLG 界面を占める割合が増える．それにも関わらず分岐無しのときに分子鎖が長くなるにつれ強度が低下したのは，本シミュレーションでは総粒子数を統一しているため，分子鎖が長くなると接着する分子鎖数が減り，MLG 近傍での絡み合い点が減少したためと考えられ
る．分岐がある場合の PP2000b は長い分子鎖にも関わらず MLG に接している割合

33

は PP500b と同程度である．PP1000b は，特に上の凹凸で広がって凹部をふさいだ
形態の分子鎖が見られる．
PP500
PP1000PP2000PP500bPP1000b
PP2000b

Maximum stress, σ zz, MPa

400

300

200

0

Fig. 4.2.2

8

100

15
Surface length, λ , nm

7.5

Change in the debonding stress with chain length and branch.

34

Local density, ρ, g/cm3

1.2
PP500
PP1000PP2000PP500bPP1000b
PP2000b

1.1

1

0.8

Fig. 4.2.3

8

0.9

15
Surface length, λ , nm

7.5

Change in the local density at the PP/MLG interface near max. force.

（a）PP500（b）PP1000（c）PP2000（d）PP500b（e）PP1000b（f）PP2000b

Fig. 4.2.4

Snapshots of 10 PP chains in VDW radius of C atoms in MLG.

35

第 5 章 結言
炭素繊維強化プラスチックの強度の鍵となる「炭素繊維−樹脂」の接着強度について原子レベルでの知見を得ることを目的として，ポリプロピレン（PP）を多層グ
ラファイト（Multi-Layer Graphite，MLG）の剛体壁で挟み込み圧着・はく離させるシミュレーションを行い，グラファイト–高分子界面の接着形態が接合強度に及ぼす影響を議論した．

2 章では，分子動力学法および PP の United Atom（UA）モデルについて概説した．また UA モデルにおける各ポテンシャルの力の導出について示した．

3 章では，アモルファス PP に剛体 MLG を圧着・はく離させるシミュレーション
を，MLG の積層方向（グラファイト層が界面に平行，垂直，斜め）および表面の凹
凸（平滑と正弦波で凹凸をつけたもの）を考慮して行った．表面凹凸がある場合，表面凹凸が小さくなるほど PP 分子鎖が入り込みにくくなり接着強度が低下した．積層方向は，分子鎖が面で接する平行配置が当然ながら高い強度を示した．PP/MLG
界面の局所密度を調べたところ，平行配置が最も密度が高く強度の傾向と一致する．
しかし，垂直および斜めの「点で接する」MLG では，表面凹凸の大きさで強度と密度が一致しなかった．これらの配置では凹凸によってグラフェンシートの平面部分が増加し接する PP 分子鎖の局所密度は上がるものの，引張に対して抵抗力を持たず容易にすべってはく離するため強度が上がらなかったと考える．

4 章では，圧着する PP ブロックを結晶 PP として分子鎖方向を圧着・はく離方向に対して 0[deg]，30[deg]，60[deg]，90[deg] と変化させたシミュレーション，ならびに分岐させた非晶 PP を用いたシミュレーションを行った．結晶 PP は界面に対して水平に配向するほど高いはく離強度を示した．これは，やはり分子鎖末端付近の36

みが接着する点での接触に比べて分子鎖全体が接着する線での接触の方が接着性が上がるためである．非晶 PP は分子鎖が長くなるほど最大応力が減少した．枝分かれがある分子鎖は，枝分かれのない分子鎖に比べて高い強度を示した．PP/MLG 界
面の局所密度を調べたところ，枝分かれのないモデルは，分子鎖長さによる差が小さく，それぞれの MLG 凹凸パターンごと同じ値をとった．このように MLG に接着している粒子数が多いにも関わらず，分子鎖が長くなるにつれ強度が低下したのは，

MLG に接着する分子鎖数が減り，MLG 近傍での絡み合い点が減少したためと考えられる．

4 章の分岐させた PP で見られるように，ケース毎に内部構造が変わって特異な挙動を生じる可能性は否定できない．したがって，本来は乱数等の条件を変えた複数の計算結果から考察することが望ましいのだが，本シミュレーションでは 150 万粒子の大規模な計算としたため，時間の都合からかなわなかった．圧着の深さや速度等による効果も考えられるため，後続の研究で検討して欲しい．

37

参考文献

(1) 中谷，中尾，複合材料と界面，40，2，pp.61 - 64，（1988）．
(2) 小柳，他，実験力学，10，4，pp.407 - 412，
（2010）．
(3) 播磨，他，科学・技術研究，5，2，pp.163 - 168，(2016)．
(4) H. Albertsen，et al. ，Com. Sci. and Tec. ，54，pp.133 - 145，
（1995）．
(5) M. Rodriguez，et al. ，Com. Sci. and Tec. ，72，15，pp.1924 - 1932，（2012）．
(6) D. Esqué-de los Ojos，et al. ，Com. Mat. Sci. ，117，pp.330 - 337，
（2016）.
(7) S. J. V. Frankland，et al. ，Com. Sci. and Tec. ，63，pp.1655 - 1661，（2003）.
(8) Y. Han，J. Eliott，Com. Mat. Sci. ，39，pp.315 - 323，
（2007）.
(9) S. Yang，et al. ，Carbon，55，pp.133 - 143，
（2013）.
(10) C. Li，et al. ，Com. A，43，pp.1293 - 1300，（2012）.
(11) A. N. Rissanou，V. Harmandaris，J. Nanopart. Res. ，15:1589，15，（2013）.
(12) K. Yashiro，et al. ，J. Soc. Mat. ，67，2，pp.242 - 248，（2018）.
(13) 上田，コンピュータシミュレーション，朝倉書店，（1990）．
(14) 洲之内，サイエンスライブラリ–理工系の数学=15，数値計算，サイエンス社，（1978）．

(15) 桑島，他，第４回計算化学シンポジウム，53，（1994）．
38

(16) M. Fukuda，S. Kuwajima，J. Chem. Phys. ，107，6，pp.2149 - 2159，（1997）．
(17) M. Fukuda，S. Kuwajima，J. Chem. Phys. ，108，7，pp.3001 - 3009，（1998）．
(18) J. L. Tsai，J. F. Tu，Mat. Des. ，31，pp.194-199，（2010）．
(19) W. G. Hoover，Phys. Rev. A，31，3，pp.1695 - 1697，（1985）．

39

謝辞

本研究を遂行するにあたり，屋代如月教授には浅学非才な著者に対し懇切丁寧に指導していただきました．ここに心より御礼申し上げます．本論文を完成させるにあたり，広い視野から研究全般に対して多くのご助言をいただきました内藤圭史助
教にも心より感謝いたします．計算機の購入に際し，科研費を充てさせて頂きました日本学術振興会に御礼申し上げます．ともに切磋琢磨し合った二村晟平氏，足立善英氏を始めとする研究室メンバーや同じ専攻の同期たち，面倒を見てくださった
先輩の方々にも御礼申し上げます．色々な方々の支えがあって今の自分があります．
最後に，6 年間の学生生活を暖かく見守り精神的にも経済的にも支えて頂いた家族に心より感謝いたします．ありがとうございました．

40