平成 30 年度

修士論文

全ての核分裂生成物核種を陽に取り扱った
原子炉動特性モデルの開発
Development of Nuclear Reactor Kinetic
Calculation Model with
the Explicit Fission Products Model

原子炉工学研究室
片桐耕司

北海道大学大学院

工学院

エネルギー環境システム専攻

目次
第1章

序論 .......................................................................................................................1

1.1

炉心安全解析.............................................................................................................1

1.2

原子炉の動特性 .........................................................................................................2

1.3

EFP モデル ................................................................................................................4

1.4

本論文の構成.............................................................................................................5

第2章

計算手法 ................................................................................................................6

2.1

中性子拡散方程式の空間離散化 ...............................................................................6

2.2

内部反復 ..................................................................................................................10

2.3

外部反復 ..................................................................................................................11

2.4

二次元体系 ..............................................................................................................122.5

動特性計算 ..............................................................................................................132.6

遅発中性子先行核密度時間変化の準解析的解法 ....................................................16

2.7

EFP モデルの適用 ...................................................................................................172.8

部分行列指数法 .......................................................................................................18

2.9

𝜃法による離散化 .....................................................................................................202.10

空間依存動特性方程式への部分行列指数法の適用 ............................................212.11行列指数の近似 .......................................................................................................212.12

加速法 ..................................................................................................................222.12.1Gauss-Seidel 法 .................................................................................................222.12.2

Successive Over Relaxation 法............................................................................222.12.3

Simplified Fixed Source Scaling Factor 法 ...........................................................23

2.13

EFP モデルで使用している核データ ......................................................................23

2.14

EFP モデルによるガス状 FP の漏洩計算 ................................................................24

2.15

本章のまとめ .......................................................................................................24

第3章

結果 .....................................................................................................................26

3.1

部分行列指数法と𝜃法による一次元動特性解析 .....................................................26

3.2

部分行列指数法と𝜃法の計算精度 ...........................................................................34

3.3

TWIGL ベンチマーク問題.......................................................................................35

3.4

ガス状 FP の漏洩事象(炉心全体) ............................................................................433.5

ガス状 FP の漏洩事象(一部領域) ............................................................................48

3.6

本章のまとめ...........................................................................................................52

第4章

結論 .....................................................................................................................53

参考文献 ................................................................................................................................54

記号表
𝒓

: 位置

𝑣

: 速度

𝑡

: 時間

Σ𝑟

: 巨視的除去断面積

𝐷: 拡散係数Σ𝑠: 巨視的散乱断面積

𝑆: 外部中性子源Σ𝑎: 巨視的吸収断面積

𝜈: 平均中性子放出数

Σ𝑓: 巨視的核分裂断面積

𝑄: 中性子源𝑘𝑒𝑓𝑓: 実効中性子増倍率

𝜒: 核分裂スペクトル𝜙: 中性子束𝜆: 崩壊定数𝛽: 遅発中性子割合𝐶: 遅発中性子先行核密度

𝑃𝑛: 遅発中性子放出確率𝑦: 独立核分裂収率𝑰: 単位行列第1章 序論
1.1 炉心安全解析
原子炉の安全な運転において、原子炉の過渡事象の解析は極めて重要な役割を占めている。特に軽水炉の安全評価においては、運転時の異常な過渡変化や事故について代表的な事象を選定した上で、これらの事象一つ一つに対して炉心を含むプラント全体の挙動が評価され、安全性の確認が行われている。以上の安全性の確認を行うための安全解析は原子炉の設置及び主要設備変更の際に実施されている。この際、運転時の異常な過渡変化と事故の二種類を評価しなければならないことが原子炉の安全設計審査指針の中で明記されている[1]。
特に運転時の異常な過渡変化は原子炉の運転中において、単一の故障もしくは誤作動ま
たは運転員の単一の誤操作、及びこれらと類似の頻度で発生すると予想される外乱によって生ずる異常な状態に至る事象を指す。具体的には、加圧水型軽水炉(PWR)、沸騰水型軽水
炉(BWR)共に以下のような事象が想定されている。
✓

炉心内の反応度または出力分布の異常な変化✓

炉心内の熱発生または熱除去の異常な変化✓

原子炉冷却材圧力または原子炉冷却材保有量の異常な変化✓

その他原子炉施設の設計により必要と認められる事象

安全解析は以上の事象が生じた場合、炉心は損傷に至ることなく、かつ原子炉施設は通常運転に復帰できる状態で事象が収束することを確認するものである。このときの判断基準は以下の通りである。
✓

最小限界熱流速比(Departure from Nuclear Boiling Ratio)または最小限界出力比(Minimum Critical Power Ratio)が許容限界値を超えないこと✓

燃料被覆管は機械的に破損しないこと✓

燃料エンタルピは許容限界値以下であること✓

原子炉冷却材圧力バウンダリにかかる圧力は、最高使用圧力の 1.1 倍以下であること

安全解析では以上の各種制限値が全て満たされることをプラント挙動の解析、伝熱計算、熱水力計算を通して確認することで、運転時の異常な過渡変化時における炉心及び燃料の健全性を担保している。
一方、事故は運転時の異常な過渡変化を超える異常な状態であり、発生する頻度はまれであるが、発生した場合は原子炉施設からの放射性物質の放出の可能性がある事象とされている。評価すべき事象としては、事故が発生した場合における工学的安全施設等、主として事故の影響を緩和する系に属する構築物、系統及び機器の設計の妥当性を確認する観点か
ら、以下が想定されている。
✓

原子炉冷却材の喪失又は炉心冷却状態の著しい変化✓

反応度の異常な投入又は原子炉出力の急激な変化
1

✓環境への放射性物質の異常な放出✓

原子炉格納容器内圧力、雰囲気等の異常な変化✓

その他原子炉施設の設計により必要と認められる事象

また、事故時の安全性に関する判断基準は以下の通りである。
✓

炉心は著しい損傷に至ることなく、かつ、十分な冷却が可能であること。

✓

燃料エンタルピは限界値を超えないこと。

✓

原子炉冷却材圧力バウンダリにかかる圧力は、最高使用圧力の 1.2 倍以下であること。

✓

原子炉格納容器バウンダリにかかる圧力は、最高使用圧力以下であること。

✓

周辺の公衆に対し、著しい放射線被ばくのリスクを与えないこと。

以上のように安全解析には様々な制約条件があり、これらの制約条件に対し炉心が健全性を担保出来ることを確認する必要がある。

1.2 原子炉の動特性
現行の安全解析では、異常な過渡変化及び事故の発生に対して、プラント特性解析、炉心動特性解析、熱流動特性解析といった複数の解析コードを複合し、原子炉全体の変化を追従している。この中でも炉心動特性は、核計算から炉心の出力及びその分布の変化を追従するものであり、重要な役割を担っている。
炉心では、U-235 などの核分裂性物質が中性子によって核分裂し、それによって発生した中性子によって連鎖的に反応が継続することで、出力を維持している。この核分裂反応時に発生する核種は核分裂生成物(Fission Products)と呼ばれ、1000 種類以上存在している。この
核分裂生成物の中には、崩壊により核分裂反応から時間遅れで中性子を放出する核種が存在する。この核種のことを遅発中性子先行核と呼び、また放出される中性子のことを遅発中性子と呼ぶ。遅発中性子の割合は、核分裂核種、入射中性子のエネルギーによって異なるが、ほとんどの場合 1%以下とわずかである。しかし、原子炉の制御という点で非常に重要である。例として、臨界定常状態に正の反応度を印加した際の中性子数変化を、遅発中性子の有無で比較した結果を Fig. 1.2.1 に示す。Fig.1.2.1 より、遅発中性子により中性子数の増加が抑えられていることが分かる。このことから遅発中性子により原子炉の制御が可能となっていることが分かる。

2

Fig. 1.2.1 遅発中性子の有無による中性子数の時間変化
遅発中性子先行核は数多く存在し、それぞれが異なる時間遅れで崩壊し遅発中性子を放出する。例えば、Br-87 は 55.6[s]と長い半減期を持ち、一方で Ga-83 は0.3[s]と短い半減期を持つ。このように半減期に幅を持つ多くの核種を簡略化して取り扱うため、現在一般的に用いられている動特性解析では、これらの遅発中性子先行核を半減期により 6 つのグループにグルーピングする 6 群モデルが採用されている。Table 1.2.1 に U-235 熱中性子核分裂における 6 群モデルのグルーピングの一例を示す[2]。
Table 1.2.1 6 群モデルグルーピングの例

先行核生成割合グループ

半減期T1⁄2 [s]崩壊定数𝜆𝑖 [s−1 ]

1

55.72

0.0124

0.052

2

22.72

0.0305

0.346

3

6.22

0.111

0.310

4

2.30

0.301

0.624

5

0.610

1.14

0.182

6

0.230

3.01

0.066

𝑎𝑖

6 群モデルは全ての先行核を、仮想的な 6 種類の先行核として扱うため、計算負荷が小さ
く、また取り扱いが容易であるなどのメリットが存在する。6 群モデルでは、それぞれの先行核は崩壊により中性子を放出し、中性子を放出しない他の核種になると仮定されている。
しかし、実際の反応では各々の核種は崩壊を繰り返し、中性子を放出する可能性がある他の核種へと変換されていく。6 群モデルではこのような物理現象を詳細に取り扱うことが出来ない。この他にも、6 群モデルでは直接的に取り扱うことの出来ない物理現象として、ガス
3

状 FP の放出や、遅発ガンマ線の放出などが挙げられる。

1.3 EFP モデル
前節で述べたように、6 群モデルは簡略化されたモデルであるため、考慮することが出来ない物理現象が存在する。そこで本研究では物理現象を詳細に取り扱うことが出来るモデルとして、EFP(Explicit Fission Products)の開発を行った。EFP モデルでは遅発中性子先行核だけでなく、1300 種類以上存在する核分裂生成物核種を全て陽に取り扱う。EFP モデルの概略図を Fig. 1.3.1 に示す。

Fig. 1.3.1 6 群モデルと EFP モデルの概略図
EFP モデルでは、各々の FP についてのデータを評価済み核データファイルから取得し、利用している。利用するデータとして、核分裂収率や半減期、崩壊分岐比などが挙げられる。
これにより、実現象と同様に、核種が崩壊し他核種へと移行する現象も再現することが出来る。また、各核種を陽に扱うため、それぞれの同位体が持つ物理的、化学的な性質を反映させることが出来る。この例として挙げられるのが、前節で述べたガス状 FP の漏洩事象である。ヨウ素やキセノンなどの FP は揮発性を持ち、事故時には出力の上昇とともに増加し、最終的に漏洩する事が想定される。しかし、6 群モデルでは核種を群としてまとめて取り扱
うため、この現象を再現することが難しい。一方、EFP モデルでは各核種を陽に取り扱うため、漏洩により核種数密度が減少していく設定をすることが出来る。このように EFP モデ
ルを用いることで、複雑な事象を正確に再現することが出来るといった利点が存在する。
先行研究では、この EFP モデルを一点炉動特性解析に適用させ、中性子数時間変化について、6 群モデルとの比較を行った。正負の反応度印加や、中性子パルス照射など、様々なケースについて解析を行ったが、いずれも有意な差は見られなかった。本研究では、空間依存解析へ EFP モデルを適用させ、6 群モデルでは扱えない事象の再現を目指す。また、そのための手法開発及び計算手法の高度化も目的の一つとする。

4

1.4 本論文の構成
本論文は全 4 章で構成される。第 2 章では、一般的な空間依存動特性解析で用いられている計算手法について、また EFP モデルを空間依存動特性解析へと適用させる際に必要な新たな手法の提案を行う。第 3 章では、一次元動特性解析において新たに提案した計算手法の有効性評価及び、体系を 2 次元に拡張し、ガス状 FP の漏洩事象といった EFP モデルで扱うことの出来る複雑な現象についての解析を行う。最後に第 4 章で結論を述べる。

5

第2章 計算手法
2.1 中性子拡散方程式の空間離散化
中性子拡散方程式の解法を示す。まず初めに、中性子拡散方程式の空間離散化手法について説明するため、固定源多群中性子拡散問題について考える。時間依存の多群中性子拡散方程式は以下の(2.1.1)式のように表される
1 𝜕𝜙𝑔 (𝒓, 𝑡)= 𝛻 ∙ 𝐷𝑔 (𝒓, 𝑡)𝛻𝜙𝑔 (𝒓, 𝑡) − 𝛴𝑟,𝑔 (𝒓, 𝑡)𝜙𝑔 (𝒓, 𝑡)𝑣𝑔
𝜕𝑡(2.1.1)𝐺+ ∑ 𝛴𝑠,𝑔′ →𝑔 (𝒓, 𝑡)𝜙𝑔′ (𝒓, 𝑡) + 𝑆𝑔 (𝒓, 𝑡) ,𝑔′ ≠𝑔

本節では定常計算での計算手法について記述する。また、簡略化のため一次元体系につい𝜕𝜙𝑔て考える。この時、固定源(定常)方程式は、 𝜕𝑡 = 0より、(2.1.2)式のようになる。
𝐺𝑑𝜙𝑔 (𝑥)𝑑− (𝐷𝑔 (𝑥)) + 𝛴𝑟,𝑔 (𝑥)𝜙𝑔 (𝑥) = ∑ 𝛴𝑠,𝑔′ →𝑔 (𝑥)𝜙𝑔′ (𝑥) + 𝑆𝑔 (𝑥) ,𝑑𝑥
𝑑𝑥′(2.1.2)𝑔 ≠𝑔

この(2.1.2)式内の各変数はそれぞれ空間的に連続な値を持っている。そのため、計算機に解かせる際には各変数を離散化して取り扱う必要がある。以下では本研究で用いられている離散化手法について述べる。
まず、体系全体を𝑁個のメッシュに分割し、あるメッシュ𝑖の幅を∆𝑥𝑖 とおく。また、メッシュ中点の座標を𝑥𝑖 とする。これより、メッシュ左側境界は−𝑥𝑖 /2、右側境界は+𝑥𝑖 /2と表される。模式図を Fig. 2.1.1 に示す。

