平成 30 年度

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

燃料デブリ取り出し戦略の構築：リスク管理と物理シミュレーションの融合

成果報告書

平成 31 年 3 月
国立大学法人 東京大学

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

東京大学が実施した平成 28-30 年度「燃料デブリ取り出し戦
略の構築：リスク管理と物理シミュレーションの融合」の成果を取りまとめたものです。

目次
1.

はじめに .................................................................. 1-1

2.

業務計画

2.1

全体計画 ............................................................ 2.1-1

2.2

平成 30 年度の成果の目標および業務の実施方法 ......................... 2.2-1

3.

業務の実施内容および成果

3.1

マルチフィジックスシミュレーション .................................. 3.1-1

3.2

リスク評価 .......................................................... 3.2-1

3.3

英国共同研究先の研究紹介 ....... エラー! ブックマークが定義されていません。

3.4

研究推進 ............................................................ 3.4-1

4.

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

1

表一覧
頁
表 3.1-1 液滴落下解析のパラメータ ・・・・・・・・・・・・・・・・・・・・・3.1-6
表 3.1-2 計算条件

・・・・・・・・・・・・・・・・・・・・・・・・・・・・3.1-11

表 3.1-3 原子数密度の条件

・・・・・・・・・・・・・・・・・・・・・・・・3.1-11

表 3.1-4 物性値 ・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・3.1-19
表 3.1-5 計算条件 ・・・・・・・・・・・・・・・・・・・・・・・・・・・・・3.1-19
表 3.1-6 切り粉漏洩評価の計算条件・・・・・・・・・・・・・・・・・・・・・3.1-23
表 3.1-7 多孔質体の物性値・・・・・・・・・・・・・・・・・・・・・・・・・3.1-24
表 3.1-8 ガスの物性値 ・・・・・・・・・・・・・・・・・・・・・・・・・・・3.1-26
表 3.1-9 装置の物性値 ・・・・・・・・・・・・・・・・・・・・・・・・・・・3.1-26
表 3.1-10 粒子の物性値・・・・・・・・・・・・・・・・・・・・・・・・・・・3.1-28
表 3-2-1 シナリオの説明 ・・・・・・・・・・・・・・・・・・・・・・・・・・3.2-17
図一覧頁
図 2.1-1 全体計画図

・・・・・・・・・・・・・・・・・・・・・・・・・・・2.1-1

図 3.1-1 液滴落下実験の概要

・・・・・・・・・・・・・・・・・・・・・・・3.1-1

図 3.1-2 実験装置の KYOWA Drop-Master DMo-501 接触角計 ・・・・・・・・・・・3.1-2
図 3.1-3 低落下速度における粒子層への液滴付着挙動（フレーム時間間隔 8 ms） ・3.1-3
図 3.1-4 高落下速度（0.4 m/s 程度）における挙動（フレーム時間間隔 8 ms） ・・3.1-4
図 3.1-5 高落下速度における液滴挙動のハイスピード撮影（フレーム間隔 1 ms）
・・3.1-5
図 3.1-6 液滴落下シミュレーションの結果

・・・・・・・・・・・・・・・・・3.1-7

図 3.1-7 固気液混相流計算結果の比較・・・・・・・・・・・・・・・・・・・・3.1-12
図 3.1-8

keff の経時変化 ・・・・・・・・・・・・・・・・・・・・・・・・・・3.1-13

図 3.1-9 液架橋の幾何学的形状

・・・・・・・・・・・・・・・・・・・・・・3.1-14

図 3.1-10 粒子-粒子間の液架橋モデル ・・・・・・・・・・・・・・・・・・・・3.1-17
図 3.1-11 回転ドラム装置図・・・・・・・・・・・・・・・・・・・・・・・・・3.1-19
図 3.1-12 湿潤粉体の挙動予測モデルによる計算結果・・・・・・・・・・・・・・3.1-20
図 3.1-13 粒子速度分布・・・・・・・・・・・・・・・・・・・・・・・・・・・3.1-21
図 3.1-14 球状粒子空隙中の浸潤液面 ・・・・・・・・・・・・・・・・・・・・3.1-24
図 3.1-15 液滴と多孔質平板の衝突および浸入 ・・・・・・・・・・・・・・・・3.1-25
図 3.1-16 気流の流速ベクトル（固体粒子なし）・・・・・・・・・・・・・・・ 3.1-27
図 3.1-17 ガス温度（固体粒子なし）・・・・・・・・・・・・・・・・・・・・ 3.1-27
図 3.1-18 気流の流速ベクトル（固体粒子あり）・・・・・・・・・・・・・・・ 3.1-29
図 3.1-19 ガス温度（固体粒子あり）・・・・・・・・・・・・・・・・・・・・ 3.1-30
図 3.1-20 粒子速度[m/s] ・・・・・・・・・・・・・・・・・・・・・・・・・3.1-31
図 3.2-1 リスクガバナンスの枠組み・・・・・・・・・・・・・・・・・・・・・3.2-1ii

図 3.2-2 燃料デブリ取出しのエンドポイントと枠組み・・・・・・・・・・・・・3.2-1
図 3.2-3 炉内安全確保に関するリスク構造 ・・・・・・・・・・・・・・・・・・3.2-2
図 3.2-4 燃料デブリ切断に関するリスク構造 ・・・・・・・・・・・・・・・・・3.2-2
図 3.2-5 リスク評価とリスク管理のための 6 つのステップ ・・・・・・・・・・・3.2-4
図 3.2-6 脅威の体系的分類 ・・・・・・・・・・・・・・・・・・・・・・・・・3.2-5
図 3.2-7 起因事象に関するリスク評価指標 ・・・・・・・・・・・・・・・・・・3.2-5
図 3.2-8 失敗シナリオに関するリスク評価指標 ・・・・・・・・・・・・・・・・3.2-6
図 3.2-9 シナリオの優先度付けに関するリスク評価指標 ・・・・・・・・・・・・3.2-6
図 3.2-10 デブリ取り出し時のリスク要因（液相部の閉じ込めについて）・・・・3.2-7

図 3.2-11 デブリ取り出し時のリスク要因（気相部の閉じ込めについて）・・・・・3.2-8
図 3.2-12 格納容器内部と原子炉圧力容器内部の解析モデル・・・・・・・・・・・3.2-9
図 3.2-13 格納容器内部の流動状況・・・・・・・・・・・・・・・・・・・・・・3.2-10
図 3.2-14 格納容器内の温度分布・・・・・・・・・・・・・・・・・・・・・・・3.2-11
図 3-2-15 出口境界近傍の温度過渡変化・・・・・・・・・・・・・・・・・・・・3.2-11
図 3.2-16 温度計１と温度計２の温度過渡変化 ・・・・・・・・・・・・・・・・3.2-12
図 3.2-17 温度計３と温度計４の温度過渡変化 ・・・・・・・・・・・・・・・・3.2-12
図 3.2-18 各温度計の温度過渡変化より求めた異常と診断される確率 ・・・・・・3.2-13
図 3.2-19 温度計 1 と 2 を組み合わせて異常と診断される確率 ・・・・・・・・・3.2-13
図 3.2-20 温度計 3 と 4 を組み合わせて異常と診断される確率 ・・・・・・・・・3.2-14
図 3.2-21 温度計を３つあるいは４つを組み合わせて異常を診断する確率・・・・ 3.2-14
図 3.2-22 数値解析で得られた温度と各シナリオで推定される温度の比較・・・・・3.2-15
図 3.2-23 各温度計のシナリオごとの信頼度

・・・・・・・・・・・・・・・・・3.2-16

図 3.2-24 1 カ所の温度を用いた場合と 4 カ所の温度を用いた場合の異常ありシナ
リオの信頼度 ・・・・・・・・・・・・・・・・・・・・・・・・・・3.2-17
図 3.2-25 4 カ所の温度を用いてシナリオ推定した結果・・・・・・・・・・・・・3.2-18
図 3.2-26 4 カ所の温度と流速を用いたシナリオ推定の結果・・・・・・・・・・・3.2-18
図 3.3-1 蒸気の体積分率

時刻 = 0.04, 0.05 s 摂動温度104 ℃・・・・・・・・・3.3-6

図 3.4-1 意見交換会後の集合写真・・・・・・・・・・・・・・・・・・・・・・ 3.4-1
図 3.4-2 山口彰研究代表者による講演・・・・・・・・・・・・・・・・・・・・ 3.4-1

iii

概略
福島第一原子力発電所の廃炉の加速に資するため、世界でも例のない事故炉からの燃料デブリの安全かつ円滑な取出しに貢献する研究について、英国 Imperial College London（ICL）との協力により推進する。燃料デブリの取り出しにおいて、デブリの冷却や放射線の遮蔽において優れている冠水工法の実施が有力であると考えられているが、本工法には格納容器の止水、耐震性、臨界管理等に問題がある。他方、冠水工法の課題をクリアした気中工法も燃料デブリの取り出しにおいて有力な技術であると考えられているが、これまでに技術検討がほとんどなされていなかった。そこで、本研究では、ICL との共同研究により、気中工法に関する安全な燃料デブリの取り出しに係わるマルチフィジックスシミュレーションとその結果に基づくリスク評価により、燃料デブリを長期管理しつつ気中取り出し方法の問題点克服を検討する。

iv

1.

はじめに

燃料デブリの取り出しにおいては、燃料デブリの冷却、従業員の被ばく防止、環境への放射性物質の放出防止が、重要な留意点である。本プロジェクトでは気中工法に注目する。本プロジェ
クトは、マルチフィジックスシミュレーションに関する研究とリスク評価に関する研究より構成される。本研究は、東京大学および ICL で実行したマルチフィジックスシミュレーション結果に基づいてリスク評価を行い、燃料デブリを長期管理しつつ気中取り出し方法の問題点克服を検討する。
マルチフィジックスシミュレーションに関する研究の内容は、①マルチフィジックスシミュレータを使用して、デブリ切断時における切り粉飛散防止スプレーの保守的な条件を決定する、②固気液三相流れの数値シミュレーションと連続エネルギーモンテカルロ法コードを接続した新しい数値解析手法を開発して、気中工法が臨界の観点から問題ないことを確認する、③気中工法における燃料デブリの冷却性能評価および飛散したデブリの切り粉の外部への漏洩の可能性の検討を行う、である。リスク評価に関する研究は、マルチフィジックスシミュレーション結果に基づいており、①燃料デブリ取り出しの技術リスクを見積もるための簡易システムの基本的方針の確
立、②燃料デブリ取り出しの技術リスク、プロジェクト遂行リスク、社会的意思決定リスクをバ
ランス良く評価するための考え方の提示、および③日英で取り組んだマルチフジックスシミュレーション結果に基づき、リスクモニタリングとリスク管理のための考え方の提案を行う。

1-1

2.

業務計画
全体計画
本業務の全体計画図を図 2.1-1 に示す。
図 2.1-1 全体計画図

実施項目

平成 28 年度

平成 29 年度

平成 30 年度

(1) マルチフィジックスシミュレーション①デブリ切断時の切り粉拡散抑制評価②固気液混相流-臨界連成手法の開発③デブリ冷却性能および切り粉漏洩評
価(2) リスク評価①技術リスク
②技術リスク、プロジェクト遂行リス
ク、社会的意思決定リスク③リスクモニタリング・リスク管理(3)研究推進平成 28 年度
【マルチフィジックスシミュレーション】①燃料デブリ取り出し時の燃料切り粉の飛散防止のためのスプレーの条件を検討するため、固気液三相流シミュレーションを専用のワークステーションで実施する。②シミュレーション結果の妥当性を検証するため、接触角計および既存の装置を組み合わせて実験を行いその結果を整理する。③ICL と計算結果の妥当性の検証方法について確認すると共に実際の気中工法への応用方
法について検討しこの結果を整理する。
【リスク評価】

