汎用炉物理解析コードシステム CBZ の概要 1
2020/9/21 千葉 豪
概 要
汎用炉物理解析コードシステム CBZ は、核分裂性物質を含む系の核分裂連鎖反応、核燃料の燃焼、放射線の輸送と減衰といった、原子炉内とその周囲における中性子、ガンマ線の輸送、それらと原子核、原子との相互作用に関わる物理現象を数値的に模擬するためのコードシステムである。日本原子力研究開発機構で開発されたコードシステム CBG をベースとして、2012 年 4 月より北大・原子炉工学研究室で開発が行われている。
プログラムはコンピュータ言語 C++で記述されており、原子炉物理の計算に関連する情報（原子炉や燃料集合体の構造を表現するための幾何形状や、原子炉を構成する媒質といったものに加えて、計算手法や収束条件なども含む）は全て C++の「クラス」として定義される。中性子、ガンマ線の輸送方程式、拡散方程式を解くソルバーが複数実装されており、それらソルバーはそれぞれ関連する複数のクラスで構成されている。また、ソルバー間のデータ（例えば中性子反応断面積データ）のやり取りはクラスのインスタンスを通して行われることか
ら、複数のソルバーを組み合わせた多種多様な炉物理計算が容易に実現可能であることが特徴である。

1多群断面積ライブラリ CBZLIBCBZ は決定論的手法に基づく計算コードシステムであるため、エネルギーを離散化して取り扱う必要がある。
従って、共鳴自己遮蔽効果を考慮した実効断面積を計算することになるが、その際に用いる多群ライブラリとして
CBZLIB が実装されている。CBZLIB はいわゆる Bondarenko 型の炉定数であり、無限希釈断面積、非等方散乱の高次のルジャンドル係数を含む散乱マトリクスに加えて、共鳴自己遮蔽因子を有している。ライブラリは NJOY-99
コード等を用いて評価済み核データファイルを処理することにより作成され、エネルギー群構造は任意である。

CBZ を用いた一般的な解析では、軽水炉を対象とした場合には 107 群（SRAC-2K6 コードのライブラリと同一の群構造）、高速炉の場合には 70 群（JAERI Fast-set と同一の群構造）、もしくは 175 群のライブラリ 2 を使用する。
CBZLIB を用いた実効断面積計算は自己遮蔽因子を用いて行う。一般的には自己遮蔽因子は NR 近似 3 により得た中性子束を用いて計算される。ただし、軽水炉解析用の CBZLIB には、中性子の減速をより厳密に考慮するため、共鳴自己遮蔽因子の計算では、荷重関数として一点炉詳細群減速方程式を解いて得られる中性子束エネルギースペクトルを用いているものもある。また、核種間の共鳴干渉効果を考慮するため、重核種に対しては多重R 因子法 [1] 4 を適用し、FP に対しては U-238 との共鳴干渉効果を考慮した荷重関数（U-238 の巨大共鳴による中性子束エネルギースペクトルの歪みを一意的に考慮したもの）を無限希釈断面積作成時に用いる工夫 [2] を行っているものもある。さらに、ピンセル体系におけるカレント重みの実効全断面積の計算精度を向上させるための補正因子を、U-238 と Pu-242 のいくつかのエネルギー群に対して導入している。これらはまとめて「AdvancedBondarenko method」と呼称され、文献 [3] に詳細が述べられている。
CBZLIB では、個々の核種のデータがそれぞれ単一のテキストファイルとして与えられている。この具体的な解説は別テキスト「CBZLIB の概要と利用マニュアル」に記載されているので、そちらを参照のこと。