Fig. 2.1.1 メッシュ分割(2.1.2)式をメッシュ𝑖内で積分することで、以下の(2.1.3)式を得ることが出来る。

6

−𝐷𝑔 (𝑥)𝑑𝜙𝑔 (𝑥)𝑑𝜙𝑔 (𝑥)|+ −𝐷𝑔 (𝑥)|+ 𝛴𝑟,𝑔,𝑖 𝜙𝑔,𝑖 ∆𝑥𝑖
𝑑𝑥 𝑥=∆𝑥 /2
𝑑𝑥 𝑥=−∆𝑥 /2𝑖𝑖(2.1.3)𝐺= ∑ 𝛴𝑠,𝑔′ →𝑔,𝑖 𝜙𝑔′ ,𝑖 ∆𝑥𝑖 + 𝑆𝑔,𝑖 ∆𝑥𝑖 ,𝑔′ ≠𝑔

ここで、𝐽𝑔 (𝑥) = −𝐷𝑔 (𝑥)𝑑𝜙𝑔 (𝑥)𝑑𝑥で表される中性子流を定義する。また、𝜙𝑔,𝑖 はメッシュ𝑖での平均中性子束を表し、(2.1.4)式で定義される。
𝜙𝑔,𝑖 =

1 ∆𝑥𝑖 /2∫
𝜙 (𝑥)𝑑𝑥 ,∆𝑥𝑖 −∆𝑥𝑖/2 𝑔(2.1.4)

また、各種断面積もメッシュ内の反応率を保存するよう、(2.1.5)式のように定義される。
∆𝑥 /2

Σ𝑔,𝑖 =∫−∆𝑥𝑖 /2 Σ𝑔 (𝑥)𝜙𝑔 (𝑥)𝑑𝑥
𝑖

∆𝑥 /2

∫−∆𝑥𝑖 /2 𝜙𝑔 (𝑥)𝑑𝑥(2.1.5),𝑖

中性子流を用いると、(2.1.3)式は(2.1.6)式のように書き換えることが出来る。
𝐽𝑔 (𝑥𝑖 +∆𝑥𝑖∆𝑥𝑖) − 𝐽𝑔 (𝑥𝑖 −) + 𝛴𝑟,𝑔,𝑖 𝜙𝑔,𝑖 ∆𝑥𝑖 = 𝑆𝑔,𝑖 ∆𝑥𝑖 ,2
2(2.1.6)

さらに、Fig. 2.1.2 で示されるような、メッシュ境界での中性子流を(2.1.7)式のように表すことで、(2.1.6)式は(2.1.8)式のように変形することが出来る。

Fig. 2.1.2 中性子流の定義

𝑥−𝐽𝑔,𝑖= 𝐽𝑔 (−

∆𝑥𝑖
),
2

∆𝑥𝑖
𝑥+
𝐽𝑔,𝑖
= 𝐽𝑔 ( ) ,2(2.1.7)

𝐺
𝑥+
𝑥−
𝐽𝑔,𝑖
− 𝐽𝑔,𝑖+ 𝛴𝑟,𝑔,𝑖 𝜙𝑔,𝑖 ∆𝑥𝑖 = ∑ 𝛴𝑠,𝑔′ →𝑔,𝑖 𝜙𝑔′ ,𝑖 ∆𝑥𝑖 + 𝑆𝑔,𝑖 ∆𝑥𝑖 ,(2.1.8)𝑔′ ≠𝑔𝑑𝜙

次に、 𝑑𝑥 の計算方法について考える。着目メッシュ𝑖と隣接するメッシュ(𝑖 + 1)の境界に∆𝑥𝑥+
おける中性子流は、前述したとおり𝐽𝑔,𝑖= −𝐷𝑔 (𝑥𝑖 + 2 𝑖 )7𝑑𝜙𝑔 (𝑥)𝑑𝑥|

𝑥=𝑥𝑖 +∆𝑥𝑖2で表される。ここで、

𝑥+
メッシュ境界における中性子束を𝜙𝑏𝑥+ 、メッシュ𝑖の媒質の拡散係数を𝐷𝑔,𝑖 とすると、𝐽𝑔,𝑖はメ

ッシュ𝑖側の差分近似により、(2.1.9)式のように表すことが出来る。
−

𝑥+
(𝐽𝑔,𝑖
) ≈ −𝐷𝑔,𝑖𝑥+𝜙𝑔,𝑏− 𝜙𝑔,𝑖,∆𝑥𝑖 /2(2.1.9)𝑥+
同様に、𝐽𝑔,𝑖はメッシュ(𝑖 + 1)側の差分近似でも(2.1.10)式のように表すことが出来る。
+

𝑥+
(𝐽𝑔,𝑖
) ≈ −𝐷𝑔,𝑖+1𝑥+
𝜙𝑔,𝑖+1 − 𝜙𝑔,𝑏,∆𝑥𝑖+1 /2−

(2.1.10)+𝑥+𝑥+
メッシュ境界において中性子流は連続であるため((𝐽𝑔,𝑖) = (𝐽𝑔,𝑖) )、(2.1.9)式、(2.1.10)式を用いることで、メッシュ𝑖とメッシュ(𝑖 + 1)境界の中性子束𝜙𝑏𝑥+ は以下の(2.1.11)式のように
表される。
𝑥+
𝜙𝑔,𝑏
=

∆𝑥𝑖 𝐷𝑔,𝑖+1 𝜙𝑔,𝑖+1 + ∆𝑥𝑖+1 𝐷𝑔,𝑖 𝜙𝑔,𝑖,∆𝑥𝑖 𝐷𝑔,𝑖+1 + ∆𝑥𝑖+1 𝐷𝑔,𝑖(2.1.11)

同様の差分近似をメッシュ(𝑖 − 1)とメッシュ𝑖境界にも適用することで、メッシュ(𝑖 − 1)とメッシュ𝑖境界の中性子束𝜙𝑏𝑥− は(2.1.12)式のように表される。
𝑥−
𝜙𝑔,𝑏
=

∆𝑥𝑖−1 𝐷𝑔,𝑖 𝜙𝑔,𝑖 + ∆𝑥𝑖 𝐷𝑔,𝑖−1 𝜙𝑔,𝑖−1,∆𝑥𝑖−1 𝐷𝑔,𝑖 + ∆𝑥𝑖 𝐷𝑔,𝑖−1(2.1.12)(2.1.11)式を(2.1.9)式に代入することで、中性子流は隣接する 2 メッシュの中性子束を用いて計算することが出来る。
𝑥+𝐽𝑔,𝑖≈ −𝐷𝑔,𝑖𝑥+𝜙𝑔,𝑏− 𝜙𝑔,𝑖
2𝐷𝑔,𝑖 𝐷𝑔,𝑖+1=−(𝜙− 𝜙𝑔,𝑖 ) ,∆𝑥𝑖 /2𝐷𝑔,𝑖 ∆𝑥𝑖+1 + 𝐷𝑔,𝑖+1 ∆𝑥𝑖 𝑔,𝑖+1(2.1.13)

同様に(2.1.12)式より、メッシュ(𝑖 − 1)とメッシュ𝑖境界の中性子流も隣接する 2 メッシュの中性子束を用いて表すことが出来る。
𝑥−𝐽𝑔,𝑖≈ −𝐷𝑔,𝑖𝑥−𝜙𝑔,𝑖 − 𝜙𝑔,𝑏2𝐷𝑔,𝑖−1 𝐷𝑔,𝑖=−(𝜙 − 𝜙𝑔,𝑖−1 ) ,∆𝑥𝑖 /2𝐷𝑔,𝑖−1 ∆𝑥𝑖 + 𝐷𝑔,𝑖 ∆𝑥𝑖−1 𝑔,𝑖(2.1.14)

最後に(2.1.13)式及び(2.1.14)式を(2.1.8)式に代入することで、着目メッシュと隣接する 2 メッシュの中性子束からなる 3 点階差の差分式が導出出来る。
𝑥+
𝑥−
(𝐶𝑔,𝑖
+ 𝐶𝑔,𝑖+ 𝛴𝑟,𝑔,𝑖 ∆𝑥𝑖 )𝜙𝑔,𝑖𝐺𝑥+𝑥−= 𝐶𝑔,𝑖𝜙𝑔,𝑖+1 + 𝐶𝑔,𝑖𝜙𝑔,𝑖−1 + ∑ 𝛴𝑠,𝑔′ →𝑔,𝑖 𝜙𝑔′ ,𝑖 ∆𝑥𝑖 + 𝑆𝑔,𝑖 ∆𝑥𝑖 ,(2.1.15)𝑔′ ≠𝑔

但し、(2.1.15)式中の各係数は(2.1.16)式、(2.1.17)式で示す通りである。
𝑥+
𝐶𝑔,𝑖
=

2𝐷𝑔,𝑖 𝐷𝑔,𝑖+1,
𝐷𝑔,𝑖 ∆𝑥𝑖+1 + 𝐷𝑔,𝑖+1 ∆𝑥𝑖(2.1.16)

𝑥−
𝐶𝑔,𝑖
=

2𝐷𝑔,𝑖−1 𝐷𝑔,𝑖,𝐷𝑔,𝑖−1 ∆𝑥𝑖 + 𝐷𝑔,𝑖 ∆𝑥𝑖−1(2.1.17)

但し、 (2.1.15)式は炉心境界においては異なる取り扱いが必要である。ここでは、零中性子
8束境界条件及び完全反射境界条件の式を導出する。
(1)零中性子束境界条件
零中性子束境界条件は体系表面で中性子束が零となる境界条件である。この際の炉心境
界における中性子流を求める。まず、炉心境界における中性子束について、外挿した直線を考える。この際、中性子束が炉心外で0となる距離(外挿距離)を𝑑とする。炉心境界における中性子束と外挿距離のイメージを Fig. 2.1.3 に示す。

Fig. 2.1.3 炉心表面の中性子束と外挿距離
この外挿直線の式は(2.1.18)で表される。但し、炉心表面の中性子束を𝜙𝑠 とする。
𝜙≈

𝜙𝑔,0 − 𝜙𝑠
𝑥 + 𝜙𝑠 ,∆𝑥0 /2(2.1.18)

𝑥 = 𝑑において、𝜙 = 0となるため、(2.1.18)式に代入することで、外挿距離と炉心表面中性子束に関する(2.1.19)式を得ることが出来る。
𝜙𝑠 =

𝑑
𝜙 ,
∆𝑥0 /2 + 𝑑 𝑔,0(2.1.19)

よって、炉心境界における中性子流は以下の(2.1.20)式で表すことが出来る。
𝑥−𝐽𝑔,0= −𝐷𝑔,0

𝐷𝑔,0
𝑑𝜙
𝜙 ,
| ∆𝑥 ≈ −0𝑑𝑥 𝑥=∆𝑥0 /2 + 𝑑 𝑔,0(2.1.20)2

零中性子束境界条件では、外挿距離を0としたものと等価であるため、零中性子束境界条件
における炉心表面の中性子流は(2.1.21)式で表すことが出来る。
𝑥−
𝐽𝑔,0
=−

2𝐷𝑔,0
𝜙 ,
∆𝑥0 𝑔,0(2.1.21)

よって、炉心境界メッシュにおける(2.1.15)式は以下の(2.1.22)式及び(2.1.23)式の様になる。

9

(ⅰ)𝑖 = 0𝑥+(𝐶𝑔,0+ 2𝐷𝑔,0 /∆𝑥0 + 𝛴𝑟,𝑔,0 ∆𝑥0 )𝜙𝑔,0𝐺𝑥+= 𝐶𝑔,0
𝜙𝑔,1 + ∑ 𝛴𝑠,𝑔′ →𝑔,0 𝜙𝑔′ ,0 ∆𝑥0 + 𝑆𝑔,0 ∆𝑥0 ,𝑔′ ≠𝑔(2.1.22)(ⅱ)𝑖 = 𝐼𝑥−(2𝐷𝑔,𝐼 /∆𝑥𝐼 + 𝐶𝑔,𝐼
+ 𝛴𝑟,𝑔,𝐼 ∆𝑥𝐼 )𝜙𝑔,𝐼𝐺𝑥−= 𝐶𝑔,𝐼
𝜙𝑔,𝐼−1 + ∑ 𝛴𝑠,𝑔′ →𝑔,𝐼 𝜙𝑔′ ,𝐼 ∆𝑥𝐼 + 𝑆𝑔,𝐼 ∆𝑥𝐼 ,(2.1.23)𝑔′ ≠𝑔(2)完全反射境界条件
完全反射境界条件は体系から出ていく中性子が全て反射され戻ってくるという境界条件である。これは一般的に対称性のある炉心の一部を計算する際に用いられる。そのため、境界面での正味の中性子流は0となる。これは外挿距離を∞としたものと等価である。よって、完全反射条件における中性子束の式は以下の(2.1.24)式及び(2.1.25)式の様になる。
(ⅰ)𝑖 = 0
𝐺
𝑥+
𝑥+
(𝐶𝑔,0
+ 𝛴𝑟,𝑔,0 ∆𝑥0 )𝜙𝑔,0 = 𝐶𝑔,0
𝜙𝑔,1 + ∑ 𝛴𝑠,𝑔′ →𝑔,0 𝜙𝑔′ ,0∆𝑥0 + 𝑆𝑔,0 ∆𝑥0 ,(2.1.24)𝑔′ ≠𝑔(ⅱ)𝑖 = 𝐼
𝐺
𝑥−
𝑥−
(𝐶𝑔,𝐼
+ 𝛴𝑟,𝑔,𝐼 ∆𝑥𝐼 )𝜙𝑔,𝐼 = 𝐶𝑔,𝐼
𝜙𝑔,𝐼−1 + ∑ 𝛴𝑠,𝑔′ →𝑔,𝐼 𝜙𝑔′ ,𝐼 ∆𝑥𝐼 + 𝑆𝑔,𝐼 ∆𝑥𝐼 ,(2.1.25)𝑔′ ≠𝑔2.2 内部反復
上記では一次元中性子拡散方程式について、3 点階差式の導出を行った。では次に実際の計算で採用されている反復解法について述べる。(2.1.15)式で右辺第三項は散乱により、他のエネルギー群から供給される中性子を表している。上群から下群への散乱のみを仮定し、中性子束を上群から計算していった場合、この項は既知となる。そのため、(2.1.15)式の右辺第三項及び第四項をまとめて(2.2.1)式の様に扱う。
𝐺

