1

汎用炉物理解析コードシステム CBZ のチュートリアル
千葉豪

1

はじめに
原子核に中性子が衝突すると、ある確率のもとで核分裂反応が起こり、原子核は通常２つに分裂します。

その分裂した原子核（核分裂片）は莫大な運動エネルギー（約 167 MeV）を持っており、それを最終的に
電気エネルギーに変換して利用するのが原子力発電です1 。
核分裂反応を起こす核燃料は原子炉容器に格納されます。原子炉容器内には種々の機器が配置されていますが、その中の核燃料が設置される周辺は炉心（Reactor core）と呼ばれ、文字通り原子力発電所の心臓部にあたります。
原子力発電では、核分裂反応を持続的に行わせるため、反応を生じさせる中性子の数を制御しています。
すなわち、核分裂反応を適切に制御するためには、我々は炉心内での中性子の振る舞いを知る必要があります。炉心内の中性子の振る舞いは Boltzman の輸送方程式と呼ばれる式に従いますので、この方程式の解を得なければなりません。この方程式を解いて中性子の振る舞いを知ることは炉心解析、炉心の核特性解析2と呼ばれ、原子力工学の重要な分野のひとつとなっています。
炉心解析を行うための種々のコードシステム（プログラム）が世界各国で開発されています。最も有名なものとしては STUDSVIK 社の CASMO/SIMULATE が挙げられるでしょう。このコードシステムは世界各国において実機軽水炉の炉心解析に広く用いられています。また、日本でも様々なものが開発されており、例えば日本原子力研究開発機構（JAEA）の SRAC コードや、原子燃料工業、原子力エンジニアリン
グが名古屋大学と共同で開発している AEGIS/SCOPE-2 などが挙げられます。
著者は JAEA 在籍時に、炉物理解析を学ぶために、独自にコードシステム CBG を開発し、それを用いていくつか論文を書きました。この CBG の開発経験を基に、北大に移ってから開発を始めたコードシステ
ムが CBZ です（このチュートリアルでは「CBZ」と「CBG」が混在していますが、どちらも同じようなものなので、あまり気にしないで下さい）。CBZ では、臨界計算や動特性計算、放射線の遮蔽計算、核燃料の燃焼計算など、多様な炉物理解析を行うことが出来ます。本チュートリアルは CBZ の使い方を学んでもら
うためのものです。

1 厳密には、核分裂片の運動エネルギーに加えて、核分裂反応とほぼ同時に発生するガンマ線や、不安定な核分裂片がベータ崩壊
した際に発生するベータ線、ガンマ線のエネルギーも利用されます。
2 核特性というのは、中性子と原子核との反応に関連した原子炉の特性のことを示しており、例えば臨界量や出力分布、反応度係数などが挙げられます。

2. インストール

2

2

インストール
CBZ では、ユーザは CBZ の種々のモジュールを利用して自身の目的に沿った計算を行います。具体的には、ユーザは C++を用いてプログラムを記述し、それをコンパイルすることで計算を行います。このコンパイルの際には、CBZ の種々のモジュールをコンパイルすることにより作成されたライブラリファイルである「libCBG.a」をインクルードする必要があります。この関係を Fig. 1 に示します。CBZ のソースプロ
グラムが図中左側の「CBZ source」に対応し、ユーザーが作成するプログラムが右側の「Main programs

for users」になります。ここでは、CBZ のソースプログラムのコンパイルを行います。ユーザが自身の目的に応じて記述する C++プログラムのコンパイルについては、これ以降の章で行うことになります。

Fig. 1: CBZ を用いた計算の構造パッケージ cbz.xxxxyyzz.tgz を解凍して展開する（xxxxyyzz は日付を示す）3 と、ディレクトリ CBZが生成され、その下に以下のディレクトリが作成されます4 。

 CBG
 CBGCAL CBGLIB CBGLIB BURN

CBG ディレクトリの下には src ディレクトリがあるはずです（dist ディレクトリもありますが、これについては後述します）。この src ディレクトリには、CBZ のソースプログラムが置かれています。一方、

CBGCAL ディレクトリは CBZ を用いていろいろな計算を行うためのデータを置くためのもので、配布パッ
ケージには Takeda.tutorial、lanl.tutorial、mza.tutorial 等のディレクトリがあるかと思います。また、CBGLIB ディレクトリは CBZ のための断面積（核データ5 ）ライブラリ CBZLIB を置くためのもの、

CBGLIB BURN ディレクトリは燃焼計算を行う際に必要となるデータファイルを置くためのものです。
それでは CBZ をコンパイルしましょう。まずは CBG/src ディレクトリに移動して下さい6 。配布されたパッケージでは既にコンパイルされている可能性がありますので、まずは make clean と打ち込んで、既にコンパイルされている場合に作成されているファイルを消去します。その後、コンパイルを行うため、makeと打ち込んでください。コンパイルが始まり、数分後くらいに libCBG.a が生成されます。ls *.a と打ち込んで、このファイルが生成されていることを確認して下さい。
3 tar

xvfz cbz.xxxxyyzz.tgz と打ち込む。

4 CBZ ディレクトリの下には階層構造でディレクトリが構成され、さまざまなファイルが置かれています。CBZ ディレクトリに含まれる全てのデータをどこかにもっていくためには、これらをひとつのファイルにまとめるのが便利です。ディレクトリ構造をひとつのファイルにまとめることを「tar する」といい、まとめられたファイルは tar ファイルと呼ばれます。例えば、CBZ ディレクト
リを tar ファイルにまとめるためには、tar cvf cbz.tar CBZ と打ち込めば、CBZ ディレクトリが cbz.tar というひとつのファイルに変換されます。また、tar されたファイルは一般にサイズが大きくなるため、圧縮してあげます。そのコマンドが gzip cbz.tar
であり、これを行うことにより、cbz.tar.gz ファイルが作成されます。この、tar して圧縮する一連の動作を一度に行うコマンドがtar cvfz cbz.tgz CBZ であり、その逆が tar xvfz cbz.tgz となります。
5 核データとは、中性子と原子核との反応を記述するためのデータです。詳細は後ほど述べます。
6 cd CBZ/CBG/src と打ち込む。

2. インストール

3

なお、gfortran がコンパイルされていない環境ではエラーで落ちますので、この場合は、gfortran を自らの計算機環境にインストールして下さい（インストールの方法が分からない場合には周りの学生や教員に聞いて下さい）。もしくは、cp Makefile wo fortran Makefile と打ち込んで、gfortran が不要な条件でMake を行って下さい。
また、gfortran とは別な要因でコンパイルに失敗する場合があります。この場合には、src ディレクトリ
内の Numeric.h の 12 行目（以下参照）のコメントアウトを外したうえで、再びコンパイルしてみて下さい。これでもコンパイルできない場合には、千葉に相談して下さい。

Listing 1: Numeric.h の冒頭部分の抜粋
1
2
3
4
5
6
7
8
9
10
11
12
13
14

#i f n d e f PIDEF#d e f i n e PIDEF#include <s t r i n g >#include <v e c t o r >#include <math . h>
typedef double r e a l ;typedef f l o a t r e a l 4 ;typedef double r e a l 8 ;//#d e f i n e

OLD COMPILER ENV

const r e a l PI = 3 . 1 4 1 5 9 2 6 5 3 5 8 9 7 9 3 2 3 8 4 6 2 6 4 3 3 ;以上でインストールは終了です。
なお、CBZ のソース（src ディレクトリ内のプログラムの一部）を新しいものに置換（上書き）した場
合や、ある別な PC でインストールが完了した CBZ ディレクトリをそれとは異なる PC にコピーして利用する場合などには、改めてここで述べられた手順でコンパイルを行い、libCBG.a を作成する必要があります。特に、「新しいソースを使って欲しい」という連絡があった場合には、新しいソースの tgz ファイルが提供されることになります。このファイルを解凍すると src ディレクトリが生成されますので、この作業を行う前に、まずはそれまで使っていた src ディレクトリの名前を変更し（例えば src old など）、その後、CBG ディレクトリで新しいソースの tgz ファイルを解凍し、再コンパイルを行うとよいでしょう。この最後で述べた部分は、今の段階では全く関係ありませんが、今後、CBZ を使っていく際には非常に重要な事項となりますので、覚えておいて下さい。

3. Takeda ベンチマークの計算

3

4

Takeda ベンチマークの計算
それでは CBZ を用いて、中性子輸送7 のベンチマーク問題として世界的に広く用いられている Takeda

ベンチマークの計算をしてみましょう8 。Takeda ベンチマーク問題には４つの体系が用意されていますが、ここでは最も簡単な体系１を計算します。体系１の概略図を Fig. 2 に示します。

Fig. 2: Takeda ベンチマーク問題の体系１（左が XY 平面、右が XZ（YZ）平面）
Fig. 2 の左が上から見た図になります。全体の 1/4 が示されており、中央部に「炉心領域（Core）」、外周部に「反射体領域（Reﬂector）」が配置されていることが分かります。また、反射体領域の一部が黒塗りされていますが、この部分は制御棒が挿入される位置で、制御棒が挿入されていない状態では「制御棒フォ
ロワ9 」（実質的には空気）となります。
右は横から見た図になっており、制御棒（もしくはフォロワ）領域が縦に延びていることが分かります。
また、軸方向についても炉心領域が反射体領域に囲まれていることが分かります。
それでは、CBGCAL/Takeda.tutorial ディレクトリの main.takeda1.dif 1.cxx をエディタで開いて下さい。
まず、Geometry data と書かれたコメント行が（20 行目くらいに）あると思います。ここから下の部分
では、計算する体系の形状データを定義しています。今回計算するのは正方格子状の体系で、X、Y、Z 方
向をそれぞれ 5 領域（ブロック）に区切っています。各々のブロックの各辺の長さは 5cm となっており（配列 xl、yl、zl で指定）、計算ではそれら各辺をさらに 5 分割しています（配列 xm、ym、zm で指定）。各辺が 25 に分割されることになるため、この体系の全計算メッシュ数は 25 × 25 × 25 = 15, 625 となります。
また、5 × 5 × 5 = 125 ブロックのそれぞれについて、どの物質（マテリアル）が割り当てられるかを配列

