炉物理の研究 第 75 号 (2023 年 3 月)

中性子多群定数処理機能(FRENDY/MG)を開発した件1[1][2]名古屋大学 山本章夫

1.はじめに多田氏を中心とする JAEA のグループが国産核データ処理コード FRENDY の開発に取り組んでいることは炉物理業界ではよく知られている[3]。FRENDY の第 1 版(V1)は 2019 年 3月にオープンソースとして公開されている。FRENDY V1 は、核データから連続エネルギー
モンテカルロコード用の ACE 形式ファイルを生成する機能を有していたが、公開された当時は「ようやく国産で NJOY を代替できる実用的なコードが開発されたか」と感慨深かった。
私も FRENDY V1(の初期バージョン)をダウンロードし、使ってみた。
FRENDY V1 は Linux
で開発されていたが、私自身はコード開発を全て Windows の Visual Studio で行っているた
め、色々苦労して FRENDY V1 を VS でコンパイルし、動かしてみた。これで ACE ファイ
ルが簡単に作れる・・・と思ったが、最初の頃は、(Linux では顕在化しなかったが)VS では
動かないメモリ関係のバグが結構あり(例：FRENDY は C++でコーディングされているが、resize していない vector 配列にデータを書き込むなど。orz)、かなりデバッグしたことは多田氏との秘密である。ちなみに、現在の FRENDY では、VS でメモリチェックも行っており、この手のバグはほぼ残っていないはずだ(異論は認める)。
さて、FRENDY V1 で ACE ファイルが作れるようになったことは、日本にとって”giantleap”2だったように思うが、処理した断面積は残念ながら MCNP や SERPENT といった連続
エネルギーモンテカルロコードでしか利用できず、いわゆる炉設計で使用する決定論的手法用の多群断面積の作成は NJOY に頼らざるを得ない状況が続いていた。この点をなんとかしたいと思ったことが多群定数作成機能の開発に取り組んだきっかけである。
評価済み核データから多群断面積を作成するのはかなり複雑な処理が必要になる。例えば NJOY コードでは、①ポイントワイズ断面積の再構成(RECONR)、②ドップラー広がり処理(BROADR)、③非分離共鳴断面積の処理(PURR)、④熱群散乱断面積の処理(THERMR)、⑤
超多群スペクトル計算と多群断面積への縮約(GROUPR)、⑥多群断面積の出力(例：MATXSR)といった処理を行う[4]。一方、FRENDY V1 で作成可能となっていた ACE ファイルに出力されるポイントワイズ断面積においては、既に①～④の処理の多くがなされている。従って、ACE ファイルを入力として多群断面積を作成すれば合理的ではないか、また、ACE ファイ
ルを入力として多群断面積を作成するコードは存在せず、使い勝手や付加価値が高くなる

第 54 回(2021 年度)の原子力学会賞特賞・技術賞を JAEA の多田健一氏、北大の千葉豪准教授とともに受賞させていただいた。受賞対象は” 純国産次世代核データ処理システム FRENDY における中性子多群断面積作成機能の開発”である。本機能の必要性や内容などは、文献[1]などに記載されており、そちらを参照頂きたい。本稿は、本機能の開発にあたっての感想を書き連
ねたものである。
2 “That’s one small step of a man, one giant leap for mankind”, Neil Armstrong, 1969.
1炉物理の研究 第 75 号 (2023 年 3 月)のではと考えた。それと前後して、別の研究プロジェクトで共鳴計算の高度化に取り組んでおり、その関係で ACE ファイルのポイントワイズ断面積を読み込んで超多群減速計算をする簡単なコードを作成していたため、それを活用できるのではとの考えもあった3。そこで、
ACE ファイルの内容を調べてみると[5]、動特性パラメータを含め、多群定数を作成するのに十分な情報があることが分かった。そもそも、MCNP などでは、様々な核特性計算を行っているのであるから当たり前の話ではある。また、ACE ファイルの内容が詳細にドキュメント化されていることもあり、ACE ファイルを入力として中性子多群定数を作成する方針とした。なお、先に述べた超多群減速計算部分は 200 行程度の短いものであり、現在のFRENDY V2 でもほぼそのまま使われている。
開発に用いているレポジトリ4のログを見ると、レポジトリの作成が 2019 年 6 月 19 日、最初のコミットは同 6 月 22 日となっている。記念すべき最初のコミットは ACE ファイルからポイントワイズ全断面積、吸収断面積、全散乱断面積、核分裂断面積を読み込む機能を有したクラスであり、合計で 400 行程度のものであった。なお、現在の多群断面積作成機能
はおそらく 2 万行弱程度ではないかと思う。
図 1 に四半期毎のコミットの推移を示す。青色のプロットが大学 PC からのコミット、赤
色が自宅/出張先からのコミット数である。これを見ると、2020 年の 1st quarter(2020 年 1 月～4 月)にコミットのピークがあるが、この時はコロナウイルス感染症により、多くの出張・会議がキャンセルになった時期である。私の実感としても、この時期は他の仕事がキャンセルされたため、かなりの時間を FRENDY/MG の開発に費やしていた記憶がある。また、2020
年の 1st quarter 以降、大学 PC からのコミット割合が増えているが、これは出張が激減し、出張中にコーディングすることがなくなったからでもある。コロナウイルス感染症の影響がこのような形で現れていることは興味深い。
FRENDY V1 との結合のための作業を 2020 年末頃(2020 年 4th quarter)からスタートさせて
おり、多群定数作成の基本的な機能はこのあたりまでで一通り完了していた。2021 年の 1stquarter と 2022 年の 1st quarter にピークがあるが、それぞれバックグラウンド断面積の自動設定機能の追加、厳密共鳴散乱取り扱い機能の追加、によるものである。

