炉物理プログラム演習：一点炉動特性方程式
千葉

豪

核分裂で発生する中性子は即発中性子 (prompt neutron) と遅発中性子 (delayed neutron) とに分類することが出来る。
即発中性子は、核分裂反応の結果発生した非常に高い励起状態にある核分裂片から核分裂直後に放出される中性子であり、遅発中性子は、核分裂片が即発中性子やガンマ線を放出した後の核分裂生成物もしくはその娘核種の β 崩壊に伴って放出される中性子である。核分裂生成物の崩壊は核種固有の時定数で起こるため、遅発中性子は核分裂反応の後、時間遅れを伴って発生する。自身の崩壊によって遅発中性子を発生する核分裂生成物を遅発中性子先行核 (Delayed neutronprecursor) と呼ぶ。
はじめに、遅発中性子を無視し、核分裂により発生するものが全て即発中性子であるとした場合を考える。即発中性子のみの一点炉動特性方程式は以下のように記述される。

dn(t)k−1=n(t)dtl(1)

ここで、n(t) は時刻 t における中性子密度、k は中性子増倍率、l は即発中性子寿命を示す。この式は、反応度 ρ = (k−1)/k 、中性子生成時間 Λ = l/k を導入することにより、以下のように書き直される。

dn(t)ρ= n(t)dtΛ
この方程式は容易に解くことが出来、解は以下のように与えられる。
(ρ )tn(t) = n(0) exp
Λ(2)(3)

ここで、式 (2) を数値的に解く方法について考える。n(ti ) が与えられた場合に、∆t 経過後の n(ti + ∆t) を計算するとする。式 (2) の両辺を区間 [ti , ti + ∆t] で積分すると以下の式を得る。

n(ti + ∆t) − n(ti ) =∫ρ ti +∆tn(t)dt
Λ ti(4)

陽解法と呼ばれる方法では、式 (4) 右辺の積分中の n(t) を t = ti の値で代表させる。この場合、n(ti + ∆t) は以下の式で与えられる。

(ρ )n(ti + ∆t) = 1 + ∆t n(ti )(5)Λ
一方、陰解法と呼ばれる方法では、式 (4) 右辺の積分中の n(t) を t = ti + ∆t の値で代表させる。この場合、n(ti + ∆t)は以下の式で与えられる。

(ρ )−1n(ti + ∆t) = 1 − ∆tn(ti )(6)Λ
また、陽解法と陰解法を組み合わせた θ 法と呼ばれる方法もある。この方法では、式 (4) 右辺を以下のように記述する。
∫ ti +∆tn(t)dt = ((1 − θ)n(ti ) + θn(ti + ∆t)) ∆t(7)ti

このとき、n(ti + ∆t) は以下の式で与えられる。

n(ti + ∆t) =1 + (ρ/Λ)(1 − θ)∆tn(ti )1 − (ρ/Λ)θ∆t(8)θ 法は、θ = 0 のとき陽解法と、θ = 1 のとき陰解法と一致する。また、θ = 0.5 の場合を Crank-Nicolson 法と呼ぶ。

問題１：l = 0.00005 s、k = 1.001、n(0) = 1 としたとき、t = 1 s までの中性子密度 n(t) を、陽解法、陰解法、Crank-Nicolson 法を用いていくつかの異なる時間ステップ幅で計算し、解析解 (3) と比較せよ。また、これら３つの手法の精度を比較せよ。
1これまでは即発中性子のみによる中性子の増倍を考えたが、実際には遅発中性子を考える必要がある。ここでは簡単のため、遅発中性子先行核が一種類のみ存在すると仮定する。このとき、一点炉の動特性方程式は以下のように書ける。

dn(t)
dt
dC(t)
dt

==

ρ−βn(t) + λC(t),Λβn(t) − λC(t)Λ(9)(10)

ここで、β は全核分裂中性子に対する遅発中性子の割合を、C は遅発中性子先行核密度を、λ は遅発中性子の崩壊定数を、それぞれ示す。

ρ が時間に対して依存しない場合、この連立方程式は解析的に解くことが出来る。


