Journal of Computer Chemistry, Japan
Online ISSN : 1347-3824
Print ISSN : 1347-1767
ISSN-L : 1347-1767
研究論文
UO2-ZrO2固溶体の融点および熱伝導率の分子動力学解析
有馬 立身
著者情報
キーワード: 分子動力学, 熱伝導率, 融点, UO2, ZrO2
ジャーナル フリー HTML

2015 年 14 巻 4 号 p. 97-104

詳細
Abstract

原子炉過酷事故の事象進展には溶融燃料の熱的物性が大きく影響する.溶融燃料の主な物質であるUO2-ZrO2固溶体の融点を古典分子動力学法を用いて組成をパラメータとして,また熱伝導率を組成・温度を変えて固体から液体状態まで評価した.融点に関しては固液2相共存法が単相で評価するよりも実験値に近い値を与えた.これは結晶性固体の単相に対して温度を上げても融解の核が生成しにくく,過加熱状態に陥りやすいことを意味している.また,UO2とZrO2を互いに固溶させることにより融点は低下し,それが融解エンタルピーの減少と相関することを明らかにした.熱伝導率は平衡分子動力学法によりグリーン-久保の関係式から,エネルギーと電荷の流れの相互効果を考慮し,固体から液体状態までの格子振動の寄与を評価した.固体の熱伝導率には,ウムクラップ散乱による温度上昇に伴う低下,低温の固溶体で見られるフォノンの不純物散乱による低下が確認できた.一方,超高温の液体状態の熱伝導率は低く,温度および組成の依存性は小さいことが分かった.

1 はじめに

東日本大震災において福島第一原子力発電所は炉心溶融を起こした.そこでは,冷却機能を喪失した炉内にて,核燃料のUO2と被覆管のジルカロイ(Zr合金)が溶融し,構造材や制御棒を巻き込みながら,圧力容器底部に落下,一部は原子炉格納容器にまで漏れ出したと考えられている.このような過酷事故の事象進展を詳細に解析することは,今後の原子炉の安全性を確保する上で必要不可欠である.しかしながら,実規模での模擬実験・検証試験は困難であり,実験室規模であっても極めて限定的なものにならざるを得ない.それに代わるものとしてマルチスケールシミュレーション技法がある.ここで言うマルチスケールスケールシミュレーションとは,ナノスケールでの電子・原子・分子の挙動を扱う第一原理計算や分子動力学(MD: Molecular Dynamics)法,マイクロ・メソケールで原子・分子の集団を扱うことが可能なモンテカルロ法,マクロスケールを計算対象とする有限要素法等,更に実規模での原子炉の安全解析が可能なプログラム群[例えば1, 2]のことを指す.材料学的には,特に事故初期に発生した超高温のUO2とZr酸化物との混合物(固体および液体)の物性がその後の事象進展に大きく影響すると予想される.そこで本研究では,原子レベルでの材料シミュレーションが可能なMD法を用いて,UO2-ZrO2固溶体の融点および熱伝導率の評価を試みた.

UO2のMD法による研究は,融点や熱伝導率に限らず,これまでにも数多く存在する[例えば3,4,5].それに比べて純粋なZrO2のものは少ない[6].なぜならZrO2の相転移挙動,つまり大気圧で温度上昇と共に単斜晶→正方晶→立方晶と結晶構造が変化することの再現に成功していないからである.本研究では一つの試みとして,ZrO2のポテンシャル関数に室温での単斜晶構造を再現できるPedone等によって開発されたものを採用した[7].彼らはガラスのMD計算用として単純なポテンシャル関数ながら,SiO2をはじめ,NaO,Al2O3など数多くの酸化物のポテンシャルパラメータを決定している.本研究では,UO2のポテンシャルパラメータをPedone等のポテンシャル関数に合うよう決定し,MD計算に利用した.

2 計算手法

2.1 原子間ポテンシャル関数

Pedone等の開発した2体原子間ポテンシャル関数は式(1)で表される[7].   

U ( r ) = z i z j e 2 r + D i j [ { 1 exp ( a i j ( r r 0 ) ) } 2 1 ] + C i j r 12
(1)

