炉物理の研究 第74号 (2022年3月)

第 52 回炉物理夏期セミナー開催報告(3)受講者の立場から

若手研究者報告

生体内の光輸送計算における反射効果に対する中性子輸送計算法の適用性
名古屋大学 工学部 エネルギー理工学科 山本研究室
学部 4 年 伊藤 雅人

1. 序論
医用工学分野では、X 線 CT や核磁気共鳴法といった生体イメージング技術が診察において、重要な役割を果たしている。各種法の課題として、X 線 CT は被曝を避けられない点、核磁気共鳴法は検査時間が長い点があげられる[1]。本研究では、光を用いた生体内イメージングの別手法である拡散光ト
モグラフィ(diffuse optical tomography, DOT)に注目する。DOT では近赤外線を用いて、生体内イメージングを行う。このため、検査時の被曝による健康被害の回避が期待できる。以下で、DOT の原理について述べる。まず、近赤外線をある表面から入射させ、別の表面から計測する。次に、生体組織分布を仮定し、入射光の生体組織内での振る舞いを解析する。具体的には、光学特性値は生体組織によって異なることを利用する。光学特性値分布𝝁𝒂 、𝝁𝒔 を仮定し、(1.1)式の輻射輸送方程式(radiative transfer
equation, RTE)を解く。最後に、RTE より得られる解と測定結果を比較する。以上の操作を両者が一致するまで繰り返す[1]。これらの結果が一致した場合に仮定した生体組織分布は正しいといえる。
[𝛀 ∙ 𝜵 + 𝜇 (𝒓) + 𝜇 (𝒓)]𝐼(𝒓, 𝛀) = 𝜇 (𝒓)𝑑Ω 𝑝(𝛀, 𝛀 ) 𝐼(𝒓, 𝛀 ) + 𝑞(𝒓, 𝜴)(1.1)

以下に、(1.1)式の変数を示す。
𝜇 (𝒓)：位置𝒓の吸収係数 [1/m]𝜇 (𝒓)：位置𝒓の散乱係数 [1/m]𝐼(𝒓, 𝛀)：位置𝒓、角度𝛀の光強度 [1/m2/s]𝑝(𝛀, 𝛀 )：位相関数
𝑞(𝒓, 𝜴)：位置𝒓、角度𝛀の外部光源 [1/m3/s]
このように、DOT では、RTE の解を求める必要がある。しかし、RTE は多変量微積分方程式であり、一般解が得られていない。そのため、数値計算を行う必要がある。本研究では、中性子輸送計算手法の適用性について検討する。生体は、異なる生体組織が配置されているような体系である。異なる生体組織の境界面では、光は屈折や反射をする。これらの現象を考慮した場合、中性子輸送計算手法で用いる角度分点が生体内光輸送解析において適切な角度分点ではないことが考えられる。より適切な
角度分点を用いて計算を行うことで、少ない分点数で精度の良い計算を行うことが可能となるため、生体内光輸送解析において適切な角度分点を選択することは、検査時間の短縮の面で重要である。そこで簡単な体系において反射効果を考慮した場合の生体内光輸送解析に適した角度分点セットについて検討をした。

1

炉物理の研究 第74号 (2022年3月)

2. 検証計算本研究では、均質な体系を用いて角度分点に関する検討を行った。まず、外部境界条件に真空境界条件を課し、複数の角度分点を用いて角度分点数を変化させながら(2~100 点)計算を実施し、各条件において誤差評価を行った。次に、外部境界条件に角度依存性のある反射境界条件を課し、同様に誤差評価を行った。

1)

計算体系・計算条件

本検討で用いた計算体系は、Fig. 2.1 のような一端に光源が配置された均質な長さ 2 cm の一次元無限平板体系である。

0.1 cm

1.9 cm光
源
Fig. 2.1 一次元無限平板体系

Table. 2.1 に光学特性値や生体内屈折率、メッシュ数、角度分点数について述べる。
Table. 2.1 計算条件
0.15
90
1
1.4
1
2
20
2~100

𝜇
[1/cm]
𝜇 [1/cm]
𝑞 [1/cm3/s]𝑛 (生体内)𝑛 (真空)
体系長さ [cm]
メッシュ数𝑁
角度分点数 M2)

誤差評価

Fig. 2.1 の体系において外部境界条件として真空境界条件を与え、Table. 2.2 に示した 4 つの分点を用いて計算を実施した。
Table. 2.2 真空境界条件を課した際に用いた角度分点
等角度分点
等方向余弦分点
Gauss Legendre (GL)分点
Double Gauss Legendre (DGL)分点

2

炉物理の研究 第74号 (2022年3月)

