炉物理プログラム演習：随伴方程式と摂動計算
千葉豪

中性子増倍体系に微小な変動（摂動）が加わった場合（例えば、体系の温度上昇や制御棒の挿入）を考える。摂動が加わる前の系の実効増倍率を k1 、摂動が加わった後の系の実効増倍率を k2 とすると、この摂動により体系に印加された
反応度 ρ は

11k2 − k1−=k1k2k1 k2

ρ=(1)と書くことが出来る。炉物理計算では、系に与えられた摂動によりどの程度の反応度が印加されたかを正確に評価することが重要である。


問題１：幅が 50 cm の一次元平板を考える。組成は一様で、エネルギー 1 群の定数の値は D=10 cm、Σa =1[/cm]、

νΣf =1[/cm] とする。また、境界条件はゼロ中性子束とする。この体系に対して、（左端から）24 cm から 26 cm の2 cm 幅の領域について Σa =1.1、1.01、1.001、1.0001 [/cm] と変化させたときに印加される反応度を k1 、k2 から計算せよ。

数値計算により中性子増倍率を求める場合、計算される中性子増倍率には有限の有効桁が存在する。従って、摂動が大きい場合には摂動前後の実効増倍率から式 (1) を用いて反応度を計算できるが、摂動が小さい場合には有効な桁で反応度を計算することが困難となる。例えば、摂動前後の系の実効増倍率として k1 = 1.00000、k2 = 1.00001 が得られたとするならば、この摂動による反応度は 0.00001 としか得られず、反応度は 0.000014 かもしれないし、0.000005 かもしれない。数値計算での反復解法において、実効増倍率に対する収束条件を厳しくすることで、得られる実効増倍率の有効桁数を大きくすれば、反応度の有効桁数を増やすことも可能であるが、その場合、収束に必要となる反復回数が莫大なものとなる可能性がある。そのような微小な反応度を計算する方法として摂動計算がある。
中性子拡散方程式を演算子を用いて以下のように記述する。

Aϕ =1Fϕk

一次元平板体系では、A、F はそれぞれ以下のように定義される。
()∑ddϕg(Aϕ)g = −Dg+ Σr,g ϕg −Σg′ →g ϕg′ ,dxdxg ′ ̸=g∑(F ϕ)g = χgνΣf,g′ ϕg′ .

(2)

(3)
(4)g′

式 (1) に対して、その随伴方程式を以下のように定義する。

A† ϕ† =1 † †F ϕk(5)ここで ϕ† は随伴中性子束である。A† 、F † は、それぞれ A、F に対する随伴演算子であり、以下の性質を有する。

⟨ †⟩⟨⟩ϕ , Aϕ = ϕ, A† ϕ† ,⟨ †⟩⟨⟩ϕ , F ϕ = ϕ, F † ϕ†

(6)(7)

ここで、記号 ⟨⟩ は全位相空間に対する積分（例えば拡散方程式の場合には、全空間、全エネルギーでの積分）を示す 1 。
⟨⟩
今回の場合は全空間、全エネルギー群に対する積分となり、例えば ϕ† , Aϕ のうち中性子吸収が関わる項については∑∑
ϕ†m,g Σa,m,g ϕm,g Vm(8)mgと記述される。ここで、ϕm,g はメッシュm における g 群の中性子束を、Vm はメッシュm の体積（一次元の場合は幅）を示す。
⟨⟩ ⟨⟩問題２： ϕ† , Aϕ 、 ϕ† , F ϕ の具体的な表式を記述せよ。
1 随伴方程式の物理的意味については、ホームページに掲載されている「原子炉特別実験：随伴中性子束と制御棒価値」が参考になるであろう。

1式 (2) で示される系に摂動が与えられ、演算子が A′ = A + ∆A、F ′ = F + ∆F となり、実効増倍率が k ′ になったとする。摂動前の随伴方程式と、摂動後の方程式（通常の方程式は随伴方程式に対して前進方程式と呼ばれる）は以下のように書ける。

A† ϕ†=A ′ ϕ′=1 † †F ϕ ,k1 ′ ′
F ϕk′(9)(10)

式 (9) の両辺に ϕ′ 、式 (10) の両辺に ϕ† を乗じ、それぞれの式の両辺について全位相空間で積分すると、以下の式を得る。

1 ⟨ ′ † †⟩ϕ ,F ϕ ,k1 ⟨ † ′ ′⟩ϕ ,F ϕk′⟨ ′ † †⟩ϕ ,A ϕ=⟨ † ′ ′⟩ϕ ,A ϕ=(11)(12)

式 (11) に対して、式 (6)(7) で示される随伴演算子の性質を用いると、以下のように変形できる。