ここで,右辺第一項はクーロン項,第二項はモース項,第三項は反発項である.ziはi番目のイオンの電荷(形式電荷の60%)rは粒子間距離,Dijaijr0Cijはイオンi-jに対するポテンシャルパラメータである.陽イオン間では,クーロン相互作用のみ働くとしている.

本研究では計算対象をUO2,ZrO2およびそれらの固溶体としており,ZrO2のポテンシャルパラメータにはPedone等によって室温での単斜晶構造および体積弾性率を再現するように決定されたものを採用した[7].UO2のポテンシャルパラメータは蛍石構造の格子定数の温度依存性および室温での体積弾性率を再現するように決定した.計算に使用したポテンシャルパラメータをTable 1に,UO2の格子定数の温度依存性をFigure 1に示す.MD計算から得られた格子定数は全温度領域にわたって実験値をよく再現できている[8].

Table 1.  Potential parameters of U-O, Zr-O[7], and O-O[7]
Dij[eV] aij−2] R0[Å] Cij[eV Å12]
U2.4- O−1.2 0.036670 2.398825 2.897745 1.0
Zr2.4-O−1.2 0.206237 2.479675 2.436997 1.0
O−1.2- O−1.2 0.042395 1.379316 3.618701 22.0
Figure 1.

 Temperature dependence of lattice constant of UO2. Solid line was obtained from measurement.

MD計算プログラムには,熱伝導率計算の機能を追加したMXDORTOを使用した[9].

2.2 融点評価法

本研究では融点評価手法として2つの方法を採用した.一つは初期状態が結晶性固体であるスーパーセルに対して温度Tを変えながら,NTPアンサンブルのもとで計算し,融解するか否かを判定するものであり,OPS (One-Phase Simulation)と呼ばれる.この方法は,融解の核が発生しにくいため過加熱状態に陥りやすく,融点が高めに評価される傾向がある.ここでは蛍石構造のユニットセル5 × 5 × 5 (= 1500粒子)のスーパーセルに対して,P = 0.1 MPa,Tを一定速度(1.25 K/ps)で徐々に上昇させながら融点を決定した.

一方,2相共存法(TPS: Two-Phase Simulation)と呼ばれる方法は,初期状態で結晶性固体と液体を隣合せて共存させておき,NTP一定の条件下で計算を走らせ,最終的に融解するか否かを判定するものである[10, 11].この方法は初期状態で液体が存在しているため,過加熱状態に陥ることはなく,融解現象だけでなく,設定温度が融点よりも低い場合は,結晶成長が観察される等の利点がある.但し,設定温度が融点に近くなると状態変化が起きにくくなるため,長時間の計算が必要な場合がある.ここでは,結晶性固体のスーパーセルと同じ大きさ,同じ温度の液体のスーパーセルを蛍石構造の(001)面で隣合せ,固液2相の共存したスーパーセルを作成した.スーパーセルの大きさは蛍石構造のユニットセル5 × 5 × 5 × 2 (= 3000粒子)とした.その後,スーパーセルを作成した時と同じ温度でNTP一定の条件下で計算した.また,OPSとTPSの両者に対して,固溶体中のUとZrは,蛍石構造の陽イオン副格子上にランダムに配置した.

2.3 熱伝導率評価法

平衡分子動力学法(EMD: Equilibrium MD)では,熱伝導率はエネルギー流の自己相関関数の積分値から,グリーン-久保の公式を用いて計算される.   

κ 0 = 1 3 k B T V 2 0 d t J ( t ) J ( 0 )
(2)

ここでkB はボルツマン定数,Vはスーパーセル体積, J ( t ) はエネルギー流である.このエネルギー流はクーロン相互作用が働くイオン体系においては,それほど単純なものとはならない.イオン結晶のものは,Bernu [12]やPierleoni [13]が一成分プラズマに対して定式化したものを基に,Sindzingre [14]やLindan [15]によって整理された.本研究の計算対象であるUO2-ZrO2 系でもクーロン相互作用が働くことから,この手法に従った.エネルギー流は以下のように表される.   

J ( t ) = i = 1 N [ m i v i 2 2 + 1 2 j i U ( r i j ) ] v i + i = 1 N v i S i (3)

ここでmii番目の粒子の質量,vi は速度である.テンソル S i は長距離でも働くクーロン相互作用の発散を避けるためにエワルド法にならって導入されたものであり,短距離力であるモース項や分散項も同時に含まれる.本計算ではMotoyama等によって与えられた S i を採用した[16].

