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

核分裂連鎖反応で形成される中性子束の振幅と随伴中性子束
電力中央研究所 名内 泰志

1 緒言
核燃料物質内の中性子輸送で興味深い性質の一つに、「ある条件下で中性子束1分布が基本モードに漸近する」ということがある。こうした条件では中性子束分布がある固有値方程式の固有関数となっている。そして数学的にはその方程式に随伴式があり、その固有関数である随伴中性子束なる物理量が定義される。それらの中で、定常方程式の随伴中性子束は感度解析などに利用されている。
筆者は連続エネルギーモンテカルロ法（MC 法）でこれらの随伴中性子束を計算する方法を研究している[1, 2, 3]。MC 法の一般的な中性子束の計算は、中性子と原子核の衝突でおこる反応の種類、反応の結果生じる中性子のエネルギーの範囲が限定され、各々のおこる確率が決定されていることを前提としている。一方、例えば、散乱された中性子がどのような入射中性子によってもたらされたかを遡る逆問題については、エネルギーの範囲が無限大までとりえることなどがあり、候補となる反応やエネルギーをとる確率をあらかじめ定めることが難しい。積分型輸送方程式とその随伴式を対比すると、随伴中性子束の評価に、この「遡る計算」が必要と見えるところがあり、このため 2009 年まで MC 法では随伴中性子束が計算できないとされてきた。これに対し筆者は、Hurwitz や Ussachoff による、定常方程式の随伴中性子束が臨界炉心に配置した中性子源に対する”Core Power の漸近値”に比例するとした論文[4,5]に行き当たり、その比例関係が反応度の符号によらず成り立つための理論を明確にし、計算アルゴリズムを構築し、MCNP コードなどで計算した[1]。固有値方程式の随伴中性子束については Bell & Glasstone[6]が文献[5]と類似の検討をおこなっていたが、これを MCNP コードで実際に計算した。
さらに遅発中性子の影響もとりこんで Natural Mode
式=固有値方程式の随伴中性子束の物理的性質をあきらかにし、MC 法での計算を可能としてきた[2]。さらに最近、固有値方程式についての理論をとりまとめ[3]、よく使われる固有値方程式が出そろった。
炉物理部会編集小委員会ご担当から第 54 回原子力学会賞論文賞対象の論文[2]の事項を書くこと依頼されたが、第 44 回原子力学会賞論文賞対象の論文[1]のことを「炉物理の研究」
で紹介していなかったことと、せっかく論文[3]もまとめたところなので、本稿では論文[1,2, 3]の随伴中性子束に関する検討箇所に絞って、それらの比較を紹介したい。
以降、類似した数式の導出が並ぶ。数式の苦手な方や時間のない方は表 4-1 を眺めていた
だければ、筆者の意図するところをご理解いただくのに十分と思われる。

1

本稿では中性子束は angular neutron flux とする
1

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

2 固有値方程式
時間 t に依存する、位置 r, エネルギーE、飛行方向の中性子束(r,E,,t)に関する微積分
型中性子輸送方程式は以下になる[7]。巨視的全断面積t, 巨視的散乱断面積s は位置毎に組成に応じて定まるが、煩雑なので r の関数とは表記していない。
1 ∂Ψ(𝐫, E, 𝛀, t)= −𝛀 ∙ 𝛁Ψ − Σt Ψ + ∫ dEin ∫ d𝛀in Σs (Ein → E, 𝛀in → 𝛀)Ψ(Ein , 𝛀in )V∂t+ Sf (𝐫, E, 𝛀, t) + Sext (𝐫, E, 𝛀, t)(1)

式中の添字 in は incident、つまりその項が核反応を引き起こす中性子のものであることを指す。V は中性子の速さ、Sf は核分裂源、Sext は外部中性子源の発生強度である。Sf は遅発中性子の扱いでかわる。ここから、定常方程式、固有値方程式、固有値方程式、固有値方
程式を導く。よく御承知おきの読者は、図 2-1、2-2 あたりを一瞥して 2-5 節に進んでほしい。
2-1 定常方程式
本節は定常方程式について書く。ここはよく知られているので、式(7)以降の行だけ読んでほしい。時間に依存しない中性子束をs として、この定常方程式の成り立つ系の Sf は
Sf,s = ∫ dEin ∫ d𝛀in χt (E, 𝛀)νt (Ein )Σf ϕ𝑠 (Ein , 𝛀in )(2)となる。f は巨視的核分裂断面積、 t は核分裂あたりの中性子発生数で、即発中性子発生数p と、群 j 毎の遅発中性子発生数dj の和である。
νt = νp + ∑ νdj(3)jt は核分裂中性子発生スペクトルで、即発中性子と遅発中性子の関係は以下になる。
χt νt = χp νp + ∑ χdj νdj(4)jについては放出エネルギーと飛行方向で積分して 1 になるよう規格化して表記している∫ dE ∫ d𝛀χt or p or dj = 1(5)式(1)と式(2) にsmode を代入して、式(2)を式(1)に代入し、外部中性子源 Sext=0 の条件で、中性子束が時間に対して増加も減少もしないように左辺の微分項をゼロとし、各項が釣り合うように固有値 kmode で核分裂中性子生成源を除したものが定常方程式となる。
(Ein , 𝛀in )0 = −𝛀 ∙ 𝛁ϕmode− Σt ϕmode+ ∫ dEin ∫ d𝛀in Σs (Ein → E, 𝛀in → 𝛀)ϕmode
s
s
s
+

1
k mode(Ein , 𝛀in )∫ dEin ∫ d𝛀in χt (E, 𝛀)νt (Ein )Σf ϕmodes

式(6)を満たす固有値 kmode に対し中性子束smode は固有関数で、kmode とsmode の組合わせは無数にあるのだが、kmode の最大値を実効増倍率 keff という。これに相当する中性子束分布smode
は基本モードのそれ、すなわちs となる。
2(6)

炉物理の研究 第 75 号 (2023 年 3 月)𝛀 ∙ 𝛁ϕs + Σt ϕs − ∫ dEin ∫ d𝛀in Σs (Ein → E, 𝛀in → 𝛀)ϕs (Ein , 𝛀in )=1∫ dEin ∫ d𝛀in χt (E, 𝛀)νt (Ein )Σf ϕs (Ein , 𝛀in )
k eff

(7)

左辺第 1,2,3 項は中性子の漏洩、衝突、散乱源であり、左辺を方向、エネルギー、位置で積分すれば、中性子の消失数となる。右辺の積分は全核分裂中性子発生数である。Duderstadt& Hamilton は、keff 自体は核分裂中性子が次世代に生み出す核分裂中性子数の期待値、として説明しつつ、式(7)の導出で「世代」についてあまり強い関心を払っていない[7]。臨界状
態では中性子束の分布は振幅も時間に対しても世代に対しても定常なのだが、物理現象として、核分裂連鎖反応における世代は厳然と存在する。世代毎に中性子束分布をとると、十分な世代を経ていれば第 m+1 世代と第 m 世代の中性子束分布s(m+1)とs(m)は基本モードに相似で、かつ振幅の期待値の比が、臨界であれば keff=1 で一定となっている。式(7)は基本モードに中性子束が収束して以降の全ての世代の中性子束に対して成り立つ式となる。式(7)のs を基本モードに十分漸近しているという第 m 世代の中性子束s(m)で置き換えると、左辺は第 m 世代の中性子の消滅を表す式となり、右辺の積分部分はs(m)により核分裂生成す
る第 m+1 世代の核分裂中性子発生数となる。その積分値を世代間中性子発生数比である keff
で除する、という行為は、ちょうど右辺が第 m 世代の核分裂中性子発生数を示すことに等しい。したがって定常方程式は、その基本モードに形状が収束した後の任意の世代における、核分裂中性子の発生と中性子の消失のバランスを示す式である、ともいえる。この式(7)は臨界以外の状況でも成り立つ式で、中性子束の振幅が世代毎に変化しながら分布は変わらないという状況で、ある世代での中性子の生成量と消滅量のバランスを示している。ただし、未臨界であれば外部源がないと連鎖反応の継続が難しいこと、超臨界ではやがて熱的フィードバックがかかることから、式(7)を満たす中性子束を臨界以外の場で観測することは容易ではない。
2-2 固有値方程式
即発臨界未満の体系に炉心に時間的に関数状に中性子を撃ち込み、その後、遅発中性子
先行核の崩壊が顕著になるまでの時間範囲に着目すると、炉心出力が時間に対して減少し、やがて単一の指数関数に沿って減少する。中性子源の投入直後、あるいは反応度の投入直後
は、中性子束は高次モードの形状関数mode と固有値mode の積の重ね合わせとなる。
mode tΨ(𝐫, E, 𝛀, t) = ∑ ϕmode
eα
α
mode

mode e mode t の中で基本モード e t だけが他のモードに比べて減衰が遅いため、やがてこの成分が支配的になる。式(1)の Sf の遅発中性子源と外部中性子源 Sext を零にして、この時間依存中性子束を代入し、e t で除すると

3(8)