⟩ 1⟨ †⟩⟨ †ϕ , Aϕ′ =ϕ , F ϕ′k
式 (12) から式 (13) を両辺について引くと、以下の式を得る。
()⟩⟨ †⟩
11 ⟨ † ′ ′⟩ 1 ⟨ †′ϕ , ∆Aϕ =−ϕ ,F ϕ +ϕ , ∆F ϕ′′kkk(13)(14)

これより、摂動による反応度が以下の積分から計算できることが分かる。

⟩ ⟨⟩1⟨ †ϕ , ∆F ϕ′ − ϕ† , ∆Aϕ′11
ρ= − ′ = kk k⟨ϕ† , F ′ ϕ′ ⟩(15)

この式に基づくと、反応度を単純な積分計算から得ることが出来るため、値が微小な反応度であっても精度良く計算することができる。
また、摂動後の中性子束を ϕ′ = ϕ + ∆ϕ と書くと、∆F ϕ′ 、∆Aϕ′ は以下のように書き直せる。

∆F ϕ′
∆Aϕ′=

∆F (ϕ + ∆ϕ) ,(16)=∆A (ϕ + ∆ϕ)(17)

これに対して、演算子の摂動量と中性子束の変動量の積 2（すなわち ∆F ∆ϕ、∆A∆ϕ）が無視できるものと仮定すると、以下の式を得る。

⟩ ⟨⟩1⟨ †ϕ , ∆F ϕ − ϕ† , ∆Aϕk
ρ=⟨ϕ† , F ′ ϕ⟩(18)

式 (15) に基づく摂動計算を厳密摂動と呼ぶのに対して、式 (18) に基づく計算を一次摂動と呼ぶ。厳密摂動では、摂動前の随伴中性子束 ϕ† と摂動後の中性子束 ϕ′ が必要となるため、摂動に応じて ϕ′ の計算が必要となる。一方、一次摂動では、必要となる中性子束 ϕ と随伴中性子束 ϕ† はいずれも摂動前のものであるため、予めこれらを計算しておけば、任意の摂動に対して反応度を計算することができる。この一次摂動の考え方は感度係数の計算で用いられる。
エネルギー 1 群の中性子拡散固有値方程式の随伴式は次のように書ける。
()d †1dD(x) ϕ (x) + Σa (x)ϕ† (x) = νΣf (x)ϕ† (x)−dxdxk(19)これから明らかなように 1 群の随伴拡散方程式は前進方程式と同一となる。すなわち、ϕ(x) = ϕ† (x) となる。このように、随伴方程式と前進方程式が一致する場合を自己随伴と呼ぶ。


問題３：問題１の摂動による反応度を、ϕ(x)、ϕ† (x)(= ϕ(x)) を用いて一次摂動計算により求めよ。また、問題１の結果と比較せよ。


問題１での摂動は Σa に対してのみ与えられるため、∆F = 0 である。
′†

問題４：問題１の摂動による反応度を、ϕ (x)、ϕ (x)(= ϕ(x)) を用いて厳密摂動計算により求めよ。
2 これを高次成分とも呼ぶ。

2一方、多群の随伴方程式は前進方程式と異なる。多群の随伴拡散方程式を以下に示す。
()
∑∑
∑dd †
1−Dg (x) ϕg (x) + Σa,g (x)ϕ†g (x) +νΣf,gχg′ (x)ϕ†g′ (x) +Σg→g′ ϕ†g′ (x)Σg→g′ ϕ†g (x) =dxdxk′′′gg(20)g

多群随伴方程式の解法は基本的に前進方程式と同一であるが、エネルギー群について解く順番が異なる。すなわち、2 群問題の場合は、核分裂中性子源を仮定したのち、2 群の随伴中性子束分布を計算し、その後、（2 群からの散乱源を加え∑て）1 群の計算を行うこととなる。なお、随伴式における g 群の核分裂中性子源は νΣf,g
χg′ ϕg′ で与えられる。そこ′g∑で、前進式と形式を合わせるため、擬似的なスカラー核分裂源として F =χg′ ϕg′ を考え、各群の核分裂中性子源をg′νΣf,g F として計算するようにすればよいであろう。n 回目の外部反復における k n は k n = F n /(F n−1 /k n−1 ) で求めればよいが、外部反復における初期核分裂源を F 0 = 1 となるように規格化し、k 0 = 1 としてやれば、F n−1 /k n−1 = 1 が保証されるため、k n = F n として計算できる。


問題５：厚さ 100 cm の一次元体系について、各定数が以下のように与えられた場合の随伴方程式を解き、k と ϕ†gを求めよ。外部境界条件はゼロ中性子束とする。また、計算した k が、同一の体系の前進方程式を解いて得た k と一致することを確かめよ。

Energy groupD[cm]Σa [/cm]

Σg,g+1 [/cm]νΣf [/cm]

χ

11.5

0.01

0.020.005

1.0

20.4