熱伝導率の計算対象となるスーパーセルの大きさは,蛍石構造のユニットセル5 × 5 × 5とし,NTPアンサンブルのもと,P = 0.1 MPaで計算した.UO2-ZrO2固溶体のUとZrは融点評価と同様に,ランダムに蛍石構造の陽イオン副格子上に配置した.EMDでの熱伝導率は,一つの計算条件につき,緩和計算を行った後の500,000step (1 step = 2 fs)のデータをもとに算出した.

3 結果および考察

3.1 UO2の融点付近の構造変化

超高温の物性の一つとして,またポテンシャルパラメータの検証として,UO2の密度と温度の関係をFigure 2に示す.図中密度が急激に変化している温度が融点に対応しており,ここでは後で示すように,MD計算による融点と実験値が一致しなかったため,UO2の実験値の融点3120 Kとした.したがって,3120 K以上の液体の密度は,一旦5000 Kで完全に液体にし,その後冷却したものに対して求めたものである.MDの結果は,Harding等の実験値と比較しても[17],固体および液体状態で良い一致を示しており,本研究で決定したポテンシャルパラメータは密度に関して実験値とよく一致する結果を与えることが確認された.

Figure 2.

 Temperature dependence of density of UO2. Solid line was given for measurement, and dotted lines at high temperatures showed the uncertainties of measurement.

UO2液体の構造は,融点が高い,それ自体放射性物質であることから評価が難しいとされていたが,最近の浮遊式溶解法などの測定装置の開発により,少しずつ明らかになりつつある[18, 19].Skinner等は [19],高温の固体においてはUに対する酸素の配位数は8であり,室温の場合と変わらないが,融解後は6.7程度に低下することを実験的に確認している.また,温度の上昇と共にU-O間の第一近接距離は短くなり,融解後更に小さくなるが,U-U間は逆に温度の上昇と共にその距離は拡がり,融解時に大きな収縮が見られることを報告している.但し,融解後は3270 Kのみのデータしかなく,更に高温の融体の挙動はYakub等[20]の開発したポテンシャル関数を使ったMD計算により,U-O間およびU-U間距離は拡がると予想している.これに合わせた本研究のMD計算の結果をFigure 3に示す.この図は温度TにおけるU-O間およびU-U間距離を300 Kの値で規格化したものである.確かにU-O間またU-U間距離はそれぞれ温度の上昇と共に短く,また長くなり,融解直後には大きな収縮が両者に見られる.ところが,U-O間距離については,本計算では融解後も収縮する傾向が見られた.この相違はMD法のポテンシャル関数に由来するものであるが,更なる実験的な検証が必要であると思われる.

Figure 3.

 Temperature dependence of nearest neighbor distance ratio: r (T)/r (300 K). These results were obtained from U-O and U-U ions.

3.2 UO2-ZrO2の融点

UO2-ZrO2固溶体の状態図について簡単に述べておく.UO2の固体は融点(3120 K)まで蛍石構造のみをとる.一方,ZrO2は1443 Kで単斜晶から正方晶へ,2558 Kで正方晶から立方晶へと相変態し,融点は2983 Kとされる[21, 22].固溶体に対しては確定的ではないが,2558 K以上で立方晶,それ以下の温度ではUO2-richな立方晶,ZrO2-richな正方晶領域およびそれらの混合相を持ち,約1380 K以下になるとZrO2-richな相が単斜晶に変わると考えられている[21].本研究で使用したポテンシャル関数では,ZrO2の相変態挙動を再現できないため,融点近傍の結晶構造を考慮し,全てのスーパーセルを立方晶蛍石構造とした.従って,中低温およびZrO2-richな領域では,本研究のものは仮想的なものとなる.

Figure 4にZrO2 40 mol %固溶体のOPSで計算した融解前後の平均自乗変位(MSD: Mean-Square Displacement)と温度の関係を示す.ここでの昇温速度は1.25 K/psであった.この図から融点は約3980 Kと判断でき,特に陽イオン(UおよびZr)のMSDが劇的に増加している温度がそれに対応する.酸素は2500 K辺りから拡散が顕著になり,融解前も極めて大きく拡散していることが分かる.MSDで比較すると,融点直下で陽イオンに比べて3ケタ程度大きい.これは蛍石構造を持つUO2-ZrO2固溶体が酸素の超イオン伝導体であることを示しており,融解前後で酸素のMSDに大きな違いがないことから,固体状態においても液体と同じような拡散が起きていることを示唆している.