mat で定義します。ここでは、0 が燃料、1 が反射体、2 が制御棒、3 がフォロワとなります。なお Takeda
ベンチマークの体系１では、制御棒挿入状態と制御棒引き抜き状態のふたつの状態について評価することになっています。この入力ファイルでは、制御棒挿入位置に物質 3 が割り振られているため、制御棒引き抜き状態に対応します。
体系の物質配置の定義の仕方がいまいちピンとこない場合には、ホームページの「CBZ」の項に「CBZの XYZ 座標系の炉心体系情報に関するクラス」という文書があるので、その 6.1 節「媒質マップと物質配置、境界条件の関係」を参考にして下さい。
次に、Material data と書かれたコメント行を（70 行目くらいに）見つけられることと思います。ここから下の部分では、物質のデータ、具体的には中性子と原子核の反応断面積を定義しています。なお、この
7 「輸送」というと、何か物を運ぶようなことが連想されますが、媒質内の中性子の移動を英語では Neutron transport と呼ぶため、これを和訳して中性子輸送となります。いずれ慣れます。
8 「Takeda ベンチマーク」とは、元大阪大学教授（現福井大学教授）の竹田敏一先生が中性子輸送方程式を数値的に解くプログラ
ムの妥当性を検証するために作成したベンチマーク問題です。1991 年に経済協力開発機構 · 原子力機関（OECD/NEA）が世界各国の研究者を召集し、このベンチマーク問題を用いて、参加者それぞれが開発、使用している計算コードの比較が行われました。
9 制御棒の「通り道」ということで、このような名前となっています。

3. Takeda ベンチマークの計算

5問題ではエネルギーを２つに分けています。ちなみに、分割されたエネルギー領域を群（group）と呼び、この例では「エネルギー群数は 2」という言い方となります。1 群が高速中性子エネルギーに、2 群が熱中
性子エネルギーに、それぞれ対応します。
最後に、それら定義した形状データ、物質データを用いて、中性子拡散方程式を解きます。なお、中性子の振る舞いは中性子輸送方程式によって記述されますが、輸送方程式に近似を導入することで拡散方程式が得られます。輸送方程式においては、中性子束は位置、エネルギー、運動方向に依存する（角度中性子束を扱う）一方、拡散方程式では運動方向の依存性を近似的に扱う（スカラー中性子束のみを扱う）ため、方程式に中性子束の運動方向依存性が現れず、計算負荷が小さくなります。ここでは、中性子拡散方程式を数値的に解くコード PLOS を呼び出して、処理を行わせています。
実際の計算は以下の通りに行います。まず cp main.takeda1.dif 1.cxx main.cxx と打ち込んで、

main.takeda1.dif 1.cxx ファイルを main.cxx にコピーして下さい。その後、make と打ち込んでコンパイルをして下さい10 。その結果、a.out という名前の実行形式（ロードモジュール）ファイルが生成されますので、確認して下さい。このファイルを確認できたら./a.out と打ち込んで、計算を実行して下さい。
以下のように処理結果がずらずらと出力されることと思います。

Listing 2: Takeda ベンチマークの計算例
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

#*** T o t a l Mesh : 15625#*************************************************#* System CBG#*
Power i t e r a t i o n#*
S o l v e r : PLOS#*
A c c e l e r a t i o n : none#*
Forward c a l c u l a t i o n#*
C o n v e r g e n c e c o n d i t i o n ( k e f f ) : 1 e =05#*( flux ): 0.0001#*( source ) : 0.0001#*************************************************#* I t :K eff: Err i n k : Err i n: Err i n#* e r ::: flux: source#*************************************************#
0:0.906483
9 . 3 5 2 e =02
2 . 2 3 3 e+00
3 . 7 0 0 e =01#
1:0.917541
1 . 2 2 0 e =02
1 . 0 1 8 e =01
7 . 7 2 3 e =02#2:0.925198
8 . 3 4 4 e =03
6 . 6 3 0 e =02
5 . 3 4 6 e =02#3:0.927640
2 . 6 4 0 e =03
1 . 5 5 2 e =02
1 . 3 7 6 e =02#4:0.927664
2 . 5 6 5 e =05
4 . 6 9 9 e =02
2 . 7 8 7 e =02#5:0.928331
7 . 1 9 8 e =04
4 . 4 2 8 e =03
2 . 9 5 2 e =03#
6:0.927417
9 . 8 4 9 e =04
1 . 4 2 4 e =03
2 . 5 6 3 e =02#7:0.927214
2 . 1 9 2 e =04
4 . 2 0 4 e =03
1 . 6 0 6 e =02#8:0.927261
5 . 1 4 6 e =05
2 . 9 9 1 e =03
7 . 0 6 6 e =03#9:0.927328
7 . 1 2 8 e =05
3 . 1 5 5 e =03
2 . 6 9 0 e =03# 10:0.9273906 . 7 4 4 e =05
6 . 7 1 0 e =04
4 . 4 0 9 e =04# 11:0.927305
9 . 2 2 1 e =05
2 . 3 3 9 e =04
2 . 5 8 0 e =03# 12:0.927291
1 . 4 9 2 e =05
4 . 3 2 4 e =04
1 . 4 4 0 e =03# 13:0.927291
6 . 5 2 0 e =07
8 . 5 5 8 e =04
7 . 0 4 4 e =04# 14:0.927304
1 . 4 1 4 e =05
1 . 4 6 0 e =04
1 . 0 2 2 e =04# 15:
0.927295
9 . 9 2 4 e =06
4 . 2 4 9 e =04
3 . 3 0 9 e =04# 16:
0.927300
5 . 6 8 9 e =06
5 . 6 4 6 e =05
3 . 0 0 9 e =05

PLOS では拡散方程式を反復法によって解きますので、反復毎に中性子実効増倍率 keﬀ や中性子束（の空間、エネルギー）分布（の推定値）が更新されます。最終的に、実効増倍率や中性子束分布が一定値に収束したことを確認して、これらを計算結果とします。具体的には、今の反復回での結果と一つ前の反復回での結果の差異を評価して、それがある値よりも小さくなった場合に、収束と判断します。上の例に示されている情報はこの反復計算の過程であり、反復回数、その反復回で得られた実効増倍率、そして、実効増倍率、中性子束分布、核分裂中性子源分布についての前の反復回との相対差異（中性子束分布、中性子源分布についてはその最大値）が表示されています。最終的な実効増倍率の計算結果は最後の行から判断することができるでしょう。
数値計算には様々な誤差が含まれますが、代表的なものとして「空間メッシュの離散化誤差」が挙げられ
10 make により、同一のディレクトリにある main.cxx ファイルのコンパイルを行います（処理の詳細は Makefile ファイルで定義されています）。従って、計算したいファイルがある場合、本文中で述べられているように一旦それを main.cxx にコピーして makeすることになります。

3. Takeda ベンチマークの計算

6ます。今回の計算では一辺 5cm のブロックを 5 × 5 × 5 メッシュに分割していますが、この計算で得られるのはあくまで近似解であり、メッシュ数を無限大相当にしたときに、はじめて空間メッシュの離散化誤差がゼロの解を得ることが出来ます。


課題 3.1：ブロック各辺のメッシュ数を 2 倍、3 倍にした計算を行い、得られた実効増倍率と計算に要した時間を比較して下さい（計算時間の測定は、計算を実行するときに、time ./a.out というように

time を挿入してやればよいです）。メッシュ数の変更は main.cxx ファイル中の配列 xm、ym、zm の値を変更することで行えます。


さて、このベンチマーク問題には、厳密数値解としてモンテカルロ法と呼ばれる方法によって中性子輸
送方程式を解いて得られた値が与えられており、制御棒が引き抜かれた体系では 0.9780 とされています。

PLOS で得られた結果はどうだったでしょうか？かなりずれていたのではないかと思います。既に述べましたが、PLOS はあくまで中性子輸送方程式の近似式である拡散方程式の解を計算します。従って、その解は近似解に過ぎないと言えます。
それでは次に中性子輸送方程式を解くプログラムを使ってみましょう。main.takeda1.sn.cxx を実行して下さい（ここでは、5cm のブロックを 10×10×10 のメッシュに分割しています）。計算結果はどうなったでしょうか？
ここで用いたプログラムは、SNT という、中性子輸送方程式を解くためのものです。SNT を使った場合、
厳密数値解（すなわちモンテカルロ法による計算値 0.9780）にかなり近い値が得られた一方、計算時間が長くなってしまったのではないかと思います。数値計算の分野では、収束に至るまでの時間を可能な限り短縮するための手法の開発がひとつの大きなテーマとなっています。原子炉解析の分野でも同様であり、種々の方法が提案されています。ここでは例として、拡散合成加速法（Diﬀusion synthetic accelecation）を使ってみましょう。これは、輸送方程式の反復収束解の推定に、近似式である拡散方程式の解を用いる、という方法です。main.cxx（もともとは main.takeda1.sn.cxx）に sys.NoDSAAcceleration() なる行が（下から 5 行目くらいに）あるでしょうから、これを削除、もしくはコメントアウト（この行の頭に「//」を追加）して下さい。これをコンパイルして実行すると、だいぶ短い時間で収束解が得られることが分かるでしょう。
拡散合成法に加えて粗メッシュ有限差分加速法（Coarse-mesh ﬁnite diﬀerence acceleration）という方法
を用いると、さらに計算を加速することが出来ます。この方法は、最終行近くにある CalIgen(‘‘none’’);を CalIgen(‘‘cmfd’’); に変更する（‘‘none’’ を ‘‘cmfd’’ に変更する）ことで用いることが出来ます。
計算時間がさらに短縮することを確認して下さい。


課題 3.2：今までは制御棒引き抜き状態について計算していましたが、次は制御棒挿入状態について計算してみましょう。これは物質配置のデータを変更することで容易に実施可能です。拡散計算、輸送計
算でそれぞれ制御棒挿入前後の実効増倍率から制御棒反応度を求めて下さい。なお、制御棒挿入前後の実効増倍率を k1 、k2 とした場合、反応度は (1/k1 − 1/k2 ) = (k2 − k1 )/(k1 k2 ) と計算されます（単位は [∆k/kk ′ ]）。ちなみにモンテカルロ法では制御棒反応度は-0.0166∆k/kk ′ と得られていますので、この参照値と比較して下さい。


