6月1日課題

原子炉物理学特論 (課題８：最終回)
反復計算への加速法の導入
千葉豪
反復計算により収束解を得る際、計算対象が大規模になるほど長い計算時間を要する。計
算時間を短縮して効率よく収束解を得るために、反復計算に対して「加速」を行うことは極
めて効果的であり、これまでにさまざまな反復計算のための加速法が開発されている。今回
は、そのなかでも比較的簡単なものを扱い、加速法の導入効果を実感することを目的とする。
本題に入る前に、反復計算における収束（計算打ち切り）条件について簡単に触れる。反
復回数のみで反復計算を制御すると、計算値が十分に収束しているのに無駄な計算を行って
いたり、計算値が十分に収束していないのに計算を打ち切ってしまい想定した結果が得られ
ていなかったりする可能性がある。そのような観点から、収束条件の導入は反復計算におい
ては必須と言える。
中性子拡散方程式をべき乗法で解く場合の収束条件について考えよう。内部反復について
は、中性子束のある反復での計算値 φn と前回の反復での計算値 φn−1 の相対差の絶対値を全
ての空間メッシュで計算し、その最大値がある値以下となった場合に反復計算を打ち切る、
という操作を行うのが一般的である。これを式で書くと、g 群の中性子束を計算していると
するならば、
|φni,g − φn−1
i,g |
max
< inner
(1)
n−1
i
φi,g
となるだろう。一方、外側反復については、中性子源と実効増倍率について同様の条件を適
用する。すなわち、収束条件は式で
|k n − k n−1 |
k n−1
n
|S − S n−1 |
max i n−1i
i
Si
< outer,k ,
(2)
< outer,f issionsource
(3)
と書ける1 。
ここで一点、内部反復の収束条件は外側反復のそれよりも厳しくしないといけないことに
注意してもらいたい。そうでないと、いくら反復を行っても外部反復の収束条件が満足され
ない、ということがある。
1
一般的には、外部反復においても中性子束の収束を監視する。
1
問題１：厚さ 60[cm] の燃料（Fuel）の両端に厚さ 30[cm] の軽水反射体 (Water reflector）
が配置されている一次元平板原子炉を考え、境界条件は左右（すなわち軽水反射体の外
側において、）ともに中性子束をゼロとする。媒質の各定数が以下のように与えられた
場合の keff を求めよ（定数は課題７の問題２のものと同様である）。なお、空間メッシュ
の幅は１メッシュあたり 0.25[cm] とし、keff については小数点以下 5 桁までの収束を保
証するものとする。また、核分裂源分布の初期値として、燃料領域の片側 30[cm] の領
域のみに値を与え、それ以外はゼロとするものとする。
Medium
Fuel
Water
reflector
Energy
group
1
2
1
2
D
[cm]
1.58
0.271
1.41
0.117
Σa
[/cm]
0.0032
0.0930
0.0
0.0191
Σs,g→g+1
[/cm]
0.0178
0.0476
-
νΣf
[/cm]
0.0
0.168
0.0
0.0
χ
[-]
1.0
0.0
1.0
0.0
なお、あまりに計算時間がかかる、もしくはかからないようであれば、空間メッシュの幅を
調整し、この問題の計算時間が 2 分程度になるようにしてほしい。
それでは加速法として、SOR 法（Successive over-relaxation method、逐次過緩和法）を
導入しよう。
内部反復では、以下の式で中性子束の計算を行っている。
(k)
(k)
(Cn + Cn−1 + Σa,n ∆xn ) φ(k+1)
= Cn φn+1 + Cn−1 φn−1 + Sn ∆xn
n
(k)
(k)
(4)
(k+1)
ガウスザイデル法を用いた場合は、右辺の φn−1 （もしくは φn+1 ）が、φn−1 （もしくは
(k+1)
φn+1 ）となるであろう。
SOR 法では、(k + 1) 回目の反復における中性子束の推定値を以下の式で求める。
(k)
(k)
(5)
φ(k)
n )
(6)
(Cn + Cn−1 + Σa,n ∆xn ) φ̃n = Cn φn+1 + Cn−1 φn−1 + Sn ∆xn ,
φ(k+1)
n
φ(k)
n
=
+ ω(φ̃n −
ω は加速因子とよばれ、ω を 1.0 より大きくすると収束解を「先取り」しようとし、逆に 1.0
より小さくすると、「慎重に」反復を行うことになる。問題によっては ω を 1.0 未満にする
ことで収束が保証されるという場合があるが、今回は ω を 1.0 より大きくして、内部反復の
収束解をより高速に得ることを目指す。問題２：問題１について、ω として計算時間の観点から最適な値を試行計算により求
めよ。
次は、SOR 法を外側反復にも適用してみよう。計算式は以下の通りである。
X
S̃i =
νΣf,g φ(n)
g ,
(7)
g
(n+1)
Si
(n)
= Si
2
(n)
+ ωo (S̃i − Si )
(8)
ここで、外側反復のための加速因子を内側反復でのものと区別するため ωo と記述している。
問題３：問題２で作ったプログラムについて外側反復にも SOR 法を適用し、ωo として
計算時間の観点から最適な値を試行計算により求めよ。
ωo の値によっては外側反復で収束に至る回数を低減することができる一方、内側反復の収束
に要する時間が増加し、全体として計算時間が増加する、という場合もありえるであろう。
問題４：
「原子炉物理学特論」は来年度も同様の形式で実施する予定です。この講義を
改善するためのアイデア、意見を述べて下さい（忌憚のないものほど参考になります）。
なお、この回答については成績評価には反映しません。
次は、源外挿法と呼ばれる加速法を導入する。その前に、背景となる理論について簡単に
述べる（一部、課題７と重複する）。
中性子拡散方程式を演算子で記述すると以下のように書ける。
Aφ =
1
Fφ
k
(9)
べき乗法では反復 (n + 1) 回目の中性子束 φ(n+1) を、φ(n) と k (n) を用いて以下のように計
算する。
1
Aφ(n+1) = (n) F φ(n)
(10)
k
これに変形を施すと、
1
φ(n+1) = (n) A−1 F φ(n)
(11)
k
となる。ここで、議論を簡単にするため、k が既知であるとすると、φ(n+1) は以下のように
書ける。
2
n
n
1 −1
1 −1 (n)
1 −1
1
(n+1)
(n−1)
A F φ
A F
φ
= A Fφ =
= ··· =
φ(1) = n A−1 F φ(1)
k
k
k
k
(12)
さて、固有値方程式 (9) に対しては、複数の固有値 λi とそれに対応する固有ベクトル fi
が解として存在する。
1
Afi = F fi
(13)
λi
ここで、λ1 を固有値のなかで最も大きな値をとる第一固有値とする。従って、λ1 = k であ
り f1 = φ である。式 (13) を変形すると、以下の式を得る。
λi fi = A−1 F fi
(14)
さて、固有ベクトル群 fi はこの系に対して完備である（すなわち、この系のいかなるベク
トルも固有ベクトルの重ね合わせで表現できる）ため、反復の初期値 φ(1) は以下のように
固有ベクトルで展開することができる。
X
φ(1) =
aj fj
(15)
j=1
3
これを式 (12) に代入すると、以下の式を得る。
(n+1)
φ
X 1
X
n X
n
1
= n A−1 F
aj fj =
aj n A−1 F fj =
aj
k
k
j=1
j=1
j=1
λj
k
n
fj
反復がある程度進んだ場合（n がある程度の大きさになった場合）、j > 3 について (λ2 /k)n (λj /k)n となるので、以下の式が得られる。
n−1
λ2
(n)
f2 ,
(16)
φ
= a1 f1 + a2
k
n
λ2
(n+1)
φ
= a1 f1 + a2
f2
(17)
k
これを辺々引くと、以下の式が得られる。
(n+1)
φ
これを
−φ
(n)
= a2
n−1
λ2
λ2
−1
f2
k
k
n−1
λ2
φ(n+1) − φ(n) = a2
f2
λ2
k
−1
k
1
(18)
(19)
のように変形すると、式 (16) の右辺と式 (19) の右辺が同一となることから、求める中性子
束に対応する a1 f1 は
1
φ(n+1) − φ(n)
a1 f1 = φ(n) − (20)
λ2
−1
k
で推定することができる。
ただし、そのためには λ2 を知っている必要がある。そこで、式 (18) と同様の以下の式を
考える。
n−2
λ2
λ2
(n)
(n−1)
φ −φ
= a2
−1
f2
(21)
k
k
式 (18) と式 (21) より
φ(n+1) − φ(n)
λ2
= (n)
k
φ − φ(n−1)
(22)
が得られ、この式を用いて λ2 を見積もることができる。
以上では中性子束に対する加速として説明したが、実際には核分裂源に対して操作を行う
ことになる。つまり、λ2 の推定については、
S (n+1) − S (n)
λ2
= (n)
k
S − S (n−1)
を用い、核分裂源の加速については
S (n+1) = S (n) − 1
λ2
−1
k
4
S̃ − S (n)
(23)
(24)
で行うことになる2 。ここで、k としては反復中の推定値を用いることになる。なお、λ2 の
推定は各空間メッシュの中性子源を用いて行うため、推定値はメッシュ毎に得られることに
なるだろう。従って、それらから何らかの方法で体系全体としての推定値に焼き直す必要が
ある。また、λ2 がある程度精度良く求まっていなければならないので、例えばある（外側）
反復時点における λ2 の推定値が、前回の反復での推定値とそれほど違う場合にのみこの方
法を適用する、等の工夫が必要となる。
参考問題１：問題２で作ったプログラムについて源外挿法を適用し、計算時間の短縮を
試み、その結果について述べよ。
参考問題２：問題１では初期中性子源を空間的に「偏った」ものとして与えたが、空間
的に偏りのないものとした場合、収束に至る反復回数はどうなるか調べ、反復回数が変
化した理由について、固有ベクトルの観点から考察せよ。
2
見方を変えると、源外挿法は SOR 法のための最適加速因子を求める方法、と言えるかもしれない。
5

Download Report