燃焼の基礎理論
名古屋大学 工学研究科 山本 章夫

1

講義の内容

燃焼の物理とモデル化
燃焼に伴う物理現象
 燃焼チェーン
 燃焼の定式化
 燃焼度、燃焼時間と絶対中性子束

燃焼の数値解法
数値解法の誤差要因
 数値解法

原子炉寿命の予測計算

Reactor Physics, Nagoya University, Akio YAMAMOTO

2

燃料の燃焼の重要性
原子炉の燃料は、核分裂を起こすことによりエネ
ルギーを発生する。
 原子炉物理では核分裂や中性子の吸収により、燃料の組成が変化していく過程を燃焼と呼ぶ。
 原子炉の特性は、燃料の組成に依存する。従って、燃焼に伴なって原子炉の特性は変化する。
 原子炉の運転可能期間の予測、運転期間を通じた原子炉の振る舞いの正確な予測の観点から燃焼解析は極めて重要

Reactor Physics, Nagoya University, Akio YAMAMOTO

3燃焼計算の重要性

ラマーシュ、「原子炉の初等理論」, p.730
Reactor Physics, Nagoya University, Akio YAMAMOTO

4燃料の燃焼に伴う物理現象
重要(1)核分裂性核種の減少
(2)親核種から核分裂性核種への転換

原子炉の増倍率に対してプラスの効果(3)核分裂生成物の発生

原子炉の増倍率に対してマイナスの効果

核分裂生成物は中性子の吸収物質であるため、原子炉の増倍率に対してマイナスの効果

上記(1)～(3)を合計すると、軽水炉では燃焼に伴い増倍率は減少する。

Reactor Physics, Nagoya University, Akio YAMAMOTO

5

燃焼チェーン

燃料内部では、燃焼に伴って(原子炉の運転に伴って)燃料内の核種が変化する。
核分裂
 中性子の吸収 α、β崩壊 (n, 2n)反応
 その他




核種の一連の変化を燃焼チェーンとして記述する。
例えば・・・
U238 (n,γ) U239 (β-崩壊) Np239 (β-崩壊) Pu239


Reactor Physics, Nagoya University, Akio YAMAMOTO

6

燃焼チェーン重要

燃焼チェーンは以下の核種に対するものに大別される
重核種(Th, U, Pu, Am, Cm)
 核分裂生成物(Xe, Sm, Pm,......)
 可燃性毒物(B, Gd, Er…)


Reactor Physics, Nagoya University, Akio YAMAMOTO

7

燃焼チェーン










炉内で生成・消滅する核種は数千にもなるが、炉心解析ではすべてを取り扱うわけではない。
炉心解析では、核的な特性に影響を与える核種を主として取り扱う。
多数の核種を取り扱うことにより、燃焼計算の精度は向上する
が、計算時間は長くなる。
従来は、計算時間の観点から簡略な燃焼チェーン(重核種20,FP30核種程度)が用いられていたが、最近は詳細な燃焼チェーン(重核種30, FP200核種程度)が用いられる場合が多い。
簡略な燃焼チェーンでは、主要なFPを個別に扱い、残りのFPはグループ化し、疑似FP (pseudo FP)として取り扱う。疑似FP
は、詳細な燃焼チェーンとの比較により作成するが、収率や断面積の作成にノウハウが必要。

Reactor Physics, Nagoya University, Akio YAMAMOTO

8

燃焼チェーン

奥村啓介、MVP-BURN利用マニュアル、日本原子力研究所、2004

MVP-BURN重核種チェーン
Reactor Physics, Nagoya University, Akio YAMAMOTO

9燃焼チェーン

MVP-BURN重核種チェーン 奥村啓介、MVP-BURN利用マニュアル、日本原子力研究所、2004
Reactor Physics, Nagoya University, Akio YAMAMOTO

10

燃焼チェーン

MVP-BURN簡略チェーン奥村啓介、MVP-BURN利用マニュアル、日本原子力研究所、2004
Reactor Physics, Nagoya University, Akio YAMAMOTO

11

燃焼チェーン

MVP-BURN標準チェーン

Reactor Physics, Nagoya University, Akio YAMAMOTO

奥村啓介、MVP-BURN利用マニュアル、日本原子力研究所、2004

12燃焼チェーン

MVP-BURN汎用チェーン
Reactor Physics, Nagoya University, Akio YAMAMOTO

奥村啓介、MVP-BURN利用マニュアル、日本原子力研究所、2004

13

燃焼チェーン
MVP-BURN開発者用チェーン

奥村啓介、MVP-BURN利用マニュアル、日本原子力研究所、2004

Reactor Physics, Nagoya University, Akio YAMAMOTO

14

燃焼チェーン


可燃性毒物の燃焼チェーン

可燃性毒物：原子炉の寿命初期の反応度を抑える(小さくする)ために、炉心内に配置される中性
子吸収物質。
 燃焼とともにその量(原子数密度)が減るため、吸収効果が小さくなる。すなわち、反応度を抑える効果も小さくなる。
Reactor Physics, Nagoya University, Akio YAMAMOTO

15

燃焼チェーンの核種数と精度


取り扱う核種数とメモリ量の関係
PWR, 4ループ, 全炉心 燃料棒単位(詳細メッシュ計算)
 軸方向60メッシュ
 取り扱い核種数300
300[核種/燃焼領域]×300[燃焼領域/集合体径方向断面]×193[集合体/炉心]×60[軸方向メッシュ]＝1×109(byte)

約1GBのメモリが必要
 目的によって適切な核種数の燃焼チェーンを使用することが望ましい

Reactor Physics, Nagoya University, Akio YAMAMOTO

16

燃焼チェーンの選択


解析目的
核特性解析
 中性子/γ線源強度・発熱量評価
 PIE解析MVP-BURNコードで用意されている燃焼チェーン

チェーン名
標準
汎用
開発者用
簡易

重核種数
21 (234U-246Cm)
21 (234U-246Cm)
28 (232Th-246Cm)
19(235U-245Cm)

FP 核種数
49
103
193
30

Reactor Physics, Nagoya University, Akio YAMAMOTO

疑似 FP 核種数
1
1
0
4