DGL 分点を用いて角度分点数を 256 点取り、計算を実施した数値解を参照解とした。角度分点数を2~100 まで取り、これらの分点を用いて計算を実施した。得られた光子束および参照解により(2.1)式のように計算される相対値の RMSE を計算し、誤差評価を行った。各角度分点で計算した場合の角度分
点数に対する RMSE の分布を Fig. 2.2 に示す。
相対値の RMSE =

𝜙 −𝜙

1
𝑁

(2.1)

𝜙

以下に(2.1)式の変数を示す。
𝜙 ：対象の条件の i 番目のメッシュの光子束の数値解
𝜙：i 番目のメッシュの光子束の参照解

𝑁：メッシュ数

1 E+04等角度

相対値のRMSE[%]

1 E+02
1 E+00

GL分点

1 E-02等cos

1 E-04
1 E-06
1 E-08

DGL分点

1 E-101 E-12
0

20

40

60

80

100

角度分点数 [-]
Fig. 2.2 角度分点数に対する RMSE の分布 (真空境界条件)

Fig. 2.2 から、Fig. 2.1 の一次元無限平板体系の外部境界条件に真空境界条件を与えた場合の計算において、検討した 4 つの角度分点のうち最も適切な角度分点は DGL 分点である。ここで、光源が配置されたメッシュ内での光子束の角度分布を Fig. 2.3 に示す。Fig. 2.3 を見ると、このメッシュ内での光子束の角度分布は𝑐𝑜𝑠𝜃 = 0で滑らかでない。DGL 分点は−1 ≤ 𝑐𝑜𝑠𝜃 ≤ 1の角度領域を𝑐𝑜𝑠𝜃 = 0で分割し、−1 ≤ 𝑐𝑜𝑠𝜃 ≤ 0及び0 ≤ 𝑐𝑜𝑠𝜃 ≤ 1の領域それぞれに GL 分点を与えることで得られる分点である。以上を踏まえると、任意の体系・計算条件に最適な角度分点は以下のステップで得られることが期待できる。

3

炉物理の研究 第74号 (2022年3月)

1) すべてのメッシュ内の光子束の角度分布を調べ、それらの角度分布が持つ不連続点あるいは滑らかでない点の角度を特定する。
2) 調査した角度点で、角度領域(−1 ≤ 𝑐𝑜𝑠𝜃 ≤ 1)を分割し、各領域に GL 分点を与える。

領域平均角度光子束 [cm-2s-1]

2.60E-02
領域12.55E-02
2.50E-02
2.45E-02
2.40E-02
2.35E-02

滑らかでない2.30E-02
2.25E-02
2.20E-02
-1.0

-0.5

0.0

0.5

1.0

cosθ
Fig. 2.3 光源が配置されたメッシュ内の光子束の角度分布(真空境界条件)

続いて、Fig. 2.1 の体系において外部境界条件として角度依存性のある境界条件を与えた場合について同様の検討をした。
角度依存性のある境界条件は(2.2)式で与えられる。
𝜓= R[θ ]𝜓(2.2)

以下に、(2.2)式の変数を示す。
𝜓𝜓

：m 番目の角度分点の境界面から放出される光子束強度：M-m+1 番目の角度分点の境界面に入射する光子束強度

M：角度分点数
R[θ ]：反射率
θ ：入射角反射率R[θ ]は生体及び真空の屈折率が決まれば、光の入射角θ のみに依存する。今回の場合、生体内から生体外へ放出される光の反射を考えればよい。本検討では生体内の屈折率𝑛 = 1.4、真空の屈折率
𝑛 = 1であるので、生体内から生体外へ放出される光の反射率R[θ ]の角度依存性を Fig. 2.4 に示す。
（生体内の屈折率𝑛 ）>（真空の屈折率𝑛 ）より、Fig. 2.4 のようにある角度𝜃 以上で光の反射率Rが 1 となる臨界角が存在する。
4

炉物理の研究 第74号 (2022年3月)

1.2

反射率 R[-]

1

R

0.8
0.6
0.4
0.2

臨界角

0
0

0.25

0.5

0.75

1

cosθ1 [-]
Fig. 2.4 反射率の角度依存性

この場合について角度分点数 88 点とし、DGL 分点を用いて計算を実施した。体系両端のメッシュ
内の光子束の角度分布を Fig. 2.5 に示す。

光源

2.9E-04
領域1

領域平均角度光子束 [cm-2s-1]

領域平均角度光子束 [cm-2s-1]

5.0E-02

4.5E-02

4.0E-02

3.5E-02

3.0E-02
-1.0

臨界角
-0.5

0.0

2.7E-04
2.6E-04
2.5E-04
2.4E-04
2.3E-04

0.5

1.0

-1.0

cosθ

領域20

2.8E-04

臨界角
-0.5

0.0

0.5

1.0

cosθ

