平成２８年度

文部科学省 国家課題対応型研究開発推進事業
英知を結集した原子力科学技術・人材育成推進事業

シビアアクシデントにおける炉心構造物移行の高精度数値シミュレーション

成果報告書

平成２９年３月
国立大学法人 東京大学

本報告書は、文部科学省の英知を結集した
原子力科学技術・人材育成推進事業による委託業務として、国立大学法人

東京大学が実施した平成２６－２８年度「シビアアクシデ
ントにおける炉心構造物移行の高精度数値シ
ミュレーション」の成果を取りまとめたものです。

目次
概略・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・

iv

1. はじめに ・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・

1-1

2. 業務計画
2.1 全体計画

・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・ 2-1

2.2 各年度計画 ・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・ 2-2
3. 業務の実施内容及び成果
3.1 調査研究及び検証実験(H26-28)

・・・・・・・・・・・・・・・・・・・・・・ 3-1

3.1.1 実験装置および実験方法(H26-28)

・・・・・・・・・・・・・・・・・・ 3-1

3.1.2 共晶反応・溶融可視化実験(H26-28)

・・・・・・・・・・・・・・・・・ 3-1

3.1.3 まとめ(H26-28)

・・・・・・・・・・・・・・・・・・・・・・・・・・ 3-7

3.2 熱流動シミュレーション(H26-28)

・・・・・・・・・・・・・・・・・・・・・ 3-22

3.2.1 数値解析モデル(H26-28)

・・・・・・・・・・・・・・・・・・・・・・ 3-22

3.2.2 シミュレーション(H26-28)

・・・・・・・・・・・・・・・・・・・・・ 3-28

3.2.3 結果・考察(H26-28)
3.3 まとめ(H26-28)

・・・・・・・・・・・・・・・・・・・・・・・ 3-29

・・・・・・・・・・・・・・・・・・・・・・・・・・・・ 3-39

4. 結言 ・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・ 4-1

i

表一覧
項
表 2.1 全体計画 ・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・ 2-1
表 3.1-1 Average temperature and standard deviation ofeach test case ・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・ 3-8
表 3.1-2

Experimental conditions for specimens ・・・・・・・・・・・・・・・

3-8

表 3.1-3

Conditions of experiments・・・・・・・・・・・・・・・・・・・・・・

3-8

表 3.3-1

計算条件(Case 3)・・・・・・・・・・・・・・・・・・・・・・・・・・・ 3-30

表 3.3-2

物性値(Case 3)・・・・・・・・・・・・・・・・・・・・・・・・・・・・ 3-30

図一覧
項
図 3.1-1

Schematic image of experimental facility・・・・・・・・・・・・・・・ 3-9

図 3.1-2

Photo of experimental facility・・・・・・・・・・・・・・・・・・・・ 3-9

図 3.1-3

Design of specimen and positions of thermocouples ・・・・・・・・・・ 3-10

図 3.1-4

Time-resolved images of the reaction in Case 1 (500–4000 s) ・・・・・ 3-11

図 3.1-5

Reaction area relative to the square root of time, in which the lines show

the least square curve fit of the experiments ・・・・・・・・・・・・・・・・・・ 3-12
図 3.1-6

Average reaction layer thickness relative to the square root of time, inwhich the lines show the least square curve fit of the JAEA
results・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・ 3-12
図 3.1-7

Arrhenius plot of the reaction rate constants for B4C powder

experiments・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・ 3-13
図 3.1-8

Sudden increase of maximum layer thickness, in which the lines show the leastsquare curve fit of the experiments・・・・・・・・・・・・・・・・・・・・・・・・ 3-13
図 3.1-9

Local migration of nonreacted B4C powder (754–3940 s)・・・・・・・・・ 3-13

図 3.1-10

Designs of two specimen types and positions of thermocouples (All units aremm) ・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・ 3-14
図 3.1-11

Time-resolved images in Case 1・・・・・・・・・・・・・・・・・・・

3-14

図 3.1-12

Temperature history of specimen in Case 1 ・・・・・・・・・・・・・

3-15

図 3.1-13

Time-resolved images in Case 2 (600–1100 s) ・・・・・・・・・・・・

3-15

図 3.1-14

Temperature history of specimen in Case 2 ・・・・・・・・・・・・・

3-16

図 3.1-15

Schematic images showing different melt behaviors for thick and thin SS

plates ・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・
図 3.1-16

Schematic of the tube-type specimen ・・・・・・・・・・・・・・・・

ii

3-16
3-17

図 3.1-17(a)Results of experimental cases (Case1-2)・・・・・・・・・・・・・

3-17

図 3.1-17(b)Results of experimental cases (Case3-4)・・・・・・・・・・・・・

3-18

図 3.1-17(c)Results of experimental cases (Case5-6)・・・・・・・・・・・・・

3-18

図 3.1-18

Temperature history in Case 2 ・・・・・・・・・・・・・・・・・・・

3-19

図 3.1-19

Measured voltage applied to tungsten heaters in Case 2・・・・・・・

3-19

図 3.1-20

History of melt bottom edge in Case 2 ・・・・・・・・・・・・・・・

3-19

図 3.1-21

Temperature history in Case 6 ・・・・・・・・・・・・・・・・・・・

3-20

図 3.1-22

Measured voltage applied to tungsten heaters in Case 6・・・・・・・

3-20

図 3.1-23

Diagnostics in Case 6 ・・・・・・・・・・・・・・・・・・・・・・・

3-20

図 3.3-1

MPS 法における温度-エンタルピーモデル ・・・・・・・・・・・・・・・

3-31

図 3.3-2

MPS 法における相状態-エンタルピーモデル ・・・・・・・・・・・・・・

3-32

図 3.3-3

CSF モデル概略図・・・・・・・・・・・・・・・・・・・・・・・・・・

3-33

図 3.3-4

CSF モデルにおける仮想粒子・・・・・・・・・・・・・・・・・・・・・

3-34

図 3.3-5

粒子初期配置(Case 3)・・・・・・・・・・・・・・・・・・・・・・・・

3-35

図 3.3-6

シミュレーション結果図(Case 3)・・・・・・・・・・・・・・・・・・・

3-36

iii

概略
福島第一原子力発電所の原子炉では、東日本大震災によりシビアアクシデントが発生した。シビ
アアクシデントにおいて、原子炉の冷却機能が喪失して崩壊熱を十分に除去できなくなると、燃料棒、制御棒および周辺の炉心構造物の温度が上昇する。このとき、異種の金属が接触していると、共晶反応が起こるため、個々の材料はその融点よりも低い温度で溶融する。なかでも、制御棒においてボロンカーバイドとステンレス鋼は、共晶反応により、融解が 1200℃程度で始まり、シビアア
クシデントの早期に融解が始まることが知られている。この溶融物は、炉心下部で冷やされ再凝固するため、一部は原子炉圧力容器内にとどまるものの、原子力圧力容器の外部へ落下・流出するものもある。従って、共晶反応を伴う制御棒の溶融は、シビアアクシデントの進展を把握するために極めて重要である。他方、既存のシビアアクシデント解析コード（MELCOR、RELAP、など）では、制御棒に関する共晶反応を精度よく模擬することができなかった。共晶反応を伴う炉内構造物の再配置を予測するために、数値流体力学の応用が期待されている。シビアアクシデントでは、融解に伴う構造物の大変形を模擬する必要があるため、Moving Particle Semi-implicit Method（以下、MPS
法と記す）のようなラグランジュ的手法による数値流体力学の応用することにより、現象を精度よく模擬できると期待されている。このような背景から、本プロジェクトの目的は、シビアアクシデ
ント解析コードの高度化に資するため、MPS 法に共晶反応モデルを導入して、高精度熱流動解析コードを開発することである。本研究では、ボロンカーバイドとステンレス鋼の共晶反応に関する実験を行い、シミュレーションに必要な境界条件などのデータベースを作成する。数値解析においては、MPS 法を用いて、制御棒における共晶反応とこれに伴う溶融・再配置の進展を模擬できるような新しいモデルを開発する。英国側の共同研究先の Imperial College London の研究者と協力して数値解析モデリングの検討および結果の妥当性検証を行う。

iv

1.

はじめに

適切にシビアアクシデントマネージメントが実施できるように、様々な条件下での炉心損傷を模擬することは極めて重要である。既存の数値シミュレーション技術は、熱水力学的にシビアアクシ
デント事象全体を対象とするものが主であり、例えば、制御棒の溶融過程については十分に模擬できなかった。そこで、本研究では、制御棒の溶融過程および再配置を模擬するための数値解析技術を開発するとともに、検証実験を行う。本研究に従事する日本側の研究者および英国側の研究者
は、それぞれ、制御棒の溶融過程および再配置を模擬するために要求される高度な数値解析技術を有しているため、協力しながら研究を遂行する。