備考/用途
UO2/MOX 用
PIE 解析など
チェーンモデルの検証
UO2/MOX 用

17

FP収率(独立収率と積算収率)
FP収率はJNDC decay data fileなどに格納されている。(1000核種以上)
 この収率は個々のFP核種に対する「独立収率」
である。
 一方、チェーンを簡略化した場合、「積算収率」を用いなければならない。
 このため、FP収率は燃焼チェーンごとに異なり、値を設定するためには大変な労力を要する。


Reactor Physics, Nagoya University, Akio YAMAMOTO

18

FP収率(独立収率と積算収率)
MVP-BURNの燃焼チェーンにおけるRh105とPd105の収率
チェーン名
開発者
標準
簡略

105

Rh
1.70062E-08
9.63844E-03---

MVP-BURN簡略チェーン
Reactor Physics, Nagoya University, Akio YAMAMOTO

105

Pd
9.99776E-13
9.99776E-13
9.63844E-03

MVP-BURN開発者用チェーン
奥村啓介、MVP-BURN利用マニュアル、日本原子力研究所、2004

19

疑似FPの評価
燃焼チェーンにおけるFP核種の取り扱い数を少なくする場合、疑似FPが必要となる。
 疑似FPは、複数のFP核種をひとまとめにしたものである。
 疑似FPの収率や断面積は、燃焼チェーンの設計に依存する。
 疑似FPの収率や断面積は、詳細な燃焼チェーンを用いた計算結果と比較することで設定する。
Reactor Physics, Nagoya University, Akio YAMAMOTO

20

原子炉の燃焼計算
燃焼に伴う原子炉特性の変化の解析を燃焼計
算と呼ぶ。
 燃焼計算は、原子炉の設計計算において、中心的な役割を果たす。
 燃焼計算では、燃料を構成する核種の原子数密度がどのように変化していくかを精度良く計算する必要がある。
 では、核種の原子数密度の変化を支配する方程式は・・・？


Reactor Physics, Nagoya University, Akio YAMAMOTO

21

燃料の燃焼計算例
1.E-01
0

10

20

30

40

50

60

70

80

1.E-02

原子数密度(×1e24)

1.E-03
1.E-04
U235
U236
U238

1.E-05
1.E-06
1.E-07
1.E-08
1.E-09燃焼度(GWd/t)

Reactor Physics, Nagoya University, Akio YAMAMOTO

22単一核種の燃焼
重要

単一核種の個数密度は燃焼に伴って以下の方
程式(燃焼方程式)に従って変化する。
235U微視的

25
吸収断面積
(
)
dN
t
25
U235の例= − N (t ) × σ a × φ (t )dt
235U原子

数密度



中性子束

つまり、「核種が中性子を吸収した分だけ減少」
する

[復習]原子数密度×微視的断面積×中性子束=反応率

Reactor Physics, Nagoya University, Akio YAMAMOTO

23単一核種の燃焼

中性子束を一定とすると・・・

この微分方程式は簡単に解けて

dN 25 (t )= − N 25 (t ) × σ a × φdt

N 25 (t ) = N 25 (0) exp(−σ aφt )235U原子

数密度の初期値

Reactor Physics, Nagoya University, Akio YAMAMOTO

24単一核種の燃焼

例えば、典型的な軽水炉の例として
σ a = 1000 [barn]
φ = 1013[1/cm2/sec]の場合： N 25 (t ) = N 25 (0) exp(−10 −8 t ) t:[sec]1.2

U235原子数(相対値)

1.0
0.8
0.6
0.4
0.2
0.0
0

5

10

15
20
運転日数(月数)

Reactor Physics, Nagoya University, Akio YAMAMOTO

25

30

35

25

単一核種の燃焼

中性子束を一定と見なせる期間は余り長くないので、実際の計算では、時間間隔Δt毎に中性子束を再計算する必要がある。従って、原子数密度もΔt後の値を求める必要がある。
N 25 (t + ∆t ) = N 25 (t ) exp(−σ aφ∆t )Reactor Physics, Nagoya University, Akio YAMAMOTO

26

複数の核種が関係する燃焼

以下の核種の燃焼を考える：複数の核種が関係する「燃焼チェーン」
U238+n→U239 (or fission)→Np239→Pu239
 Pu239+n→Pu240 or fission
 Pu240+n→Pu241 or fission
 Pu241+n→Pu242 or fission

軽水炉で重要な転換反応

Reactor Physics, Nagoya University, Akio YAMAMOTO

27

複数の核種が関係する燃焼

U238

各核種の原子数密度に対する微分方程式(1)dN 28 (t )= − N 28 (t ) × σ a28 × φ U235と同じ考え方
dt

49
dN
(t )
Pu239
= − N 49 (t ) × σ a49 × φ + N 28 (t ) × σ c28 × φ
dt

吸収による減少分

U238の捕獲反応による増加分

U238+n→U239 (or fission)→Np239→Pu239
Pu239+n→Pu240 or fission

Reactor Physics, Nagoya University, Akio YAMAMOTO

28

複数の核種が関係する燃焼

各核種の原子数密度に対する微分方程式(2)Pu240

dN 40 (t )= − N 40 (t ) × σ a40 × φ + N 49 (t ) × σ c49 × φ
dt

Pu241

dN 41 (t )= − N 41 (t ) × σ a41 × φ + N 40 (t ) × σ c40 × φ
dt

Reactor Physics, Nagoya University, Akio YAMAMOTO

29

複数の核種が関係する燃焼

解析解は非常に煩雑になり、一般化は困難

Reactor Physics, Nagoya University, Akio YAMAMOTO

30

複数の核種が関係する燃焼
四つの連立微分方程式を解けばU238, Pu239,
Pu240, Pu241の原子数密度の変化が求まる。
 さらに複雑な燃焼チェーンに対しても同様の考え
で原子数密度の変化を記述できる。
 実際の原子炉設計では、核種の原子数密度は燃焼チェーンを記述する連立微分方程式を以下の方法で解いている。


解析的方法：Bateman法
 数値解析法：指数行列法・ルンゲ・クッタ法など