Fig. 2.5 体系の両端のメッシュ内の光子束の角度分布(角度依存性のある反射境界条件)

Fig. 2.5 のように両端のメッシュ内の光子束角度分布は、𝑐𝑜𝑠𝜃 = 0及び臨界角𝑐𝑜𝑠𝜃 = ±𝑐𝑜𝑠𝜃 で滑らかでない。真空境界条件での検討から、角度領域−1 ≤ 𝑐𝑜𝑠𝜃 ≤ 1をこの三点で分割し、各領域

5炉物理の研究 第74号 (2022年3月)−1 ≤ 𝑐𝑜𝑠𝜃 ≤ −𝑐𝑜𝑠𝜃 , −𝑐𝑜𝑠𝜃 ≤ 𝑐𝑜𝑠𝜃 ≤ 0, 0 ≤ 𝑐𝑜𝑠𝜃 ≤ 𝑐𝑜𝑠𝜃 , 𝑐𝑜𝑠𝜃 ≤ 𝑐𝑜𝑠𝜃 ≤ 1のそれぞれに GL 分点を与えた分点が Fig. 2.1 の体系の外部境界条件に角度依存性のある反射境界条件を与えた場合に適切な分点であると考えられる。本検討ではこの分点を Quadruple Gauss Legendre (QGL)分点と呼ぶこととする。
Table. 2.2 に示した 3 つの分点を用いて計算を実施した。
Table. 2.3 角度依存性のある反射境界条件を課した際に用いた角度分点
Gauss Legendre (GL)分点
Double Gauss Legendre (DGL)分点
Quadruple Gauss Legendre (QGL)分点
QGL 分点を用いて角度分点数を 256 点取り、計算を実施した数値解を参照解とした。角度分点数を2~100 まで取り、これらの分点を用いて計算を実施した。得られた光子束および参照解により(2.1)式のように計算される相対値の RMSE を計算し、誤差評価を行った。各角度分点で計算した場合の角度分
点数に対する RMSE の分布を Fig. 2.6 に示す。

1 E+03

相対値のRMSE [%]1 E+02
1 E+01
1 E+00
1 E-01

DGL分点

1 E-02

QGL分点

1 E-03
0

20

40

60

80

100

角度分点数 [-]
Fig. 2.6 角度分点数に対する RMSE の分布 (角度依存性のある反射境界条件)

Fig. 2.6 より多くの角度分点数で QGL を用いた計算での相対値の RMSE が最小である。また、角度分
点数に対する相対値の RMSE の分布をみると GL 分点あるいは DGL 分点を用いた場合は振動してい6

炉物理の研究 第74号 (2022年3月)

るのに対し、QGL 分点を用いた場合は単調減少している。以上より一次元無限平板体系の外部境界条
件に角度依存性のある反射境界条件を与えた場合の計算において、検討した 3 つの角度分点のうち最も適切な角度分点は QGL 分点であることが分かった。

3. まとめと今後の課題
1)

まとめ

本研究では、角度依存性のある反射効果を考慮した生体内光輸送計算における最適な角度分点について検討を行った。ここでは均質な一次元無限平板体系を用いて検討を行った。始めに外部境界条件
として真空境界条件与えた場合について検討を行い、全メッシュ内の光子束の角度分布が不連続ある
いは滑らかでない点で角度領域を分割し、各領域に GL 分点を与えたものが対象とする体系、条件での輸送計算に適した分点であると予想した。次に、先の検討で得た予想を基に外部境界条件として角度依存性のある反射境界条件与えた場合について検討を行い、体系の臨界角に注目することで反射を考慮した場合の対象とする体系の輸送計算に適した分点を与えることが可能であることが分かった。

2)

今後の課題

今回は、単純な体系を用いて検討した。一般に生体は複数の生体組織からなる。このため、体系内の非均質を考慮した体系を扱わなければならない。この場合、本検討では考えなかった屈折を取り扱う必要がある。MOC 計算においてこの屈折を取り扱いが課題となる。また、このような複数の臨界角を持つような体系において本検討で提案した角度分点の取り方が有効であるかを調査する必要もある。

4. 参考文献[1] 天野虎之介, “光輸送解析への中性子輸送計算手法の適用性に関する検討”, 卒業論文, 名古屋大学(2021)

7

炉物理の研究 第74号 (2022年3月)

第 52 回炉物理夏季セミナー報告
若手研究会報告
遺伝的アルゴリズムを用いたトリウム燃料集合体の転換比の検討、暗中模索の日々
大阪大学大学院 工学研究科 環境エネルギー工学専攻
原子力社会工学領域（北田研究室） 修士２年 竹石太一