この研究プロジェクトを一緒に行っていた組織からは、オープンソースコードとしてプロジ
ェクトで使用したコードの一部を使用する許可を頂いた。改めて感謝である。
4 Subversion である。当時は github ではプライベートレポジトリを無料で作成することが出来ないとの認識であったが、今調べてみると、2019 年初め頃からプライベートレポジトリを作ることが出来るようになっていたようだ(知らなかった)。
3

炉物理の研究 第 75 号 (2023 年 3 月)

図 1 FRENDY/MG 開発に関する四半期毎のレポジトリコミット数
表 1 に FRENDY/MG 開発にあたっての主要なマイルストーンをまとめた。これを見て頂
くと、いくつか開発に時間を要している項目があることが分かる。次章では、それらの点も含めて、開発時に印象に残ったことをまとめてみたい。
表 1 FRENDY/MG の主要な開発工程
日時

開発内容

2019/6/22

最初のコミット

2019/7

ACE ファイルから角度分布などの読み込み実装

2019/8弾性散乱の取り扱い、データコンテナー実装

2019/9非弾性散乱、エネルギー分布の取り扱い実装

2019/10エネルギー分布の内挿(method of corresponding energy)実装

2019/10ν 値(即発、遅発)、核分裂スペクトル(即発、遅発)の計算実装

2019/11非分離共鳴断面積処理、バンド断面積とバンド確率から多群断面積の計算実装

2019/11

Free gas の熱群散乱断面積計算実装

2019/11

超多群計算機能、超多群スペクトルを用いた多群断面積の計算実装

2019/12

Thermal cutoff の扱いなどの実装。KRAM 形式での多群実効断面積の出力。核特性計算が一応可能に。

2020/1

計算効率化、コードのリファクタリングなど

2020/1

Free gas 熱群散乱断面積の高速化

2020/2

Free gas 熱群散乱断面積のエネルギー積分高速化

2020/3g++など Visual C++以外のコンパイラーの適合性向上。高速炉解析用に Current 重み全断面積の計算機能実装

2020/3

Free gas 以外の S(α,β)の取り扱い実装

炉物理の研究 第 75 号 (2023 年 3 月)

2020/4

検証結果を受けてのデバッグなど。特に非弾性散乱処理における実験室系と重心系の座標の取り違えの修正など。LAW=44, 66 の実装。

2020/5

Free gas 以外の S(α,β)の取り扱い。Free gas に対する Short Collision Time approximationの取り扱い実装。

2020/6

検証計算の結果を受けてのテストとデバッグ

2020/7

MATXS 形式での断面積出力実装

2020/8

Free gas 断面積の計算速度向上

2020/9-12

テストとデバッグをぼちぼち

2021/2

バックグラウンド断面積の自動設定機能実装

2021/3

GENDF 形式での断面積出力実装

2021/71000 群などエネルギー群数が多いときの所要メモリ量削減と処理高速化

2022/2

厳密共鳴散乱の取り扱い実装

2.開発中に印象に残ったこと
2.1 水素の弾性散乱の取り扱い
炉物理の教科書を開いてみると、中性子のエネルギー分布に関する記述としては、核分裂
スペクトル5に関することと、減速材核種と中性子の散乱によるエネルギーの低下、共鳴吸収による落ち込みについての記述がある。減速計算の中でも、水素と中性子の弾性散乱は最も取り扱いが簡単な例として、必ず取り上げられているのではないだろうか。水素と中性子の弾性散乱については、・共鳴領域において重心系でほぼ等方散乱、実験室系では前方散乱
・中性子が水素の原子核を「かすった」場合には実験室系での散乱角(もともとの進行方
向と散乱後の進行方向のなす角度)はほぼ 0 であり、エネルギー損失もほぼ 0 である。
・実験室系での散乱角が大きくなるにつれ、エネルギー損失が大きくなり、中性子が水素原子核と「真正面から」衝突する場合、中性子のエネルギーはゼロとなり、水素原子核
が中性子の持っていたエネルギーで弾き飛ばされる。
・つまり、散乱後の中性子のエネルギーは入射エネルギーから 0 までの値をとり得る。散乱後のエネルギー分布の確率は、単位エネルギーあたりで考えたときに一様である。
・散乱角の最大値は 90°であり、実験室系で後方散乱はしない。
・水素からなる無限媒質中の中性子スペクトルは 1/E となる余談であるが、典型的な軽水炉の中性子スペクトルをプロットすると、1MeV 近辺をピークとした核分裂に起因するスペクトルの山が見える。この山は、核分裂のエネルギー分布(やはり1MeV 近辺をピークにした形状になる)によるもの、と多くの場合考えられているようであるがそうではない。(私も長い間、同様に誤解していた)。水素や酸素の断面積がエネルギーに関して一定であれば、核分裂に起因するスペクトルの山は明確には現れないはずである。
5

炉物理の研究 第 75 号 (2023 年 3 月)などの性質がある。一言で言うと、「わかりやすく解析的に扱いやすい」現象であると言える。
しかし、中性子が陽子+電子であることを考えると、中性子の方が若干であるが陽子(水素の原子核)より重たい。中性子の重さを 1 としたとき、水素原子の重さは A=0.999167 程度の値となる。水素原子核は、全核種の中で唯一、入射粒子の中性子より軽いのである。従って、中性子が水素原子と真正面から衝突した場合、・散乱後のエネルギーはゼロにならず、次式で与えられる6。
𝐸′ = (

𝐴−1 𝟐
0.999167 − 1 𝟐
) =(
) 𝐸
𝐴+1
0.999167 + 1(1)・散乱角はゼロ、つまり散乱後、入射方向と同じ方向に飛行する。
となる。ここまでの話をまとめると、水素以外の核種の弾性散乱については、散乱後のエネ
ルギーが小さくなるほど実験室系における散乱角は大きくなるが、水素に限っては、この関係が厳密には成り立たないと言うことになる。図 2 に A=0.999167 の場合、図 3 に A=1.0, 2.0の場合の弾性散乱後のエネルギーと散乱角(いずれも実験室系)を示す。いずれも入射エネル
ギーは 100 eV とした。