2.1-1

燃料デブリ取り出しの手順が決まった場合のそのプロセスの技術リスクを見積もるための簡易シ
ステムを構築するに際しての基本的な方針を検討する。
【研究推進】
研究代表者の下で、ICL および東京大学内の各研究項目担当者間における連携を密にして研究を進める。
平成 29 年度
【マルチフィジックスシミュレーション】
研究分担者の開発した手法を発展させて、固気液三相流れの数値シミュレーションと連続エネ
ルギーモンテカルロ法コードを接続して、切り粉の堆積状態と臨界安全との関係についても検討する。限られた条件下で、日本側の計算結果と英国側（ICL）の計算結果を比較し、計算結果の妥当性を検証する。
【リスク評価】
燃料デブリ取り出しに係る戦略の具体化を合理的に進めるためには、技術的リスクのみでなく、
プロジェクトの進捗に係るリスク、社会の理解醸成に係るリスクは注視すべき項目である。これ
らリスクを燃料デブリ取り出しに係る意思決定の中にどのように取り組むべきかを検討し、リスク指標として整理する。この成果により、燃料デブリ取り出しの技術リスク、プロジェクト遂行
リスク、社会的意思決定リスクをバランス良く評価するための考え方を提示する。
【研究推進】
研究代表者の下で、ICL および東京大学内の各研究項目担当者間における連携を密にして研究を進める。
平成 30 年度
【マルチフィジックスシミュレーション】
実際のジオメトリーを参考にして、気中工法における燃料デブリの冷却性能評価用の熱流動解析を行う。さらに、デブリの切り粉が飛散した場合の原子炉外部への漏洩可能性を検討する。
【リスク評価】
燃料デブリ取り出しにはさまざまな作業に関する意思決定をその都度行うために、評価されたリスクを常時モニタリングすることが有効である。最終年度には、日英で取り組んだマルチフジックスシミュレーション結果に基づき、リスクモニタリングとリスク管理のための考え方を提案する。燃料デブリを取り出す作業は、一時的にはリスクが増加する段階を経て、安定的にリスクが抑制された段階へ移行するなど、冷静にリスク分析を行なわれなければ、適正な手順を踏めない状況が生じうる。そのような作業のリスク変動の定量化とリスク低減効果を関係者が共有するために有効な、動的なリスクをモニタリングするシステムを提案する。実作業中に実測可能なパラメータを反映できる評価式を構築したい。また、燃料デブリ取り出しに係る技術リスクについてはこれまでの研究成果に基づき継続的に高度化する。
【研究推進】
研究代表者の下で、ICL および東京大学内の各研究項目担当者間における連携を密にして研究
を進める。
2.1-2

平成 30 年度の成果の目標および業務の実施方法
平成 30 年度の計画の概要は、以下に示すとおりである。
(1)マルチフィジックスシミュレーション① 固気液混相流-臨界連成手法の開発① -1 シミュレーション解析
福島第一原子力発電所の実際のジオメトリーを参考にして、気中工法における燃料
デブリの冷却性能評価用の熱流動解析を行う。
①-2 検証方法の確認等（ICL と共同）
計算結果について東京大学と ICL で議論する。
② デブリ冷却性能および切り粉漏洩評価
固気混相流の熱流動解析により、デブリの切り粉が飛散した場合の原子炉外部への漏洩可能性を評価する。
(2)リスク評価
これまでの研究成果に基づき継続的に高度化することを含め燃料デブリ取り出し作業のリスク変動を定量化し、リスク低減効果を関係者が共有するために有効なリスクをモ
ニタリングする手法を検討しまとめる。実作業中に実測可能なパラメータが反映できることを考慮する。
(3)研究推進研究代表者の下で、ICL および東京大学内の各研究項目担当者間における連携を密にして研究を進める。

2.2-1

3.

業務の実施内容および成果
マルチフィジックスシミュレーション
3.1.1 デブリ切断時の切り粉拡散抑制評価（平成 28 年度）
本研究では、接触角計を用いて、液滴落下実験を行った。また、それに相当した計算体系に、DEM-VOF-DNS 計算手法を適用し解析した。本研究の試験・検証段階で得られた実験結果と計算結果をまとめて報告する。
(1) 実験目的と体系
本研究では、自由落下する液滴が球状粒子層に接触する際の挙動を観察することを目的とする。一般の平坦基板と違い、固体粒子層に形成される多孔質が複雑な表面構造を成し、液滴の挙動に大きな影響を与える。このような現象を理解する上に、実験観察によるアプロー
チが不可欠である。そのため、現段階の定性的な観察を成功させることによって、実験手法の確立と現象の把握に重要である。
本実験の概要図を下に示す。ある高さから液滴を落下させる。下方に固体粒子層を設置し、液滴との接触挙動を観察し記録する。

図 3.1-1 液滴落下実験の概要(2) 実験設備と条件
本実験では、協和界面科学株式会社製接触角計 Drop-Master シリーズ DMo-501（以降、接触角計と記する）を使用する。その装置の写真は図 3.1-2 に示す。写真のように、接触角計に平面試料台とガラス製シリンジが装備されている。固体基板（または本研究で扱う固体粒子）を設置した後、上方のシリンジより液滴を作成し落下させる。接触角計本体にカメラと照明ライトがあり、基板上の液滴の動きを直接に撮影することができる。なお、接触角計本
3.1-1

体をパソコンに接続し、同社製の測定解析統合システムソフトウェア FAMAS（interFAce
Measurement & Analysis System）を用いた接触角測定も可能である。

図 3.1-2 実験装置の KYOWA Drop-Master DMo-501 接触角計
液体に水（室温で密度 1000 kg/m3、粘性 1.0016 Pa s）を使用した。ガラスビーズ（株式会社ユニオン製 UB-2325LN、粒径分布 1.70 ～ 2.36 mm）をシャーレに入れ固体粒子層を作成した。粒子層の表面を平らにするために、シャーレの蓋で軽くプレスした。界面物性値として、水の表面張力が 0.072 N/m、水とガラスの間の接触角が 0°に近い。
液滴を作成する際、接触角計の自動押し出し機能を利用した。回転スクリューによって連続的にシリンジをプッシュし、約 16 ± 1 μL 体積の液滴が周期的に落下する。液滴形状を球形だと仮定したところ、その直径は 3.1 mm に相当する。さらに、シリンジの針先位置を調整することによって、落下速度を変えて実験した。
(3) 実験結果
落下速度の違いに着目し、二パターンの実験結果を報告する。
① 低速度における付着挙動
まず、液滴が十分に低い速度で粒子層表面に接触するケースを観察する。シリンジの針
先を粒子層に近づけ、液滴の離脱とともに粒子表面に接触させた。このプロセスにおける液滴挙動を接触角計本体に取り付けたカメラで撮影した。液滴が変形する挙動は、図 3.1-3 のように捉えられた。画像に示されたように、液滴が不均一な粒子層表面に接近した際、表面張力と微視的接触角（水‐ガラスの場合約 0°）の作用で最初に接触した粒子の方に付着し、液滴形状に大きな偏りが生じた。また、粒子表面を迅速に濡れ始め、針先から引き離され、粒子層内部へ完全に浸入した。この液滴変形挙動は表面張力に支配され、約 20 ms の時間で
発生した。

3.1-2

(a) 接触直前：t = 0 ms(b) 接触後：t = +8 ms(c) t = +16 ms(d) t = +24 ms

図 3.1-3 低落下速度における粒子層への液滴付着挙動（フレーム時間間隔 8 ms）②

より高い落下速度における衝突挙動

次に、より高い落下速度を有する液滴の挙動を考察する。接触角計器材の制限上、シリン
ジ針先位置は数センチまで調整が可能である。微小液滴にかかる空気抵抗のため、正確に落下速度を計算するのは困難であるが、その推定値は約 0.4 m/s 程度である。これは画像における液滴中心位置の差分計算によく一致する。今後、画像処理により落下速度を算出する手法を考案する予定である。
図 3.1-4 に示された落下過程に注目したところ、液滴が粒子層表面に接触した後、速やかに粒子層内部へ浸入したことが分かる。また、次の図 3.1-5 の高精度なハイスピードカメラ（1024 * 1024 px，1000 fps）を用いた撮影はこの挙動の詳細を示した。固体層表面を通過した時間に着目すると、高速液滴の場合、接触から完全浸入まで約 8 ms しか経たなかった。
この点について前節の低速付着挙動に異なる。また、衝突過程の最中に、顕著な液滴変形も観察されなかった。
このように、液滴と固体粒子層の相互作用には、落下速度が重要なパラメータであることが示された。本研究において開発された実験手法を用いて、数値モデルの開発・検証および液滴問題の現象解明を進めたい。

3.1-3

図 3.1-4 高落下速度（0.4 m/s 程度）における挙動（フレーム時間間隔 8 ms）

3.1-4

図 3.1-5 高落下速度における液滴挙動のハイスピード撮影（フレーム間隔 1 ms）(4) 数値解析結果
実験結果との比較考察を行うために、高速度での液滴落下シミュレーションを DEM-VOFDNS 手法で実行した。ここで、液体と周囲大気の物性に水と空気の値を使用した。また、パ

3.1-5

ッキングされた固体粒子層を再現するために、分子動力学的なモンテカルロ手法でランダムな球体充填を生成した。全体 8 mm × 6 mm × 8 mm の矩形空間に、直径 2 mm の球状粒子が52 個配置された（固体体積分率約 55%）。直径 3.2 mm の液滴を粒子層の上面に配置し、初期速度を 0.4 m/s と設定した。他の解析条件の詳細は表 3.1-1 に示す。
表 3.1-1 液滴落下解析のパラメータ
Physical properties
Liquid/gas density

1000/2 kg/m3

Liquid/gas viscosity

0.001/0.00001 Pa s

Surface tension

0.07 N/mContact angle

10°

Dimensions and sizes
Domain size

8 × 10 × 8 mm

Particle diameter

2.0 mmParticle number

52

Solid layer height

6.0 mmSolid volume fraction

55%

Liquid droplet diameter

3.2 mmComputational parameters
Grid size

0.125 mmTime step

1.0 μs

解析開始後 8 ms までのスナップショットを図 3.1-6 に示す。球状液滴が固体粒子に接触
し固体層に吸収される様子が再現されている。接触過程の初期段階（約 t = 3 ～ 5 ms）では、液滴の上部がほぼ変形しなかった。また、液滴の接触から固体層表面に広がるまでの挙動は、約 10 ms の時間スケールで確認できた。これらの点については、実験による観察と定性的に一致している。一方、液滴の拡散速度は実験よりやや遅くなっている。これは、実際の落下速度と粒子分布の違いによるものだと考えられる。

3.1-6

図 3.1-6 液滴落下シミュレーションの結果

3.1-7