2実効群断面積の計算（共鳴自己遮蔽効果の処理）
均質媒質における実効群断面積の計算は、個々の共鳴核種について背景断面積を求め、背景断面積、温度と R因子をパラメータとした自己遮蔽因子テーブルから自己遮蔽因子を内挿することで行う。一方、非均質格子系における共鳴自己遮蔽効果の計算は等価原理に基づいて行う。非均質格子系での背景断面積の評価は、燃料ピンセル体系に対しては以下の方法で行う。まず燃料からの脱出確率は Bell 因子を用いた一項有理近似で記述する [4]。ここで、Bell 因子は軽水炉の燃料ピン系に対してエネルギー群毎に最適化された値を用いている [1]。また、格子効果（複数の燃料棒の配列効果）はエネルギー群依存の Dancoﬀ 係数により考慮する [4]。Dancoﬀ 係数は杉村、山本の方法を用いて燃料集合体での評価も行うことができる [5]。高速炉解析で馴染みの深い東捻の方法 [6] も利用可能である。
1 /Document/CBG Manual/Summary2 この 175 群ライブラリの群構造は、主に中性子遮蔽計算のために開発された VITAMIN-J と呼ばれるものを僅かに変更したものであり、超小型高速炉の解析など、中性子の漏洩寄与が大きく 70 群構造では不十分となる場合にこの群構造を用いている。
3 中性子束 ϕ(E) が全断面積 Σ(E) に反比例するという近似である。
4 ある共鳴核種の実効断面積の計算では、通常、それ以外の核種は非共鳴核種であると仮定される。しかし実際には種々の異なる共鳴核種が混在しており、その干渉効果（共鳴干渉効果）が重要となる。R 因子は着目共鳴核種と共鳴干渉を考慮する核種の数密度の比として定義される。R 因子について、いくつかの値で予め共鳴干渉効果を考慮した自己遮蔽因子を計算しておくことで、実効断面積計算において特定の核種との共鳴干渉効果を近似的に考慮することが可能となる。

1

3

中性子、ガンマ線の輸送・拡散ソルバー
CBZ には中性子、ガンマ線の輸送方程式、拡散方程式を解くための多様なソルバーが実装されている。これらのうち、ただひとつのソルバーを除いて、全てが同一の基底クラス（GeneralSystem クラス）を継承した構造となっている。これは、それらソルバーがべき乗法による固有値計算や固定源計算の枠組みを共有できるからである。また、拡散方程式、輸送方程式いずれにおいても現れる物理量であるスカラー中性子束（ガンマ線束）を用いた種々の計算（反応率分布の計算や摂動計算の一部）はこの基底クラスに実装されている。随伴方程式の計算も殆どのソルバーで可能となっている。
拡散方程式を解くためのソルバーとしては、PLOS、DHEX がある。これらは有限体積法に基づいたコードであり、PLOS はデカルト座標系、円筒体系、球体系に、DHEX は Hexagonal-Z 体系に適用することが出来る。DHEX は、六角メッシュに加えて、六角メッシュを 6 等分した三角メッシュ体系での計算も可能となっている。内部反復の加速には SOR 法を、外部反復の加速には粗メッシュ有限差分（CMFD）加速 [7] を採用している。これらに加えて、階層領域分割型境界要素法 [8] に基づくソルバー ABEMIE も実装されている。ABEMIE は任意の線分で構成される二次元体系の計算が可能である。
輸送方程式を解くためのソルバーとしては、SNR、SNRZ、SNT がある。これらはいずれも離散座標法に基づいており、それぞれ球体系、円筒体系、デカルト座標系に適用する。離散座標法は中性子の角度方向を離散的に扱う方法であり、Sn 法とも呼ばれる [9, 10]。いずれのソルバーも、内側反復には拡散合成法 [11] を、外側反復にはCMFD 法を採用し、計算の高速化を図っている。
また、燃料ピンセルや燃料集合体など、二次元平面上の複雑形状における輸送方程式を扱うためのソルバーとしては PJI、MEC がある。PJI は衝突確率法 [10]、MEC は特性曲線法（Characteristic 法）に基づくソルバーであり、Ray-tracing[12, 13] により衝突確率やパスラインの計算を行う。なお、Ray-tracing の計算は PJI、MEC ともに同一のモジュールを利用している。

