原子炉動特性の基礎 (2)
一点炉遅発中性子六群動特性方程式
千葉豪
平成 28 年 11 月 3 日

文献 [1] に基づいて、一点炉、遅発中性子六群の動特性方程式の解析解を計算する。
一点炉、遅発中性子六群の動特性方程式（以下、「遅発中性子六群」の記述は省略）は、以下のように書ける。

dn(t)
dt=6
∑
ρ−βn(t) +λi Ci (t),
Λi=1

(1)

dCi (t)dt=

βin(t) − λi Ci (t),Λ(2)(i = 1, ..., 6)

この連立常微分方程式の一般解は、以下のように指数関数の重ね合わせとなることが分かっている。

n(t) =7∑

Aj exp(sj t)

(3)

j=1

Ci (t) =7∑

Cij exp(sj t),(i = 1, ..., 6)

(4)

j=1

また、指数関数の肩 sj は、以下の s に対する方程式（逆時間方程式）を満足する。

ρ=6sl1 ∑sβi+sl + 1 sl + 1 i=1 s + λi(5)式 (5) は s についての 7 次方程式となるため、ρ が与えられた場合、解が 7 つ（s1 、· · ·、s7 ）求まることになる。
ここで、式 (3)、(4) 中の係数 Aj 、Cij について考えよう。式 (4) を式 (2) に代入することにより、以下の関係式が得られる。
(7∑βiexp (sj t) (sj + λi ) Cij − Aj
Λ
j=1)

= 0,

(i = 1, ..., 6)(6)

この式が恒等的に成り立つことから、以下の関係式が得られる。

Cij =βiAj ,Λ(sj + λi )(i = 1, ..., 6, j = 1, ..., 7)(7)

従って、式 (2) は以下のように書き直せる。
7dCi (t)βi ∑Aj=exp (sj t) ,dtΛ j=1 sj + λi

1(i = 1, ..., 6)(8)係数 Aj については、n と Ci に対する初期条件 n(0)、Ci (0) を用いて、以下の式を満足するように一意的に決定される。

n(0) =

7∑Aj ,(9)j=1
7βi ∑Aj,Λ j=1 sj + λi

Ci (0) =(i = 1, ..., 6)(10)

仮に初期条件として臨界定常状態を仮定する場合には、

Ci (0) =βin(0),Λλi(i = 1, ..., 6)(11)より、式 (10) は以下のように書き直せる。

n(0) =7∑λiAj ,s + λij=1 j(i = 1, ..., 6)(12)

Keepin の U-235 熱中性子核分裂のデータを用いて計算を行った。定数データを Table 1 に示す。

Table 1: Constants for numerical calculations
Group123456
l[s]βeﬀ

λi [1/s]
0.01240.0305
0.111
0.301
1.14
3.01

Relative abundance ai
0.033
0.219
0.196
0.395
0.115
0.042
0.0000492
0.007611

ai /λi
2.66
7.18
1.77
1.31
0.100.01

定常臨界状態に対して t = 0 において ρ = 0.6 [$] の反応度を印加したときの中性子数と遅発中性子先行核密度の時間変化を計算した結果を Fig. 1 に示す。中性子数の図から明らかなように、反応度が挿入された直後に瞬間的な跳躍があり、その後は緩やかに増加していく様子が分かる。このときの逆時間方程式の根 si と係数 Ai を Table 2 に示す。負の根に対応する係数が負となっていることが分かる。
中性子数の根毎の寄与、及び 1 番目の根の寄与に対する各根の比を Fig. 2 に示す。1 番目の正の根が漸近状態の振る舞いを与えること、即発跳躍は 7 番目の根が大きく寄与していることが分かる。また、時定数が長く、Ai が比較的大きな値となる 3 番目の根の成分が比較的長く残ることが分かる。ただし、反応度印加 15 秒後程度でその寄与も無視できる程度となる。
また、印加反応度を即発臨界未満の範囲内で変えたときの、A2 から A7 の反応度に対する依存性を計算したものを Fig. 3 に示す。全ての係数が常に負の値となっていることが分かる 1 。
1