𝑆𝑔,𝑖 ∆𝑥𝑖 = ∑ 𝛴𝑠,𝑔′ →𝑔,𝑖 𝜙𝑔′ ,𝑖 ∆𝑥𝑖 + 𝑆𝑔,𝑖 ∆𝑥𝑖 ,(2.2.1)𝑔′ ≠𝑔

これにより、(2.1.15)式は行列形式で(2.2.2)式の様に表記される。
𝑨𝝓 = 𝑺 ,10(2.2.2)(2.2.2)式は𝑨の逆行列を計算することで解を得ることが出来る。しかし、行列𝑨のサイズが大きい場合、ガウスの消去法等で逆行列を求めることが困難である。このような場合、逆行列を用いた手法ではなく、反復解法が用いられる。反復解法では、はじめに中性子束分布を仮定し、各メッシュについて(2.1.15)式を使って𝜙𝑖 を計算する。この時、𝜙𝑖−1、𝜙𝑖+1 については仮定された値を用いる。これを全てのメッシュについて繰り返し行う。例えば、𝑘回目、𝑘 +1回目の反復における中性子束をそれぞれ𝜙 (𝑘) 、𝜙 (𝑘+1) とすると𝜙 (𝑘+1) は(2.2.3)式の様に記述
される。
(𝑘+1)

𝑥+
𝑥−
(𝐶𝑔,𝑖
+ 𝐶𝑔,𝑖+ 𝛴𝑟,𝑔,𝑖 ∆𝑥𝑖 )𝜙𝑔,𝑖

(𝑘)

(𝑘)

𝑥+
𝑥−
= 𝐶𝑔,𝑖𝜙𝑔,𝑖+1 + 𝐶𝑔,𝑖𝜙𝑔,𝑖−1 + 𝑆𝑔,𝑖 ∆𝑥𝑖 ,(2.2.3)

この反復を解が収束するまで繰り返す。なお収束の判定には各メッシュにおける𝜙 (𝑘) 、𝜙 (𝑘+1) の相対差の最大値を用いる。
(𝑘+1)

𝜙
𝜖 = max | 𝑖
𝑖

(𝑘)

− 𝜙𝑖

(𝑘)

𝜙𝑖

(2.2.4)|,散乱による中性子源については、ある群の内部反復が収束した段階で、その群から下群への散乱を計算し、中性子源に加える操作を行う。また、第一群の散乱源については0として扱
い、中性子源は外部中性子源のみとする。これを全てのエネルギー源について適用させることで中性子束分布を求めることが出来る。

2.3 外部反復
次に核分裂中性子源がある場合の数値解法について考える。この場合、多群中性子拡散固有値方程式は(2.3.1)式で記述される。
−

𝑑𝜙𝑔 (𝑥)𝑑(𝐷𝑔 (𝑥)) + 𝛴𝑟,𝑔 (𝑥)𝜙𝑔 (𝑥)
𝑑𝑥
𝑑𝑥
𝐺

𝐺

𝑔 ≠𝑔

𝑔

𝜒𝑔
= ∑ Σ𝑠,𝑔′ →𝑔 (𝑥)𝜙𝑔′ (𝑥) +∑ 𝜈Σ ′ (𝑥)𝜙𝑔′ (𝑥) ,𝑘𝑒𝑓𝑓 ′ 𝑓,𝑔′(2.3.1)

この方程式を数値的に解く場合、右辺第二項の核分裂中性子源項に関する反復計算が必要となる。計算手順としては、まず核分裂源の初期値として、適当な分布を仮定する。この仮定値を中性子源として、内部反復を行う。得られた収束解を基に、再び核分裂源の分布を求める。また固有値𝑘𝑒𝑓𝑓 についても同様に、初期値として適当な値を仮定し、反復を行う。固
有値の推定値は体系全体の核分裂源積分値から(2.3.2)式のように求めることが出来る。
′
𝑘𝑒𝑓𝑓= 𝑆 ′ /(𝑆/𝑘𝑒𝑓𝑓 ) ,(2.3.2)′

但し、𝑘𝑒𝑓𝑓 及び𝑆 は現在のステップにおける固有値の推定値及び、内部反復により得られた中性子束分布を基に計算された核分裂源の積分値であり、(2.3.3)式で計算される。
𝑆 = ∑ (∑ 𝜈𝛴𝑓,𝑔,𝑖 𝜙𝑔,𝑖 ) 𝑉𝑖 ,
𝑖

𝑔

11

(2.3.3)

但し、𝑉𝑖 はメッシュ𝑖の体積を示す。こうして得られた核分裂源分布及び、固有値についても相対差の最大値をとり収束条件と比較し収束を判定する。この核分裂源に関する反復計算
を外部反復と呼ぶ。一般的に外部反復の収束条件を内部反復の条件より厳しくすると、外部反復の収束が達成されないことがある。そのため、本研究では内部反復の収束条件は外部反復の条件より一桁厳しく設定している。

2.4 二次元体系
前節まででは、一次元体系での空間離散化手法について述べた。本節では体系を二次元に拡張した場合の空間離散化について述べる。なお簡略化のため 2.1 節と同様に固定源多群中性子拡散問題について考える。二次元体系の場合、固定源多群中性子拡散方程式は以下の(2.4.1)式で表される。
−

𝜕𝜙𝑔 (𝑥, 𝑦)𝜕𝜙𝑔 (𝑥, 𝑦)𝜕𝜕(𝐷𝑔 (𝑥, 𝑦))−(𝐷𝑔 (𝑥, 𝑦)) + 𝛴𝑟,𝑔 (𝑥, 𝑦)𝜙𝑔 (𝑥, 𝑦)
𝜕𝑥
𝜕𝑥
𝜕𝑦
𝜕𝑦
(2.4.1)𝐺= ∑ 𝛴𝑠,𝑔′ →𝑔 (𝑥, 𝑦)𝜙𝑔′ (𝑥, 𝑦) + 𝑆𝑔 (𝑥, 𝑦) ,𝑔′ ≠𝑔(2.4.1)式の左辺第一項及び左辺第二項は 2.1 節と同様の有限差分法を用いて離散化することが出来る。有限差分法を用いることにより、中性子束について着目メッシュと𝑥方向及び𝑦
方向それぞれに隣接するメッシュの 5 つのメッシュを用いた 5 点階差式を得ることが出来る。二次元体系におけるメッシュ分割の定義を Fig. 2.4.1 に、5 点階差式を(2.4.2)式に示す。

Fig. 2.4.1 メッシュ分割(二次元)
12

𝑦−

𝑦+

𝑥+
𝑥−
(𝐶𝑔,𝑖,𝑗+ 𝐶𝑔,𝑖,𝑗+ 𝐶𝑔,𝑖,𝑗 + 𝐶𝑔,𝑖,𝑗 + 𝛴𝑟,𝑔,𝑖,𝑗 ∆𝑥𝑖 ∆𝑦𝑗 )𝜙𝑔,𝑖,𝑗
𝐺𝑥+𝑥−= 𝐶𝑔,𝑖𝜙𝑔,𝑖+1,𝑗 + 𝐶𝑔,𝑖,𝑗𝜙𝑔,𝑖−1,𝑗 + ∑ 𝛴𝑠,𝑔′ →𝑔,𝑖,𝑗 𝜙𝑔′ ,𝑖,𝑗 ∆𝑥𝑖 ∆𝑦𝑗
𝑔′ ≠𝑔(2.4.2)+ 𝑆𝑔,𝑖,𝑗 ∆𝑥𝑖 ∆𝑦𝑗 ,但し、(2.4.2)式中の各変数は以下で定義される。
𝑥+
𝐶𝑔,𝑖,𝑗=2𝐷𝑔,𝑖,𝑗 𝐷𝑔,𝑖+1,𝑗 ∆𝑦𝑗,
𝐷𝑔,𝑖,𝑗 ∆𝑥𝑖+1 + 𝐷𝑔,𝑖+1,𝑗 ∆𝑥𝑖(2.4.3)𝑥−
𝐶𝑔,𝑖,𝑗=2𝐷𝑔,𝑖−1,𝑗 𝐷𝑔,𝑖,𝑗 ∆𝑦𝑗,
𝐷𝑔,𝑖−1,𝑗 ∆𝑥𝑖 + 𝐷𝑔,𝑖,𝑗 ∆𝑥𝑖−1(2.4.4)2𝐷𝑔,𝑖,𝑗 𝐷𝑔,𝑖,𝑗+1 ∆𝑥𝑖,
𝐷𝑔,𝑖,𝑗 ∆𝑦𝑗+1 + 𝐷𝑔,𝑖,𝑗+1 ∆𝑦𝑗(2.4.5)2𝐷𝑔,𝑖,𝑗−1 𝐷𝑔,𝑖,𝑗 ∆𝑥𝑖,
𝐷𝑔,𝑖,𝑗−1 ∆𝑦𝑗 + 𝐷𝑔,𝑖,𝑗 ∆𝑦𝑗−1(2.4.6)𝑦+𝐶𝑔,𝑖,𝑗 =𝑦−𝐶𝑔,𝑖,𝑗 =また、炉心境界での階差式も、一次元の場合と同様に外挿距離を考えることで導出することが出来る。参考に、炉心西側境界メッシュで零中性子束境界条件の場合の階差式を(2.4.7)式に示す。
𝑥+(𝐶𝑔,0,𝑗+2𝐷𝑔,0,𝑗 ∆𝑦𝑗𝑦−𝑦++ 𝐶𝑔,0,𝑗 + 𝐶𝑔,0,𝑗 + 𝛴𝑟,𝑔,0,𝑗 ∆𝑥0 ∆𝑦𝑗 ) 𝜙𝑔,0,𝑗∆𝑥0
𝐺𝑥+= 𝐶𝑔,0𝜙𝑔,1,𝑗 + ∑ 𝛴𝑠,𝑔′ →𝑔,0,𝑗 𝜙𝑔′ ,0,𝑗 ∆𝑥0 ∆𝑦𝑗
𝑔′ ≠𝑔(2.4.7)(𝑗 ≠ 0 , 𝑗 ≠ 𝐽),+ 𝑆𝑔,0,𝑗 ∆𝑥0 ∆𝑦𝑗(2.4.2)式を用いて一次元と同様に内部反復、外部反復を行うことで、二次元体系でも中性子束分布を計算することが出来る。

2.5 動特性計算
前節までの内容は全て、時間的に定常な方程式についての解法であった。ここからは中性子束分布が時間的に変化する場合の中性子拡散方程式及びその計算手法について述べる。
中性子拡散動特性方程式は(2.5.1)式及び(2.5.2)式で表される。簡略化のため、一次元体系について記述する。

13

1 𝜕𝜙𝑔 (𝑥, 𝑡)= ∇ ∙ 𝐷𝑔 (𝑥, 𝑡) ∙ ∇𝜙𝑔 (𝑥, 𝑡) − Σ𝑟,𝑔 (𝑥, 𝑡)𝜙𝑔 (𝑥, 𝑡)𝑣𝑔
𝜕𝑡𝐺+ ∑ Σ𝑠,𝑔′ →𝑔 (𝑥, 𝑡)𝜙𝑔′ (𝑥, 𝑡)𝑔′ ≠𝑔(2.5.1)+ 𝜒𝑝,𝑔 (1 − 𝛽) ∑ 𝜈Σ𝑓,𝑔′ (𝑥, 𝑡)𝜙𝑔′ (𝑥, 𝑡)𝑔′+ ∑ 𝜒𝑑,𝑚,𝑔 𝜆𝑚 𝐶𝑚 (𝑥, 𝑡) ,𝑚

𝜕𝐶𝑚 (𝑥, 𝑡)= 𝛽𝑚 ∑ 𝜈Σ𝑓,𝑔′ (𝑥, 𝑡)𝜙𝑔′ (𝑥, 𝑡) – 𝜆𝑚 𝐶𝑚 (𝑥, 𝑡) ,𝜕𝑡′(2.5.2)𝑔(2.5.1)式の左辺は時間微分項であり、数値計算の際に直接取り扱うのは困難である。そのため、𝜃法などの手法を用いて離散化する必要がある。まず、簡略化のために(2.5.1)式を(2.5.3)
式のように記述する。
1 𝜕𝜙𝑔 (𝑥, 𝑡)= 𝑅𝑔 (𝑥, 𝑡) ,𝑣𝑔
𝜕𝑡

(2.5.3)

但し、式中の𝑅は(2.5.4)式で表される。
𝑅𝑔 (𝑥, 𝑡) = ∇ ∙ 𝐷𝑔 (𝑥, 𝑡) ∙ ∇𝜙𝑔 (𝑥, 𝑡) − Σ𝑟,𝑔 (𝑥, 𝑡)𝜙𝑔 (𝑥, 𝑡)𝐺+ ∑ Σ𝑠,𝑔′ →𝑔 (𝑥, 𝑡)𝜙𝑔′ (𝑥, 𝑡)𝑔′ ≠𝑔(2.5.4)+ 𝜒𝑝,𝑔 (1 − 𝛽) ∑ 𝜈Σ𝑓,𝑔′ (𝑥, 𝑡)𝜙𝑔′ (𝑥, 𝑡)𝑔′+ ∑ 𝜒𝑑,𝑚,𝑔 𝜆𝑚 𝐶𝑚 (𝑥, 𝑡) ,𝑚