問題２：式 (9)、(10) の解を、ラプラス変換、ラプラス逆変換を用いて求めよ（この際、千葉のホームページの項目「原子炉動特性」のテキスト「原子炉動特性の基礎 (1) 一点炉遅発中性子一群動特性方程式」の 3 節を参考にするとよい）。なお、t = 0 で臨界である（ρ = 0）とし、t = 0 以降に反応度 ρ がステップ状に印加されるとする。中性子密度の初期値は n(0) = 1 とする（t = 0 で臨界であることから、C(0) は dC/dt = 0 から求められるであろう）。
また、n(t)、C(t) はいずれも指数関数の和で表現されることに留意せよ。得られた解を用いて、l = 0.000015 s、λ = 0.15 /s、β = 0.0065 の臨界体系に、t = 0 で ρ = 0.001 がステップ状に印加された場合の n(t)、C(t) を計算せよ。なお、Λ ≈ l としてよい。


式 (9)、(10) を陰解法で解く場合には、以下の式が得られるであろう。

n(t + ∆t) − n(t) =C(t + ∆t) − C(t) =ρ−βn(t + ∆t)∆t + λC(t + ∆t)∆t,
Λβn(t + ∆t)∆t − λC(t + ∆t)∆tΛ(11)(12)n(t + ∆t)、C(t + ∆t) を得るためには、この連立方程式を解く必要がある。θ 法の場合も同様である。
問題３：問題２と同様の問題について、n(t) と C(t) を、陽解法、陰解法、Crank-Nicolson 法で計算し解析解と比較することにより、それらの精度を比較せよ。ただし、t = 10(s) までとし、いくつかの時間ステップ幅で計算を行うこととする。


問題４：問題２と同様の体系について t = 0 以降も臨界であるとする。n(0) = 1、C(0) = 0 とした場合、直感的には系が臨界であることから n(t) = 1 と考えてしまうが、これは正しくない。そのことを実際に数値計算により確かめ、n(t) = 1 とならない理由を述べよ。


動特性方程式を数値的に解く方法としては以下に述べるものもあり、空間依存の動特性方程式を解く場合は一般的に用いられているようである。
式 (10) の両辺に exp(λt) を乗じることにより以下の式を得る。

βdC(t)
exp(λt) = n(t) exp(λt) − λC(t) exp(λt)
dtΛ

(13)

この式は以下のように書き直せる。

dβ(C(t) exp(λt)) = n(t) exp(λt)
dtΛ
この式の両辺を [ti , ti + ∆t] で積分すると以下を得る。
∫β ti +∆t ′n(t ) exp(λt′ )dt′
C(ti + ∆t) exp (λ(ti + ∆t)) − C(ti ) exp(λti ) =Λ ti(14)(15)

この式は以下のように整理される。

∫β ti +∆t ′n(t ) exp (−λ(ti + ∆t − t′ )) dt′
C(ti + ∆t) = C(ti ) exp(−λ∆t) +Λ ti
この式の右辺第二項の積分について、t̂ = ti + ∆t − t′ とおいて、以下のように書き直す。
∫ ti +∆t∫ ∆t′′′n(t ) exp (−λ(ti + ∆t − t )) dt =n(ti + ∆t − t̂) exp(−λt̂)dt̂
ti

02(16)(17)ここで、[0, ∆t] での n(ti + ∆t − t̂) での振る舞いを、以下のように一次関数で近似するものとする。なお、精度を向上させるために、より高次の関数で近似する場合もあるようである。
()()n(ti ) − n(ti + ∆t)t̂t̂n(ti + ∆t − t̂) ≈ n(ti + ∆t) +t̂ = 1 −n(ti + ∆t) +n(ti )∆t∆t∆t(18)

これにより、最終的に式 (16) は以下のように書けることになる。

C(ti + ∆t) = µC(ti ) + ηn(ti ) + ξn(ti + ∆t)(19)µ=(20)η=ξ=

ここで、exp(−λ∆t),()β1 − exp(−λ∆t)− exp(−λ∆t) ,Λλλ∆t()
1 − exp(−λ∆t)
β1−Λλλ∆t(21)(22)である。
式 (9) を [ti , ti + ∆t] で積分し、θ 法を適用すると以下のように書ける。

n(ti + ∆t) − n(ti ) =ρ−β((1 − θ)n(ti ) + θn(ti + ∆t)) ∆t + λ ((1 − θ)C(ti ) + θC(ti + ∆t)) ∆tΛ

