Reactor physics programming exercise:
Adjoint equation and perturbation calculations
Go CHIBA

Let us consider that a small perturbation, such as temperature increase and control rod insertion, is given to a neutron
multiplication system. When neutron multiplication factors before and after the perturbation are denoted to as k1 andk2 , a reactivity inserted by this perturbation ρ can be presented as1k2 − k1
1−=.
k1
k2
k1 k2

ρ=(1)

In reactor physics calculations, it is quite important to accurately evaluate how much reactivity is inserted by perturbation.
✓✏

Problem 1: Let us consider a one-dimensional slab reactor with width of 50[cm]. This system is homogeneous, andone-group constants are given as D = 10[cm], Σa = 1[/cm] and νΣf = 1[/cm]. External boundary conditions are zeroneutron flux. Calculate reactivity from k1 and k2 when Σa in a region from 24[cm] to 26[cm] is changed to 1.1, 1.01,1.001 and 1.0001.
✒✑
When neutron multiplication factor is numerically calculated, effective digit of calculated k should be finite. When aperturbation is large, inserted reactivity can be easily calculated from two ks before and after the perturbation as shownin Eq. (1), but when a perturbation is small, it is difficult to obtain reactivity with high accuracy. For example, let
us consider that we obtain k1 = 1.00000 and k2 = 1.00001 in a perturbation. In this case, we can calculate reactivity
only as 0.00001, and this can be 0.000014 and 0.000005. By adopting strict (or rigorous) convergence criteria to iterative
calculations, we can obtain numerical results with higher accuracy and then reactivity is numerically calculated with highaccuracy, but it requires long computation time in the iterative calculations. To tackle this problem, the perturbation
calculations are generally used.
Neutron diffusion equation can be represented by operators as 1
Aφ =

1F φ.
k(2)

Adjoint equation to this equation is defined as1 † †F φ ,kwhere φ† is adjoint neutron flux. A† and F † are adjoint operators to A and F , and the following relations hold:A† φ† =φ† , Aφ†

φ ,Fφ(3)=φ, A† φ† ,(4)=† †(5)

φ, F φ,where hi is an integration over whole phase spaces. In the present case, this is an integration over whole space and wholeenergy group, and neutron absorption component in φ† , Aφ is explicitly represented asXX
φ†m,g Σa,m,g φm,g Vm ,(6)mgwhere φm,g is the gth group neutron flux in spatial mesh m and Vm is volume (or width in one-dimension) of mesh m.
✓✏

Problem 2: Derive the explicit expressions of φ† , Aφ and φ† , F φ .
✒✑Let us consider that a perturbation is given to a system presented as Eq. (2), operators are changed as A′ = A + ∆Aand F ′ = F + ∆F , and neutron multiplication factor changes to k ′ . The adjoint equation before the perturbation and theforward equation after the perturbation can be written as follows:A† φ†=A′ φ′=1 † †F φ ,k1 ′ ′F φ.
k′(7)
(8)

1 As shown in this equation, 1/k is multiplied to the operator F . Thus the operator F is a neutron generation operator by fissions and canbe explicitly presented as F φ = χgX

νΣf,g′ φg′ . Scattering is included to the operator A.

g′1

By multiplying φ′ to the both sides of Eq. (7) and φ† to the both sides of Eq. (8), and integrating them over whole phasespaces, the following equation can be derived:
φ′ , A† φ†=φ† , A′ φ′=

1 ′ † †φ ,F φ ,k1
φ† , F ′ φ′ .
k′(9)
(10)

When we consider the nature of the adjoint operators as shown in Eqs. (4) and (5), Eq. (9) can be rewritten asφ† , Aφ′ =1 †φ , F φ′ .
k(11)By substituting Eq. (11) from Eq. (10), the following equation is derived:
11 †1†′φ , ∆Aφ =φ† , F ′ φ′ +−φ , ∆F φ′ .
k′kk(12)From this equation, it can be found that reactivity inserted by this perturbation can be calculated by the followingequation:
1 †φ , ∆F φ′ − φ† , ∆Aφ′
11k.
(13)
ρ= − ′ =k khφ† , F ′ φ′ iThis equation suggests that small reactivity can be accurately calculated because reactivity is obtained from simple integralcalculations.
Furthermore, when we introduce an approximation φ′ ≈ φ to Eq. (13), the following equation is derived:
1 †φ , ∆F φ − φ† , ∆Aφk.
ρ=hφ† , F ′ φi(14)] The perturbation calculation presented in Eq. (13) is referred to as the rigorous perturbation and that in Eq. (14) is as
the first-order perturbation. In the rigorous perturbation calculations, adjoint neutron flux before the perturbation φ† andforward neutron flux after the perturbation φ′ are required, so calculations of φ′ should be done for each perturbation.
On the other hand, in the first-order perturbation calculations, forward neutron flux φ and adjoint neutron flux φ∗ before
the perturbation are required, so reactivity given by arbitrary perturbation can be calculated from these quantities. This
theory is generally used in sensitivity calculations.
One-group neutron diffusion eigenvalue adjoint equation can be written asdd1D(x) φ† (x) + Σa (x)φ† (x) = νΣf (x)φ† (x).
−(15)dxdxk
This is the exactly same as the one-group forward equation; φ(x) = φ† (x).
✓✏Problem 3: Calculate reactivity given by a perturbation in the problem 1 by using the first-order perturbation theory;
φ(x) and φ† (x)(= φ(x)). Compare results with those in problem 1.
✒✑
Perturbation is given only to Σa in the problem 1, so ∆F = 0.
✓✏