Reactor Physics, Nagoya University, Akio YAMAMOTO

31

燃焼方程式の一般化


核種の生成・消滅は、以下の式により記述できる。
dN i= 生成率 − 変換率 − 崩壊率dt

一般的に書くと以下のようになる

dN i= ∑ γ jiσ f , j N jφ +σ c ,i −1 N i −1φ + ∑ λk N k − σ a ,i N iφ − λi N i
dt
j
k
核種jの核
分裂によ
る生成

Reactor Physics, Nagoya University, Akio YAMAMOTO

核種i-1の
捕獲によ
る生成

核種kの
崩壊によ
る生成

核種iの 核種iの吸収によ 崩壊る変換

32

燃焼方程式の一般化


前出の燃焼方程式は連立の一階常微分方程式
である。核種の原子数密度をベクトルと見ると、一般的に以下のように書ける。

dN i= ∑ γ jiσ f , j N jφ +σ c ,i −1 N i −1φ + ∑ λk N k − σ a ,i N iφ − λi N idtjk
核種の原子数密度をベクトルとみて・・



dN
= AN
dt
Reactor Physics, Nagoya University, Akio YAMAMOTO

33

燃焼方程式の一般化
「燃焼マトリックス」Aは、以下の構造を持つ。
1

1

A

150

100

50

221

核分裂によるFP発生N

1FP（原子量小)

50

50

100

100

150

150

221

…

×

1

50

100
1

吸収・崩壊による消滅(対角要素)
捕獲による生成(対角要素より一段左)*)β崩壊による生成

1

150

10

20

221

FP（原子量大)
重核種（原子量小)…
重核種（原子量大)

221
28

1

10

10

20

20

α崩壊による生成
Reactor Physics, Nagoya University,
Akio YAMAMOTO
28

1

10

20

28

28

*)核種ベクトルの並び方によっては、二段以上左に現れる場合がある。

34

燃焼方程式の一般解
dN 25 (t )= −σ a × φ × N 25 (t )dt

N (t ) = exp(−σ aφt ) N (0)
25

25

単一核種

Reactor Physics, Nagoya University, Akio YAMAMOTO



dN
= AN
dt


N (t ) = exp(At )N (0)
複数核種

35

燃焼方程式の一般解
exp(At)は指数行列と呼ばれる。2x2など、小さい行列に対しては、解析解を計算することも可能。
 一般的な燃焼行列(行列Aのサイズが数十～数百)では、何らかの仮定を用いないと解析解を求めることは困難。従って、数値解法を用いることが一般的。


Reactor Physics, Nagoya University, Akio YAMAMOTO

36

燃焼度の定義
燃焼度は初期の重金属量1トンあたりに発生した積算熱量で与えられる。
 一般的には、MWd/t, GWd/tが使用される。
 1 MWd/tは、重金属1 tonあたり、1MWの出力を1day継続して出した熱量に相当する。
 重金属量は、燃焼前の値を用いる。


Reactor Physics, Nagoya University, Akio YAMAMOTO

37

燃焼度の定義


ウラン235 1gが全て核分裂したときの燃焼度は?
1[g]/235.04[g/mol]×6.02×1023[個/mol]
×200[MeV]×1.6×10-13[J/MeV]=8.2×1010[J]1gあたりこれだけの熱量を発生するので、8.2×1010[J/g]
一方、
1[MWd]=106[W・d]×3600×24[s/d]=8.6×1010[W・s]=8.6×1010[J/s・s]=8.6×1010[J]
これより、
1[J]=1/8.6×1010[MWd]=1.16×10-11[MWd]以上より、8.2×1010[J/g]=8.2×1010[J/g]×1.16×10-11[MWd/J]=0.95[MWd/g]= 0.95[MWd/g]×106[g/t]=9.5×105 [MWd/t]

Reactor Physics, Nagoya University, Akio YAMAMOTO

38

燃焼度の定義


ウラン235 1tが全て核分裂したときの燃焼度は?

Reactor Physics, Nagoya University, Akio YAMAMOTO

39

燃焼時間の算出
炉心の実効増倍率や相対出力分布のみを求める計算では、中性子束の絶対値は不要である。
 これに対し、燃焼計算では、中性子束の絶対値と燃焼時間が不可欠である。
 中性子束の絶対値と燃焼時間は以下の手順で算出する。


核分裂による発生熱量を体系全体で積算し、炉心熱出力(入力値)になるように中性子束を規格化する。
 燃焼ステップをMWd/t単位で入力する場合、日数(day)=MWd/t÷炉心熱出力(MW)×重金属量(t)から、燃焼日数を算出。
 燃焼ステップをday単位で入力する場合、MWd/t=炉心熱出力(MW)×日数(day)÷重金属量(t)から、MWd/tを算出。
 重金属量は、燃焼前(寿命初期)の値を使用する。


Reactor Physics, Nagoya University, Akio YAMAMOTO

40

燃焼時間の算出


中性子束の規格化と燃焼時間の計算は、非常に間違えやすいので、十分にチェックする必要がある。具体的には：入力する出力の単位(W/kW/MW) 対象とする体系(炉心全体、集合体一体、高さ1cmあたりの集合体
一体、など)
 重金属量(g/kg/t/st) 炉心の体積などは要注意

コード内で計算された燃焼日数(燃焼度で入力した場合)、
燃焼度(時間で入力した場合)などをチェックすることで、かなりの間違いは防ぐことができる。

Reactor Physics, Nagoya University, Akio YAMAMOTO

41

燃焼計算の誤差要因

反応率(生成率、変換率、崩壊率)の誤差
→反応率の計算精度に依存

燃焼チェーンの簡略化に伴う誤差
 燃焼方程式の数値計算法に伴う誤差

例：燃焼チェーンの線形化に伴う誤差

燃焼方程式の時間離散化に伴う誤差

燃焼方程式の時間離散化誤差
数値計算の丸め誤差
反応率の離散化誤差

熱出力の規格化に伴う誤差→数値解法に依存→実機の解析では重要

初期組成の誤差

照射後試験(PIE)の解析では重要

Reactor Physics, Nagoya University, Akio YAMAMOTO

42