1. 序論
トリウムは、その豊富さと熱中性子領域での高い中性子再生率と転換率から高転換軽水炉の燃料として考えられている。しかし、トリウム燃料に関する商業的な製造・再処理イン
フラはまだ整っておらず、実用化はされていない。トリウム燃料のメリットをより高めるよ
うな条件が提案できれば、トリウム燃料を用いる選択肢が商業的に魅力的なものになっていき、実用化に向けた市場の変化が生まれる可能性もある[1]。
トリウム燃料の実用可能性を広げるための一つの観点として、より高い転換比を得られる条件をまとめることが挙げられる。本研究ではトリウム燃料を一般的な PWR 集合体の条件で燃焼させた際の、燃料棒配置と転換比の関係を明らかにする検討を行った。
PWR 集合体は燃料棒の組成とその配置によってその核特性が変わり、転換比も変化する。
具体的な燃料棒の組み合わせは非常に多様なケースが考えられ、すべてについて網羅的に検討するのは莫大な時間が掛かってしまう。そこで、遺伝的アルゴリズムを用いることで、これらの多様な選択肢の中から効率的かつ網羅的に転換比の高い配置を探索することができるのではないかと考えた。
若手研究会発表では、遺伝的アルゴリズムを用いて有意義な結果を得られるまでの奮闘の記録と、得られた結果について報告した。

2. 遺伝的アルゴリズムの条件設定
まず、今検討で考慮する PWR 集合体の条件設定を説明する。
燃料棒が 17×17 本で構成される PWR 集合体の 1/8 対称燃料棒配置を対象に考える。用いる燃料棒は、U233 と O16 からなる UO2 を 0〜4%、それ以外は TH232 と O16 からな
る ThO2、という 5 種の組成を用意して、これらを 1/8 集合体に自由に配置する条件を設定した。

炉物理の研究 第74号 (2022年3月)

Fig.1 集合体の条件設定
遺伝的アルゴリズムについては python の DEAP というライブラリを用いて実装した。
ピン配置データを遺伝子として用いており、遺伝子を基に SCALE[2]中の Polaris モジュール（2DMOC 計算コード）を用いて燃焼計算を行い、初期 U233 数密度に対する 50GWd/t
での U233 の数密度の比で算出した転換比を評価関数として採用した。
遺伝子を変化させる設定としては、交叉と突然変異を実装した。交叉は 2 点交叉を用いており、世代中の個体を順に 2 つ 1 組ずつ選出し 0.6 の確率で交叉を行う。突然変異は世代中の個体の 4 割を選び出し、それぞれの遺伝子に対して 0.05 の確率でランダムに組成を変更する。フローチャートを以下に示す。

Fig.2 遺伝的アルゴリズムのフローチャート
以上の条件で 1 世代の個体数 300・50 世代分の計算を行ったところ、すべてのピンが 0%になるピンが最適解として得られた。

炉物理の研究 第74号 (2022年3月)

Fig.3 得られた最適解
評価関数は初期 U233 数密度を分母に持つ関数なので、これが小さくなるほど転換比は向上する。そのため遺伝的アルゴリズムによって初期 U233 数密度の最小化が行われてしまった、という結果に終わってしまった。

3.数密度を固定するような条件の再設定
自由度が高すぎる設定では有意義な結果が得られなかった。本研究では集合体中の燃料配置に関する考察を行いたいので、燃焼初期の集合体全体の U233 の数密度が一定になる
ように条件の再設定を行った。
初期世代の個体生成・交叉・突然変異の段階において、集合体全体の U233 数密度が変化しないように再実装を行った。個体生成についてはランダムに生成した集合体中の U233 数
密度の総和がある値になった時だけ採用し、交叉・突然変異についてもその操作を行った後でそれぞれの個体の集合体中の U233 数密度の総和が変わっていない場合だけ採用し、そうでない場合はやり直すという形式をとっている。
以上の条件で 1 世代の個体数 300・50 世代分の計算を行ったところ、以下の図のようなピン配置が得られた。

Fig.4 条件再設定後に得られた最適解

炉物理の研究 第74号 (2022年3月)

0%のピンと 4%のピンが多い結果が得られた。また、集合体の角に U233 の多い燃料棒が集中し、集合体の中央に U233 のない燃料棒が入る結果が得られた。
この結果から、一様な集合体ではなく、U233 と Th232 を分離して集合体中に装荷した方が良い転換比を示していることが分かる。
原因の考察としては、断面積による影響が考えられる。熱中性子領域での断面積に着目す
ると、U233 の吸収断面積が Th232 の 100 倍近くある。そのため、U233 と Th232 が同じ燃料棒に均質的に配置されていると、U233 が熱中性子を捕獲してしまう分 Th232 の捕獲反応が抑えられ、転換比が低くなるのではないかと考えられる。
しかし、この考察だけでは中性子源が集中する理由の説明にはならない。今後、詳細なデータでの検討を進めていく。