Problem 3: Calculate reactivity in the problem 1 by the rigorous perturbation theory; φ′ (x) and φ† (x)(= φ(x)).
✒✑
On the other hand, multi-group adjoint equation is different from multi-group forward equation. Multi-group adjoint
equation is presented as−ddxDg (x)d †φ (x)dx g

+ Σa,g (x)φ†g (x) +X

Σg→g′ φ†g (x)g′=XX1
Σg→g′ φ†g′ (x).
χg′ (x)φ†g′ (x) +νΣf,gk′′g(16)g

Numerical procedure to solve multi-group adjoint equation is basically the same as that for the forward equation, but anorder of energy group in the procedure is different. In two-group problems, fission neutron source is initially assumed,and then adjoint neutron flux is calculatedfor group 2, and then calculation is performed for group 1. Note that fission
X′′χg φg in the adjoint equation. Since computer program for the forward equation
neutron source in group g is νΣf,gg′2

might be extended to adjoint equation treatment, it is better that scalar fission source is defined as F =X

χg′ φg′ andg′fission source in each group is defined as νΣf,g F . Neutron multiplication factor at the nth outer iteration k n can beguessed by k n = F n /(F n−1 /k n−1 ), but is the initial fission source of the outer iteration is set as F 0 = 1 and k 0 = 1 isalso assumed, k n = F n since F n−1 /k n−1 = 1.
✓✏
Problem 5: Let us consider a one-dimensional slab with width of 100[cm] and constants are given in the table below.
External boundary conditions are zero neutron flux. Obtain k and φ†g by solving the adjoint equation, and confirm
that the calculated k coincides with k obtained by solving the forward equation of the same system.
Energy group
12

D
1.5
0.4

Σa0.010.1

Σg,g+10.02-νΣf0.005
0.141χ1.00.0

✒
✓

✑
✏

Problem 6: Obtain reactivity by the rigorous perturbation calculations when σa,1 is increased by 5% in a range from48[cm] to 52[cm] in the slab reactor given in the problem 5. Also calculate reactivity induced by different perturbation
which are given to νΣf,2 and Σ1→2 in the same spatial range.
✒✑✓✏
Problem 7: On the slab reactor given the problem 6, calculate reactivity when D1 and D2 are increased by 1% in awhole region.
✒✑
Reactivity calculations for diffusion coefficients perturbation are a bit complicated, so descriptions are provided below.
When a perturbation is given to diffusion coefficients, numerator of perturbation equation ρnume is calculated asfollows:
Zdφ′nume† d∆D(17)ρ=dx φdx
dx





Zdφ′
dφ†
dφ′
=φ† ∆D− dx
∆D.
(18)
dx
dx
dx
Since an integration in Eq. (17) is done over whole spatial region, so quantities of the first term of the right hand side ofEq. (18) should be on external boundaries, and those are zero. 2 Finally ρnume can be presented as † 
 ′ Zdφdφρnume = − dx
∆D(19)dx
dx
Neutron current J is represented as follows in the diffusion approximation:
J = −Ddφ.
dx

Thus Eq. (19) can be rewritten by using neutron current as † 
 ′ ZJJ
nume
ρ
= − dx∆D.
DD′(20)(21)

In general numerical procedure to diffusion equation, neutron current j on the mesh surface is defined from neutron flux
at mesh surface and mesh center, so J should be constant in left-half and right-half in each mesh. When these are denoted
to as Jlef t and Jright , ρnume can be written asρnume = −∆D ∆xi  ††′′JJ+JJ
leftrightleftrightDD′ 2

2 In the zero neutron flux condition, φ† = 0 on the boundary. In the reflective boundary condition, dφ′ /dx = 0 on the boundary.

3(22)