燃焼方程式の数値解法


燃焼方程式を以下のように記述する。


dN
= AN
dt


N : 原子数密度ベクトル
A : 生成、変換、崩壊率を表す行列

反応率(σφ)は燃焼ステップごとに与えられるとする。燃焼ステップ間で反応率は一定と仮定。

Reactor Physics, Nagoya University, Akio YAMAMOTO

43

主要な数値解法
オイラー法
 ルンゲ・クッタ法
 指数行列法
 行列分解法 Bateman法
 Padé近似法
 Krylov部分空間法
 Chebychev Rational Approximation法(CRAM)
 Mini-max polynomial approximation (MMPA)

Reactor Physics, Nagoya University, Akio YAMAMOTO

44

オイラー法
有限差分法に基づく、もっとも単純な数値解法
 単純ではあるが、精度が良くないため、燃焼計
算に用いられることは少ない。

 l +1  l
ldN N − N
≈
= AN
∆t
dt
 l +1  ll
N = N + AN ∆t l +1  l
 l +1dN N − N
≈
= ΑN
dt
∆t
 l +1  l(I − Α∆t )N = N
Reactor Physics, Nagoya University, Akio YAMAMOTO

陽的解法

陰的解法
45

オイラー法燃焼(タイム)ステップを細かくすると、方程式の離散化誤差は減少する。
 一方、燃焼ステップ数が多くなるため、数値計算の丸め誤差が蓄積する。
 燃焼ステップ数が多くなると計算時間は増加する。


上記の観点から、適正な燃焼ステップを決める必要がある。

Reactor Physics, Nagoya University, Akio YAMAMOTO

46

ルンゲ・クッタ法
常微分方程式の数値解法としてよく用いられている手法。
 2次、4次のルンゲ・クッタ法が有名
 2次のルンゲ・クッタ法 dN l  l
d1 =t , N ∆tdt
 dN  l ∆t  l d  t + , N + 1 ∆t
d2 =
dt 
2
2 
 l +1  l 
N = N + d2()

Reactor Physics, Nagoya University, Akio YAMAMOTO

47

ルンゲ・クッタ法
4次のルンゲ・クッタ法
 計算精度が比較的良好で、計算時間が短いた
め、燃焼計算に用いられることがある dN l  l
d1 =t , N ∆tdt
 dN  l ∆t  l d  t + , N + 1  ∆t
d2 =
dt 
2
2 


 dN  l ∆t  l d  t + , N + 2  ∆t
d3 =
dt 
2
2 

 dN ll 
d4 =t + ∆t , N + d 3 ∆t
dt


 
 l +1  l d1 + 2d 2 + 2d 3 + d 4
N =N +
6

(

(

)

)

Reactor Physics, Nagoya University, Akio YAMAMOTO

48

指数行列法


指数行列法では、燃焼方程式の解析解を以下のように
与える。



N (t + ∆t ) = exp(A∆t )N (t )

指数行列は、Taylor級数で与えることができる。

(
A∆t ) (A∆t )exp(A∆t ) = I + (A∆t ) +
+
+ ...
2

2!



3

3!

一見、理想的な解法に見えるが・・・

Reactor Physics, Nagoya University, Akio YAMAMOTO

49

指数行列法

簡単のため、exp(-100)をTaylor級数展開で計算することを考える。
2
3(− 100 ) (− 100 )exp(− 100 ) = 1 + (− 100 ) +
+
+ ...

2!

3!

exp(-100)～4×10-44である。
 上記のTaylor級数は100項目近辺で最大となり、(-100)100/100!～1×1042 である。
 これより、exp(-100)を5桁の精度で計算するた
めには、44+42+5=92桁以上の精度が必要。
 このような計算は非常に困難。
Reactor Physics, Nagoya University, Akio YAMAMOTO

50

指数行列法以下の級数を数値計算で使用するためには、行
列AΔtのノルムが十分小さくなければならない。
32(A∆t ) (A∆t )exp(A∆t ) = I + (A∆t ) +
+
+ ...

2!

3!

Δtが大きいと、AΔtのノルムは大きくなる。
 Aに含まれている核種の最小の半減期の10程度
がΔtの目安になる。→短半減期核種が含まれている場合、詳細なタイムステップが必要Reactor Physics, Nagoya University, Akio YAMAMOTO

51

指数行列法
例えばMVP-BURNの開発者用チェーンの場合、Rh106の半減期が最も短い。(29.8sec)
 この場合、Δtは約300sec以下にする必要がある。
 一方、通常の燃料の燃焼期間は数年。(~108sec)
 これより、燃焼ステップを108/300~3×105と非常に多数とる必要が生じる→計算時間の点から困難 従って、短半減期核種は平衡解、解析解(Bateman法)などの方法で計算する必要がある。
Reactor Physics, Nagoya University, Akio YAMAMOTO52

指数行列法

指数行列の単純なTaylor展開ではΔtに制限が生じるため、以下のような変形を行った上、級数を計算する場合もある。適切なmを選ぶことにより、級数の収束を早めることができる。（Scaling and Squaring)  A∆t 
exp(− A∆t ) = exp −  m m

この方法でも、単寿命の核種が存在する場合は取り扱いが難しい。
 行列のべき乗の操作を行うため、丸め誤差が発生しやすい。
Reactor Physics, Nagoya University, Akio YAMAMOTO

53

指数行列法とルンゲ・クッタ法の関係


指数行列法の展開次数とオイラー法、ルンゲ・クッ
タ法の関係
2項目(一次の項)まで考慮： オイラー法
3項目(二次の項)まで考慮： ルンゲ・クッタ法(2次)
5項目(四次の項)まで考慮： ルンゲ・クッタ法(4次)(A∆t ) (A∆t )exp(A∆t ) = I + (A∆t ) +
+
+ ...
2

2!

Reactor Physics, Nagoya University, Akio YAMAMOTO

3

3!

54

行列分解法
前述の指数行列法では、どのような燃焼チェー
ンでも取り扱うことができる。言い換えると、行列
Aの形に制限はない。
 しかし、指数行列の計算に時間を要する。
 そこで、Aを対角化することで計算時間の短縮を図る→燃焼チェーンの「線形化」