これについては、数値的に計算しなくとも、解析的に証明することが可能であると考えられる。

2

(b) Delayed neutron precursors
Number of delayed neutron precusors(a) Number of neutrons
6

Number of neutrons

5432
1000.5

1Time [s]

1.5

2

8000

C1C2C3C4C5C67000
6000
5000
4000
3000
2000
1000
002

46810Time [s]Fig. 1: Time dependences of numbers of neutrons and delayed neutron precursors (ρ = 0.6 [$])
Table 2: Roots of inhour equation and coeﬃcients (ρ = 0.6 [$])
i
123
4
567

si [1/s]
0.299-0.013-0.043-0.153-0.908-2.768-63.182Ai
3.308-0.031-0.283-0.212-0.233-0.097-1.452(a) Absolute contribution(b) Relative contribution to the 1st Comp.
0Ratio to the first component

5Number of neutrons

43

i=1
i=2
i=3
i=4
i=5
i=6
i=7

2
1
0-1
-2

-0.02
-0.04
-0.06

i=2
i=3
i=4
i=5
i=6
i=7-0.08
-0.1

00.05

0.1

0.15

0.2 0.25
Time [s]

0.3

0.35

0.4

02

4

68
Time [s]101214

Fig. 2: Time dependence of component-wise neutrons (ρ = 0.6 [$])

参考文献
[1] 外池幸太郎、内山軍蔵、「一点炉動特性モデルを適用した臨界実験装置シミュレータ（CASIM）」、JAEA-Data/Code 2009-013、日本原子力研究開発機構、(2013).
[2] 三澤毅、他、「原子炉物理実験」、京都大学学術出版会、(2010).

30Coefficients

-0.2
-0.4
-0.6
i=2
i=3i=4i=5i=6i=7-0.8-1-1.200.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008
Reactivity

Fig. 3: Reactivity dependence of coeﬃcients of exponential functions[3] C. Meyer, Matrix analysis and applied linear algebra, Society for industrial and Applied
Mathematics, Philadelphia (2000).
[4] 深谷裕司、「次世代炉心解析システム MARBLE 用一点炉動特性ソルバー PointKinetics の開発」、JAEA-Data/Code 2011-014 (2011).
[5] 伴雄一郎、「空間依存動特性方程式の統一的解法の開発」、名古屋大学大学院工学研究科、修士論文、(2011).
[6] S. Yamoah, et al., ‘An accurate solution of point reactor neutron kinetics equations ofmulti-group of delayed neutrons,’ Ann. Nucl. Energy, 54, p.104-108 (2013).

4A

行列指数を用いた動特性方程式の解析解の導出

ここでは、行列指数を用いた動特性方程式の解析解の導出について、文献 [3, 4] の記述に基づいて説明する。
式 (1)、(2) に示した動特性方程式は以下のように行列形式で書ける。

dn= Andt(13)

ここで、n = (n, C1 ..., C6 )T であり、
(ρ − β)/Λ λ1λ2λ3λ4λ5λ6β1 /Λ
−λ100000 β2 /Λ0−λ20000 

A=
β3 /Λ
00−λ3000 



β/Λ
000−λ0044
β5 /Λ
0000−λ50 


β6 /Λ
00000−λ6(14)である。
式 (13) の解は行列指数を用いて以下のように書ける。

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

行列指数 exp(At) は、指数関数のテーラー展開と同一の形式で以下のように定義される。

exp(At) = I + At +(At)2 (At)3++ ...
2!
3!

(16)

−1
ここで、行列 A が対角化可能である場合、si を行列 A の固有値とすると、
At ==() PDP
k −1−1kkk−1Pdiag(s1 t, ..., sn t)P と書けるので、(At) = PD P = Pdiag (s1 t) , ..., (sn t) P となり、

eAt