Figure 4.

 Mean-square displacements (MSD) of O, U, and Zr ions.

Figure 5にZrO2 40 mol %固溶体のTPSで計算した融点前後のスーパーセルの様子を示す.Figure 5 (a)は3200 Kの計算で使用した初期状態であり,(b)はその温度で終状態が固体,(c)は3300 Kで終状態が液体になったことを示している.これらのスーパーセルは3次元可視化プログラムVESTAを使って描画した[23].この結果からTPSでの融点は3250 Kと判断でき,同じ組成のOPSから求めたものより低い値となった.固溶体に対するOPSおよびTPSから求めた融点と組成の関係をFigure 6 (a)に示す.評価法のところで述べた通り,全組成にわたってOPSはTPSよりも750 K程度融点が高く見積もられている.これはOPSでは融解の起点となる核形成が起きにくく過加熱状態になりやすいこと,TPSではすでに固液相が共存しており,融解の障壁が低いためと考えられる.しかしながら,TPSから求まる融点はOPSより実験値に近いもののUO2-ZrO2の全組成領域に亘って300 K-600 K程度高くなった[24].また,いずれの方法でもUO2とZrO2が互いに固溶することにより,融点が低下する傾向が表れている. Figure 6 (b) にTPSから求めた融解エンタルピーを示す[25].つまり,固溶の効果により,融解の障壁となる融解エンタルピーが小さくなり,融点が低下したと解釈できる.

Figure 5.

 (a) Supercell of (U0.6Zr0.4) O2 at 3200 K as the initial configuration of TPS, (b) crystalline solid at 3200 K as the final configuration, (c) liquid at 3300 K as the final configuration. Blue: U; light blue: Zr; green: O.

Figure 6.

 (a) Melting points calculated by OPS and TPS. Dotted and solid lines stand for liquidus and solidus ones, respectively. These lines were drawn by the equations used in SCDAP/RELAP5-3D© code. (b) Enthalpies of fusion of UO2-ZrO2 solid solutions. Solid symbols were estimated by MD calculations, and the blank one was the experimental value recommended for UO2 by Konings[25].

3.3 UO2-ZrO2の熱伝導率

UO2やZrO2は超イオン伝導体であり,酸素の拡散が大きくなるような高温においては,それらの熱伝導率は式(2)だけからは正確には決まらない.すでに,その他の超イオン伝導体(NaCl,KCl等)に対して,現象論的方程式をもとに,エネルギーの流れ J E ( t ) と電荷(ここではイオン電荷)の流れ J Z ( t ) から熱伝導率が求められている[26].それぞれ以下のように与えられる.   

J E ( t ) = L EZ T 2 ( μ Z T ) L EE ( T )
(4)   
J Z ( t ) = L ZZ T ( μ Z T ) L ZE ( T )
(5)

ここでμZは電気化学ポテンシャルである.これらの式には,それぞれの流れが互いに影響を及ぼしあう相互効果が含まれており,オンサーガーの相反定理(Lαβ = Lβα)を適用すると,最終的に熱伝導率は式(6)で表される.本研究でも相互効果が無視できないT ≥ 2500 Kの高温ではこの式を使って熱伝導率を評価した.   

κ = L EE L EZ 2 L ZZ T
(6)

ここで輸送係数Lαβはグリーン-久保の公式から,それぞれ以下のように求まる.   

L αβ = 1 V 0 d t C αβ ( t )
(7)   
C ZZ ( t ) = 1 3 k B T j Z ( t ) j Z ( 0 )
(8)
  
C EZ ( t ) = 1 3 k B T 2 j E ( t ) j Z ( 0 )
(9)
  
C EE ( t ) = 1 3 k B T 2 j E ( t ) j E ( 0 )
(10)

結局MD計算では輸送係数の定義式中の微視的な流れ j α を求めることになる. j E は式(3)と同じであり,イオン電荷の流れは次式で定義される.   

j Z ( t ) = i = 1 N z i e v i
(11)