参考文献
[1] C. . Hirt and B. . Nichols, “Volume of fluid (VOF) method for the dynamics
of free boundaries,” J. Comput. Phys., vol. 39, no. 1, pp. 201–225, 1981.
[2] P. A. Cundall and O. D. L. Strack, “A discrete numerical model for granularassemblies,” Géotechnique, vol. 29, no. 1, pp. 47–65, 1979.
[3] J. N. Israelachvili, Intermolecular and Surface Forces, 3rd Editio. Academic
Press, 2010.
[4] T. Kajishima, S. Takiguchi, H. Hamasaki, and Y. Miyake, “Turbulence Structureof Particle-Laden Flow in a Vertical Plane Channel Due to Vortex Shedding,” JSMEInt. J. Ser. B, vol. 44, no. 4, pp. 526–535, Sep. 2001.
[5] J. H. Ferziger and M. Peric, Computational Methods for Fluid Dynamics, 3rd ed.
Berlin, Heidelberg: Springer Berlin Heidelberg, 2001.
[6] D. Gueyffier, J. Li, A. Nadim, R. Scardovelli, and S. Zaleski, “Volume-ofFluid Interface Tracking with Smoothed Surface Stress Methods for Three-Dimensional
Flows,” J. Comput. Phys., vol. 152, no. 2, pp. 423–456, 1999.
[7] W. Aniszewski, T. Ménard, and M. Marek, “Volume of Fluid (VOF) type advection
methods in two-phase flow: A comparative study,” Comput. Fluids, vol. 97, pp. 52–
73, Jun. 2014.
[8] K. Yokoi, “Efficient implementation of THINC scheme: A simple and practicalsmoothed VOF algorithm,” J. Comput. Phys., vol. 226, no. 2, pp. 1985–2002, 2007.
[9] G. Tryggvason, R. Scardovelli, and S. Zaleski, Direct Numerical Simulationsof Gas-Liquid Multiphase Flows, 1st Editio. Cambridge University Press, 2011.
[10] R. Peyret and T. D. Taylor, Computational Methods for Fluid Flow. Berlin,
Heidelberg: Springer Berlin Heidelberg, 1983.
[11] M. M. Francois, S. J. Cummins, E. D. Dendy, D. B. Kothe, J. M. Sicilian, andM. W. Williams, “A balanced-force algorithm for continuous and sharp interfacialsurface tension models within a volume tracking framework,” J. Comput. Phys., vol.
213, no. 1, pp. 141–173, Mar. 2006.
[12] J. . Brackbill, D. . Kothe, and C. Zemach, “A continuum method for modelingsurface tension,” J. Comput. Phys., vol. 100, no. 2, pp. 335–354, Jun. 1992.
[13] T. Yabe, K. Chinda, and T. Hiraishi, “Computation of surface tension andcontact angle and its application to water strider,” Comput. Fluids, vol. 36, no.
1, pp. 184–190, Jan. 2007.
[14] M. Sakai et al., “Lagrangian–Lagrangian modeling for a solid–liquid flow ina cylindrical tank,” Chem. Eng. J., vol. 200–202, no. null, pp. 663–672, Aug. 2012.
[15] T. Kajishima and S. Takiguchi, “Interaction between particle clusters andparticle-induced turbulence,” Int. J. Heat Fluid Flow, vol. 23, no. 5, pp. 639–
646, Oct. 2002.
[16] B. Lafaurie, C. Nardone, R. Scardovelli, S. Zaleski, and G. Zanetti,
3.1-8

“Modelling Merging and Fragmentation in Multiphase Flows with SURFER,” J. Comput.

Phys., vol. 113, no. 1, pp. 134–147, Jul. 1994.
[17] M. Uhlmann, “An immersed boundary method with direct forcing for thesimulation of particulate flows,” J. Comput. Phys., vol. 209, no. 2, pp. 448–476,
Nov. 2005.
[18] T. R. Auton, J. C. R. Hunt, and M. Prud’Homme, “The force exerted on a bodyin inviscid unsteady non-uniform rotational flow,” J. Fluid Mech., vol. 197, pp.
241–257, Apr. 1988.
[19] M. Sussman, “An Adaptive Mesh Algorithm for Free Surface Flows in General
Geometries,” in Adaptive Method of Lines, A. Vande Wouwer, P. Saucez, and W.
Schiesser, Eds. Chapman and Hall/{CRC}, 2001.
[20] X. Sun and M. Sakai, “Numerical simulation of two-phase flows in complex
geometries by using the volume-of-fluid/immersed-boundary method,” Chem. Eng. Sci.,
vol. 139, pp. 221–240, Jan. 2016.
[21] J. N. Israelachvili, Intermolecular and Surface Forces: Revised Third Edition,
3rd editio. Academic Press, 2011.
[22] Y. Tsunazawa, D. Fujihashi, S. Fukui, M. Sakai, and C. Tokoro, “Contact
force model including the liquid-bridge force for wet-particle simulation using thediscrete element method,” Adv. Powder Technol., vol. 27, no. 2, pp. 652–660, Feb.
2016.
[23] K. A. Brakke, “The Surface Evolver,” Exp. Math., vol. 1, no. 2, pp. 141–
165, Jan. 1992.
[24] D. J. Broesch and J. Frechette, “From concave to convex: capillary bridges
in slit pore geometry.,” Langmuir, vol. 28, no. 44, pp. 15548–54, Nov. 2012.
[25] Y. I. Rabinovich, M. S. Esayanur, and B. M. Moudgil, “Capillary forces
between two spheres with a fixed volume liquid bridge: theory and experiment.,”

Langmuir, vol. 21, no. 24, pp. 10992–7, Nov. 2005.
[26] P. Lambert, A. Chau, A. Delchambre, and S. Régnier, “Comparison between two
capillary forces models.,” Langmuir, vol. 24, no. 7, pp. 3157–63, Apr. 2008.
[27] K. Washino, H. S. Tan, A. D. Salman, and M. J. Hounslow, “Direct numericalsimulation of solid–liquid–gas three-phase flow: Fluid–solid interaction,” Powder

Technol., vol. 206, no. 1–2, pp. 161–169, Jan. 2011.
[28] K. Yokoi, D. Vadillo, J. Hinch, and I. Hutchings, “Numerical studies of theinfluence of the dynamic contact angle on a droplet impacting on a dry surface,”

Phys. Fluids, vol. 21, no. 7, p. 72102, 2009.
[29] A. S. Joshi and Y. Sun, “Wetting dynamics and particle deposition for anevaporating colloidal drop: A lattice Boltzmann study,” Phys. Rev. E, vol. 82, no.
4, p. 41401, Oct. 2010.

3.1-9

3.1.2 固気液混相流-臨界連成手法の開発（平成 29・30 年度）(1) 平成 29 年度①-1 シミュレーション解析
1) 解析条件
固気液三相流れの数値シミュレーションと連続エネルギーモンテカルロ法コードを接続して、切り粉の堆積状態と臨界安全との関係についても検討した。計算体系は前節の崩壊水柱の多孔質
壁透過問題である。ただし、水柱の初期高さが 0.40m、多孔質材の平均粒径が 1.0mm に設定した。
三つの計算ケースには、多孔質壁の空隙率をそれぞれ 0.9、0.7、0.5 に設定した（表 3.1-2）。
2) 解析結果・考察
図 3.1-7 に計算結果を示す。Case 1 では porosity が 0.9 であり、多孔質壁で流体が減速されず壁面を通過する様子が模擬された。Case 2 では porosity が 0.7 で、Case 1 と比較して流体の速度が減衰されている様子が示された。本体系では、多孔質体における固液相互作用力の計算に用いられる代表粒径は全て 1.0 mm と同じ値を用いているため、空隙率によって流体が減速され、挙動が変わることが示唆された。Case 3 では porosity が 0.5 で、流体速度の減衰が大きくなる
ことが示された。
これらの結果を基に、MCNP を用いて実効増倍率 keff の経時変化を計算した。表 3.1-3 に、MCNPのインプットに用いた核種の原子数密度を示す。多孔質壁(固相)は 4.4%濃縮ウラン燃料、気相は空気、液相は軽水とした。MCNP セルにおける各原子の粒子数密度は、VOF 関数の値と多孔質体の空隙率を基に計算した。keff の経時変化をプロットした結果が図 3.1-8 である。Porosity の値が大きくなるにつれて、keff の値が小さくなっている。今回の計算体系では、中性子の平均自由行程に対して計算領域が狭く、中性子がすぐに計算領域の外へ出てしまうために、水の通過に伴って中性子の増加は起こらないことが示唆された。
以上の結果から、MCNP と DEM-VOF 法を連成させることで、連続的な keff の計算および臨界評価
手法の開発に成功した。デブリ取り出し時のリスク評価を連続的に行えるため、安全性の向上に資する手法であることが示された。

3.1-10表 3.1-2 計算条件
Case 1Case 2

Case 3

Domain size

0.89 m × 0.58 m

Cell size

5 mm
Water (1000 kg/m3, 0.001 Pa s)

Liquid phase

Air (1 kg/m3, 10-5 Pa s)

Gas phase
Water dam: Height

0.40 mSolid phase
Porous media: Porosity

4.4wt% enriched Uranium fuel
0.9

0.7

Porous media: Grain size

0.5

1.0 mm
表 3.1-3 原子数密度の条件

Material

4.4wt% enriched Uranium fuelAir

NuclideAtomic density [atoms/b/cm]

235

U

1.03467E-3

238

U

2.21965E-2

16O

4.64624E-2

14N

4.39292e-5

16O

1.02232e-5

H

6.67382e-2

16

3.33698e-2

1

Light waterO

3.1-11Case 1Case 2
図 3.1-7 固気液混相流計算結果の比較

3.1-12

Case 3

0.01
50%

70%

90%

0.008

keff

0.006

0.004

0.002

0

0

0.5

1

1.5

2

2.5

Time [s]図 3.1-8

keff の経時変化

3.1-133

3.5

4

(2) 平成 30 年度①-1 シミュレーション解析
1) 液架橋の幾何学的形状
図 3.1-9 (a)に、粒子-粒子間に構成される液架橋の幾何学的形状を示す。 、 、 、 、は、それぞれ粒子半径、粒子表面距離、偏角、接触角、接触半径を表しており、下付きの 1 および 2 はそれぞれ粒子 1、2 を示している。図 3.1-9 (b)は、粒子-壁面間に構成される液架橋を示しており、右の粒子が固体平面に置き換えられている。粒子-壁面間の液架橋に対しては、偏角を定義することはできないが、接触半径 は粒子-粒子間の液架橋と同様に定めることができる。
重力効果が支配的である場合を除き、静的液架橋の界面は一定の曲率を有している。軸対称な液架橋では、曲率 は 2 つの曲率半径（半径方向の内曲率半径率半径）を用いて次式で表される。
=

1

−

1

と経線方向の外曲(1)

次に、液架橋力の具体的な計算方法について記述する。粒子-粒子間、粒子-壁面間にはたらく液架橋力は、ラプラス圧と表面張力の和で表される
=−

sin

+ 2!

sin sin" + #(2)

ここで、 はラプラス圧、!は表面張力を示す。また、液架橋内におけるラプラス圧は表
面張力と曲率から求めることができる。

= σκ

(3)

(a) Axisymmetric liquid bridge(b) Axisymmetric liquid bridge

between two spheres.

between sphere and plane.

図 3.1-9 液架橋の幾何学的形状
2) 既存の液架橋力モデル
液架橋力の分析モデルはこれまでに数多く提案されており、Israelachvili（2010)やRabinovich et al,（2005）のモデルは、その単純さから複数の研究で使用されている(1)(2)。しかし、これらの分析モデルが妥当な解を導き出すのは、液量、粒子表面距離、接触角が比較的小さい場合に限定されている。他にも、厳密解として得られた点データをフィッティングすることで、モデル方程式を導出するという手法があるが、元々のデータへの依存性が強いため、一般条件下に対して正常に拡張されない可能性が考えられる。

3.1-14