(2.5.3)式を、𝜃法を用いて近似することで、(2.5.5)式を得る。
1 𝜙𝑔 (𝑥, 𝑡 + ∆𝑡) − 𝜙𝑔 (𝑥, 𝑡)≈ 𝜃𝑅𝑔 (𝑥, 𝑡 + ∆𝑡) + (1 − 𝜃)𝑅𝑔 (𝑥, 𝑡) ,𝑣𝑔
∆𝑡

(2.5.5)(2.5.5)式を変形することで、(2.5.6)式を得る。
𝜙𝑔 (𝑥, 𝑡 + ∆𝑡)
𝜙𝑔 (𝑥, 𝑡)
1− 𝑅𝑔 (𝑥, 𝑡 + ∆𝑡) = ( − 1) 𝑅𝑔 (𝑥, 𝑡) +,𝜃𝑣𝑔 ∆𝑡𝜃𝜃𝑣𝑔 ∆𝑡(2.5.4)式を(2.5.6)式に代入すると、以下の(2.5.7)式を得ることが出来る。

14(2.5.6)−∇ ∙ 𝐷𝑔 (𝑥, 𝑡 + ∆𝑡) ∙ ∇𝜙𝑔 (𝑥, 𝑡 + ∆𝑡) + {Σ𝑟,𝑔 (𝑥, 𝑡 + ∆𝑡) +1} 𝜙 (𝑥, 𝑡 + ∆𝑡)𝜃𝑣𝑔 ∆𝑡 𝑔𝐺− ∑ Σ𝑠,𝑔′ →𝑔 (𝑥, 𝑡 + ∆𝑡)𝜙𝑔′ (𝑥, 𝑡 + ∆𝑡)𝑔′ ≠𝑔(2.5.7)− 𝜒𝑝,𝑔 (1 − 𝛽) ∑ 𝜈Σ𝑓,𝑔′ (𝑥, 𝑡 + ∆𝑡)𝜙𝑔′ (𝑥, 𝑡 + ∆𝑡)𝑔′

𝜙𝑔 (𝑥, 𝑡)1− ∑ 𝜒𝑑,𝑚,𝑔 𝜆𝑚 𝐶𝑚 (𝑥, 𝑡 + ∆𝑡) = ( − 1) 𝑅𝑔 (𝑥, 𝑡) +,𝜃𝜃𝑣𝑔 ∆𝑡𝑚(2.5.7)式中の核分裂中性子源項及び次のタイムステップにおける遅発中性子先行核密度
𝐶𝑚 (𝑥, 𝑡 + ∆𝑡)について、以下の(2.5.8)式及び(2.5.9)式のように記述する。(2.5.9)式の導出については後述する。
𝑃(𝑥, 𝑡) = ∑ 𝜈Σ𝑓,𝑔′ (𝑥, 𝑡)𝜙𝑔′ (𝑥, 𝑡) ,(2.5.8)𝑔′

𝐶𝑚 (𝑥, 𝑡 + ∆𝑡) = 𝜇𝑚 𝐶𝑚 (𝑥, 𝑡) + 𝜂𝑚 𝑃(𝑥, 𝑡) + 𝜉𝑚 𝑃(𝑥, 𝑡 + ∆𝑡) ,(2.5.9)(2.5.8)式及び(2.5.9)式を用いると、(2.5.7)式は(2.5.10)式に書き換えることが出来る。
′ (𝑥,−∇ ∙ 𝐷𝑔 (𝑥, 𝑡 + ∆𝑡) ∙ ∇𝜙𝑔 (𝑥, 𝑡 + ∆𝑡) + Σ𝑟,𝑔𝑡 + ∆𝑡)𝜙𝑔 (𝑥, 𝑡 + ∆𝑡)𝐺′= 𝜒𝑝,𝑔𝑃(𝑥, 𝑡 + ∆𝑡) + ∑ Σ𝑠,𝑔′ →𝑔 (𝑥, 𝑡 + ∆𝑡)𝜙𝑔′ (𝑥, 𝑡 + ∆𝑡)(2.5.10)𝑔′ ≠𝑔+ 𝑄𝑔 (𝑥, 𝑡) ,但し、(2.5.10)式中のそれぞれの変数は以下の(2.5.11)式、(2.5.12)式、(2.5.13)式で表される。
′ (𝑥,
Σ𝑟,𝑔
𝑡) = Σ𝑟,𝑔 (𝑥, 𝑡) +1,
𝜃𝑣𝑔 ∆𝑡′𝜒𝑝,𝑔= 𝜒𝑝,𝑔 (1 − 𝛽) + ∑ 𝜒𝑑,𝑚,𝑔 𝜆𝑚 𝜉𝑚 ,(2.5.11)(2.5.12)𝑚

1
𝑄𝑔 (𝑥, 𝑡) = ∑ 𝜒𝑑,𝑚,𝑔 𝜆𝑚 (𝜇𝑚 𝐶𝑚 (𝑥, 𝑡) + 𝜂𝑚 𝑃(𝑥, 𝑡)) + ( − 1) 𝑅𝑔 (𝑥, 𝑡)𝜃𝑚(2.5.13)𝜙𝑔 (𝑥, 𝑡)+,𝜃𝑣𝑔 ∆𝑡(2.5.13)式中には中性子束の空間微分項を含む変数である𝑅が含まれている。これを他の変数を用いて表記する。(2.5.10)式より以下の(2.5.14)式を得る。
′ (𝑥,−∇ ∙ 𝐷𝑔 (𝑥, 𝑡) ∙ ∇𝜙𝑔 (𝑥, 𝑡) + Σ𝑟,𝑔𝑡)𝜙𝑔 (𝑥, 𝑡)𝐺′= 𝜒𝑝,𝑔𝑃(𝑥, 𝑡) + ∑ Σ𝑠,𝑔′ →𝑔 (𝑥, 𝑡)𝜙𝑔′ (𝑥, 𝑡) + 𝑄𝑔 (𝑥, 𝑡 − ∆𝑡) ,
𝑔′ ≠𝑔

15

(2.5.14)(2.5.14)式を(2.5.4)式に代入することで、以下の(2.5.15)式を得ることが出来る。
𝑅𝑔 (𝑥, 𝑡) =𝜙𝑔 (𝑥, 𝑡)+ ∑ 𝜒𝑑,𝑚,𝑔 𝜆𝑚 (𝐶𝑚 (𝑥, 𝑡) − 𝜉𝑚 𝑃(𝑥, 𝑡)) − 𝑄𝑔 (𝑥, 𝑡 − ∆𝑡) ,𝜃𝑣𝑔 ∆𝑡(2.5.15)𝑚(2.5.15)式を用いることで、(2.5.10)式中にある中性子束の空間微分項を消去することが出来た。これにより、(2.5.10)式右辺第一項は核分裂中性子源項として取り扱うことが出来、また第三項は既知の値として外部中性子源として扱うことが出来る。これにより(2.5.10)式は核分裂中性子源を持つ固定源問題と同様の解法を用いることが出来ることが分かる。計算の流れとしては、ある時間の中性子束分布を内部反復及び外部反復により求め、算出された中性子束を用いて𝑅及び𝑄を計算し次の時間ステップの計算に用いるという繰り返しで、時間依存の中性子束を計算することが出来る。但し、初期状態として臨界定常状態を仮定する時、𝑅の初期値は以下のように0となる。
𝑅𝑔 (𝑥, 0) =1 𝜕𝜙𝑔 (𝑥, 0)
=0,
𝑣𝑔
𝜕𝑡

(2.5.16)

2.6 遅発中性子先行核密度時間変化の準解析的解法
次のタイムステップの遅発中性子先行核密度分布𝐶𝑚 (𝑡 + ∆𝑡)は、(2.5.6)式で表されると述べた。ここで、準解析的に(2.5.6)式を導出する方法を説明する。遅発中性子先行核密度は(2.5.2)式の方程式で表現される。
𝜕𝐶𝑚 (𝑥, 𝑡)= 𝛽𝑚 ∑ 𝜈Σ𝑓,𝑔′ (𝑥, 𝑡)𝜙𝑔′ (𝑥, 𝑡) – 𝜆𝑚 𝐶𝑚 (𝑥, 𝑡) ,𝜕𝑡′(再掲)(2.5.2)

𝑔(2.5.2)式の両辺にexp(𝜆𝑚 𝑡)を乗じる。
exp(𝜆𝑚 𝑡)

𝜕𝐶𝑚 (𝑥, 𝑡)= exp(𝜆𝑚 𝑡) 𝛽𝑚 𝑃(𝑥, 𝑡)– exp(𝜆𝑚 𝑡) 𝜆𝑚 𝐶𝑚 (𝑥, 𝑡) ,𝜕𝑡

(2.6.1)(2.6.1)式は(2.6.2)式のように変形することが出来る。
𝜕(exp(𝜆𝑚 𝑡) 𝐶𝑚 (𝑥, 𝑡)) = exp(𝜆𝑚 𝑡) 𝛽𝑚 𝑃(𝑥, 𝑡) ,𝜕𝑡

(2.6.2)(2.6.2)式の両辺を[𝑡, 𝑡 + ∆𝑡]で積分する。
𝑡+∆𝑡

∫
𝑡

𝑑𝑡 ′