=∞
∑
(At)kk=0k!

=∞∑
P D k P −1k=0k!

=P(∞)∑ Dkk!
k=0

P −1 = P diag(es1 t , ..., esn t )P −1(17)と書ける。このことは、行列 A が対角化可能であれば、行列指数は行列 A の固有値を肩にもつ指数関数の重ね合わせで記述できることを示している。
それでは次に、動特性方程式に現れる行列 A が対角化可能であるかについて考えよう。
一般に、n × n 行列が n 個の distinct な固有値をもつ場合、その行列は対角化可能であると言える。行列 A の固有値は以下の条件式より求められる。

det |A − sI| = 0
行列 (A − sI) は(18)(ρ − β)/Λ − s
λ1λ2λ3λ4λ5λ6
β1 /Λ
−λ1 − s00000
β2 /Λ
0−λ2 − s0000

A − sI = β3 /Λ
00−λ3 − s000β/Λ000−λ−s0044β5 /Λ
0000−λ5 − s0β6 /Λ00000−λ6 − s(19)5と書けるが、特徴的な疎な構造をもっているため、その行列式は以下のように書くことが出来る。
(det|A − sI| =)∏∑ βi ∏ρ−β−s(−λj − s) −λi(−λj − s)ΛΛ j̸=iji

これより(det|A − sI| = 0 =(20))

∑βi
ρ−βλi
−s +Λ
Λ(λi + s)i(21)が得られ、Λ = (1 − ρ)l を代入すると、逆時間方程式 (5) が得られる。式 (5) は s について 7 つの特異点をもち、その負側、正側の極限接近値に対応する ρ の値がそれぞれ −∞、∞ となること、
s → ∞ および s → −∞ の極限で ρ → 1 となることから [5]、式 (5) は、ρ ̸= 1 の場合、7 つの解を持つことが言える。行列 A が 7 つの distinct な解をもつことから、行列 A は対角化可能であり、式(13) の解は 7 つの指数関数の重ね合わせで記述されることが言える。
以上の議論から、n は次のような形で書けることが分かる。

n(t) = Pdiag(es1 t , ..., es7 t )P−1 n(0)(22)

これを変形すると、以下の式を得る。

P−1 n(t) = diag(es1 t , ..., es7 t )P−1 n(0)(23)

ここで、n を P−1 によって線形変換したベクトル n̄ = P−1 n を考えると、式 (23) は以下のように書ける。
n̄(t) = diag(es1 t , ..., es7 t )n̄(0)
(24)
すなわち、

n̄i (t) = n̄i (0) exp(si t)(25)と書け、n̄i は単一の指数関数で記述されることが分かる。
ここで、si は値が大きい順に並んでいるものとする（s1 が安定ペリオドに対応）。n = Pn̄ と書けること、n̄ が単一の指数関数であることから、Pij は j 番目の指数関数成分の ni への寄与を示していることが分かる。例えば列ベクトル P∗1（安定ペリオドの固有値に対応する固有ベクトル）
は、安定ペリオドに対応する指数関数成分の、各々の n に対する寄与を示していると言える。
固有値 si に対応する固有ベクトル ui は以下のように容易に計算することが可能である [6]。固()T
有ベクトル ui を ui = 1, ui1 , ui2 , ..., ui6 と書くと、以下の式が成り立つ。
 ρ−β

− siΛ
β1 /Λ
β2 /Λ
..
.
β6 /Λ

λ1···

λ2λ6

−λ1 − si0···00−λ2 − si · · ·0..
..
..
..
.
.
.
.
00· · · −λ6 − si

これより、uij =βjΛ(λj + si )







10

  ui  0 
 1 
 i   u2 
= 0 
  ..   ..  .   . ui6(26)0(27)が得られる。