ここまでは実効増倍率の観察を行ってきましたが、次は中性子束の空間分布を見てみましょう。

X 軸方向に沿った中性子束分布を出力するメソッド ShowNeutronFluxAlongXAxis がありますので、それを使います。main.takeda1.dif 1.cxx のファイルの最後の方にこのメソッドがコメントアウトしてあるはずなので、実行させるときはコメントアウトを外して下さい。なお、このメソッドの引数ですが、一つ目、二つ目がそれぞれ中性子束分布を出力させる Y、Z 方向の位置を示します。これがふたつともゼロだと、Y、Z 方向それぞれについて一つ目のメッシュ位置の X 軸方向の中性子束分布を出力することになります。また、三つ目の引数は出力させる中性子束分布のエネルギー群を示します（なお、この数値を-1 にす3. Takeda ベンチマークの計算

7

ると全てのエネルギー群が出力されます）。ここで、メッシュ位置、エネルギー群ともに 0 から値が始まることに注意して下さい。


課題 3.3：制御棒引き抜き状態について、制御棒挿入位置を横断するように、X 軸に沿った中性子束分
布を 1 群、2 群ともに出力させて、Gnuplot や Excel 等で図示して下さい。また、制御棒を挿入した状
態でも同様の図を作成し、制御棒挿入前後で中性子束の空間分布がどのように変化するか、観察して
下さい。




なお、中性子束の単位は [/cm2 /s] ですが、PLOS 等の計算では全核分裂中性子源が 1.0 となるように中性子束を規格化しているため、中性子束の絶対値自体は意味を持ちません。また、1 群の中性子束は「高速中性子」、2 群の中性子束は「熱中性子」にそれぞれ対応します。
計算結果は画面に出力されるため、これを編集するのは少々面倒臭いです。そこで、./a.out > outputのようにしてやると、画面への出力がファイル output に書き込まれていくので、あとは output を編集す
ればよいことになります。このような方法をとるのがよいでしょう。

4. コマンド cbz の活用

4

8

コマンド cbz の活用
ここで述べる内容は Mac には対応していないことが分かっていますので、Mac ユーザーは読み飛ばして下さい。
これまでに解説したように、CBZ で計算を行う際には、ユーザが準備したファイルを main.cxx にコピーし、make を行ったのち、./a.out によって実行するという流れとなります。このような方式は、慣れてくると面倒臭く感じるでしょう。また、入力ファイルを何度も修正、変更するような場合には、オリジナルの入力ファイル mainXXXXXX.cxx を変更しているのか、main.cxx を変更しているのか、作業中に混乱してしまうような場合もあり、作業の効率性を著しく損ねてしまう可能性があります。
このような問題を解決するため、2020 年 9 月に修士を卒業した Qin Tian 君が、当研究室に在籍中にコ
マンド cbz を開発しました。このコマンドは、例えば計算を行いたい入力ファイルとして mainXXXXX.cxxがあった場合に、

cbz mainXXXXX.cxxと入力することにより、上記の一連の動作を自動的に行ってくれるという優れものです。
また、画面に表示される出力を output というファイルに出力したい場合は、

cbz mainXXXXX.cxx outputと入力すればよいです。
このコマンドを使うためにはいくつかの簡単な設定を行う必要がありますが、それについては CBZ/CBG
ディレクトリの ReadMe.txt ファイルを参照して下さい。それほど難しくはありません。
このコマンドは活用してみると、かなり有効であることに気付くと思いますので、面倒臭がらずにこの設定を行うことを強く推奨します。

5. 京都大学臨界集合体 KUCA の軽水減速架台の計算

59京都大学臨界集合体 KUCA の軽水減速架台の計算
核分裂エネルギーから電気を作る原子炉を実機原子炉、出力炉と呼んだりしますが、一方で研究用の臨界集合体というものがあります。これは原子炉内の中性子の振る舞いを実験的に検証するための装置であり、核分裂反応は持続しますが（臨界）、熱出力は極めて低く抑えられます（つまり中性子束レベルが極めて低いということ）。臨界集合体は日本にもいくつかあり、例えば、京都大学の臨界集合体 KUCA（Kyoto

University Critical Assembly）が挙げられます。
また、原子炉は中性子のエネルギー分布から分類することも可能であり、大部分が高速中性子で構成される原子炉を高速中性子炉と呼び、中性子の減速材（例えば水やポリエチレン、炭素）を用いることにより熱中性子の寄与を支配的にする原子炉を熱中性子炉と呼びます。
ここでは、京都大学臨界集合体 KUCA の軽水減速架台の計算を行います。軽水減速架台とは、軽水で満たされたプール中に、アルミニウムで被覆された薄い金属ウラン板が複数配置されるような体系です。軽水がふんだんに使われますから、熱中性子炉となります。
今回計算する対象である KUCA の体系図を Fig. 3 に示します。Z 軸方向の長さは固定されています。今回の計算では 57 cm として下さい。また、Y 軸方向の長さですが、今回の計算では 30 cm として下さい。
燃料は X 軸方向に追加されていくので、X 軸方向の長さが可変となります。

Fig. 3: KUCA の体系図
燃料領域では薄い燃料板が何枚も並べられますが、板と板との間には一定の間隙があり、そこには軽水が満たされることになります。ある燃料板の中心位置と隣の燃料板の中心位置の間の距離を燃料配列ピッチと呼びます。燃料板の厚さは一定ですので、燃料配列ピッチの値が大きくなると、間隙に入り込む軽水の割合が増加することになります。KUCA では、このピッチの長さがおよそ 3.0 mm の燃料を「C30」、3.5 mmのものを「C35」、4.5 mm のものを「C45」と、それぞれ呼称します。
燃料領域は、燃料板と間隙内の軽水、さらにはその他の構造物で構成されます。これらの構造を正確に模擬して数値計算するのは手間を要するので、ここでは「等価な均質の媒質」に置換したデータを用います。その値を Table 1 に示します。表から明らかなように、エネルギー群は 2 群として扱うことになります。また、この表には燃料領域を取り囲む軽水領域のデータも与えられています。


課題 5.1：C30、C35、C45 の 2 群断面積データを観察して、各断面積の燃料間の違いと、その違いが生じる理由について説明して下さい。

それでは実際に計算を始めましょう。CBGCAL/KUCA ディレクトリに移動して下さい。そこに main.keff 1.cxx
なるファイルがあると思いますので、それをエディタで開いて下さい。
この入力には、C30、C35、C45、水領域の断面積データが既に与えられています11 。入力下部に体系を
11 なお、燃料領域の 2 群核分裂生成断面積（νΣf,2 ）に、高速核分裂因子 ϵ を乗ずることによって、高速中性子による核分裂反応

5. 京都大学臨界集合体 KUCA の軽水減速架台の計算

10Table 1: KUCA の 2 群断面積データ構築するところがありますので、皆さんは適宜そこを変更して計算して下さい。なお、計算メッシュのサイズに留意する必要があるので、体系のサイズを変更する場合は、計算メッシュも適切に設定して下さい（１メッシュあたり 1cm 等、目安を決めておくことが必要です）。


課題 5.2：燃料領域を取り囲む水領域を無視した場合、燃料をどの程度装荷すれば（どの程度の X[cm]とすれば）臨界になるか、燃料として C30、C35、C45 を用いた場合それぞれについて計算して下さい。なお、計算メッシュは 1cm あたり 1 メッシュとすればよいでしょう。


以降の問題では水反射体を燃料の外側に設置しますので、main.keff 2.cxx を利用して下さい。


課題 5.3：燃料領域を厚さ 30cm の水領域が取り囲むとした場合、燃料をどの程度装荷すれば（どの程度の X[cm] とすれば）臨界になるか、燃料として C30、C35、C45 を用いた場合それぞれについて計算して下さい。


課題 5.4：中性子実効増倍率の観点から燃料領域を取り囲む水領域の厚さが無限相当と見做せる最小値
を、燃料として C30、C35、C45 を用いた場合それぞれについて計算して下さい。






課題 5.5：課題 5.3 で計算した臨界体系において、炉中心位置（YZ 平面でのほぼ中心位置）における X
軸方向の 1 群、2 群の中性子束分布を計算し、図示して下さい（C30、C35、C45 全てについて）。中性
子束分布は前節で説明した ShowNeutronFluxAlongXAxis メソッドで取り出して下さい。また、C30、

C35、C45 間の結果の違いについて説明して下さい。
を考慮しています。

6. 高速中性子炉の炉心燃料の中性子束エネルギースペクトルの計算

611

高速中性子炉の炉心燃料の中性子束エネルギースペクトルの計算
これまでの Takeda ベンチマークの計算や KUCA の計算では、中性子束のエネルギーは１群もしくは２群で取り扱っていました。これらの計算は一連の炉物理計算の手続きの中で「下流」にあたるわけですが、「上流」に近いところでは中性子束のエネルギー依存性をもう少し詳細に取り扱っています。その例として、
ここでは高速中性子炉の炉心燃料に関する中性子束のエネルギースペクトル12 計算を行いたいと思います。
ここで対象とするのは、中型高速中性子炉の炉心燃料集合体について、燃料ペレットや被覆管、冷却材、ラッパ管といった非均質構造を均質（ごちゃごちゃに混ぜあわせたもの）としてモデル化したものとなります。さらに、このような媒質が無限に広がると仮定するので、中性子束は空間的に一様となり、所謂「一点炉」となります。この場合、中性子束の空間的な振る舞いを考慮する必要がないため、中性子束のエネル
ギー的な振る舞いのみ、すなわち減速過程のみを考慮することになります。減速といっても高速中性子炉を対象としていますので、軽水のような減速材は媒質には含まれておらず、減速の程度は非常に小さいものとなります。
それではパッケージの monju ディレクトリ中の main.innercore.cxx を走らせてみて下さい。このプログラムは、対象媒質に対して中性子の減速計算を行い、中性子無限増倍率 k∞ と、エネルギー群毎の中性子束、中性子吸収断面積 Σa 、核分裂による中性子生成断面積 νΣf 、(n,2n) 反応断面積 Σ(n,2n) を出力します。
以下にその例を示します。