𝑡+∆𝑡
𝜕
′)
′ ))
(𝑥,
(exp(𝜆𝑡𝐶𝑡=𝛽∫𝑑𝑡 ′ exp(𝜆𝑚 𝑡 ′ ) 𝑃(𝑥, 𝑡 ′ ) ,
𝑚
𝑚
𝑚
𝜕𝑡 ′
𝑡

(2.6.3)(2.6.3)式より、(2.6.4)式を得る。
𝐶𝑚 (𝑥, 𝑡 + ∆𝑡) = 𝐶𝑚 (𝑥, 𝑡) exp(−𝜆𝑚 ∆𝑡)𝑡+∆𝑡+ 𝛽𝑚 ∫

𝑑𝑡 ′ exp(−𝜆𝑚 (𝑡 + ∆𝑡 − 𝑡 ′ )) 𝑃(𝑥, 𝑡 ′ ) ,(2.6.4)𝑡

ここで、(2.6.4)式中の核分裂源項について、以下の線形近似を用いる。
𝑃(𝑥, 𝑡 ′ ) = 𝑃(𝑥, 𝑡) +1
{𝑃(𝑡 + ∆𝑡) − 𝑃(𝑡)}(𝑡 − 𝑡 ′ ) ,∆𝑡
16(2.6.5)(2.6.5)式を(2.6.4)式に代入し積分することで、以下の(2.6.6)式を得ることが出来る。
𝐶𝑚 (𝑥, 𝑡 + ∆𝑡) = 𝐶𝑚 (𝑥, 𝑡) exp(−𝜆𝑚 ∆𝑡)+𝛽𝑚 1 − exp(−𝜆𝑚 ∆𝑡){− exp(−𝜆𝑚 ∆𝑡)} 𝑃(𝑥, 𝑡)
𝜆𝑚
𝜆𝑚 ∆𝑡

+

𝛽𝑚
1 − exp(−𝜆𝑚 ∆𝑡){1 −} 𝑃(𝑥, 𝑡 + ∆𝑡) ,
𝜆𝑚
𝜆𝑚 ∆𝑡

(2.6.6)

但し、(2.6.6)式中の各変数は以下の(2.6.7)、(2.6.8)式、(2.6.9)式で表される。
𝜇𝑚 = exp(−𝜆𝑚 ∆𝑡) ,𝜂𝑚 =(2.6.7)𝛽𝑚 1 − exp(−𝜆𝑚 ∆𝑡){− exp(−𝜆𝑚 ∆𝑡)} ,
𝜆𝑚
𝜆𝑚 ∆𝑡

(2.6.8)𝛽𝑚
1 − exp(−𝜆𝑚 ∆𝑡)
{1 −
},
𝜆𝑚
𝜆𝑚 ∆𝑡

(2.6.9)𝜉𝑚 =(2.6.7)式、(2.6.8)式、(2.6.9)式の各パラメータは(2.5.9)式で用いられているパラメータ(𝜇𝑚 、
𝜂𝑚 、𝜉𝑚 )と同一のものである。また、初期状態として臨界定常状態を仮定したとき、遅発中性子先行核密度の初期値は以下の(2.6.10)式で表される。
𝐶𝑚 (𝑥, 0) =𝛽𝑚 𝑃(𝑥, 0),𝜆𝑚(2.6.10)

2.7 EFP モデルの適用
ここまで記述してきた計算手法は、全て従来の 6 群モデルを適用した場合のものである。
では、中性子拡散動特性方程式に、EFP モデルを適用した場合の計算手法について考える。
EFP モデルを適用した中性子拡散動特性方程式は以下の(2.7.1)式及び(2.7.2)式で表される。
1 𝜕𝜙𝑔 (𝑥, 𝑡)= ∇ ∙ 𝐷𝑔 (𝑥, 𝑡) ∙ ∇𝜙𝑔 (𝑥, 𝑡) − Σ𝑟,𝑔 (𝑥, 𝑡)𝜙𝑔 (𝑥, 𝑡)𝑣𝑔
𝜕𝑡𝐺+ ∑ Σ𝑠,𝑔′ →𝑔 (𝑥, 𝑡)𝜙𝑔′ (𝑥, 𝑡)𝑔′ ≠𝑔(2.7.1)+ 𝜒𝑝,𝑔 (1 − 𝛽) ∑ 𝜈Σ𝑓,𝑔′ (𝑥, 𝑡)𝜙𝑔′ (𝑥, 𝑡)𝑔′+ ∑ 𝜒𝑑,𝑚,𝑔 𝜆𝑚 𝐶𝑚 (𝑥, 𝑡)𝑃𝑛𝑚 ,𝑚𝜕𝐶𝑚 (𝑥, 𝑡)= 𝑦𝑚 ∑ Σ𝑓,𝑔′ (𝑥, 𝑡)𝜙𝑔′ (𝑥, 𝑡) – 𝜆𝑚 𝐶𝑚 (𝑥, 𝑡)𝜕𝑡′𝑔(2.7.2)+ ∑ 𝜆𝑚′ 𝐶𝑚′ (𝑥, 𝑡)𝑃𝑚′ →𝑚 ,𝑚′

一般的な 6 群モデルでは核分裂により発生した核種は、全て崩壊により中性子を放出し、中17性子を放出しない他の核種になると仮定されている。そのため、(2.5.1)式及び(2.5.2)式中の変数𝐶は遅発中性子先行核密度を表している。一方、EFP モデルでは崩壊により中性子を放出しない核種についても陽に取り扱う。そのため、変数𝐶は正確には核分裂生成物核種の数
密度を表すことになる。そのため、今後 EFP モデルで記述する際には変数𝐶は核分裂生成物(FP)数密度と呼称する。また、𝑃𝑛𝑚 は核種の崩壊により遅発中性子が発生する確率を表しており、また𝑃𝑚′ →𝑚 は核種𝑚′の崩壊により核種𝑚が生成する確率を表している。(2.7.2)式右辺
第三項は核種の崩壊による他核種の生成を表す項となっている。この項には FP 数密度の時間変化を示す𝐶が含まれているため、これらの連立微分方程式は互いが独立ではない。そのため、FP 数密度に関する(2.7.2)式の解法について、従来の 6 群モデルと同様の準解析的な手法を用いることが出来ない。そこで、本研究では(2.7.1)式及び(2.7.2)式を行列形式で表記
し、離散化及び数値解を得る方法を提案する。

2.8 部分行列指数法(2.7.1)式及び(2.7.2)式を行列形式で表記すると、(2.8.1)式で表される。
𝑑𝒏(𝑡)= 𝑨𝒏(𝑡),𝑑𝑡(2.8.1)

ここで、𝒏は各メッシュの中性子束及び𝑀核種の FP 数密度を要素として持つベクトルである。𝒏の形式を(2.8.2)式に示す。
𝒏(𝑡)= (𝜙1(1)

(1)

… 𝜙𝐺

(1)

𝐶1

(1)

… 𝐶𝑀

(2)

𝜙1

(2)

… 𝐶𝑀

(𝑁)

… 𝜙𝑔

𝑇

(𝑁)
… 𝐶𝑀 ) ,(2.8.2)(2.8.2)式中の添え字については、Fig. 2.8.1 で表される。

Fig. 2.8.1 (2.8.2)式中の添え字の意味
また、𝑨は正方係数行列である。二次元体系の場合、𝑨は Fig. 2.8.2 で示されるような対角区間とその隣接する区間のみ値を持つ疎行列である。なお、Fig. 2.8.2 の網掛け部分と、×部分
は𝑨が値を持つ部分を表している。

18

Fig. 2.8.2 行列𝑨の構成(2.8.1)式のような行列形式の微分方程式は、一般的には形式的に(2.8.3)式のように解を得ることが出来る。
𝒏(𝑡) = exp(𝑨𝑡) 𝒏(0),(2.8.3)

ここで、exp(𝑨𝑡)は行列指数関数と呼ばれ、以下の(2.8.4)式で定義される。
1
exp(𝑨𝑡) = 𝑰 + 𝑨𝑡 + (𝑨𝑡)2 + ⋯2(2.8.4)(2.8.3)式は行列指数関数をなんらかの方法で近似することで、数値解を得ることが出来る。
しかし、空間依存動特性方程式の場合、係数行列𝑨に対して行列指数を精度良く求めることは困難である。そこで、本研究では行列𝑨を分割し部分的に行列指数法を適用する手法と、𝜃法により離散化する手法の二種類を提案する。
最初に、(2.8.1)式のあるメッシュに対応した部分を取り出したものを(2.8.5)式のように分
割する。
𝑑 𝒏1
𝑨
(𝒏 ) = ( 11
𝑨
𝑑𝑡 2
21

𝑨12 𝒏1)( ),𝑨22 𝒏2(2.8.5)

これらは連立形式では(2.8.6)式及び(2.8.7)式で示される。
𝑑𝒏1 (𝑡)= 𝑨11 𝒏1 (𝑡) + 𝑨12 𝒏2 (𝑡) ,𝑑𝑡(2.8.6)𝑑𝒏2 (𝑡)= 𝑨21 𝒏1 (𝑡) + 𝑨22 𝒏2 (𝑡) ,𝑑𝑡(2.8.7)

ここで、分割によって行列𝑨22 について行列指数を計算することが可能になったと仮定する。
この場合、部分的に行列指数法を適用する手法を説明する。 (2.8.7)式の両辺に左からexp(−𝑨22 𝑡)を乗じる。
exp(−𝑨22 𝑡)

𝑑𝒏2 (𝑡)= exp(−𝑨22 𝑡) 𝑨21 𝒏1 (𝑡) + exp(−𝑨22 𝑡) 𝑨22 𝒏2 (𝑡) ,𝑑𝑡(2.8.8)行列指数関数の定義式((2.8.4)式)を用いて(2.8.8)式を変形すると、(2.8.9)式を得る。
𝑑(exp(−𝑨22 𝑡) 𝒏2 (𝑡)) = exp(−𝑨22 𝑡) 𝑨21 𝒏1 (𝑡) ,𝑑𝑡(2.8.9)式を[𝑡, 𝑡 + ∆𝑡]で積分すると(2.8.10)式を得る。
19(2.8.9)exp(−𝑨22 (𝑡 + ∆𝑡)) 𝒏2 (𝑡 + ∆𝑡) − exp(−𝑨22 𝑡) 𝒏2 (𝑡)𝑡+∆𝑡exp(−𝑨22 𝑡 ′ )𝑨21 𝒏1 (𝑡 ′ ) 𝑑𝑡 ′ ,=∫(2.8.10)𝑡(2.8.10)式に左からexp(𝑨22 (𝑡 + ∆𝑡))を乗じ、整理すると以下の(2.8.11)式を得る。
𝑡+∆𝑡𝒏2 (𝑡 + ∆𝑡) = exp(𝑨22 ∆𝑡) 𝒏2 (𝑡) + ∫exp(𝑨22 (𝑡 − 𝑡 ′ + ∆𝑡))𝑨21 𝒏1 (𝑡 ′ ) 𝑑𝑡 ′ ,(2.8.11)𝑡

ここで、(2.8.11)式中の𝒏1 (𝑡 ′ )に対して(2.6.5)式と同様の線形近似を適用し、積分をすること
で、以下の(2.8.12)式を得ることが出来る。
𝒏2 (𝑡 + ∆𝑡) = 𝑴𝒏2 (𝑡) + 𝑬𝒏1 (𝑡) + 𝑿𝒏1 (𝑡 + ∆𝑡) ,(2.8.12)

但し、(2.8.12)式中の𝑴、𝑬、𝑿はそれぞれ(2.8.13)式、(2.8.14)式、(2.8.15)式で表される変数である。
𝑴 = exp(𝑨22 ∆𝑡) ,𝑬 = 𝑨−122 {(2.8.13)𝑨−1
22(𝑰 − 𝑴) − 𝑴} 𝑨21 ,∆𝑡(2.8.14)𝑨−1
22(𝑰 − 𝑴)} 𝑨21 ,∆𝑡(2.8.15)

𝑿 = 𝑨−1
22 {𝑰 −

また、(2.8.6)式に𝜃法を適用すると、以下の(2.8.16)式を得る。
𝒏1 (𝑡 + ∆𝑡) − 𝒏1 (𝑡)= {𝑨11 𝒏1 (𝑡 + ∆𝑡) + 𝑨12 𝒏2 (𝑡 + ∆𝑡)}𝜃
∆𝑡(2.8.16)+ {𝑨11 𝒏1 (𝑡) + 𝑨12 𝒏2 (𝑡)}(1 − 𝜃) ,(2.8.16)式を𝒏1 (𝑡 + ∆𝑡)について整理すると、(2.8.17)式を得ることが出来る。
𝑵0 𝒏1 (𝑡 + ∆𝑡) = 𝑵1 𝒏1 (𝑡) + 𝑵2 𝒏2 (𝑡) ,(2.8.17)

但し、(2.8.17)式中の𝑵0、𝑵1、𝑵2はそれぞれ(2.8.18)式、(2.8.19)式、(2.8.20)式で表される変数
である。
𝑵0 = 𝑰 − (𝑨11 + 𝑨12 𝑿)𝜃∆𝑡 ,(2.8.18)

𝑵1 = 𝑰 + 𝑨12 𝑬𝜃∆𝑡 + 𝑨11 (1 − 𝜃)∆𝑡 ,(2.8.19)

𝑵2 = 𝑨12 𝑴𝜃∆𝑡 + 𝑨12 (1 − 𝜃)∆𝑡 ,(2.8.20)(2.8.17)式、(2.8.18)式、(2.8.19)式、(2.8.20)式を用いることで、𝒏1 (𝑡)及び𝒏2 (𝑡)が既知であれば、𝒏1 (𝑡 + ∆𝑡)が計算できることが分かる。また、算出された𝒏1 (𝑡 + ∆𝑡)と(2.8.12)式を用いる
ことで、𝒏2 (𝑡 + ∆𝑡)も算出することが出来る。このようにして、部分的に行列指数法を適用
し、行列形式の微分方程式を解く手法を部分行列指数法(Partially Matrix Exponential method)と呼ぶこととする。

2.9 𝜃法による離散化
前節では(2.8.6)式及び(2.8.7)式を、部分行列指数法を用いることで数値的に解いた。本節では(2.8.7)式を𝜃法により離散化する方法について述べる。(2.8.7)式に𝜃法を適用すると20

(2.9.1)式を得る。
𝒏2 (𝑡 + ∆𝑡) − 𝒏2 (𝑡) = 𝑨21 (𝜃∆ 𝑡𝒏1 (𝑡 + ∆𝑡) + (1 − 𝜃)∆𝑡𝒏1 (𝑡))+ 𝑨22 (𝜃∆ 𝑡𝒏2 (𝑡 + ∆𝑡) + (1 − 𝜃)∆𝑡𝒏2 (𝑡)),(2.9.1)(2.9.1)式を整理すると、以下の(2.9.2)式を得る。
(𝑰 − 𝜃∆𝑡𝑨22 )𝒏2 (𝑡 + ∆𝑡) = (𝑰 + (1 − 𝜃)∆𝑡𝑨22 )𝒏2 (𝑡)+ (1 − 𝜃)∆𝑡𝑨21 𝒏1 (𝑡) + 𝜃∆ 𝑡𝑨21 𝒏1 (𝑡 + ∆𝑡),(2.9.2)(2.9.2)式より、𝒏2 (𝑡 + ∆𝑡)についての(2.9.3)式を得ることが出来る。
𝒏2 (𝑡 + ∆𝑡) = 𝑴𝜃 𝒏2 (𝑡) + 𝑬𝜃 𝒏1 (𝑡) + 𝑿𝜃 𝒏1 (𝑡 + ∆𝑡),(2.9.3)

但し、(2.9.3)式中の各変数は以下の(2.9.4)式、(2.9.5)式、(2.9.6)式で定義される。
𝑴𝜃 = (𝑰 − 𝜃∆𝑡𝑨22 )−1 (𝑰 + (1 − 𝜃)∆𝑡𝑨22 ) ,(2.9.4)𝑬𝜃 = (1 − 𝜃)∆𝑡(𝑰 − 𝜃∆𝑡𝑨22 )−1 𝑨21 ,(2.9.5)𝑿𝜃 = 𝜃∆ 𝑡(𝑰 − 𝜃∆𝑡𝑨22

)−1

𝑨21 ,

(2.9.6)2.10 空間依存動特性方程式への部分行列指数法の適用
空間依存動特性方程式へ部分行列指数法を適用する。空間依存動特性方程式を(1.7.5)式のように行列表記し分割した場合、𝒏1は中性子束に対応し、𝒏2は FP 数密度に対応する。この場合、𝑨22 は FP の崩壊と他の核種への移行を表す行列となる。これは燃焼行列とも呼ばれる行列で、行列指数関数が計算可能である。分割された行列を用いて、部分行列指数法を適用した場合、FP 数密度分布は以下の(2.10.1)式で表される。
𝑪(𝑡 + ∆𝑡) = 𝝁𝑪(𝑡) + 𝑃(𝑡)𝜼 + 𝑃(𝑡 + ∆𝑡)𝝃 ,(2.10.1)

但し、(2.10.1)式中の各パラメータは(2.10.2)式、(2.10.3)式、(2.10.4)式、(2.10.5)式及び(2.10.6)式で定義される。
𝝁 = exp(−𝝀∆𝑡),(2.10.2)1
𝜼 = 𝝀−1 { 𝝀−1(𝑰 − 𝝁) − 𝝁} 𝒚,∆𝑡(2.10.3)

1 −1
𝝀 (𝑰 − 𝝁)} 𝒚,∆𝑡(2.10.4)

𝝃 = 𝝀−1 {𝑰 −𝒚 = (𝑦1 , 𝑦2 , … , 𝑦𝑀 )𝑇 ,𝜆1
⋯ −𝑃𝑀→1 𝜆𝑀
⋮
⋱
⋮
𝝀=(
)
−𝑃1→𝑀 𝜆1 ⋯𝜆𝑀

(2.10.5)(2.10.6)

2.11 行列指数の近似
前節では、部分行列指数法を空間依存動特性方程式へと適用させた。ここで、(2.10.1)式で表される次のタイムステップでの FP 核種数密度を算出するためには、行列指数を数値的に求める必要がある。行列指数の近似法には CRAM や MMPA 法などの手法がある[4]が、本研究では Padé 近似法を用いている。Padé 近似法では、指数関数exp(𝑥)は以下の式で近似す21

ることが出来る。
exp(𝑥) ≈𝑁𝑝𝑞 (𝑥),𝐷𝑝𝑞 (𝑥)(2.11.1)

但し、(2.11.1)式中の各変数は以下の(2.11.2)式、(2.11.3)式で表される。
𝑝(𝑝 + 𝑞 − 𝑘)! 𝑝!
𝑥𝑘,(𝑝 + 𝑞)! 𝑘! (𝑝 − 𝑘)!

(2.11.2)(𝑝 + 𝑞 − 𝑘)! 𝑞!
(−𝑥)𝑘 ,
(𝑝 + 𝑞)! 𝑘! (𝑞 − 𝑘)!

(2.11.3)

𝑁𝑝𝑞 (𝑥) = ∑𝑘=0𝑞

𝐷𝑝𝑞 (𝑥) = ∑𝑘=0

ここで、(2.11.2)式及び(2.11.3)式中の𝑝、𝑞は任意の整数であり、Padé 近似の次数である。一
般的には𝑝 = 𝑞となる値が用いられる。(2.11.1)式はスカラーの指数に対する式であるが、同様に行列指数に対しても適用させることが出来る。行列指数exp(𝑨∆𝑡)に対して Padé 近似を適用した場合、以下の(2.11.4)式で近似することが出来る。
−1 (𝑨
exp(𝑨∆𝑡) ≈ 𝐷𝑝𝑞∆𝑡) ∙ 𝑁𝑝𝑞 (𝑨 ∆𝑡) ,(2.11.4)

但し、(2.11.4)式中の各変数は以下の(2.11.5)式及び(2.11.6)式で表される。
𝑝(𝑝 + 𝑞 − 𝑘)! 𝑝!
(𝑨∆𝑡)𝑘 ,(𝑝 + 𝑞)! 𝑘! (𝑝 − 𝑘)!

(2.11.5)(𝑝 + 𝑞 − 𝑘)! 𝑞!
(−𝑨∆𝑡)𝑘 ,(𝑝 + 𝑞)! 𝑘! (𝑞 − 𝑘)!

(2.11.6)

𝑁𝑝𝑞 (𝑥) = ∑𝑘=0𝑞

𝐷𝑝𝑞 (𝑥) = ∑𝑘=0

本研究では特に断りが無い場合、6 次の Padé 近似(𝑝 = 𝑞 = 6)を用いて行列指数の近似を行っている。

2.12 加速法
本節では、2.2 節及び 2.3 節で述べた内部反復と外部反復の収束を速めるための手法について記述する。

2.12.1 Gauss-Seidel 法
内部反復では(2.2.3)式を用いて中性子束の更新を収束まで行う。この時、𝑖 = 1,2,3, …
(𝑘+1)という順番で計算を行うとすると、𝑖番目のメッシュの中性子束(𝜙𝑔,𝑖)を更新する際には(𝑖 − 1)番目のメッシュの中性子束は既に更新済みである。従って、反復計算に以下の(2.12.1)式を用いることが出来る。
(𝑘+1)

𝑥+
𝑥−
(𝐶𝑔,𝑖
+ 𝐶𝑔,𝑖+ 𝛴𝑟,𝑔,𝑖 ∆𝑥𝑖 )𝜙𝑔,𝑖

(𝑘)

(𝑘+1)

𝑥+
𝑥−
= 𝐶𝑔,𝑖𝜙𝑔,𝑖+1 + 𝐶𝑔,𝑖𝜙𝑔,𝑖−1 + 𝑆𝑔,𝑖 ∆𝑥𝑖 ,(2.12.1)

この反復解法は Gauss-Seidel 法と呼ばれ、一般的に収束が速い。

2.12.2 Successive Over Relaxation 法
Successive Over Relaxation(SOR)法は、外挿を用いた加速法である。SOR 法では中性子束の更新を以下の(2.12.2)式を用いて行う。

22

𝜙 (𝑘+1) = 𝜙 (𝑘) + 𝜔 (𝜙 ′(𝑘+1)

− 𝜙 (𝑘) ) ,(2.12.2)

但し、(2.12.2)式中の各変数の意味は以下の通りである。
𝜙 (𝑘+1)：SOR 法適用後の中性子束𝜙 (𝑘)：𝑘回目の反復によって得られた中性子束𝜙𝜔′ (𝑘+1)：𝑘 + 1回目の反復によって得られた SOR 法適用前の中性子束：加速パラメータ(𝜔 = 1で通常の反復計算と等価)

この時、加速パラメータの値域は1 < 𝜔 < 2である。SOR 法では、Fig. 2.12.1 で示されるよ
うに通常の反復計算で得られた更新値を加速パラメータによって強調することで収束を速める手法である。

Fig. 2.12.1 SOR 法のイメージ

2.12.3 Simplified Fixed Source Scaling Factor 法
Simplified Fixed Source Scaling Factor 法[[6]は核分裂源を含む固定中性子源問題を解く際に有効な加速法である。まず、正確な中性子束の値(𝜙𝑒𝑥 )が既知の場合、以下の(2.12.3)式が成立する。
(2.12.3)

𝐋𝜙𝑒𝑥 = 𝐌𝜙𝑒𝑥 + 𝐒 ,但し、𝐋、𝐌、𝐒はそれぞれ中性子の消失、中性子の生成、外部中性子源を表す演算子である。
しかし、(2.12.3)式は中性子束が十分に収束していない場合、成立しない。そのため、Simplified
Fixed Source Scaling Factor 法では、𝑘回目の反復での中性子束を用いた補正係数が用いられる。この補正係数は(2.12.4)式で表される。
𝑓=〈𝐒〉〈𝐋𝜙 (𝑘) 〉 − 〈𝐌𝜙 (𝑘) 〉,(2.12.4)

但し、(2.12.4)式中の〈 〉は空間、エネルギーの積分を表す演算子である。Simplified Fixed
Source Scaling Factor 法では、補正係数を計算し、式のように𝑘回目の反復に適用させること
で、中性子束の収束を速めることが可能となる。
𝜙̃ (𝑘) = 𝑓𝜙 (𝑘) ,(2.12.5)

2.13 EFP モデルで使用している核データ
EFP モデルでは各々の核種を陽に取り扱うため、核種が持つ定数を評価済み核データラ
23イブラリから取得している。各核種が持つ定数として、核分裂収率、半減期、崩壊分岐比、崩壊により放出される粒子、崩壊後の原子核の情報などがある。このうち、核分裂収率については JENDL/FPY-2011 から取得しており、その他の情報については JENDL/FPD-2011 から取得している。この核データライブラリでは、崩壊による中性子放出確率𝑃𝑛については直接的に定義されていない。そのため、𝑃𝑛を以下の(2.13.1)式で定義する。
𝑐𝑖

𝑃𝑛𝑖 = ∑(𝑅𝑗 ∙ 𝑁𝑗 ) ,(2.13.1)𝑗

但し、𝑐、𝑅、𝑁はそれぞれ崩壊チャンネルの総数、崩壊チャンネルに対する分岐比及び崩壊により放出される中性子数を表す。また、核分裂収率には独立核分裂収率と累積核分裂収率の 2 種類が存在するが、本研究では全ての核種の崩壊による他核種への移行を陽に取り扱うため、独立核分裂収率を用いている。

2.14 EFP モデルによるガス状 FP の漏洩計算
EFP モデルでは、事故時に炉心からガス状 FP が漏洩する事象を取り扱うことが出来る。
本節ではこの事象を再現する際に必要な計算について述べる。
まずガス状 FP の漏洩を再現するため、漏洩する速さを定義する。本研究では単位時間当たりに漏洩する割合を漏洩定数𝑙[/s]として定義した。また、(2.14.1)式のように定義すること
で、漏洩により核種数密度が半減する時間(漏洩半減期𝐿1/2 [𝑠])としても定義することが出来る。
𝐿1/2 =

ln 2,
𝑙(2.14.1)

実際の動特性計算時には、まず(2.10.1)式を用いて次タイムステップの FP 核種数密度を算出する。その後、漏洩を仮定した核種𝑚′及び漏洩を仮定した領域のメッシュ𝑖 ′ 、については以′下の式で表される計算を行い、漏洩を含めた次タイムステップの FP 核種数密度𝐶𝑚′ ,𝑖 ′ (𝑡 +∆𝑡)を算出する。
′
𝐶𝑚′ ,𝑖 ′ (𝑡 + ∆𝑡) = 𝐶𝑚′ ,𝑖 ′ (𝑡 + ∆𝑡)(1 − exp(𝑙∆𝑡)) ,(2.14.2)

但し、(2.14.2)式中の𝐶𝑚′ ,𝑖 ′ (𝑡 + ∆𝑡)は、(2.10.1)式を用いて算出した FP 核種数密度である。

2.15 本章のまとめ
本章ではまず定常問題における拡散動特性方程式の解法について、固定源問題(2.1 節)から固有値問題(2.3 節)へと拡張し述べた。また、2.4 節において体系を二次元に拡張した場合の計算について記述した。
2.5 節以降では本研究のメインテーマである動特性計算に着目し、その計算手法について述べた。2.6 節までは従来一般的に用いられている数値解法について、
また 2.7 節以降では本研究で提案した EFP モデルで用いられる数値解法について記述した。
特に 2.8 節では EFP モデルによる空間依存動特性解析へ行列指数法を適用させるために新24

たに開発した PME 法について詳しく説明した。次章では、本章で記述した計算手法を基に実装した動特性解析コードを用いて行った解析の結果を示す。

25

第3章 結果
3.1 部分行列指数法と𝜃法による一次元動特性解析
前章では空間依存動特性解析の計算手法と、EFP モデルを適用させるための計算手法について述べた。この EFP モデルを適用させるための手法として、部分行列指数(PME)法及び𝜃法の 2 種類を紹介した。本節では一次元動特性方程式へ EFP モデルを適用させ、2 種類の計算手法の比較を行った。解析対象として、一次元軽水反射体付きの体系を設定した。体系図を Fig. 3.1.1 に示す。また、各定数について Table 3.1.1 に示す。

Fig. 3.1.1 一次元軽水反射体付き体系体系図
Table 3.1.1 各領域の定数

RegionEnergyGroupΣ𝑠 [/cm]𝐷[cm]

Σ𝑎 [/cm]𝜈𝛴𝑓 [/cm]
𝑔→1

𝑔→2

1

1.58

0.0032000.0178

2

0.271

0.093

0.124

0

0

1

1.41

0

0

0

0.04762

0.117

0.0191

0

0

0

Fuel

Reflector

𝜈𝑑𝑣1 [cm/s]

𝑣2 [cm/s]

𝛽0.01691.00 × 107

2.00 × 105

0.007526

二種類の手法を比較するため、各々のタイムステップを∆𝑡 = 0.1[s]とした。また、Reference
として、PME 法で∆𝑡 = 0.001[s]とした結果を用いた。Fig. 3.1.2 にそれぞれの手法による出力の変化を示す。

Fig. 3.1.2

出力の変化

また、𝑡 = 0[s]近傍の結果を Fig. 3.1.3 に示す。

27

Fig. 3.1.3

出力の変化(𝑡 = 0[s]近傍)

Fig. 3.1.2 及び Fig. 3.1.3 から、PME 法、𝜃法のどちらも𝑡 = 0[s]近傍で振動していることが分かる。しかし、∆𝑡 = 0.1[s]とした時の手法による差はほとんど見られず、最大となる点
でも 0.042%だった。Fig. 3.1.2 及び Fig. 3.1.3 では、PME 法と𝜃法のいずれも全体的に過小評価の傾向が見られる。そこで、Fig. 3.1.4 には∆𝑡 = 0.01[s]の条件で PME 法及び𝜃法で計算した結果を示す。

28

Fig. 3.1.4 出力の変化(∆𝑡 = 0.01[s])Fig. 3.1.4 より、即発跳躍領域でわずかに差異が見られるものの、それ以降の時間領域での過小評価が改善されていることが分かる。また、即発跳躍領域での振動も見られなくなったことが分かる。
次に、この時の核種数密度の変化を比較する。比較する核種として、非常に短い半減期を持つ FP である Co-78(𝑇1/2 = 4.84[ms])、比較的短い半減期を持つ FP である Rh121(𝑇1/2 = 0.121[s])及び長い半減期を持つ FP である Br-87(𝑇1/2 = 55.7[s])を選出した。Fig.
3.1.5 に Co-78 の核種数密度の変化を、また Fig. 3.1.6 にその𝑡 = 0[s]近傍の結果を示す。
Fig. 3.1.7 及び Fig. 3.1.8 には同様に Rh-121 の結果を示し、Fig. 3.1.9、Fig. 3.1.10 には Br-87の結果を示す。いずれの図も、PME 法、𝜃法、Reference の 3 つの結果をプロットしている。

29

Fig. 3.1.5 Co-78 炉内平均数密度変化

Fig. 3.1.6 Co-78 炉内平均数密度変化(𝑡 = 0[s]近傍)

30

Fig. 3.1.7 Rh-121 炉内平均数密度変化

Fig. 3.1.8 Rh-121 炉内平均数密度変化(𝑡 = 0[s]近傍)

31

Fig. 3.1.9 Br-87 炉内平均数密度変化

Fig. 3.1.10 Br-87 炉内平均数密度変化(𝑡 = 0[s]近傍)

32

Fig. 3.1.5 及び Fig. 3.1.6 から、Co-78 の核種数密度について PME 法、𝜃法のどちらも最初の数ステップで結果が振動していることが分かる。Reference との相対差は、どちらの手法
も 1 ステップ目で最大となり、PME 法で-10.86%、𝜃法で-12.50%となった。どちらの手法も∆𝑡 = 0.1[s]という厳しい条件では、十分な精度を得ることが出来なかった。Fig. 3.1.7 及び Fig.
3.1.8 は比較的短い半減期を持つ Rh-121 の結果であるが、この核種についても 1 ステップ目
で Reference との相対差が最大となり、PME 法で-11.08%、𝜃法で-11.92%となった。この 2 核種については、即発跳躍部分で Reference との差が大きくなるという結果となった。一方、
長い半減期を持つ Br-87 については、即発跳躍部分でも Reference との相対差は PME 法、𝜃
法ともに 0.1%以下となった。また、この核種については二つの手法間の差も10−3%以下となり、他の二核種の結果と異なる傾向が見られた。次節ではこの結果について、行列指数の精度の点から考察を行う。
加えて、いずれの核種の結果でも出力と同様に PME 法及び𝜃法が過小評価する傾向が見られた。∆𝑡 = 0.01[s]で計算した際の Co-78 数密度の結果を Fig. 3.1.11 に示す。

Fig. 3.1.11 Co-78 炉内平均数密度変化(∆𝑡 = 0.01[s])Fig. 3.1.11 より、出力の結果と同様即発跳躍領域で過小評価する傾向は見られるが、それ以降の時間領域での過小評価が改善した。また、𝑡 = 0.0[s]近傍での振動も見られなくなった
ことが分かる。

33

3.2 部分行列指数法と𝜃法の計算精度
2.8 節で述べたように、部分行列指数法では行列指数を用いて FP 核種数密度の算出を行っている((2.8.12)式、(2.8.13)式)。
𝒏2 (𝑡 + ∆𝑡) = 𝑴𝒏2 (𝑡) + 𝑬𝒏1 (𝑡) + 𝑿𝒏1 (𝑡 + ∆𝑡) ,(再掲)(2.8.12)

𝑴 = exp(𝑨22 ∆𝑡) ,(再掲)(2.8.13)

また、行列指数は Padé 近似を用いて近似を行っている((2.11.4)式、(2.11.5)式、(2.11.6)式)。
−1 (𝑨
exp(𝑨∆𝑡) ≈ 𝐷𝑝𝑞∆𝑡) ∙ 𝑁𝑝𝑞 (𝑨 ∆𝑡) ,𝑝(𝑝 + 𝑞 − 𝑘)! 𝑝!
(𝑨∆𝑡)𝑘 ,(𝑝 + 𝑞)! 𝑘! (𝑝 − 𝑘)!

(再掲)(2.11.5)(𝑝 + 𝑞 − 𝑘)! 𝑞!
(−𝑨∆𝑡)𝑘 ,(𝑝 + 𝑞)! 𝑘! (𝑞 − 𝑘)!

(再掲)(2.11.6)

𝑁𝑝𝑞 (𝑥) = ∑𝑘=0𝑞

𝐷𝑝𝑞 (𝑥) = ∑𝑘=0(再掲)(2.11.4)

ここで、Padé 近似の次数を 1 次とした場合(𝑝 = 𝑞 = 1)、以下の式を得る。
𝑴 = exp(𝑨22 ∆𝑡) = (𝑰 −

−1
∆𝑡
∆𝑡
𝑨22 ) (𝑰 + 𝑨22 ) ,22(3.2.1)

また、部分行列に分割した微分方程式に𝜃法を適用した場合は、前述の(2.9.3)式、(2.9.4)式で表される。
𝒏2 (𝑡 + ∆𝑡) = 𝑴𝜃 𝒏2 (𝑡) + 𝑬𝜃 𝒏1 (𝑡) + 𝑿𝜃 𝒏1 (𝑡 + ∆𝑡),

(再掲)(2.9.3))−1 (𝑰(再掲)(2.9.4)

𝑴𝜃 = (𝑰 − 𝜃∆𝑡𝑨22+ (1 − 𝜃)∆𝑡𝑨22 ) ,ここで、(3.2.1)式は(2.9.4)式で𝜃 = 0.5とした場合(Crank Nicolson 法)と等価であることが分かる。つまり、部分行列指数法と比較した場合、行列指数に関して、𝜃法は一次の Padé 近似相
当の精度しか持たないと言い換えることが出来る。そのため、𝑨の対角成分が非常に大きくなる Co-78 などの超短半減期核種を計算する場合、二次以降の項の影響が大きくなるため、𝜃法での近似精度が悪くなったと考えられる。一方、𝑨の対角成分が小さくなる Br-87 などの長半減期核種については、二次以降の項が非常に小さくなるため、𝜃法でもほとんど差が無い結果となったと考えられる。加えて、PME 法と𝜃法の近似精度及び、PME 法における行列指数計算精度が結果に与える影響を検証した。PME 法で Padé 近似の次数を1とした場
合((2.11.5)式、(2.11.6)式において、𝑝 = 𝑞 = 1)と𝜃法で Co-78 の炉内平均数密度を計算した結果を Fig. 3.2.1 に示す。なお、どちらも∆𝑡 = 0.1[s]の条件で計算を行った。

34

Fig. 3.2.1 Co-78 炉内平均数密度変化(𝑝 = 𝑞 = 1)
Fig. 3.2.1 より、一次の Padé 近似を用いた PME 法の結果が、𝜃法の結果と非常に近いことが分かる。このことから、𝜃法が一次の近似と等価であることを確認することが出来た。また、PME 法における近似精度についても、行列指数の計算精度が支配的であることが分かった。
前節と本節の検証より、粗いタイムステップを用いた PME 法と𝜃法は、どちらも Referenceに対する近似精度は悪いものの、短半減期核種については核種数密度の算出に差が出ることが分かった。しかし、Fig. 3.1.2 及び Fig. 3.1.3 から分かるように、出力の結果に対しては二手法間にほとんど差が見られなかった。このことから、Co-78 に代表されるような、超短半減期核種は出力に対してほとんど影響を持たないということが分かった。

3.3 TWIGL ベンチマーク問題
本節からは体系を二次元に拡張し、動特性解析を行う。解析体系として、TWIGL ベンチ
マーク問題を設定した。体系図及び、各定数を Fig. 3.3.1、Table 3.3.1 に示す。

35

Fig. 3.3.1 TWIGL ベンチマーク問題体系図
Table 3.3.1 定数の定義

Material

EnergyGroupΣ𝑠 [/cm]𝐷[cm]

Σ𝑎 [/cm]νΣf [/cm]

𝜒
𝑔→1

𝑔→2

1

1.4

0.010

0.007

1.0

0.0

0.01

2

0.4

0.150

0.200

0.0

0.0

0

1

1.4

0.010

0.007

1.0

0.0

0.01

2

0.4

0.150

0.200

0.0

0.0

0.0

1

1.3

0.008

0.003

1.0

0.0

0.01

2

0.5

0.050

0.060

0.0

0.0

0.0

1

2

3

𝜈

𝑣1 [cm/s]

𝑣2 [cm/s]

𝛽2.43

1.00 × 107

2.00 × 105

0.0075

この体系に𝑡 = 0.0[s]でステップ状反応度を印加する場合、Material 1 の吸収断面積を以下36

の(3.3.1)式のように変化させる。
Σa,2 (𝑡) = Σ𝑎,2 (0) − 0.0035(3.3.1)(𝑡 ≥ 0) ,Fig. 3.3.1 で、Material 1 及び Material 2 はいずれも燃料領域を構成しており、また Material 3はブランケット領域を構成している。Material 1 と Material 2 は同じものであるが、Material
1 は摂動対象である。体系の西側、及び南側は完全反射境界条件となっている。そのため、
TWIGL ベンチマーク問題は、燃料領域とブランケット領域から成る1/4炉心体系であるということが出来る。また、オリジナルの TWIGL ベンチマーク問題では、遅発中性子先行核を 1 群(𝜆 = 0.08[/s])で取り扱っている。Fig. 3.3.2 にステップ状反応度印加後の炉心相対出力の変化を示す。Fig. 3.3.2 には Reference(先行核１群)、EFP モデル、6 群モデルの 3 ケースの結果を示す。但し、この 6 群モデルはフィッティングにより EFP モデルと整合を持った
6 群定数を使用しており、Table 1.2.1 の Keepin の 6 群モデルとは異なる[8]。使用した定数を Table 3.3.2 に示す。

Table 3.3.2 EFP モデルに整合を持つ 6 群モデル

先行核生成割合グループ

半減期T1⁄2 [s]崩壊定数𝜆𝑖 [s−1 ]

1

47.8

0.01450.05532

19.7

0.0351

0.191

3

4.65

0.149

0.215

4

1.91

0.363

0.363

5

0.499

1.39

0.119

6

0.169

4.09

0.058737𝑎𝑖Fig. 3.3.2 炉心相対出力の変化
Reference とその他 2 ケースでは、遅発中性子先行核の取り扱いが異なるため、結果が大きく異なっているが、EFP モデルと EFP モデルに整合を持つ 6 群モデルの結果にはほとんど差が見られなかった。これは、先行研究における一点炉での検証結果と同様の傾向を示している。
次に、EFP モデルを用いて計算した中性子束分布の結果を示す。Fig. 3.3.3、Fig. 3.3.4 はそ
れぞれ高速群中性子束及び熱群中性子束の初期分布を表している。これらは、固有値計算で𝜈を規格化し𝑘 = 1.0として計算した際の中性子束分布である。

38

Fig. 3.3.3 高速群中性子束の初期分布

Fig. 3.3.4 熱群中性子束の初期分布

39

Fig. 3.3.5(a)~(c)はそれぞれ、反応度印加から𝑡 = 0.01[s]、0.1[s]、1.0[s]経過後の高速群中性子束分布を表している。同様に、Fig. 3.3.6(a)~(c)に熱群中性子束の結果を示す。

(a)𝑡 = 0.01[s](b)𝑡 = 0.1[s](c)𝑡 = 1.0[s]Fig. 3.3.5 高速群中性子束分布の変化(a)𝑡 = 0.01[s](b)𝑡 = 0.1[s](c)𝑡 = 1.0[s]Fig. 3.3.6 熱群中性子束分布の変化
次に、EFP モデルによって計算した FP の炉内平均数密度変化及び空間分布の結果について記述する。核種として、短半減期核種である Co-78(𝑇1/2 = 4.84[ms])及び、ガス状 FP である I-139(𝑇1/2 = 2.28[s])の 2 核種についての結果を示す。Fig. 3.3.7 は Co-78 の炉内平均数密度変化である。

40

Fig. 3.3.7 反応度印加後の Co-78 炉内平均数密度変化
また、Co-78 の数密度分布変化をカラーマップ形式で Fig. 3.3.8(a)~(c)に示す。但し、(a)~(c)はそれぞれ反応度印加から𝑡 = 0.01[s]、0.1[s]、1.0[s]経過後の数密度分布を示している。

(a)𝑡 = 0.01[s](b)𝑡 = 0.1[s](c)𝑡 = 1.0[s]Fig. 3.3.8 Co-78 数密度分布の変化
Co-78 と同様の結果を I-139 についても示す。Fig. 3.3.9 は、反応度印加後の I-139 炉内平均数密度変化を示した図である。

41

Fig. 3.3.9 反応度印加後の I-139 炉内平均数密度変化
また、I-139 の数密度分布変化をカラーマップ形式で (a)~(c)に示す。但し、(a)~(c)は Co-78と同様にそれぞれ反応度印加から𝑡 = 0.01[s]、0.1[s]、1.0[s]経過後の数密度分布を示している。

(a)𝑡 = 0.01[s](b)𝑡 = 0.1[s](c)𝑡 = 1.0[s]Fig. 3.3.10 I-139 数密度分布の変化
Fig. 3.3.8 及び Fig. 3.3.10 から、Co-78 と I-139 の数密度は異なる傾向で変化していることが分かる。これは、各々の核種の生成過程に依存すると考えられる。Co-78 の炉内平均数密度
変化は Fig. 3.3.2 の炉心出力変化と似た傾向を示している。これは、Co-78 が主に核分裂反応によって生成されるため、核分裂反応量(炉心出力)の増加に伴い Co-78 も増加するためだ
42

と考えられる。一方、I-139 は他核種の崩壊による生成が支配的であるため、出力の上昇か
ら時間遅れで数密度が増加する。そのため、Co-78 とは異なる傾向で数密度が増加したと考えられる。また、空間分布についてどちらの核種も Material 間で非連続的に変化していることが分かる。一方、Fig. 3.3.5、Fig. 3.3.6 から分かるように、中性子束は Material の境界でも連続的に変化している。これは、中性子は空間メッシュ間で拡散するが、FP 核種は空間メッシュ間での拡散は無いためである。そのため、FP 核種は隣接する空間メッシュでの増減の影響を直接は受けないため、Material 境界部では非連続的な変化をする。また、これらの核種数密度分布は熱群中性子束分布と似た傾向を持つことが分かる。これは、FP が核分裂反応によって生成されることに起因する。Table 3.3.1 の核分裂断面積からも分かるように、
核分裂反応は主に熱群の中性子によって発生する。そのため、熱群中性子束分布と FP 核種
数密度分布が似た形となる。このように、空間依存動特性解析へ EFP モデルを適用させる
ことで、各々の核種を陽に取り扱った解析が実現可能となった。

3.4 ガス状 FP の漏洩事象(炉心全体)
1.3 節で述べたように、
EFP モデルでは事故時に発生した揮発性の FP 核種が炉心の全体、または一部領域から漏洩する事象を取り扱うことが出来る。本節及び次節では、TWIGL ベ
ンチマーク問題を対象として、この事象を考慮した場合の炉心相対出力変化、中性子束分布
変化、及び FP 核種数密度、分布の変化を比較する。なお、漏洩核種として Kr、Xe、I の全ての同位体を仮定した。また、漏洩定数も全ての核種で𝑙 = 2.0[/s](𝐿1/2 = 0.35[s])とした。
まず、炉心全体から漏洩が発生したと仮定した場合の結果を示す。は炉心相対出力の変化を示している。なお、Fig. 3.4.1 では、ⅰ漏洩を考慮しない場合、ⅱ反応度印加直後(𝑡 = 0.0[s])に漏洩が発生した場合、ⅲ反応度印加から𝑡 = 0.1[s]後に漏洩が発生した場合の 3 ケースの結果を示す。

43

Fig. 3.4.1 FP 漏洩事象での炉心相対出力変化
Fig. 3.4.1 より、即発跳躍部分では漏洩による出力の変化がほとんど発生しないことが分かる。これは、出力の上昇における即発中性子の割合が大きいため、遅発中性子発生量の減少が影響しないためと考えられる。一方、即発跳躍を超えた時間帯になると出力上昇が急激に緩やかになっていることが分かる。これは漏洩核種の中性子源としての寄与が大きいため、漏洩により中性子源が減少し中性子束が Reference と比べて低下したためである。
次に、Fig. 3.4.2(a)~(c)では反応度印加直後に漏洩が発生したケースにおける高速群中性子束分布の変化を示す。但し、Fig. 3.4.2(a)~(c)はそれぞれ、反応度印加から𝑡 = 0.01[s]、0.1[s]、1.0[s]経過後の中性子束分布を表している。

(a)𝑡 = 0.01[s](b)𝑡 = 0.1[s](c)𝑡 = 1.0[s]Fig. 3.4.2 高速群中性子束分布の変化(炉心全体漏洩、𝑡 = 0.0[s])44Fig. 3.3.5 の漏洩無の事象と比べ、特にピークの中性子束が低くなっていることが分かる。
𝑡 = 1.0[s]での中性子束分布を漏洩の有無で比較するため、カラーマップで表示したものをFig. 3.4.3 に示す。

Fig. 3.4.3 𝑡 = 1.0[s]での高速群中性子束分布(左：漏洩無、右：漏洩有)
また、Fig. 3.4.4(a)~(c)では、同様に熱群中性子束分布の変化を示す。

(a)𝑡 = 0.01[s](b)𝑡 = 0.1[s](c)𝑡 = 1.0[s]Fig. 3.4.4 熱群中性子束分布の変化(炉心全体漏洩、𝑡 = 0.0[s])熱群中性子束についても、高速群同様に Fig. 3.3.6 に比べて上昇が緩やかであることが分かる。熱群中性子束についても𝑡 = 1.0[s]での中性子束分布を Fig. 3.4.5 にカラーマップで示す。

45

Fig. 3.4.5 𝑡 = 1.0[s]での熱速群中性子束分布(左：漏洩無、右：漏洩有)熱群についても漏洩により全体的な中性子束が低下したことが分かる。また特に中性子束がピークとなる炉心中央部で顕著に見られた。
次に、FP の炉内平均数密度変化及び空間分布の結果を示す。また、3.3 節と同様に Co-78と I-139 の結果を示す。Fig. 3.4.6 は漏洩の有無による Co-78 数密度の変化である。

Fig. 3.4.6 FP 漏洩事象での Co-78 炉内平均数密度変化
Fig. 3.3.8 の漏洩無の事象と同様に、Co-78 の数密度は出力を追従する形で変化していることが分かる。そのため、即発跳躍領域では漏洩の影響をほとんど受けないことが分かる。また、46

Co-78 の数密度分布変化をカラーマップ形式で Fig. 3.4.7(a)~(c)に示す。

(a)𝑡 = 0.01[s](b)𝑡 = 0.1[s](c)𝑡 = 1.0[s]Fig. 3.4.7 Co-78 数密度分布の変化(炉心全体漏洩、𝑡 = 0.0[s])同様に、I-139 の炉内平均数密度変化を Fig. 3.4.8 に示す。

Fig. 3.4.8 FP 漏洩事象での I-139 炉内平均数密度変化
I-139 は Co-78 と異なり揮発性を持ち、本解析において漏洩を仮定している。そのため、漏洩事象発生直後から指数関数的に減少していることが分かる。数密度分布の変化についても同様に Fig. 3.4.9(a)~(c)にカラーマップ形式で示す。

47

(a)𝑡 = 0.01[s](b)𝑡 = 0.1[s](c)𝑡 = 1.0[s]Fig. 3.4.9 I-139 数密度分布の変化(炉心全体漏洩、𝑡 = 0.0[s])Fig. 3.4.9 より、全領域から I-139 が漏洩し減少していく様子が見て取れる。以上の解析より、漏洩事象での炉内平均数密度変化には、漏洩による直接的な減少と間接的な減少の二種類があることが分かった。直接的な減少は、漏洩により核種数密度が減少することで、漏洩を仮定した核種に限られる。一方、間接的な減少は漏洩により核分裂量(出力)が減少し、それにより核種生成量が減少する効果や、漏洩核種を親核種に持つ核種の生成量が減少する効
果が考えられる。次節では漏洩事象をより詳細に再現するため、部分的な漏洩を模擬した解析を行う。

3.5 ガス状 FP の漏洩事象(一部領域)
前節では炉心全体からガス状 FP が漏洩する事象を解析した。本節では漏洩領域を Fig.
3.3.1 の Material1 部分のみと仮定し、出力変化、中性子束分布、FP 数密度について結果を示す。は反応度印加後𝑡 = 0.1[s]後に漏洩定数𝑙 = 2.0で漏洩した事象での炉心相対出力変化を示したものである。

48

Fig. 3.5.1 FP 漏洩事象での炉心相対出力変化(一部領域)
同時刻から同じ漏洩定数で漏洩が開始したが、一部領域のみからの方が、出力が高くなるというのは自明である。次に、漏洩後の中性子束分布の比較を行う。Fig. 3.5.2 は、(a)漏洩しな
かった場合、 (b)𝑡 = 0.1[s]に炉心全体から漏洩した場合、(c)𝑡 = 0.1[s]に Material 1 から漏洩が発生した場合の 3 ケースについて、反応度印加から𝑡 = 1.0[s]後の高速群中性子束分布を比較したものである。

(a)漏洩無(b)全体漏洩(c)一部漏洩

Fig. 3.5.2 𝑡 = 1.0[s]での高速群中性子束分布(漏洩領域)

Fig. 3.5.2 で、漏洩無の場合と比較すると、全炉心漏洩事象の場合、炉心全体で中性子束分布が低下していることが分かる。一方、一部漏洩事象の場合、漏洩を仮定した部分とその周辺
49

のみ中性子束分布が低くなっていることが分かる。このように、漏洩を仮定した部分付近で中性子束分布の減少が顕著に見られた。
次に、3.4 節と同様に Co-78、I-139 の二核種について炉内平均数密度変化及び数密度分布を示す。Fig. 3.5.3 FP 漏洩事象での Co-78 炉内平均数密度変化(漏洩領域)は Co-78 の炉内平均数密度変化を示したものである。

Fig. 3.5.3 FP 漏洩事象での Co-78 炉内平均数密度変化(漏洩領域)
Co-78 については、炉心全体漏洩事象と同様に出力の変化を追従する様な挙動を示している
ことが分かる。では、核種数密度分布を比較するため、反応度印加から𝑡 = 1.0[s]後の核種数密度分布を Fig. 3.5.4(a)~(c)に示す。但し(a)~(c)はそれぞれ、(a)漏洩しなかった場合、 (b)𝑡 =0.1[s]に炉心全体から漏洩した場合、(c)𝑡 = 0.1[s]に Material 1 から漏洩が発生した場合の 3ケースである。

50

(a)漏洩無(b)全体漏洩(c)一部漏洩

Fig. 3.5.4 𝑡 = 1.0[s]での Co-78 数密度分布(漏洩領域)
Fig. 3.5.4 から、漏洩無の事象と比べて、一部漏洩事象では Material 1 の領域でわずかに数密度が減少していることが分かる。漏洩しない核種については、部分的な漏洩の場合、数密度分布に対して漏洩の影響は少ないことが分かる。次に、I-139 の結果を示す。Fig. 3.5.5 は I139 の炉内平均数密度変化を示したものである。

Fig. 3.5.5 FP 漏洩事象での I-139 炉内平均数密度変化(漏洩領域)
漏洩領域を限定した場合、自明ではあるが炉内の核種減少量は低下し、時間経過に伴い増加に変わった。次に、核種数密度分布の比較を行う。反応度印加から𝑡 = 1.0[s]後の核種数密度51分布を Fig. 3.5.6(a)~(c)に示す。(a)~(c)のケースは Fig. 3.5.4 と同様である。

(a)漏洩無(b)全体漏洩(c)一部漏洩

Fig. 3.5.6 𝑡 = 1.0[s]での I-139 数密度分布(漏洩領域)
Fig. 3.5.6 より、部分的な漏洩により、Material 1 の領域の核種数密度が著しく低下している

ことが分かる。一方、そのほかの領域では、現象は見られるものの大きな変化とはなっていない。これは、FP については空間メッシュ間での拡散を扱っていないため、Material 1 の領域での核種数密度減少が他の領域へは直接的には影響を与えないためである。

3.6 本章のまとめ
本章では、第 2 章で記述した計算手法を基に、EFP モデルによる空間依存動特性解析を行い、その結果について記述した。まず 3.1 節、3.2 節では一次元軽水反射体付きの体系に対して、EFP モデルによる動特性解析を行い、PME 法と𝜃法の性能を評価した。また、それを受け行列指数の近似精度に関する検討を行った。これらの検討により、超短半減期核種が出力変化へほとんど影響を与えないという知見を得た。また、3.3 節以降では体系を二次元に拡張し、TWIGL ベンチマーク問題を対象に解析を行った。これにより、従来の６群モデル
では再現が難しい解析を行うことが可能になった。また、ガス状 FP の漏洩事象が出力変化に与える影響についても、EFP モデルによって解析が可能となった。

52

第4章 結論
本研究では、原子炉動特性解析に着目し、従来の 6 群モデルに対してより物理現象を詳細に扱うことが出来る EFP モデルの開発を行った。この EFP モデルは先行研究において一点
炉動特性解析には適用が完了している。その結果では、6 群モデルとの間に重要な差異はほとんど見られなかった。そのため、本研究では空間依存動特性解析へ EFP モデルを適用させることで、6 群モデルでは再現することが困難な事象の解析を実現することを目的とした。
先行研究の一点炉解析では、行列形式で記述された微分方程式に対して、行列全体に行列指数法を適用することが可能であった。しかし、本研究で EFP モデルを適用させることを目的としていた空間依存動特性方程式については、係数行列全体の行列指数関数を精度良く計算することが困難であった。そのため、空間依存動特性方程式に行列指数法を適用させる手法として、部分行列指数法という新たな手法を開発した。空間依存動特性解析に適用可能なその他の計算手法としては、𝜃法などが考えられるが、本研究で行列指数法の適用を目指した。これは、他研究において、燃焼計算では超短半減期核種の計算精度が重要であるということが示唆されていたためである[4]。そのため、本研究でも低次の近似手法である𝜃法と比較する形で、超短半減期核種の解析を行った。その結果、低次の近似では超短半減期核種の近似精度が悪化するという結果を得ることが出来た。しかし、動特性解析の主たる結果である炉心出力には、近似の次数に対する差異は発生しなかった。これより、動特性解析において、超短半減期核種は炉心出力の結果に大きな影響を与えないという知見を得ることが出来た。
また、体系を二次元に拡張し、TWIGL ベンチマーク問題を対象に動特性解析を行った。
この解析によって、各々の核種に着目した解析やガス状 FP の漏洩事象など、従来の 6 群モ
デルでは再現が難しい事象について再現が可能となった。特に、ガス状 FP の漏洩事象では、漏洩により炉心出力に顕著な影響が出るということが分かった。
本研究により、空間依存動特性解析に対して EFP モデルを適用させることが可能となった。今後の取り組みとしては、遅発ガンマ線が減速材へ直接エネルギーを付与する効果を取り扱った動特性解析などが考えられる。また、PME 法では(2.6.5)式で表される線形近似を用いている部分が存在する。この近似精度を二次、三次と高めることで、粗いタイムメッシュ
における核種数密度算出精度の改善が見込めるなどの検討が考えられる。

53

参考文献
[1] 原子力安全委員会、「発電用軽水炉原子炉施設の安全評価に関する審査指針」、原子力安全委員会(2001)[2] G. R. Keepin, T. F. Wimett, and R. K. Zeigler, “Delayed Neutrons from Fissionable Isotopes of
Uranium, Plutonium, and Thorium”, PHYSICAL REVIEW JOURNALS ARCHIEVE, August 15,
1957[3] ジェームス J.ドゥデルシュタット、ルイス J.ハミルトン、「原子炉の理論と解析」、現代
工学社、(1980)[4] Yosuke Kawamoto, Go Chiba, Masashi Tsuji, Tadashi Narabayashi, “Numerical solution of matrix
exponential in burn-up equation using mini-max polynomial approximation”, Annals of Nuclear
Energy, Volume 80, June 2015, Pages 219-224[5] 辻田浩介、「Characteristics 法を用いた動特性計算に関する研究」、2012 年度名古屋大学
修士論文(2013)[6]Yuichiro Ban, Tomohiro Endo, Akio Yamamoto, Yoshihiro Yamane, “Explicit Time Integration
Scheme Using Krylov Subspace Method for Reactor Kinetics Equation”, Journal of NUCLEARSCIENCE and TECHNOLOGY, Vol. 48, No. 2, p.243-255(2011)[7] 片桐耕司、「EFP モデルによる原子炉動特性解析の高度化に関する研究」、2016 年度北海道大学卒業論文(2017)[8]Koji Katagiri, Go Chiba, “Nuclear Reactor One-Point Kinetic Calculations with Explicit Fission
Product Model”, Proc. RPHA2017, Chengdu, China, Aug. 24-25, 2017, (2017) (CD-ROM)54

謝辞
本研究に際し、千葉豪准教授、山本泰功助教授からは様々なご指導を頂きました。特に担当教官である千葉豪准教授には、学部時代を含め 3 年間非常に多くのことを教えて頂きました。ここに心から感謝申し上げます。また、炉工の同期や後輩、修了された先輩方を含め、研究生活を支えて頂いた全ての皆様にお礼を申し上げます。ありがとうございました。

55