行列 A の固有値、固有ベクトルを用いた一点炉動特性方程式の解法については文献 [6] において試みられている。実用性が高いとは言えないが、例えば、いくつかの反応度 ρ に対して固有値、6固有ベクトルを用意しておき、任意の ρ の値に対する固有値、固有ベクトルはそれらから内外挿により計算し、逐次、適切な線形変換を施して n の時間挙動を指数関数より求める、という方法が考えられるであろう。
最後に、式 (3)(4) における Aj 、Cij と P の関係について考えよう。式 (3)(4) に基づくと、n(t)はn(t) = Gf (t)と書ける。ここで、A1 A2
 C11 C12G =  C21 C22
 .
..
 .
.
 .
C61 C62(28)· · · A7
· · · C17 · · · C27 .. ..
.
. · · · C67(29)であり、f (t) = (es1 t , ..., es7 t )T である。一方、P を用いた場合、n(t) はn(t) = Pdiag(es1 t , ..., es7 t )n̄(0) = Pdiag(n̄1 (0), ..., n̄7 (0))f (t)(30)とも書くことが出来る。従って、

G = Pdiag(n̄1 (0), ..., n̄7 (0))の関係が成り立つことが分かる。

7(31)

B

原子炉安定ペリオドと反応度の関係

式 (5) をペリオド T = 1/s で記述すると、

ρ=

6lT ∑βi+T + l T + l i=1 1 + λi T(32)と書ける。これが原子炉安定ペリオドと反応度の関係を表す式となる。また、|T | ≫ l の場合には簡略化して
6∑βi
ρ≈
(33)
1+λi Ti=1と書け、この場合にはドル単位の反応度は66∑
βi /βaiρ ∑≈=β1 + λi T1 + λi Ti=1i=1(34)と書ける。
次に、与えられたペリオドに対して遅発中性子 6 群定数を用いて計算される反応度を再現する
1 群に丸めた遅発中性子の崩壊定数 λ̄ を求めてみよう。
λ̄ は遅発中性子先行核の「平均的な」崩壊定数に対応する。ペリオドを T 、再現すべき反応度をρ とすると、λ̄ は以下の式を満足すべきである。

ρ=lTβ+·T + l T + l 1 + λ̄T(35)

この式より、λ̄ が以下のように与えられる。

λ̄ =1β−(T + l)ρ − l T(36)

この式に基づいて計算した λ̄ を Fig. 4 に示す。値は 0.08 から 0.2 程度を示しており、6 群定数では 2、3 群に対応している。これより、ペリオドが長くなる、すなわち印加される反応度が小さくなるほど、遅発中性子先行核の平均的な崩壊定数は小さくなっていくことが分かる。

Equivalent decay constant [1/s]0.190.180.170.160.150.140.130.120.110.100.090.08020406080 100 120 140 160 180 200
Period [s]Fig. 4: One-group equivalent decay constant8

C

原子炉安定ペリオドから決定される反応度の遅発中性子先行核定数に対する感度

原子炉安定ペリオドから決定される反応度 ρ の λi 、βi に対する微係数は、式 (33) より、以下のように書ける。

∂ρ∂λi /λi∂ρ∂βi /βi

βi λi T,
(1 + λi T )2βi
1 + λi T= −(37)=(38)

ペリオドが正の場合における微係数 ∂ρ/(∂λi /λi ) の計算結果を Fig. 5 と Fig. 6 に、∂ρ/(∂βi /βi )の計算結果を Fig. 7 と Fig. 8 に、それぞれ示す。ペリオドが比較的長い場合には 2 群の感度が大きいことが分かる。
では、T → ∞ の極限をとった場合には、各群の定数の感度の大小はどのようになるであろうか？
上式から明らかなように、T → ∞ では相対感度はゼロとなるため、ここでは感度に T を乗じて考える。すると、∂ρTT→∞ ∂λi /λi∂ρlimTT→∞ ∂βi /βilim= −=βi,λiβiλi(39)(40)が得られる。Table 1 に U-235 熱中性子核分裂に対する Keepin の評価した遅発中性子データを示しているが、ai /λi （すなわち βi /λi ）は 2 群が最も大きいことが分かる。従って、T → ∞ の極限をとった場合であっても 2 群の感度が支配的となることが分かる。
一方、ペリオドが負の場合には、逆時間方程式で最も値が大きい（絶対値が小さい）解 s は
s > −λi でなければならないため、ペリオドがとり得る値の範囲は T = 1/s < −1/λi となる（−∞の反応度、すなわち、keﬀ → 0 となるような大きな負の反応度を与えた場合に T → −1/λi となる）。
ペリオドが負の場合における微係数 ∂ρ/(∂λi /λi )、∂ρ/(∂βi /βi ) の計算結果を Fig. 9 と Fig. 10に、それぞれ示す。反応度が負側に大きくなるほど（ペリオドが −1/λ1 に近づくほど）、最も崩壊定数が小さい（半減期が長い）1 群の感度が大きくなることが分かる。