4、結言
本稿では、トリウム燃料を用いた PWR の転換比を遺伝的アルゴリズムで最適化した結果について報告した。
遺伝的アルゴリズム中に、集合体中の U233 の数密度を固定するという観点を組み込むことで、集合体中の U233 の配置が非均質な集合体の方が良い転換比を持つという結果が得られた。
U233 と Th232 の断面積を考慮すると、U233 と Th232 を分離した燃料棒配置の場合にTh232 の捕獲反応が効率よく進むため、転換比が高くなると考えられる。
参考文献[1] INTERNATIONAL ATOMIC ENERGY AGENCY, Role of Thorium to Supplement Fuel
Cycles of Future Nuclear Energy Systems, Nuclear Energy Series, 2012[2] B. T. Rearden and M. A. Jessee, Eds., SCALE Code System, ORNL/TM-2005/39, Version
6.2.3, Oak Ridge National Laboratory, Oak Ridge, Tennessee (2018). Available from
Radiation Safety Information Computational Center as CCC-834.

炉物理の研究 第74号 (2022年3月)

反応度事故を想定した核熱結合計算モジュールの開発
北海道大学大学院

1.

工学院

原子炉工学研究室

原口太志

緒論

1.1 核分裂生成物の計算における取り扱いと EFP モデル
核分裂反応が起こると、反応と同時に励起状態の核分裂片（fission fragments）から中性子が発生する。一連の流れで発生する中性子の 99%以上の大部分はこのとき発生し、これを即発中性子（prompt neutron）という。即発中性子、即発ガンマ線を放出し基底状態になったものが核分裂生成物（fission products, FP）であり、この中には崩壊して中性子を放出するものが存在する。この核種のことを遅発中性子先行核（delayed neutron precursor）と言い、この遅発中性子
先行核の崩壊によって核分裂反応から少し時間をおいて発生する中性子を遅発中性子（delayed
neutron）という。例えば 235𝑈 熱中性子核分裂において即発中性子に対する遅発中性子の割合は
約 0.65%程度と割合はかなり小さいが、この遅発中性子によって原子炉の過渡変化が緩やかになることで原子炉の制御にとって極めて重要な役割を果たしている[1]。
遅発中性子先行核としては、87𝐵𝑟 、88𝐵𝑟 、137𝐼 、138𝐼 、139𝐼 などの核種が存在し、その数はおよそ 300 核種である。またその各々が異なる半減期でベータ崩壊し遅発中性子を放出する。先行
核がそれぞれのタイミングで中性子を放出すればそのタイミングで原子炉の出力が変化することになるので原子炉動特性解析においてはこれらすべてを別々に考慮すべきであるが、計算コストがとても大きくなってしまうため、遅発中性子先行核をその崩壊定数𝜆によって６つにグルーピングし、６種の仮想的な先行核として取り扱う６群モデルが一般的に用いられている[2]。
このような近似的な取り扱いは十分に高い精度を有していると考えており、またその精度をさらに向上させるために仮想的なグループの数を 8 つに増加させる検討なども行われている[3]。
一方、核分裂生成核種の崩壊、生成過程においては、遅発ガンマ線の生成とその媒質内での輸送、事故時における揮発性 FP 核種のシステムからの漏洩など、従来動特性モデルでは取り扱いが困
難な物理現象が存在する。また６群モデルは、従来の原子炉設計、安全解析への適用を目的として開発されたため、想定内の対象、物理現象、物理パラメータを精度よく評価することは可能であるが、想定を超える対象や物理現象パラメータについて対応することはできない。
そのため本研究室では先行研究において全ての核分裂生成物を陽に取り扱う原子炉解析モデ
ル EFP モデル(Explicit Fission Products モデル)の開発を行った[4]。EFP モデルでは核分裂生成物の崩壊による他の核分裂生成物の生成および中性子放出を詳細に解析し追跡している。これにより EFP モデルでは核分裂生成物の崩壊チェーンを厳密に再現することが可能になり、より実現象に近い再現ができる。また EFP モデルでは遅発ガンマ線の放出も同様に計算すること

炉物理の研究 第74号 (2022年3月)ができる。６群モデルと EFP モデルの概略図を図１に示す。

図１：６群モデルと EFP モデルの概略図