炉物理の研究 第 75 号 (2023 年 3 月)αϕ = −𝛀 ∙ 𝛁ϕα − Σt ϕα + ∫ dEin ∫ d𝛀in Σs (Ein → E, 𝛀in → 𝛀)ϕα (Ein , 𝛀in )V α+ ∫ dEin ∫ d𝛀in χp (E, 𝛀)νp (Ein )Σf ϕα (Ein , 𝛀in )(9)という、所謂固有値方程式が導出される[6]。論文や教科書によっては、測定できる中性子
束の時間変化割合が基本的に負値のため、を正値・中性子束の時間変化を e-t として表現し、式(9)左辺に負号がつく場合がある。またいくつかの論文では固有値方程式という用語を後述の natural mode 式に使っている例があり、そうした場合式(9)のをp で書くことがある。は高次モードの影響は出るものの、実効増倍率と異なり、体系が未臨界から遅発
臨界程度まで実測出来る量である。の測定にはパルス中性子法や炉雑音法等が利用される。
2-3 Natural Mode Equation あるいは固有値方程式
外部中性子源がなく、遅発中性子の発生が有意な時間領域での中性子輸送を考える際は、
固有値方程式で無視していた遅発中性子の生成項を Sf に加える。
Sf,nat = ∫ dEin ∫ d𝛀in χp (E, 𝛀)νp (Ein )Σf Ψ(𝐫, E, 𝛀, t) + ∑ χdj λj Cj

(10)

j

C は遅発中性子先行核密度、は先行核崩壊定数となる。式(1)は1 ∂Ψ(𝐸, 𝛀, t)= −𝛀 ∙ 𝛁Ψ − Σt Ψ + ∫ dEin ∫ d𝛀in Σs (Ein → E, 𝛀in → 𝛀)Ψ(Ein , 𝛀in , t)V∂t+ ∫ dEin ∫ d𝛀in χp (E, 𝛀)νp (Ein )Σf Ψ(Ein , 𝛀in , t) + ∑ χdj λj 𝐶j

(11)

j

となる。一方先行核の密度に関する式は以下となる。
dCj= ∫ dEin ∫ d𝛀in νdj (Ein )Σf Ψ(Ein , 𝛀in , t) − λj Cj
dt
Natural Mode 式は中性子束と遅発中性子先行核濃度の振幅が時間に対して同じ指数関数で(12)変化する「漸近条件」に現れる固有値方程式であり、その際の中性子束分布w 及び先行核濃度分布 cj が時間に対して不変となる。
Ψ(E, 𝛀, t) → ϕ𝜔 eωt(13)

Cj (t) → cj eωt(14)(12)に代入し et で除すると、

式(13)、(14)を式(11)、ω
ϕ = −𝛀 ∙ 𝛁ϕ𝜔 − Σt ϕ𝜔 + ∫ dEin ∫ d𝛀in Σs (Ein → E, 𝛀in → 𝛀)ϕ𝜔 (Ein , 𝛀in )V 𝜔+ ∫ dEin ∫ d𝛀in χp(E, 𝛀)νp (Ein )Σf ϕ𝜔 (Ein , 𝛀in ) + ∑ χdj λj cj

(15)

j

と
ω𝑐j = ∫ dEin ∫ d𝛀in νdj (Ein )Σf ϕ𝜔 (Ein , 𝛀in ) − λj cjが得られる。
式(16)から cj の式を作って式(15)に代入すると
4(16)

炉物理の研究 第 75 号 (2023 年 3 月)ω
ϕ = −𝛀 ∙ 𝛁ϕω − Σt ϕω + ∫ dEin ∫ d𝛀in Σs (Ein → E, 𝛀in → 𝛀)ϕω (Ein , 𝛀in )V ω+ ∫ dEin ∫ d𝛀in χp(E, 𝛀)νp (Ein )Σf ϕω (Ein , 𝛀in )

(17)

λj
+∑
∫ dEin ∫ d𝛀in χdj νdj (Ein )Σf ϕω (Ein , 𝛀in )ω + λjj

これが Natural Mode Equation あるいは固有値方程式と呼ばれる[8]。このは炉周期の逆数であり、反応度が投入された後の出力の時間に対する指数関数状の変化から得ることができる。
2-1~2-3 節の方程式の評価するところの違いを明示するため、図 2-1 を用意した。数字は世代を表し、折点は中性子の衝突点、分岐は核分裂、横破線は遅発中性子先行核を指す。世代と時間は混同されがちなのだが、例えば世代 3 の中性子の生成にかかわる連鎖反応に遅発中性子が現れるかどうかで時間的には相当「散在」することがイメージできるかと思う。
2-1 節の定常方程式では時間的な散らばりは関係なく、③の数と②の数の比、④の数の③の数に対する比で世代毎中性子発生率比を計算する。一方 2-2 節の固有値方程式は遅発中性
子が出る前の状況を示す式なので、ソースをおいてから間もない時間領域で、中性子の世代には関係なく、時間に対する中性子束の変化を扱う。また 2-3 節の固有値方程式では、遅発中性子発生が平衡状態、つまり相当時間が進んだ状況で、中性子の世代には関係なく、時間に対する中性子束の変化を扱う。
0
3
2 3
2

1

0

4

2

2
2

0

4

3

1

1
1

4

5
5

3

3
3
3

固有値方程式の成立する時間領域

4
3

2
2

2
2

5
5

番号: 世代折点: 衝突分岐: 核分裂
横破線: 遅発中性子先行核

固有値方程式の成立する時間領域
時間軸

図 2-1 核分裂連鎖反応での世代と時間の模式図
2-4 固有値方程式
従前の PWR 燃料集合体は軸方向の形状と組成が概ね一様で、無限増倍率が大きいながら、
5炉物理の研究 第 75 号 (2023 年 3 月)

水平方向の中性子漏洩率を大きくとるようにすることで、燃料が 1 体で保管されている際の未臨界性を担保している。鉛直方向2に一様で水平方向のバックリングが物質バックリン
グよりも大きい未臨界体系において、ある高さ(位置 r0 の z 成分)に中性子源を配置すると、
位置 r0 から高さ方向に十分に離れると中性子束が鉛直方向に単一の指数関数的に従って減衰することが知られる。図 2-2 に模式図を示す。この条件での水平方向中性子束分布が固有値方程式であらわされる。
まず便宜的に時間的に変動しない外部中性子源を配置するとして、中性子の輸送方程式を簡略化していくと、時間定常なので式(1)の左辺の時間微分が 0 となり、核分裂源項は定常方程式の式(2)と同形状となる。ただし keff での除算は行わない。
0 = −𝛀 ∙ 𝛁Ψ − Σt Ψ + ∫ dEin ∫ d𝛀in Σs (Ein → E, 𝛀in → 𝛀)Ψ(Ein , 𝛀in )+ ∫ dEin ∫ d𝛀in χt (E, 𝛀)νt (Ein )Σf Ψ(Ein , 𝛀in ) + S𝑒𝑥𝑡 (E, 𝛀)𝛿(𝐫, 𝐫0 )(18)

位置 r0 から離れた位置では、外部中性子源強度が 0 となる。高さ方向に十分に離れると中性子束が鉛直方向に単一の指数関数的に従って減衰することを以下で表現する。
Ψ(𝐫, E, 𝛀) → ϕ𝛾 (𝑥, 𝑦)𝑒 −γz(19)

式(18)に式(19)を代入して e-z で除すると以下が得られる。
0 = −𝛀 ∙ {(𝒆𝑥

𝜕ϕ𝛾
𝜕ϕ𝛾
+ 𝒆𝑦
)} + γ𝛀 ∙ 𝒆𝑧 ϕ𝛾 − Σt ϕ𝛾
𝜕𝑥𝜕𝑦+ ∫ dEin ∫ d𝛀in Σs (Ein → E, 𝛀in → 𝛀)ϕ𝛾

(20)

+ ∫ dEin ∫ d𝛀in χt (E, 𝛀)νt (Ein )Σf ϕ𝛾
少し整理した式、γΩ𝑧 ϕ𝛾 = 𝛀 ∙ {(𝒆𝑥

𝜕ϕ𝛾
𝜕ϕ𝛾
+ 𝒆𝑦
)} + Σt ϕ𝛾 − ∫ dEin ∫ d𝛀in Σs (Ein → E, 𝛀in → 𝛀)ϕ𝛾
𝜕𝑥𝜕𝑦

− ∫ dEin ∫ d𝛀in χt (E, 𝛀)νt (Ein )Σf ϕ𝛾
が固有値方程式と呼ばれる[9]。値を測定する実験は指数実験とよばれる[10]。

2

本稿では便宜的に鉛直方向を z 方向とし、中性子束は中性子源より上方のみを考慮する
6(21)

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

z-axis= verticalz-axisfissionhorizontal
plane

captureneutronexternal
sourcelog(flux)