dN i= ∑ γ jiσ f , j N jφ +σ c ,i −1 N i −1φ + λi −1 N i −1 − σ a ,i N iφ − λi N idtj= Fi + (σ c ,i −1φ + λi −1 )N i −1 − (σ a ,iφ + λi )N i

Reactor Physics, Nagoya University, Akio YAMAMOTO55

行列分解法

線形化した燃焼チェーンは、一般的に以下のように書くことができる。
 dN= − AN + F
dt

0
0

0 
 σ a ,1+λ1
−σ −λ σ + λ 


2
a,2

 c ,1 1



 
A= 0



0

0
σ
λ
+
a , n −1
n −1


 0

0 −σ c ,n −1−λn −1 σ a ,n + λn 




解は

[

]

 l +1l−1−1N = exp(− A∆t ) N − A F + A F

Reactor Physics, Nagoya University, Akio YAMAMOTO

56

行列分解法

行列Aを対角化する行列Lを求める。ここでDは対角化行列。

L−1AL = D



すると



となる。この変形には、以下の関係式を使用している。

()

exp(− A∆t ) = exp − LDL−1∆t = L exp(− D∆t )L−1

2
3(
− A∆t ) (− A∆t )exp(− A∆t ) = I + (− A∆t ) +
+
+ ...

2!

3!

チェーンの「線形化」を行うと、行列の対角化が容易になる。
 Dは対角行列なので、exp(-DΔt)は容易に計算できる。


Reactor Physics, Nagoya University, Akio YAMAMOTO

57

行列分解法
線形化チェーンを用いた方法は、短い時間で精度よい計算が可能である。
 しかし、以下のような問題点もある。


α崩壊、β崩壊、(n, 2n)など、複雑な分岐および合流が存在するチェーンをそのまま扱うことができない。
 すべての燃焼チェーンを線形化する必要がある。(分岐するチェーンの扱いが煩雑) 変換率(=吸収率+崩壊率)が全く同じ核種が二つあると数値的に(ゼロ割が発生するため)計算できない。
 循環形式のチェーンをそのまま扱うことが出来ない。
Reactor Physics, Nagoya University, Akio YAMAMOTO

58

行列分解法

分岐するチェーンの取り扱いの例


パスごとに線形化したチェーンを定義する必要がある

Pu241→Pu242→Am243→Cm244
Pu241→Am241→Am242m→Am243→Cm244
Pu241→Am241→Am242g→Pu242→Am243→Cm244
Pu241→Am241→Am242g→Cm242→Cm243→Cm244
Pu241→Am241→Am242g→Am243→Cm244
Reactor Physics, Nagoya University, Akio YAMAMOTO

59

Bateman法


線形な一階の常微分方程式に対する解析解を直接用いる方法。行列分解法で用いた「線形化された」燃焼方程式に対する解析解は：N i (t + ∆t ) = N i (t ) exp(−d i ∆t )
 i −1
i −1
exp(
)
exp(
)
d
t
d
t
−
∆
−
−
∆
a
j
i
n +1, n 
a j +1, j ∏+ ∑ N k (t ) ∑
 j =kdi − d j
k =1
n=k d n − d j


n≠ j
i −1

ここで、aijは、核種iからjへの変換率[行列Aの(i,j)の要素]、
diは核種iの消滅率[=aii]。
 解析解を用いるので、計算精度は高い。
 短半減期の核種を含む場合に有効。
 行列分解法と同じく、燃焼チェーンの「線形化」が必要。
 SRAC, MVP-BURNなどで使用されている。
 行列分解法の「解析解」版
Reactor Physics, Nagoya University, Akio YAMAMOTO

60

Padé近似法


𝒆𝒆−𝒙𝒙

𝟏𝟏
𝟏𝟏 + 𝒙𝒙𝟏𝟏 − 𝒙𝒙

指数行列の近似として、Taylor展開ではなく、
Padé近似を使用する。
N pq ( A∆t )exp( A∆t ) ≅D pq ( A∆t )2
N 22 ( A∆t ) 12 + 6 A∆t + ( A∆t )e.g .
=D22 ( A∆t ) 12 − 6 A∆t + ( A∆t )2Padé近似は、一般にTaylor級数より収束が早い
 MathematicaやMATLABなどで用いられている指数行列の一般的な数値解法


Reactor Physics, Nagoya University, Akio YAMAMOTO

61

Krylov部分空間法


指数行列の第m項までの級数展開の一般形は






1
2
m −1
N (t + ∆t ) ≅ c0 N (t ) + c1 (A∆t ) N (t ) + c 2 (A∆t ) × N (t ) + ⋅ ⋅ ⋅ + c m −1 (A∆t ) × N (t )

1
c
=
k

k!

のとき、Taylor展開になる Krylov部分空間法では、ckをN (t + ∆t )の最小自乗残差が最も小さくなるように選ぶ→行列Aに応じた最適な展開係数ckを求める
 計算理論はやや難解

Reactor Physics, Nagoya University, Akio YAMAMOTO

62

Krylov部分空間法
4.1wt% UO2ピンセル計算時の核種個数密度の最大誤差(1MWd/t~2300sec時)、MVP-BURNの開発者用チェーンを使用
Absolute maximum percent error of nuclide
number density

1.E+23
1.E+21
1.E+191.E+17
1.E+15
1.E+13
1.E+11
1.E+09
1.E+07
1.E+05
1.E+03
1.E+01
1.E-01
1.E-03
1.E-05Taylor expansion
Krylov subspace

0

5

10

15

20

Order of expansion/Krylov dimension
Reactor Physics, Nagoya University, Akio YAMAMOTO

25

30

63

CRAM


Chebyshev rational approximationにおいては、指数行列は以下のように近似される。
α jI
exp[ A∆t ] = α 0 I + ∑j =1 A∆t − θ j Ik k / 2 α jI 
= α 0 I + 2 Re ∑

A
t
I
θ
∆
−
j
=
1

j 




従って、