UO2の熱伝導率の温度依存性の計算結果をFigure 7に示す.ウムクラップ過程(フォノン-フォノン相互作用)による温度上昇にともなう熱伝導率の低下が再現されている.比較のためにRonchi等によって整理された実験値とその内訳(格子振動およびスモールポーラロンの寄与)も示す[27].超高温では輻射による寄与もあるとされるが,未だによく分かっていない.MD法では格子振動による寄与しか計算されないため,1500 K付近から増加するスモールポーラロンの寄与は再現されない.更に,MD計算では実験と異なり,計算対象に不純物や転位,粒界等の格子欠陥が含まれないことから,特に低温では実験値に比べて過大評価となっている.一方,3120 K以上の液体の熱伝導率に対しては,年々実験値が低下しており,現状では測定誤差が大きいと言わざるを得ない[28,29,30,31].

Figure 7.

 Thermal conductivity of UO2. Circle symbols were given by EMD calculations. Solid line was the experimentally recommended value, which was attributed to lattice and small polaron components. Other symbols were measured values for liquid UO2.

Figure 8 (a)にUO2-ZrO2固溶体の熱伝導率の組成依存性を示す.白抜きのデータは状態図には存在しない仮想的な結晶性固体の熱伝導率である.また,固溶体と液体は簡単のため3000 K以下を固相,3500 K以上を液相とした.低温においては,UO2およびZrO2で熱伝導率は高く,お互い混合すると熱伝導率は低下する.これは,UO2-ZrO2固溶体の系で,不純物によるフォノン散乱が原因で熱伝導率が低下することが再現されていることを意味する.温度が高くなると,不純物効果の寄与は小さくなり,ウムクラップ効果の寄与が大きくなる.その結果,固溶体の熱伝導率の組成依存性は低くなる.Figure 8 (b)にUO2-ZrO2液体の熱伝導率を示す.温度依存性も組成依存性も小さくなるという結果となった.最後に,MD計算で求めた熱伝導率の誤差について述べる(Figure 7およびFigure 8参照).MD計算による熱伝導率は,輸送係数を決定する際の積分値の収束性,スーパーセルの大きさ,固溶体の陽イオン配置,計算時間(ステップ数)等に依存する誤差を持つ.今回十分な誤差評価を行ってはいないが,これまでの経験から計算誤差は熱伝導率の10%程度と予想される[32].

Figure 8.

 Thermal conductivities of UO2-ZrO2 (a) solid and (b) liquid states.

4 まとめ

UO2-ZrO2固溶体の融点を組成をパラメータとして,熱伝導率を組成・温度をパラメータとして,古典分子動力学法を用いて評価した.ポテンシャル関数には,将来の研究(ガラス固化体や溶融燃料など)を見据えて,まだ改良の余地があるものの,Pedone等によって多成分ガラス用に開発されたものを適用した.融点に関しては,2相共存法から求めた融点が実験値より300 K-600 K程高くなったものの,固溶の効果によりZrO2 50 mol %付近で最も低くなり,それが融解エンタルピーの減少と関係することを明らかにした.熱伝導率を平衡分子動力学法により,グリーン-久保の公式から算出し,イオン性固体で無視できないエネルギーと電荷の流れの相互効果を取り入れることにより,超高温領域まで格子振動の寄与を評価した.UO2固体に対して500 K-3000 Kの範囲では,おおよそ実験値を再現する結果となった.また,UO2-ZrO2固溶体に対しても,温度の上昇とともに熱伝導率が低下するウムクラップ過程および低温で顕著に表れる不純物効果による熱伝導率の低下が確認できた.一方,液体状態における熱伝導率は低く,格子振動の寄与は温度および組成依存性は小さいことが分かった.