A=0.999167

A=0.999167
1.0

0.8
0.4
0.2

0.0
-0.2 0

20

40

60

-0.4
-0.6

-0.8
-1.0

散乱後エネルギー(eV)

80

100

散乱角の余弦散乱角の余弦

0.6

1.0
0.8
0.6
0.4
0.2
0.0
-0.2
-0.4
-0.6
-0.8
-1.0
1.0E-05 1.0E-04 1.0E-03 1.0E-02 1.0E-01 1.0E+00 1.0E+01 1.0E+02
散乱後エネルギー(eV)

図 2 水素(A=0.999167)における弾性散乱後のエネルギーと散乱角の方向余弦。いずれも実験室系。左図は変化をわかりやすくするため、マークを追加。右図は、散乱後エネルギー0eV近辺での変化を拡大するため、横軸を対数でプロット。

厳密には、「原子核」の質量を用いるべきかもしれない。中性子の質量を 1 としたとき、陽子の質量は 0.99862 程度である。MCNP などのモンテカルロコードでは、A の値として、中性子の質量を 1 としたときの(電子の質量込みの)「原子」の質量を用いている(と思う)。ENDF で定義されている AWR も同じである。
6

炉物理の研究 第 75 号 (2023 年 3 月)

A=2.0
1.0

0.8

0.8

0.6

0.6

0.4

0.4

0.2
0.0
-0.2 0

20

40

60

-0.4

80

100

散乱角の余弦散乱角の余弦

A=1.0
1.0

0.2
0.0
-0.2 0

-0.6

-0.6

-0.8

-0.8

-1.0

20

40

60

80

100

-0.4

-1.0

散乱後エネルギー(eV)

散乱後エネルギー(eV)

図 3 A=1.0, 2.0 における弾性散乱後のエネルギーと散乱角の余弦。いずれも実験室系。
これより、水素については、低エネルギー側で散乱が極端に前方に偏ることが分かる7。
また、水素以外の核種については、実験室系での散乱角が 0～180°(散乱角の余弦は+1～-1)の間で分布するが、A=1.0 の場合、散乱角の最大値は 90°(散乱角の余弦は 0)であり、水素は散乱角の最大値は 90°より小さく、散乱角の余弦は常に正の値を取る。
多群の弾性散乱断面積を求めるためには、以下の二重積分を実行する必要がある[1]。
𝜎𝑠,𝑙,𝑔→𝑔′ =∫𝐸∈𝑔 𝜙𝑙 (𝐸) ∫𝐸′ ∈𝑔′ 𝜎𝑠 (𝐸, 𝐸 ′ , 𝜇(𝐸, 𝐸 ′ ))𝑃𝑙 (𝜇(𝐸, 𝐸 ′ ))𝑑𝐸 ′ 𝑑𝐸∫𝐸∈𝑔 𝜙𝑙 (𝐸)𝑑𝐸,(2)

where E, E’, 𝜎𝑠,𝑙,𝑔→𝑔′ , 𝜙𝑙 (𝐸) , 𝜎𝑠 (𝐸, 𝐸 ′ , 𝜇) , and 𝑃𝑙 (𝜇) are the incident energy (LAB), outgoingenergy (LAB), multi-group scattering cross section (l-th order anisotropic), l-th order angular flux
moment (Pl flux), scattering cross section from E to E’, scattering angle 𝜇 (LAB), and l-th order
Legendre polynomial 𝑃𝑙 (𝜇). Note the scattering angle is a function of E and E’.
ただし、散乱角の余弦𝜇が入射と散乱後エネルギーの関数であり、先ほど説明したように、
水素の場合は、𝜇(𝐸, 𝐸 ′ )が散乱後エネルギー𝐸 ′ ≈ 0近辺で急激に変化することを考慮しなけ
れば、精度良い数値積分が出来ない。しかし、当初はこのような事情を理解していなかった
ため、(2)式の二重積分をたいした工夫もなく8実行していた。計算結果を詳細に確認していく過程で、水素については、𝐸 ′ ≈ 0近辺で急激に散乱が前方に偏る現象があること、また、この現象が弾性散乱断面積の非等方成分として精度良く再現できていないことに気づき、どのように対応するかを考えた。その結果、中性子より軽い原子核については、以下の手順
で散乱後エネルギー𝐸 ′ についての数値積分を実施することにした。
・実験室系における散乱角の最大値を解析的に計算
・最大の散乱角に対応する散乱後のエネルギーを逆算
・上記のエネルギーを積分区間の設定の際に考慮

もっとも、低エネルギー側で散乱が前方に偏るエネルギー範囲は狭いので、A=1.0 として計算しても、核特性には大きく影響しないだろう。
8 エネルギー区間を log で等分割、Gauss–Lobatto 積分の採用などは行っていた。
7

炉物理の研究 第 75 号 (2023 年 3 月)

実験室系における散乱角の最大値は、(3)、(4)式から計算することが出来る9。詳細について
は、FRENDY Version 2 のマニュアルに記載されるはずである。
𝜇𝐶𝑀,𝑚𝑖𝑛 = −𝐴
𝜇𝐿𝐴𝐵,𝑚𝑖𝑛 =(3)

1 + 𝐴𝜇𝐶𝑀,𝑚𝑖𝑛(4)√1 + 𝐴2 + 2𝐴𝜇𝐶𝑀,𝑚𝑖𝑛