Simplified P3(SP3) 方程式 [14, 15] に対する計算も可能である。デカルト座標系については PLOS の固定源機能を組み合わせて用いることにより実現している。また、エネルギー１群に限定しているが、ABEMIE ソルバーもSP3 方程式を解くことが出来る [16]。
また、通常はベクトル形式で考慮する核分裂スペクトルに対して、入射中性子エネルギー依存性を陽に取り込んだ計算（核分裂スペクトルを行列形式で考慮する計算）も可能である [17]。

4燃焼計算
核燃料に含まれる種々の原子核は、原子炉の運転に伴いその数密度が時々刻々変化する。この計算はいわゆる「燃焼方程式」を解くことで実現されるが、CBZ にはこのためのモジュールとして Burnup クラスが実装されている。燃焼方程式の数値解法としては行列指数法を採用し、行列指数は Mini-Max Polynomial Approximation(MMPA)[18, 19]、もしくは Chebyshev 有理近似法（Chebyshev Rational Approximation Method, CRAM）[20] により評価している（これらの方法では Krylov 部分空間法 [21] で解けない短半減期核種を含むチェーンでも対応可能である）。また、燃焼チェーンを作成するモジュール BurnupChainGenerator が実装されており、評価済みの崩壊データファイル、核分裂収率データファイルから、ユーザーが必要と判断した核種で構成される任意の燃焼チェーンを作成することが可能である。なお、CBZ を用いた一般的な軽水炉の燃焼計算では、燃焼中の反応度が精度良く計算できるように独自のアルゴリズムにより選定された 138 の FP からなるチェーン [22] を用いている。Burnup クラス
は入力された核種数密度、中性子束情報（エネルギースペクトル、レベル）、燃焼期間から、燃焼後の数密度を計算する。中性子束の情報は、中性子輸送・拡散ソルバーから引き渡すことになる。

Burnup モジュールは一群断面積に基づく燃焼計算モジュールであるが、一般的な燃焼計算では、燃焼領域の中性子束のエネルギースペクトルを計算し、それに基づいて一群断面積を求め燃焼方程式を解く、という手続きがとられる。CBZ では、単一のピンセル体系もしくは複数のピンセルからなる体系の燃焼計算を行うための基底クラス
GeneralBurner が実装されており、それを継承する形で単ピンセル体系向けのモジュール Burner と複数ピンセル体系向けのモジュール MulticellBurner が実装されている。Burner は燃料、被覆管、冷却材の三領域からなる正方及び六角格子体系について、MulticellBurner は正方格子体系の複数のピンセルからなる体系について、それぞれ燃焼計算を行う。中性子束計算は Burner は衝突確率法に基づく PJI モジュール、MulticellBurner は MOCに基づく MEC モジュールでそれぞれ行う。いずれのモジュールでも、出力運転中の計算とともに、原子炉運転後の2冷却計算も行うことが可能であり、燃焼（冷却）期間に依存した各種原子核の数密度や放射能、崩壊熱といった、燃焼に伴う諸量を計算することが出来る。MulticellBurner モジュールでは、可燃性毒物を含んだ体系の計算において時間に対する離散化誤差を低減する目的で、Predictor-Corrector 法（PC 法）、さらには Optimally-Weighted

PC 法（OWPC 法）[23, 24] が実装されている。
また、高速炉や加速器駆動システム（Accelerator-Driven System、ADS）の全炉心を二次元円筒体系に模擬した燃焼計算も FRBurnerRZ モジュールにより可能となっている [25]。このモジュールでは燃料交換も簡略的に考慮する。

5