1-1

2.

業務計画

2.1 全体計画
本研究において、日本側の体制はシビアアクシデントおよび流体実験の専門家と計算力学の専門家で構成され、英国側はシビアアクシデントの専門家と計算力学の専門家で構成される。
本研究は、「調査研究及び検証実験」と「熱流動シミュレーション」で構成される。各研究項目は研究代表者と研究分担者が取りまとめる。研究の実務は、それぞれの研究者が指導する大学院生または研究員の若手研究者が担当する。表 2.1 に示すように、本研究は３年で構成されている。
表 2.1 全体計画
年度
平成 26 年度

平成 27 年度

平成 28 年度

項目
フェーズ

要素研究フェーズ
（特性解明、要素開発、理論検討）

検証フェーズ（制御効果、適用性評価）(1)調査研究及び検証実験①調査と予備実験

予備検討
実験実施（数値解析条件への反映を確認）②実験実施とﾌｨ—ﾄﾞﾊﾞｯｸ

実験実施③実験データによる評価(2)熱流動ｼﾐｭﾚ—ｼｮﾝ

数値解析・実験意見交換①基盤技術開発
数値解析・実験意見交換

シンプルな体系における②検証

数値解析の実行東大と IC の同一条件での（Verification）

数値解析の実行と結果比較③妥当性確認

数値解析の実行（実験結果を参考）（Validation）(3)研究推進①まとめ

まとめ・評価まとめ・評価②まとめ（２年目）③まとめ（３年目）

まとめ・評価初年度においては、炭化ホウ素とステンレス鋼の共晶反応に関する調査及び検証実験を行うとともに、MPS 法を用いて溶融のようなシビアアクシデントでは基本となる伝熱流動現象を模擬できることを示す。２年目に、平成 26 年度に製作した共晶反応実験装置を用いて可視化手法に関する実験を行
い、実験データを取得する。数値シミュレーションに関しては、平成 26 年度に開発した粒子法を使用した熱流動シミュレーションコードの妥当性を検証する。３年目は、共晶反応実験装置を用いて、2-1

共晶反応について検討を行う。また、MPS 法を用いて、これまでの実験成果を参考にして共晶反応を考慮した固体溶融に関する数値シミュレーションを実行する。数値シミュレーションにより、制御棒の溶融に係わる実験で観察される典型的な現象のキャンドリングを模擬できることを確認する。

2.2 各年度計画
本業務では、シビアアクシデント時の事故進展を模擬するため、東京大学及び英国 Imperial
College London（以下、「IC」という。）が協力して、制御棒の共晶過程及び再配置を模擬するための詳細な熱流動シミュレーション技術をそれぞれ開発して、検証（Verification）と妥当性確認（Validation）に関する研究を行った。
(1) 調査研究及び検証実験
平成 26 年度
①調査と予備実験
シビアアクシデントの熱流動シミュレーションで使用する、ホウ素と鉄の共晶状態に関する相図を調査する。また、シミュレーションの検証データを取得するための共晶反応実験装置を製作し、可視化手法に関する予備実験を行う。
平成 27 年度
②実験実施とフィードバック
熱流動シミュレーションコードの妥当性確認（Validation）のため、実験条件が適切に数値解析に反映できることを検討する。前年度に製作した共晶反応実験装置を用いて、可視化手法に関する実験を行い、実験データを取得する。
平成 28 年度
③実験データによる評価
熱流動シミュレーションコードの妥当性確認（Validation）のため、平成 26 年度に製作した共晶反応実験装置を用いて、実験を行い、実験データを取得するとともに、共晶反応について検討を行う。
(2) 熱流動シミュレーション
平成 26 年度
①基盤技術開発
MPS 法を使用して熱流動シミュレーションのコードを開発する。基礎的な数値実験として、シンプ
ルな体系において高粘性流体の数値シミュレーションを実行する。IC と協力して、数値シミュレー
ション結果の妥当性について考察する。
平成 27 年度
②検証（Verification）

2-2

平成 26 年度に開発した粒子法を使用した熱流動シミュレーションコードについて、妥当性を検証する。検証実験を参考にして、計算体系および計算条件を選定する。IC と協力して、数値シミュレ
ーション結果の妥当性について考察するとともに、手法間の特徴を整理する。
平成 28 年度
③妥当性確認（Validation）
これまでに開発した MPS 法の数値解析技術を用いて、これまでの実験成果を参考にして共晶反応を考慮した固体溶融に関する数値シミュレーションを実行する。数値シミュレーションにより、制御棒の溶融に係わる実験で観察される典型的な現象のキャンドリングを模擬できることを確認する。
(3) 研究推進
平成 26 年度
研究代表者の下で各研究項目間における連携を密にして研究を進める。
平成 27 年度
研究代表者の下で各研究項目間における連携を密にして研究を進める。研究成果は、国際会議などで発表する。
平成 28 年度
研究代表者の下で各研究項目間における連携を密にして研究を進める。研究成果は、国際会議などで発表する。

2-3

3.

業務の実施内容及び成果

3.1 調査研究及び検証実験(H26-H28)
炉心構造物移行の高精度シミュレーションの V&V を目的として、BWR 制御棒材料における溶融進展と溶融物落下挙動データベースを構築している。実験では、大きく分けて 2 種類の試験体を用いた可視化実験を行った。1 つは板状試験体を用いた実験、もう 1 つは実際の制御棒と同様な形状であるチューブ状試験体を用いた実験である。実験で取得した画像データを解析することで、定性的・定量的なデータを得た。

3.1.1 実験装置および実験方法(H26-H28)
実験装置は、試験体をタングステンヒーターからの輻射伝熱で加熱する構造となっている。実験中は、共晶反応の進展や溶融物の挙動をカメラで撮影すると同時に、熱電対とデータロガーを用いて高精度に試験体の温度を記録した。図 3.1-1 に実験装置の構成を示す。
実験装置は、異種材料が互いに接することを避ける構造とした。超高温では金属材料が大きく熱膨張し、異種材料は異なる熱膨張率を持つため、接合面が存在すると接合面で大きな応力が生じて実験装置が破断する。したがって、異種材料が接する面積を極力小さくして熱応力が逃げやすい構造とした。例えば、試験体とタングステンヒーターは接面をなくし、試験体は電極部分から吊り下げるようにして保持した。実験中は、タングステンヒーターが 2300–3300 K 程度の高温になるため、ヒーター両端を固定支持にしてしまうと熱膨張によって座屈してしまう。そこで、下端をローラー支持とし、おもりでヒーターを下方向へ引張ることでヒーターが熱膨張で歪むことがないようにした（図 3.1-1：Traction system）。試験体に電流が流れこむと試験体が直接ジュール加熱されてしまい、溶融によるジオメトリ変化に伴って発熱位置が変わってしまう。このような境界条件はシミュレーションには適さないため、試験体を電気的に絶縁する必要がある。そこで、試験体と電極部分の接面にセラミック板を挟みこむことで、試験体-電極間の絶縁性と断熱性を確保した。
1473 K 以上の超高温環境下では、酸化物が非常に生じやすい。そこで、実験装置に Ar ガスを実験装置下部から流入することで酸化を抑制した。装置内酸素濃度を 0.1%以下にした後に、実験を開始した。試験体の温度が所望の温度に達した後は、入力電力を適切に設定して温度を制御した。図 3.1-2
は、図 3.1-1 に示した実験装置の写真である。

3.1.2 共晶反応・溶融可視化実験(H26-H28)(1) 板状試験体を用いた可視化実験
実験では、粉末のボロンカーバイド（以下、B4C と呼称）をステンレス鋼 （SUS304、以下 SS と呼称）ではさんだ試験体を用いた。粉末 B4C の平均粒径は 20–30 μm である。図 3.1-3 に試験体とタングステンの概略図を示す。左図は正面から見た試験体、右図は試験体を側面から見た図である。
図中には熱電対の挿入位置もしめしている。熱電対を挿入するための穴は，正面図における SS の横
方向の中心に穴あけし、熱電対挿入深さは奥行きの半分 （2.5 mm）とした。また、熱電対には K タ
イプを用いた。

3-1