y-axisx-axisnuclear material
図 2-2 指数実験条件の模式図
2-1~2-４節で 4 種類の固有値方程式を示した。少し煩雑なので、ここでは消滅演算子を使って式を簡略にして少し見通しをよくする。関数 f に対して以下の演算子を用意する。
Lf = 𝛀 ∙ 𝛁f + Σt f − ∫ dEin ∫ d𝛀in Σs (Ein → E, 𝛀in → 𝛀)f(Ein , 𝛀in )L𝑥𝑦 f = 𝛀 ∙ (𝒆𝑥𝜕f𝜕f+ 𝒆𝑦 ) + Σt f − ∫ dEin ∫ d𝛀in Σs (Ein → E, 𝛀in → 𝛀)f(Ein , 𝛀in )
𝜕𝑥
𝜕𝑦

(22)
(23)

Fp f = ∫ dEin ∫ d𝛀in χp (E, 𝛀)νp (Ein )Σf f(Ein , 𝛀in )(24)

Fdj f = ∫ dEin ∫ d𝛀in χdj νdj (Ein )Σf f(Ein , 𝛀in )(25)

Ff = Fp f + ∑ Fdj fj(26)

これらをつかうと各固有値方程式が以下であらわされる。
定常方程式
0 = −Lϕs +固有値方程式
71
Fϕk eff s(27)

炉物理の研究 第 75 号 (2023 年 3 月)αϕ = −Lϕα + Fp ϕαV αNatural mode 式あるいは固有値方程式λjωϕ𝜔 = −Lϕ𝜔 + Fp ϕ𝜔 + ∑F ϕV
ω + λj 𝑑𝑗 𝜔
j

(28)

(29)

固有値方程式γΩ𝑧 ϕ𝛾 = L𝑥𝑦 ϕ𝛾 − Fϕ𝛾

(30)

となる。
4 つの方程式をそれぞれ満たす 4 種の中性子束は当然異なった形状となる。一方で、
式(27)は核分裂中性子発生数と中性子消失数を外部源なしで釣り合わせたもので、式(28)や(29)は核分裂中性子発生数が中性子消失数を超過すれば時定数が正、つまり中性子束の振幅が時間的に増加することを意味し、また式(30)はある高さでの水平面上で核分裂中性子数より消失数が多い分を鉛直方向の中性子流で補償する式となっている。要はどれも中性子数のバランスを示している式である。
2-5 随伴方程式
式(27)~(30)の固有値方程式には共役な式が存在し、それを満たす随伴式が知られる。一般に演算子 A と関数 f に対する随伴演算子 A†と随伴関数 f†は以下の関係を満たすものとされる。
〈f † 𝐴𝑓〉 = 〈𝑓A† f † 〉(31)

かぎ括弧はエネルギー、飛行方向、空間全体での積分値となる。
式(22), (23)に現れるt は自己随伴演算子である。
続いて散乱演算子や式(24)(25)等の核分裂演算子について考える。以下添え字の out は反応で発生する中性子に関するものを表す。
∫ d𝐫 3 ∫ dEout ∫ d𝛀out f † (Eout , 𝛀out ) ∫ dEin ∫ d𝛀in Σs (Ein → Eout , 𝛀in → 𝛀out )f(Ein , 𝛀in )(32)

で積分の順番を入れ替えれば∫ d𝐫 3 ∫ dEin ∫ d𝛀in f(Ein , 𝛀in ) ∫ dEout ∫ d𝛀out Σs (Ein → Eout , 𝛀in → 𝛀out )f † (Eout , 𝛀out )(33)となる。すなわち、s と f†の積を、反応で発生する中性子のエネルギー、飛行方向で積分する、という演算子が随伴演算子となる。
次に、輸送項に現れる∙∇について述べる。スカラー関数 f とその随伴関数 f†の積に∙∇を作用させる。
𝛀 ∙ 𝛁ff † = f𝛀 ∙ 𝛁f † + f † 𝛀 ∙ 𝛁f

(34)

f † 𝛀 ∙ 𝛁f = −f𝛀 ∙ 𝛁f † + 𝛀 ∙ 𝛁ff †

(35)

移項する。
全体をエネルギー、飛行方向、位置で積分する。
〈f † 𝛀 ∙ 𝛁f〉 = −〈f𝛀 ∙ 𝛁f † 〉 + 〈𝛀 ∙ 𝛁ff † 〉8(36)

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

ここで右辺第 2 項でガウスの発散定理を利用。
〈f † 𝛀 ∙ 𝛁f〉 = −〈f𝛀 ∙ 𝛁f † 〉 + ∫ dE ∫ d𝛀 𝛀 ∙ ∫ ff † 𝐧ds(37)

ここで右辺第 2 項の dS は体系の外表面上のある微小面積で、n はその法線ベクトルである。
f に中性子束、f†に随伴中性子束をあてる場合を考える。計算境界の外からの中性子流入はゼロである。
f = 0 on dS, if 𝛀 ∙ 𝐧 < 0(38)となる。一方、その外表面から外に抜けていく中性子はその後原子炉に何も影響を与えることはないため、外向きの中性子にとっての随伴中性子束はゼロとなる。
f † = 0 on dS, if 𝛀 ∙ 𝐧 > 0(39)となる。したがって、f に中性子束、f†に随伴中性子束をあてる場合、式(37)の第 2 項がゼロとなり〈f † 𝛀 ∙ 𝛁f〉 = −〈f𝛀 ∙ 𝛁f † 〉(40)となる。随伴演算子の性質から、∙∇の随伴演算子は- ∙∇となる[11]。
これらを使って随伴演算子の定義は以下となる。
L† f † = −𝛀 ∙ 𝛁f † + Σt f † − ∫ dEout ∫ d𝛀out Σs (E → Eout , 𝛀 → 𝛀out )f † (Eout , 𝛀out )L†xy f † = −𝛀 ∙ (𝒆𝑥

𝜕f †
𝜕f †
+ 𝒆𝑦
) + Σt f †𝜕𝑥𝜕𝑦− ∫ dEout ∫ d𝛀out Σs (E → Eout , 𝛀 → 𝛀out )f † (Eout , 𝛀out )

(41)

(42)

Fp† f † = ∫ dEout ∫ d𝛀out χp (E𝑜𝑢𝑡 , 𝛀out )νp Σf f † (Eout , 𝛀out )

(43)

† †
Fdj
f = ∫ dEout ∫ d𝛀out χdj (E𝑜𝑢𝑡 , 𝛀out )νdj Σf f † (Eout , 𝛀out )

(44)

† †
F † f † = Fp† f † + ∑ Fdj
f
j

(45)

これらを使って書き下すと、以下の随伴式が得られる。
定常方程式の随伴式
0 = −L† ϕ†𝑠 +

1 † †
F ϕ𝑠
k eff

(46)

α固有値方程式の随伴式
α †ϕ = −L† ϕ†α + Fp† ϕ†αV αNatural mode 式あるいは固有値方程式の随伴式

9

(47)

炉物理の研究 第 75 号 (2023 年 3 月)λj
ω †ϕ = −L† ϕ†ω + Fp† ϕ†ω + ∑F † ϕ†V ω
ω + λj dj ω(48)

γΩz ϕ†γ = L†xy ϕ†γ − F † ϕ†γ

(49)

j

固有値方程式の随伴式

3. 中性子連鎖反応の順序に沿った重要度関数のバランス式
本章では核分裂連鎖反応の経過に沿った重要度関数の振幅について述べる。2 章では各固有値方程式の固有関数の振幅は任意の正値であったが、説明のしやすさから、基本モードについては積分値を 1 に規格化する。
〈ϕ𝑠 〉 = 1(50)〈ϕ𝛼 〉 = 1

(51)

〈ϕ𝜔 〉 = 1

(52)〈ϕγ 〉xy = ∫ dx ∫ dy ∫ dE ∫ d𝛀ϕγ = 1(53)

以降、式が並び、煩雑であるが、以降の下線①②③を考慮して式(58)(67)(84)(97)が立式され、式(65)(73)(94)(105)が同様の手順で導出されている。
3.1 十分な世代経過後の中性子束の振幅
核物質の体系があり、中性子が未だ一つも飛んでいないものとする。ここで位置 r0 にエネルギーE0、飛行方向0 の中性子を配置する。中性子核分裂連鎖反応を核分裂から核分裂までの世代で区切ることとし、最初の核分裂までの中性子の飛行を第 0 世代とし、以降核分裂で世代が進むたびに世代数を 1 加えていく。第 n 世代の中性子数の第 n-1 世代に対する中性子発生化数比を k(r0E00),n とする。世代間中性子発生率比は初期の中性子源位置に依存する。しかし、連鎖反応が繰り返され、この n の数が大きくなると、k(r0E00),n は keff に収束する。
k (𝐫0,E0,𝛀0),n → k eff
その際の中性子束の分布はs に相似となる。
Ψ(𝐫, E, 𝛀, t) → IFP (𝐫0 , E0 , 𝛀0 )ϕ𝑠 (𝐫, E, 𝛀)
右辺の(r,E,)は中性子束を評価している位相空間座標で、(r0,E0,0)は中性子源を配置している位相空間座標である。興味深い点として、中性子源をどこに配置しているかによらず、十分な世代後の子孫の中性子束分布は同じ形状になる。一方で、振幅 IFP は中性子源の位置に依存する。また、ここで中性子源としたが、外部源、核分裂源、散乱源だけでなく、飛行中の中性子すべてに IFP は定義できる。この点を考慮し図 3-１のように位置 r0 から r0+0dに
10