Listing 3: 高速中性子炉の炉心燃料の減速計算の出力例
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

# K= i n f: 1.40722# E top
1 . 0 0 0 0 0 e+07
7 . 7 8 8 0 1 e+06
6 . 0 6 5 3 1 e+06
4 . 7 2 3 6 7 e+06

Flux
1 . 2 8 6 5 0 e =01
3 . 2 8 0 0 3 e =01
7 . 2 1 1 6 7 e =01
1 . 1 8 8 5 0 e+00

Sigma a
Nu=S i g m a f
1 . 4 7 7 4 7 e =02 3 . 3 7 5 9 6 e =02
1 . 1 4 0 2 2 e =02 2 . 6 6 4 0 0 e =02
7 . 6 9 6 2 0 e =03 1 . 8 1 6 2 8 e =02
7 . 1 9 2 7 1 e =03 1 . 7 5 1 7 9 e =02(n ,2 n)7 . 0 6 9 5 3 e =03
1 . 7 9 6 6 6 e =03
1 . 0 2 1 8 5 e =05
0 . 0 0 0 0 0 e+00

( omitted )8 . 7 6 4 2 5 e =01
6 . 8 2 5 6 0 e =01
5 . 3 1 5 7 9 e =01
4 . 1 3 9 9 4 e =01
3 . 2 2 4 1 9 e =01

5 . 1 9 1 3 6 e =10
5 . 1 5 2 4 4 e =10
2 . 4 0 0 7 2 e =10
4 . 6 9 7 2 0 e =11
4 . 0 0 8 7 8 e =10

3 . 7 1 9 2 4 e =01
2 . 8 5 2 4 4 e =01
5 . 0 1 5 8 4 e =01
1 . 7 8 3 8 7 e+00
1 . 5 2 9 6 3 e+00

1 . 7 4 1 7 6 e =01
3 . 0 8 2 4 9 e =01
7 . 5 1 4 1 7 e =01
2 . 9 7 9 8 2 e+00
2 . 6 8 2 1 1 e+00

0 . 0 0 0 0 0 e+00
0 . 0 0 0 0 0 e+00
0 . 0 0 0 0 0 e+00
0 . 0 0 0 0 0 e+00
0 . 0 0 0 0 0 e+00

以下の課題は、main.innercore.cxx を使って解いて下さい。


課題 6.1：高速中性子炉の炉心燃料について、中性子束エネルギースペクトルを図示して下さい。なお、
横軸は対数軸とし、縦軸は線形軸、対数軸の２通りとすること。


中性子束に加えて、反応率 Σϕ を観察することで、どのようなエネルギー領域で反応が起こっているかを知ることが出来ます。


課題 6.2：高速中性子炉の炉心燃料について、吸収反応率 Σa ϕ、核分裂生成反応率 νΣf ϕ、(n,2n) 反応率 Σ(n,2n) ϕ のエネルギースペクトルを図示して下さい。なお、横軸は対数軸とし、縦軸は線形軸、対数軸の２通りとすること。
エネルギー群依存で定義される反応率 Σg ϕg を全ての g について和をとることで、エネルギー積分された反応率を得ることが出来ます。エネルギー積分された吸収反応率を Ra 、生成反応率を Rf 、(n,2n) 反応率をR(n,2n) と記述した場合、中性子無限増倍率は Rf /(Ra − R(n,2n) ) と計算することが出来ます（(n,2n) 反応
は中性子が 1 個増加する反応ですので、負の吸収反応と見做されます）。

12 中性子束のエネルギースペクトルとは、中性子束のエネルギー分布のことです。略して、「中性子スペクトル」と言う場合もあります。中性子束のエネルギー分布が高い側に偏っている場合は「硬いスペクトル」、低い側に偏っている場合は「軟らかいスペクトル」と言います。

6. 高速中性子炉の炉心燃料の中性子束エネルギースペクトルの計算

12





課題 6.3：エネルギー積分された反応率を用いて無限中性子増倍率を計算し、プログラムの出力に表示されているものと一致することを確認して下さい。

なお、課題 6.2、6.3 については、プログラムに自分で変更を加えて実施するのもよいでしょう。

7. 超小型臨界集合体 Godiva、Flattop-25 の計算

13超小型臨界集合体 Godiva、Flattop-25 の計算

7

ここでは、米国ロスアラモス国立研究所において、核兵器開発のために構築された超小型の臨界集合体

Godiva、及び Flattop-25 の計算を行います。Godiva は金属ウラン（U-235 の濃縮度が 90%以上）で構成される裸の球、Flattop-25 は金属ウラン球の外側に天然ウランを巻いた球体系であり、いずれも高速中性子炉に分類されます。臨界集合体 Flattop の体系図を Fig. 4 に示します。中心に核分裂性物質の「core」が配置され、その周囲を、漏れ出た中性子を反射して戻すための天然ウランが取り囲んでいます。

Fig. 4: Flattop の体系図
パッケージでは、CBGCAL/lanl.tutorial ディレクトリに、関連するデータが置かれています。

7.1

媒質の多群断面積の計算

単一もしくは複数の媒質から構成される原子炉の多群中性子輸送（もしくは拡散）方程式を解く際には、各媒質の多群断面積を予め計算しておく必要があります。前節の Takeda ベンチマークの問題では多群（２群）断面積が予め与えられていましたが、ここでは CBZ により多群断面積を計算する方法について説明します。なお、高速中性子炉と熱中性子炉とでは計算の手続きが僅かに異なることに注意が必要です。

lanl.tutorial ディレクトリにおいて、多群断面積の計算を行うファイルは main.xs.cxx なので、適宜それを横目で眺めながら本稿を読み進めて下さい。

7.1.1

多群断面積ライブラリ CBZLIB の読み込み

媒質の多群断面積を計算するためには、個々の原子核の、中性子に対する反応断面積の情報（多群断面積セット）が必要となります。CBZ では CBZLIB と呼ばれる多群断面積セットが用意されており、それを事前に呼び出しておく必要があります。CBZLIB は、基になっている評価済み核データファイル、エネル
ギー群数に応じて様々存在します。高速中性子炉を計算する際の CBZLIB の読み込み例を以下に示します。

Listing 4: CBZLIB の読み込み（高速中性子炉計算時）
1
2
3
4
5
6
7
8

// x s l i b r a r ys t r i n g l i b d i r ( ” . . / . . / CBGLIB/ j 4 . 7 0 g . i w t 7 / ” ) ;XSLibrary x s l i b ( l i b d i r , ”N=ENERGY” ) ;s t r i n g filename []={
”U233” , ”U234” , ”U235” , ”U238” , ” Pu239 ” ,” Pu240 ” , ” Pu241 ” , ” Pu242 ” , ” Ga000” ,};x s l i b . ReadFile (9 , l i b d i r , filename ) ;7. 超小型臨界集合体 Godiva、Flattop-25 の計算

14

2 行目の string 型の変数 libdir には、用いる CBZLIB のデータが置かれているディレクトリ名を指定します。CBZLIB は核種毎に異なったテキストファイルで同一のディレクトリにデータが置かれています。
この例ではディレクトリ名として j4.70g.iwt7 が指定されていますが、ここで使う CBZLIB は評価済み核データファイル JENDL-4.0 に基づく（j4 に対応）70 群の（70g に対応）ライブラリであることになります（iwt7 は多群断面積の作成に用いた重み関数を示しますが、ここでは分からなくてもよいです）。

3 行目では CBZLIB の実体である XSLibrary クラスのインスタンス xslib を作成しています。一つ目の引数ではライブラリデータが置かれているディレクトリの位置を指定しています。二つ目の引数ではエネル
ギー群のエネルギー境界の情報が格納されているテキストファイルの名前を指定しており、この部分は変更する必要はありません。

4 行目から 6 行目では、読み込む核種の指定を行っています。中性子との反応断面積が評価されている原子核の種類数は 400 を越えますが、いちいちこれら全ての原子核のデータを読み込むのは無駄なので、計算に必要な（計算対象の媒質に含まれる）核種のデータのみをここで指定します。
最終行では、XSLibrary クラスに実装されている ReadFile メソッドで核種毎のテキストファイルを読み込みます。引数は、一つ目が読み込む核種数に、二つ目が CBZLIB のディレクトリに、三つ目が読み込
むファイルの名前の配列に、それぞれ対応しています。
高速中性子炉を計算する場合には、一般的に 70 群、もしくは 175 群のライブラリを用いますが、熱中性子炉を計算する場合には 107 群のライブラリを用います。参考までに、熱中性子炉を計算する際の CBZLIBの読み込み例を以下に示します。

Listing 5: CBZLIB の読み込み（熱中性子炉計算時）
1
2
3
4
5
6
7
8
9
10
11
12
13
14

// +++ x s l i b r a r ys t r i n g l i b d i r ( ” /home/ c h i b a /CBGLIB/ j 4 . 1 0 7 g . i w t 7 / ” ) ;XSLibrary x s l i b ( l i b d i r , ”N=ENERGY” ) ;s t r i n g filename []={
”U234” , ”U235” , ”U236” , ”U238” , ” F019 ” ,”H001” , ”O016” , ” Al027 ” , ” S i 0 2 8 ” , ” S i 0 2 9 ” ,” S i 0 3 0 ” , ” Cu063 ” , ”Cu065 ” , ” Zn000 ” , ”Mn055”};
x s l i b . ReadFile (15 , l i b d i r , filename ) ;// +++ t h e r m a l s c a t t e r i n g d a t a reads t r i n g l i b t h ( ” /home/ c h i b a /CBGLIB/ j 4 . 1 0 7 g . i w t 7 / Thermal / ” ) ;x s l i b . GetLibData ( 1 0 0 1 0 ) . GetThScat ( ) . R e a d F i l e ( l i b t h , ”H. H2O” ) ;x s l i b . GetLibData ( 8 0 1 6 0 ) . GetThScat ( ) . R e a d F i l e ( l i b t h , ”O016” ) ;ここで注意すべきなのは、高速中性子炉用の CBZLIB を読み込む場合と異なり、熱中性子散乱データ（Ther-