感度計算、不確かさ評価
実効増倍率の核データに対する感度係数は中性子束と随伴中性子束から計算出来るため、随伴中性子束を計算できるソルバーでは実効増倍率の感度係数の計算が可能である（拡散/輸送、ピンセル/燃料集合体/全炉心）。反応度の感度係数は実効増倍率の感度係数の差をとることにより求める。反応率比の感度係数では一般化随伴中性子束の計算 [26] が必要となるが、基本的には固定源随伴問題となるため、GeneralSystem クラスに実装されているメソッドで対応可能である。実効遅発中性子割合、中性子生成時間等、定義に随伴中性子束が現れるパラメー
タの感度については、Modified k-ratio 法に基づく方法で計算を行う [27]。
また、単一のピンセル体系および燃料集合体を含む複数のピンセル体系における燃焼後の核特性（数密度、実効増倍率、反応率比）の核データに対する感度を計算するメソッドも実装されている [2, 28, 29, 30]。PC 法、OWPC法を用いた燃焼計算についても対応が可能となっている [31]。
感度係数は SensitivityData クラスで定義する。感度の画面への出力など、いくつかの便利なメソッドがSensitivityData クラスに実装されている。
核データに起因する不確かさの定量化は、NJOY-99 の ERRORR モジュールが計算した多群共分散ファイルをCBZ 固有のフォーマットに変換した多群共分散ライブラリと SensitivityData クラスのインスタンスを用いて、UNC モジュールが行う。断面積の不確かさのみならず、原子核の崩壊データや核分裂収率データの不確かさも考慮することが可能である。

6

中性子・ガンマ線結合計算
ガンマ線の輸送計算を行う場合、中性子核反応によるガンマ線生成データ、ガンマ線輸送データ（光子と原子の反応確率データ）については、NJOY-99 で処理を行った CBZ 固有の多群ライブラリを用いる。ガンマ線の輸送計算は、はじめに中性子輸送計算を行い、それに基づいてガンマ線源を求め、ガンマ線に対する固定源計算を行うという手続きをとる。燃焼計算と組み合わせた遅発ガンマ線放出の時間依存性を正確に取り扱った計算も実現可能である [32]。

参考文献[1] G. Chiba, “CBGLIB: a multi-group neutron library for precise neutronics simulations,” Proc. of the 2007
symposium on Nuclear Data, JAEA-Conf 2008-008, p.114 (2008).
[2] G. Chiba, K. Okumura, A. Oizumi, M. Saito, “Sensitivity analysis of fission product concentrations forlight water reactor burned fuel,” J. Nucl. Sci. Technol., 47, p.652 (2010).
[3] G. Chiba, T. Narabayashi, “Advanced Bondarenko method for resonance self-shielding calculations indeterministic reactor physics code system CBZ,” Ann. Nucl. Energy, 96, p.277 (2016).
[4] R. J. J. Stamm’ler, M. J. Abbate, Methods of steady state reactor physics in nuclear design, Academic,London, New York (1983).
[5] N. Sugimura, A. Yamamoto, “Evaluation of Dancoﬀ factors in complicated geometry using the method ofcharacteristics,” J. Nucl. Sci. Technol., 43, p.1182 (2006).
3