幾何学的近似によって液架橋の構造を単純化することは、方程式の導出を容易にする上で有効性が高く、中でもトロイダル近似は液体橋モデリングによく使用されている手法である。Huppmann and Riegger（1975）の液架橋力モデルは、トロイダル近似に基づいた、おそらく最初の液架橋力モデルであり、妥当な精度と汎用性をもって液架橋力を予測できる(3)。トロイダル近似を用いることで、軸対称な液架橋のメリディオナル形状は円形にフィッティングされる。これにより、液架橋力の近似解を求めるにあたり、更なる幾何学的関係を導出することが可能となる。
Huppmann-Riegger モデルでは、液架橋の形状に対してトロイダル近似を行うことにより、偏角 と曲率半径、の間に次の関係が成り立っている。

=

sin

=

"1 − sin" + ##−"1 − cos # + /2
cos" + #更に偏角 は、液量)* に対して次の関係式を持っている。
)* = 2π,cos" + # − " /2 −− #-" .

#+

+

(4)

cos" + #(5)

液架橋力を求めるためには、上の陰関数方程式を解き、偏角 を求めなければならない。
その複雑さにも関わらず、Huppmann-Riegger モデルは広範囲条件の液架橋に適用可能である。しかし、Huppmann-Riegger モデルは、同じパラメータの粒子間において発生する
液架橋について記述したものであるため、図 3.1-9 のように粒子径や接触角といった粒子のパラメータが異なる場合は、Derjaguin 近似に基づいた有効方程式を用いることで対応
1 1
1⎧ 4 + 6
1
2 5
=
1
⎨
/00
⎩
2 5

できる。

sphere − spheresphere − plane1
cos /00 = "cos 5 + cos
2

(6)

#

3) 液架橋力の数値解析手法
本研究では、Sun and Sakai(2018)による液架橋力モデルを使用した(4)。この液架橋力
モデルはトロイダル近似を使用しており、広範囲の液量および接触角に対して適用可能である。加えて、異なる粒子半径と接触角を有する液架橋に対しても直接適用できるため、非常に汎用性が高いモデルである。
図 3.1-10 は、粒子-粒子間の液架橋モデルを表している。粒子半径 5 、、粒子表面距

離 、接触角 5 、 、液量)* については既に値が定められていると仮定すると、変数は偏
角 5、

と曲率半径

、

である。これらの変数に対して方程式を導出し、変数を 5のみに減らすことで、最終的に液架橋力を算出する。
曲率半径、は、次式のように偏角 5 、

3.1-15で表すことができる。

5 "1 − >?@ 5 # +"1 − >?@ # +>?@" 5 + 5 # + >?@" + #

=

=

5 sin

=

5−

sin

−

式(8)、式(9)より、偏角 5 、
−

5 sin

5−

sin
sin" 5 + 5 # − sin" +これらの式を解くことで、
−

5 sin

5−

sin
sin" 5 + 5 # − sin" +上記の方程式より、変数"1 − sin" 5 +"1 − sin"

+

5 ##

##

(7)
(8)
(9)

の間に次の関係式が得られる。

は

#
#

、

5 "1 − cos 5 # +"1 − cos # +cos" 5 + 5 # + cos" + #

=

(10)

の関数として計算できる。
5 "1 − cos 5 # +"1 − cos # +cos" 5 + 5 # + cos" + #

=

(11)

はすべて 5 の関数として計算することが可能と、なる。最後に、液量)* が一定であることを用いて方程式を導出する。液架橋の境界面を回転させることによって得られる体積)= は、
)=

=

""

+
−

−

# +

+

"
"

+
+

. ">?@ .

3#">?@ ∅5 + >?@ ∅ ##"∅5 + ∅ − ##"@BC ∅5 >?@ ∅5 + @BC ∅ >?@ ∅ #∅5 + >?@ . ∅ #

∅5 =

∅ =

5+

+

5

したがって、液量)* は次式を満たしていなければならない。
)* = )" 5 ,,# = )=,(12)− )5 − )(13)ここで、)5 と) は、液体中に浸された粒子の球冠部分の体積であり、それぞれ次のように求められる。

)5 =
)5 =

3
3

.
5 "2 − 3 cos
5 "2 − 3 cos.

5 + >?@
5 + >?@

.
.

5#

(14)

5#

液架橋が対称性を有しているとき、式(12)には Huppmann-Riegger モデルの全ての項が含まれており、このモデルの妥当性を裏付けている。
粒子-壁面間の液架橋モデルについては、右側の偏角
ここでは

= 0、

と球冠体積) が存在しないため、= 0とおく。式(11)は除外されることになるが、式(7)、式(8)、式(13)の関係は維持されているため、液架橋力の計算が可能である。

液架橋の変曲点の位置によって、液架橋力の計算方法は少し異なる。図 3.1-10 (a)に示される液架橋では、変曲点は内曲率半径に対応して実際の液架橋上に存続する。その変曲点上において、液架橋力はラプラス圧の定義式（3）を用いて次のように計算される。
=−

= !

+2 !

"1 +

/

#

(15)

例外的ではあるが、図 3.1-10 (b)のように、変曲点の位置が粒子内部に位置する場合が3.1-16

存在する。このような液架橋は、粒子表面の濡れ性が大きく異なる場合などに、実際に構成される。したがって、この液架橋力モデルでは、より汎用性を高めるべく、この特殊な場合に対しても液架橋力の計算式を拡張している。円形表面の中心は、
EF =

5 cos

5+

cos" 5 +5#(16)

これは変曲点の x 座標と等値である。したがって、次の条件が成立する。
5 cos

5 G EF G

5++"1 − cos#(17)

この条件が満たされない場合は、実際の境界面に変曲点が存在しないことを意味する。
この場合、内曲率半径を用いている式（15）を用いることは適切ではない。ここでは、は以下の式を使用して推定することができる。

5 =

IIII
H =
5

sin" 5 +

1
" +
2 5

5#

,

=

#

sin"

+

(18)

#

ここで、 5 および は固体粒子上の接触円の半径である。

式（1）の

=

5 =

sin" 5 +

5 sin

5

5# −

sin"

+

#

(19)

をIIIIと置き換えることで、式（3）よりラプラス圧を求めることができHる。そして、式（2）より、2 つの接触部における液架橋力がそれぞれ計算される。
5 =−

=−

5

+ 2!

+ 2!

5 sin" 5 +

sin"

+

5#

#

(20)

これらの平均から、変曲点の存在しない液架橋における液架橋力が求められる。
=

1
" +
2 5

#

(21)

この処理は、多くの接触角の組み合わせを有する capillary 液架橋にとっては重要なものとなり得る。また、粒子壁面間の場合についても、
上記の式を変形することなく計算することができる。

(a) Toroidal approximationof a liquid bridge.

= 0、= 0と定めることで、(b) Liquid bridge without
real neck point.

図 3.1-10 粒子-粒子間の液架橋モデル

3.1-17

4) モデル応用
1. 計算条件
燃料デブリ取り出し時には、湿潤粉体の挙動を予測できるモデルが必要となる。本研究では、DEM における液架橋力モデルを用いて、湿潤粉体の挙動が定性的に模擬可能となることを示す。
表 3.1-4 に物性値を示す。固相にはガラスビーズの物性を用いており、密度は 2,500kg/m3 、バネ定数は 1,000 N/m、反発係数は 0.9、摩擦係数は 0.3、表面張力係数は72.75x10-3 N/m、粒子の直径は 1.8 mm とした。固体粒子の接触角は 60°、壁面の接触角
は 100°とし、粒子壁面間の有効接触角は次式に基づき 80.6°とした。
1
cos /00 = ,cos J + cos K 2(22)

表 3.1-5 に計算条件を示す。計算粒子数は 70,000、時間刻みは 2.5x10-5 sec、液量は粒子総体積の 5%で 1.07x10-5 m3 とした。
図 3.1-11 に、シミュレーションに用いた回転ドラムを示す。ドラムの直径は 120 mm、奥行きは 100 mm で、回転速度は 30 rpm とした。
2. 結果・考察
図 3.1-12 に、計算結果を示す。図 3.1-12 (b)は 2.0 sec(1 回転)時の粒子配置を示して
おり、液架橋によってドラム側面に粒子が付着する挙動が確認できた。また、粒子層は楕円形状を形成し、ドラムとの間にはギャップが生じていることが図 3.1-12 (b)-(d)で確認できた。これは、粒子間の液架橋が引力としてはたらいた結果、凝集体を形成したためと考えられる。その後、さらにドラムを回転させ、図 3.1-12 (e)のように 15 回転させたあたりから、ドラムと粒子層の両端との間のギャップがほとんどなくなり、平坦な表面を形成することがシミュレーションによって模擬された。このとき、ドラム側面に付着した粒子は図 3.1-12 (b)-(d)と比較して数が多くなっており、粒子層の両端がドラム側面側へと引っ張られたため、ギャップがなくなったと考えられる。
図 3.1-13 に、固体粒子の速度分布を示す。図 3.1-13 (b)-(d)では、壁面に付着した粒子の速度は凝集体下部のところで他の部分よりも遅いことが示された。本モデルでは、配位数に応じて液量を等分しているため、粒子-粒子間の液架橋力の方が粒子-壁面間の液架橋力よりも支配的となり、その結果として、凝集体下部で速度がドラムの回転速度よりも遅くなったと考えられる。時間の経過とともに、壁面に付着する粒子層が厚くなってくると、凝集体下部における粒子の速度はドラムの回転速度に近づいており、図3.1-13 (f)に示すように、30 回転時には、ドラム側面付近の固体粒子の速度はほぼ一様となった。回転数が少ない時は凝集体として回転していたため、粒子層の表面において粒子の速度は 0.1 m/s 程度とドラムの回転速度と比べて遅かったが、回転数が増加するに従い、ドラムによって持ち上げられた粒子が雪崩のように粒子層の表面を流れ落ちる挙動へと遷移しており、その結果として、粒子層表面における固体粒子の速度もドラムの回転速度とほぼ同等の 0.18 m/s 程度まで上昇した。以上の結果から、湿潤粉体の挙動をDEM によって定性的に模擬可能であることが示唆された。
3.1-18

表 3.1-4 物性値
Item

Unit

Value

Density

kg/m3

2,500

Spring constant

N/m

1,000

Restitution coefficient-

0.9

Friction coefficient-

0.3

Surface tension coefficient

N/m

72.75x10-3

Contact angle (particle-particle)

deg

60

Contact angle (particle-wall)

deg

80.6

表 3.1-5 計算条件
Item

Unit

Value

Time stepsec

2.5x10-5

Number of particles-

70,000Liquid volume

m3

1.07x10-5

図 3.1-11 回転ドラム装置図

3.1-19(a)0.0 sec(b)2.0 sec(0 revolution)(1 revolution)(c) 6.0 sec(d) 10.0 sec(3 revolutions)(5 revolutions)(e) 30.0 sec(f) 60.0 sec(15 revolutions)(30 revolutions)

図 3.1-12 湿潤粉体の挙動予測モデルによる計算結果

3.1-20(a) 0.0 sec(b) 2.0 sec(0 revolution)(1 revolution)(c) 6.0 sec(d) 10.0 sec(3 revolutions)(5 revolutions)(e) 30.0 sec(f) 60.0 sec(15 revolutions)(30 revolutions)
図 3.1-13 粒子速度分布3.1-21

3.1.3 デブリ冷却性能および切り粉漏洩評価（平成 29・30 年度）(1) 平成 29 年度
多孔質構造の内部へ浸入した液面を完全に再現することが不可能であるため、既存手法を適用できない。実際に多孔質材内部において、細孔構造（Pore structure）によって、連続的な界面が存在できず、固体表面に付着した微小界面に分断される。前述の VOF 関数の全局的界面曲率がもはや無意味となり、局所曲率を真の値に置き換えなければならない。この場合、空孔中の浸漬液面に対するモデリングが必須となる。本手法における計算格子スケールに比べて、多孔質材の空孔サイズがより小さいので、これは界面曲率の Sub-grid Scale モデルに相当する。
ここで、球状粒子充填によって形成された多孔質固体に限定して議論する。図 3.1-14 に示されたように、液面が距離 h と接触角 θ の固体表面を浸潤し、ある位置で（重力などの外部作用と）つりあい状態に達した。界面の形状が半径 r の円としてよく近似できるとし、その曲率が半径の逆数より概算できる。