mal scattering data）も読み込む必要があるという点です。それがこの例での 12 行目から 14 行目に対応しており、ここでは水分子中の水素原子と、（自由気体と仮定した）酸素原子の熱中性子散乱データを外部ファイルから読み込んでいます。このデータを読み込んで CBZLIB に含めることにより、熱中性子エネル
ギー領域における散乱での分子束縛効果を考慮することが可能となります。媒質に水を含む熱中性子炉を解析する際には、このデータを含めることは必須となります。

7.1.2

媒質データの作成

ここでは媒質データの作成方法について述べます。媒質データとしては、基本的には、含まれる原子核の種類とその数密度、媒質の温度が挙げられます。CBZ では Medium クラスが媒質に対応し、これらのデータを保持します。Medium クラスのインスタンスの作成例を以下に示します。

Listing 6: 媒質データの作成
1
2
3
4
5
6
7
8

// *** Medium d a t a
Medium god med ; // g o d i v a
god med . PutImax ( group ) ;god med . PutPL ( p l , 1 ) ;i n t god mat [ ] = { 9 2 2 3 4 0 , 9 2 2 3 5 0 , 9 2 2 3 8 0 } ;r e a l g o d d e n [ ] = { 4 . 9 1 8 4 e = 4 , 4 . 4 9 9 4 e = 2 , 2 . 4 9 8 4 e = 3};
god med . P u t N u c l i d e ( 3 , god mat , g o d d e n ) ;god med . P u t T e m p e r a t u r e F o r A l l N u c l i d e ( 3 0 0 . ) ;7. 超小型臨界集合体 Godiva、Flattop-25 の計算

15

2 行目では god med という名前の Medium クラスのインスタンスを作成しています。
3 行目では Medium クラスの PutImax メソッドにより、媒質データのエネルギー群数を指定しています。
この例では int 型の変数 group が引数として与えられています（この例では事前に group に値が与えられています）。なお、ここで与えるエネルギー群数は、用いる CBZLIB のエネルギー群数と一致させる必要があります。

4 行目では PutPL メソッドにより、その媒質データに対して散乱の非等方性をどの程度考慮するかを指定します。一つ目の引数が非等方散乱断面積で考慮する最大の（ルジャンドル多項式に展開した際の）次数に対応します。この例では int 型変数 pl にて指定しています。二つ目の引数については説明を省略します
が、通常は 1 を指定して下さい。

7 行目では PutNuclide メソッドにより、媒質に含まれる原子核の種類と数密度の情報を与えています。
引数の一つ目が原子核の種類数、二つ目が原子核の ID 13 、三つ目が数密度に対応しています。この例では、含まれる原子核の ID は 5 行目で、数密度は 6 行目でそれぞれ与えられています。
最下行では PutTemperatureForAllNuclide メソッドにより、媒質に含まれる全ての原子核の温度を指定しています。温度（単位は K）は引数で与えます。
解析する体系によっては、空気の媒質データを必要とする場合があります。そのようなときのために、空気の媒質データの作成例を以下に示します。MAT 番号の 70140、80160 はそれぞれ N-14、O-16 に対応しています。

Listing 7: 空気の媒質データの作成
1
2
3
4
5
6
7
8

// ( a i r )Medium a i r ;a i r . PutImax ( group ) ;a i r . PutPL ( 1 ) ;i n t mat0 [ ] = { 7 0 1 4 0 , 8 0 1 6 0 } ;r e a l den0 [ ] = { 3 . 2 2 6 9 e = 5 , 8 . 6 5 6 9 e = 6};
a i r . P u t N u c l i d e ( 2 , mat0 , den0 ) ;a i r . PutTemperatureForAllNuclide ( 2 9 3 . ) ;7.1.3

媒質の実効微視的断面積と巨視的多群断面積の計算、およびファイル出力媒質に含まれる核種の断面積が中性子の入射エネルギーに対して強い依存性を持っている場合（共鳴構造を有している場合）、多群（すなわちエネルギー平均）の微視的断面積は、その核種が媒質にどの程度含まれているか（換言すると、その核種が他の核種にどの程度希釈されているか）によって、大きく値が変わります14 。そして、媒質に応じて計算される多群の微視的断面積を多群実効微視的断面積と呼ぶこともあります。
また、媒質に含まれる全ての核種の実効微視的断面積 σ が得られれば、数密度 N を用いて媒質の巨視的∑
断面積 Σ（=Ni σi 、i は核種を示す index) を計算することができます。この巨視的断面積が中性子輸
i

送（拡散）方程式に現れるパラメータとなります。実効微視的断面積と巨視的断面積の計算例を以下に示します。

Listing 8: 実効微視的断面積、巨視的多群断面積の計算とファイル出力
1
2
3
4// *** s e l f = s h i e l d i n g t r e a t m e n t
O n e P o i n t C a l c u l a t o r opc ;opc . C a l S e l f S h i e l d i n g I n f i n i t e S y s t e m ( god med , x s l i b ) ;opc . C a l I n c i d e n t E n e r g y D e p e n d e n t C h i ( god med , x s l i b ) ;13 CBZ における原子核に対する ID の与え方ですが、原子核の原子番号を Z 、質量数を A、励起レベルを L としたとき、その原子核の ID を Z × 10000 + A × 10 + L と定義します。基本的に大部分の原子核は基底状態にあるため L = 0 となりますが、一部の原子核で励起状態で準安定となるものがある（例えば Am-242 など）ため、その場合は L > 0 となります。また、天然同位体組成の原子核については、質量数を 0 として定義します。例えば天然同位体組成の鉄については、ID は「260000」として与えられます。
14 一般に、断面積が大きい値を示すエネルギーでは中性子束は歪みます（小さくなります）。多群断面積は中性子束荷重で平均化しますので、断面積が大きい値を示すエネルギー点での荷重が小さくなるほど、エネルギー平均の断面積は小さくなります。この効果
を共鳴自己遮蔽効果と呼びます。ただし、断面積が共鳴構造を有する核種が他の核種に強く希釈されている場合には、共鳴による中性子束の歪みは小さくなるため、エネルギー平均の断面積は大きくなります。

7. 超小型臨界集合体 Godiva、Flattop-25 の計算
5616

god med . W r i t e F i l e ( x s d i r , ” god ” , m i c x s w r i t e ) ;微視的多群断面積、巨視的断面積の計算には CBZ の OnePointCalculator クラスを用いるので、そのインスタンス opc を 2 行目で作成しています。

3 行目では OnePointCalculator クラスのメソッド CalSelfShieldingInfiniteSystem により、Medium
クラスのインスタンスに含まれる核種の微視的多群断面積を計算し、さらに媒質の巨視的断面積の計算を行っています。このメソッドでは、一つ目の引数に計算する Medium クラスのインスタンスを、二つ目の引数に計算に用いる XSLibrary クラスのインスタンスを、それぞれ指定します。

4 行目では OnePointCalculator クラスのメソッド CalIncidentEnergyDependentChi により、核分裂
スペクトルの計算を行っています15 。核分裂性物質を含む媒質に対してはこのメソッドを行っておいたほうが確実ですが、中性子エネルギーが比較的高い超小型原子炉を除けば、このメソッドは使わなくてもよいです。
ここまでで中性子輸送計算に必要な巨視的多群断面積が作成され、Medium クラスのインスタンス god medに情報が格納されました。引き続いて中性子輸送計算を行う場合には、このインスタンスをそのまま用いればよいわけですが、一般的に、引き続いて行う中性子輸送計算は異なるファイルで行うため、ここで作成した god med のデータを一度外部ファイルに保存させる必要があります。ファイルへの保存は最下行で示されているように、Medium クラスの WriteFile メソッドで行います。このメソッドの引数は、一つ目がデータファイルを書き出すディレクトリ、二つ目がファイル名、三つ目が微視的断面積もファイルに出力するかを指定する boolean 変数となっています。
なお、熱中性子炉の媒質を計算する場合には、熱中性子散乱を考慮する必要があります。その場合の巨視的断面積の計算例を以下に示します。

Listing 9: 実効微視的断面積、巨視的多群断面積の計算とファイル出力（熱中性子炉の場合）
1
2
3
4
5
6

// s e l f = s h i e l d i n g t r e a t m e n t
O n e P o i n t C a l c u l a t o r opc ;opc . C a l S e l f S h i e l d i n g I n f i n i t e S y s t e m ( s o l 1 , x s l i b ) ;opc . C a l T h e r m a l S c a t t e r i n g M a t r i x ( s o l 1 , x s l i b , 3 . 9 3 ) ;s o l 1 . W r i t e F i l e ( ” . / XS/CBGXS/ ” , ”SOL1” ) ;この例では Medium クラスのインスタンス sol1 の巨視的断面積を計算しています。熱中性子散乱を考慮させるため、OnePointCalculator クラスの CalThermalScatteringMatrix メソッドを用いています。このメソッドの三つ目の引数は熱中性子散乱を考慮する上限のエネルギーに対応しており、この例のように3.93 eV を指定してやれば十分です。なお、この例では CalIncidentEnergyDependentChi メソッドを省略しています。
それでは main.xs.cxx を実行させて、Godiva、Flattop-25 を構成する媒質の断面積データを作成してみて下さい。XS/CBGXS ディレクトリに god、fu、fb という三つのファイルが生成されたことと思います。このうち god は Godiva の燃料領域、fu、fb はそれぞれ Flattop-25 の燃料、ブランケット領域の断面積デー
タに対応しています。

7.2

原子炉の計算

原子炉の各媒質の多群断面積が計算されたので、次は原子炉の中性子の振る舞いを支配する多群中性子輸送（拡散）方程式を数値的に解き、原子炉内の中性子の挙動を評価します。ここでは、超小型臨界集合体

Flattop-25 を、離散座標法に基づく中性子輸送ソルバー SNR で計算する場合を例として示します。
15 CBZLIB には入射エネルギー依存の核分裂スペクトルが与えられていますが、中性子輸送計算ソルバーでは入射エネルギーに依存しない核分裂スペクトルを用います。そこで、このメソッドでは、近似的に媒質の中性子エネルギースペクトルを計算し、それを用いて入射エネルギーに依存しない核分裂スペクトルの計算を行っています。