(2) 結果と考察
画像解析によって得られた定性的・定量的なデータを報告する。表 3.1-1 に実験ケースの条件を示す。σは測定温度の標準偏差を表している。Average temp.は 4 つの熱電対の値を平均したものである。ただし、所望の温度に到達した後の温度データを平均化している。それぞれのケースにおいて試験体の左右の SS 部分で測定された温度が異なったため、解析は左右別々に行った。L と R はそれ
ぞれ Left と Right を表す。また、実験中に温度の大きな変化が確認された Case 1 R については、
実験開始から温度変化が確認された時間までと温度変化後から実験終了時間までの 2 つの期間に分けて解析を行った。
図 3.1-4 に Case 1 の可視化画像を示す。本実験の設定温度は SS が溶融する温度ではないため、B4Cと SS の境界部分から生じている溶融物は全て共晶溶融物である。試験体の上下端には少量の酸化タ
ングステンが見られるが、鮮明な可視化画像を得ることができている。2000 s には溶融物が B4C 領域に侵入する様子が見られた。これと同様な現象を示唆する結果が、過去の金属組織学的アプロー
チによる実験でも得られている[1]。本研究によって、B4C 領域へ溶融物が侵入していく過程がリアルタイムに確認された。これは、B4C が SS に融けていくと同時に、粉末に存在した小さな空隙がつながって、大きな空隙が生じることに起因する現象であると考えられる。3500 s には、未反応 B4C
粉末の塊が溶融物内を未反応 SS 側へ移流していく様子が見られた。反応領域では活発な対流が見られる低粘性部分と移流が確認されない高粘性部分が共存していた。これは、共晶反応中は非定常状態であり、溶融物の組成が均一ではなかったことに起因すると考えられる。また、実験中に溶融物が試験体から流れ出ることはなかった。溶融物が試験体下部の未反応部分に支えられたこと、溶融物自体が非常に大きな表面張力を持っていたこと、固相と液相が共存していたことが原因だと考えられる。
図 3.1-5 は共晶反応面積が時間経過と共に拡大していく様子を示すグラフである。図中に示された
線は、最小二乗法によって作成された近似曲線である。横軸は実験時間の平方根、縦軸は共晶反応
面積をそれぞれしめしている。反応面積の増加と実験時間の平方根は線形な関係を持つことが確認された。B4C と SS の共晶反応による反応面積の増加が放物線律速によって表現できるということを示している。
図 3.1-5 の時間分解測定データを、過去に Nagase ら[1]によって取得されたデータと比較したもの
が図 3.1-6 である。時間分解データについて Nagase らは B4C ペレットを用いた実験結果のみを報告しているため、図 3.1-6 の参考データは全てペレットについてのものである。Nagase らは加熱実験
後の試験体を観察して、共晶反応部分の平均反応層厚さを測定した。したがって、ここでは平均反応層厚さについて本実験結果と比較した。図中に示された直線は Nagase らが最小二乗法で作成した
線形近似曲線である。
図 3.1-6 を見てわかるように、本研究で得られたほぼ全てのデータは、Nagase らの 1473 K と 1498
K の線形近似曲線の間にプロットされている。本研究の試験体温度である 1473 K 程度では、Nagase
らの実験結果よりも本実験の方が厚い反応層を形成している。また、Nagase らの手法は加熱実験後に試験体を観察するという静的アプローチであるため、実験結果の時間分解能が悪い。一方、本研
究ではリアルタイム可視化という動的アプローチを採用したため、良好な時間分解能で定量的なデータ得ることができている。

3-2

図 3.1-7 は共晶反応のアレーニウスプロットである。横軸は温度の逆数、縦軸は共晶反応による反
応層成長率である。図中に示された直線は、それぞれ Nagase ら[1]、Hofmann ら[2]および Shibata
ら[3]が報告したアレーニウスプロットの近似曲線であり、全てのプロットが粉末 B4C を用いた実験の結果である。本研究でも粉末の B4C を用いているため、粉末 B4C のアレーニウスプロットを比較することで本研究手法の妥当性を検証することができる。本研究と既存研究の実験温度が同じ領域では反応層成長率の差異はほとんどなく、その差は 1 オーダー以下である。1 オーダー程度の反応層成
長率の差異は Shibata ら[3] によっても報告されており、主な原因は粉末の充填率の差異であると考えられている。粉末の充填率が低いと表面エネルギーは高くなるため、充填率が比較的低い今回の実験では液相出現温度が低下する。既存研究に比べて低温で生じる液相がホウ素の移流拡散を促進したために本研究は既存研究よりも少し大きな反応層成長率を示したと考えられる。以上より、
本研究手法は妥当であると判断できる。既存研究では 1473 K 付近で不連続な部分が存在し、この温度付近におけるアレーニウスプロットが欠落している。この不連続部分は、B4C と SS が共晶反応を開始し、液相が出現する温度帯である。したがって、反応層厚さの測定では高い時間分解能が要求され、過去の静的なアプローチでは所望のデータを得ることは困難であった。本研究の動的可視化
手法は、既存研究の静的アプローチでは確認することができなかった領域もカバーすることができている。
図 3.1-8 は Case 1 について最大反応層厚さをプロットしたものである。Case 1 R は Case 1 L に比べて線形近似曲線の傾きが急であり、図 3.1-5 と同様な傾向を示している。また、ある時刻までは両結果とも図 3.1-5 と同様に放物線律速に従っている。しかし、Case 1 L では 3600 s 付近で、Case
1 R では 4225 s 付近で急激な最大反応厚さの増加が見られた。このような急激な変化は図 3.1-5 では見られなかったことから B4C の移流に関係した局所的な現象が原因になっていると考えられる。
最大反応厚さの急激な増加が見られた時間帯に生じた現象は局所的な未反応 B4C 粉末の移流であった（図 3.1-4）。図 3.1-9 は未反応 B4C 粉末を詳しく確認するために図 3.1-4 の一部分に着目した拡大画像である。図 3.1-9 では移流を確認しやすくするために 1 マスあたり 0.5 × 0.5 mm のグリッドを加えた。図 3.1-9 を見ると、未反応 B4C 粉末が試験体の内側から外側へ向かって融液内を流れていることがわかる。未反応 B4C 粉末の移流メカニズムとして融液における自然対流やマランゴニ
効果が考えられる。しかし、本現象には局所的な温度分布と組成分布が関係していると考えられ、本画像解析のみから未反応 B4C 粉末の移流メカニズムを特定することは難しいと考えられる。

(3)(ア) 板状試験体を変更した実験（追加実験）

上述の実験では、共晶溶融物の落下や流動といったリロケーション挙動を確認することはできなかった。板状試験体を用いた実験の考察から、試験体の垂直方向温度分布、溶融物の表面積と体積の比率および試験体の SS 部分の厚さが溶融物の挙動に影響するのではないかという仮説を立てた。本
仮説を検証すべく下記の追加実験を行った。
図 3.1-10 は試験体の概略図である。本実験では、これまでの実験で使用した板状試験体と同じ
Front-uncovered specimen と板状試験体の正面に SS 薄板を取り付けた Front-covered specimen の2 種類を用いて実験を行った。粉末 B4C の充填率は約 45%であり、空隙は空気で満たされている。試験体温度は K タイプ熱電対で直接測定した。表 3.1-2 は実験条件である。Case 1 では、Front-

3-3uncovered specimen を上下逆さまに設置し、上端が最高温度になり下方へ行くに従って温度が低くなるような温度勾配を生じさせた。試験体温度が 1473K に達した後は、最高温度 1523 K まで 0.08
K/s で昇温した。Case 2 では、「 (ア) 板状試験体を用いた可視化実験」と同様に試験体を配置した。また、実験中は試験体温度を 1493 K で維持した。Case 2 が「(ア) 板状試験体を用いた可視化実験」と異なる点は、正面を SS 薄板で覆っている Front-covered specimen を用いている点と試験体の長さのみである。