その結果、水素について、𝐸 ′ ≈ 0近辺で非等方散乱断面積成分が急激に大きくなる状況を精度良く計算出来るようになった。後で NJOY の水素の多群弾性散乱断面積を確認してみると、やはり𝐸 ′ ≈ 0近辺で非等方散乱断面積成分が急激に大きくなっており、「真面目に」散乱断面積を計算していることが確認できたことは印象的であった。
なお、二体の質点が弾性衝突する、いわゆる二体問題は、大学の初等力学で出てくるが、実際のコード開発に使おうと思うと式の導出はそれなりに面倒である。非弾性散乱については、さらに関係式が複雑になり、一筋縄ではいかなかった。多分、私が力学の勉強をサボっていたからであろう。NJOY のマニュアルを読むと、使用している式は記載されているが、導出が十分でなく、そのフォローにはそれなりの時間を要した。JAEA で実施されていた
NJOY 勉強会の資料も参考にさせていただきつつ、式の導出やフォローを行ったというのが実情である。この辺の式の導出については、文献[1]の Supplemental material として公開されている10。
2.2 ENDF/ACE ファイルにおけるエネルギーと角度分布の取り扱い
ENDF ファイルと ACE ファイルにおいては、反応後の中性子のエネルギーと角度分布は、
テーブルあるいは様々な式で与えられている。エネルギーと角度分布が独立に与えられている場合もあれば、これらが切り離せないものとして同時に与えられている場合もある。
ACE ファイルにおいては、角度は主として AND ブロックで、エネルギー分布は主としてDLW ブロックで与えられるが、整理すると表 2 となる。
表 2 ACE ファイルにおける角度及びエネルギー分布の記述
LAW(Elastic scattering, mt=2)

Angular distribution
AND blockAND blockEnergy distribution
Analytically calculated fromangular distribution
Analytically calculated fromangular distribution
DLW block

3 (Inelastic scattering, mt=5190)
4, 7, 9, 11, 66 (Fission, (n,2n)reactions, etc.)
44, 61 (Continuum inelastic
scattering, mt=91)

AND block

DLW block

DLW block

ここで、LAW は、DLW ブロックにおいて散乱後のエネルギーや角度を与えるフォーマット
9

𝜇𝐶𝑀,𝑚𝑖𝑛 は、CM 系における𝜇の最小値ではない(𝜇𝐶𝑀 の最小値は-1)。これは、「実験室系における散乱角の余弦𝜇𝐿𝐴𝐵 の最小値」を与える𝜇𝐶𝑀 である。
10 https://www.tandfonline.com/doi/suppl/10.1080/00223131.2021.1921631?scroll=top&role=tab

炉物理の研究 第 75 号 (2023 年 3 月)である。表 2 に示されているように様々なフォーマットがある。ACE ファイルの説明書によると以下の通りである[5]。
LAW=3—Level Scattering
LAW=4—Continuous Tabular Distribution
LAW=5—General Evaporation Spectrum
LAW=7—Simple Maxwellian Fission Spectrum
LAW=9—Evaporation Spectrum
LAW=11—Energy Dependent Watt Spectrum
LAW=22—Tabular Linear Functions of Incident Energy
LAW=24—Tabular Energy Multipliers
LAW=44—Kalbach-87 Formalism
LAW=61—Like LAW=44, but tabular angular distribution instead of Kalbach-87
LAW=66—N-body phase space distribution
LAW=67—Laboratory Angle-Energy Law
これらのうち、LAW=5, 22, 24, 67 は少なくとも最近の核データ(ENDF/B-VII.0 以降、
JENDL-4 以降、JEFF-3.3、TENDL-2019)を処理した ACE ファイルには現れないようだ(核データファイルにもあまり現れないのかもしれないが確認したことがないので分からない)。
いずれにせよ、角度とエネルギーの組み合わせに加え、これだけ多数のフォーマットを処理する必要があるのは面倒だ。重要なのでもう一度繰り返す。面倒である！！なお、ACE フ
ァイルについては、多くの場合これらの LAW を ENDF からコピーしており、LAW を角度
分布やエネルギー分布に再構成することは基本的に行っていない(と思う)。再構成の処理は、連続エネルギーモンテカルロコード内で実施されている。
私の理解では、これらのフォーマットは、エネルギーや角度分布をできるだけ少ないデー
タ量で表すためのものであり、計算機資源が乏しかった過去では不可欠であったであろう
ことは容易に想像できる。また、核データの評価者が評価しやすいフォーマット、という側面もあると推測している。ただし、LAW=66 のように、ENDF/B の重水素でしか使われていないものも存在する。結局のところ、評価者による評価から核データ処理のプロセスの中で、どこで詳細なエネルギー分布や角度分布を再構成するか、という問題であり、それが現時点
では核データ処理コードの所掌である、という形になっていると理解している。しかしながら、先ほどの LAW=66 のようなマイナーなものは核データ評価としては Obsolete としてはどうだろうか。エネルギー分布についても、現在の計算機資源を考えると、LAW=4 のみで良く、エネルギーと角度分布が切り離せない場合についても LAW=61 のみで良いのではないだろうか。核データ評価をやったことがない人間の勝手な感想ではあるが、このようにすることで、核データ処理コードの品質保証は大幅にやりやすくなる気がする。

炉物理の研究 第 75 号 (2023 年 3 月)

2.3 熱群散乱断面積の計算
表 1 を見ると、free gas という単語が何回も出てくる。これはすなわち、free gas の散乱断面積(incoherent inelastic scattering)の計算にかなり苦労したことを示している。ACE ファイルにおいて、熱群の散乱断面積は他の原子との化学結合や結晶構造を考えない free gas の取り扱いと化学結合や結晶構造を考慮した取り扱いの二つに大別される。後者については、散乱則𝑆(𝛼, 𝛽)を格納した専用の ACE ファイルがあり、MCNP 等ではこれを使用することになる。
例えば、H2O 内の H, 黒鉛の C 等について熱群の中性子散乱を計算する際、化学結合や結晶構造を考慮することは、熱中性子炉の計算では不可欠である。一方、化学結合・結晶構造を考慮した散乱則が与えられているのはごく一部であり、ほとんどの核種については熱群の散乱を free gas モデルで取り扱うことになる。このモデルでは、散乱断面積は解析的に以下の式で与えられる[4]。