≈ M N5(23)ℎ = 2M ⋅ >?@(24)

さらに、球体粒子に比べて、空隙のサイズが十分に小さいと仮定すれば、単純な幾何学的関係か
ら次式が導かれる。

つまり、多孔質内部空孔に付着した界面の曲率は次式のように与えられる。
≈

2>?@
ℎ

(25)

このように、界面曲率の計算問題は、多孔質内部の代表的な空孔径の見積もりと同等であることが分かった。
連結的かつ均質的な多孔質材料に対して、その空孔径 h が空隙率 ε の関数となるのは自然であろう。特に、単純球充填による多孔質固体の場合、空孔径と粒子径の比 h/dp は空隙率に関して
単調増加挙動を示すべきである。その漸近的極限として、ε=0 でゼロとなり（空隙なし）、ε=1
で無限大となる（完全開放）。理論的・実験的にその関係式を導出できれば、固体内部の界面曲率を計算することが可能となる。ただし、この問題に対するアプローチは本研究の取り組みの対象外であり、実際の多孔質デブリ構造を解明する将来の課題に帰結すべきだと考えられる。その代わりに、簡便性のため、球充填多孔質の空孔‐粒子比を妥当に直接設定し、この修正表面張力
計算手法を本数値解析フレームワークに導入することによって、固体内部の毛細管効果を模擬できることを示す。
計算例として、液滴と多孔質平面の衝突問題を解析した。本計算における流体の物性値は、前節で取り上げた低速液滴‐粒子衝突問題と一致した。詳細は表 3.1-6 に参照されたい。液滴直径が 1 mm で、初期速度が 0.1 m/s である。また、水平面に設置した多孔質材は、均一的な空隙率0.5 を有し、直径 dp = 0.1 mm の球状粒子によって構成されるものとする。この空隙率の下で、単純立方格子構造を考慮した粒子間の最小距離が約 0.05dp、最大距離が約 0.8dp と推定できる。
そのため、多孔質固体の相当空孔径をその中間値に近い h = 0.33dp と設定した。現段階において、この設定にはある程度の任意性が存在するが、実測定値に基づいたモデルを用いて定式化す
ることが可能である。また、接触角をゼロ度（完全濡れ）と与えた。これらの多孔質平面の物性値を表 3.1-7 にまとめた。
3.1-22

計算結果の代表的なスナップショットを図 3.1-15 に示す。液滴表面にある赤い線は液滴と多孔質が交差する境界を示す。このケースでは、外部液滴より多孔質体内部のミクロ的な曲率が遥かに大きいため、強い負圧と吸引力が発生した。多孔質領域を通過する際の抵抗力に比べて、この毛管作用が支配的となり、液滴全体が迅速に多孔質体内部へ吸い込まれた（＜ 5 ms）。濡れ性の高い液体の浸潤挙動を定性的によく再現したことが分かった。
表 3.1-6 切り粉漏洩評価の計算条件
Case 1 (low speed)

Case 2 (high speed)

Domain size

1.6 mm × 1.2 mm × 1.6 mm

Cell size

12.5 μm

Liquid density

1000 kg/m3

Liquid viscosity

0.001 Pa s

Gas density

10 kg/m3

Gas viscosity

10-5 Pa s

Surface tension

0.03 N/mContact angle (floor)

90 degree

Contact angle (particle)

10 degree

Particle density

2500 kg/m3

Particle diameter

100 μmDroplet diameter1000 μm

Time step

0.1 μs0.01

Droplet velocity

0.1 m/s

4.0

Re number

100

4000

We number

0.333

533

3.1-23図 3.1-14 球状粒子空隙中の浸潤液面
表 3.1-7 多孔質体の物性値
Parameters
Porosity

0.5

Contact angle

0°Particle diameter

100 μmEquivalent gap distance

33.3 μm

3.1-24図 3.1-15 液滴と多孔質平板の衝突および浸入（VOF 関数が 0.5 となる等値面を黄色で示した）

3.1-25

(2) 平成 30 年度
1) 炉内の気流評価
1. 解析条件計算に用いたガスの物性値は、密度：1.0 kg/m3、粘度：1.8×10-5 Pa・s、比熱：1000 J/(kg・K)、熱伝導率：0.024 W/(m・K)、温度：27 ℃とした。ガスの流量は 7.85×10-4 m3/s（流出部径0.1 m、流出速度 0.1 m/s 相当）とした。また熱源を考慮するケースでは、装置底部の円盤部表面よりガスへ発熱量 35 kW を与えた。円盤部以外の装置底面において、環境温度: 27 ℃、熱伝導係数: 167 W/(m2・K)の境界条件にてガスと壁間の熱交換を考慮する。ガスの物性値および計算条件をそれぞれ表 3.1-8 および表 3.1-9 に示す。
2. 解析結果
以下に、粒子が含まれないガスの挙動について示す。熱源の有無で区別されたケース（102.96
秒経過後）を示す。気流の流速ベクトル（図 3.1-16）に関して、熱源なしのケースでは気流の上方移動が見られなかったのに対して、熱源ありのケースでは気流が垂直上方に移動した。また上昇した気流はその後、側壁に沿って下方に気流が移動したことも示された。熱源ありのケースに関して、炉内の温度分布を図 3.1-17 に示す。最高温度は 68.5 ℃だった。

表 3.1-8 ガスの物性値
Item

Unit

Value

Density

kg/m3

1.0

Viscosity

Pa・s

1.8 x 10-5Specific heatJ/kg・K

1,000

Thermal conductivity

W/m・K

0.024

Temperature℃

27

Flow ratem3/s

7.85×10-4

表 3.1-9 装置の物性値
Item

Unit

Value

Gas outlet diameterm

0.1

Thermal conductivity

W/m2・K167

3.1-26

熱源なし

熱源あり

図 3.1-16 気流の流速ベクトル（固体粒子なし）
矢印の向きと長さは表示断面（鉛直断面）上の速度の大きさおよび方向を示し、矢印の色は断面に垂直な成分も含めた速度絶対値を示す。

熱源あり
図 3.1-17 ガス温度（固体粒子なし）
2) 炉内の切り粉漏洩評価
1. 解析条件
今度はガス内に粒子を入れて解析を行った。ガスおよび計算領域に関する物性値・計算条件は、粒子なしのケース（熱源あり）と同じくそれぞれ表 3.1-8 および表 3.1-9 の通りである。粒子の物性値は、密度: 10,970 kg/m3、実粒子径: 0.025 mm、代表粒子径: 0.003 mm、反発係数: 0.9、
摩擦係数: 0.3 である（表 3.1-10）。
2. 解析結果[実粒径: 1.0 µm または 10.0 µm]、[熱源位置: 中央または壁面近傍]の組み合わせである 4ケースならびに[実粒径 25.0 µm、熱源位置: 壁近傍]の 1 ケース（計 5 ケース）を計算した。以下に、これらのケースに関して、図 3.1-18 に炉内におけるガス流速ベクトル、図 3.1-19 に温度

3.1-27

分布、図 3.1-20 に粒子の速度ベクトルを示す。
ガス流速ベクトルに関して、熱源位置が中央にあるケースの方が壁面近傍のケースよりも大きな速度が炉中央付近で観察された。粒子のサイズ依存性は見られなかった。
温度分布に関して、熱源の位置に応じて偏りが見られ、中央にあるケースでは対称な分布が得られた。粒子サイズ依存性は見られなかった。最高温度は全ケースで 67〜68 ℃であった。
粒子の速度ベクトルに関しては、ガス流速ベクトル、温度分布に比べ、粒子径 1.0 µm のケー
スと 10 µm および 25 µm のケースで異なる傾向が見られた。熱源が壁面近傍の場合、粒子径が10.0 µm 以上では最高速度が等しく(0.8 m/s)、粒子径 1.0 µm ではより低い最高速度が見られた（0.75 m/s）。最小速度については、1.0 µm のケースで 7.1 x 10-5 m/s であったのに対して、10.0 µm および 25.0 µm のケースでは 1.3 x 10-7 m/s だったことから、100 倍程度の違いが見られた。小さい粒子（1.0 µm）の方が全体的に吹き上がる傾向が見られ、炉底で滞留している粒子も静止せず低い位置を維持したまま移動を続けていることが分かる。他方、大きい粒子（10.0µm および 25.0 µm）は上方に吹き上がる粒子の割合は小粒子に比べると少なく、炉底で沈降して静止している粒子が多数見られる。また熱源が中央付近にある場合でも、粒子速度の最小値は粒子径 1.0 µm の方が粒子径 10.0 µm のケースよりも 9 倍程度大きかった。これは、粒子の流体に対する追従性の違いで説明できる。小粒子は追従性が高く、炉内全体におよぶガス移動に合わせてほとんどの粒子が移動を続けるのに対して、大粒子は追従性がそれほど高くないために、比較的流速が高い熱源付近のガス移動だけに巻き込まれ、そこから離れた箇所での流体の移動に対しては引きずられなくなるからである。
表 3.1-10 粒子の物性値
Item

Unit

Value

Density

kg/m3

10970

Actual diameterµm

25.0 / 10.0 / 1.0

Representative diameter

µm

3.0

Restitution coefficient-

0.9

Friction coefficient-

0.3

3.1-28熱源位置: 中央

熱源位置: 壁近傍
実粒径: 1.0 µm

熱源位置: 中央

熱源位置: 壁近傍
実粒径: 10.0 µm

熱源位置: 壁近傍、実粒径: 25.0 µm
図 3.1-18 気流の流速ベクトル（固体粒子あり）
矢印の向きと長さは表示断面（鉛直断面）上の速度の大きさおよび方向を示し、矢印の色は断面に垂直な成分も含めた速度絶対値を示す。

3.1-29熱源位置: 中央

熱源位置: 壁近傍
実粒径: 1.0 µm

熱源位置: 中央

熱源位置: 壁近傍
実粒径: 10.0 µm

熱源位置: 壁近傍、実粒径: 25.0 µm
図 3.1-19 ガス温度（固体粒子あり）

3.1-30熱源位置: 中央

熱源位置: 壁近傍
実粒径: 1.0 µm

熱源位置: 中央

熱源位置: 壁近傍
実粒径: 10.0 µm

熱源位置: 壁近傍、実粒径: 25.0 µm
図 3.1-20 粒子速度[m/s]3.1-31

参考文献(1) J. N. Israelachvili, “Intermolecular and Surface Forces,” Academic Press,
2010.
(2) Y. I. Rabinovich, M. S. Esayanur, B. M. Moudgil, “Capillary forces between
twospheres with a fixed volume liquid bridge: theory and experiment,”
Langmuir 21, 10992–10997, 2005.
(3) W. Huppmann, H. Riegger, “Modelling of rearrangement processes in liquid
phase sintering,” Acta Metall, 23, 965–971, 1975.
(4) X. Sun, M. Sakai, “A liquid bridge model for spherical particles applicable
to asymmetric configurations” Chem. Eng. Sci, 182, 28-43, 2018.

3.1-32