0.1-0.1410.0


問題６：問題５の組成をもつ厚さ 100 cm の一次元平板について、（左端から）48 cm から 52 cm までの領域の Σa,1を 5%増加させたときの反応度を摂動計算により求めよ。また、νΣf,2 、Σ1→2 に対しても同様の摂動を与えた場合の反応度も計算せよ。


問題７：問題６と同様の体系について、炉心全体の D1 、D2 を 1%増加させたときの反応度を摂動計算により求めよ。

拡散係数の摂動に対する反応度の計算は多少複雑となるので、以下に解説を加える。
拡散係数に摂動が与えられた場合の、摂動理論に基づく反応度計算式 (15) の分子 ρnume は以下のように与えられる。
()}∫ {dϕ′nume† dρ=ϕ∆Ddx(21)
dxdx[]{()()}∫dϕ′
dϕ†
dϕ′=ϕ† ∆D−
∆Ddx
(22)
dxdxdx式 (21) の積分は全空間に対するものであるため、式 (22) の右辺第一項は外部境界上での値となるが、この項はゼロとなる 3 。従って、最終的に ρnume は次のように書ける。
( ′ )}∫ {( † )dϕdϕ∆Ddx
ρnume = −dxdx(23)

拡散方程式では中性子流 J は以下の式で記述される。

J = −Ddϕdx

従って、式 (23) は中性子流を用いて以下のように書ける。
( ′ )}∫ {( † )JJnumeρ=−∆DdxDD′(24)(25)

一般的な拡散方程式の数値解法では、メッシュ表面の中性子流 J をメッシュ表面と中心の中性子束の値から定義するため、中性子流 J はメッシュの右半分、左半分で一定値をとると考えてよい。これらを JL 、JR のように記述したとき、
ρnume は以下のように記述される。
)∆D ∆xi ( † ′† ′JJ+J
ρnume = −J(26)LRLRDD′ 2
3 ゼロ中性子束境界条件が課された場合には境界上で ϕ† = 0 となるためこの項はゼロとなり、反射境界条件が課された場合も境界上で dϕ′ /dx = 0となるため同様となる。

3次に、核分裂反応が無い、固定源方程式の場合について考えよう。
このような場合の解くべき中性子拡散方程式は以下となる。

Aϕ = S(27)

ここで、中性子束分布によって決まる応答パラメータ R として以下を考える。

R = ⟨Σϕ⟩(28)

例えば、位置 r⃗′ での反応率が応答パラメータである場合は、Σ は以下のように定義できるであろう。
()
Σ = δ ⃗r − r⃗′ Σ(⃗r)(29)

また、以下のように Σ を定義すると、R は位置 r⃗′ の g ′ 群の中性子束そのものとなる。
()
Σ = δ ⃗r − r⃗′ δg,g′(30)ここで式 (27) に対して、以下の随伴方程式を定義する。

A† ϕ† = Σ(31)

この式の両辺には ϕ を乗じ、全位相空間で積分すると以下の式を得る。

⟨⟩ϕ, A† ϕ† = ⟨Σϕ⟩ = R(32)⟨⟩ ⟨⟩ ⟨⟩
この式の左辺は ϕ, A† ϕ† = ϕ† , Aϕ = ϕ† S のように変形できるので、以下の式が得られる。
⟨⟩
R = ϕ† S(33)

これより、応答パラメータ R を計算する方法として、式 (27) を解いて中性子束分布を求め、それを用いるというものに加えて、式 (31) を解いて随伴中性子束分布を求め、それと外部源分布 S を用いるというものがあることが分かる。
さて、ある基準状態に対する随伴方程式が式 (31) のように与えられる一方、それに摂動が加えられた系（A′ = A + ∆A）での拡散方程式が以下のように与えられるものとする。

A′ ϕ′ = S(34)

式 (31) の両辺に ϕ′ 、式 (34) の両辺に ϕ† を乗じ、全位相空間について積分すると以下を得る。

⟨ ′ † †⟩ϕ ,A ϕ=⟨ † ′ ′⟩ϕ ,A ϕ=式 (35) は以下のように変形できる。

⟨Σϕ′ ⟩ = ⟨Σ (ϕ + ∆ϕ)⟩ = R + ∆R,⟨ †⟩Sϕ = R,

⟨ ′ † †⟩ ⟨ †⟩ϕ , A ϕ = ϕ , Aϕ′ = R + ∆R(35)(36)(37)

従って、式 (37) から式 (36) を辺々引くことにより、以下が得られる。

⟨⟩⟨⟩∆R = − ϕ† , ∆Aϕ′ ≈ − ϕ† , ∆Aϕ(38)すなわち、ϕ† を一度計算しておけば、任意の摂動 ∆A に対する応答パラメータの変動 ∆R を容易に計算出来ることが分かる。

4