𝜎 𝑖𝑛𝑐,𝑖𝑛𝑒𝑙𝑎 (𝐸, 𝐸 ′ , 𝜇) =𝛼=(𝛼 + 𝛽)2𝜎𝑏 𝐸 ′ 1√
exp (−),
2𝑘𝑇 𝐸 √4𝜋𝛼
4𝛼𝐸 ′ + 𝐸 − 2𝜇√𝐸𝐸 ′
,
𝐴𝑘𝑇
𝛽=

(5)

𝐸′ − 𝐸,𝑘𝑇

where
𝜎𝑏 = 𝜎𝑓𝑟𝑒𝑒(𝐴+1)2
𝐴2: the incoherent inelastic scattering cross section for the bounded state,

𝜎𝑓𝑟𝑒𝑒 : elastic scattering cross section for unbounded (free) state obtained by pointwise elastic
scattering cross section,
𝐴: relative mass of a nucleus to a neutron,𝑘: Boltzmann constant [eV/K],
𝑇: Temperature [K],𝐸: incident energy in the LAB system,
𝐸 ′ : outgoing energy in the LAB system,
𝜇: scattering cosine angle in the LAB system.
結局のところ、散乱断面積は、温度、入射/散乱後エネルギー、原子核の質量、散乱角の余弦などの関数であることから、(5)式を入射エネルギー𝐸、散乱後エネルギー𝐸 ′ 、散乱角の余弦𝜇で積分することにより、多群の散乱断面積を得ることが出来る。なお、非等方成分については、(5)式にルジャンドル関数を乗じた上で積分すれば良い。
このように、原理的には簡単な話であるが、数値計算を実際に行ってみるとそんなに簡単
ではないことが分かる。ごく簡単に言うと、熱群散乱断面積の計算は、数値積分の精度向上、

炉物理の研究 第 75 号 (2023 年 3 月)

桁落ち、オーバーフロー、アンダーフローとの戦いである。
熱群散乱断面積は、入射エネルギーと散乱後エネルギーが近い場合𝐸 ≈ 𝐸′で大きくなるが、散乱余弦が𝜇 ≈ 1の場合、(5)式から𝛼 ≈ 0となる。(5)式を見ると、𝛼が分母として二カ所に現れてきており、数値計算的に積分の精度を高めることが困難な形になっていることが分かる。実際、当初は(5)式を元に𝐸, 𝐸 ′ , 𝜇で三重数値積分を行うコードを作成したが、全く計算精度が出なかった。(5)式の𝜇依存性の一例を図 4 に示す。このように、特に𝐸 ≈ 𝐸′の場合、𝜇 ≈1近辺でデルタ関数状に値が大きくなっており、これを精度良く数値積分する必要がある。
単純な台形公式を用いた積分では、𝜇について、数万点の分割が必要になるレベルであるが、多群断面積を計算するには、𝐸, 𝐸 ′ についても数値積分を行う必要があり、𝜇についてこのような詳細な積分を実施していては日が暮れてしまう(日が暮れても計算が終わらないかも)。