(54)
(55)

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

動く中性子の飛跡に着目して IFP のバランス式を作る。
位置r0 + 0 d 位置
r0 + 0 d 
衝突

透過中性子
(Eout3 ,out3)

エネルギーE0方向 0(Eout2,out2)
核分裂中性子位置(即発+ 遅発)(Eout1,out1 )r0
散乱中性子
図 3-1 IFP、I、Iのバランスを検討する中性子の飛跡

IFP を便宜的に式(54-55)条件の成り立つある世代(仮にＮとする)の振幅とする。ここでは、源にある中性子にとっての IFP は、この飛跡を通過した中性子、あるいは飛跡上での反応で生じる中性子（散乱による方向・エネルギー変換も含む）にとっての IFP の和と等しいという関係を数式化する。
中性子が微小距離 dの飛跡を通過する確率は e-td≒1- tdとなる。この飛跡上で中性子が散乱しエネルギーEout, 飛行方向 out に中性子を生じる確率は微小距離であるという条件を考慮すると
Σs (E0 → Eout , 𝛀0 → 𝛀out )dξ

(56)

この飛跡上で中性子が核分裂をおこしエネルギーEout、方向 out に中性子を生じる確率は
χt (Eout , 𝛀out )νt (E0 )Σf (E0 )dξ(57)

それぞれ位置での IFP を考慮し、バランス式を作るとIFP (𝐫0 , E0 , 𝛀0 ) = (1 − Σt dξ)IFP (𝐫0 + 𝛀0 dξ, E0 , 𝛀0 )+ ∫ dEout ∫ d𝛀out Σs (E0 → Eout , 𝛀0 → 𝛀out )dξI(𝐫0 + 𝛀0 εdξ, Eout , 𝛀out )+1∫ dEout ∫ d𝛀out χt (Eout , 𝛀out )νt (E0 )Σf (E0 )dξI(𝐫0 + 𝛀0 εdξ, Eout , 𝛀out )
k eff

(58)

ここでは
0≤ε≤1
であり、飛跡の上で反応がおきていることを指す。ここで式(58)の keff について説明が必要となる。IFP という量を、便宜的に第 N 世代後の中性子束の振幅とする。①左辺は(r0,E0,0)におかれた中性子にとっての第 N 世代後の中性子束の振幅で、右辺第１項は透過中性子にとっての第 N 世代後の中性子束、右辺第 2 項は散乱中性子束にとっての第 N 世代後の中性子束の振幅である。第 3 項については核分裂で世代が一つ進む。式(57)に IFP を乗じただけでは、第 N+1 世代の中性子束の振幅を示すこととなる。
「十分な世代」以降は振幅が世代毎にkeff 倍となる。このため、keff で除すことで第 N 世代の中性子束の振幅が得られる。すなわ
11

(59)

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

ち、核分裂源に含まれる 1/keff は第 N 世代後の中性子束の振幅のバランス式を作成するために必要となる。
続いて、式(58)の左辺をさしひく。
0 = (IFP (𝐫0 + 𝛀0 dξ, E0 , 𝛀0 ) − IFP (𝐫0 , E0 , 𝛀0 )) − dξΣt IFP (𝐫0 + 𝛀0 dξ, E0 , 𝛀0 )+ ∫ dEout ∫ d𝛀out Σs (E0 → Eout , 𝛀0 → 𝛀out )dξI(𝐫0 + 𝛀0 εdξ, Eout , 𝛀out )1+∫ dEout ∫ d𝛀out χt (Eout , 𝛀out )νt (E0 )Σf (E0 )dξI(𝐫0 + 𝛀0 εdξ, Eout , 𝛀out )
k eff

(60)

右辺第 1 項は Tayler 展開で以下になる。
IFP (𝐫0 + 𝛀0 dξ, E0 , 𝛀0 ) = IFP (𝐫0 , E0 , 𝛀0 ) + dξ𝜕IFP (𝐫0 , E0 , 𝛀0 )
𝜕ξ

(61)

右辺第 2 項の偏微分は、0 方向に変分をとることに相当する。これは方向ベクトル0 とナ
ブラ演算子の内積をとることで表現でき、
IFP (𝐫0 + 𝛀0 dξ, E0 , 𝛀0 ) = IFP (𝐫0 , E0 , 𝛀0 ) + dξ𝛀0 ∙ 𝛁IFP (𝐫0 , E0 , 𝛀0 )(62)となる。これを式(60)に代入する
0 = dξ𝛀0 ∙ 𝛁IFP (𝐫0 , E0 , 𝛀0 ) − dξΣt IFP (𝐫0 + 𝛀0 dξ, E0 , 𝛀0 )+ ∫ dEout ∫ d𝛀out Σs (E0 → Eout , 𝛀0 → 𝛀out )dξI(𝐫0 + 𝛀0 εdξ, Eout , 𝛀out )+1∫ dEout ∫ d𝛀out χt (Eout , 𝛀out )νt (E0 )Σf (E0 )dξI(𝐫0 + 𝛀0 εdξ, Eout , 𝛀out )
k eff

(63)

続いて全体を dで除す。
0 = 𝛀0 ∙ 𝛁IFP (𝐫0 , E0 , 𝛀0 ) − Σt IFP (𝐫0 + 𝛀0 dξ, E0 , 𝛀0 )+ ∫ dEout ∫ d𝛀out Σs (E0 → Eout , 𝛀0 → 𝛀out )I(𝐫0 + 𝛀0 εdξ, Eout , 𝛀out )+1∫ dEout ∫ d𝛀out χt (Eout , 𝛀out )νt (E0 )Σf (E0 )I(𝐫0 + 𝛀0 εdξ, Eout , 𝛀out )
k eff

(64)

さらに d→0 をとる。
0 = 𝛀0 ∙ 𝛁IFP (𝐫0 , E0 , 𝛀0 ) − Σt IFP (𝐫0 , E0 , 𝛀0 )+ ∫ dEout ∫ d𝛀out Σs (E0 → Eout , 𝛀0 → 𝛀out )I(𝐫0 , Eout , 𝛀out )+1∫ dEout ∫ d𝛀out χt (Eout , 𝛀out )νt (E0 )Σf (E0 )I(𝐫0 , Eout , 𝛀out )k eff

式(46)と式(65)は、IFP がs†が同一であれば同じ式となる。従って IFP が正で意味のある解をもつならば、IFP が定常方程式にとっての随伴中性子束に比例することとなる。
この式の導出は Ussachoff によって 1956 年になされている[4]。しかしながら、Ussachoffの論文では前記下線部 ① で述べた、核分裂項を keff で除した理由がかかれておらず、Importance を無限大時間後の Core Power としている[4]。図の 2-1 から類推されるように、「十分な時間後」と「十分な世代後」では意味合いが異なる。しかし Ussachoff や Hurwitz はその違いを認識していなかった、もしくは意味合いが実効的に同等となる臨界条件の随伴
中性子束にしか興味がなかったのだろうと思う。筆者はこの keff の役割を明確にするととも
12

(65)

炉物理の研究 第 75 号 (2023 年 3 月)に、連続エネルギーモンテカルロ法での随伴中性子束の計算法を明確にしている[1]。
なお、本稿では IFP について「中性子束の振幅」として説明したが、基本モードの中性子束に比例する物理量であれば何でも IFP に比例する。MCNP コードでは現状、基本モードの中性子束によって誘導される核分裂による中性子発生数を体系全体で積分したものを IFP に相当する物理量として評価している模様である。
3.2 遅発中性子発生の影響がない時間範囲で十分な時間経過後の中性子束の振幅
本節も中性子が未だ一つも飛んでいない核物質の体系で、時刻 0 に位置 r0 にエネルギー
E0、飛行方向0 の中性子を配置する。遅発中性子の発生が無視できる時間領域では、時間が十分に経過すると中性子束分布が時間の経過に対して一定のに漸近し、その際の中性子束の変化が et に従うようになる。a に式(51)で規格化をかけてしまっているので振幅を別に Iととる。ここでは Iを、中性子源を(r0,E0,0)に配置してから十分な時間を経過した後のある時刻 T 以降で、中性子束が以下の式を満たす場合の振幅とする。
Ψ(𝐫, E, 𝛀, t) → Iα (𝐫0, E0 , 𝛀0 )ϕα (𝐫, E, 𝛀)eα(t−T)(66)