1.2 核熱結合計算
原子炉の安全評価では、いくつかの過渡・事故条件を想定し、原子炉システムへの影響を定量的に評価して、その安全性を担保する。この想定の一つとして、反応度事故が挙げられる。反応度事故の具体的な例としては、制御棒が誤って引き抜かれた際に、即発臨界を越えるような正の反応度が瞬間的に印加されるようなものが考えられる。このような場合、即発臨界を越えるために原子炉出力は急激に上昇するが、それに伴う燃料温度の急激な上昇により負のドップラー効
果が印加され、原子炉は短時間で即発臨界未満となり、出力は急激に減少を始める。さらに、燃料ペレットから被覆管、冷却材への熱が伝わり、時間遅れで冷却材温度が上昇し、冷却材密度上
昇による負の反応度が印加され、さらに反応度が低下する。この例で明らかなように、反応度事
故においては負の反応度を印加する燃料ドップラー反応度が重要なフィードバック機構となるが、即発減速材加熱効果（Prompt Moderator Heating effect）も瞬時の負の反応度効果として寄与することが分かっている。PMH 効果とは、即発臨界を越えたことにより中性子束レベルが瞬時に上奏し、冷却材領域での弾性散乱反応率が急激に増加して、それにより冷却材温度が上昇して、負の反応度効果を与えるというものである。燃料ペレットからの熱伝導・熱伝達による熱移動とは異なり、反応度印加とほぼ同時に冷却材の温度上昇に寄与するため、燃料ドップラー反応度と同様に瞬時のフィードバック効果を与えるものである。実際の安全評価での PMH 効果の扱いは現時点では我々にとって不明であるが、無視されているか、発生エネルギーの一定割合が冷却材に瞬時に伝わると言った過程で考慮されているものと考えられる。
以上のように、反応度事故においては、種々の気候によって負の反応度が体系に印加されることになるが、核分裂反応の結果生成する核分裂生成物が崩壊した際に発生する遅発ガンマ線に

炉物理の研究 第74号 (2022年3月)

よるエネルギー輸送も、同様に負の反応度を体系に印加する機構の一つとして考えられる。遅発
ガンマ線は核分裂生成物の半減期に応じて発生率が変わるが、半減期が極端に短い核分裂生成物からの遅発ガンマ線は、PMH 効果と同様に、反応度事故直後に負の反応度を印加する機構になりうる。
以上を踏まえて、本研究では EFP モデルに基づく空間依存動特性解析を熱水力計算と結合した、いわゆる核熱結合計算を行うことで反応度事故を想定した核分裂生成物、中性子、ガンマ線の挙動を直接的に模擬した数値計算を行い、種々の負のフィードバック機構の影響を定量的に評価することを最終的な目的とする。具体的には、日本原子力研究所で開発された核熱水力結合
動特性解析コード EUREKA 相当の解析手法の核計算部分を EFP モデルに置換することを行い、反応度事故の解析を実施する。この最終的な目的のために、今回の研究では 6 群モデルでの核熱結合計算モジュールの開発を行った。

2．熱計算モジュールの開発
核熱結合計算に用いることを目的とした、熱計算モジュールの開発を行った。燃料棒の計算において必要となる被覆表面温度は、冷却材計算から導かれるパラメータである。同様に冷却材計
算において必要となる被覆管表面熱流束は、燃料棒の計算から導かれるパラメータである。燃料棒の計算モジュールと冷却材の計算モジュールそれぞれ単体では、被覆管表面温度と被覆管表
面熱流束をそれぞれ任意に設定していた。このモジュールでは燃料棒の計算と冷却材の計算を交互に行い、それぞれに必要な境界条件を受け渡しして計算を進めていった。概略図を図 2 に示す。

図 2：結合熱計算の概略図

そして過渡解析においての計算プロセスの概略図を図 3 に示す。𝑡 = 𝑡𝑖+1 における燃料棒と冷却材の計算は、𝑡 = 𝑡𝑖 での燃料棒と冷却材の計算結果を用いて、各計算に必要となる境界条件を渡した。そして作成したモジュールにおいて、反応度事故相当の問題を解いた。燃料棒と冷却材それぞれが定常の状態から、出力を 1 タイムステップである 0.1 秒間大きくして、その後元炉物理の研究 第74号 (2022年3月)の定常状態の出力に戻すような問題での燃料と冷却材の挙動の時間変化を調べた。つまり、ある時間に大きな反応度が印加されて、フィードバックによって直後に反応度が抑えられるような現象を再現した出力変化を設定した。

図 3：過渡解析の概略図

そして、結合熱計算モジュールを用いて反応度事故相当の出力上昇を考慮した問題を解いた。冷却材計算と燃料棒計算共に定常状態から、図 3 で説明したように出力を 1 タイムステップ分である 0.1 秒間大きくし、その後元の定常状態の出力に戻すような問題を解いた。解析体系を示す。
＊解析体系
・流路高さ：4.0m
・水力等価直径：0.01103m
・入口温度：553.15K
・入口圧力：14.0Mpa
・入口流量：0.3472kg/s
・軸方向ノード数：50
この解析体系における、結合熱計算での冷却材温度の過渡変化を図 4、燃料棒温度の過渡変化を図 5 に示す。冷却材に関しては出口温度、燃料棒に関しては軸方向と径方向の中心点の温度をそれぞれ図示した。出力変化の幅による結果の差を確認するために、出力を10 倍、30 倍、50 倍と増加させた 3 つのパターンを結果として掲載した。