7. 超小型臨界集合体 Godiva、Flattop-25 の計算

17パッケージでは main.flatu.cxx がこの部分に該当していますので、以降はそのファイルを横目で眺めながら読み進んでいって下さい（Godiva の場合は main.godiva.cxx が対応します）。

7.2.1

体系情報の作成

球形の原子炉の体系情報は CartMeshInfo クラスにより定義します。Flattop-25 の体系情報の作成例を以下に示します。

Listing 10: 球形の原子炉の体系情報の作成
1
2
3
4
5
6
7

// CartMeshInfo
C a r t M e s h I n f o cmi ;i n t xm[ ] = { 3 0 , 1 6 , 3 0 } ;// f o r f l a t t o p =u
r e a l x l [ ] = { 6 . 1 1 5 6 , 3 . , 1 5 . 0 0 8 6 } ; // f o r f l a t t o p =u
i n t mat [ ] = { 0 , 1 , 1 } ;cmi . PutMeshInfo ( 3 ,xm, x l , mat ) ;cmi . PutBoundaryCondition ( ” R e f l e c t i v e ” , ”Vacuum” ) ;2 行目では CartMeshInfo クラスのインスタンス cmi を作成しています。
6 行目では cmi に対して、領域の分割、物質配置、メッシュの分割情報を PutMeshInfo メソッドにより与えています。引数はそれぞれ、一つ目が領域分割数を、二つ目が各領域のメッシュ分割数を、三つ目が各領域の幅（単位は cm）を、四つ目が各領域の物質番号を示しています。この例では、領域のメッシュ分割
数、幅、物質番号が 3 行目から 5 行目に与えられています。球体系を扱っているため、一つ目の領域の幅
は球中心位置からの幅となります。すなわち、この例では、半径 24.1cm の球の内側に、半径 6.1156cm の球があると考えればよいわけです。なお、各領域の幅を領域の外側位置（球の中心を 0 としたときの各領域の右端位置）で指定することも可能で、その場合は PutMeshInfo メソッドの五つ目の引数として文字列‘‘cumulative’’ を追加すればよいです。
最終行では境界条件を指定しており、球形状の原子炉の場合はこの例の通り、左側に Reflective（完全反射）、右端に Vacuum（真空）を指定します。

7.2.2

中性子輸送計算のための角度求積セットの指定

離散座標法に基づく中性子輸送ソルバーでは、いくつかの特定の方向を向いた中性子束についてのみ輸送計算を行い、（スカラー）中性子束などの角度積分値は求積法で求めます。そのため、この求積法に用いる角度離散点とその重みからなる求積セットを計算の前に指定する必要があります。角度求積セットの指定の例を以下に示します。

Listing 11: 角度求積セットの指定
1
2
3// Quadrature
SNRQuadrature quad1 ( p l ) ;quad1 . PutDoubleGaussian ( 2 4 ) ;2 行目では、SNRQuadrature クラスのインスタンス quad1 を作成しています。引数としては非等方散乱
断面積に対して考慮する最大次数を与えます。

3 行目では PutDoubleGaussian メソッドにより、二重ガウス求積セット（Double Gaussian angularquadrature set）を指定しています。引数として与えられているのはその離散点数となります（離散点数がn の場合、次数は n/2 − 1 となります）。
非等方散乱断面積に対する次数、求積セットに与える離散点数はそれぞれ十分大きな値とする必要があります。すなわち、これらをパラメータとした計算を実施し、これらの値を増加しても計算結果に有意に影響を与えない程度（例えば実効増倍率では 0.0001 以内程度）まで大きな値を設定する必要があります。

7. 超小型臨界集合体 Godiva、Flattop-25 の計算
7.2.3

18

媒質情報の読み込み

ここでは計算対象としている原子炉に含まれる媒質の情報を外部ファイルから読み込みます。すでに体系情報の作成の箇所で説明したように、この例で対象としている原子炉 Flattop-25 は二種類の媒質（濃縮ウ
ランの燃料と、天然ウランのブランケット）から構成されています。媒質情報のファイルからの読み込みの例を以下に示します。

Listing 12: 媒質情報の読み込み
1
2
3
4
5

// Medium
Medium med [ mednum ] ;s t r i n g mdir ( ” . / XS/CBGXS/ ” ) ;med [ 0 ] . R e a d F i l e ( mdir , ” f u ” , 5 ) ;med [ 1 ] . R e a d F i l e ( mdir , ” f b ” , 5 ) ;

// o
// o

まず、Medium クラスのインスタンスの配列 med を作成します。この例では int 型の変数 mednum が用いられていますが、この値は 2 となっています。Medium クラスにはファイルからデータを読み込むメソッド
として ReadFile が実装されており、引数としてファイルの位置、ファイル名を指定します。

7.2.4

計算オプションの設定

中性子輸送計算には反復計算の収束条件、最大反復回数などの計算オプションが必要となります。その情報を持つのが GeneralOption クラスであり、通常の輸送方程式（前進式 forward equation）を解くか、随伴式 adjoint equation を解くかもこのクラスで定義します。GeneralOption クラスの使用例を以下に示します。

Listing 13: 計算オプションの設定
1
2// GeneralOption
GeneralOption opt ;この例では GeneralOption クラスのインスタンス opt を作成しています。インスタンス作成時には種々の計算条件はデフォルト値に設定されます。解くべき方程式については、デフォルトでは前進式に設定されます。随伴式を解く場合には、PutAdjointCal メソッドを用います。

7.2.5

中性子輸送方程式の数値計算

以上により中性子輸送計算に必要な情報は全て準備されました。ここでは中性子輸送計算について説明します。
中性子輸送方程式を解くためのプログラムは CBZ にいくつか備えられていますが、球体系を離散座標法
を用いて計算する場合は SNR を用います。SNR を用いた輸送方程式の計算例を以下に示します。

Listing 14: SNR による計算
1
2
3
4
5
6
7
8
9
10

// SNR f o r a d j o i n t
SNRSystem s n r ( dim , group , mednum ) ; // group & mediums n r . AddMedium ( med [ 0 ] ) ;s n r . AddMedium ( med [ 1 ] ) ;s n r . PutCartMeshInfo ( cmi , ” S p h e r e ” ) ;s n r . PutPL ( p l ) ;s n r . S e t Q u a d r a t u r e (&quad1 ) ;snr . PutGeneralOption ( opt ) ;s n r . PutWriteFlux ( ) ;r e a l k1=s n r . C a l I g e n ( ) ;2 行目では SNR の実体となるクラス SNRSystem のインスタンス snr を生成しています。引数は一つ目が次元数（球体系の場合は 1）、二つ目がエネルギー群数、三つ目が原子炉に含まれる媒質の数に対応しています。

7. 超小型臨界集合体 Godiva、Flattop-25 の計算

19

3、4 行目では媒質を示す Medium クラスのインスタンスの情報を追加していきます。追加された順番で
SNR は Medium クラスのインスタンスに物質 ID を付与していくため、この例では SNR は med[0] の物質
ID を 0、med[1] の物質 ID を 1 と認識しています。
5 行目では体系情報を示す CartMeshInfo クラスのインスタンスの情報を SNR のインスタンスに与えています。引数の二つ目は体系の形状を定義する文字列であり、球体系を扱う場合はこの例の通りにして下さい。

6 行目では非等方散乱断面積で考慮する最大次数を定義します。
7 行目では作成した角度求積セットを定義する SNRQuadrature クラスのインスタンスを、8 行目ではGeneralOption クラスのインスタンスを snr に与えます。
9 行目は角度中性子束の情報を SNR に持たせておくという指定であり、通常は不要です（感度計算、摂動計算を行う際には角度中性子束の情報が必要となるため、このメソッドを行う必要があります）。
最後の 10 行目で snr に中性子輸送計算を行わせ、得られた実効増倍率を real 型の変数 k1 に代入して
います。




課題 7.1：それでは実際に main.godiva.cxx、main.flatu.cxx を実行させてみて下さい。実効増倍率
としてどの程度の値が得られましたか？この体系は臨界ですので、実効増倍率の実験値はほぼ 1.0 となります。計算値がどの程度実験値と一致しているか確認してみて下さい。また、計算値と実験値のず
れは、燃料領域の厚さにしてどの程度か評価してみて下さい。


7.2.6



炉中心反応率比の計算

以上では、超小型臨界集合体 Godiva、Flattop-25 の中性子実効増倍率を計算する方法について説明してきましたが、炉物理パラメータは中性子増倍率以外にもいろいろあります。ここでは臨界集合体を用いた臨界実験において測定される「反応率比（スペクトルインデックス）」の計算方法について説明します。
反応率比の計算には、反応断面積に加えて着目位置での中性子束のエネルギースペクトルの情報が必要となります。従って、中性子輸送計算を行ったあとに初めて計算可能となります。
体系の任意のメッシュにおける微視的反応率比を計算するメソッドとして、GeneralSystem クラスにはCalPointReactionRateRatio が実装されています。その使用例を以下に示します（サンプルの main.godiva.cxxの下側にコメントアウトされています）。

Listing 15: CalPointReactionRateRatio メソッドの使用例
12
34// ( Re action r a t e r a t i o c a l c u l a t i o n )i n t r r r p o s =0;
r e a l tmp1=s n r . C a l P o i n t R e a c t i o n R a t e R a t i o ( r r r p o s , 9 2 2 3 8 0 , s i g f , 9 2 2 3 5 0 , s i g f ) ;cout<<”#\n# R e a c t i o n r a t e r a t i o o f U= 238(n , f ) t o U= 235(n , f ) : ”<<tmp1<<” \n” ;引数ですが、一つ目が反応率比を計算するメッシュ位置となります。この例では 0 が指定されていますので、最も中心寄りのメッシュで計算されることになります。二つ目と三つ目の引数は反応率比の分子の反応率を定義するためのもので、それぞれ核種 ID と反応名（核分裂反応だと sigf、捕獲反応だと sigc を指定）となっています。この例では U-238 の核分裂反応となっています。四つ目と五つ目の引数は分母の反応率を定義するためのものです。
注意しなければならないのが、微視的反応率比を計算する場合は、各媒質データ（Medium クラスのインスタンス）が核種毎の微視的断面積データを持っている必要がある点です。多群断面積を作成するファイル