3.1 節と同様の議論となるが、中性子源をどこに配置しているかによらず、十分な時間経過
後の中性子束分布形状は同一になる。ここで(r0,E0,0)に中性子を配置してから前記ある時
刻 T で形成される中性子束について、その振幅に関するバランス式を作成する。3.1 節同様、
図 3-1 でいうところの r0 から r0 + 0dにかけた飛跡での透過と衝突事象に着目する。3.1 節と異なる点は、②透過中性子や散乱中性子は(r0,E0,0)におかれた中性子に比べ発生時点に時間遅れが生じている点である。例えば r0+ 0dで衝突が起きるなら、その時刻は(r0,E0,0)に中性子を配置した時刻に比べd /V だけの時間遅れを生じる。時刻 T での中性子束の振幅に関するバランス式を作るなら、この時間遅れを補正する必要がある。この時刻では中性子束の振幅は et で変化していくので、バランス式の導出では e-・遅れ時間を乗じて補正すればよい。これを念頭に、遅発中性子発生を無視してバランス式をつくると以下となる。
αIα (𝐫0 , E0 , 𝛀0 ) = e−Vdξ (1 − Σt dξ)Iα (𝐫0 + 𝛀0 dξ, E0 , 𝛀0 )α+ e−Vεdξ ∫ dEout ∫ d𝛀out {Σs (E0 → Eout , 𝛀0 → 𝛀out )(67)+ χp (Eout , 𝛀out )νp (E0 )Σf (E0 )}dξIα (𝐫0 + 𝛀0 εdξ, Eout , 𝛀out )右辺第 1 項目を Tayler 展開し一次の項をとる。
αIα (𝐫0 , E0 , 𝛀0 ) = (1 − dξ) (1 − Σt dξ)Iα (𝐫0 + 𝛀0 dξ, E0 , 𝛀0 )Vα+ e−Vεdξ ∫ dEout ∫ d𝛀out {Σs (E0 → Eout , 𝛀0 → 𝛀out )+ χp (Eout , 𝛀out )νp (E0 )Σf (E0 )}dξIα (𝐫0 + 𝛀0 εdξ, Eout , 𝛀out )左辺を差し引く。この際、右辺第一項を展開した際の 2 次の微小項は無視する。

13

(68)

炉物理の研究 第 75 号 (2023 年 3 月)α
0 = Iα (𝐫0 + 𝛀0 dξ, E0 , 𝛀0 ) − Iα (𝐫0 , E0 , 𝛀0 ) − (Σt + ) dξIα (𝐫0 + 𝛀0 dξ, E0 , 𝛀0 )Vα+ e−Vεdξ ∫ dEout ∫ d𝛀out {Σs (E0 → Eout , 𝛀0 → 𝛀out )(69)+ χp (Eout , 𝛀out )νp (E0 )Σf (E0 )}dξIα (𝐫0 + 𝛀0 εdξ, Eout , 𝛀out )右辺第 1 項を Tayler 展開する。0∙∇の使い方は 3.1 節と同じである。
α0 = Iα (𝐫0 , E0 , 𝛀0 ) + dξ𝛀0 ∙ 𝛁Iα (𝐫0 , E0 , 𝛀0 ) − Iα (𝐫0 , E0 , 𝛀0 ) − (Σt + ) dξIα (𝐫0 + 𝛀0 dξ, E0 , 𝛀0 )Vα+ e−Vεdξ ∫ dEout ∫ d𝛀out {Σs (E0 → Eout , 𝛀0 → 𝛀out )(70)+ χp (Eout , 𝛀out )νp (E0 )Σf (E0 )}dξIα (𝐫0 + 𝛀0 εdξ, Eout , 𝛀out )全体を dで除する。
α0 = 𝛀0 ∙ 𝛁Iα (𝐫0 , E0 , 𝛀0 ) − (Σt + ) Iα (𝐫0 + 𝛀0 dξ, E0 , 𝛀0 )Vα+ e−Vεdξ ∫ dEout ∫ d𝛀out {Σs (E0 → Eout , 𝛀0 → 𝛀out )(71)+ χp (Eout , 𝛀out )νp (E0 )Σf (E0 )}Iα (𝐫0 + 𝛀0 εdξ, Eout , 𝛀out )さらに d→0 をとる。
α0 = 𝛀0 ∙ 𝛁Iα (𝐫0 , E0 , 𝛀0 ) − (Σt + ) Iα (𝐫0 , E0 , 𝛀0 )V
+ ∫ dEout ∫ d𝛀out {Σs (E0 → Eout , 𝛀0 → 𝛀out )(72)+ χp (Eout , 𝛀out )νp (E0 )Σf (E0 )}Iα (𝐫0 , Eout , 𝛀out )少し整理してαI (𝐫 , E , 𝛀 ) = 𝛀0 ∙ 𝛁Iα (𝐫0 , E0 , 𝛀0 ) − Σt Iα (𝐫0 , E0 , 𝛀0 )V α 0 0 0
+ ∫ dEout ∫ d𝛀out Σs (E0 → Eout , 𝛀0 → 𝛀out )Iα (𝐫0 , Eout , 𝛀out )(73)+ ∫ dEout ∫ d𝛀out χp (Eout , 𝛀out )νp (E0 )Σf (E0 )Iα (𝐫0 , Eout , 𝛀out )この式は式(47)と一致する。従って固有値方程式の随伴式を満たす随伴中性子束は、評
価点に中性子源をおいた際に十分な時間経過後に形成される中性子束の振幅に比例する[2]。
この式の導出は Bell & Glasstone によってなされているが[6]、上記下線部②の「時間遅れを指数関数で補正する」という考え方は入っておらず、時間依存の importance 関数 A(r,E,,t)の距離 Vt での変分をあっさり
lim [

Δt→0

A(𝐫0 + 𝛀0 VΔt, E0 , 𝛀0 , t + Δt) − A(𝐫0 , E0 , 𝛀0 , t)∂A
] = 𝛀0 ∙ 𝛁A +
VΔt
V ∂t

と書いている。
3.3 遅発中性子の発生が有意な時間領域で十分な時間経過後の中性子束の振幅
本節では、時刻ゼロで中性子が一つもない炉心の位置 r0 にエネルギーE0、飛行方向0 の
14

(74)

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

中性子源を配置した後、遅発中性子発生が所謂平衡条件に漸近する程度の時間 Tを経過した後、中性子束がΨ(𝐫, E, 𝛀, t) → Iω (𝐫0 , E0 , 𝛀0 )ϕω (𝐫, E, 𝛀)eω(t−Tω )(75)という漸近条件となる際の、中性子束の振幅 Iに関するバランス式を作成する。このバランス式の導出は、透過中性子、散乱中性子、核分裂即発中性子に関しては 3.2 節と同様なので遅発中性子について詳説する。r0 と r0+ 0 dの間の微小距離 dでの核分裂発生確率は
Σf (E0 )dξ

(76)

であり、j 群の遅発中性子先行核の発生確率は
νdj (E0 )Σf (E0 )dξ(77)

j 群の遅発中性子が時刻から+dに発生する確率は、先行核発生時点d/V からまでの間での先行核密度の減衰と崩壊定数を考慮して
dτλj eεdξ−λj (τ−)V ν (E )Σ (E )dξ
dj
0 f 0

(78)

時刻から+dに発生に発生する j 群の遅発中性子の飛行方向がout でエネルギーが Eout となる確率は−λj (τ−dτλj eεdξ)V χ (Edj
out , 𝛀out )νdj (E0 )Σf (E0 )dξ

(79)

この遅発中性子発生に起因して生じる中性子にとってある時間 Tにつくる中性子束の振③幅の期待値は、遅発中性子の発生時点がで、ゼロ秒に比べて秒遅れていることを補正して−λj (τ−dτλj eεdξ)−ωτV χ (E
dj out , 𝛀out )νdj (E0 )Σf (E0 )dξIω (𝐫0 + 𝛀0 εdξ, Eout , 𝛀out )e(80)

指数関数の項をまとめる。
εdξ

dτλj eλj V e−(ω+λj )τ χdj (Eout , 𝛀out )νdj (E0 )Σf (E0 )dξIω (𝐫0 + 𝛀0 εdξ, Eout , 𝛀out )(81)

遅発中性子の発生事象を時間的に積分する。先行核は発生時点d/V から Tまでとする3
Tω

∑∫
j

εdξ
V

εdξ

dτλj eλj V e−(ω+λj )τ ∫ dEout ∫ d 𝛀out {χdj (Eout , 𝛀out )νdj (E0 )Σf (E0 )dξIω (𝐫0(82)+ 𝛀0 εdξ, Eout , 𝛀out )}煩雑なのでの積分だけここですませておく。
εdξ

∑ eλj Vjεdξλj(e−(ω+λj ) Vω + λj− e−(ω+λj )Tω ) ∫ dEout ∫ d 𝛀out {χdj (Eout , 𝛀out )νdj (E0 )Σf (E0 )dξIω (𝐫0+ 𝛀0 εdξ, Eout , 𝛀out )}遅発中性子も含めた Iのバランス式は以下となる。

本来は遅発中性子が生じてから Tまで十分な時間をとる必要があるが、実質的には
1/min{j}の数倍の時間を過ぎれば遅発中性子発生がゼロとなり、十分な時間はとれる。
15
3

(83)

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

