平成 28 年度

卒業論文

EFP モデルによる原子炉動特性解析モデルの高度化に関する研究

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

北海道大学

工学部

機械知能工学科

目次
第1章

序論

1.1

緒論 ..................................................................................................................... 1

1.2

即発中性子と遅発中性子 ................................................................................. 1

1.3

6 群モデルと EFP モデル ................................................................................. 2

1.4

研究目的 ............................................................................................................. 4

第2章

解析手法

2.1

解析仕様 ............................................................................................................. 52.2

一点炉動特性方程式 ......................................................................................... 5

2.3

Pade 近似法 ........................................................................................................ 6

2.4

Marquardt 法による 6 群定数のフィッティング ........................................... 72.5

解析体系、解析条件 ......................................................................................... 9

第3章

解析結果

3.1

Marquardt 法によるフィッティング ............................................................. 10

3.2

フィッティングによる差異軽減 ................................................................... 11

3.3

正の反応度印加 ............................................................................................... 133.4

負の反応度印加 ............................................................................................... 143.5

重み付けによる再現性の向上 ....................................................................... 163.5

負の反応度における再現性向上に対する考察 ........................................... 22

3.6

未臨界パルス照射体系 ................................................................................... 23

第4章

結論 ............................................................................................. 25

第 1 章 序論
1.1

緒論

原子力発電所では、その中心部である炉心において中性子による核分裂連鎖反応が起こり、膨大なエネルギーが発生している。中性子の挙動を、より正確に計算機上で再現することで、信頼性の高い炉心設計解析・安全解析を行うことが出来る。
原子炉施設の安全設計の妥当性を確認するためには、異常状態すなわち「運転時の異常な過渡変化」及び「事故」について解析と評価を行う必要がある。このうち「運転時の異常な過渡変化」とは、発生頻度は高い（原子炉の寿命期間中に 1 回以上発生する可能性のあるもの）が、燃料及び冷却材圧力バウンダリの破損は起こらず、放射線影響の観点からは無視できる事象である。一方、「事故」とは「運転時の異常な過渡変化」を超える異常な状態であって、発生する頻度は低いものの、発生した場合は原子炉施設及び周辺に重大な影響を及ぼす恐れのある想定事象のことである。運転時の異常な過渡変化では「原子炉起動時における制御棒の異常な引き抜き」や「出力運転中の制御棒の異常な引き抜き」、
事故においては、「制御棒飛び出し」（PWR）、「制御棒落下」（BWR）といった事象が想定される。これらの事象により原子炉に反応度が投入されると、炉心での中性子の挙動が変化し、中性子数が増加する。中性子数つまり出力が急激に増加すると、燃料破損、更にはメルトダウン等のシビアアクシデントへと発展してしまう。原子炉の安全性を高め、より正確な設計解析を行うためには、このような反応度印加時の中性子変化を追跡する正確なモデルが必要不可欠である。

1.2

即発中性子と遅発中性子

核分裂反応が発生すると、反応と同時に励起状態の核分裂片（Fission fragments）から中性子が発生する。この中性子を「即発中性子」（Prompt neutron）と言う。即発中性子、即発ガンマ線を放出し基底状態になったものが核分裂生成物(Fission Products、FP)であり、この中には崩壊により中性子を放出する核種が存在する。この核種のことを「遅発中性子先行核」(Delayed neutron precursor)と言い、遅発中性子先行核の崩壊により、核分裂反応から時間遅れで発生する中性子のことを「遅発中性子」(Delayed neutron)と言う。
235U 熱中性子核分裂において、即発中性子に対する遅発中性子の割合は約 0.65%程度であるが、この遅発中性子によって原子炉の制御が可能となっている。
例として始めに、遅発中性子が存在しない体系を考える。この体系では、時刻tにおける中性子数n(t)は以下で示される単一の指数関数で記述される。
−1( ) = (0) exp(1.2.1)= (0) exp
:即発中性子平均寿命[s]:中性子増倍係数= ( − 1)⁄ ：反応度
= ⁄ :中性子生成時間[s]= 0.000015[s]とし、反応度 = 0.0001を印加すると、1 秒後の中性子数は約 785 倍にまで増加する。
次に、1 種類の遅発中性子先行核が存在する体系を考える。この体系で時刻tにおける中性子数n(t)及び遅発中性子先行核濃度C(t)が従う微分方程式は以下のように記述される。

1

( )

=

"( )

=

−