N (t + ∆t ) = exp[ A∆t ]N (t )  k / 2 α jI
N (t )
= α 0 N (t ) + 2 Re ∑ j =1 A∆t − θ j Ik /2   −1= α 0 N (t ) + 2 Re ∑ x j (t ) x j (t ) = (A∆t − θ j I ) α j N (t )
 j =1

Reactor Physics, Nagoya University, Akio YAMAMOTO

64

CRAM


係数の一例

M. Pusa, “Rational Approximations to the Matrix Exponential in Burnup calculations,”
Nucl. Sci. Eng.. 169, 155 (2011).

Reactor Physics, Nagoya University, Akio YAMAMOTO

65

CRAM


計算精度


極短半減期の核種に対しても良い精度
M. Pusa, “Rational Approximations to the Matrix Exponential in Burnup calculations,”
Nucl. Sci. Eng.. 169, 155 (2011).



計算時間


比較的短い

Reactor Physics, Nagoya University, Akio YAMAMOTO

66

MMPA
燃焼方程式の解析解を以下のように行列の多項式で近似する。
 係数は、最大誤差が最も小さくなるように(minimax)決定する𝑛𝑛

𝒏𝒏 𝑡𝑡 = 𝑎𝑎0 𝒏𝒏𝟎𝟎 + � 𝑎𝑎𝑖𝑖 𝐀𝐀𝑡𝑡 + 𝑐𝑐𝐈𝐈 𝐀𝐀𝑡𝑡 − 𝑐𝑐𝐈𝐈 −1 𝑖𝑖 𝒏𝒏𝟎𝟎𝑖𝑖=1

展開次数が同じであれば、CRAMと同等の精度であるが、計算の安定性及び計算時間の観点か
ら有利

Reactor Physics, Nagoya University, Akio YAMAMOTO

67

各手法の比較

ルンゲ・クッタ法








複雑なチェーンをそのまま取り扱い可能、4次は計算精度よい、計算時間短い短半減期核種が存在する場合は、タイムステップを細かくする必要あり。場合によっては、計算が困難。
短半減期核種は平衡解やBateman法などで評価する必要あり。

指数行列法(Taylor展開)







複雑なチェーンをそのまま取り扱い可能、計算精度よい
短半減期核種が存在する場合は、タイムステップを細かくする必要あり。場合によっては、計算が困難。
短半減期核種は平衡解やBateman法などで評価する必要あり。
計算時間長いORIGEN2コードなど

Reactor Physics, Nagoya University, Akio YAMAMOTO

68

各手法の比較

行列分解法(線形化チェーン)

計算精度よい複雑な燃焼チェーンの取り扱いは煩雑

Bateman法(線形化チェーン)

計算精度よい複雑な燃焼チェーンの取り扱いは煩雑
DCHAIN, SRAC, MVP-BURN

Reactor Physics, Nagoya University, Akio YAMAMOTO

69

各手法の比較Padé近似法


指数行列法(Taylor級数展開)より収束早い
燃焼計算への適用実績少ない

Krylov部分空間法


多数の核種が関係する非常に複雑な燃焼チェーンを取り扱い可能
短半減期核種を直接取り扱い可能
計算時間短いMARBLE等

Reactor Physics, Nagoya University, Akio YAMAMOTO

70

各手法の比較


CRAM







多数の核種が関係する非常に複雑な燃焼チェーンを取り扱い可能
極短半減期核種を直接取り扱い可能
計算時間短いSERPENT2コードなど

MMPA





多数の核種が関係する非常に複雑な燃焼チェーンを取り扱い可能
極短半減期核種を直接取り扱い可能
計算時間短い数値安定性高い

Reactor Physics, Nagoya University, Akio YAMAMOTO

71

予測子・修正子法
燃焼計算の誤差要因として、時間の離散化(タイプステップ)のほかに、反応率の離散化がある。
 燃焼方程式を解く際、タイムステップ間では反応率を一定にする。これが誤差要因となる。特に、Gd燃料のように燃焼とともに中性子スペクトルが大きく変化する場合は、誤差が顕著に表れる。
σ (T + ∆t )φ (T + ∆t )σ (T )φ (T )T

T + ∆t燃焼に伴う吸収率の変化の例
Reactor Physics, Nagoya University, Akio YAMAMOTO

72

予測子・修正子法
反応率の誤差はタイムステップを小さくすることで低減することが可能である。
 しかし、実機の解析においては、反応率の計算はもっとも時間を要するため、単純にタイムステップを増加させるのは得策ではない。(実効断面積の作成→スペクトル計
算に時間がかかる)

そのため、予測子・修正子法が一般的に用いられる。

Reactor Physics, Nagoya University, Akio YAMAMOTO73

予測子・修正子法予測子・修正子法では、以下の手順で計算を行う。
初期状態の原子数密度を元に反応率を計算する。
 この反応率で燃焼計算を行い、1ステップだけ燃焼を進める。(predictor計算)
 1ステップ燃焼を進めた時点の原子数密度を元に反応率を計算する。
 この反応率で初期状態からの燃焼計算をもう一度行う。(corrector計算)
 Predictor計算とcorrector計算の原子数密度を平均し、最終的な原子数密度の計算値とする。


Reactor Physics, Nagoya University, Akio YAMAMOTO

74

予測子・修正子法
予測子・修正子法では、反応率の計算回数は2倍になるが、計算精度はタイムステップ(反応率
計算の間隔)を1/10程度にしたものに相当する。
従って、トータルの計算時間は短くなる。
 なお、予測子・修正子法にはいくつかのバリエー
ションがあり、predictorとcorrectorの反応率の平均値を用いて燃焼計算をやり直す方法もある。


Reactor Physics, Nagoya University, Akio YAMAMOTO

75

Projected PC法(PPC法)


燃焼計算の主要な数値計算誤差(1)時間離散化(2)燃焼ステップ内における反応率一定の仮定



特にGd燃料においては、燃焼に伴い反応率が急激に変化するため、(2)が大きな誤差要因になる
そのため、予測子修正子法(PC法)の使用が不可欠
 しかし、PC法を使用しても、燃焼ステップは0.2GWd/t
程度が限界
 結果的に、集合体燃焼計算には、200-300MWd/t程
度の燃焼ステップが必要となる。
76