ωdξIω (𝐫0 , E0 , 𝛀0 ) = e− V (1 − Σt dξ)Iω (𝐫0 + 𝛀0 dξ, E0 , 𝛀0 )ωεdξ+ e− V ∫ dEout ∫ d𝛀out {Σs (E0 → Eout , 𝛀0 → 𝛀out )+ χp (Eout , 𝛀out )νp (E0 )Σf (E0 )}dξIω (𝐫0 + 𝛀0 εdξ, Eout , 𝛀out )
εdξ
εdξ
λj
+ ∑ eλj V(e−(ω+λj )ε V
ω + λj(84)j− e−(ω+λj )Tω ) ∫ dEout ∫ d 𝛀out {χdj (Eout , 𝛀out )νdj (E0 )Σf (E0 )dξIω (𝐫0+ 𝛀0 εdξ, Eout , 𝛀out )}右辺第 1 項の指数関数を Tayler 展開して第 1 項までとる。
ωIω (𝐫0 , E0 , 𝛀0 ) = (1 − (Σt + ) dξ) Iω (𝐫0 + 𝛀0 dξ, E0 , 𝛀0 )Vωεdξ+ e− V ∫ dEout ∫ d𝛀out {Σs (E0 → Eout , 𝛀0 → 𝛀out )+ χp (Eout , 𝛀out )νp (E0 )Σf (E0 )}dξIω (𝐫0 + 𝛀0 εdξ, Eout , 𝛀out )
εdξ
εdξ
λj
+ ∑ eλj V(e−(ω+λj ) Vω + λj(85)j− e−(ω+λj )Tω ) ∫ dEout ∫ d 𝛀out {χdj (Eout , 𝛀out )νdj (E0 )Σf (E0 )dξIω (𝐫0+ 𝛀0 εdξ, Eout , 𝛀out )}左辺を差し引くω
0 = Iω (𝐫0 + 𝛀0 dξ, E0 , 𝛀0 ) − Iω (𝐫0 , E0 , 𝛀0 ) − (Σt + ) dξIω (𝐫0 + 𝛀0 dξ, E0 , 𝛀0 )Vωεdξ+ e− V ∫ dEout ∫ d𝛀out {Σs (E0 → Eout , 𝛀0 → 𝛀out )+ χp (Eout , 𝛀out )νp (E0 )Σf (E0 )}dξIω (𝐫0 + 𝛀0 εdξ, Eout , 𝛀out )
εdξ
εdξ
λj
+ ∑ eλj V(e−(ω+λj ) Vω + λj(86)j− e−(ω+λj )Tω ) ∫ dEout ∫ d 𝛀out {χdj (Eout , 𝛀out )νdj (E0 )Σf (E0 )dξIω (𝐫0+ 𝛀0 εdξ, Eout , 𝛀out )}右辺第 1 項を 1 次まで Tayler 展開して第 2 項をさしひく。
ω
0 = dξ𝛀0 ∙ 𝛁Iω (𝐫0 , E0 , 𝛀0 ) − (Σt + ) dξIω (𝐫0 + 𝛀0 dξ, E0 , 𝛀0 )Vωεdξ+ e− V ∫ dEout ∫ d𝛀out {Σs (E0 → Eout , 𝛀0 → 𝛀out )+ χp (Eout , 𝛀out )νp (E0 )Σf (E0 )}dξIω (𝐫0 + 𝛀0 εdξ, Eout , 𝛀out )
εdξ
εdξ
λj
+ ∑ eλj V(e−(ω+λj ) Vω + λjj−(ω+λj )Tω−e) ∫ dEout ∫ d 𝛀out {χdj (Eout , 𝛀out )νdj (E0 )Σf (E0 )dξIω (𝐫0+ 𝛀0 εdξ, Eout , 𝛀out )}全体を dで除す。
16

(87)

炉物理の研究 第 75 号 (2023 年 3 月)ω
0 = 𝛀0 ∙ 𝛁Iω (𝐫0 , E0 , 𝛀0 ) − (Σt + ) Iω (𝐫0 + 𝛀0 dξ, E0 , 𝛀0 )Vωεdξ+ e− V ∫ dEout ∫ d𝛀out {Σs (E0 → Eout , 𝛀0 → 𝛀out )+ χp (Eout , 𝛀out )νp (E0 )Σf (E0 )}Iω (𝐫0 + 𝛀0 εdξ, Eout , 𝛀out )
εdξ
εdξ
λj
+ ∑ eλj V(e−(ω+λj ) Vω + λj(88)j− e−(ω+λj )Tω ) ∫ dEout ∫ d 𝛀out {χdj (Eout , 𝛀out )νdj (E0 )Σf (E0 )Iω (𝐫0+ 𝛀0 εdξ, Eout , 𝛀out )}さらに d→0 をとる。
0𝜔= 𝜴0 ∙ 𝜵𝐼𝜔 (𝒓0 , 𝐸0 , 𝜴0 ) − (𝛴𝑡 + ) 𝐼𝜔 (𝒓0 , 𝐸0 , 𝜴0 )𝑉
+ ∫ 𝑑𝐸𝑜𝑢𝑡 ∫ 𝑑𝜴𝑜𝑢𝑡 {𝛴𝑠 (𝐸0 → 𝐸𝑜𝑢𝑡 , 𝜴0 → 𝜴𝑜𝑢𝑡 ) + 𝜒𝑝 (𝐸𝑜𝑢𝑡 , 𝜴𝑜𝑢𝑡 )𝜈𝑝 (𝐸0 )𝛴𝑓 (𝐸0 )}𝐼𝜔 (𝒓0 , 𝐸𝑜𝑢𝑡 , 𝜴𝑜𝑢𝑡 )
+∑
𝑗

𝜆𝑗
(1
𝜔 + 𝜆𝑗(89)− 𝑒 −(𝜔+𝜆𝑗 )𝑇𝜔 ) ∫ 𝑑𝐸𝑜𝑢𝑡 ∫ 𝑑 𝜴𝑜𝑢𝑡 {𝜒𝑑𝑗 (𝐸𝑜𝑢𝑡 , 𝜴𝑜𝑢𝑡 )𝜈𝑑𝑗 (𝐸0 )𝛴𝑓 (𝐸0 )𝐼𝜔 (𝒓0 , 𝐸𝑜𝑢𝑡 , 𝜴𝑜𝑢𝑡 )}対比のために変形する。
ωI (𝐫 , E , 𝛀 ) = 𝛀0 ∙ 𝛁Iω (𝐫0 , E0 , 𝛀0 ) − Σt Iω (𝐫0 , E0 , 𝛀0 )V ω 0 0 0
+ ∫ dEout ∫ d𝛀out {Σs (E0 → Eout , 𝛀0 → 𝛀out )Iω (𝐫0 , Eout , 𝛀out )}+ ∫ dEout ∫ d𝛀out {χp (Eout , 𝛀out )νp (E0 )Σf (E0 )Iω (𝐫0 , Eout , 𝛀out )}λj+∑(1
ω + λj(90)j− e−(ω+λj )Tω ) ∫ dEout ∫ d 𝛀out {χdj (Eout , 𝛀out )νdj (E0 )Σf (E0 )Iω (𝐫0 , Eout , 𝛀out )}ここで遅発中性子先行核の性質を考慮する。遅発中性子を含む炉心の中性子束の相対時間変化率と、遅発中性子先行核の減衰定数には以下の関係があることが知られる[7]。
−min{λj } < ω(91)

ω + λj > 0(92)

この関係からとなるので、Tが十分に大きければ必ず以下の漸近性が成り立つ。
e−(ω+λj )Tω → 0

式(93)の成り立つ、つまり十分な時間が経過した後は

17

(93)

炉物理の研究 第 75 号 (2023 年 3 月)ωI (𝐫 , E , 𝛀 ) = 𝛀0 ∙ 𝛁Iω (𝐫0 , E0 , 𝛀0 ) − Σt Iω (𝐫0 , E0 , 𝛀0 )V ω 0 0 0
+ ∫ dEout ∫ d𝛀out {Σs (E0 → Eout , 𝛀0 → 𝛀out )Iω (𝐫0 , Eout , 𝛀out )}+ ∫ dEout ∫ d𝛀out {χp (Eout , 𝛀out )νp (E0 )Σf (E0 )Iω (𝐫0 , Eout , 𝛀out )}
+∑
j

(94)

λj
∫ dEout ∫ d 𝛀out {χdj (Eout , 𝛀out )νdj (E0 )Σf (E0 )Iω (𝐫0 , Eout , 𝛀out )}ω + λj