炉物理の研究 第74号 (2022年3月)

図 4：燃料棒の中心点の温度の過渡変化

図 5：冷却材出口温度の過渡変化
燃料棒、冷却材共に出力上昇に伴って温度が急上昇し、その後ピークを迎えて元の定常状態に収束していく結果を得た。それぞれの温度は出力と境界条件によって決定されるので、出力をもとの値に戻せば温度も最終的に元の定常状態に収束していくが、その様子を
図 4 と図 5 から確認することが出来た。出力上昇に対応してまず燃料棒の温度が急上昇す
るが、冷却材に関しては熱の移動によるタイムラグがあるので、比較的に緩やかに温度が上昇する結果となった。また燃料棒中心の温度が出力をもとに戻した後も上昇し続けるのは、燃料棒の熱伝導率の変化が関係している。出力上昇に伴って燃料温度が急上昇するが、その際に燃料棒の熱伝導率が急減する。それによって表面からの熱流束が急減するた
め、その分使われなかった熱量が燃料棒の温度上昇に使われるために出力減少後も温度はある程度上昇し続ける。冷却材の温度上昇のピークのタイミングよりも燃料棒の温度上昇

炉物理の研究 第74号 (2022年3月)のピークのタイミングが遅いのはそのためである。以上の結果から、冷却材計算と燃料棒計算を結合した熱計算モジュールが適切に計算できていることを確認した。なお過渡変化の妥当性についての確認は行っていないので、その点に関してはさらなる検討が必要である。

3．6 群モデルによる核熱結合計算
遅発中性子 6 群の動特性計算モジュールと前章で述べた熱計算モジュールを結合した核
熱結合計算モジュールを作成し、反応度フィードバックを考慮して反応度をステップ入力
した問題に対する各パラメータの過渡変化を調べた。まず、核熱結合計算において、核計算
は熱計算によって求められる燃料棒温度や冷却材温度の反応度フィードバックを用いて計算され、熱計算は核計算によって求められる出力を用いて計算される。反応度フィードバッ
クを考慮した反応度は、燃料温度𝑇𝑓 と冷却材温度𝑇𝑐 を用いて計算した。各モジュール単体の計算では任意に設定していた出力履歴や反応度履歴等のパラメータを、各計算から受け渡して過渡計算を進めていった。概略図を図 6 に示す。

図 6：核熱結合計算の概略図
そして過渡解析における計算プロセスの概略図を図 7 に示す。前章で説明したように、𝑡 =𝑡𝑖+1 における燃料棒と冷却材の計算は、𝑡 = 𝑡𝑖 での燃料棒と冷却材の計算結果を用いて行った。
そして得られた𝑡 = 𝑡𝑖+1 における熱計算の結果を持ちいて𝑡 = 𝑡𝑖+1 における核計算を行い、得られた結果を𝑡 = 𝑡𝑖+2の熱計算に用いた。そして作成したモジュールを用いて、反応度事故
相当の問題を解いた。核計算が臨界定常状態かつ熱計算が定常上の状態から、反応度をステップ状に印加する問題での反応度や線出力、燃料や冷却材の過渡変化を調べた。つまり、それぞれのモジュール単体の計算で任意に設定していた反応度や出力履歴を一切設定せず、反応度をステップ状に印加するだけの問題を解いた。

炉物理の研究 第74号 (2022年3月)

図 7：核熱結合計算の計算プロセスの概略図
また、反応度フィードバックを考慮した核熱結合計算においては、燃料温度及び冷却材温度の反応度へのフィードバックを定量化する必要があった。この反応度係数を評価するた
めに、資料[5]より反応度係数を計算した。今回の計算で使用した反応度係数は、JENDL-3.3に基づく ORIGEN コード用のライブラリの作成のために準備された PWR ピンセル[5]に準じたモデルで計算することとした。資料をもとに計算した燃料と冷却材の反応度係数は式(1)、式(2)のとおりである。
・冷却材の反応度係数
Tcoolant1=585 K、Tfuel1=1200 K のとき kinf1=1.31966
Tcoolant2=595 K、Tfuel2=1200 K のとき kinf2=1.31526
1
1
(
−
)
𝑘𝑖𝑛𝑓1 𝑘𝑖𝑛𝑓2𝑑𝜌== −0.00025
𝑑𝑇𝑐 𝑇𝑐𝑜𝑜𝑙𝑎𝑛𝑡2 − 𝑇𝑐𝑜𝑜𝑙𝑎𝑛𝑡1(1)・燃料の反応度係数
Tcoolant1=585 K、Tfuel1=1200 K のとき kinf1=1.31966
Tcoolant2=585 K、Tfuel2=1300 K のとき kinf2=1.31672
1
1
(
−
)
𝑘𝑖𝑛𝑓1 𝑘𝑖𝑛𝑓2𝑑𝜌== −0.0000169
𝑑𝑇𝑓
𝑇𝑓𝑢𝑒𝑙2 − 𝑇𝑓𝑢𝑒𝑙1(2)