[6] T. Tone, “A numerical study of heterogeneity eﬀects in fast reactor critical assemblies,” J. Nucl. Sci.
Technol., 12, p.467 (1975).
[7] K. S. Smith, “Nodal method storage reduction by non-linear iteration,” Trans. Am. Nucl. Soc., 44, p.265(1984).
[8] G. Chiba, M. Tsuji, Y. Shimazu, “A hierarchical domain decomposition boundary element method with ahigher order polynomial expansion for solving 2-D multi-region neutron diﬀusion equations,” Ann. Nucl.
Energy, 28, p.895 (2001).
[9] E. E. Lewis, W. F. Miller Jr., Computational methods of neutron transport, American Nuclear Society(1993).
[10] 小林啓祐、「原子炉物理」、コロナ社（1996）.
[11] M. L. Adams, E. W. Larsen, “Fast iterative methods for discrete-ordinates particle transport calculations,”
Prog. Nucl. Energy, 40, p.3 (2002).
[12] K. Tsuchihashi, “Collision probability in two-dimensional lattice by ray-trace method and its applications
to cell calculations,” JAERI-M 85-034 (1985).
[13] R. Sanchez, L. Mao, S. Santandrea, “Treatment of boundary conditions in trajectory-based deterministictransport methods,” Nucl. Sci. Eng., 140, p.23 (2002).
[14] E. W. Larsen, J. E. Morel, J. M. McGhee, “Asymptotic derivation of the multigroup P1 and simplified PNequations with anisotropic scattering,” Nucl. Sci. Eng., 123, p.328 (1996).
[15] P. S. Brantley, E. W. Larsen, “The simplified P3 approximation,” Nucl. Sci. Eng., 134, p.1 (2000).
[16] G. Chiba, “Application of the hierarchical domain decomposition boundary element method to the simplified P3 equation,” Ann. Nucl. Energy, 38, p.1033 (2011).
[17] G. Chiba, “Criticality calculations with fission spectrum matrix,” Proc. of the 2008 Symp. on Nuclear
Data, JAEA-Conf 2009-004, p.107 (2009).
[18] Y. Kawamoto, G. Chiba, M. Tsuji, T. Narabayashi, “Numerical solution of matrix exponential in burn-up
equation using mini-max polynomial approximation,” Ann. Nucl. Energy, 80, p.219 (2015).
[19] G. Chiba, Y. Ohoka, K. Yamamoto, H. Nagano, “Revisiting mini-max polynomial approximation method
for nuclear fuel depletion calculations,” Proc. Int. Conf. on Math. Comp. Methods Applied to Nucl. Sci.
Eng., 2019 Aug 25-29, Portland (OR), (2019).
[20] M. Pusa, J. Leppanen, “Computing the matrix exponential in burnup calculations,” Nucl. Sci. Eng., 164,p.140 (2010).
[21] T. Yamamoto, et al., “Numerical Solution of stiﬀ burnup equation with short half lived nuclides by theKrylov subspace method,” J. Nucl. Sci. Technol., 44, p.147 (2007).
[22] G. Chiba, M. Tsuji, T. Narabayashi, Y. Ohoka, T. Ushio, “Important fission products identification method
for simplified burnup chain construction,” J. Nucl. Sci. Technol., 52, p.953 (2015).
[23] S. Okumura, G. Chiba, “Development of nuclear fuel depletion calculation capability for LWR fuel assemblyin reactor physics code system CBZ,” Proc. of Reactor Physics Asia 2017 Conference (RPHA2017), (2017).
[24] J. Sasuga, G. Chiba, Y. Ohoka, K. Yamamoto, H. Nagano, “Improvement of optimally-weighted predictorcorrector methjod for nuclear fuel burnup calculations,” Proc. of Reactor Physics Asia 2019 Conference(RPHA2019), (2019).

4

[25] J. Fan, G. Chiba, “Verification of FRBurner module of CBZ code system based on OECD/NEA benchmark
report,” Proc. of Reactor Physics Asia 2019 Conference (RPHA2019), (2019).
[26] 三谷浩、黒井英雄、「積分量の感度係数と一般化摂動法」、JAERI-M 4760 (1972).
[27] G. Chiba, M. Tsuji, K. Sugiyama, T. Narabayashi, “JENDL-4.0 benchmarking for eﬀective delayed neutron
fraction of fast neutron systems,” J. Nucl. Sci. Technol., 48, p.1471 (2011).
[28] G. Chiba, M. Tsuji, T. Narabayashi, “Uncertainty quantification of neutronics parameters of light water
reactor fuel cells with JENDL-4.0 covariance data,” J. Nucl. Sci. Technol., 50, p.751 (2013).
[29] G. Chiba, T. Narabayashi, “Development of a fuel depletion sensitivity calculation module for multi-cellproblems in a deterministic reactor physics code system CBZ,” Ann. Nucl. Energy, 96, p.313 (2016).
[30] G. Chiba, “Burnup sensitivity calculations with CBZ for light water reactor assembly problems,” Proc. ofInt. Conf. on Physics of Reactors (PHYSOR2020), Cambridge, UK (2020, cencelled).
[31] G. Chiba, “Perturbation theory for nuclear fuel depletion calculations with predictor-corrector method,”J. Nucl. Sci. Technol., 55, p.290 (2018).
[32] G. Chiba, M. Tsuji, T. Narabayashi, “Photon transport eﬀect on intra subassembly thermal power distribution in fast reactor,” Ann. Nucl. Energy, 65, p.41 (2013).

5