導出した式は随伴式(46)に一致する。従って Natural Mode 式、固有値方程式の随伴式を満たす随伴中性子束は、評価したい地点においた中性子源が十分な時間経過後になす中性子束の振幅に比例する。この式の導出が文献[2]の主たる成果である。この式を導出した後から式(73)を導出したところ、Bell & Glasstone の文献[6]に気付いて少々気落ちしたが、式(94)の導出には前記下線部③を考慮してのバランス式の想起が必要で、Tの要件も明示しており、それなりの差別化を図れたと自負している。
3.4 指数実験が成立する条件での中性子束の振幅指数実験が成立する条件を考える。便宜的に細長い、長手方向が一様な核物質の長手方向
を z 軸方向とし、x,y 方向の幾何学バックリングが物質バックリングよりも大きい条件を問題として設定する。ここで𝐫0 = x0 𝐞x + y0 𝐞y + z0 𝐞z(95)という位置にエネルギーE0、飛行方向0 の中性子源を配置する。中性子源の配置の仕方として 3.1 節は第 0 世代のみ、3.2, 3.3 節は時刻ゼロに瞬間的に配置としたが、この節では中性子源はおいたままとし、「世代」や「時間」に対して定常状態なものを検討する。r0 から
z 方向に十分な距離 Z 離れた位置の中性子束分布は基本モード(x,y)に漸近し、振幅は Z 離れた位置で I (x0,y0,E0, 0)になるとする。3.1~3.3 節と異なり、Iが源地点の z 座標に依存しないとしているが、そもそも指数実験条件では z 方向は長く、かつ組成が一様であることによる。
指数法条件では高さが z0+Z のあたりから中性子束が z 方向に指数関数的に減少する。
z≠z0+Z の位置まで含めて中性子束を表すと以下となる。
Ψ(𝐫, E, 𝛀, t) → Iγ (𝑥0 , 𝑦0 , E0 , 𝛀0 )ϕγ (x, y, E, 𝛀)𝑒 −γ(z−(𝑧0 +𝑍))
ここで、z=z0 上の地点 r0 から r0+ 0dに飛行する中性子に対する z=z0+Z 面上での中性子束の振幅に関するバランス式を作成する。ここで図 3-1 と類似だが、空間的な状況を想像しやすいように図 3-2 を示す。

18

(96)

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

鉛直方向

中性子束
初期中性子起因
透過中性子起因
核分裂中性子起因
散乱中性子起因
z=z0+Z面Z
透過
核分裂
r 0+ d  0d  0 e z
r 0 + d  0

d

0e z
初期
散乱
r 0 =(x 0 ,y 0 ,z 0 )

核物質
図 3-2 指数実験条件での Iのバランス式についての模式図
3.1~3.3 節と異なるのは、散乱や核分裂反応地点、透過した地点と、地点 r0 で z 方向高さが異なることに着目する点である。例えば r0+0d地点は 0∙ezdだけ高さが異なる。例によって z=z0+Z 面での振幅のバランスをとる際は、この領域で中性子束が指数関数 e-z に従って減少していくことを考慮して、補正が必要となる。その補正を加味してバランス式は以下になる。
Iγ (𝐫0 , E0 , 𝛀0 ) = eγ𝐞z∙𝛀0dξ (1 − Σt dξ)Iγ (𝐫0 + 𝛀0 dξ, E0 , 𝛀0 )+ eγ𝐞z∙𝛀0εdξ ∫ dEout ∫ d𝛀out {Σs (E0 → Eout , 𝛀0 → 𝛀out )(97)+ χt (Eout , 𝛀out )νt (E0 )Σf (E0 )}dξIγ (𝐫0 + 𝛀0 εdξ, Eout , 𝛀out )右辺第 1 項の指数部を Tayler 展開の一次の項までとる。
Iγ (𝐫0 , E0 , 𝛀0 ) = (1 + γ𝐞z ∙ 𝛀0 dξ)(1 − Σt dξ)Iγ (𝐫0 + 𝛀0 dξ, E0 , 𝛀0 )+ eγ𝐞z∙𝛀0εdξ ∫ dEout ∫ d𝛀out {Σs (E0 → Eout , 𝛀0 → 𝛀out )(98)+ χt (Eout , 𝛀out )νt (E0 )Σf (E0 )}dξIγ (𝐫0 + 𝛀0 εdξ, Eout , 𝛀out )右辺第 1 項目の二つの丸括弧の積で二次の微小量を無視する。
Iγ (𝐫0 , E0 , 𝛀0 ) = Iγ (𝐫0 + 𝛀0 dξ, E0 , 𝛀0 ) + (γ𝐞z ∙ 𝛀0 − Σt )dξIγ (𝐫0 + 𝛀0 dξ, E0 , 𝛀0 )+ eγ𝐞z∙𝛀0 εdξ ∫ dEout ∫ d𝛀out {Σs (E0 → Eout , 𝛀0 → 𝛀out )+ χt (Eout , 𝛀out )νt (E0 )Σf (E0 )}dξIγ (𝐫0 + 𝛀0 εdξ, Eout , 𝛀out )右辺第 1 項を Tayler 展開して一次の項までとる。

19

(99)