( ) + !"( )( ) − !"( )(1.2.2)(1.2.3)β:遅発中性子割合
λ:遅発中性子先行核の崩壊定数[% &' ]二元連立微分方程式(1.2.2)式、(1.2.3)式はラプラス変換によって解析的に解くことが出来、以下の(1.2.4)式を得る。
( ) = ( exp() ) + * exp(+ )(1.2.4)
ここで、(、*、)、+は 、 、 、λから定まる定数である。
この場合、同条件( = 0.0001、 = 0.000015、! = 0.078)での 1 秒後の中性子数増加は1.69%にとどまる。
(1.2.1)式及び(1.2.4)式で表される中性子数変化を、Fig.１-1 に示す。
このように遅発中性子の存在により、人間が原子炉を制御することが可能になっている
ことが分かる。即発中性子に比べ、僅かしか存在しない遅発中性子であるが、原子炉の制
御・動特性解析を行う上で極めて重要な役割を持つ。

Fig.1-1: 中性子数の時間変化

1.3 6 群モデルと EFP モデル
遅発中性子先行核としては、87Br、88Br、137I、138I、139I などの核種が存在し、その数はおよそ 300 核種である。またその各々が異なる半減期で崩壊し遅発中性子を放出する。
そのため計算コストなどの面から、原子炉動特性解析では遅発中性子先行核をその崩壊定数に従って 6 つにグルーピングし、6 種の仮想的な先行核として取り扱う 6 群モデルが一般的に用いられている。Table 1-1 に 6 群モデルのグループ分けの例を示す。

2

Table 1-1: 235U 熱中性子核分裂による遅発中性子先行核の 6 群データ(JENDL-4.0 より)
グループ

半減期T'⁄/ [s]崩壊定数λ0 [s &' ]先行核生成割合(1

1

55.72

0.01240.033072

22.72

0.03050.21901

3

6.22

0.111

0.19594

4

2.30

0.301

0.39496

5

0.610

1.14

0.11504

6

0.230

3.01

0.04199

現在用いられている 6 群モデルは、実際の原子炉における挙動を再現するように作られており、原子力発電所の設計・解析に広く用いられている。しかし、この 6 群モデルは多くの物理現象を簡略化しモデル化している。そのため 6 群モデルによる計算は、計算時間
等の点で極めて実用的である一方、直接考慮することが出来ないいくつかの物理現象が存在する。また、6 群モデルは、従来の原子炉設計・安全解析への適用を目的として開発されたため、想定内の対象・物理現象・物理パラメータを精度良く評価することは可能であるが、想定を超える対象や物理現象、パラメータについて対応する事は出来ない。
そのため、本研究室では先行研究において、全ての核分裂生成物(FP)を陽に取り扱う原子炉解析モデル EFP(Explicit representation of Fission Products)モデルの開発を行った[1]。EFP モデルでは、FP の崩壊による他の FP の生成及び中性子放出を詳細に解析し追跡している。これにより EFP モデルでは FP の崩壊チェーンを厳密に再現することが可能になり、より実現象に近い体系を再現することが出来る。これを陽に扱うことを、これ以降チェーン効果(Chain-effect)と呼称する。また、EFP モデルでは遅発ガンマ線も陽に計算することが出来る。6 群モデルと EFP モデルの概略図を Fig.1-2 に示す。

Fig. 1-2: 6 群モデルと EFP モデル概略図

3

Table1-1 で示したような、一般的に用いられている 6 群モデルは実験データを基に作成されている。これは具体的には、瞬時核分裂後の時間依存の中性子放出数に関する実験を再現するように 6 群定数が決められている。一方、EFP モデルで用いる核データは理論計算に依るところが大きい。そのため、この 6 群モデルを EFP モデルとの比較に用いると、実験データと EFP モデルでの基になる核データの違いに起因する差異が発生してしまう。
この差異が先行研究において顕著に現れてしまい本来の EFP モデルと 6 群モデル間の差を比較するのが困難になってしまっていた。本研究はこのデータに起因する差異を取り除くことを目的としコードの開発及び検討を行った。
EFP モデルではチェーン効果を無視して計算を行うことも出来る。チェーン効果を無視
した場合、中性子を放出しない核種は無視することが出来るため、扱く核種数を減らすことが出来る。本研究で用いた JENDL/FP のデータの場合、扱う核種数はチェーン効果を考慮した場合には 1284 核種、無視した場合には 315 核種となる。

1.4 研究目的
1.3 節で述べたように EFP モデルでは先行研究において、基となる核データの違いによって、解析結果に差異が生じてしまう問題が発生していた。そこで本研究では、最小二乗法の一種である Marquardt 法を用いた非線形フィッティングにより、基となる核データの整合を取ることで、それに起因する差異を排除することを目指した。またそれによってモデルの違いにのみ起因する差異のみを抽出し、6 群モデルによる EFP モデルの再現精度について確かめ、検討を行う。
また、チェーン効果を考慮するモデル及び考慮しないモデルそれぞれで、同条件の現象を模擬した際、それぞれの EFP モデルと 6 群モデル間における差異を比較することによって、チェーン効果を考慮することによる再現性の違いを検討する。
これらの比較について、更に 6 群モデルによる EFP モデルの再現性を向上させるため、
Marquardt 法における重み付けを導入する。この重み付けとは、基となるデータに重みを付けることにより、特定のデータにおけるフィッティング誤差を低減することが出来る。
それにより更なる再現性の向上が見込まれる。また重み付けの条件を変更させることにより、重み付け条件に伴う再現性の変化を確認する。

4

第 2 章 解析手法
2.1 解析仕様
本研究では 235U の熱中性子核分裂時の空間依存性を無視した、一点炉における過渡挙動
を EFP モデルと 6 群モデルで計算し比較した。そのために一点炉動特性方程式を導出し、それを解くために C++によるプログラムの開発を行った。データは JENDL /FPD-2011 及び FPY-2011 を用いた。

2.2 一点炉動特性方程式
まず 6 群モデルの方程式について記述する。空間依存性を無視した体系において、時刻tにおける中性子数 ( )及び、i番目の家系に属する遅発中性子先行核密度"1 ( )は以下の(2.2.1)式,(2.2.2)式で表される連立微分方程式に従う[2]。
"1 ( )

( )

=

=
1

Λ

−

4

( ) + 3 !1 "1 ( ),15'( ) − !1 "1 ( ), (i = 1,2, … ,6)(2.2.1)(2.2.2)= − 1⁄ :反応度
:遅発中性子割合:中性子生成時間[s]次に EFP モデルの場合、動特性方程式は以下の(2.2.3)式、(2.2.4)式で表すことが出来る。
( )"1 ( )

=

=1
>

−
Λ

<

( ) + 3 !1 "1 ( ); 1 ,(2.2.3)( ) − !1 "1 ( ) + 3 !? "? ( );?1 , (@ = 1,2, … , A)

(2.2.4)

=

<

?5'

15'

; 1 :核種@の崩壊の際の中性子放出確率= :核分裂収率>:核分裂あたりの平均放出中性子数;?1 :核種Bの崩壊により核種@が生成する確率
A:考慮する核種数(2.2.4)式は、チェーン効果を考慮した場合の式である。(2.2.4)式の右辺第三項は、FP の崩壊による娘核の生成を表現している。これらを数値的に解くことにより、中性子数の時間変化を解析することが出来る。
(2.2.4)式で、核分裂収率を用いているが、核分裂収率には「独立核分裂収率」と「累積核分裂収率」の 2 種類が存在する。前者は核分裂反応によりその核種が直接発生する量を示すのに対して、後者は核分裂反応で発生した様々な各種の崩壊の結果、その核種が発生する量を示している。そのため、チェーン効果を無視した場合は累積核分裂収率を、チェーン効果を考慮した場合では独立核分裂収率を用いる。また、チェーン効果を無視した場
合、(2.2.4)式の右辺第三項は無視される。
次に、遅発中性子割合 は以下の式で表すことが出来る。
>C :核分裂あたりの平均放出遅発中性子数

=

5

>C
>

(2.2.5)

>は核分裂あたりの平均放出即発中性子数>D と核分裂あたりの平均放出遅発中性子数>C の和で表すことが出来る。また、>C は以下のように表すことが出来る。
>C = 3 =E@ ; @(2.2.6)

1=E:累積核分裂収率F
EFP モデルで解析を行う際、中性子増倍係数として、即発中性子増倍係数 D を用いるのが便利である。即発中性子増倍係数とは、即発中性子のみによる中性子増倍係数であり、以下の(2.2.7)式で表される。
(2.2.7)
D = (1 − )これらより、(2.2.3)式は以下のように書き換えることが出来、(2.2.8)式を得る。
( )

=
=
=
=

<

−

−1

( ) + 3 !1 "1 ( ); 1
15'−

(1 − ) − 1
D−1<( ) + 3 !1 "1 ( ); 1
15'<( ) + 3 !1 "1 ( ); 1
15'<( ) + 3 !1 "1 ( ); 1 ,(2.2.8)15'

また、(2.2.4)式も(2.2.9)式のように書き換えが出来る。
"1 ( )

=1
=
>
=
=
=

1

<

( ) − !1 "1 ( ) + 3 !? "? ( );?1

∙
D

∙

D

1−
∙

H1 −

D =1

>D

?5'

<

=1
( ) − !1 "1 ( ) + 3 !? "? ( );?1
>
?5'

=1

>C
I J>C + >D K
>C + >D<<( ) − !1 "1 ( ) + 3 !? "? ( );?1

( ) − !1 "1 ( ) + 3 !? "? ( );?1 ,
?5'

?5'

(@ = 1,2, … , A)

(2.2.9)

次の節ではこれらの式の数値的解法について述べる。

2.3 Pade 近似法(2.2.8)式、(2.2.9)式を行列表記に直すとL( )= ML( )

(2.3.1)と書ける。なお、L( ) = ( ( ), "' ( ), "/ ( ), … , "< ( )) である。ここで() は転置行列を表している。
(2.3.1)式の解は形式的に以下のように書くことが出来る。
L( ) = exp(M ) L(0),(2.3.2))は行列指数関数(Matrix
ここでexp(Mexponential)である。行列指数関数を数値的に計算する方法はいくつあるが、本研究では Pade 近似法を用いた。スカラー変数Oに対する指数N

6N

exp (O)の[p/q]次の Pade 近似は次の(2.3.3)式の様に書くことが出来る。
ADS (O)
exp(O) ≈,TDS (O)
ここで
D

ADS (O) = 3
S

V5W

TDS (O) = 3

V5W() + + − )! )!
OV,() + +)! ! () − )!