これに対して式 (19) を代入することで、n(ti )、C(ti ) から n(ti + ∆t) が以下のように計算できる。
) }{(ρ−β(1 − θ) + λθη ∆t n(ti ) + (λ(1 − θ) + λµθ) C(ti )∆t1+Λ()n(ti + ∆t) =ρ−β1−θ∆t − λθξ∆tΛ(23)(24)

また、C(ti + ∆t) については、式 (19) から計算できる。
この方法では、n(ti + ∆t) と C(ti + ∆t) を一括して求めないことから、式 (11)、(12) のような連立方程式を解く必要がなく、コーディングが容易になるという利点があると考えられる。


問題５：問題２と同様の問題について、n(t) と C(t) を上記で述べた方法で計算し、Crank-Nicolson 法と計算精度を比較せよ。ただし、t = 10(s) までとし、いくつかの時間ステップ幅で計算を行うこととする。


以上では遅発中性子先行核を一種類として扱っているが、実際は 100 を超える種類の遅発中性子先行核が存在する。一般的な計算では、先行核の半減期に応じて複数種をグループ化（群化）して扱う（通常は 6 群で扱っている）。この場合、動特性方程式はdn(t)dt=∑ρ−βn(t) +λm Cm (t),ΛmdCm (t)dt=βmn(t) − λm Cm (t),Λ(25)(m = 1, 2, ..., M )(26)と書ける。ここで、M は遅発中性子先行核のグループ数（先行核の家系数）を示す。また、これを行列表記に直すとと書ける。なお、β =∑dn= Andt(27)

T

βm であり、n = (n C1 C2 ... CM ) である。

m

これを陰解法で解く場合にはn(t + ∆t) − n(t) = An(t + ∆t)∆t(28)

よりn(t + ∆t) = (I − A∆t)

−1n(t)(29)(I + (1 − θ)∆tA) n(t)(30)と書くことができ、θ 法で解く場合にはn(t + ∆t) = (I − θ∆tA)−13と書くことが出来る。
式 (29)、(30) のいずれにおいても、逆行列を求める必要が生じる。行列のサイズが非常に小さいため、逆行列の計算にはガウスの消去法を用いればよいであろう。
一方、逆行列を計算しない方法として、先行核１群のときと同様のものがある。この場合は、m 群の先行核密度は以下のように記述される。

Cm (ti + ∆t) = µm Cm (ti ) + ηm n(ti ) + ξm n(ti + ∆t)(31)µm(32)

ここで、=exp(−λm ∆t),()βm1 − exp(−λm ∆t)=− exp(−λm ∆t) ,Λλmλm ∆t()βm
1 − exp(−λm ∆t)=1−Λλmλm ∆tηmξmである。これらのパラメータを用いることにより、n(ti + ∆t) は以下で計算できる。
) }{(∑∑
ρ−β1+(1 − θ) +λm θηm ∆t n(ti ) +(λm (1 − θ) + λm µm θ) Cm (ti )∆tΛmm()n(ti + ∆t) =∑ρ−β1−
θ∆t −λm θξm ∆t
Λm(33)(34)(35)



問題６：問題２と同様の問題について n(t) を計算し、問題２で得られた n(t) と比較せよ。ただし、遅発中性子先行核を６群で扱い、遅発中性子データとして以下を用いることとする。

Groupλm [1/s]Relative abundance am = βm /β10.0124

0.033

20.0305

0.219

30.111

0.196

40.3010.395

51.140.115

63.010.042



式 (27) の解は形式的に次のように書くことが出来る。

n(t) = exp(At)n(0)(36)

ここで exp(At) は行列指数（matrix exponential）といい、これが計算できれば n(t) は式 (36) より直接計算することができる。
行列指数を数値的に計算する方法はいくつかあるが、ここではパーデ近似法を説明する [3]。スカラー変数 x に対する指数 exp(x) の [p/q] 次のパーデ近似は次のように書ける。

Npq (x),Dpq (x)(37)p∑(p + q − k)!p! kx ,(p + q)!k!(p − k)!

(38)q∑(p + q − k)!q!
(−x)k .
Dpq (x) =(p + q)!k!(q − k)!

(39)

exp(x) ≈ここでNpq (x) =k=0

k=0である。
同様に、行列指数に対しても以下のパーデ近似を適用することが出来る。

exp (A∆t) ≈4

Npq (A∆t)
Dpq (A∆t)(40)一般的には p = q とした近似が用いられるようである。


問題７：問題６と同じ問題をパーデ近似法を用いて計算し、その精度を陽解法、陰解法、Crank-Nicolson 法と比較せよ。
参考文献[1] 伴雄一郎、「空間依存動特性方程式の統一的解法の開発」、名古屋大学大学院工学研究科、修士論文、(2011).
[2] 深谷裕司、「次世代炉心解析システム MARBLE 用一点炉動特性ソルバー Pointkinetics の開発」、JAEA-Data/Code2011-014 (2011).
[3] 山本章夫、「燃焼の基礎理論」、第 38 回炉物理夏期セミナーテキスト、(2006).

5