(4) 結果と考察
図 3.1-11 は Case 1 の可視化結果、図 3.1-12 は温度の測定結果である。図 3.1-11 では最初に試験
体上部で共晶反応が生じていることがわかる。これは温度が最も高い試験体の上端で B4C が SS へ最も速く拡散したことによると考えられる。490 s には共晶溶融物が B4C 領域に侵入しており(図 3.111 内(1))、1090 s にはさらに多くの共晶溶融物が粉末 B4C に吸収され混合している(図 3.1-11 内(2))。690 s から確認される白く光っている部分はハレーションである(図 3.1-11 内(3))。これは共晶反応によって生じた溶融物が表面張力によって丸みを帯び、タングステンからより多くの光をカ
メラ側へ反射するようになったためだと考えられる。690 s までは共晶反応が SS 部分の表面まで到達していない。690 s からは、共晶反応が SS 部分の表面まで到達して、共晶溶融物が下部へ移行していく挙動を観察することができた(図 3.1-11 内(4))。これは、リロケーション挙動の一端を観察することができたということである。「 (ア) 板状試験体を用いた可視化実験」との主な違いは試
験体に生じた温度勾配である。「 (ア) 板状試験体を用いた可視化実験」では、高温の中心部から低温の両上下端へと温度勾配が生じていたが、本追加実験では、高温の上部から低温の下部への温度勾配が生じた。つまり、本追加実験と同様な温度勾配が存在する時の方が流動を伴うリロケーシ
ョンが生じやすいと考えるられる。これは、上部から溶融が開始すると溶融物の大部分が雰囲気にさらされて表面積が大きくなり、未反応の固体 SS に囲まれている場合よりも流れやすくなるからだと考えられる。そのため、共晶反応の進行方向に未反応の固体 SS が存在する間は溶融物が流れ落ちなかったと考えられる。
図 3.1-13 は Case 2 の可視化結果であり、図 3.1-14 は測定された温度の記録である。正面の SS 薄
板右に取り付けられた熱電対（Surface/Right）は途中で試験体から外れてしまったため、他の熱電対と異なる挙動を示しており、信頼性が低いと考えられる。図 3.1-13 では、最初に試験体の中心部から溶融が生じていることがわかる。これは試験体に中心部から両端へ向かって温度が低くなるよ
うな垂直方向の温度勾配が生じていたためである。さらにタングステンヒーターで試験体の両端から輻射加熱しているため水平方向の温度勾配は中心部から両端へ向かって温度が高くなるような勾配であるため、試験体の中心部でも特に SS 薄板の B4C と接している範囲内の最も外側から溶融が開始している(図 3.1-13、600 s、(1))。800 s には試験体内の B4C が Ar 雰囲気へ露出され始める(図3.1-13、(3))。この時、共晶反応によって溶融した SS 薄板は中心部から試験体の右端へと移動する挙動を示した。900 s には右端へ移動した溶融物は試験体に沿って下方へ移動し液滴を形成した(図3.1-13、(4))。900-1100 s にかけて、この液滴に溶融物が流れていくことで液滴の大きさは大きくなっていった。Case 2 では Case 1 で確認された溶融物の B4C 領域への侵入・吸収は確認されなかった。むしろ、Case 2 では溶融物は粉末 B4C にはじかれるような挙動を示した。

3-4

図 3.1-15 は異なる SS 板厚が持つ溶融物挙動への影響を定性的に説明した模式図である。左右には板厚が厚い場合と薄い場合がそれぞれ示されている。板厚が厚い場合（本実験では 1 mm）では、共晶反応開始後しばらくは溶融物表面が固相 SS によって囲まれるため、溶融物が試験体内部に比較的長い時間保持される。この間、溶融物は重力によって試験体下方へ流れないため、溶融物が粉末 B4Cの空隙へ侵入していく時間的余裕が生じる。その結果、溶融物は粉末 B4C と混合しながら溶け落ちたと考えられる。一方、板厚が薄い場合（本実験では 0.1 mm）では、溶融物界面の多くは気液界面となる。固液界面の面積が小さいため、共晶反応後に溶融物は容易に移動することが可能となり、溶融物の流れが生じて液滴形成につながったと考えられる。溶融物の粘性は非常に小さいため、溶
融物形成後は重力によって速やかに下方へ流れたと考えられる。その結果、溶融物は粉末 B4C に吸収される前に下方へ移動していったと考えられる。すなわち、Case 1 と Case 2 の溶融物挙動の違い
は溶融物体積と固液界面（試験体-溶融物）のバランスに由来すると考えることができる。
本実験で試験体の垂直方向温度勾配と SS 厚さが流れを伴うリロケーションを引き起こす重要なパラメータとなり得ることが示唆された。上述の実験結果およびその考察から得られた知見を元に以下にしめす「(イ)チューブ状試験体を用いた可視化実験」を行った。

(5) (イ)チューブ状試験体を用いた可視化実験
図 3.1-16 はチューブ状試験体の概略図である。本実験では、試験体を吊り下げるのではなく、下部
電極部の上に固定した。一部の実験ケースでは試験体下部にセラミックホルダーを設置することで試験体下部からの熱逃げを低減した。これにより、試験体を約 1573 K 超の高温まで加熱することが可能となった。実験ではセラミックホルダーを使用した場合と SS チューブのみを使用した場合の両方を行った。

(6) 実験結果
表 3.1-3 に各ケースの実験条件を示した。D は SS チューブの外径、d は SS チューブの内径、Length
は SS チューブの高さをそれぞれ表している。Packing ratio (P. R.)は粉末 B4C の充填率を表している。B4C 粉末の空隙は空気で満たされている。SS: B4C は試験体の SS と B4C の質量比を表している。Top と Mid.は、それぞれ試験体上端からの熱電対位置を表している。Temp.は Top 位置における試験体の平均温度である。
図 3.1-17 は実験結果画像である。Case 1 では、共晶溶融物が未反応の B4C にまとわりついて膜を形成した。このとき、SS チューブ内の多くの B4C が未反応のまま残っており、試験体は崩落することなく元の形状を維持した。溶融によって上方の一部の SS が下方へ移動することで、チューブ内に充填した B4C が Ar 雰囲気へ露出した。Case 2 では SS チューブの下端にセラミックホルダーを設置して、Case 1 よりもさらに高い温度で実験を行ったが、試験体表面に形成された膜は流れ落ちなかった。また、セラミックホルダーを設置したことによって SS チューブの温度勾配が Case 1 よりも緩やかになったが、その影響は確認されなかった。Case 3 では、Case 1 と 2 と同じ断面寸法の SS チ
ューブを使用し、SS チューブ内にセメントをコーティングしたφ3.0 mm の棒を挿入した。SS チュー
ブ内に棒を挿入して、充填できる B4C の量を減らすことで SS: B4C を調整した。溶融物は試験体の真3-5

中部分で上下方向に分裂して、それぞれが縮む様にして移動していった。Case 4 では、Case1 と 2と同様な溶融物の膜を形成すると同時に溶融物の液滴が形成された。溶融物の膜は下方へ流れなかったが、液滴は時間の経過とともに下方へ流れていった。Case 5 では溶融物の膜は形成されず溶融物は全て液滴を形成した。試験体の形状を維持した Case 1–4 と異なり試験体は崩落していった。
Case 6 では SS チューブの下端にセラミックホルダーを設置したが、セラミックホルダーを設置したことに起因する現象の変化は観察されなかった。また、Case 5 と同様に溶融物は全て液滴を形成した。
「(イ) チューブ状試験体を用いた可視化実験」では大きく分けて溶融物の液滴形成と膜形成という2 種の現象が観察された。試験体の組成比が Fe-B の共晶点から離れており、SS 部分の断面積が小さいほど膜を形成する傾向があった。逆に、試験体の組成比が Fe-B の共晶点に近く、SS 部分の断面積が大きいほど液滴を形成する傾向があった。試験体の組成比が Fe-B の共晶点から離れているとき、試験体の B4C 領域と SS 領域の境界付近にある少量の B4C のみでほぼ全ての SS が溶融する。したがって、多くの未反応 B4C が試験体に残留する。多くの残留した未反応 B4C は溶融物に対してメソ多孔質体のように振る舞い、Case 1 や Case 2 のように溶融物を表面に吸着すると考えられる。溶融物が表面に吸着した後、等温凝固によって溶融物の粘性が上がっていった可能性が考えられる。一方で、試験体の組成比が Fe-B の共晶点に近いとき、SS との共晶反応によってほとんどの B4C が SS に融け込む形で消費される。また、共晶溶融物の組成比が共晶点に近い値となって粘性が小さくなり、Case 3 と Case 5 および Case 6 のように溶融物が流れやすくなると考えられる。SS 部分の断面積が大きいと、共晶反応によって生じる溶融物の体積が大きくなる。溶融物の体積が大きくなると溶融物に生じる重力が表面力に勝り、Case 5 と Case 6 のように溶融物は液滴を形成しながら試験体下方へ流れたと考えられる。
本実験で得られた知見を元にフクシマ事故における制御棒溶融過程を考察する。福島第一原子力発
電所の BWR 制御ブレードにおける組成比 (SS: B4C)と SS 部分の断面積のみを考慮すると、制御ブレードが SS の融点に達するまでの間は Case 2 のように膜と液滴がどちらも形成されていたと考えられる。しかし、本実験と実際の過酷事故では境界条件が異なるため、慎重な議論が必要である。例
えば、本実験とフクシマ事故の大きな違いのひとつとして環境雰囲気の違いが挙げられる。本実験ではアルゴンガス雰囲気である一方で、フクシマ事故では大量の水蒸気が炉内に存在する酸化雰囲
気であったと考えられている。酸化雰囲気下では、共晶溶融物表面に酸化ホウ素 (B2O3)が生成する
ことが知られている[4]。酸素分圧が高い環境で生成した B2O3 は対流混合を抑制し、融液内の組成偏析を引き起こす [5]。その結果、融液の見かけの粘性が高くなると考えられ、フクシマ事故における制御ブレード溶融移行は本研究で観察された溶融移行よりも進行が緩慢であったと予想できる。