() + + − )! +!
(−O)V ,() + +)! ! (+ − )!

(2.3.3)(2.3.4)(2.3.5)である。
同様に、行列指数関数に対しても以下の Pade 近似を適用することが出来る。
ADS (M ∆ )exp(M∆ ) ≈(2.3.6)
TDS (M ∆ )一般的には p=q とした近似が用いられている。本研究では p=q=6 の 6 次近似を用いた。

2.4 Marquardt 法による 6 群定数のフィッティング
1 章で述べたように、本研究では 6 群モデルと EFP モデル間で基になる核データの整合を取ることにより、基となる核データの違いに起因する誤差の排除を目指した。この節ではフィッティング法について記述する。
まず 235U が一個、中性子によって核分裂（瞬時核分裂）した後の遅発中性子放出率
C ( )に着目する。初めに先行核種数密度の初期値に核分裂収率を与える。このときの収率には、チェーン効果を考慮する場合は独立核分裂収率を、無視する場合は累積核分裂収率を用いた。
"1 (0) = =1
(2.4.1)
即発中性子は考えないため、k Z = 0とおく。すると、C0 (t)が従う微分方程式は(2.4.2)の様に記述される。
"1 ( )<= −!1 "1 ( ) + 3 !? "? ( );?1 ,時刻 の遅発中性子放出率は以下の様になる。
<?5'

C ( ) = 3 ; 1 !1 "1 ( )15'(2.4.2)(2.4.3)

これを Pade 近似により数値的に解く。出力は(2.4.3)式で求めた C ( )に時間と 100 を乗じて、 C ( ) × × 100とした。これはx軸をログスケールとしてグラフに描画する場合、積分量をリニアスケールの場合と等しくするためである。 C ( )をプロットした際、ある区間の積分量はその区間で発生した遅発中性子量を示す。そのためログスケールでもこれを保存することで発生している遅発中性子量を分かりやすく判断することが出来る。
この遅発中性子放出率時間変化を 6 群モデルで表すと以下の(2.4.4)式となる。
4>C 3 ]? !? expJ−!? K,
C ( ) = \\\]? :遅発中性子先行核生成割合!? :遅発中性子先行核崩壊定数[1/s]?

但し、J∑? ]? = 1K

7

(2.4.4)