参考文献
  • 1   K. R. Robb, M. W. Francis, L. J. Ott, Nucl. Technol., 186, 145 (2014).
  • 2   Y. Tobita, S. Kondo, H. Yamano, K. Morita, W. Maschek, P. Coste, T. Cadiou, Nucl. Technol., 153, 2245 (2006).
  • 3   K. Govers, S. Lemehov, M. Verwerft, J. Nucl. Mater., 366, 161 (2007).
  • 4   K. Govers, S. Lemehov, M. Hou, M. Verwerft, J. Nucl. Mater., 376, 66 (2008).
  • 5   T. Arima, K. Idemitsu, Y. Inagaki, Y. Tsujita, M. Kinoshita, E. Yakub, J. Nucl. Mater., 389, 149 (2009).
  • 6   P. K. Schelling, S. R. Phillpot, D. Wolf, J. Am. Ceram. Soc., 84, 1609 (2001).
  • 7   A. Pedone, G. Malavasi, M. C. Menziani, A. N. Cormack, U. Segre, J. Phys. Chem. B, 110, 11780 (2006).
  • 8   D. G. Martin, J. Nucl. Mater., 152, 94 (1988).
  • 9   K. Hirao, K. Kawamura, Material Design using Personal Computer, Shokabo, Tokyo (1994).
  • 10   A. B. Belonoshko, L. S. Dubrovinsky, Geochim. Cosmochim. Acta, 59, 1883 (1995).
  • 11   K. Harafuji, T. Tsuchiya, K. Kawamura, J. Appl. Phys., 96, 2501 (2004).
  • 12   B. Bernu, P. Vieillefosse, Phys. Rev. A, 18, 2345 (1978).
  • 13   C. Pierleoni, G. Ciccotti, J. Phys. Condens. Matter, 2, 1315 (1990).
  • 14   P. Sindzingre, M. J. Gillan, J. Phys. Condens. Matter, 2, 7033 (1990).
  • 15   P. J. D. Lindan, M. J. Gillan, J. Phys. Condens. Matter, 3, 3929 (1991).
  • 16   S. Motoyama, Y. Ichikawa, Y. Hiwatari, A. Oe, Phys. Rev. B, 60, 292 (1999).
  • 17   J.H. Harding, D.G. Martin, P.E. Potter, Thermophysical and Thermochemical Properties of Fast reactor Materials, Commission of European Communities Report EUR 12402 (1989).
  • 18   A. Navrotsky, Science, 346, 916 (2014).
  • 19   L. B. Skinner, C. J. Benmore, J. K. R. Weber, M. A. Williamson, A. Tamalonis, A. Hebden, T. Wiencek, O. L. G. Alderman, M. Guthrie, L. Leibowitz, J. B. Parise, Science, 346, 984 (2014).
  • 20   E. Yakub, C. Ronchi, D. Staicu, J. Chem. Phys., 127, 094508 (2007).
  • 21   K. A. Romberger, C. F. Baes, H. H. Stone, J. Inorg. Nucl. Chem., 29, 1619 (1967).
  • 22   T. Massalski, Binary Alloy Phase Diagrams, ASM International, Materials Park, Ohio (1990).
  • 23   K. Momma, F. Izumi, J. Appl. Cryst., 44, 1272 (2011).
  • 24   SCDAP/RELAP5–3D© Code Development Team, SCDAP/RELAP5–3D© Code Manual, INEEL/EXT-02–00589, Idaho National Engineering and Environmental Laboratory (2003).
  • 25   R. J. M. Konings, O. Benes, A. Kovacs, D. Manara, D. Sedmidubsky, L. Gorokhov, V. S. Iorish, V. Yungman, E. Shenyavskaya, E. Osina, J. Phys. Chem. Ref. Data, 43, 013101 (2014).
  • 26   N. Galamba, C. A. Nieto de Castro, J. F. Ely, J. Chem. Phys., 120, 8676 (2004).
  • 27   C. Ronchi, M. Sheindlin, M. Musella, G. J. Hyland, J. Appl. Phys., 85, 776 (1999).
  • 28   C. S. Kim, R. A. Harley, J. Fischer, M. G. Chasanov, L. Leibowitz, Proceedings of the 7th Symposium on Thermophysical Properties, 338 (1977).
  • 29   C. Otter, D. Damien, High Temp. High Press., 16, 1 (1984).
  • 30   J. K. Fink, L. Leibowitz, High Temp. High Press., 17, 17 (1985).
  • 31   M. Sheindlin, D. Staicu, C. Ronchi, L. Game-Arnaud, B. Remy, A. Degiovanni, J. Appl. Phys., 101, 093508 (2007).
  • 32   T. Arima, K. Yoshida, T. Matsumoto, Y. Inagaki, K. Idemitsu, J. Nucl. Mater., 445, 175 (2014).
 
© 2015 日本コンピュータ化学会
feedback
Top