(7) Validation data
上述のチューブ状試験体の実験 Case 1–6 の内、膜を形成した Case 2 と液滴を形成した Case 6 について時間分析 Validation data を作成した。はじめに Case 2 の Validation data を説明する。図3.1-18 は測定した温度データである。Top と Mid.はそれぞれ Top 位置と Mid.位置における試験体温度、Gas は出口ガス温度、Window は実験装置の可視化窓温度、Wall は実験装置の外壁温度である。
図 3.1-19 はタングステンヒータに印加されていた電圧の測定データである。図 3.1-19 の青色部分

3-6で電圧と電流ともに一定にしたため、1000 s よりも前のデータについて数値計算結果と比較するのが適切である。また、この電圧と電流が一定である間のタングステンヒータ温度を電気抵抗率から推定した[6]。その結果、タングステンヒータ温度は約 2028 K であると推定することができた。図3.1-20 は共晶溶融物の下端位置のプロットである。縦軸の”Relative position from top”とは試
験体の上端から溶融物下端位置までの鉛直方向の距離を表している。
次に Case 6 の Validation data について説明する。図 3.1-21 は測定した温度データ、図 3.1-22 はタングステンヒータに印加されていた電圧の測定データ、図 3.1-23 は数値計算と比較するための診断量である。Case 2 のときと同様にタングステンヒータの温度を推定したところ、約 2421 K であった[6]。図 3.1-23 の左図は溶融物で形成された液滴の横幅、右図は試験体を正面から見たときの溶融部分の面積をそれぞれ表している。
以上の Validation data を用いて数値計算結果を検討することで数値計算モデルの妥当性を定性的・定量的に評価することが可能である。

3.1.3 まとめ(H26-H28)


BWR 制御棒材料を対象とした超高温環境下の共晶反応可視化手法を開発した。



本研究で開発した共晶反応可視化手法を用いて B4C /SS の溶融進展と溶融物落下挙動を可視化した。



溶融進展については、主に以下の２つが分かった。

 未反応 B4C 粉末が融液内を移流していくことで、局所的に溶融進展が加速されることが分かった。
 溶融進展の最中は、反応領域において固相と液相が共存することが分かった。


溶融物落下挙動については、主に以下の３つが分かった。

 落下挙動には、いくつかのモードが存在することが明らかになった。
 落下挙動のモードは大きく分けて、膜を形成して元の形状を維持するモードと液滴を形成して大きく崩落するモードが確認された。
 本研究の実験条件では、上記モードの決定因子は SS 部分の断面積と B4C /SS の組成比であった。
以上の実験について、シミュレーションで必要となる実験中の熱的な境界条件が得られた。

3-7

表 3.1-1 Average temperature and standard deviation of each test case

Time
(s)

Average temp.
(K)σWidth of B4C region(mm)Case 1 L

700–3000

1477

16.21

4.0

Case 1 R

700–1500

1458

2.51

4.0

1500–4500

1472

6.11

Case 2 L

250–1500

1455

19.37

2.0

Case 2 R

250–1500

1464

15.52

2.0

表 3.1-2

Name
Case 1Case 2

Experimental conditions for specimens

Front

Specimencondition

direction

Uncovered

Upside down

CoveredNormalTemperature [K]Height [mm]1523 *

35

**

65

1493

*

: Maximum temperature, **: Constant temperature

表 3.1-3
Case(D, d) [mm]Length [mm]Conditions of experiments
P. R.

SS:B4C

Top, Mid.

Temp. [K][%]1

4.2, 3.64

6550

68.2:31.8

1.0, 35.0

1498

2*

4.20, 3.64

70 (25)**38

73.2:26.8

0.0, 10.0

1544

3

4.2, 3.64*

65

38

90.8:9.2

1.0, 35.0

1530

45.0, 4.0704280.1:19.9

2.0, 40.0

1473

56.0, 4.0

7055

85.8:14.2

1.0, 35.01540

65.0, 3.0

65 (35)41

92.8:7.2

0.0, 25.0

1486* Case 2 では SS チューブ寸法を変えないで SS:B4C を変化させるために、セメントでコーティング
したφ3.0 mm のチューブを中心に挿入した (図 3.1-16、Double tube)。
** Length の（）内の数字は SS 部分の長さである。
** 実験では K タイプ熱電対を用いているため、これらの温度は参考温度である。

3-8

(a)側断面図(b) 正断面図

図 3.1-1 Schematic image of experimental facility

図 3.1-2 Photo of experimental facility

3-9図 3.1-3 Design of specimen and positions of thermocouples (Units are [mm])3-10図 3.1-4 Time-resolved images of the reaction in Case 1 (500–4000 s)

3-11図 3.1-5 Reaction area relative to the square root of time, in which the lines show
the least square curve fit of the experiments

1.6
Case 1 L

Average rection layer thickness [mm]1.4Case 1 R
Case 2 L

1.2

Case 2 R

1★ Nagase et al. (pellet)1473 K

Nagase et al. (pellet)1498 K

0.8
0.6
0.4
0.2
0
0

20

40

60
80
Square root of time [s1/2]

100

120

図 3.1-6 Average reaction layer thickness relative to the square root of time, inwhich the lines show the least square curve fit of the JAEA results

3-12図 3.1-7 Arrhenius plot of the reaction rate constants for B4C powder experiments

図 3.1-8 Sudden increase of maximum layer thickness, in which the lines show theleast square curve fit of the experiments

図 3.1-9 Local migration of nonreacted B4C powder (754–3940 s)

3-13

(a) Front-uncovered specimen (b) Front-covered specimen (c) Front view of both
specimen types
図 3.1-10 Designs of two specimen types and positions of thermocouples (All units
are mm)

図 3.1-11 Time-resolved images in Case 13-14

図 3.1-12 Temperature history of specimen in Case 1

図 3.1-13 Time-resolved images in Case 2 (600–1100 s)

3-15図 3.1-14 Temperature history of specimen in Case 2

図 3.1-15 Schematic images showing different melt behaviors for thick and thin SS
plates3-16図 3.1-16 Schematic of the tube-type specimen

Case 1 (0–1487 s)

Case 2 (0–1554 s)

図 3.1-17(a) Results of experimental cases (Case1-2)

3-17Case 3 (0–505 s)

Case 4 (0–1100 s)

図 3.1-17(b) Results of experimental cases (Case3-4)

Case 5 (0–1250 s)

Case 6 (0–1863 s)

図 3.1-17(c) Results of experimental cases (Case5-6)

3-18図 3.1-18 Temperature history in Case 2

図 3.1-19 Measured voltage applied to tungsten heaters in Case 2

Relative position form top [mm]

20
18

Right

16

Left

14

Middle

12
10
8
6
4
2
0
0

200

400

600

800
1000
Time [s]

1200

1400

1600

図 3.1-20 History of melt bottom edge in Case 2

3-19

Temperature [K]

1600
1400
1200
1000
800
600
400
200
0

Top
Mid.
Wall

0

300

600

900
1200
Time [s]

1500

1800

図 3.1-21 Temperature history in Case 6

図 3.1-22 Measured voltage applied to tungsten heaters in Case 6

Increase of droplet width
図 3.1-23 Diagnostics in Case 6

3-20

Development of reaction area

参考文献[1] F. Nagase, H. Uetsuka and T. Otomo, 1997. Chemical Interactions between B4C andstainless steel at high temperatures. Journal of Nuclear Materials, 245(1), pp. 52–
59.
[2] P. Hofmann, M.E. Markiewicz and J.L. Spino, 1994. Reaction behavior of B4C absorbermaterial with stainless steel and zircaloy in severe light water reactor accidents.
Nuclear Technology, 90(2), pp. 226–244.
[3] H. Shibata, K. Sakamoto, A. Ouchi and M. Kurata, 2015. Chemical interaction between
granular B4C and 304L-type stainless steel materials used in BWRs in Japan. J. Nucl.
Sci. Technol., 52(10), pp. 1313–1317.
[4] A. Erdemir, C. Bindal, G.R. Fenske, 1996. Formation of ultralow friction surface
films on boron carbide. Applied Physics Letters, 68, pp. 1637.
[5] M. Furuya and S. Morooka, 2015. Visual Observation of Natural Convection due toEutectic Melting of Stainless Steel and Boron Carbide and Three-Dimensional Raman
Spectrometry. 第 52 回日本伝熱シンポジウム講演論文集.
[6] H. T., Evans Jr., 1994. Handbook of Chemistry and Physics 75th Ed. Boston: CRC Press.

3-21

3.2 熱流動シミュレーション(H26-H28)
3.2.1 数値解析モデル(H26-H28)(1) MPS 法①

支配方程式

流体の支配方程式は、微圧縮を仮定した連続の式および Navier-Stokes 方程式を用いる。

⋅

(3.1)

0

1

(3.2)

⋅

ここで、、t、u、P、および G はそれぞれ、密度、時間、流速、圧力、粘性応力テンソルおよび外力による加速度である。

②

微分演算子の離散化モデル

ｃおよびラプラシアンは次式のように離散化される。
Φ
Φ
〈 Φ〉