90.00-0.01

Relative sensitivity

-0.02
-0.03
-0.04
-0.05
-0.06-0.07

i=1i=2i=3i=4i=5i=6-0.08-0.09-0.10024681012
Period [s]14161820Fig. 5: Relative sensitivity of reactivity to decay constant (positive period 1)

0.00
-0.01

Relative sensitivity

-0.02
-0.03
-0.04
-0.05
-0.06
-0.07

i=1
i=2
i=3
i=4
i=5
i=6-0.08
-0.09
-0.10
020

40

60

80

100 120
Period [s]140160180200Fig. 6: Relative sensitivity of reactivity to decay constant (positive period 2)

10

0.40
i=1
i=2
i=3
i=4
i=5
i=60.35

Relative sensitivity

0.30
0.25
0.20
0.15
0.10
0.05
0.00
02

46810
12
Period [s]14161820

Fig. 7: Relative sensitivity of reactivity to relative abundance (positive period 1)

0.40
i=1
i=2
i=3
i=4
i=5
i=60.35

Relative sensitivity

0.30
0.25
0.20
0.15
0.10
0.05
0.00
02040

60

80

100 120
Period [s]140160180200Fig. 8: Relative sensitivity of reactivity to relative abundance (positive period 2)

11

0.40
0.35

Relative sensitivity

0.30

i=1
i=2
i=3
i=4
i=5
i=60.25
0.20
0.15
0.10
0.05
0.00
-280 -260 -240 -220 -200 -180 -160 -140 -120 -100
Period [s]-80

Fig. 9: Relative sensitivity of reactivity to decay constant (negative period)

0.00
-0.05

Relative sensitivity

-0.10
-0.15
-0.20
-0.25
-0.30
-0.35

i=1i=2i=3i=4i=5i=6-0.40-280 -260 -240 -220 -200 -180 -160 -140 -120 -100
Period [s]-80

Fig. 10: Relative sensitivity of reactivity to relative abundance (negative period)12D

逆動特性法に関するメモ

このメモではペリオド法による反応度の測定を説明したが、同一の問題、すなわちステップ状に反応度が印加され、中性子数が漸近的に増加している領域（漸近領域）での反応度の測定に対しては、逆動特性法を用いても同様の値が算出される。以下、逆動特性法について説明する。
式 (2) を積分型に書き直すと以下のように書ける。
∫βi ∞Ci (t) =n(t − u) exp(−λi u)du
Λ 0(41)

これを式 (1) に代入すると以下の逆動特性法の式を得る。
∫ρ=β+6∞ n(t − u)
Λ dn ∑−λi βiexp(−λi u)dun(t) dt i=1n(t)0(42)

ここで、漸近領域の中性子数変動から逆動特性法を用いて得られる反応度を計算する。漸近領域
では中性子数は指数関数に従い変動するため、n(t) = N exp(st) と書ける。これを式 (42) に代入
し、Λ = l/k = l(1 − ρ) の関係式を用いると、逆時間方程式と同一の式を得ることができる。すな
わち、逆動特性法で計算される反応度はペリオド法によるものと一致する。

13