図 4 𝜎 𝑖𝑛𝑐,𝑖𝑛𝑒𝑙𝑎 (𝐸, 𝐸 ′ , 𝜇)のプロット。E=0.01 [eV], E’=0.00999 [eV], A=1, T=300 [K], 𝜎𝑏 =1.0.
いろいろ調べていたところ、文献[6]においては、以下の式で熱群散乱断面積の非等方性分を計算していることが分かった11。
√𝑎+𝑏 2𝑐1
𝛽 2
𝑎 − 𝑢2
𝑖𝑛𝑐,𝑖𝑛𝑒𝑙𝑎
′)
(𝐸,
𝜎𝑙
𝐸 =∫
exp (− (𝑢 + ) ) 𝑃𝑙 (
) 𝑑𝑢,
4
𝑢
𝑏
√𝑎−𝑏 𝑏𝑎=𝐸′ + 𝐸,
𝐴𝑘𝑇

𝑏=

2√𝐸𝐸 ′
,
𝐴𝑘𝑇

𝑐=

(6)

𝜎𝑏 1
𝐸′
√ ,
4𝑘𝑇 √4𝜋 𝐸
𝑢2 = 𝛼11

文献[6]において、(6)式に相当する式においては、𝑃𝑙 (𝑎−𝑢2𝑏)の a, b が逆に表記されている。これ以外にも、熱中性子散乱則の計算を行う論文(厳密共鳴散乱の論文含む)については、式の誤りが散見される。

炉物理の研究 第 75 号 (2023 年 3 月)

2𝑐

1

𝛽 2

𝑎−𝑢2

ここで、 𝑏 exp (− 4 (𝑢 + 𝑢 ) ) 𝑃0 ( 𝑏 )の u 依存性を図 5 に示す。図 5 より、𝑢 ≈ 0近辺で値が急激に変化しているが、数値積分は図 4 に比べると容易であろうことが見て取れる。そのため、FRENDY/MG では、(6)式を用いて数値積分を行うこととした。さらに、足し算(数値積分)の回数を減らすために、数値積分のための離散化条件を動的に変化させ、被積分関数が急激に変化している部分では細かい区間で積分する adaptive integration を使用している。

図 5 (6)式の被積分関数(l=0)の u 依存性。E=0.01 [eV], E’=0.00999 [eV], A=1, T=300 [K], 𝜎𝑏 =1.0.
このように、数値積分について色々工夫を凝らしたが、それでもなお三重積分の結果得られる多群散乱断面積の精度に不安が残った。
そこで、最も重要な P0 散乱断面積については、
文献[6]に記載されている(7)式を利用して規格化することとした。高次成分は P0 散乱断面積と同じ規格化を行っている。(7)式は、(5)式を解析的に𝜇で積分して得られるものである。(7)式の計算は、数値積分でもなく、パラメータを代入すれば正確な値が得られるように思えるが、実はこの式をそのまま数値計算に利用すると、桁落ちのため適切な結果が得られないことがある。
𝑖𝑛𝑐,𝑖𝑛𝑒𝑙𝑎(𝐸, 𝐸 ′ ) =𝜎0,𝑎𝑛𝑎𝜎𝑓 𝜂 2𝐸
𝐸′
𝐸
𝐸′[exp ( − ) {erf (𝜂√ − 𝜌√ )
𝐸 2
𝑘𝑇 𝑘𝑇
𝑘𝑇
𝑘𝑇

𝐸
𝐸′
± erf (𝜂√ + 𝜌√ )}
𝑘𝑇
𝑘𝑇

(7)

𝐸′
𝐸
𝐸′
𝐸
+ {erf (𝜂√ − 𝜌√ ) ∓ erf (𝜂√ + 𝜌√ )}],
𝑘𝑇
𝑘𝑇
𝑘𝑇
𝑘𝑇
𝐴+1

𝐴−1

where erf is the error function, 𝜂 = 2√𝐴, 𝜌 = 2√𝐴, and upper/lower signs correspond to 𝐸 < 𝐸′ or
𝐸′ < 𝐸, respectively.

炉物理の研究 第 75 号 (2023 年 3 月)

例えば、T=300 K、A=16，E=1 eV のとき、0 ≤ 𝐸′ ≤ 1.0 eV に対する(7)式の[ ]内を計算す
ると、図 5 のようになる。

𝐸
𝐸
𝐸
𝐸
𝐸
𝐸′
𝐸
𝐸′
𝐸
) {erf (𝜂√ − 𝜌√ ) − erf (𝜂√ + 𝜌√ )} + {erf (𝜂√ − 𝜌√ ) + erf (𝜂√ + 𝜌√ )}]の値
図 5 (7)式の[exp (𝑘𝑇𝐸 − 𝑘𝑇
𝑘𝑇
𝑘𝑇
𝑘𝑇
𝑘𝑇
𝑘𝑇
𝑘𝑇
𝑘𝑇
𝑘𝑇
′

′

′

一方、実際に数値計算してみると、例えば𝐸 ′ = 0.6 eV のとき、以下のようになる。
[exp (

𝐸
𝐸′
𝐸
𝐸′
𝐸
𝐸′
− ) {erf (𝜂√ − 𝜌√ ) − erf (𝜂√ + 𝜌√ )}
𝑘𝑇 𝑘𝑇
𝑘𝑇
𝑘𝑇
𝑘𝑇
𝑘𝑇

𝐸′
𝐸
𝐸′
𝐸
+ {erf (𝜂√ − 𝜌√ ) + erf (𝜂√ + 𝜌√ )}]
𝑘𝑇
𝑘𝑇
𝑘𝑇
𝑘𝑇
(8)
= [5244501.9003447441114
× {0.99999999670593152823 − 1.0000000000000000000}+ {−0.95599531167630235159 + 1.0000000000000000000}]= [−0.017275748360090493309 + 0.044004688323697648411]= 0.026728939963607155102
つまりは、特に赤字部分の計算において、非常に桁落ちしやすい計算となっている。
そのため、erfc(𝑥) ≡ 1 − erf(𝑥)を用いて、(7)式に現れるerf(𝑥) + erf(𝑦)については(9)式、
erf(𝑥) − erf(𝑦)については、(10)式の形でerfc(𝑥)に置き換えて数値計算を行っている。
𝑥 × 𝑦 > 0, erf(𝑥) + erf(𝑦),𝑥 × 𝑦 < 0, erf(𝑥) + erf(𝑦) = erf(𝑥) − erf(−𝑦) = 1 − erfc(𝑥) − (1 − erfc(−𝑦))(9)= erfc(−𝑦) − erfc(𝑥),𝑥 × 𝑦 > 0, erf(𝑥) − erf(𝑦) = 1 − erfc(𝑥) − (1 − erfc(𝑦))= erfc(𝑦) − erfc(𝑥),(10)炉物理の研究 第 75 号 (2023 年 3 月)

𝑥 × 𝑦 < 0, erf(𝑥) − erf(𝑦) ,erfc(𝑥)を用いると、(11)式のように計算出来る。
[exp (

𝐸
𝐸′
𝐸
𝐸′
𝐸
𝐸′
− ) {− erfc (𝜂√ − 𝜌√ ) + erfc (𝜂√ + 𝜌√ )}
𝑘𝑇 𝑘𝑇
𝑘𝑇
𝑘𝑇
𝑘𝑇
𝑘𝑇

𝐸′
𝐸
𝐸′
𝐸
+ {erf (𝜂√ − 𝜌√ ) + erf (𝜂√ + 𝜌√ )}]
𝑘𝑇
𝑘𝑇
𝑘𝑇
𝑘𝑇
= [5244501.9003447441114(11)× {3.2940684717751522982 × 10−9
+ 2.5942036981195424107 × 10−217 }+ {−0.95599531167630235159 + 1.0000000000000000000}]= [−0.017275748360090493309 + 0.044004688323697648411]= 0.026728939963607155102
erf(𝑥)に加えてerfc(𝑥)という関数が C++で用意されている理由を分かっていなかったが、今回その必要性を理解できた。なお、桁落ちの問題については、free gas の取り扱いだけではなく、resonance up-scattering の取り扱いを実装した際[7]にも相当悩まされたのが実態である。
ここまでやればなんとかなるだろうと思っていたところ、resonance up-scattering の取り扱いを実装する際、0 K に近い温度で計算を実施する必要があり、10-10 K など、絶対零度に近い温度で計算すると断面積の計算値が異常になることが判明した。調べてみると、図 6 のよ
うに、絶対零度に近い場合、(6)式の被積分関数が、デルタ関数状になることが分かった。そのため、積分の精度が出ず、場合によっては散乱断面積が 0 として計算されていたことが判明した。そこで、(6)式の被積分関数が最大値をとる点(𝑢 = |𝛽|)を数値積分時の区間に追加することで計算精度が落ちないようにした。

図 6 (6)式の被積分関数(l=0)の u 依存性。
E=0.01 [eV], E’=0.00999 [eV], A=1, T=0.01 [K], 𝜎𝑏 =1.0.

炉物理の研究 第 75 号 (2023 年 3 月)

NJOY の THERMR では、(5)式の評価を入射エネルギー118 点に対して実施している。
GROUPR では、この結果を内挿することにより、多群の散乱断面積を計算している。しかしながら、このエネルギー入射点数は、H in ZrH など、複雑な散乱断面積を示すものがあることを考えると、不足していると考えられる。そこで、FRENDY/MG では、熱エネルギー領域に対して約 20,000 点程度の入射エネルギー点(NJOY の 100 倍以上)に対して散乱断面積の計算を行っている12。そのため、FRENDY/MG では、(6)式の被積分関数の評価を 108～109 回程度実施しており、それなりの計算時間を要しているが、熱群散乱断面積の計算精度は
NJOY よりも高いものと考えている。一例として、図 7 に Xe135 の熱群散乱断面積の非等方成分の比較を示す。入射エネルギー1eV の中性子に対する P0, P1, P2, P3 全散乱断面積の比を比較したものである。Mathematica は、(6)式を数値積分して計算したものであり、参照値となる。
図 7 より、
NJOY の結果は高次散乱成分について Mathematica との差異が大きいが、FRENDY/MG の結果は、Mathematica の結果をほぼ再現していることが分かる。
6.0E-03

4.0E-032.0E-03
P2/P0
0.0E+00-2.0E-03

1
P1/P0

2

3
P3/P0

-4.0E-03
FRENDYNJOY

Mathematica

図 7 Xe135 熱群断面積の非等方性分の比較。600K, 入射エネルギー1eV。散乱後エネルギーについて積分した値。
以上、free gas モデルを用いた多群散乱断面積の計算についての紆余曲折を示した。現在の計算アルゴリズムで計算精度は得られていると考えているが、さらに計算効率の高い手法について研究の余地があると感じている。この問題は、「(5)式を数値積分する」という単純な問題ではあるが、ここまでに示したように、「理論式は導けるが、数値計算が難しい」ということである。今後、さらに実践的な高効率アルゴリズムの開発が望まれる。
2.4 MATXSR の不具合
NJOY で処理した多群断面積は、いくつかのフォーマットで出力することが可能である。
MATXS 形式はその一つであり、汎用的なフォーマットであるためよく使用されている。
MATXS 形式の多群断面積は NJOY の MATXSR モジュールで出力されるが、特定の計算条

12

代表的な入射エネルギー点に対して散乱後のエネルギーと角度分布のテーブルを用意しておき、内挿する方式などもテストしたが、内挿では十分精度が出ず、結局詳細な積分を行っている。

炉物理の研究 第 75 号 (2023 年 3 月)

件において、多群の核分裂マトリックスが適正に出力されない不具合がある。
MATXS 形式には、散乱マトリックスなどの行列形式の断面積をそのまま格納しているわ
けではなく 0 のデータを省くなどの処理でデータの圧縮を図っている。その一つのやり方
として、一部の核分裂スペクトルのように、入射エネルギーにかかわらず放出エネルギー分
布が一定である場合、行列形式のデータを放出エネルギーと反応断面積に分離して格納する constant sub-block という方法がとられる。MATXS 形式では、この形式は 10d ブロックと呼ばれている。
MATXS 形式では、複数の温度点やバックグラウンド断面積が同一の MATXS ファイル内に含まれている場合、代表的な断面積を base の値として一式格納するが、base 以外の温度やバックグラウンド断面積点に対する断面積は、base の値からの差分を格納する。このようにすることで、温度やバックグランド断面積が変化しても変化しない断面積については、差分は 0 となり、データとして圧縮することが可能となる。
さて、10d ブロックについては、複数の温度点あるいはバックグランド断面積点が一つの
MATXS ファイルに含まれている場合、MATXSR は正しい差分を出力しない。10d ブロック
は、実態としては核分裂マトリックスに対して適用されているが、NJOY99 までは、核分裂
マトリックスが無限希釈断面積のみ(一つ)に対して出力されていたため、問題が顕在化しなかった。しかし、NJOY2016 では核分裂マトリックスがバックグランド断面積依存で出力されるようになったため、不具合が顕在化した。この不具合を回避するためには、行列形式の断面積を圧縮する手段として 10d を使用しなければ良い。この回避方法については、文献[1]の Appendix に記載したので、参考にしていただきたい。
もう一点、核データにデータが与えられているにもかかわらずデータが全てゼロの場合、
MATXS 形式で正しくデータが出力されない不具合がある。一般的に、全てがゼロの核デー
タを評価済み核データファイルに与えることはしないと思われるため、ある意味 NJOY の仕様かもしれない。この影響で、JENDL-4 の Cs-133, Cs-135, Cs-137, JEF-3.3 の Cs-135 については、MATXS 形式の多群断面積を用いると、吸収反応率に差異が生じる。この問題については、JENDL-5 において、評価済み核データに修正がなされている。
2.5 NJOY の不具合
FRENDY/MG の検証にあたっては、処理した多群断面積を用いて 1/E、軽水炉、高速炉のスペクトルを仮定して一群の全反応率、吸収反応率、散乱反応率、核分裂反応率、生成反応率を計算し、NJOY2016 を用いて同様に処理した結果と比較した。対象としたのは、ENDF/BVII.1, VIII.0, JEFF-3.3, TENDL-2019, JENDL-4.0, -4.0u, -5 の全核種である。先の MATXS 形式の不具合もその過程で発見されたものである。なお、文献[1]を執筆した際には JENDL-5 の正式版が公開されていなかったため、JENDL-5𝛼4を対象としている。なお、後に JENDL-5の正式リリース版に対して同様の検証を行い、問題が無いことを確認している。
多群定数の検証では、NJOY と開発した断面積処理コードでそれぞれ処理した多群断面積

炉物理の研究 第 75 号 (2023 年 3 月)を用いた核特性解析結果を比較するという形で実施されていることが多い。しかし、この方法では炉心解析などに現れない、あるいは影響が小さい核種について、適正に処理されているかどうかが明らかにならない問題がある。そこで、上述のように、最近の評価済み核デー
タ全てを対象とした検証を行うこととした。このような形で多群定数処理コード(NJOY を含む)が検証されたことはなく、初めての試みであったと考えている。この比較は、FRENDY/MG 開発の過程で何度も行っているが、最初に実施したときは FRENDY/MG の不具合で正常に処理が終了しなかったケースも多かったし、また、NJOY との差異が大きな核種も多かった。差異が大きな核種の処理について検討した結果、FRENDY/MG(あるいは、上流側の FRENDY V1)の処理に誤りが見つかり、修正を行ったケースも多い。
このような検討を進めていった結果、NJOY の不具合に起因すると考えられる差異が何点か明らかになっている。詳細については、文献[1]に示してあるが、概要だけ以下にまとめる。
・Eu-156 of ENDF/B-VII.1, -VIII.0、Pm-148g of JEFF-3.3、Cm-247 of JEFF-3.3：スムース断面積の内挿が粗いことに起因する差異
・Be-7 of ENDF/B-VII.1：172 XMAS 群構造の 7 群目で全断面積がゼロになる
・Ca-40 of ENDF/B-VIII.0：mt=600 の 1 MeV 近辺の断面積が適切に再構成されていない(クーロン関数の数値計算上の問題)
・Sc-45 (mt = 1) in ENDF/B-VII.1, Sn-124 (mt = 102), Po-210 (mt = 102) in ENDF/ B-VIII.0：非分
離共鳴と分離共鳴のエネルギー境界が適切に設定されない
・Ra-226 of JEFF-3.3：ν 値が適正に設定されない
NJOY のように長く使われてきている断面積処理コードでもこのように不具合がある、と見ることも出来るだろうし、あるいは断面積処理の複雑さを考えると、この程度の差異しかないのか、と考えることも出来るであろう。
2.6 投稿論文の proofreading の件
FRENDY/MG の開発そのものではないが、投稿した論文の校正原稿の出来映えが非常によろしくなかったことは、今回もっとも印象に残ったことの一つである。例えば、原稿では
𝜎𝑔,𝑥,𝑙 =∫𝐸∈𝐸 𝜎𝑥 (𝐸)𝜙𝑙 (𝐸)𝑑𝐸
𝑔∫𝐸∈𝐸 𝜙𝑙 (𝐸)𝑑𝐸

,

(12)

𝑔

という表記が proofreading の初版では以下のようになっていた。

この論文には、約 50 個の式があるが、後は推して知るべし、である。おそらくであるが、炉物理の研究 第 75 号 (2023 年 3 月)

今後、これ以下の出来映えの proofreading に出会うことはないだろうと予想している。

3. まとめFRENDY/MG の開発の経緯、印象に残ったことについて思いつくままにまとめた。エンターテイメントとして読んでいただければと考えている。
FRENDY/MG の開発は多数の方のご協力があって前進した。論文の共著の皆さんはもちろんのこと、論文発刊後に BWR において検証計算を実施してくださった GNF-J の東條氏、小野氏、共同研究の成果の利用を快諾していただいた組織、FRENDY の開発に関する
定期的な打ち合わせの場での JAEA 今野氏のアドバイス、JAEA 内で実施されていた
NJOY(GROUPR)の勉強会の資料、などが助けとなった。改めて感謝申し上げたい。
今後、FRENDY Version2 が国内はもとより海外でも広く使われ、使用経験を受けて改善されるとともに日本の炉物理/核データ処理のプレゼンス向上に貢献できることを期待して
いる。

参考文献
[1] A. Yamamoto, K. Tada, G. Chiba, T. Endo, "Multi-group neutron cross section generation
capability for FRENDY nuclear data processing code," J. Nucl. Sci. Technol., 58, pp.1165-1183(2021).
[2] A. Yamamoto, T. Endo, K. Tada, "Adaptive setting of background cross sections for generationof effective multi-group cross sections in FRENDY nuclear data processing code," J. Nucl. Sci.
Technol., 58, pp.1343-1350 (2021).
[3] K. Tada, Y. Nagaya, A. Kunieda, K. Suyama, T. Fukahoari, “Development and Verification of a
New Nuclear Data Processing System FRENDY,” J. Nucl Sci Technol., 54, 806 (2017).
[4] R. E. MacFarlane, D. W. Muir, “NJOY99.0 Code System for Producing Pointwise andMultigroup Neutron and Photon Cross Sections from ENDF/B Data.” PSR-480/NJOY99.00, Los
Alamos National Laboratory, (2000).
[5] J. Conlin [ed], ACE—A Compact ENDF—Format Specification,
https://github.com/NuclearData/ACEFormat[6] W. Zeng, R. McClarren, “Semi-analytic benchmark for multi-group free-gas Legendre moments
and the application of Gauss quadrature in generating thermal scattering Legendre moments,”
Ann. Nucl. Energy, 85, 1131 (2015).
[7] A. Yamamoto, T. Endo, G. Chiba, K. Tada, "Implementation of Resonance Up-scattering
Treatment in FRENDY Nuclear Data Processing System," Nucl. Sci. Eng., (2022).
https://doi.org/10.1080/00295639.2022.2087833