(3.3)

〈 ⋅

(3.4)

〈

〉

Φ〉

2

Φ

(3.5)

Φ

ここで、d、n0、、r、w、および はそれぞれ、次元数、初期粒子数密度、スカラー量、位置ベクトル、重み関数、ベクトル量およびラプラシアン定数である。2 粒子の中点で発散を計算する場
合は、式(3.4)の代わりに次式を用いる。

〈 ⋅

〉

2

(3.6)

初期粒子数密度およびラプラシアン定数は粒子が格子状に並んだ状態を仮定し、次式によって求めた値を使い続ける。

(3.7)
1

(3.8)

MPS 法における重み関数は、影響半径内でのみ正の値をもち、原点において無限大に発散する関数を用いている。本研究では次式に示す重み関数を用いる[1], [2]。

2

0

(3.9)

0

3-22

0

(3.10)

0
ここで、R、Reff および wgrad はそれぞれ、粒子間距離、影響半径および圧力勾配項の計算にのみ用いる重み関数である。圧力勾配項の計算には次式に示すように、式(3.10)の重み関数により計算した
初期粒子数密度を用いる。

(3.11)
影響半径は、粒子径の定数倍で設定される。本研究では、影響半径を粒子径の 2 倍に設定した。

③

圧力勾配項
Explicit-MPS 法では、圧力を粒子数密度の関数として次式のように計算する[1], [2]。

17

(3.12)0
ここで、 c および n は音速および粒子の粒子数密度である。粒子数密度は式(3.7)と同様にして、粒子ごとに毎ステップ計算する。Explicit-MPS 法において、音速 c は流体の圧縮性を決定するパラメータとして、ユーザー側で決定する。マッハ数| |/ が 1 よりも十分に小さい場合、すなわ

ち音速を十分に大きく設定すれば流れは非圧縮とみなすことができる。一方、マッハ数が 1 に近い
場合は微圧縮の仮定が成り立たないため、解析精度が十分確保できない可能性がある。圧力勾配項
は、式(3.12)で求めた圧力を用いて、次式で計算される。

(3.13)
本研究では、圧力勾配項は陽的に、粘性応力項は陰的に計算を行っている。式(3.13)によって圧力勾配を計算した後、次式によって中間速度を求める。
∗

1(3.14)

Δ④

粘性応力項
粘性応力項は、次式のように離散化される。
22⋅(3.15)

メッシュフリー粒子法である SPH 法および MPS 法において、式(3.15)の離散化はよく用いられている[3]–[6]。式(3.14)で求めた中間速度 u*を用いて、次式のように速度を陰的に計算する。本研究では、ソルバーに前処理付き共役勾配法を用いた。
∗

Δ(3.16)⋅3-23

式(3.14)および式(3.16)は、圧力勾配項および粘性応力項を同時に陰的に計算する Fractional
Time Step 法と同様に陽的に更新する部分と陰的に更新する部分を分けて計算している。本研究では
圧力項は陽的に計算しているため、Kim & Moin[7] とは異なりポテンシャル関数のポアソン方程式を解く必要がない。

(2) 熱伝導モデル①

支配方程式
MPS 法を用いて熱伝導問題を扱う際、連続の式と Navier-Stokes 方程式に加えて、次式に示すエ

ネルギー保存則を計算する。

(3.17)⋅

ここで、h、k、T および Q はそれぞれ、エンタルピー、熱伝導率、温度およびソース項である。式(3.17)の右辺第一項は、粘性応力項と同様に次式のように離散化される。
2
2
⋅

②

(3.18)

温度-エンタルピーモデル
粒子の温度は、既存の手法[3]と同様に、エンタルピーの関数として次式のように計算される。

(3.19)

ここで、Tm、hs、hl および cp はそれぞれ、融点、固相のエンタルピー、液相のエンタルピーおよび比熱である。式(3.19)による温度-エンタルピーモデルを図 3.3-1 に示す。このモデルでは、相変化過程における粒子温度を融点で一定であると仮定している。固相および液相のエンタルピーは次式で求められる。

(3.20)(3.21)
ここで、L は潜熱である。

③

相状態-エンタルピーモデル
粒子の相状態についても温度と同様に、エンタルピーの関数として次式のように計算される。

1(3.22)
03-24

ここで、 は粒子の相状態係数である。 が 1 のとき粒子の相状態は固相を、0 のとき液相を示し、相変化過程における相状態係数は線形補間した。式(3.22)による相状態-エンタルピーモデルを図3.3-2 に示す。

④

粘性係数-相状態モデル
固液間の相変化では、密度変化は気液間相変化と比較して無視できるほど小さい。しかし、粘性

係数は固相および液相で非常に大きく変わるため、固液間相変化を計算する際には粘性係数のモデル化が必要となる。本研究では、粘性係数を相状態係数の関数として次式に示す Arrhenius の式を用いた。

(3.23)

exp 2.5
ここで、r および C はそれぞれ、液相の粘性係数に対する比粘度 (/ ) および経験定数である。経験定数 C は 4 から 8 の間の値を取ることが知られている[8]。本研究では、固相を超高粘性の流体として模擬するため、C の値を一律に 8 に設定した。このとき、固相における粘性係数は液相のおよそ 108 倍程度まで大きくなる。

⑤

境界条件

(ア) Dirichlet 境界条件
Dirichlet 境界条件では、境界面で温度を固定する。従って、壁面からの熱輸送は式(3.20)に壁面の温度を代入すればよいため、ソース項 Q は(3.24)

0として計算される。壁面を粒子ではなく SDF を用いて設定した場合、壁面からの熱輸送は式(3.18)の代わりに次式で離散化される。

2

⋅

2

(3.25)(イ) Neumann 境界条件
Neumann 境界条件では、境界面で熱流束を固定する。MPS 法のようなメッシュフリー粒子法では、自由液面の面積を正確に計算することが困難なため、輻射伝熱も環境温度固定として Dirichlet
境界条件として設定されてきた[3], [4], [9], [10]。本研究では、熱流束の発散を計算することによって熱流入量を求める発散モデル[11]を用いた。発散モデルを用いると、ソース項 Q は⋅(3.26)⋅

で求められる。ここで q は計算領域内部に設定した熱流束である。Neumann 境界条件を設定するには、熱流束ベクトルの方向に応じて表面を判定する必要がある。表面粒子の判定は粒子数密度によって行われ、次式を満たす粒子を表面粒子として扱う。

(3.27)

Γ3-25

ここでは表面判定のための係数であり、本研究では先行研究[12], [13]と同じ 0.97 に設定した。
次に熱流束境界条件のための面法線ベクトル N h e a t を次式によって求める。

∑(3.28)∑
式(3.28)で求めた面法線ベクトルを用いて、次式を満たす粒子が熱流束 q によって境界面となる。

⋅(3.29)

0(3) 表面張力モデル
本研究では、表面張力を Continuum Surface Force (CSF)モデルで計算する。CSF モデルでは、表面張力を境界条件として与えるのではなく界面における力として計算する[14]。MPS 法では、式(3.2)の外力項は次のように表される[15], [16]。

(3.30)
ここで、g、、、および N s u r f はそれぞれ、重力加速度、表面張力係数、曲率、デルタ関数および表面張力項における面法線ベクトルである。
表面張力を界面における力として計算するため、まず式(3.27)によって表面粒子の判定を行う。
表面粒子において、表面張力項の計算には、次式に示す粒子数密度および重み付け関数を用いる。
,

,

1
0

,

,

(3.31)

0(3.32),(3.33),,

1 0(3.34)
0
式(3.27)によって表面と判定される粒子は必ずしも 1 層とはならず、図 3.3-3 に示すように多層と,なることがある。式(3.31)および式(3.32)の式では、表面と判定された粒子のうち、内側に位置す
る粒子は外側の粒子よりも大きな値をもつことになり、過大評価されてしまう可能性がある。そこ
で、式(3.33)および式(3.34)によって粒子数密度を再評価することによって、表面判定の精度を向上させている。これら再計算された値を用いて、曲率を次式で求める。

2 cos(3.35)

また、θは,2(3.36),

3-26

で求められる。nsurf,0 は表面張力の計算における初期粒子密度であり、曲率が 0 であるような平らな状態に粒子を配置し、計算の最初に求めておく。
表面張力項における面法線ベクトルを計算するために、図 3.3-4 に示すように粒子 i の周りに仮想的に粒子を配置する。周囲の粒子は、それぞれ粒子径だけ軸方向に離れた位置にする。これらの位置に実際の粒子は存在しないため、粒子数密度を次式で計算する。
,,(3.37)

周囲の粒子について、式(3.37)を用いて得られた粒子数密度をそれぞれ nx+、nx-、ny+、ny-、nz+および nz-として、面法線ベクトルを次式で求める。

22(3.38)