炉物理の研究 第 75 号 (2023 年 3 月)∂∂∂
+ 𝐞𝐲 + 𝐞𝐳 ) Iγ (𝐫0 , E0 , 𝛀0 )∂x∂y∂z
+ (γ𝐞z ∙ 𝛀0 − Σt )dξIγ (𝐫0 + 𝛀0 dξ, E0 , 𝛀0 )Iγ (𝐫0 , E0 , 𝛀0 ) = Iγ (𝐫0 , E0 , 𝛀0 ) + dξ𝛀0 ∙ (𝐞x+ eγ𝐞z ∙𝛀0εdξ ∫ dEout ∫ d𝛀out {Σs (E0 → Eout , 𝛀0 → 𝛀out )(100)+ χt (Eout , 𝛀out )νt (E0 )Σf (E0 )}dξIγ (𝐫0 + 𝛀0 εdξ, Eout , 𝛀out )左辺を差し引く0 = dξ𝛀0 ∙ (𝐞x∂∂∂+ 𝐞𝐲 + 𝐞𝐳 ) Iγ (𝐫0 , E0 , 𝛀0 ) + (γ𝐞z ∙ 𝛀0 − Σt )dξIγ (𝐫0 + 𝛀0 dξ, E0 , 𝛀0 )∂x∂y∂z
+ eγ𝐞z ∙𝛀0εdξ ∫ dEout ∫ d𝛀out {Σs (E0 → Eout , 𝛀0 → 𝛀out )(101)+ χt (Eout , 𝛀out )νt (E0 )Σf (E0 )}dξIγ (𝐫0 + 𝛀0 εdξ, Eout , 𝛀out )全体を dで除する。
0 = 𝛀0 ∙ (𝐞x∂∂∂+ 𝐞𝐲 + 𝐞𝐳 ) Iγ (𝐫0 , E0 , 𝛀0 ) + (γ𝐞z ∙ 𝛀0 − Σt )Iγ (𝐫0 + 𝛀0 dξ, E0 , 𝛀0 )∂x∂y∂z
+ eγ𝐞z ∙𝛀0εdξ ∫ dEout ∫ d𝛀out {Σs (E0 → Eout , 𝛀0 → 𝛀out )(102)+ χt (Eout , 𝛀out )νt (E0 )Σf (E0 )}Iγ (𝐫0 + 𝛀0 εdξ, Eout , 𝛀out )さらに d→0 をとる
0 = 𝛀0 ∙ (𝐞x∂∂∂+ 𝐞𝐲 + 𝐞𝐳 ) Iγ (𝐫0 , E0 , 𝛀0 ) + (γ𝐞z ∙ 𝛀0 − Σt )Iγ (𝐫0 , E0 , 𝛀0 )∂x∂y∂z
+ ∫ dEout ∫ d𝛀out {Σs (E0 → Eout , 𝛀0 → 𝛀out )(103)+ χt (Eout , 𝛀out )νt (E0 )Σf (E0 )}Iγ (𝐫0 , Eout , 𝛀out )ここで I(r0,E0,0)という量は、位置 r0 から z 方向に Z だけ離れた面での中性子束を示すもので、中性子源を配置した位置の z 座標が z2 なら、I(r0,E0, 0)は z=z2+Z 面の振幅を指すだけであること、また、指数実験法の条件が成り立つのは z 方向の組成が一様であることを思料すると、z 方向の偏微分はゼロとなり、r0 の関数としていたものは全て(x0,y0)の関数に置き換わる。
0 = 𝛀0 ∙ (𝐞x∂∂+ 𝐞𝐲 ) Iγ (x0 , y0 , E0 , 𝛀0 ) + (γ𝐞z ∙ 𝛀0 − Σt )Iγ (x0 , y0 , E0 , 𝛀0 )∂x∂y
+ ∫ dEout ∫ d𝛀out {Σs (E0 → Eout , 𝛀0 → 𝛀out )+ χt (Eout , 𝛀out )νt (E0 )Σf (E0 )}Iγ (x0 , y0 , E0 , 𝛀0 )比較のために整理をかける。

20

(104)

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

γ𝐞z ∙ 𝛀0 Iγ (x0 , y0 , E0 , 𝛀0 )= −𝛀0 ∙ (𝐞x∂∂+ 𝐞𝐲 ) Iγ (x0 , y0 , E0 , 𝛀0 ) + Σt Iγ (x0 , y0 , E0 , 𝛀0 )∂x∂y

− ∫ dEout ∫ d𝛀out Σs (E0 → Eout , 𝛀0 → 𝛀out )Iγ (x0 , y0 , E0 , 𝛀0 )(105)

− ∫ dEout ∫ d𝛀out χt (Eout , 𝛀out )νt (E0 )Σf (E0 )Iγ (x0 , y0 , E0 , 𝛀0 )これは式(49)と同形状となっている。したがって固有値方程式の随伴式の固有関数である
随伴中性子束(x,y,E,0)は、位置(x, y)に配置した中性子を起源として、z 方向に十分な距離をはなれた地点に形成する中性子束の振幅に比例するといえる。

4. まとめ本稿で述べてきたことをまとめると表 4-1 になる。そもそも固有値方程式は、時間依存中性子輸送方程式が、世代の経過や時間の経過、空間位置の方向に沿って同じ形状で変化する
状況、すなわち形状が基本モードに漸近して振幅が変化する状況を示すものである。これに対し、随伴中性子は、評価地点に源をおいた際に形成される基本モードの中性子束の振幅に比例する、ということとなる。また、「中性子束の振幅」としたが、その中性子束によって誘導される反応の反応率は随伴中性子束と基本的に比例する。当然、総核分裂数、総核分裂
中性子発生数、発熱量も比例する。
随伴中性子束と、基本モード中性子束の振幅の比例性を使うと、随伴中性子束は、連鎖反応の経過に沿った計算で評価できる。これにより、随伴中性子束の計算に連続エネルギーモ
ンテカルロ法が利用できる。定常方程式の随伴中性子束は、評価点に中性子源をおいて固有値計算を行い、世代毎中性子発生数比 kn をとり、それの総積をとることで得られる。固有
値方程式の随伴中性子束については、遅発中性子発生を強制的に零として、評価点からの固定源計算を行い、時間に対して中性子束を tally することで評価できる。固有値方程式も同様である。ただし遅発中性子発生を含む固定源計算を行うこと、さらには遅発中性子が平衡条件に達するまで連鎖反応をシミュレートする必要がある。筆者は文献[12]でべき乗法を利用した動特性計算で 15000~20000 世代の計算を行うことで随伴中性子束を評価しているが、計算資源を相当に要している。固有値方程式の随伴中性子束は、指数実験法が成立する体系で、評価点に点源を配置して、z 軸方向で遠方での中性子束の指数関数的減衰を計算する
ことで評価可能である。
なお本稿で説明する随伴中性子束計算はどの点に対しても計算することができるものの、
1 回の計算で 1 点のデータしか得られず、同時に広い領域を計算することはできない4。この点は、固有値と固有関数の分布を同時に計算できることに比べて煩雑なところで、効率的

4複数のソースに対して同時に計算するには「標識」をつけて、子孫の中性子に「標識」を相続させ、評価段階で標識付き中性子による中性子を計数するコーディングが必要

21

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

計算には工夫が必要となる。
筆者は従前の随伴中性子束に関する教科書[7]の説明がどうも腑に落ちなかったことがあり、3 章の検討を行って初めて随伴中性子束の意味を認識できた5。しかし、従来の随伴中性子束に関する説明をよく理解されておられる読者の方には、ぜひ本稿と見比べて新な展開をご検討いただきたい。また本稿で述べている随伴中性子束の計算法は、特に固有値問題の随伴中性子束計算等は、かなり負荷が大きい。なんでも MC 法で計算するべきだなどと考えずに、目的と差し迫った納期に応じて使いやすい道具を使っていただきたい。
表 4-1. 固有値方程式と随伴中性子束の比較
固有値方程式

定常方程式

固有値方程式

Natural Mode 式=固有値方程式

固有値方程式
固有値方程式が成

核分裂連鎖反応の即発中性子による遅発中性子と即発指数実験が可能な立する条件

世代が進み、中性

核分裂連鎖反応の中性子による連鎖

条件で、水平方向

子束分布形状が世時間範囲で、時間

反応の時間範囲での中性子束分布の代によって変わらに対して中性子束

時間に対して中性

形状が鉛直方向位

ず、振幅の世代間

形状が変わらず、

子束形状変わら置に対して変わら比が keff に収束す

振幅が時間の指数

ず、振幅が時間のず、振幅が鉛直方る条件

関数で変化する条指数関数で変化す

向位置の指数関数件る条件で変化する条件

固有値

固有値

固有値

即発中性子減衰定

炉周期の逆数指数減衰定数

固有値

実効増倍率 keff数固有値の意味固有値の測定可否

即発中性子連鎖反遅発中性子を含む

指数実験条件での毎の中性子束の振応での時間に対す連鎖反応での時間

鉛直方向位置に対幅比る中性子束の相対に対する中性子束する中性子束の相変化率の相対変化率

対変化率

未臨界～遅発臨界

即発未臨界条件で未臨界・かつ指数

条件程度まで可可実験成立条件で可遅発臨界でのみ可中性子束

s(r,E,) (r,E,)(r,E,)(x,y,E,)

計算法

べき乗法

遅発中性子を含ま

遅発中性子を含む

固有値計算、固定ない固定源時間依

時間依存動特性計

源計算

存動特性計算、-k

算、-k 法、weight法他

balance 法、他、† (r,E,)†(r,E,)

随伴中性子束5連鎖反応での世代s†(r,E,)式のフォローをされた読者の方にはお世話様と申し上げたい22†(x,y,E,)炉物理の研究 第 75 号 (2023 年 3 月)

随伴中性子束の物

評 価 す る (r,E,)

即発中性子による遅発中性子と即発指数実験が成立す理的意味に配置した中性子

核分裂連鎖反応の中性子による連鎖

る燃料において、源に誘導される、

時間範囲で、評価

反応の時間範囲

評 価 す る (x,y,E,)

十分な世代を経たする(r,E,)に配置

で 、 評 価 す るに配置した中性子後のある世代に形した中性子源に誘(r,E,)に配置した源をよって誘導さ成される定常方程導される十分な時

中性子源に誘導される、鉛直方向に式を満たす基本モ

間を経たある時間れる十分な時間を十分な距離で形成

ード中性子束の振に形成される 固経たある時間に形される γ 固有値方幅

有値方程式の基本成される固有値

程式の基本モード

モード中性子束の方程式の基本モー中性子束の振幅

振幅

ド中性子束の振幅。十分な時間の必 要 条 件 は exp((j+)t)→0

随伴中性子束の計

評価点(r,E,)に中評価点(r,E,)に中評価点(r,E,)に中

評 価 点 (x,y,E,) に算法

性子源を配置して性子源を配置し性子源を配置して中性子源を配置しべき乗法で世代間て、遅発中性子発

遅発中性子先行核て遅発中性子も含中性子発生数比 km

生を off にして固の崩壊も含めた固

めて固定源計算をを計算し、ある世定源計算である時定源計算を行い、実施し、ある高さ

間後の中性子束をある時間後の中性での中性子束を

tally する

子束を tally するtally する

Mk

代 M までm=1
を計算する

m

参考文献
[1] Yasushi Nauchi, Takanori Kameyama,"Development of calculation technique for iterated fission
probability and reactor kinetic parameters using continuous-energy Monte Carlo method", Journal of
Nuclear Science and Technology 47(11), 977–990, 2010.
[2] Yasushi Nauchi, “Adjoint flux calculation of natural mode equation by time dependent neutrontransport”, Journal of Nuclear Science and Technology, 57(8), 1000–1013, 2020.
[3] Yasushi Nauchi and Tetsuo Matsumura,” Adjoint flux and perturbation theory of γ-mode
eigenvalue problem appearing in exponential experiments”, Nuclear Science and Engineering 162(11)1306-1322, 2022.
[4] H. Hurwitz, Naval reactor physics handbook, Vol. 1, A. Radkowsky (Ed.), U. S. Atomic Energy
Commission, 864–869 (1964).
[5] L. N. Ussachoff “Equation for the importance of neutrons, reactor kinetics and the theory of
perturbations,” Proc. Int. Conf. Peaceful Uses of Atomic Energy, Geneva, Switzerland, August 8–21,
23

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

1955, Vol. 5, p.503 (1956).
[6] G. I. Bell and S. Glasstone, Nuclear reactor theory, Robert E. Krieger Publishing Company,
Malabar, Florida (1970).
[7] J. J. Duderstadt and L. J. Hamilton, Nuclear reactor analysis, John Wiley & Sons Inc., New York(1976).
[8] W. M. Stacey, Jr. Space-Time Nuclear Reactor Kinetics, Academic Press Inc, NY, US and London,
UK, (1969).
[9] T. Yamamoto and Y. Miyoshi, “An algorithm of - and -mode eigenvalue calculations by Monte
Carlo method,” Proc. 7th Int. Conf. on Nuclear Criticality Safety, (ICNC 2003), Tokai, Japan, October
20–24, 2003, JAERI-Conf 2003–019, Japan Atomic Energy Research Institute, 2003.
[10] T. Suzaki, “Subcriticality Determination of Low-Enriched UO2 Lattices in Water by Exponential
Experiment”, Journal of Nuclear Science and Technology, 28(12),1067, 1991.
[11] 小林啓祐,「原子炉物理」, コロナ社, (1996).
[12] Yasushi Nauchi, "Estimation of time-dependent neutron transport from point source based onMonte Carlo power iteration", Journal of Nuclear Science and Technology, 56(11), 996-1005, 2019.

24