Gd燃料の反応率の誤差
Gd同位体 微視的反応率

Step
n

Step
n+1

Burnup77

従来のPC法
Gd同位体微視的反応率

Step
n

Step
n+1

Burnup78

PC法の離散化誤差
Gd同位体微視的
過小評価
反応率(燃焼の遅れ)

Step
n

Step
n+1

Burnup79

Effective micro reaction rate of Gd155

PPC法のコンセプト
5.4E-02
5.2E-02
5.0E-02
4.8E-02
4.6E-02
4.4E-02
4.2E-02
4.0E-02
2.2E-04

Step
n+1

Step
n
2.4E-04

Number density of Gd155[1/barn/cm]

2.6E-0480

PPC法の計算フロー(1)従来のPC法の計算を実施(2)「実効的な」微視的反応率を計算
N np+1 = N n exp(− R p )N nc+1 = N n exp(− R c )R p = − ln( N np+1 / N n )R c = − ln( N nc+1 / N n )R p − Rc
p
c
R=
N
N
+
R
(
−
)
n +1
p
N n − N n +1

Effective micro reaction rate of Gd155(3)原子数密度と実効的な微視的反応率の関係を評価
5.4E-02
5.2E-02
5.0E-02
4.8E-02
4.6E-02
4.4E-02
4.2E-02

R cp
R

4.0E-02
2.2E-04

N np+1 N n
2.4E-04

Number density of Gd155[1/barn/cm]

81

2.6E-04

PPC法の計算フロー(4)Corrector計算のために、より正確な反応率を外挿して計算↓ N 従来のPC法で評価
p
c
R −R
p
c
p
c
c ,mod
N
N
N
R
R
+
−
+
=
[(
)
/
2
]
n +1
n +1
n +1
p
N n − N n+1(5)Correctorステップの原子数密度を再計算
pn +1

c , mod
N nc+,mod
=
N
−
R
exp(
)
1
n

(6)Predictorステップとcorrectorステップの原子数
密度を平均し、次ステップの最終的な原子数密
度とする

(

)

p
N n +1 = N nc+,mod+N1
n +1 / 2

82

100

1.30

0

1.25
PPC, 0.25GWd/t
PPC, 0.5GWd/t
PPC, 1.0GWd/t
PC, 0.25GWd/t
PC, 0.5GWd/t
PC, 1.0GWd/t
Reference k-infinity

-100
-200
-300

1.20
1.15
1.10
1.05
1.00

-400

k-infinity

Difference of k-infinity(pcm)

PPC法の精度 (6wt%Gd2O3)

0.95

-500

0.90

-600

0.85

-700

0.80
0

10

20

30

40

Assembly burnup (GWd/t)

50

60
83

サブステップ
燃焼方程式を解く際には、タイムステップを細分化することで計算精度を向上することができる。
 しかし、反応率の再計算には時間がかかる。
 そこで、反応率の再計算は行わずに、(燃焼方程式を解く)タイムステップのみを細分化する手法をサブステップ法という。


出力

Predictor計算

Corrector計算

Pn
Pn −1

t m t m +1

tn
Reactor Physics, Nagoya University, Akio YAMAMOTO

Pn +1

t n +1

時間

84

サブステップ法

サブステップ法を用いることは、以下の点からメリットがある。


時間に関する離散化誤差の低減



燃焼ステップ内での出力の規格化

出力

出力

入力

入力
計算

計算

tn

t n +1

時間

ステップ内で再規格化なし
Reactor Physics, Nagoya University, Akio YAMAMOTO

tn

t n +1

時間

ステップ内で再規格化あり
85

冷却計算


崩壊による核種の変換は、炉心特性に大きな影響を及ぼすので、実機の解析ではこの影響(冷却効果)を考慮する必要がある。軽水炉の解析では以下の崩壊が重要

I135→Xe135(半減期6.7hr)

Np239→Pu239(半減期2.4d)

再使用燃料の反応度に影響

Pu241→Am241(半減期14.4y)


再使用燃料の反応度に影響

Eu155→Gd155(半減期4.7y)

通常の起動試験時の炉心特性に影響

Pm148→Sm148(半減期2.6y)

通常の起動試験時の炉心特性に影響

Pm148m→Sm148(半減期41d)

スクラム後の再起動時の臨界点に影響（Pm149→Sm149の壊変と反応度的にキャンセルする傾向)

Pm149→Sm149(半減期54hr)

スクラム後の再起動時の臨界点に影響

再使用燃料の反応度に影響(特にMOX燃料で重要)冷却計算は、中性子束レベルを非常に小さくした燃焼計算を行うことで模擬することができる。

Reactor Physics, Nagoya University, Akio YAMAMOTO

86

原子炉の燃焼計算
重要

燃焼計算は以下の手順で行う。




ある時点での原子炉内の巨視的断面積の計算
Σ = Nσ いろいろな核種について合計ある時点での原子炉内の中性子束分布の計算(拡散理論)

1− D∇ φ + Σ aφ = νΣ f φk 原子炉出力から中性子束の絶対値を計算
2

P=∫

Core


 
E R Σ f ( r )φ ( r )dr

燃料を構成する核種の反応率(核分裂、吸収・・・)を

計算 R = φ (r )σ a
 ある時間後における核種の原子数密度を予測Reactor Physics, Nagoya University, Akio YAMAMOTO

87

原子炉寿命の予測

原子炉の寿命を予測するためには、炉心内の全ての燃料について繰り返し燃焼計算を実施する必要がある。
Σ = σN
巨視的断面積の更新
12
 中性子束の計算：中性子束の更新 − D∇ φ + Σ aφ = νΣ f φ
kdN 燃焼計算：原子数密度の更新= −λN + γΣ f φ


繰り返
し計算

dt



上記の計算を繰り返し、炉心の実効増倍率が1.0になる時点を原子炉の寿命とする。

Reactor Physics, Nagoya University, Akio YAMAMOTO

88

原子炉寿命の予測

実際の原子炉における計算は非常に複雑である
ため、大型の計算プログラム(数万～数十万行)を用いて計算を実行する。
典型的なPWRにおける燃焼に伴う反応度変化
0.14