2(3.39)| |最後に表面張力を規格化する必要がある。これは既に述べたとおり、界面において表面と判定される粒子は何層かの厚みを持つためである。先行研究[15]ではΓ

0.97とした場合、厚みと粒子径の比を 1.55 として表面張力項をこの値で除算して規格化しており、本研究でもこの値を採用した。

(4) 共晶反応モデル
福島第一原子力発電所を含む BWR では、燃料棒に B4C を用いている。B4C(融点：2,620 K)はステ
ンレス鋼との間で物質拡散が進むことによって、ステンレスの融点である 1,673～1,723 K より低い
温度で溶融が始まることが知られている(共晶反応)[17]。そのため BWR において、シビアアクシデ
ント発生時には最初に制御棒が落下すると考えられている。シビアアクシデント発生時における事
故進展を詳細に模擬するためには、共晶反応のモデル化が必要となる。
共晶反応のモデル化にあたって、元素の存在比をもとに相図を用いて共晶点を決定する必要があるが、今回対象としている B4C とステンレス鋼では拡散係数のデータが存在しないため、物質拡散に基づくモデル化ができない。本研究では、影響半径内部に異なる種類の粒子がいる場合に融点が低くなるという単純なモデル化を行った。MPS 法では、粒子の温度や相状態はエンタルピーによって管理される。影響半径内部に異なる種類の粒子が存在する場合、固相のエンタルピーは式(3.20)の代わりに次式で計算される。

(3.40)
ここで Te は共晶点である。

(5) 壁面モデル(5-1) 符号付距離関数法(SDF)
MPS 法では、粒子を並べて壁面との相互作用を計算するが、この場合曲率をもつような複雑な体系で計算することが困難という問題点があった。本研究では任意形状を模擬することのできる符号付距離関数法[18]を用いて、壁面との距離と法線ベクトルを求める。

3-27

SDF では、計算領域を立方格子に分割し、格子点ごとに壁面までの符号付距離を計算する。符号付距離関数は次式で定義される。

(3.41)⋅

ここで、、d、s および x はそれぞれ、符号付距離関数、壁面からの距離、計算領域内外を示す符号および座標である。符号は次式で定義される。

11(3.42)

計算領域内の任意の点において、壁面までの距離は各格子における距離の線形補間で求める。法線ベクトルは次式で計算される。

|(3.43)|(5-2) 壁重み関数
壁面内部における重み関数は壁内部に粒子が充填されているとして計算する。重み関数は一定の距離ごとにあらかじめ次式で計算しておき、計算領域内部の任意の点における壁面からの重み関数の値は、壁面からの距離に応じて線形補間して求める。

(3.44)

cos(3.45)

ここで、W、Wgrad、Rwall および  はそれぞれ、壁重み関数、圧力勾配項に用いる壁重み関数、壁面からの距離および壁面の法線と粒子方向のなす角である。壁内部では、粒子が単純立方格子状に充填していると仮定しており、本研究でも同じ手法を用いた。

3.2.2 シミュレーション(H26-H28)(1) 解析条件
図 3.3-5 に示すように粒子を立方体状に並べて構造物を作り、計算領域内部に鉛直下向きの熱流束場を設定した。構造物の一辺の長さは 0.10 m とした。図 3.3-5 に示したように、一部異なる粒子を配置した。計算条件を表 3.3-1 に示す。時間刻み、重力加速度、音速、粒子径、影響半径、粘性
経験定数、熱流束はそれぞれ、10-5 s、9.8 m/s2、100 m/s、5.0 mm、10.0 mm、8.0、40,000 W/m2 である。粒子は x、y および z 軸方向にそれぞれ、20 個ずつ並べ、計 8,000 個の粒子で計算を行った。
粒子の物性値を表 3.3-2 に示す。Material1,2 ともに同じ物性値であり、その密度、動粘性係数、比熱、熱伝導率、融点、共晶点および潜熱はそれぞれ、1,000 kg/m3、10-3 m2/s、100 J/(kg K)、50
W/(m K)、10,000 K、320 K および 10,000 J/kg である。壁面は全て 300 K の Dirichlet 境界条件を与えた。高粘性の溶融物がロウソクのように壁面を伝って流れ落ちるキャンドリング現象を模擬す

3-28

るため、Case 3 では表面張力モデルを導入した。表面張力係数は水と同じ 0.073 N/m とした。壁面
における境界条件はすべて 300 K の温度固定境界とした。

3.2.3 結果・考察(H26-H28)
図 3.3-6 に結果を示す。Material1,2 ともに融点は 10,000 K と、設定した熱流束場では溶融が始まらない温度である。本研究では共晶反応を模擬するために、影響半径内部に異なる種類の粒子が存在するとき融点は共晶点まで下がるようにしている。図が示すように、角にある異物質の境界でのみ溶融が起こっており、簡易的にではあるが共晶反応による溶融現象の模擬に成功した。表面張力モデルの導入により、1.5 sec の図が示すように、溶融物が壁面を伝うように流れ落ちている。本
解析では底部壁面に温度固定境界を設定している。溶融物は底部壁面まで流れ落ちた後再凝固しており、本モデルによって、共晶反応を考慮した炉内構造物の溶融および再配置の模擬が可能となった。
共晶反応は、境界面において物質拡散が進むことによって起こることがわかっている。詳細計算のためには、物質拡散の評価が必要である。しかしながら共晶反応発生時における物質拡散の定量的な評価は未だ報告されていないため、本研究では単純なモデルを用いて検証を行った。シビアア
クシデント発生時には炉内の温度は非常に高温となるため、溶融現象は短時間で進展することが知られている。従ってシビアアクシデント解析では、今回のような単純なモデルでも精度よく計算できると考えられる。残念ながら、3.1 節で実施した実験で得られたデータと同様の体系での詳細計算はできなかったが、定性的な比較としては共晶の状態を評価することができた。

3-29

表 3.3-1 計算条件(Case 3)

Item

Unit

Value

Time steps

10-5

Gravitational acceleration

m/s2

9.8

Sound speed

m/s

100

Particle diameter

mm

5.0

Effective radius

mm

10.0

Empirical constant for viscosity

-

8.0
2

Intensity of heat flux

W/m

40,000

表 3.3-2 物性値(Case 3)

Item

Unit

Material 1

Material 2

Density

kg/m3

1,000

1,000

2

-3

Kinematic viscosity

m /s

10

10-3

Specific heatJ/(kg K)

10

10

Thermal conductivity

W/(m K)

50

50

Melting point

K

10,00010,000

Eutectic temperature

K

320

320

Latent heat

J/kg

10,00010,0003-30

図 3.3-1 MPS 法における温度-エンタルピーモデル

3-31

図 3.3-2 MPS 法における相状態-エンタルピーモデル

3-32

図 3.3-3 CSF モデル概略図

3-33

図 3.3-4 CSF モデルにおける仮想粒子3-34

図 3.3-5 粒子初期配置(Case 3)

3-35

0 sec

1.5 sec

3.0 sec

5.0 sec

図 3.3-6 シミュレーション結果図(Case 3)

3-36

参考文献
[1]

Y. Yamada and M. Sakai, “Lagrangian-Lagrangian simulations of solid-liquid flows
in a bead mill,” Powder Technol., vol. 239, pp. 105–114, 2013.

[2]

X. Sun, M. Sakai, M. T. Sakai, and Y. Yamada, “A Lagrangian-Lagrangian coupled
method for three-dimensional solid-liquid flows involving free surfaces in arotating cylindrical tank,” Chem. Eng. J., vol. 246, pp. 122–141, 2014.

[3]

G. Li, Y. Oka, and M. Furuya, “Experimental and numerical study of stratificationand solidification/melting behaviors,” Nucl. Eng. Des., vol. 272, pp. 109–117,
2014.

[4]

R. Chen, Y. Oka, G. Li, and T. Matsuura, “Numerical investigation on melt freezing
behavior in a tube by MPS method,” Nucl. Eng. Des., vol. 273, pp. 440–448, 2014.

[5]

P. W. Cleary, “Extension of SPH to predict feeding , freezing and defect creation
in low pressure die casting,” Appl. Math. Model., vol. 34, no. 11, pp. 3189–3201,
2010.

[6]

M. Tong and D. J. Browne, “International Journal of Heat and Mass Transfer An
incompressible multi-phase smoothed particle hydrodynamics ( SPH ) method formodelling thermocapillary flow,” Int. J. Heat Mass Transf., vol. 73, pp. 284–292,
2014.

[7]

J. Kim and P. Moin, “Application of a fractional-step method to incompressible
Navier-Stokes equations,” J. Comput. Phys., vol. 59, no. 2, pp. 308–323, 1985.

[8]

M. Ramacciotti, C. Journeau, F. Sudreau, and G. Cognet, “Viscosity models forcorium melts,” Nucl. Eng. Des., vol. 204, no. 1–3, pp. 377–389, 2001.