リスク評価
3.2.1 技術リスク（平成 28 年度）
燃料デブリ取り出しの手順が決まった場合のそのプロセスの技術リスクを見積もるための簡易システムを構築するに際しての基本的な方針を検討する。
さまざまな活動に対するリスク管理は、その活動を円滑に合理的に遂行し、目的を達成するために重要なプロセスである。図 3.2-1 に、国際リスクガバナンス協議会[1]の示すリスクガバナ
ンスの枠組みを示す。リスクガバナンスとは、リスク評価とリスク管理を継続的に実施するための構造と意志を示すものであり、それをコミュニケーションによりアウトリーチする仕組みである。

図 3.2-1 リスクガバナンスの枠組み
リスク管理においては理に適う意思決定を行うためにリスク評価が行われ、それの適切性が示されることが大切である。よって、リスク管理では、リスクの事前評価を行うとともに、リスクを評定する（定量的にリスクを求め、それを分析する）ことが必要である。以下では、リスクの事前評価に相当する分析について述べる。
福島第一原子力発電所の燃料デブリの取出しに係る方針を決める時期にある。その作業は、以下の手順（図 3.2-2）に従うこととされている。

図 3.2-2 燃料デブリ取出しのエンドポイントと枠組み
廃炉作業を進めていく過程で留意すべき放射線リスク管理上の課題は、①再臨界防止 ②冷却機能確保 ③放射性物質の漏洩防止 ④水素爆発防止 ⑤原子炉建屋の耐久性確保 ⑥作業員の安全3.2-1

確保の６点であろう。それらに焦点を当て、燃料デブリ取出しに関するリスク管理の基本的な方針を検討する。
国際リスクガバナンス協会によれば、リスクとは私たちが価値を認めるものに関する事象や活動がもたらす不確かな影響である。福島第一原子力発電所の廃止措置等において、放射線安全、
労働安全、環境の保全、プロジェクト進捗が達成すべき使命であり、私たちの考える価値である。
それの円滑な遂行・達成への影響がリスクであり、これらがリスク管理の対象となる。このうち、
本研究では、技術リスクすなわち放射線リスクを対象とする。
リスクとは、どのようなことが起きうるのか（シナリオ）、それはどれくらい現実的なのか（確からしさ）、その結果どのような被害あるいは影響があるのか（影響度）の三要素からなる。
定量的リスク評価のプロセスは６ステップからなる。これは、あらゆるリスクを伴う行為のリス
クマネジメントに適用できる。
ステップ 1：対象システムを明確に定義する（成功するプロセスを記述すること）
ステップ 2：成功するプロセスに影響を及ぼしうる危険要因（hazards）を列挙し特性化する（危険要因がどのような影響を与えるか）ステップ 3：成功するプロセスが失敗に至るあらゆるシナリオを策定するとともに、それらをカテゴリー化する（リスクトリプレットの第一要素）
ステップ 4：各シナリオの尤度を定量化する（そのシナリオの確からしさでありリスクトリプ
レットの第二要素）ステップ 5：各シナリオの被害度を定量化するとともに、適切に分類してリスクプロファイルを記述する。（リスクトリプレットの第三要素）
ステップ 6：リスクマネジメントの指針となるべく結果の分析と解釈を行い、これによりリスク管理に関する意思決定（対策）の優先度を定める
本研究においてもこの枠組みが適用できると考えられる。
図 3.2-3 と図 3.2-4 に炉内安全確保と燃料デブリ切断（図 3.2-2 参照）に関する作業とリスク要因、それに対する対応策を示す。これをリスク構造と呼んでいる。

図 3.2-3 炉内安全確保に関するリスク構造

図 3.2-4 燃料デブリ切断に関するリスク構造3.2-2

以上の通り、燃料デブリ取出しに係るエンドポイントと枠組み、一般的に適用可能なリスク管理プロセス、燃料デブリ取出しの二つのタスクを対象としたリスク構造の明示を実施した。これをさらに具体化、定量化することにより燃料デブリの取出しのための簡易システムを構築できる。
参考文献[1] An introduction to the IRGC Risk Governance Framework, International Risk
Governance Council, 2012 ISBN 978-2-9700772-2-0
3.2.2 技術リスク、プロジェクト遂行リスク、社会的意思決定リスク（平成 29 年度）
燃料デブリ取り出しに係る戦略の具体化を合理的に進めるためには、技術的リスクのみでなく、プロジェクトの進捗に係るリスク、社会の理解醸成に係るリスクは注視すべき項目である。これらリスクを燃料デブリ取り出しに係る意思決定の中にどのように取り組むべきかを検討し、リスク指標として整理した。この成果により、燃料デブリ取り出しの技術リスク、プロジェクト遂行リスク、社会的意思決定リスクをバランス良く評価するための考え方を提示する。
技術戦略プラン 2017[文献(7)]では、廃止措置の目標は、放射線リスクを速継続的にかつ速やかに低減させることとしている。また、5 つの基本的考え方として、(1)安全（放射性物質に起因するリスクの低減及び労働安全の確保）、(2)確実（信頼性が高く、柔軟性のある技
術）、(3)合理的（リソース(ヒト、モノ、カネ、スペース等)の有効活用）、(4)迅速（時間軸の意識）、(5)現場指向（徹底した三現(現場、現物、現実)主義）を指摘している。リスク評価に当たっては、この目的ならびに基本方針を念頭におくことが適切である。
廃止措置においては、いかに最善と考えられる方策をとったとしても残留リスクが存在する。原子力発電所の運転において残留リスクは発電という便益と対比しつつその正当性が議論される。操作や保全に関する手順を一時的にでも変更する場合、その変更による一時的なリスクの増加は、実質的な安全が確保されているならば許容することによって円滑な廃止措置の進捗がなされる。福島第一原子力発電所の廃炉においては、原子力発電所の運転の場合と異なり、発電のような積極的な便益はない。しかし、現存するリスクを低減することが、総体としての安全とセキュリティにつながる。したがって、長期的視点をもって着実にリスクが低減されるならば、公益につながる。廃炉は技術的なチャレンジを伴う作業であり、エンドポイントを明確に定義し、長期的観点からリスクを制御しマネジメントすることが重要である。
リスク評価と管理の戦略として、6 つのステップを提案する。これを図 3.2-5 に示す。ステップ 1 では、エンドポイントに到達するためのシステム要素を定義する。システムとはエンドポイントへの成功パスであり、それを構成する要素はそれに必要なタスクである。ステップ 2
は、成功パスの構成要素に対する脅威や障害を同定する。タスクと脅威あるいは障害に組み合わせはひとつの起因事象を定義する。起因事象の全リストを提示し、それらの尤度を概括評価する作業はステップ 3 である。ここまでで、システムに関する情報が必要となり、これを用いて起因事象の進展が与えられる。これは成功パスの当該構成要素に脅威が作用した時の3.2-3

失敗シナリオとなる。
安全のための設備などが、起因事象の発生防止ならびにその影響緩和のために効果を発揮するべく設置されている。これら安全設備等を考慮し、失敗のシナリオとその影響度をそれぞれステップ 4 とステップ 5 で定量評価する。失敗シナリオの尤度はステップ 3、4、5 でそれぞれ評価されるが、段階が進むにつれて詳細度を高めていく。最後のステップ 6 では、リスク活用によるリスク管理にかかる意思決定の段階である。この時点においていずれの失敗シナ
リオに対してどのような体躯を施すことが適切かを優先度とともに決定することができる。
これらのリスク要因については画一的ではなく、その特性に応じて分類して評価する必要がある。
図 3.2-6 に一例を示す。使用済燃料プールからの燃料取出しをサブエンドポイントとするとき、燃料取扱機により燃料をプールから取り出すタスクを考える。ここで考えられる多様なリスク要因はタスクに対する脅威として定義され、その構造が図中に示されている。このタ
スクについて、考えられる脅威を、システム・機器に関する要因、社会的な要因、マネジメ
ント要因の３種類に分類した。システム・機器に関する要因は、ランダム故障、自然事象、
ヒューマンファクターにさらに分類する。社会的要因は、社会からの信頼失墜と悪意の行動に分類する。
図 3.2-6 に示されるように、タスクと脅威の組み合わせのそれぞれについて、一つのシナリ
オが定義される。

Step 1
Success Path
Description

DefineEndpoint

Step 2Threats
Identification

Identify
Threats

Step 3
Initiating Eventand Likelihood

Identify
Initiating
Event

Step 4
FailureScenario and
Likelihood

Describe
FailureScenario

Step 5Consequencesand Likelihood

Evaluate
Consequences

Describe
Provisions
Describe
Provisions

Identify Tasks

Order Tasks

Estimate
Possibility

Interpret Risk
and Risk
ProfileValue Impact
Analysis
Quantify
Likelihood

Estimate
Likelihood
DecidePossible
AffectStep 6
Risk-Informed
Decision
Making

Evaluate Risk
ImportanceComparealternatives

Judge Priority

Decisionmaking

図 3.2-5 リスク評価とリスク管理のための６つのステップ

3.2-4

Threat

System and Equipment Factor
RandomfailureElement of success path

NaturalhazardSocietal Factor
Human error

Public trust

Management
FactorMalicious
activity

Projectmanagement

Loss ofpower supply

Sub-Endpoint

Task

Fuel Retrieval
from Spent Fuel
Pool

Transport Fuelusing FuelHandling
MachineEarthquake(Small)Earthquake(Large)Typhoon /Strong windMiss
operation
Report minor
incident
Sabotage
Lack ofworkers

図 3.2-6 脅威の体系的分類
さらに、ここで定義されるそれぞれのシナリオに関してリスク評価指標が定義される。図3.2-7、図 3.2-8、図 3.2-9 は、それぞれ、起因事象（ステップ３）、失敗シナリオ（ステップ４）、シナリオの優先度付け（ステップ５）についてリスク評価指標を定義し、定量化の考え方を示すものである。

Step 3: Initiating Event (IE Likelihood)
Initiating Event
Description

Countermeasure and Reasoning of
Likelihood

Likelihood(H/M/L/VL)

Combinationof the task
and threat,

Built-in countermeasures are
considered.

H (high),
M (medium),
L (low) andVL (very low).

likelihood is evaluated with reason
possible
initiating
event

The reason will be used in the peer Initiating events
review that follows
VL are screenedout

図 3.2-7 起因事象に関するリスク評価指標

3.2-5Step 4: Failure Scenario (Scenario Likelihood and Risk Importance)
FailureScenario

Scenario Quantification and Likelihood
Reasoning of Likelihood(H/M/L)

Risk Importance

System failure
scenario
analysis

Considering mitigation
measures event
development (scenario) islisted up

H(high/unknown)
L (low).

Systemanalysis
methods(ET/FT, Graph)

ConsiderLikelihood of the scenario is mitigation
evaluated with reason
measures

H (high)
M (medium)L (low)

Failure scenarios
od rank L arescreened out

The reason will be used in
the peer review that
follows

図 3.2-8 失敗シナリオに関するリスク評価指標

Step 5: Scenario Prioritization
Consequence
Likelihood (H/M/L)
Priority
Consequence ofLikelihood of failure With the
failure scenario is
scenario andconsequence
evaluated in detailconsequence isand likelihood,evaluated in detailpriority(1-10)
Possible measures
is determined