算出した反応度係数によるフィードバックを加えた時の核熱結合計算の解析体系と各パラメータの計算結果の図 8～図 11 を以下に示す。
〇解析体系

炉物理の研究 第74号 (2022年3月)

・冷却材入口温度：Tin = 280 ℃
・冷却材入口圧力：Pin = 1.4 × 107 Pa・冷却材入口流量： Gin = 0.5361 kg/s
・初期最大線出力： q ini = 364 W/cm
・即発中性子寿命： l = 1.50 × 10−6 𝑠・遅発中性子割合： 𝛽 = 0.0006
・燃料温度の反応度係数： 𝜕𝜌/𝜕𝑇𝑓 = −1.69 × 10−5 (Δ𝑘/𝑘)/K
・冷却材温度の反応度係数： 𝜕𝜌/𝜕𝑇𝑐 = −0.00025 (Δ𝑘/𝑘)/K・印加反応度：0.1𝛽 (Δ𝑘/𝑘)〇解析結果

図 8：反応度の時間変化

炉物理の研究 第74号 (2022年3月)

図 9：線出力の時間変化

図 10：冷却材出口温度の時間変化

炉物理の研究 第74号 (2022年3月)

図 11：燃料中心温度の時間変化
まず反応度については、ステップ状に印加された反応度に対して燃料温度と冷却材温度のフィードバックが効いて、ネットリアクテビティ―が減少していく様子を再現した。反応度係数は冷却材温度の方が絶対値が大きいが、燃料温度の方が温度変化が大きいためフィードバックは燃料温度の方が大きくなった。線出力に関しては、熱計算モジュール単体では出力が一時的に上昇してすぐ出力が下がる形に自分で設定していた。核熱結合計算でのフィードバックを考慮した計算による反応度の急減によって線出力も急減するか、または増加が抑えられる様子を再現することが出来た。ここで燃料中心の温度上昇と冷却材出口温
度の温度上昇を比較した際に、燃料中心の温度上昇の方が緩やかであった。燃料から冷却材に熱が伝わることを考えると、これは一見実際の現象と異なっているようにみえるが、熱伝導率の変化から説明できる。反応度の投入によって中性子数が急激に増加することで、燃料温度が急激に増加する。この時、燃料の熱伝導率は逆に急減するので、燃料中で熱が伝わりにくくなる。燃料中の発生熱量は表面からの熱の放出と燃料自身の温度上昇に使われるが、燃料中の熱伝導が悪くなることで熱が燃料自身の温度上昇により使われる。これによって燃料の径方向温度分布がより急になるように推移していく。この推移にタイムラグがある
ので、冷却材が温度上昇のピークを迎えるより遅く、燃料中心が温度上昇のピークをとるという現象になった。
以上の結果より、6 群モデルの動特性計算モジュールによる核熱結合計算が適切に行うことが出来ていると結論付けた。なおそれぞれのパラメータの過渡変化の妥当性については検証を行っていないので、それについての検証はこれから行う必要がある。

炉物理の研究 第74号 (2022年3月)

4．結論
核熱結合計算に用いる熱計算モジュールを開発し、6 群モデルの動特性計算モジュールとの核熱結合計算モジュールで反応度事故相当の問題の解析を行った。得られた結果より、6群モデルの動特性計算モジュールによる核熱結合計算が適切に行うことが出来ていると結論付けた。なおそれぞれのパラメータの過渡変化の妥当性については検証を行っていない
ので、それについての検証は今後行う必要がある。

参考文献[1] W.M.Stacy, Jr., “Space-time nuclear reactor kinetics,” Academic Press, New York andLondon, (1969).
[2] G. R. Keepin, T. F. Wimett, R. K. Ziegler, “Delayed neutrons from fissionable isotopesof uranium, plutonium, and thorium,” Phys. Rev., 107, 1044-1049 (1957).
[3] G.D.Spriggs, J.M.Campbell, “An 8-group delayed neutron model based on a consistent
set of half-lives,” Prog. Nucl. Energy, 41, 223-251 (2002).
[4] K.Katagiri, G.Chiba, “Spatially-dependent nuclear reactor kinetics calculations with theexplicit fission product model," Ann. Nucl. Energy, 133, 202-208 (2019).
[5] J. Katakura, et al., JAERI-Data/Code 2004-015 (2004).