main.xs.cxx では、微視的断面積の情報はファイルに書き出さないようになっていますので、これを指定する boolean 変数 micxs write を true にして、多群断面積を再作成する必要があります。

8. 臨界実験装置 ZEBRA に構築された炉心 MZA の計算

820臨界実験装置 ZEBRA に構築された炉心 MZA の計算
ここではナトリウム冷却高速中性子炉の特性を把握するために英国の臨界実験装置 ZEBRA で構築された炉心 MZA の解析を行います。パッケージでは、CBGCAL/mza.tutorial ディレクトリに、関連するデー
タが置かれています。
ここで、繰り返しになるかもしれませんが、簡単に原子炉、臨界集合体の核特性解析方法について説明します。
実機原子炉では、複数本の燃料ピンが規則的に配置された燃料集合体が炉心内に規則的に配置されます。
一方、臨界集合体では、さまざまな物質板が規則正しく並べられた入れ物が、炉心内に規則的に配置されます。臨界実験装置 ZEBRA では、鉛直方向に延びた入れ物に薄い物質プレートが積み上げられる形で配置されます。これは実機原子炉の燃料集合体に該当するもので、燃料要素と呼ばれます。燃料要素の概念図をFig. 5 に示します。物質板の一辺の長さは 2 インチ（5.08 cm）程度となっています。

Fig. 5: 臨界実験装置 ZEBRA における燃料要素
また、原子炉の中央部に配置される燃料要素内には、燃料領域を構成する部分、燃料領域を上下に挟む
ブランケット領域を構成する部分、さらにブランケット領域を挟む中性子反射体（もしくはガスプレナム）
領域を構成する部分があり、その各々は、規則正しく配置された種々の物質プレートで構成されます。この
燃料要素内の燃料領域部分、ブランケット部分などを単位格子と呼びます。Fig. 6 に MZA 炉心で燃料領域を構成する単位格子の物質配置例を示します。Pu プレート（図中「Pu METAL」）やナトリウムプレー
ト（同「Na」）で構成されており、実機高速炉の燃料集合体の組成を模擬しています。

8. 臨界実験装置 ZEBRA に構築された炉心 MZA の計算

21

Fig. 6: MZA 炉心の単位燃料格子の物質配置例
原子炉や臨界集合体の全炉心体系を計算する際に、燃料集合体、燃料要素内の詳細な構造を考慮することは非常に大変です。そこで、一般的には、単位格子の計算と全炉心の計算を分離する取り扱いをします。
前者は格子計算、後者は炉心計算と呼ばれます。
格子計算では、単位格子のみを取り扱い、それに周期境界条件を課すことにより、単位格子が無限に配列されるような系を考え、その系での中性子輸送方程式を解きます。その後、単位格子と等価となるような均質の媒質を作成し、それを炉心計算で用いるようにします。すなわち、炉心計算においては単位格子は均質な媒質として扱われ、内部の詳細な構造は全く考慮する必要が無くなります。この単位格子に対する操作は均質化と呼ばれ、炉心解析での重要な手続きとなっています。
なお、このチュートリアルで用意されている入力データ main.xs.cxx では、単位格子に含まれる種々の原子核の数密度はすでに体積重みで均質化（平均化）されたものとして与えられていますので、この入力を用いた場合、単位格子は単一の物質で構成されていると見做されます。

MZA の炉心体系図を Fig. 7 に示します。左は上から見た図で、アルファベット A から D で示されている箇所に燃料を含んだ燃料要素が装荷されています。その外側で何も記載されていないところにはブランケットを含んだ燃料要素が装荷されます。また、右は横から見た図（全炉心の 1/4）であり、軸方向中心に
燃料単位格子（C11-1A 等、「Core」と記載されている領域）が、その上下に軸方向ブランケット（B11-1
等）が配置されていることが分かります。

Fig. 7: MZA の炉心体系図（左：径方向、右：軸方向）

8. 臨界実験装置 ZEBRA に構築された炉心 MZA の計算

8.1

22

格子計算（多群断面積の計算）

MZA は 10 種を超える単位格子（単位媒質）から構成されますので、その各々について多群断面積を計算する必要があります。入力ファイル main.xs.cxx を用いて多群断面積の計算を行ってみましょう。XS/CBGXS
ディレクトリに多群断面積データが生成されることと思います16 。
ここで、main.xs.cxx を眺めてみて下さい。この入力ファイルでは 15 個の媒質の多群断面積を作成しますが、「C11」で始まるものは炉心燃料に、「B11」で始まるものはブランケットに、それぞれ対応します。
また、「C11-1AV」「B11-1V」のように末尾に「V」がついているものは、「V」がついていない媒質からナトリウムが除去されたことを意味します。つまり、「C11-1AV」は「C11-1A」から冷却材模擬物質であるナトリウムが除去されたものとなります。これらは後述する Na ボイド反応度を実験によって測定する際に使われます。なお、実験は室温で行われているため、ナトリウムは固体状のプレートとなっています。実際の実験では、ナトリウムがボイド化したことを模擬するためステンレスで覆われた中空のプレートをナトリ
ウムプレートの代わりに設置しています。

8.2

炉心計算

8.2.1

中性子実効増倍率の計算

はじめに、MZA の XYZ 体系の固有値（実効増倍率）を計算してみましょう。main.xyz.cxx を走らせてみて下さい。固有値計算が始まって、ちょっとしたら終わるはずです。

main.xyz.cxx の中身を眺めてみて下さい。はじめに燃料要素の情報を定義し（Assembly、AssemblySetのあたり）、次に、それが XY 平面上でどのように配置されるかを指定しています（CartCore のあたり）。
その後、断面積データを読み込んで、固有値計算を行います。
次に、いろいろ計算オプションを変えて計算してみましょう。main.cxx を開いて下さい。
ファイルの頭部分で string solver=‘PLOS’ という箇所があります。ここでは、使用するソルバー名を指定しますが、さきほどの計算は拡散ソルバー PLOS で行っています。では、輸送ソルバー SNT で計算してみましょう。PLOS を SNT と書き換えて実行してみて下さい。始めに、PLOSE というソルバーが起動しますが、これは初期値を求めるための端メッシュ差分拡散ソルバーです（PLOS は中心メッシュ差分）。な
お、PLOSE は別に、SNT 内部の拡散合成加速のソルバーとしても使われています。PLOSE による計算が終了して初期値が推定されたのち、SNT による計算が始まります。

PLOS 及び SNT で得た実効増倍率はどの程度になったでしょうか？この体系の実効増倍率の実験値としては 1.0100 が与えられています。大きく異なる結果が得られたと思いますが、これは格子計算（多群断面積の計算）において、単位格子の非均質な構造を無視し、体積重みで均質化された核種数密度データを用いていることに起因しています。実は、パッケージには燃料単位格子 C11-1A、C11-1B、C11-1C、C11-1D について、その非均質構造を一次元平板にモデル化した格子計算のための入力データ main.xs het.cxx が含
まれています。それではこの入力を走らせてみましょう。この入力データでは、Fig. 6 のような体系に対して、中性子輸送方程式を衝突確率法により解き、中性子束分布を求め、それを用いて均質断面積を計算させています。均質断面積を定義するファイルは XS/CBGXS ディレクトリに同一の名前で上書きされることにな
ります。このようにして格子計算で非均質構造を考慮して得た均質断面積を用いて、再度、main.xyz.cxx
により体系の実効増倍率を求めてみましょう。だいぶ実験値（1.0100）に近い値が得られるものと思います。

16 さきほどは Flattop-25 の計算を行っていましたが、ここでは MZA が計算対象となっていますので、新たに多群断面積の計算が必要になります。

8. 臨界実験装置 ZEBRA に構築された炉心 MZA の計算

23





課題 8.1：非均質構造を考慮した格子計算により均質断面積を求めることにより実効増倍率が増加する
結果が得られます。その要因は複数ありますが、そのひとつとして格子平均の中性子吸収断面積が低下することが挙げられます。格子平均の中性子吸収断面積 Σ̃a は、非均質体系における領域 i の中性子束分布を ϕi 、体積を Vi としたとき、以下の式で求められます。
I∑Σi ϕ i V iΣ̃a = i=1I∑(1)
ϕ i Vii=1

非均質構造を無視した格子計算（main.xs.cxx に対応）は、中性子束分布 ϕi を平坦と仮定していることとほぼ同一と考えてよいです。Σ̃a を計算する際、ϕi の分布を考慮することによって、考慮しない場合と比較して Σ̃a が小さくなる理由を述べて下さい。


8.2.2



摂動理論によるナトリウムボイド反応度の計算

それでは次に、高速炉の炉心設計で非常に重要なパラメータであるナトリウムボイド反応度の計算をしてみましょう。ナトリウムボイド反応度とは、冷却材であるナトリウムが沸騰してボイド化した場合に炉心に添加される反応度です。
反応度 ρ はふたつの状態の炉心の実効増倍率 k1 、k2 から、ρ = 1/k2 − 1/k1 として計算されますが、k1と k2 との差が小さい場合には反応度の有効数字は小さくなってしまいます。そこで、通常、反応度計算は摂動理論を用いて行います。摂動理論については、炉物理の教科書や、千葉のホームページの「炉物理プロ
グラム演習：随伴方程式と摂動計算」を読んで勉強して下さい。
それではナトリウムボイド反応度を計算してみましょう。なお、ここでは格子計算を均質モデルで揃える必要があるので、main.xs.cxx を再度流して、全ての格子の均質断面積を再計算して下さい。