図 3.2-9 シナリオの優先度付けに関するリスク評価指標
Kaplan と Garrick[文献(8)]はリスクトリプレットの概念を提案した。これによれば、リス
クは、1) what can go wrong?、 2) How likely is it?、 3) What are the consequences if
it happens?といった 3 つの質問い答えることにより特性づけられる。ここに示した 6 つのス
テップはリスクトリプレットの 3 つの質問に対応している。ステップ１とステップ 2 は問題のフレーミングでありすなわちシステムと脅威の定義である。ステップ 3 では起因事象一式が定義され、ステップ 4 では “what can go wrong”シナリオを与える。ステップ 5 は、損害の大
きさと影響度を評価するもので、 “what are the consequences?”に答える。“How likely
is it?”はステップ 3 から 5 で、それぞれ異なる詳細度で答えている。主要な評価項目は、成功パス、脅威、起因事象、シナリオ、影響度ということになる。シナリオの影響度と損害の様相に応じてカテゴリー化される。これによりリスクプロファイルとリスク曲線が得られる。
上記の考え方に従えば、デブリ取り出しについても実施すべきタスクを明確に定義できればそのリスク評価ならびにリスク管理をすることが可能である。一方、現時点では具体的な取り出し方法やデブリの位置と量が同定されていない。現時点で必要と考える作業を、液相部の閉じ込めについて図 3.2-10 に示す。また、気相部の閉じ込めについて図 3.2-11 に示す。
液相部の閉じ込めは、デブリの切断や切断時のデブリと液相の閉じ込めについて示している。気相部の閉じ込めは、切断後に原子炉容器内および格納容器内に浮遊すると考える切粉を含む気相を閉じ込めることと気相のサンプリングによって臨界を防止することを示してい3.2-6

る。今後、作業の具体化が進むにつれて、本節に示したリスク評価ならびにリスク管理手法
適用の具体化を検討する。

図 3.2-10 デブリ取り出し時のリスク要因（液相部の閉じ込めについて）3.2-7

図 3.2-11 デブリ取り出し時のリスク要因（気相部の閉じ込めについて）
3.2.3 リスクモニタリング、リスク管理（平成 30 年度）(1) デブリ取り出しにおける格納容器の状況分析
燃料デブリ取り出しにはさまざまな作業に関する意思決定をその都度行う。リスクモニタリ
ングとリスク管理のための考え方を提案する。燃料デブリを取り出す作業は、一時的にはリスクが増加する段階を経て、安定的にリスクが抑制された段階へ移行するなど、冷静にリスク分析を行なわれなければ、適正な手順を踏めない状況が生じうる。そこで、システムの状態異常を早期に確実に検出することが重要である。そのような作業のリスク変動の定量化とリスク低減効果を関係者が共有するために有効な、動的なリスクをモニタリングするシステムを提案する。
デブリ取り出し状況における格納容器内部と原子炉圧力容器内部の熱流動解析を実施した。
解析モデルを図 3.2-12 に示す。格納容器上部に面積 0.1m2 の流入口を、格納容器下部に面積0.1m2 の出口を両側部に設けた。原子炉圧力容器の底部は解放端とした。入口は 1m/s の流速一定
境界条件とした。格納容器底部に燃料デブリを模擬した熱源を置き、熱流束は崩壊熱相当（400J/m2•s）とした。図 3.2-12 左図は概略形状を、右図は実形状にあわせて詳細化したメッシュを示す。右図には、４箇所に想定した温度計の位置を示す。

3.2-8

図 3.2-12 格納容器内部と原子炉圧力容器内部の解析モデル
定常状態を得るため、40,000 秒の計算を行った。その時点での格納容器内部の流動状況を図
3.2-13 の(1)上部、(2)下部右側、(3)下部左側に示す。入口からの流れは遅いため、片側（ここでは右側）に偏流し、格納容器壁面に沿って流下する。そのまま格納容器底部の燃料デブリに沿って流れ、それを冷却するとともに左側の壁に沿って上昇流が生じている。
図 3.2-13 に格納容器内の温度分布（40,000 秒時点）を示す。格納容器内部はほぼ均一温度に安定冷却されている。燃料デブリ左端近傍に最高温度点があり、80℃程度に達している。なお、
図 3.2-14 に示すように、出口境界近傍（温度計１）の温度は 40000 秒後には定常状態が得られ
ている。

3.2-9

(1)上部(2)下部右側(3)下部左側
図 3.2-13 格納容器内部の流動状況

3.2-10図 3.2-14 格納容器内の温度分布

図 3.2-15 出口境界近傍の温度過渡変化(2) リスクモニタリング手法
その異常を検知することによりリスクをモニタリングする手法を提案する。40,000 秒において、
過熱事象と流量喪失事象が起きることを想定する。
① 過熱事象図 3.2-16 と図 3.2-17 に燃料デブリの熱流束が１０倍（崩壊熱が１０倍）になった条件での温度計１と２、３と４（図 3.2-12 参照）の温度変化をそれぞれ示す。温度計１と２は、主流の方向であり、温度変動が大きいことがわかる。一方、温度上昇も大きい。温度計３と４は低温側であり温度変動も温度上昇も小さい。Bayes 推定により、各温度計の計算値から、温度上昇を検出する確率を図 3.2-18 に示す。温度計３と４は温度変動が小さいため、いち早く異常を検知できている。温度計３あるいは４では、およそ 1100 秒後に異常と判断する確率が１に達している。
一方、温度計２では 1200 秒、温度計１では 1700 秒で異常と確信する結果となっている。温度計
3.2-11４は異常の診断は早期に可能であるが、600 秒程度から確率が上昇し、誤診断する可能性もあると考えられる。本来は、誤診断することなくかつ、早期に確実に異常を検知するアルゴリズムが求められる。

図 3.2-16 温度計１と温度計２の温度過渡変化

図 3.2-17 温度計３と温度計４の温度過渡変化

3.2-12図 3.2-18 各温度計の温度過渡変化より求めた異常と診断される確率
次に、複数の温度データを組み合わせて異常を診断する場合を考える。図 3.2-19 には温度計１と２の両方のデータを使った場合を、温度計１あるいは２を単独で使用した場合と比較して示す。複数の温度計データを用いることによって、誤診断することなく、早期に異常を検知できている。図 3.2-20 に示すのは温度計３と４を組み合わせて用いた場合である。この場合は、誤診断する確率は大きく低減され、かつ早期の診断が可能である結果が得られている。

図 3.2-19 温度計１と２を組み合わせて異常と診断される確率

3.2-13図 3.2-20 温度計３と４を組み合わせて異常と診断される確率
温度計を３つあるいは４つを組み合わせて異常を診断する確率を図 3.2-21 に示す。
このように３つの温度計を用いれば異常を診断する能力は大幅に向上することが示される。

図 3.2-21 温度計を３つあるいは４つを組み合わせて異常を診断する確率② 流量喪失事象
次に流量喪失事象の分析を行う。41000 秒まで熱流動計算を行った。それを初期条件として、41000 秒で入口流速が 1.0m/s から 0.1m/s に変化した場合に数値解析で得られた温度、流速を用いる。以下の過渡変化図では、40000 秒を原点（0 秒）として示す。図 3.222 に温度計１、２、３と４の温度変化を示す。まずは４カ所の温度を用いて温度異常推定を行う。図 3.2-23 は温度計ごとの数値解析で得られた温度と、異常なしシナリオで推定される温度、異常ありシナリオで推定される温度の比較である。
3.2-14図 3.2-22 数値解析で得られた温度と各シナリオで推定される温度の比較

3.2-15

図 3.2-23 各温度計のシナリオごとの信頼度
図 3.2-23 は温度計 1～4 の異常ありのシナリオが正しい確率を示した図である。図 3.223 の灰色の線は異常なしのシナリオの正しい確率、赤色の線は異常ありのシナリオの正しい確率を示している。温度計２，３，４では、異常発生(1000 秒)から 200 秒以内に異常を検知できている。温度計１では、異常の検知に 300 秒以上かかっている。これは、図3.2-23 の数値解析で得られた温度を見て分かるように 1000 秒から 1200 秒あたりまで温度が減少しているためと考えられる。
図 3.2-24 は 1 カ所の温度を用いた温度異常検知と、4 カ所の温度を用いた温度異常検知の結果で、異常ありシナリオの信頼度を示している。

3.2-16

図 3.2-24 １カ所の温度を用いた場合と 4 カ所の温度を用いた場合の異常ありシナリオの信頼度
表 3.2-1 シナリオの説明

シナリオ１

異常なし

シナリオ２

熱流束800J/m2

シナリオ３

熱流束2800J/m2

シナリオ４

流入速度0.5m/s

シナリオ５

流入速度0m/s

本項では、前項で行った温度異常検知の結果から、t 異常=1170[s]としてシナリオ推定を行う。表 3.2-1 のように 5 つのシナリオを考えることする。シナリオ１，２は熱源から発生する熱流束の上昇という異常である。シナリオ 1 は熱源から発生する熱流束が 800Jになった場合で、燃料デブリの一部が再臨界を起こした場合を想定している。シナリオ 2は熱源から発生する熱流束が 2800J になった場合で、シナリオ 1 よりも多くの部分で再臨界を起こした場合を想定している。シナリオ３，４は流入気体の減少という異常である。
シナリオ３は冷却材である空気の流入速度が 0.5m/s になった場合で、冷却機能が半分喪失した場合を想定している。シナリオ 4 は冷却材である空気の流入速度が 0m/s となった
場合で、冷却機能が完全に喪失した場合を想定している。
図 3.2-25 が 4 カ所の温度を用いてシナリオ推定した結果である。温度異常検知からある程度の時間が経つまでシナリオ２の信頼度が大きく、途中からシナリオ５の信頼度が大きい。この傾向はケース１、２と同じである。

3.2-17

図 3.2-25 ４カ所の温度を用いてシナリオ推定した結果

図 3.2-26 ４カ所の温度と流速を用いたシナリオ推定の結果
図 3.2-26 が 4 カ所の温度と出口流速を用いてシナリオ推定した結果である。この結果
は、４カ所の温度を用いた場合と似た結果となっている。異なる点としては、シナリオ
３の信頼度が最大でも 0.8 を超えない点と、シナリオ５の信頼度が大きくなる時間が少し早いことである。このことから流速を用いることで、シナリオの推定の精度が改善されているといえる。以上に述べた通り、本研究にて開発した手法は、複数のランダムな温度データから異常を確実に診断できるアルゴリズムであり、これを用いれば、燃料デブリ取り出し時のさまざまな異常事象に対してリスクモニタリングが可能であることが確認された。
3.2-18

英国共同研究先の研究紹介
3.3.1 序論

デブリ・ベッドが存在する反応コアに水を注入する冷却過程は、複雑な多相流れ、熱移動、沸騰などが関わるため、過程の理解およびモデル化は困難である。本研究では、
この複雑なマルチ・フィジクスな現象について効率的な数値計算フレームワークを提案する。
多相流れの系は最近数十年の間に数値流体力学(CFD)の文脈で多数の研究例が見られる。この流れには、数多くのメカニズムが関与しており、例えば、いくつかの相どうしの相互作用、相変化、熱移動、密度コントラストによる不安定性などがある。
気体、液体、固体タイプの流れについては、全ての相は連続体として考えることができる。相境界では、質量、運動量、エネルギー交換が連続的に生じ、それぞれの相は異なる温度と速度を持っている。大規模スケールでこの流れを解くには、混合モデルを使うことが考えられる。ここで、液体の水、蒸気、固体相を異なる準連続的相と捉え、それぞれ固定オイラーモデルにおいて物理的な相分布から独立していると考える。
このとき、ミクロレベルでの相境界の場所を解くよりも、流体流れと熱移動のマクロな特徴を予測することができる。混合モデルは、混合システムの運動量、連続性、エネルギー方程式を解くことで、n 個の相（流体または粉体）をモデル化することができる。マクロなモデルと対照的に、相境界のトラッキングが、レベルセット法、VOF 法などで可能となる。
相間での運動量交換については、いくつかの考慮すべき力がある。主なものはドラッグで、固体相あるいは気体相から流体相にはたらく。それ以外の揚力などの力も考慮する必要がある。蒸着、沸騰過程のケースでは、熱あるいは質量交換についても考える必要がある。
3.3.2 連続体モデリング
各相（ni）の流れはそれぞれ異なる特性、密度、速度、圧力を持ち、いずれもオイラー記法で記述される。各相についてナビエ・ストークス式を考える。
3.3.2.1 質量保存ΓB = Γ"Q, #を相Bの単位体積あたりの生産率とする（場所Q，時間 ）。トータルでの質量は一定とし、相の間でのみ移動を考える。
T