この、]? 及び!? を変数としてフィッティングを行う。
フィッティングには最小二乗法の一種である Marquardt 法を用いた[4][5]。まず、フィッティングに用いるデータを抽出した。このデータは EFP モデルによって計算した 235U
瞬時核分裂後の遅発中性子放出率のx軸を対数的に等間隔に区切り抽出したものである。次に変数に初期値を与えた。Marquardt 法では、初期値の値によって収束性が大きく変化し、適切な初期値を与えなければ良いフィッティング結果を得ることが出来ない。
この初期値を_W = ((W' , (W / , … , (W < )N とおく。ここで、Aはフィッティングパラメータの数を示す。また、(2.4.4)式を_W と の関数とし、`(_W , 1 )とする。フィッティング誤差の評価として、相対二乗残差a W を以下の様に与える。
e/c1 − `(_W , 1 )W
a = 3bd(2.4.5)`(_W , 1 )15':フィッティングデータ数
c1 :基となるデータの値(EFP モデルによる C ( )の計算値)
次に、誤差行列Mを以下の様に定義する。
ef1? = J1 + ! g1? K 3 hV
V5'i`(_W , V ) i`(_W , V )
∙
,
i(W 1
i(W?

(2.4.6)（但し、g1? :クロネッカーのデルタ）

hV :重み係数! :対角係数
対角係数!の初期値には 1 を与えた。
更に、残差に偏微分を乗じた N 行 1 列の行列を作成する。
ei`(_p , 1 )3 h1 nc1 − `(_W , 1 )oi(W' tm
15'
l esi`(_p , 1 )slW
3 h1 nc1 − `(_ , 1 )oj= li(W / s
l 15'
s
⋮
l
s
l e
p
i`(_ , 1 )s
3 h1 nc1 − `(_W , 1 )oi(W < rk15'

(2.4.7)

この M と j を用いて、新しいパラメータ_' = (('' , (' / , … , (' < )N を以下の様に求める。
_' = _W + M&' j(2.4.8)
次に新しい_' を用いて、更新二乗残差a' を計算する。もし、a' < aであれば、_' を新し
い初期値_W とし、またa' を新しいaとして、更に! を 0.1 倍する。また、もしa' > aであれば、_' 及びa' を破棄し! を 10 倍する。この作業を行った後、再び(2.4.6)式の計算まで戻w
り、それ以降の計算を繰り返す。適当な収束条件( ≈ 1)で計算を終了し、_W を取り出すこwxとでフィッティングが完了する。

8

2.5 解析体系、解析条件
様々な条件による検証を行うため、まず初めに、EFP モデルで 235U 瞬時核分裂後の遅発中性子放出率を計算し、そのデータを基に Marquardt 法によるフィッティングを行う。
これによって重み付けの条件を変えたいくつかの 6 群定数を作成する。この際、チェーン
効果を無視した体系で用いている 6 群定数は、チェーン効果を無視した EFP モデルで計算した 235U 瞬時核分裂後の遅発中性子放出率を基に作成しているため、チェーン効果を考慮した体系の比較で用いている 6 群定数とは異なる。
作成した 6 群定数を用いて、以下の 2 つの検証を行う。
ⅰ:臨界定常状態から正/負の反応度を印加した際の中性子数の変化
ⅱ:未臨界体系に中性子を打ち込んだ際の中性子数の変化ⅰで計算を行う際の、臨界定常状態時の遅発中性子先行核密度（初期値）は、EFP モデ
ルでは k=1 で計算を行い、十分時間が経過した後の組成を用いた。6 群モデルでは"( )⁄ = 0より、以下の(2.5.1)式で算出した。
"1 (0)=

1(0) − !1 "1 (0) = 0,

"1 (0) =(0),(2.5.1)!1
また、ⅱの検証においては遅発中性子先行核密度の初期値は全て 0 とし、印加反応度は−25￠とした。
以降の計算では、EFP モデルとそれに整合を持つ 6 群モデルで様々な体系を模擬し、EFP モデル、6 群モデル間の差異を算出することにより、6 群モデルによる EFP モデルの再現性を検証する。また、 >D = 2.42048、>C = 0.0169281、 = 0.000015、= 0.00694512を用いた。
19

第 3 章 解析結果
3.1 Marquardt 法によるフィッティング
まず 2.4 節で述べた Marquardt 法によるフィッティングの結果を示す。(2.4.5)式中の_Wに初期値として Keepin の 235U 熱中性子核分裂時の 6 群定数を代入する。この定数をTable 3-1 に示す。
Table 3-1:Keepin の 6 群定数
α{
!?

1

2

3

4

5

6

3.31e-2

2.19e-1

1.96e-1

3.95e-1

1.15e-1

4.20e-2

1.24e-2

3.05e-2

1.11e-1

3.01e-1

1.14

3.01

また(2.4.5)式中のc1 (EFP モデルによる計算値)として、Fig.3-1 の点で示されたデータを用いた。このデータ数は 57 個となっている。
Keepin の 6 群定数を用いて計算した 235U 瞬時核分裂後の遅発中性子放出率と、EFP モ
デルで計算した値の比較を Fig.3-1 に示す。これより、フィッティングを行う前では大きな乖離が生じている事が分かる。

Fig. 3-1: フィッティング前の遅発中性子放出率
次に Fig. 3-2 にフィッティング後の 6 群定数を用いて計算した 235U 瞬時核分裂後の遅発
中性子放出率とフィッティングに用いた EFP モデルでの計算値を示す。また Table 3-2 にはフィッティングにより作成した新たな 6 群定数の値を示す。

10

Fig. 3--2: フィッティング後の遅発中性子放出率

Table 3-2:Marquardt
3 2:Marquardt 法によるフィッティング後の 6 群定数
α{
!?

1

2

3

4

5

6

5.35e-2

1.91e--1

2.1
2.15e-1

3.63e-1

1.19e-1

5.87ee-2

1.45e-2

3.51e--2

1.49e
1.49e-1

3.63e-1

1.39

4.09

Fig. 3-2 より、フィッティングにより大きな乖離が解消していることが分かる。次節からはこのこの整合を取った
整合を取った 6 群モデルを用いて反応度を印加する計算を行う。
群モデルを用いて反応度を印加する計算を行う。

3.2 フィッティングによる差異
フィッティングによる 異軽減
この節では反応度印加時の中性子数変化を、
この節では反応度印加時の中性子数変化を、Keepin の 6 群定数を用いて解析した
反応度印加時の中性子数変化を、Keepin
群定数を用いて解析した場合
場合と、EFP
EFP モデルと整合を取った 6 群定数を用いた場合で比較する。
群定数を用いた場合で比較する。Fig. 3
3-3 に正の反応度(ρ = 0.001, +14.4
4￠)を印加した際の変化を、を印加した際の変化を、を印加した際の変化を、Fig.
3--4 に負の反応度(ρに負の反応度 = −00.001、 − 14.4￠)14を印加した際の変化を示す。
を印加した際の変化を示す。Fig.3-3
及び Fig.3-4 内で用いられている EFP
EFP-consistentconsistent とはフィッティングにより EFP モデルと整合をとった 6 群モデルを用いて計算を行ったも
のである。

11

Fig. 3-3: 印加反応度| = p. pp}での中性子数時間変化

Fig. 3-4: 印加反応度 = −0.001での中性子数時間変化
Fig.3-3 及び Fig.3-4 より、基となる核データの整合を取った事により、EFP モデルと 6群モデル間で、核データの違いに起因する大きな乖離が解消した。これにより、EFP モデルと 6 群モデルとの比較がモデルに起因する差異のみを見ることが出来るようになった。
次節以降では主に、EFP モデルと 6 群モデル間の相対誤差を比較する。

12

3.3 正の反応度印加
Fig.3-5 に、臨界定常状態に正の反応度(+25￠、 + 50￠、 + 60￠、 + 70￠)を印加した
際の中性子数を示す。なお、Fig.3-5 における中性子数の計算値は EFP モデルを用いて計算したものである。また縦軸は対数軸としている。

Fig. 3-5:反応度印加時の中性子数変化(正の反応度印加時)
Fig. 3-6(a)～(d)に、臨界定常状態から正の反応度(+25￠、 + 50￠、 + 60￠、 + 70￠)を印加した際の、EFP モデルと 6 群モデル間における中性子数の差異を示す。中性子数の相対差は以下の(3.3.1)式で定義される。

a =e~•€•‚ƒ & e„…†
e„…†

× 100

(3.3.1)(3.3.1)式において ‡ˆ‰ 及び、 4Š‹CŒ• はそれぞれ中性子数の EFP モデルによる計算値と 6群モデルによる計算値を表している。

13

(a)+25￠投入(b) +50￠投入(c) +60￠投入(d) +70￠投入

Fig. 3-6: EFP モデルと 6 群モデル間相対誤差(正の反応度印加時)
Fig. 3-6 より、25￠の反応度印加では差異が 2%程度と小さい。また、反応度が大きくなるにつれて差異が拡大するが、Fig.3-5 から分かるように中性子数の絶対値が極端に大きいため 30%程度の差異が問題とならないと考えられる。

3.4 負の反応度印加
Fig.3-7 に、臨界定常状態から負の反応度(−25￠、 − 50￠、 − 60￠、 − 70￠)を印加した際の中性子数を示す。なお、Fig.3-7 においても中性子数の計算値は EFP モデルを用いて計算したものである。また縦軸は同様に対数軸としている。

14

Fig. 3-7 :反応度印加時の中性子数変化(負の反応度印加時)
Fig. 3-8(a)～(d)に、臨界定常状態から負の反応度(−25￠, −50￠, −60￠, −70￠)を印加した際の、EFP モデルと 6 群モデル間中性子数誤差時間変化を示す。

(b) −50￠投入(a)−25￠投入(c) −60￠投入(d) −70￠投入
Fig.3-8:チェーン効果を考慮した体系での EFP モデルと 6 群モデル間相対誤差15Fig. 3-8 より、どのケースにおいても差異は数%程度と小さいことが分かった。チェー
ン効果を考慮した体系では、およそ 150[s]を境に 6 群モデルが中性子数を過大評価する傾向へと変化している。時間経過とともに相対誤差が正の方向へ拡大している傾向は、正の反応度印加の場合と逆の結果となった。また、相対誤差の絶対値は負の反応度印加の場合の方が小さくなっている。
チェーン効果を無視した体系では、中性子数の時間変化は複数の指数関数の和のみで表すことが出来る。そのため、6 群モデルによる EFP モデルの再現性は、チェーン効果を考慮した体系より小さくなるのではないかと予測していた。しかし、3.3 節及び 3.4 節の結果から分かるように、検証した全ての反応度印加条件において、チェーン効果を無視した体系の再現性が悪化するという傾向が見られた。

3.5 重み付けによる再現性の向上
2 章で述べたように、Marquardt 法ではデータに重み付けをすることで特定の区間でのフィッティング精度を向上させることが出来る。そこで本節では特定の区間への重み付けによる 6 群モデルでの EFP モデルの再現性の変化を検証する。
まず重み付け区間による再現性の変化を調べるため、フィッティングに用いている 235U
瞬時核分裂後の遅発中性子放出率グラフの x 軸を Fig.3-9 に示すように、6 つの区間に分割する。なお以降の検証において EFP モデルでの計算はチェーン効果を考慮した体系で行っている。

Fig. 3-9: 重み付け区間の分割

16

それぞれの区間に 10 倍の重みをつけて作成した 6 群モデルと EFP モデルの中性子数相
対誤差を以下に示す。なお Fig.3-10 では正の反応度(+60￠)印加を、Fig.3-11 には負の反応度(−60￠)印加を示す。

Fig. 3-10: 重み付け区間による再現性の変化(+60￠投入)

Fig. 3-11: 重み付け区間による再現性の変化(-60￠投入)

17

Fig.3-10 より正の反応度の場合、1 番の区間への重み付けにより最も再現性が向上する
事が分かった。一方、Fig.3-11 より負の反応度印加の場合は 1 番に加えて 6 番への重み付
けでも再現性が向上する事が分かった。
次に重さによる再現性の変化を検証する。先ほどの検証で再現性が向上する重み付け区間が分かったので、その区間へ 2 倍から 20 倍の重み付けをして 6 群定数を作成する。
Fig.3-12 に正の反応度印加の場合を、Fig.3-13,Fig.3-14 に負の反応度の場合の 1 番への重み付け及び 6 番区間への重み付けをそれぞれ示す。

Fig. 3-12: 1 番に重み付けした時の重さによる再現性の変化(+60￠投入)

18

Fig. 3-13: 1 番に重み付けした時の重さによる再現性の変化(-60￠投入)

Fig. 3-14: 6 番に重み付けした時の重さによる再現性の変化(-60￠投入)

19

Fig. 3-12 より 10 倍の重みを付けた際に最も再現性が良くなる事が分かった。また Fig.
3-13 及び Fig. 3-14 より負の反応度印加の際は 1 番、6 番のどちらの区間でも 20 倍の重み付けにより時間経過時の再現性が最も向上した。
以上の検証より、正の反応度の場合は 1 番の区間に 10 倍の重み付けを、また負の反応
度の場合には 1 番及び 6 番の区間に 20 倍の重みを付けることで再現性が向上する事が分かった。最後にこれらの重み付けの印加反応度による普遍性を検証する。
Fig. 3-15、Fig. 3-16、Fig. 3-17 には、それぞれ(正の反応度、1 番、10 倍)、(負の反応度、1 番、20 倍)、(負の反応度、6 番、20 倍)の条件で印加反応度を変化させた際の再現性変化を示す。なおそれぞれの条件で、(a)±25￠、(b)±50￠、(c)±60￠、(d)±70￠の反応度
投入時の結果である。

(a)+25￠投入(b) +50￠投入(c) +60￠投入(d) +70￠投入

Fig. 3-15: 1 番区間,10 倍の条件で反応度による再現性の変化

20

(b) −50￠投入(a)−25￠投入(c) −60￠投入(d) −70￠投入
Fig. 3-16: 1 番区間、20 倍の条件で反応度による再現性の変化(b) −50￠投入(a)−25￠投入(c) −60￠投入(d) −70￠投入

Fig. 3-17: 6 番区間、20 倍の条件で反応度による再現性の変化

21

Fig. 3-15、Fig. 3-16、Fig. 3-17 より、全ての検証で全ての反応度において再現性の向上を確認する事が出来た。

3.5 負の反応度における再現性向上に対する考察
3.4 節より正の反応度では短時間区間に、負の反応度では加えて長時間区間への重み付
けも再現性向上に効果があると分かった。本節ではなぜ負の反応度の場合のみ長時間区間
への重み付けに効果があるのかについて、遅発中性子放出率に着目して考察を行う。
Fig. 3-18 には 235U 瞬時核分裂後の遅発中性子放出率を点線で示している。また実線では反応度印加後の遅発中性子先行核密度を基にした遅発中性子放出率を示している。これについて説明を述べる。まず 3.3 節及び 3.4 節の検証と同様に臨界定常状態から反応度を印加し、中性子数 ( )及び@番目の家系に属する遅発中性子先行核密度"1 ( )を計算する。十分な時間(600[s])が経過した後の遅発中性子先行核密度"1 ( = 600)を初期値として(2.4.3)式に代入し、 C ( )を計算する。これをピークの高さが等しくなるよう規格化したものが実線で示したものである。

Fig. 3-18: 反応度印加後の遅発中性子放出率
この図より、正の反応度を印加した場合は瞬時核分裂後の遅発中性子放出率に近い形となっているが、負の反応度印加後の遅発中性子放出率は大きく長時間区間へとピークがずれている事が分かる。そのため、負の反応度印加の場合は図の10/ < < 10• 区間への重み付けにより再現性が向上したと考えられる。

22

3.6 未臨界パルス照射体系
未臨界体系にt = 0[s]で中性子を打ち込んだ際の中性子数変化を EFP モデルと 6 群モデルで比較する。Fig.3-19 に、未臨界体系(−25￠)に中性子を打ち込んだ際の中性子数変化
を EFP モデルで計算したものと、6 群モデルで計算したものを示す。この図では両者の中性子数の相対差異も点線で示す。また拡大図(0 ≤ t ≤ 0.1)を Fig.3-20 に示す。

Fig. 3-19: 未臨界パルス照射体系における中性子数変化

Fig. 3-20: 未臨界パルス照射体系における中性子数変化(拡大図)23

Fig.3-19 より、t < 0.02[%]部分の即発中性子による減衰時に1.5[%]程度の差異が生じているが、中性子数が遅発中性子によるもののみになると差異は 0.5[%]以下となっている事が分かる。これより、未臨界パルス照射体系では EFP モデルと 6 群モデルの間で、パルス打ち込みの瞬間に 1[%]程度の差異が生じるが、時間経過とともに差異は減少していくということが分かった。

24

第 4 章 結論
本研究では、先行研究で開発された EFP モデルを引き継ぎ、6 群モデルとの比較精度向上を目指し検証を行った。先行研究で比較に用いられていた 6 群定数は実験データから作成されたものであった。そのため、実験データと EFP モデルの計算結果の差異が反応度を印加した際の EFP モデルと 6 群モデル間の結果にも現れてしまい、モデルの違いによる差異のみの比較が出来ていなかった。そこで Marquardt 法を用いることにより、基となる核データ起因の大きな乖離を解消することが出来た。これにより、本研究の目的である EFP モデルと 6 群モデルの、モデルに起因する差異を検証する事が可能になった。検証結果より、EFP モデルと 6 群モデル間には中性子数評価に以下の傾向があるという知見を得た。
(ⅰ)正の反応度を印加した場合、6 群モデルは EFP モデルと比べて中性子数を過小評価する傾向がある。この傾向は印加反応度が大きいほどまた時間が長くなるほど顕著となる。
(ⅱ)負の反応度を印加した場合、およそ 150 秒を境に、6 群モデルによる中性子数の評価が過小評価から過大評価へと変化する。それ以降は、誤差は増大する傾向がありこれは印
加反応度の絶対値が大きいほど顕著となる。
(ⅲ)チェーン効果を無視した体系ではチェーン効果を考慮した体系に比べて、検証した全ての反応度において 6 群モデルによる EFP モデルの再現性が悪化した。
(ⅳ)未臨界体系にパルス状に中性子を打ち込んだ場合、打ち込んだ直後(t < 0.02[%])に 1%程度の差異が生じるが、その後は差異が減少し大きな差は見られない。
また Marquardt 法でのフィッティング時における重み付けによって、6 群モデルによる
EFP モデルの再現性向上を図った。検証を行った結果、再現性が向上する重み付けの条件
は正負の反応度で異なっていた。しかし検証した全ての反応度において符号が一致している反応度では同じ条件で再現性が向上した。この正負の反応度によって再現性が向上する
重み付け区間が異なる現象を遅発中性子放出率に着目して説明を行った。
今後の展望としては、扱う物理現象の幅を広げ、より実現象に近いモデルを計算機上で再現することがあげられる。FP 崩壊に伴うガンマ線の放出をモデル化し、集合体体系において解析を行うことを検討中である。最終的には全炉心体系への拡張を行い、従来の 6群モデルでは対応できなかった事象、物理現象の数値シミュレーションを目指す。

25

参考文献
[1]松浦健太, 原子炉反応度事故解析の高度化に関する研究,北海道大学 平成 27 年度修士論文[2] Go Chiba, Masashi Tsuji, Tadashi Narabayashi , Sensitivity and uncertainty analysis for reactor
stable period induced by positive reactivity using one-point adjoint kinetics equation, Journal of
Nuclear Science and Technology, 50[12], p.1150-1160(2013)[3]瀬戸雄介, Seto’s Home Page マルカール(Marquardt)法 ,http://pmsl.planet.sci.kobe-u.ac.jp/~seto/?page_id=61(2016/12/2 アクセス)[4]若狭智嗣 ,九州大学理学研究院実験核物理研究室ホームページ ,教育・講義ノート,
数値計算法第 14 章 非線形関数の最小 2 乗法
http://ne.phys.kyusyu-u.ac.jp/~wakasa/cp/cp14.pdf(2016/12/2 アクセス)[5]三沢毅, 宇根崎博信, 卞哲浩,(2010) 原子炉物理実験, 京都大学学術出版会

謝辞
本研究を進めるにあたり、指導教官である千葉豪准教授、奈良林直教授に深く感謝を申し上げます。特に千葉豪准教授には、数多くのご指導ならびにご指摘を頂き、心より感謝いたします。また、たくさんのアドバイスを下さった先輩方、二平氏をはじめ同期の皆様にも厚く御礼申し上げます。

26

付録
EFP モデル（チェーン効果有）ソースコード(使用言語：C++)#include<iostream>#include<fstream>#include<vector>#include<sstream>#include<cmath>using namespace std;//階乗関数
double make_factorial(int fact){
double fa=1.;
for(int i=1;i<=fact;i++){fa*=i;}
return fa;}//行列乗算関数
void multi(vector< vector<double> >&A1,vector< vector<double> >&A2,vector<
vector<double> >&A3,int N){for(int i=0;i<N;i++){for(int j=0;j<N;j++){A3[i][j]=0;}}
for(int i=0;i<N;i++){for(int j=0;j<N;j++){if(A1[i][j]!=0){
double temp=A1[i][j];for(int k=0;k<N;k++){A3[i][k]+=temp*A2[j][k];
}
}
}
}
}
//逆行列関数void inverse(vector< vector<double> >&B1,vector< vector<double> >&B2,int N){for(int i=0;i<N;i++){for(int j=0;j<N;j++){if(i==j){B2[i][j]=1.;}}}
double w;//記憶用変数
for(int i=0;i<N;i++){w=1/B1[i][i];for(int j=0;j<N;j++){27

B1[i][j]*=w;
B2[i][j]*=w;}
for(int j=0;j<N;j++){if(i!=j){w=B1[j][i];for(int k=0;k<N;k++){B1[j][k]-=B1[i][k]*w;
B2[j][k]-=B2[i][k]*w;
}
}
}
}
}
//単位行列関数
void unit(vector< vector<double> >&I,int N){for(int i=0;i<N;i++){for(int j=0;j<N;j++){if(i==j){I[i][j]=1.;}else{
I[i][j]=0.;}}}}//Pade 近似
void pade(int N,vector< vector<double> >&A,vector<double>&ne,int p,int q,double
dt,double time,double rho){//累乗計算用行列 B,C (B=A^k)vector< vector<double> >B(N);for(int i=0;i<N;i++){B[i].resize(N);}vector< vector<double> >C(N);
for(int i=0;i<N;i++){C[i].resize(N);}unit(B,N);//fa1=p!,fa2=q!,fa3=(p+q)!
double fa1=make_factorial(p);double fa2=make_factorial(q);double fa3=make_factorial(p+q);//Npq,Dpq ベクトルを定義
vector< vector<double> >Npq(N);for(int i=0;i<N;i++){
Npq[i].resize(N);}vector< vector<double> >Dpq(N);for(int i=0;i<N;i++){28

Dpq[i].resize(N+1);}//Npq,Dpq を 0 に更新
for(int i=0;i<N;i++){for(int j=0;j<N;j++){Npq[i][j]=0.;
Dpq[i][j]=0.;}}//以下、Npq,Dpq の計算
for(int k=0;k<=p;k++){//B=A^kif(k==0){
unit(C,N);}else{
multi(B,A,C,N);}
B=C;//fa4=k!,fa5=(p+q-k)!,fa6=(p-k)! or (q-k)!
double fa4=make_factorial(k);double fa5=make_factorial(p+q-k);
double fa6=make_factorial(p-k);//Npq の計算
double temp=fa5*fa1/(fa3*fa4*fa6);
for(int i=0;i<N;i++){for(int j=0;j<N;j++){Npq[i][j]+=temp*pow(dt,k)*B[i][j];Dpq[i][j]+=temp*pow((-1.*dt),k)*B[i][j];}}}//Dpq の逆行列 x を求める
vector< vector<double> >x(N);for(int i=0;i<N;i++){x[i].resize(N);}inverse(Dpq,x,N);//行列 E=exp(AΔt)を求める
vector< vector<double> >E(N);for(int i=0;i<N;i++){E[i].resize(N);}multi(Npq,x,E,N);
vector<double> z(N);//記憶用配列
for(int k=0;k<=(time/dt);k++){if(k==0){cout<<k<<" "<<ne[0]<<"\n";}else{
for(int i=0;i<N;i++){29z[i]=0.;
for(int j=0;j<N;j++){z[i]+=E[i][j]*ne[j];}}for(int i=0;i<N;i++){ne[i]=z[i];}cout<<k*dt<<" "<<ne[0]<<"\n";}if(k==(time/dt)){
cout<<-1<<"\n";cout<<"##r= "<<rho<<"\n";}/*//臨界定常状態における初期組成出力用
if(k==(time/dt)){for(int i=0;i<N;i++){cout<<ne[i]<<"\n";
}
}
*/
}
}

int main(){//FissionYield ファイルの読み込み
ifstream fy("./NuclearData/FissionYieldU235_2016IndependentText.dat");//エネルギー点数
int le1;
fy>>le1;//データ格納マトリックス
vector<double> e(le1);//入射エネルギー
vector<int> nfp(le1);//核分裂生成物核種の数
vector< vector<int> > iz(le1);//原子番号vector< vector<int> > ia(le1);//質量数vector< vector<int> > ist(le1);//励起準位
vector< vector<double> > yc(le1);//核分裂収率
vector< vector<double> > dyc(le1);//絶対標準偏差//データの収納for(int i=0;i<le1;i++){fy>>e[i];fy>>nfp[i];iz[i].resize(nfp[i]);ia[i].resize(nfp[i]);ist[i].resize(nfp[i]);yc[i].resize(nfp[i]);dyc[i].resize(nfp[i]);for(int j=0;j<nfp[i];j++){30

fy>>iz[i][j];fy>>ia[i][j];fy>>ist[i][j];fy>>yc[i][j];fy>>dyc[i][j];}}//DecayData ファイルの読み込み
ifstream dd("./NuclearData/DecayDataText.dat");
int max_size=1285;
int channel_max=3;//データ格納マトリックス
vector<int>iz_d(max_size);//核種原子番号vector<int>ia_d(max_size);//核種質量数vector<int>liso(max_size);//励起準位
vector<int>lis(max_size);//vector<int>ndk(max_size);//崩壊チャンネル数
vector<double>awr(max_size);//vector<double>t(max_size);//半減期[s]vector<double>dt(max_size);//絶対標準偏差
vector<double>elp(max_size);//崩壊エネルギー
vector<double>delp(max_size);//絶対標準偏差
vector<double>eem(max_size);//vector<double>deem(max_size);//vector<double>ehp(max_size);//vector<double>dehp(max_size);//vector<vector <double> >br(channel_max);//各チャンネルへの分岐比
vector<vector <double> >dbr(channel_max);//絶対標準偏差
vector<vector <int> >iz2_d(channel_max);//崩壊後の原子番号
vector<vector <int> >ia2_d(channel_max);//崩壊後の質量数vector<vector <int> >irfs(channel_max);//崩壊後の励起準位
vector<vector <int> >idn(channel_max);//発生する中性子数
for(int i=0;i<channel_max;i++){br[i].resize(max_size);
dbr[i].resize(max_size);
iz2_d[i].resize(max_size);
ia2_d[i].resize(max_size);
irfs[i].resize(max_size);
idn[i].resize(max_size);}//データの収納
int loop=0;
while(1){if(loop>=max_size){
cout<<"error1!!\n";return 0;}dd>>iz_d[loop];if(iz_d[loop]==-1){
break;}
dd>>ia_d[loop];

31

dd>>liso[loop];dd>>awr[loop];dd>>t[loop];dd>>dt[loop];dd>>elp[loop];dd>>delp[loop];dd>>eem[loop];dd>>deem[loop];dd>>ehp[loop];dd>>dehp[loop];dd>>lis[loop];dd>>ndk[loop];if(ndk[loop]>channel_max){cout<<"error2!!\n";return 0;}
for(int i=0;i<ndk[loop];i++){dd>>br[i][loop];dd>>dbr[i][loop];dd>>iz2_d[i][loop];dd>>ia2_d[i][loop];dd>>irfs[i][loop];dd>>idn[i][loop];}loop+=1;}//loop=(核種数)vector< vector<int> > iz_fy(le1);
vector< vector<int> > ia_fy(le1);
vector< vector<int> > ist_fy(le1);
vector< vector<double> > yc_fy(le1);
vector< vector<double> > dyc_fy(le1);//DecayData と FissionYieldData の番号整合用配列
for(int i=0;i<le1;i++){iz_fy[i].resize(loop);
ia_fy[i].resize(loop);
ist_fy[i].resize(loop);
yc_fy[i].resize(loop);
dyc_fy[i].resize(loop);}//Decay データ数のカウント
int num_d=0;//DecayData と FissionYieldData の番号整合
for(int i=0;i<loop;i++){for(int j=0;j<nfp[0];j++){if((iz_d[i]==iz[0][j])&&(ia_d[i]==ia[0][j])&&(liso[i]==ist[0][j])){for(int k=0;k<le1;k++){iz_fy[k][i]=iz[k][j];ia_fy[k][i]=ia[k][j];ist_fy[k][i]=ist[k][j];yc_fy[k][i]=yc[k][j];32dyc_fy[k][i]=dyc[k][j];}num_d+=1;
break;}}}/*//番号の整合確認用
for(int i=0;i<loop;i++){if((iz_fy[0][i]!=iz_d[i])||(ia_fy[0][i]!=ia_d[i])||(ist_fy[0][i]!=liso[i])){cout<<i<<" failure\n";return 0;}
cout<<iz_fy[0][i]<<"\n";}*///emission Probability per decay
vector<double>Pn(num_d);for(int i=0;i<num_d;i++){for(int j=0;j<channel_max;j++){Pn[i]+=br[j][i]*idn[j][i];}}//Probability i→j(i 行核種が j 列核種になる確率 matrix)
vector< vector<double> >Pij(num_d);for(int i=0;i<num_d;i++){Pij[i].resize(num_d);}for(int i=0;i<num_d;i++){for(int j=0;j<channel_max;j++){for(int k=0;k<num_d;k++){if((iz_d[i]==iz2_d[j][k])&&(ia_d[i]==ia2_d[j][k])&&(liso[i]==irfs[j][k])){Pij[k][i]+=br[j][k];}}}}//即発中性子項を num_d に追加
num_d+=1;//A matrix を定義 (num_d)行(num_d)列 matrix
vector<vector <double> >A(num_d);for(int i=0;i<num_d;i++){A[i].resize(num_d,0);}//遅発中性子平均放出数 nu_d
ifstream data("./Data/EFP_nu_d.dat");double nu_d;33

data>>nu_d;//即発中性子寿命 l,遅発中性子割合 beta,即発中性子増倍率 kp,即発中性子平均放出数 np,
反応度 rhodouble l=0.000015;
double nu_p=2.42048;
double rho=0.002;
double beta=nu_d/(nu_d+nu_p);double kp=(1.-beta)/(1.-rho);double ln2=log(2);//中性子 matrix nevector<double>ne(num_d);ne[0]=1.;//A matrix 作成
A[0][0]=(kp-1)/l;for(int i=1;i<num_d;i++){A[i][0]=yc_fy[0][i-1]*kp/(l*nu_p);if(t[i-1]!=0){
A[0][i]=Pn[i-1]*ln2/t[i-1];}
for(int j=1;j<num_d;j++){if((i==j)&&(t[i-1]!=0)){
A[i][j]=-ln2/t[i-1];}
if((i!=j)&&(t[j-1]!=0)){
A[i][j]=Pij[j-1][i-1]*ln2/t[j-1];}}
if(t[i-1]!=0){
ne[i]=yc_fy[0][i-1]*kp*ne[0]/(l*nu_p*ln2/t[i-1]);}}//中性子の初期値を代入(定常状態）
ifstream data_neutron("./Data/num_neutron_EFP_steady.dat");for(int i=0;i<num_d;i++){
data_neutron>>ne[i];}double temp=ne[0];for(int i=0;i<num_d;i++){ne[i]=ne[i]/temp;}
int p=6;
int q=p;//Pade 近似次数
double delta=0.001;//タイムステップ
double TrialTime=600.;//試行時間
pade(num_d,A,ne,p,q,delta,TrialTime,rho);//Pade 近似及び出力
return 0;}34