4)分子動力学計算の多成分系ガラスの機械的性質への応用

特 集
ガラスとコンピューターシミュレーション
分子動力学計算の多成分系ガラスの機械的性質への応用
旭硝子株式会社
谷
口
中央研究所
健
英
Application of Molecular Dynamics Simulation to Mechanical Property of
Commercial Soda-lime-silica Glass.
TANIGUCHI, Taketoshi
Asahi Glass Co.ltd.Research Center
1.はじめに
を組み合わせてかなり構造情報が得られるが,
ガラス物性に大きく関わっている中・長距離構
ガラスは窓やビンに代表されるように生活に
造の知見を得るのは難しい。さらに破壊のメカ
密接な材料の一つである。さらに近年における
ニズムは動的な観察が必要であるが AFM によ
光通信・情報処理産業の発達によりフォトニク
る観察[1]が報告されている程度で研究数は
スやエレクトロニクス材料としても注目されて
少ない。
いる。このように時代に対応したガラス材料の
分子動力学(Molecular Dynamics,MD)計
開発は目覚しいが,ガラスに備わる“脆さ”が
算はガラスの構造モデルを提供する手法として
問題として伴いその解明が急務である。
古くから用いられてきたが,温度・圧力の変化
ガラスの強度については多くの研究報告があ
に伴う構造変化を動的に推測する強力な研究方
るが,試験方法は押し込みや引っ掻き試験など
法である。この計算の詳細は前掲の河村先生の
多岐に亘り,試験方法により破壊のメカニズム
総説を参照されたいが,原子間のポテンシャル
が変わるので強度も異なる。そのためガラスの
を仮定して原子の移動についてニュートンの運
変形・破壊を本質的に知るために原子レベルで
動方程式を解くだけの素朴な手法である。この
の解明が必要であるが実験的には困難である。
手法を用いてすでに Soules[2]が SiO2―Na2O
なぜなら Zachariasen が示したようにガラス構
系のガラスの変形計算を Swiler[3]が歪み
造には原子の配列に周期性がなく,得られる情
速度一定での SiO2 ガラスの変形計算を報告し
報は平均的な構造を反映したもので正確な原子
ている。一方,筆者らは多成分系ケイ酸塩ガラ
配列の記載は不可能だからである。それでも X
スの変形について詳細な記載を試みている。本
線解析や赤外,ラマン,NMR などの実験手法
稿では機械的性質のための MD 計算方法を述
〒2
2
1―8
7
5
5 神奈川県横浜市神奈川区羽沢町1
1
5
0番地
7
4―7
3
4
6
TEL 0
4
5―3
7
4―8
8
6
7
FAX 0
4
5―3
Email : taketoshi-taniguchi@agc.co.jp
べたあと,MD 計算で明らかにしてきた多成分
2
8
系ソーダライムガラスの構造とその変形[4―
6]
,同じ組成の冷却過程が異なるガラスの変
NEW GLASS Vol.
22 No.
12
00
7
形[7]
,さらに低脆性ガラスの変形[8]に
て近似する。
ついて記す。
4
$ t %5
γ=γ0+γ∞81−exp(− )9
τ
2.計算方法
先に述べたように MD 計算は温度・圧力を
変化させて時間の関数として原子を動かすシミ
ュ レ ー シ ョ ン で あ る。プ ロ グ ラ ム は
ここで γ0 は弾性変形量で静水圧の場合は体
積,一軸応力の場合は歪みとなる。また γ∞は t
=∞の時の変形量,τ は遅延時間である。
体 積 弾 性 率 K0 は Birch―Murnaghan の 状 態
MXDORTO[9]を使用して行っている。系
方程式に従って求める。
は融液やガラスの場合,物性は等方性なので立
7
5
2
−
− 3
3 6$ V % 3 $ V % 37
P= K0 ( ) −( )
2 : V0
V0
;
方体のセルが基本形状となる。独立な粒子数は
計算結果に依存しないように通常は数百から数
千のオーダーである。MD 計算の結果は原子間
のポテンシャルに強く依存する。たとえば,溶
融温度の場合,井戸が深いポテンシャルであれ
ば原子同士の結合が強く初期配置から平衡状態
1
!
2
32
,
−
3
6$ V % 3−17/
.1+ (K0’
−4)
4
;1
:(V0)
0
2
!
ここで V は圧 力 P の 体 積,V0 は 1 気 圧 の
体積,K’
は体積弾性率の圧力微分である。
に移行するのに時間がかかる。従ってポテンシ
またヤング率については弾性変形量を3次式
ャルが実際よりも硬いものを使用すると通常で
に近似して0GPa における変化量として算出
は考えられない1
0
0
0
0K 以上の温度を与えて原
した。
子を動かして平衡状態に移行させる必要がある
なお粘弾性の挙動は通常のスケーリングによ
(例えば[1
0]
)
。筆者らは松井が開発した 2 体
る温度・圧力の制御では生じないので本研究で
間ポテンシャル[1
1]を使用している。このポ
の温度と圧力の制御はそれぞれ Nose[1
3]と
テンシャルは緩やかなので溶融温度は2
0
0
0K
Andersen[14]の方法で行った。
以上で十分であり,結晶構造や体積弾性率にフ
ィットするように経験的に求められているので
物質の機械的性質を調べるのに適している。
シミュレーションの結果がどの程度実際の状
態を再現しているのか確かめることは重要であ
る。この確認には X 線動径分布のほかに機械
的性質に関しては体積弾性率やヤング率が適当
圧力の設定については,設定圧力を Pset とす
るとテンソル
"P
&P
*P
xx
Pxy
xy
Pyy
xz
Pyz
Pxz#
Pyz'
Pzz +
で表現される圧力成分のうち,静水圧の場合
は以下のように設定した。
である。結晶や融液の体積弾性率の MD 計算
Pxx = Pyy = Pzz = Pset
が行われている([1
2]など)が,これらは力
Pxy = Pxz = Pyz=0
を加えると弾性変形するのに対しガラスは高分
子のような粘弾性として振舞う。例として図1
に引っ張り応力を加えたときのアルミノケイ酸
塩ガラスの歪変化を示す。応力を加えた直後は
一軸応力の場合は Pxx,Pyy,Pzz のうち 1 成
分のみ Pset とした。
3.多成分ソーダライムガラスの構造計算
弾性変形で歪みが急激に大きくなり,その後は
最初に多成分ソーダライム(SL)ガラスの
遅延弾性変形と粘性変形により緩やかに変化す
構造について述べる。ガラス組成は 6 章で述
る。この変化は Maxwell モデルと Voigt のモ
べる低脆性(LB)ガラスの組成とともに表1
デルを組み合わせた粘弾性モデルにより以下の
に掲げる。このガラスの冷却過程は図2に示す
式のように示す時間 t における変形量 γ によっ
が,原子をランダムに配置させたあと3
0
0
0K
2
9
NEW GLASS Vol.
22 No.
12
0
0
7
単位:mol%
4.多成分ソーダライムガラスの変形計算
次に SL ガラスの変形について述べる。この
ガラスに一軸応力の圧縮と引っ張り応力を加え
る(図5)と,それぞれ3GPa 以上で歪みが大
図1 引っ張り応力によるガラスの変形
きくなり応力解放後も変形は戻らなくなった。
よって3GPa 付近が降伏応力と考えられる。
で十分に平衡状態になるまで計算し徐々に3
0
0
Marsh の理論[1
6]から得られるソーダライ
K まで温度を下げた。温度に対するネットワー
ムガラスの降伏応力は3.
3GPa であるので MD
クのリング数は図3に示されるように変化し
計算は応力に対する変形をよく再現していると
た。このガラスの網目形成イオンを Si と Al と
いえる。また弾性変形領域の歪からヤング率を
したときに架橋酸素数は3.
3で組成から算出さ
算出すると6
8GPa であった。さらに体積弾性
れる値に等しかった。
率については静水圧をかけたときの弾性変形領
図4に3
0
0K のガラス構造を示す。ネ ッ ト
域の体積変化から4
7GPa と得られた。これら
ワークは様々なサイズのリングで構成されてお
の値も実験値(ヤング率:7
4GPa,体積弾性
り,河村ら[1
5]が報告しているリング数のカ
率:4
6GPa[1
7]
)に近い。
ウント法に従うとサイズが大きくなるほど数は
SL ガラスに 5 GPa の引っ張り応力を加える
増える。また図から修飾イオンは不均質に分布
と破壊した。図6に時間に対する歪と体積の変
していることがわかる。
図2 SL ガラスと急冷(SL―q)ガラスの冷却過程
3
0
図3 冷却過程における SL ガラスのネットワークリン
グの変化
NEW GLASS Vol.
2
2 No.
12
0
07
図4 SL ガラスの構造(3x3x1nm)とネットワークリング
図6 SL ガラスの引っ張り応力(5 GPa)下の体積と
歪みの時間変化
図5 SL ガラスの一軸応力による歪み変化
化を示す。本研究での破壊は急激な体積膨張と
こると緩和されて狭くなった。Onb は時間とと
定義しており SL ガラスの破壊時間は9
5ps で
もに徐々に増えていくが変動は大きい。これは
あった。破壊に至るまでに弾性変形(0∼2.
5
ネットワークの切断と再結合が繰り返されてい
ps)
,
体積が変わらず歪が変化する流動変形(2.
5
るためで,このようなネットワークの繋ぎ換え
∼8
2ps)
,
そして流動の伴った膨張(8
2∼9
5ps)
が流動の原因と考えられる。ネットワークリン
が生じた。このときに起こる主な構造変化とし
グは 3・4 員環が時間とともに増加するが 5
て 図7に Si―O―Si 角 度,非 架 橋 酸 素(Onb)数
員環以上の大きなリングは膨張の開始とともに
およびネットワークリングの変化を示す。Si―
減少し骨格構造が崩壊していくことがわかる。
O―Si 角度は弾性変形時には広がるが流動が起
図8に原子配列の変化と半径 1 nm の空隙を
3
1
NEW GLASS Vol.
22 No.
12
0
0
7
図7 引っ張り応力(5 GPa)下の SL ガラスの!Si―O―Si 角度,"酸素数1
0
0あたりの非架橋酸素,#網目形成(Si
と Al)イオン100あたりのネットワークリングの変化
図8 引っ張り応力(5 GPa)下の SL ガラスの原子配列(上図)と空隙(下図)の変化
最小とした時のガラス中の空隙変化を示す。8
2
のリング数分 布 は ガ ラ ス を 徐 冷 し た と き の
ps までに形成された空隙が95ps までに膨張す
1
0
0
0K の分布に近く,SL―q ガラスはこの温度
る様子が伺える。この空隙の成長は Celarie ら
の構造に近いと考えられる。SL―q ガラスのヤ
[1]が AFM でガラスに引っ張り応力を与え
ング率は5
6GPa と得られ SL ガラスより小さ
たときに観察された空隙変化と対応すると思わ
れる。
5.冷却過程が異なる多成分ソーダライム
ガラスの変形
次に SL ガラスと同組成で急冷した(SL―q)
ガラスについての構造と変形について述べる。
ガラスを急冷すると脆性が低くなることが知ら
れている[1
8]
。MD 計算では SL ガラスより
さらに急冷したガラスを得るために3
0
0
0K で
0
0
0
0K/fs で冷
緩和させた融液を3
0
0K まで―1
3で SL
却した。SL―q ガラスの架橋酸素数は3.
ガラスの値と等しい。しかしネットワークリン
グは 3・4 員環が多く 5 員環以上の大きいリ
ングが少なくネットワークが SL ガラスに比べ
余り発達していないことが伺える(図 9)
。こ
3
2
図9 3
00K における SL―q と SL ガラスのネットワー
クリング
NEW GLASS Vol.
2
2 No.
120
0
7
図10 5 GPa の引っ張り応力による SL―q と SL ガラ
スの歪み変化
図1
1 5 GPa の引っ張り応力下の SL―q と SL ガラス
の非架橋酸素(Onb)の変化
くなった。これは大きなリングが少ないのでネ
スの変化と共に示す。SL ガラスと同様に破壊
ットワークが柔らかくなったことを意味してい
までに大きく増減しているものの全体としては
る。
増加する傾向にない。さらに破壊直前の同じ歪
SL―q ガラスは SL ガラス同様に 5 GPa の引
みの構造を比較すると SL ガラスでは空隙が生
っ張り応力を加えると破壊した(図1
0)が,
じているのに対し SL―q ガラスでは空隙がほと
破壊に至るまでの時間は短くなり歪みは大きく
んど生じなかった(図1
2)
。これらの結果は SL
なった。図1
1に SL―q の Onb の変化を SL ガラ
―q ガラスが SL ガラスに比べ流動しやすい構造
図1
2 歪み ε=2.
0における SL―q と SL ガラスの構造
3
3
NEW GLASS Vol.
22 No.
12
0
0
7
図13 6 GPa の!圧縮と"引っ張り応力における LB
と SL ガラスの歪みと体積変化
をしていることを示している。
図1
4 LB ガラスのネットワークリングの変化
れ4
1GPa と6
1GPa(実 測 は3
9GPa と6
8GPa
6.低脆性ガラスの構造と物性
最 後 に Na2O―MgO―CaO―Al2O3―SiO2 系 で 開
[1
7]
)となった。LB ガラスは packing density
が小さいので柔らかい構造を持つと考えられ
る。
発された低脆性(LB)ガラスの構造と変形を
LB ガラスは 6 GPa 以上加えたときに破壊
記す。このガラスは SL ガラスより SiO2 量が
した。6 GPa の引っ張り応力を掛けた時の歪
多い組成で(表 1)ブリトルネスが低い[1
7]
。
みと体積の変形を SL ガラスの変形とともに示
MD 計算で構造を描くと SiO2 に富むためにネ
す(図1
3)
。LB ガラスでは破壊時間が長く歪
ットワークリング数は SL ガラスに比べ多くネ
みも大きくなった。SL ガラスの歪みや体積は
ットワークが発達している。架橋酸素数は Si
時間とともに大きくなるのに対し,LB ガラス
と Al を網目形成 イ オ ン と す る と3.
9で あ っ
では歪みは大きくなるものの体積は5
0ps 付近
た。体積弾性率とヤング率を計算するとそれぞ
から膨張したあと7
0ps 付近で体積が収縮し始
図1
5 破壊直前の LB ガラスの構造変化
3
4
NEW GLASS Vol.
22 No.
12
00
7
め9
0ps で再び膨張して破壊に至った。このと
きのネットワークリングの変化を図1
4に示
す。時間とともに 3・4 員環が増加し 5 員環
以上の大きなリングが減少したが7
0ps を越え
ると変化が小さくなった。これはネットワーク
の崩壊が緩やかになったことを示している。図
1
5に膨張した時(7
1ps)の構造と収縮した時
(8
2ps)の構造を示す。膨張時では空隙が形成
されているが収縮時では Na が空隙を占めてい
ることがわかる。LB ガラスでは SL ガラスよ
り非架橋酸素が少なく Na は SL ガラスよりさ
らに動きやすくなっているために空隙形成によ
って発生した非架橋酸素へ移動して空隙を塞
ぎ,体積が収縮したものと考えられる。
7.おわりに
本稿では MD 計算を用いて多成分系ソーダ
ライムガラスの構造とその変形,および同組成
で冷却過程が異なるガラスと組成が異なるガラ
スの変形計算結果を紹介した。MD 計算ではガ
ラス構造は原子レベルでは不均質であることを
示した。また引っ張り応力をかけると空隙が生
じて破壊するが,仮想温度の高いガラスではネ
ットワークが融液構造に近く流動を起こしやす
いので空隙が生じにくいこと,さらにネット
ワークが発達していても非架橋酸素が少ない構
造では Na は動きやすく発生した空隙の自己修
復して破壊時間を延ばすという様々な破壊過程
を推測した。ヤング率などは実験でも求められ
るが,ガラスの不均質さや極短い時間で生じる
参考文献
[1]F.Celarie,et al.,Phys.Rev.Let.,9
0,07
5
5
04
(2
0
0
3).
[2]T.F.Soules,R.
F.
Busbey,J.Chem.Phys.
78,
6
3
07
(1
9
8
3).
[3]T.S.Swiler,GlassResearcher1
1,1
0(2
0
0
2)
.
[4]T.Taniguchi,S.Ito Phys.Chem.Glasses,43
C,49
3(2
0
0
2)
.
[5]T.Taniguchi,S.Ito Phys.Chem.Glasses,
4
5,
18
3
(2
0
0
4)
.
[6]谷口健英,伊藤節郎ガラス産業連合会共済企画第
1 回ガラスシンポジウム(2
00
5)
.
[7]S.Ito,T.Taniguchi,J.Non―Cryst.Solids,
34
9,
1
7
3
(2
0
04)
.
[8]S.Ito,T.Taniguchi,Proceedings of ICG Annual
Meeting 20
0
2(2
0
02Glass ODYSSEY)Montepellier
France,
SC2,C2―6.
[9]K.Kawamura MXDORTO.JCPE Program#29
(1
99
6)
.
[1
0]T.F.Soules,
“Stochastic and Molecular Dynamic
Models of Glass Structure”; in Glass Science and
Technology 4 A.Edited by Kreidl and Uhlmann,
Academic Press,New York(199
0)
.
[1
1]M.Matsui,Geophysical Monograpf1
01.1
45
(1
99
8)
.
[1
2]M.Matsui,Geophys.Res.Letters,
2
3,
3
95
(1
9
96)
.
[1
3]S.Nose,J.Chem.Phys.
,
8
4
[1]
5
11―5
19
(1
9
84)
.
[14]H.C.Andersen,J.Chem.Phys.,7
2[4]
2
3
84―2
39
3(1
9
80)
.
[1
5]X.Qiang,K.Kawamura,T.Yokokawa,J.
Non―Cryst.Solids,1
0
4,2
61(1
98
8)
.
[16]D.M.Marsh,Proc.Roy.Soc.A,2
79,4
20
(1
9
6
4)
.
[1
7]J.Sehgal,S.Ito,J.Am.Ceram.Soc.,8
1,2
48
5
(1
9
98)
.及び J.Sehgal 私信.
[1
8]S.Ito,J.Sehgal,S.Deutschbein,Proceedings
ICG Annual Meeting20
00.
[1
9]J.Sehgal,S.Ito,J.Non―Cryst.Solids,2
5
3,1
2
6
(1
9
98)
.
破壊過程を実験的に証明するのは難しい。しか
しこれらは MD 計算でしか導き出せない利点
でもある。MD 計算が機械的性質のみならず
様々なガラスの特性を明らかに出来れば構造的
観点からの主要な材料開発の手法となる可能性
は大きい。そのためにより速い計算処理を行う
研究環境の構築とともに計算結果に左右する原
子間のポテンシャルを充実させる必要である。
3
5