RΓ = 0(26)

相Bの体積をΩ とし，動きを追跡する．対応する相Bの質量保存式はで記述できる。体積分について微分をとると、
W
U V "
W
YZ

#+∇ ∙"

3.3-1\ #−Γ] = 0(27)

Ω は任意の開集合なので，
W\+ "\ ∙ ∇#\ = −∇ + ∇ ∙ ^_ + `W(28)となる。ΓB は質量交換に相当するソース（またはシンク）である。
3.3.2.2 線形運動量保存式
W
"
W

#+∇∙"

\ # = Γ , B = 1,2,3(29)

ただし、! ′ はコーシー応力の非対角成分、 は圧力である。複数流れ問題では、
W
"
YZ W

U V

#+∇ ∙"\ #−Γ] = 0(30)

ただし、aB は相間連成項で、ドラッグ、付加質量、相間圧力などからなる（詳細は 4.3 節）。

3.3.2.3 エネルギー保存bc+ ∇ ∙ d − ∇ ∙ "^\# = M
b

(31)

We
+ fJ "g ∙ ∇#e − ∇ ∙ ∇e = M + ∇ ∙ "^\#
W

(32)

ただし、cは内部エネルギー、dは熱流ベクトル、Mは放射による熱吸収。dはd =− ∇eで表される（ は熱伝導度）。すると、エネルギー方程式はfJとなる。コーシー応力の中で粘性項を無視して、連続の式を入れると、∇ ∙ "^\# = − ∇ ∙ \ =熱方程式は、
fJとなる。

が得られる。

(33)

We
1b+ fJ "\ ∙ ∇#e − ∇ ∙ ∇e −=MWb

多数相のケースでは、各相Bについて
fJh

1b
b

We
+
W

fJ "\ ∙ ∇#e − ∇ ∙b∇e − 4
6 + ai =
b

(34)

M

(35)

3.3.2.4 状態方程式(EOS)
状態方程式(EOS)は、測定値を使って、純物質の圧力/体積/温度（PVT）挙動を記述する。圧縮性流体を仮定したとき、
3.3-2

fJh

We
+
W

fJ "\ ∙ ∇#e − ∇ ∙b∇e − 4
6 + ai =b

ただし、>は音速、j/kJ は熱膨張係数である。

M

(36)

3.3.3 相変化、連成そして相関3.3.3.1 熱移動由来の質量移動に関する境界 Stefan 条件

気化および凝結について、運動エネルギーおよび粘性作用項を無視すると、境界での質量転移が記述できる。体積平均した質量転移率Γは
Γ=

境界への純熱移動(37)

気化の際の潜熱として記述できる。質量保存でのソース（シンク）・タームはΓ liquid phasel = m −Γ gas phase
0 solid phase(38)で与えられる。気化および凝結が生じる条件は、熱流束のバランスで決まる。境
界間質量転移は
Γ=−

vtJ
st
∑vtJ
− eu # − ∑w ,ew − e st u "e

xy

(39)

でパラメータ化可能である。ただし、xℎ = ℎz − ℎ{ は相変化のエンタルピーに相当|}する。熱伝導係数Σ{
|}

Σ{
=

⎧max m ‚{ 12 „e…
⎪
@} ,{
ƒ
⎨
⎪
⎩

min ŽR

ただし、j = 1、|}およびΣz は次式で計算できる。
{

x
z y

max
u

f {

j,

‚{

ƒ

42 + 0.74 ‡ { c{ ˆŠ xy w u w
,
•
u− w

Σw

vtJ

=

‰.Š

3.6 × 10’ “
ƒ“

6‹

3.6

if e

ƒ

@}

if e

> e{

@}

G e{

(40)

N5

= max" z , 10 #である。

最終的な液体および気体の熱伝導係数の最終値は|} ,maxΣ{=⎧ min ”• z , { { – xy
⎪ Δ
max ™|e st − e |
u

Mc˜

⎨
• z z xy⎪st
⎩Δ Mc˜ max ™|e − eu |
3.3-3

if e

@}

if e

@}

G e{

> e{

(41)

⎧ min›• w , u u œ xy
⎪ Δ =/0 max ™|e st − ew |vtJ,max
Σw=• w w xy⎨⎪st
⎩Δ =/0 max ™|e − ew |if e st G ew

if e st > ew

で制約を受ける。ただし、™ = 10N5‰ 、• = 10、Δ =/0 = 0.01 sである。
3.3.3.2 相間運動量移動(1) ドラッグ・フォース

相間で作用するドラッグ・フォースは相間の相互作用として最も重要である。ドラッグ・フォースa1 は，\z , \u , \s を未知数として

Σuw ,\w − \u - + Σsu "\s − \u # 液相

a5 = •Σuw ,\u − \w - + Σsw ,\s − \w - 気相(42)

Σsu "\u − \s # + Σsw ,\w − \s - 固相

液相および固相を連続相とし、気相は粉体相として仮定すると、相間ドラッ
グ・フォースは Ergun の式から
Σ@z = 150
Σ@{ = 150

Σuw =

3
f
4 Ÿ

2

ž

z@ z
z@
2

2

ƒ

ž

{@ z
@{

2

ƒ

+ 1.75z@ z+ 1.75z@ {uw wu u |\w − \u |ƒJ|\z − \@ |ƒ|\{ − \@ |

,

ƒ

(43)

N . Š
uw

で記述できる。流体中のバブルについては、流体-気体間の抗力係数を使い
24‰. ¢£¡1 + 0.15, uw cuw ¤ if uw cuw G 1000 粘性領域
fŸ = • uw cuw
0.44
if uw cuw > 1000 乱流領域

ただし、 c{z は

cuw =平均バブル直径ƒ はƒ“ =u ¥\w − \u ¥ƒ“

žu

¦/ !

u ,\w − \u -で計算される。正規化された体積分率は

3.3-4

(44)

(45)

(46)

で求められる。ただし、

§ =

+

§

(47)

§ は相¨における相Bの体積分率である。

(2) 仮想質量力

粒子またはバブル（質量 m）が流体中を移動するとき、流体から粒子に効力が作用する。これは付加質量力あるいは仮想質量力として知られている。付加され
る力は

ƒg
ƒg
⎧− w u © u − w ª
液相
2
ƒ
ƒ
⎪
ƒgw
a =
w u ƒgu気相⎨ 2 ©ƒ − ƒ ª
⎪
0
固相
⎩

で表される。

(48)(3) 境界圧

特に気化または凝結が生じる流れにおいて、相間の質量転移は、次のようにモデル化される。
a. =

0

⎧
⎪1

液相

w u ,\w − \u - ∇ w

⎨4
⎪
⎩

0

(49)

気相

固相

3.3.3.3 相間熱転移のパラメータ化

温度方程式（eu , ew , es が各相に対応する未知の温度変数)に付加される温度交換は⎧Σuw ,ew − eu - + Σus "es − eu #iii
˜5 = Σuw,eu − ew - + Σsw,es − ew ⎨ ii⎩Σsu "eu − es # + Σsw ,ew − es i

i

液相

気相

(50)

固相

で記述できる。沸騰による温度交換は
Σu

˜ =•
i

vtJ

"e st − eu # + Γu "ℎu − fJu eu + x/ #vtJ
Σw ,e st0

− ew - + Γw "ℎw − fJw ew #

液相
気相

(51)

固相

伝熱係数は、無次元の Nusselt 数で表現できる。
«g =

Σƒ

= 2 + 0.6 c 5/ ¬M5/.

(52)« 、 c、¬Mは無次元の数で、すべり速度およびバブル直径ƒ“ 、粒子直径ƒJ を用いて次のように書くことができる。

3.3-5

csu =

u |\u − \s |ƒJ

žu

csw =

w ¥\w − \s ¥ƒJ

žw

cuw =

u ¥\w − \u ¥ƒ“

žu

(53)

3.3.4 数値計算スキーム
複相モデルは次のような構成として要約できる．
未知数: \ ,,方程式:, e および 全領域質量保存式
 各相における体積分率 各相におけるナビエ・ストークス運動量方程式
 各相での温度方程式
 状態方程式（EOS) 相間相互作用および相変化に関する関係式
3.3.5 数値計算結果

空隙 率 s = 0.4で、粒子 径は0.001 mを想定 した。計算 空間の寸法 は、 0.09 × 0.09 ×

0.25 mとした。充填層の密度は8,000 kg/m3 、熱容量は0.5 kJ/kg K、ゼロ熱伝導度であった。
重力z = 9.81 m/s2 とした。液相中に発生させた摂動を沸騰の誘導に用いた。摂動領域は球形で、中心座標は(0.045, 0.045 0.1)に置き、0.02 mの半径とした。サイドの 4 面にはすべり境界をもうけ、上側壁面は開放して大気圧条件とした。計算領域のメッシュは各軸について20 × 40 × 40ノードの三角メッシュをもうけた。シミュレーションの様子は図3.3-1 に示す通りである。

図 3.3-1 蒸気の体積分率

時刻 = 0.04, 0.05 s 摂動温度104 ℃

3.3.6 結論
本研究は、オイラー的複相モデルを用いてデブリ層の 3 次元複相流れをシミュレーショ
ンした。実装をするにあたり、既存のオープン・ソース・ソフトウェアである FLUIDITYを用いた。各相での連続体モデルを連携させながら、それぞれの相での保存則に関わる輸送方程式を各相で解いた。相間連成モデルは、大規模スケール流れのマクロな挙動をシミュレーションするのに最適である。

3.3-6

研究推進
東京大学側で行った固気液三相流の数値シミュレーションおよびを検証実験について、ICLと廃炉事業への適用性について議論した。DEM-DNS-VOF 法が気中工法における燃料デブリから発生する切粉の飛散防止検討に役立てられることを確認した。さらに、混相流の数値シミュレーションと臨界計算との連成手法についても検討事項を確認した。
共同研究打合せおよび意見交換会（セミナーを含む）も行った（図 3.4-1 および図 3.4-2）。
平成 29 年度は、平成 29（2017） 年 7 月 18 日に Imperial College London に訪問し、共同研
究打合せおよび意見交換会（セミナーを含む）を行った。平成 30 年度は、平成 30（2018） 年
10 月 19 日および平成 31（2019） 年 3 月 8 日・平成 31（2019） 年 3 月 11 日に Imperial
College London を訪問し、共同研究打合せおよび意見交換会（セミナーを含む）を行った。
以上のように、本事業は、研究代表者の下で、ICL および東京大学内の各研究項目担当者間における連携を密にして研究を進められた。

図 3.4-1 意見交換会後の集合写真

図 3.4-2 山口彰研究代表者による講演

3.4-1

4.

結言
この研究では、以下の内容について行った。
・福島第一原子力発電所の実際のジオメトリーを参考にして、気中工法における燃料デブリの冷却性能評価用の熱流動解析を行った。
・上記の結果について、Imperial College London と議論した。
・これまでの研究成果に基づき継続的に高度化することを含め燃料デブリ取り出し作業のリ
スク変動を定量化し、リスク低減効果を関係者が共有するために有効なリスクをモニタリングする手法を検討しまとめた。実作業中に実測可能なパラメータが反映できることを考
慮した。

4-1

添付資料