main.xyz.vid.cxx を眺めて下さい。燃料要素として 13 種類が定義されていますが、最後の５つでは、
燃料領域、ブランケット領域（の一部）で冷却材がボイド化しています（「c11-1av」「b11-1v」という名前で区別されています）。摂動計算では、基準体系の随伴中性子束と、摂動体系（すなわち冷却材がボイド化した体系）の中性子束が必要になりますので、基準体系と摂動体系について中性子輸送（拡散）方程式を解きます。CBZ では収束条件等の計算条件は GeneralOption クラスで定義され、そのインスタンスがソルバーに渡されます。随伴系の計算をするためには、GeneralOption クラスで PutAdjointCal メソッドを使います。main.xyz.vid.cxx の中を確認してみてください。
また、摂動体系の体系情報は、基準体系の集合体配置を定義する CartCore クラスのインスタンスに対してChangeAssembly メソッドを実行させることにより体系を変更させて得ています。ChangeAssembly(0,0,8)とした場合は「X=0、Y=0 の集合体を No.8 集合体に置換せよ」というメソッドになります。この例では炉中心位置にある 9 つの燃料要素が、ナトリウムが一部ボイド化した燃料要素に置換されます。
ふたつの体系の固有値を計算したあと、CalReactivity メソッドにより反応度が計算されます。なお、随伴系では CMFD 加速は有効ではないことに注意して下さい。
それでは main.xyz.vid.cxx を実行してみましょう。固有値計算が二回行われ、最後に反応度が出力されたかと思います17 。
摂動計算によるナトリウムボイド反応度の計算例を以下に示します。

Listing 16: ナトリウムボイド反応度の計算例
12# +++ R e s u l t s +++#17 CBZ では計算結果は全て画面上に出力されるので、ファイルに保存したい場合には、（cbz コマンドを利用していないならば）./a.out > output のようにしてください。

8. 臨界実験装置 ZEBRA に構築された炉心 MZA の計算
3
4
5
6
7
8
9
10
11
12
13
14

24

# Yield: = 1.46472 e =05# A b s o r p t i o n : 1 . 6 0 0 0 3 e =04# S c a t t e r i n g : 3 . 4 4 8 6 8 e =04# N2N: = 1.68593 e =10# ( Non=Leak ) : 4 . 9 0 2 2 4 e =04# Leakage =x: = 1.91749 e =05# Leakage =y: = 1.91785 e =05# Leakage =z: = 4.31141 e =05# ( Leakage ): = 8.14674 e =05## ** P e r t u r b a t i o n Cal .
: 0.00040876# ** D i r e c t Cal .
: 0.00039815

摂動計算では、反応度は次のように要因別に、さらにはエネルギー群別に分類され出力されています。

 Yield 成分：摂動によって核分裂断面積が変化し、核分裂による中性子生成率が変化したことにより印加される反応度成分

 Absorption 成分：摂動によって吸収断面積が変化し、中性子消滅率が変化したことにより印加される
反応度成分

 Scattering 成分：摂動によって散乱断面積が変化し、中性子スペクトルが変化したことにより印加される反応度成分 (n,2n) 成分：摂動によって (n,2n) 断面積が変化し、(n,2n) 反応による中性子生成率が変化したことにより印加される反応度成分

 Leakage 成分：摂動によって全断面積が変化し、中性子漏洩率が変化したことにより印加される反応
度成分
例えば、ナトリウムがボイド化した場合、媒質の吸収断面積はナトリウムの分だけ小さくなるはずです。それによる反応度効果は上記の Absorption 成分に対応します。
これら各成分の総和が Perturbation Cal. として出力されます。なお、Leakage 成分を除いた成分の和は「Non-leakage 成分」としてまとめられます。また、Leakage 成分は方向毎にも与えられます。

Direct Cal. はふたつの実効増倍率からの反応度計算値を示します。両者は収束精度内でおおよそ一致していることを確認して下さい（なお、固有値計算の実効増倍率に対する収束条件は 10−5 となっています）。

課題 8.2：main.xyz.vid.cxx を流して得られるナトリウムボイド反応度の計算結果について、要因毎の反応度を横軸をエネルギーとして図示して下さい。


ナトリウムボイド反応度の Scattering 成分は高エネルギー領域では負の値をとる。このことは、η 値（νΣf /Σa ）のエネルギー依存性で説明することが出来る。


課題 8.3：エネルギー依存の η 値をプロットとし、その情報を用いて、課題 8.2 で得られたナトリウム
ボイド反応度の Scattering 成分が高エネルギー領域で負の値をとることを説明して下さい。なお、核分裂生成断面積 νΣf 、吸収断面積 Σa は main.xyz.cxx の下部にコメントアウトしてある部分を流すことによって得ることが出来ます。






課題 8.4：main.xyz.vid.cxx には、ナトリウムの一部がボイド化している燃料要素が５種類あります（Assembly クラスのインスタンスの asmi[8] から asmi[12]）。デフォルトでは炉中心位置の燃料要素のいくつかを asmi[8] に置換した場合の反応度を計算させていますが、これを asmi[9] から asmi[12]の各々に置換した場合の反応度を要因毎に計算し、結果を表にまとめて下さい。なお、それぞれの燃料要素ではどの領域のナトリウムがボイド化されているのか、しっかり把握して下さい（asmi がどのように定義されているかを見れば分かると思います）。


8. 臨界実験装置 ZEBRA に構築された炉心 MZA の計算
8.2.3

25

実効遅発中性子割合の計算

原子炉の動特性計算や、反応度の dk/kk ′ 単位からドル単位への換算に用いられる実効遅発中性子割合

βeﬀ も重要な炉物理パラメータです。入力ファイル main.xyz.beff.cxx はその計算を行うためのもので、詳細は省きますが、とりあえず流してみましょう。以下のような結果が得られる筈です。

Listing 17: 実効遅発中性子割合の計算例
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

################################################################ Delayed n e u t r o n p a r a m e t e r c a l c u l a t i o n## B e t a e f f = 3 . 4 8 5 7 2 5 e =03## [ N u c l i d e =w i s e ( r e l a t i v e c o n t r i b u t i o n ) ]## 922340 : 1 . 1 2 8 8 0 2 e =07( 3 . 2 3 8 3 5 7 e = 05)# 922350 : 1 . 7 0 5 1 2 0 e =04( 4 . 8 9 1 7 2 4 e = 02)# 922380 : 1 . 5 4 0 0 0 3 e =03( 4 . 4 1 8 0 2 8 e = 01)# 942390 : 1 . 4 5 5 3 9 1 e =03( 4 . 1 7 5 2 8 9 e = 01)# 942400 : 1 . 2 5 5 7 1 4 e =04( 3 . 6 0 2 4 4 7 e = 02)# 942410 : 1 . 8 9 9 3 1 5 e =04( 5 . 4 4 8 8 3 7 e = 02)# 942420 : 4 . 2 0 3 1 4 2 e =06( 1 . 2 0 5 8 1 6 e = 03)## [ N u c l i d e&f a m i l y =w i s e ]
#
#
1
2
3
4
5
6
# 922340 6 . 1 1 8 e =09 2 . 2 0 5 e =08 2 . 0 2 6 e =08 4 . 4 0 1 e =08 1 . 4 9 8 e =08 5 . 4 6 3 e =09# 922350 5 . 8 4 7 e =06 3 . 6 8 1 e =05 3 . 5 8 3 e =05 6 . 6 2 4 e =05 2 . 1 5 1 e =05 4 . 2 7 6 e =06# 922380 1 . 9 8 3 e =05 2 . 1 0 8 e =04 2 . 4 7 1 e =04 6 . 0 0 1 e =04 3 . 4 6 2 e =04 1 . 1 5 9 e =04# 942390 4 . 7 6 4 e =05 4 . 1 6 3 e =04 3 . 1 6 7 e =04 4 . 7 5 3 e =04 1 . 2 5 0 e =04 7 . 4 4 8 e =05# 942400 3 . 9 8 7 e =06 3 . 1 8 8 e =05 1 . 8 8 1 e =05 4 . 1 4 9 e =05 2 . 2 5 1 e =05 6 . 8 9 1 e =06# 942410 3 . 4 1 8 e =06 4 . 2 8 1 e =05 2 . 6 8 7 e =05 6 . 6 3 8 e =05 3 . 7 4 7 e =05 1 . 2 9 9 e =05# 942420 8 . 2 2 9 e =08 9 . 7 6 4 e =07 5 . 2 3 7 e =07 1 . 3 7 2 e =06 9 . 4 6 9 e =07 3 . 0 1 9 e =07#( sum ) 8 . 0 8 1 e =05 7 . 3 9 6 e =04 6 . 4 5 8 e =04 1 . 2 5 1 e =03 5 . 5 3 7 e =04 2 . 1 4 9 e =04## [ Averaged d ec ay c o n s t a n t ]#( beta =w e i g h t )( beta =w e i g h t = i n v e r s e )# 1
1 . 2 9 2 e =02
1 . 2 9 1 e =02# 2
3 . 0 7 2 e =02
3 . 0 7 0 e =02# 3
1 . 2 8 2 e =01
1 . 2 7 6 e =01# 4
3 . 3 8 4 e =01
3 . 3 7 0 e =01# 5
1 . 2 7 7 e+00
1 . 2 4 2 e+00
# 6
3 . 4 4 7 e+00
3 . 3 2 5 e+00###############################################################

まずはこの体系の βeﬀ として 0.00349 が得られています。これはすなわち全核分裂中性子のうち遅発中性子が占める割合がおよそ 0.35%であることを示しています。

βeﬀ は核分裂核種毎に定義される βeﬀ,i の総和として計算されるため、核分裂性核種毎の寄与も得ることができます。その結果が上記の例における Nuclide-wise (relative contribution) に対応しており、

βeﬀ,i が相対寄与（括弧内）とともに示されています。この体系では U-238 と Pu-239 からの遅発中性子が主であることが分かるでしょう。
さらに、βeﬀ は遅発中性子先行核家系（もしくは先行核群）毎にも定義されることから、核分裂性核種毎
でかつ先行核家系毎の数値もそのあとに示されています。
最後の Average decay constant は、先行核家系毎の崩壊定数の核種平均値 λi の結果を示しています（先行核の崩壊定数は核分裂性核種毎に異なる値が与えられるため、炉心全体としての値を得るためには平均化の操作が必要となります）。平均のとり方にいくつか異なる方法があり、ここでは２通りの方法に基づく結果が示されています。
これらの βi や λi といったパラメータが、原子炉の動特性解析に用いられることになります。

9. おわりに

9

26

おわりに
ということで、おおよその流れは掴めたものと思います。お疲れさまでした。それでは、それぞれ、目的に応じて CBZ を使ってみて下さい。