PWRにおける炉心寿命計算例

0.12

炉心反応度

0.10
0.08
0.06
0.04
0.02
0.00
0

2

4

6

Reactor Physics, Nagoya University, Akio YAMAMOTO810
運転日数(ヶ月)

12

14

16

89

燃焼行列と固有値・固有ベクトル

燃焼行列の固有値と固有ベクトル

燃焼方程式の解析解

初期組成を固有ベクトルで展開Au n = λn u nN (t ) = exp( Αt ) N (0)
N (0) = ∑ an u nnReactor Physics, Nagoya University, Akio YAMAMOTO

90

燃焼行列と固有値・固有ベクトル

解析解に代入N (t ) = exp( Αt )∑ an u nn( At ) 2= I + At ++ ... ∑ an u n
2!

 n
2


(
λn t )= ∑ an 1 + λn t ++ ...u n
2!
n



= ∑ an exp(λn t )u nnReactor Physics, Nagoya University, Akio YAMAMOTO

91

燃焼行列と固有値・固有ベクトル 組成 u n は exp(λn t )で減少する。
絶対値が小さい固有値ほど減少が少ない。
 固有値と固有ベクトルは燃焼行列Aに固有の量

核変換処理の指標として使える?

Reactor Physics, Nagoya University, Akio YAMAMOTO

92

燃焼行列と固有値・固有ベクトルMVP-BURN重核種(Th232-Cm246)チェーンに対する燃焼マトリックスの固有値・固有ベクトル
Half life

Eigenvalue-2.60E-10 -7.80E-10 -1.03E-09 -2.31E-09
1.4E+10 YEARS
TH232
2.4E-11 -1.0E+00
1.4E-12-5.2E-07
3.3E+04 YEARS
PA231
2.4E-15-1.0E-04
1.4E-16-6.0E-11
2.7E+01
DAYS
PA233
6.2E-14-2.5E-03
3.5E-15-1.3E-09
6.9E+01 YEARS
U-232
8.8E-10-2.5E-04-8.9E-13
3.0E-06
1.6E+05 YEARS
U-233
3.0E-10-5.2E-02-3.2E-12
7.4E-072.5E+05 YEARS
U-234
2.1E-05-2.1E-02-3.5E-07
2.9E-02
7.0E+08 YEARS
U-235
9.8E-06-1.0E-02-1.7E-07
1.6E-02
2.3E+07 YEARS
U-236
1.4E-05-1.5E-02-3.7E-07
8.4E-01
6.8E+00
DAYS
U-237
1.0E-06-2.9E-05-7.0E-10
1.6E-03
4.5E+09 YEARS
U-238
1.0E+00
5.1E-04
5.6E-08-7.1E-03
1.5E+05 YEARS
NP236
1.5E-10-4.6E-09-1.3E-13
3.2E-07
2.1E+06 YEARS
NP237
1.4E-04-4.1E-03-1.2E-07
2.8E-01
2.4E+00
DAYS
NP239
6.8E-05
3.4E-08-3.8E-10 -4.5E-07
2.9E+00 YEARS
PU236
6.8E-10-2.1E-08-6.0E-13
1.5E-06
8.8E+01 YEARS
PU238
4.6E-04-4.8E-03-6.6E-06
4.0E-01
2.4E+04 YEARS
PU239
5.4E-03-8.7E-04-1.2E-06
7.6E-02
6.6E+03 YEARS
PU240
1.5E-03-2.4E-04-2.6E-05
2.4E-02
1.4E+01 YEARS
PU241
2.4E-03-3.9E-04-4.3E-05
4.0E-02
3.7E+05 YEARS
PU242
2.6E-03-4.6E-04-6.2E-04
5.8E-02
4.3E+02 YEARS
AM241
1.3E-04-2.1E-05-2.3E-06
2.3E-03
1.6E+01 HOURS
AM242
2.7E-07-4.5E-08-4.9E-09
4.9E-06
1.4E+02 YEARS AM242M 2.2E-06-3.8E-07-4.2E-08
4.1E-05
7.4E+03 YEARS
AM243
1.8E-03-3.2E-04-4.4E-04
4.6E-02
1.6E+02
DAYS
CM242
5.3E-05-8.9E-06-9.9E-07
1.0E-03
2.9E+01 YEARS
CM243
2.4E-06-4.2E-07-6.3E-08
5.1E-05
1.8E+01 YEARS
CM244
4.5E-03-9.1E-04-1.3E-03
1.9E-01
Reactor Physics, Nagoya University,
Akio YAMAMOTO
8.5E+03
YEARS
CM245
5.5E-04-1.1E-04-1.7E-04
2.4E-02
4.7E+03 YEARS
CM246
3.3E-03-2.1E-03-1.0E+00 -8.7E-02
Nuclide93

[復習]原子炉の出力と燃料の消費(1)P MWの運転中の原子炉における核分裂率：

P[MW] × 10 6 [= 5.40 × 10 23


J
1
MeV
1 fission
sec
]×
[
]
[
]
86400
[
]
×
×
−13
MW ⋅ sec 1.6 × 10
Jday
E R MeVP fission[]
E R day

核分裂性物質の原子量をAとすると：燃料消費 = 5.40 × 10

23g
P fission
A
[
]×
[
]
23
6.02 × 10 fission
E R day

PA g[]= 0.895
E R day
Reactor Physics, Nagoya University, Akio YAMAMOTO

94

[復習]原子炉の出力と燃料の消費(2)

U235が燃料の場合、ER=200[MeV], A=235。
従って、燃料消費 = 1.05P[

g
]
day

3000MWの典型的な軽水炉なら、約3kg/day。
 ちなみに、3000MW・dayの電気料金はだいた
い：


Yen
] = 14.4億円
3000 × 10 × 24[kWh ] × 20[
kWh3

Reactor Physics, Nagoya University, Akio YAMAMOTO

95

まとめ
燃焼時に起こっている物理現象
 燃焼チェーン
 燃焼方程式の数値解法
 原子炉における寿命計算
 燃焼行列の固有値と固有ベクトル

Reactor Physics, Nagoya University, Akio YAMAMOTO

96