[9]

G. Li, M. Liu, G. Duan, D. Chong, and J. Yan, “Numerical investigation of erosionand heat transfer characteristics of molten jet impinging onto solid plate with
MPS-LES method,” Int. J. Heat Mass Transf., vol. 99, pp. 44–52, 2016.

[10]

Y. Liang, Z. Sun, G. Xi, and L. Liu, “Numerical models for heat conduction andnatural convection with symmetry boundary condition based on particle method,”Int. J. Heat Mass Transf., vol. 88, pp. 433–444, 2015.
[11]

K. Takabatake, X. Sun, M. Sakai, D. Pavlidis, J. Xiang, and C. C. Pain,“International Journal of Heat and Mass Transfer Numerical study on a heattransfer model in a Lagrangian fluid dynamics simulation,” Int. J. Heat Mass

Transf., vol. 103, pp. 635–645, 2016.
[12]

S. Koshizuka, A. Nobe, and Y. Oka, “Moving-particle semi-implicit method forfragmentation of incompressible fluid,” Nucl. Sci. Eng., vol. 123, pp. 421–434,
1996.

3-37

[13]

X. Sun, M. Sakai, K. Shibata, Y. Tochigi, and H. Fujiwara, “Numerical modeling on
the discharged fluid flow from a glass melter by a Lagrangian approach,” Nucl.

Eng. Des., vol. 248, pp. 14–21, 2012.
[14]

J. U. Brackbill, D. B. Kothe, and C. Zemach, “A continuum method for modelingsurface tension,” J. Comput. Phys., vol. 100, no. 2, pp. 335–354, 1992.

[15]

K. Nomura, S. Koshizuka, Y. Oka, and H. Obata, “Numerical Analysis of Droplet
Breakup Behavior using Particle Method,” J. Nucl. Sci. Technol., vol. 38, no. 12,
pp. 1057–1064, 2001.

[16]

Z. Sun, G. Xi, and X. Chen, “Numerical simulation of binary collisions using amodified surface tension model with particle method,” Nucl. Eng. Des., vol. 239,
no. 4, pp. 619–627, 2009.

[17]

A. HIDAKA, “Effect of B4C Absorber Material on Melt Progression and Chemical Formsof Iodine or Cesium under Severe Accident Conditions,” Trans. At. Energy Soc.

Japan, vol. 14, no. 1, pp. 51–61, 2015.
[18]

Y. Shigeto and M. Sakai, “Arbitrary-shaped wall boundary modeling based on signed
distance functions for granular flow simulations,” Chem. Eng. J., vol. 231, pp.
464–476, 2013.

3-38

3.3 まとめ(H26-H28)
本研究では、シビアアクシデント時の事故進展を模擬するため、東京大学及び英国 Imperial
College London（以下、「IC」という。）が協力して、制御棒の共晶過程及び再配置を模擬するための詳細な熱流動シミュレーション技術をそれぞれ開発して、検証（Verification）と妥当性確認（Validation）に関する研究を行った。

(1) 調査研究及び検証実験
平成 26 年度
①調査と予備実験

シビアアクシデントの熱流動シミュレーションで使用する、ホウ素と鉄の共晶状態に関する相図を調査した。また、シミュレーションの検証データを取得するための共晶反応実験装置を製作し、可視化手法に関する予備実験を行った。
平成 27 年度
②実験実施とフィードバック

熱流動シミュレーションコードの妥当性確認（Validation）のため、実験条件が適切に数値解析に反映できることを検討した。平成 26 年度に製作した共晶反応実験装置を用いて、可視化手法に関する実験を行い、実験データを取得した。
平成 28 年度
③実験データによる評価

熱流動シミュレーションコードの妥当性確認（Validation）のため、平成 26 年度に製作した共晶反応実験装置を用いて、実験を行い、実験データを取得するとともに、共晶反応について検討を行った。共晶実験結果を整理し、validation 用のデータを整備した。

(2) 熱流動シミュレーション
平成 26 年度
① 基盤技術開発

MPS 法を使用して熱流動シミュレーションのコードを開発した。基礎的な数値実験として、
シンプルな体系において高粘性流体の数値シミュレーションを実行した。IC と協力して、数値シミュレーション結果の妥当性について考察した。
平成 27 年度
②検証（Verification）

3-39

平成 26 年度に開発した粒子法を使用した熱流動シミュレーションコードについて、妥当性を検証した。検証実験を参考にして、計算体系および計算条件を選定した。IC と協力して、数値シミュレーション結果の妥当性について考察するとともに、手法間の特徴を整理した。
平成 28 年度
③妥当性確認（Validation）

MPS 法を用いて、これまでの実験成果を参考にして共晶反応を考慮した固体溶融に関する数
値シミュレーションを実行した。数値シミュレーションにより、制御棒の溶融に係わる実験で観察される典型的な現象のキャンドリングを模擬できることを確認した。
(3) 研究推進
平成 26 年度

研究代表者の下で各研究項目間における連携を密にして研究を進めた。
平成 27 年度

研究代表者の下で各研究項目間における連携を密にして研究を進めた。研究成果は、国際会議などで発表した。
平成 28 年度

研究代表者の下で各研究項目間における連携を密にして研究を進めた。研究成果は、国際会議などで発表した。また、J. Nucl. Sci. Technol.および Int. J. Heat Mass Transf.において、成果を発表した。
以上、３年目として本年度の業務項目を実施し、所期の目標を達成した。
今後の展望として、実際の原子炉構造物を用いた数値シミュレーションの実行、熱流動解析に必要となる詳細な物性データの取得などの研究があげられる。

3-40

4.

結言
シビアアクシデント発生時には、B4C とステンレス鋼の間で共晶反応が進むことによって本来の融点よりも低い温度で溶融が始まり、最初に制御棒が落下すると考えられている。溶融現象は自由液面の大変形を伴うため、伝熱モデルとともに自由液面流れのモデル化が必要となる。本研究においては、実験によって、共晶現象を可視化し、数多くの実験データを取得するとともに、シミュレ
ーションコードの開発を実施した。
実験においては、BWR 制御棒材料を対象とした超高温環境下の共晶反応可視化手法を開発した。
開発した共晶反応可視化手法を用いて B4C /SS の溶融進展と溶融物落下挙動を可視化した。その結果として、溶融進展については、未反応 B4C 粉末が融液内を移流していくことで、局所的に溶融進展が加速されることが分かった。また、溶融進展の最中は、反応領域において固相と液相が共存することを実験的に確認した。次に、溶融物落下挙動については、複数のモードが確認されたが、大きく分けて、膜を形成して元の形状を維持するモードと液滴を形成して大きく崩落するモードに分類できる。本研究の実験条件では、上記モードの決定因子は SS 部分の断面積と B4C /SS の組成比であった。さらに、以上の実験データについて，シミュレーションで必要となる実験中の熱的な境界条件を得ることができた。
シミュレーションについては、自由液面の模擬に優れており、炉内構造物の溶融および再配置の計算に適した MPS 法を用いた。本研究では MPS 法における Neumann 境界条件の設定手法を開発した。既存の MPS 法では、Dirichlet 境界条件しか設定できなかったが、本モデルを用いることで、自由液面に直接 Neumann 境界条件を設定することが可能となり、熱流束による構造物の溶融および溶融物の再凝固を模擬できた。溶融現象は自由液面の大変形を伴うため、シビアアクシデント時の炉内構造物の溶融および再配置の計算が既存の手法よりも精度よく計算できるようになった。最後に、シビアアクシデント発生時に問題となる共晶反応のモデル化を行った。影響半径内部に異なる種類の粒子が存在していれば融点が低くなるという単純なモデルであるが、境界面から溶融が始まる様子が模擬できており、共晶反応による溶融現象を計算することが可能となった。
本事業においては、インペリアルカレッジロンドンと、相互に出張をすることでミーティングを実施したり、遠隔テレビ会議システムを有効に使って、密な交流を行いつつ事業を進めてきた。実験結果やシミュレーション手法に関する詳細な議論の結果が、本事業の成果に反映されている。な
お、本事業の範囲では、実験体系をシミュレーションで完全に模擬することは目的としていなかったため、実験において取得したデータを、シミュレーションによって定量的に比較するまでには至らなかったが、現象を定性的に比較し、同等の流動状況が得られていることを確認している。
シビアアクシデントは高温・高放射線量という過酷な環境で、かつ、非線形性が極めて大きな問題であり、現象の解明には数値シミュレーション技術の導入が期待されている。本研究で開発したモデルによって、シビアアクシデントの一つである共晶反応の評価が可能となった。これら詳細計

4-1

算の結果を、シビアアクシデント解析コードに反映させることによって、事故進展のより詳細な解析が期待できる。

4-2