- Email: [email protected]
The optimum design of power distribution for pressurized water reactor
Annals of Nuclear Energy 50 (2012) 126–132
Contents lists available at SciVerse ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
The optimum design of power distribution for pressurized water reactor Chunhui Dai ⇑, Xinyu Wei, Yun Tai, Fuyu Zhao School of Energy and Power Engineering, Xi’an Jiaotong University, No. 28, Xianning West Road, Xi’an, Shaanxi 710049, PR China
a r t i c l e
i n f o
Article history: Received 13 September 2011 Received in revised form 25 July 2012 Accepted 29 July 2012 Available online 21 September 2012 Keywords: PWR Power distribution Two-level optimization
a b s t r a c t The aim of this work is to develop a two-level optimization method for designing the optimum initial fuel loading pattern and burnable poison placement in pressurized water reactors. At the lower level, based on the fuel loading pattern (LP) optimized by backward diffusion calculation theory, Pontryagin’s maximum principle is employed to investigate the optimum arrangement of burnable poison (BP) that can generate the lowest radial power peaking factor (PPF). At the upper level a multi-objective problem (MOP), with LP and BP as two objective functions, is proposed by coordinate the interrelationship of LP and BP, and optimized by non-dominated sorting genetic algorithm (NSGA-II). The results of optimum designs called ‘Pareto optimum solutions’ are a set of multiple optimum solutions. After sensitivity analysis is performed, the final optimum solution which is chosen based on a typical VVER-1000 reactor reveals that the method could not only save the fuel consumption but also reduce the PPF in comparison to published data. Ó 2012 Elsevier Ltd. All rights reserved.
1. Introduction The average output power of a reactor with a given volume is inversely proportional to the radial power peaking factor and directly proportional to the maximum power density. To achieve the minimum PPF, the fuel assemblies are usually placed in different areas of the core, with BPs reasonably arranged. The optimization of LP and BP is difficult in reactor design, but several techniques have been developed to facilitate this process. Haling (1964), Kim et al. (1987), and Suh and Levine (1990) decoupled fuel and BP to determine their placement; stochastic algorithms such as genetic algorithm, simulated annealing were used by Chapot et al. (1999) and Kropaczek and Turinsky (1991) to implement the global optimized calculation for searching the optimum arrangement of fuel LP and BP directly. Most of these methods require long computation time. Chao et al. (1986) obtained the theoretically power distributions with the backward diffusion calculation, which could reduce the calculation time. However, it is only applicable to the LP optimization without paying specific attention to BP. In this paper the typical Russian PWR VVER-1000 is taken as the subject to study the one-dimensional optimum initial design parameters of the reactor core. A two-level optimum method is developed and implemented successfully for the multi-objective optimization of power distribution, based on minimizing the total fuel potential required for certain reactivity and obtaining the lowest PPF subject to the constraint condition of the maximum power density and marginal neutron flux of the reactor core. At the lower ⇑ Corresponding author. E-mail address: [email protected] (C. Dai). 0306-4549/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.anucene.2012.07.017
level optimization, backward diffusion calculation theory (Chao et al., 1986) and Pontryagin’s maximum principle (Pontryagin et al., 1962) are individually applied to get the optimum solutions of LP and BP divided from the optimum power distribution problem. On this basis, the MOP with LP and BP treated as two objective functions is proposed and computed by NSGA-II method (Deb et al., 2002) to obtain the optimum designs. Because the results of the optimum designs ‘Pareto optimum solutions are a set of multiple optimum solutions, sensitivity analysis is performed so that the final core design scheme could be chosen and reported. The optimum results can provide theoretical references for improving the initial designs of LP and BP in reactor core. 2. The VVER-1000 reactor core The reactor core of a typical VVER-1000 consists of 163 hexagonal fuel assemblies with equal dimensions. The reactor core includes six different enrichments (FSAR, 2003). Table 1 and Fig. 1 show specifications of fuel assemblies and their typical arrangement in first cycle respectively. 3. Optimization of core fuel LP 3.1. Backward diffusion calculation theory The two group diffusion equations are written as follows:
D1 r2 /1 þ R1 /1 ¼
1 ðm1 Rf 1 /1 þ m2 Rf 2 /2 Þ keff
ð1Þ
127
C. Dai et al. / Annals of Nuclear Energy 50 (2012) 126–132 Table 1 Main description of fuel assemblies. Assembly type
Number of fuel rods (enrichment/%)
Number of BP rods
Density of BP (kg/m3)
16 24 36 24B20 24B36 36B36
311(1.6) 311(2.4) 245(3.7) + 66(3.3) 311(2.4) 311(2.4) 245(3.7) + 66(3.3)
0 0 0 18 18 18
0 0 0 20 36 36
where r is made a dimensionless quantity by dividing actual radial distance (R) by L1. Using Green’s function (Chao et al., 1983), K⁄ could be obtained
PðrÞ
K ðrÞ ¼
/0 II00ðrðrÞ0 Þ þ A1 ðrÞ þ A2 ðrÞ A3 ðrÞ
A1 ðrÞ ¼ K 0 ðrÞ
Z
r
Pðr 0 ÞI0 ðr 0 Þr 0 dr
0
ð8Þ
ð9Þ
0
A2 ðrÞ ¼ I0 ðrÞ
Z
r0
Pðr 0 ÞK 0 ðr 0 Þr 0 dr
0
ð10Þ
r
36
36B36
24B36
24B20
16
16
24
24
16
24
16
24
24B36 16
16
36
24B36
24 16
36
r0 ¼
36B36
v 1 Rf 1 v 2 Rf 2 R12 keff R1
þ
keff R2 R1
k
e ¼ v 2 Rf 21R
ð2Þ
ð3Þ
ð4Þ
12
keff R2 R1
Referring to the previous studies (Tran et al., 2009), the ‘‘effective’’ fast neutron diffusion equation can be written as:
L21 r2 /1 /1 ¼ K /1
ð5Þ
where L1 is the fast diffusion length, which is assumed to be constant as it vary mildly over the entire core; depending only on the group constants, K⁄ is defined as the fuel potential for evaluating the fuel utilization (Chao et al., 1986) L22
K ¼ k1
1 þ eL 2 1
L2
0
ð11Þ
0
Pðr0 ÞI0 ðr 0 Þr 0 dr h i ð1b1 ÞL1 0Þ r 0 I0 ðr 0 Þ II10 ðr þ 2ð1þb ðr 0 Þ ÞD 1 1 0
ð12Þ
R0 L1
ð13Þ
3.2. Objective functions and constrains conditions
where D1, D2 are the fast and thermal diffusion coefficient; u1, u2 are the fast and thermal neutron flux; v1Rf1, v2Rf2 are the fast and thermal fission cross section; R1, R2 are the fast and thermal removal cross section; R1,2 is the transfer cross section from fast to thermal neutron; keff is the effective multiplication factor which is assumed to be unity. In the design of initial fuel LP, keff value should be larger than 1 to guarantee adequate excess reactivity. The infinite multiplication factor of the whole reactor core k1 and the fast fission factor e can be written as:
k1 ¼
Pðr 0 ÞI0 ðr 0 Þr 0 dr
0
where R0 is the equivalent radius of the reactor core.
Fig. 1. Typical fuel assembly arrangement of VVER-1000.
D2 r2 /2 þ R2 /2 ¼ R1;2 /1
r0
where P(r) is one-dimensional radial power distribution; K0 is the zeroth order modified Bessel function of the first kind; I0 is the zeroth order modified Bessel function of the second kind; b1 is the albedo of fast neutron; r0 is
36
24B36
16
Z
R r0
/0 ¼
36
24B36 16
24B36
16 24
16
I0 ðrÞK 0 ðr 0 Þ I0 ðr 0 Þ
36
24
16
A3 ðrÞ ¼
ð6Þ
1 þ eL22 k1
As shown in Eq. (8), the fuel potential K⁄ is a function of the power distribution. Integrating K⁄ over the space and divide it by the area of core, we get the average total fuel potential, denoted as Fp
R r0 Fp ¼
0
2prK ðrÞdr pr20
ð14Þ
Since Fp is the average total potential relevant to keff, the LP optimization can be obtained by seeking proper P(r) to minimize Fp under the power peaking constraint condition 0 6 PðrÞ 6 Plim : For simplicity, the power distribution P(r) is set as a third order polynomial and normalized by Eq. (16)
PðrÞ ¼ a0 þ a1 r þ a2 r 2 þ a3 r 3 R r0 2prpðrÞdr a0 a2 r 20 a3 r 30 0 ¼ 2 þ þ ¼1 2 4 5 pr20
ð15Þ ð16Þ
where a0, a1, a2, a3 are coefficients which indicate the power distribution. In the center of core, the neutron is total reflection, so
P0ðrÞ jr¼0 ¼ a1 ¼ 0
ð17Þ
P(r)should not be negative in the overall core, and suppose P(r) get maximum at rm = m r0,
P0 ðrÞjr¼rm ¼ a1 þ 2a2 þ 3rm a3 ¼ 0
ð18Þ
a0 þ r 2m a1 þ r 2m a2 þ r 3m a3 ¼ Plim
ð19Þ
3.3. Calculation and analysis
1
where L2 is the thermal diffusion length. In cylindrical coordinates, Eq. (5) could be written as 2
d /1 ðrÞ dr
2
þ
1 d / ðrÞ /1 ðrÞ ¼ K ðrÞ/1 ðrÞ r dr 1
ð7Þ
The admissible value of Plim is around 1.3 (Boroushaki et al., 2003; Mazrou and Hamadouche, 2006). Now P(r) could be obtained through calculating Eq. (20) which is combined by Eqs. (16), (18), and (19).
CA¼J
ð20Þ
128
C. Dai et al. / Annals of Nuclear Energy 50 (2012) 126–132
8 _ wðtÞ ¼ @H > @X > < _ XðtÞ ¼ @H @w > > _ _ : _ _ _ HðX ; U ; w; tÞ P HðX ; U; w; tÞ
ð25Þ
where ^ denotes optimum quantities in admissible region. In the g l ðXÞ 6 0 Eq. (23) is replaced case X(t) and U(t) are all constrained qj ðUÞ 6 0 by more general condition (Zhao and Ren, 2009):
@H ¼0 @U
ð26Þ
with
H ¼ H þ
n s X X ki ðtÞpi ðX; UÞ þ lj ðtÞqj ðUÞ i¼1
where ki ðtÞ and lj(t) are Lagrange factors, pi(X, U) is the lowest order derivative of gi(X) with respect to t in the derivatives of containing explicit functions U(t). If gi(X) < 0, ki ðtÞ ¼ 0 satisfy, optimum trajectories do not be in the constrained boundary gi(X) = 0; if qj(U) < 0, lj(t) satisfy, admissible control U(t) do not lie in the con strained boundary qj(U) = 0; if qj(U) = 0 and lj(t) – 0, admissible control U(t) lie in qj(U) = 0.
Fig. 2. The relationship between Fp and m.
where J ¼ ½1 P lim 0T , A ¼ ½a0 a2 a3 T
2
1 r 20 =4 r 30 =5
6 C ¼ 41
r 2m
0
2
3
7 r 2m 5 3r m
4.2. Modeling and calculation
m is the most important parameter as it determines the peaking point of the power distribution, and A is computed by giving m in its allowable domain (0 6 m 6 1). The relationship of Fp and m is shown in Fig. 2. Fp reaches to the minimum value when the value of m is between 0.5 and 0.6, neither m = 0 nor m = 1.0. Some values of Fp are calculated with different m, listed in Table 2. 4. Optimization of BP Burnable poisons should be arranged in the core to compensate the excess reactivity and flatten the power distribution shape on the basis of the optimum fuel LP obtained in Section 3. Using the single group diffusion theory, the optimization of BP is carried out under the constraint condition of the maximum power density. 4.1. Pontryagin’s maximum principle (Pontryagin et al., 1962) Let the plant
X_ ¼ f ðX; U; tÞ
ð21Þ
Single group diffusion equation is
r Dr/ þ Ra / þ Rc / ¼
Z
mRf /
ð28Þ
keff
where Ra is the absorption cross section, vRf is the fission cross section, keff is the effective multiplication factor and D is the diffusion coefficient. A variable denoted as Rc is introduced to represent the burnable poisons. Assuming that there is no leakage in the axial direction of core, we can get
D 0 mRf / D/00 þ D0 þ Ra / Rc / ¼ 0 / þ R keff
ð29Þ
Let x1(R) = u(R), x2(R) = u0 (R), x3(R) = R, Eq. (28) could be rewritten as
8 dx1 ðRÞ > ¼ x2 ðRÞ > > < dR 0 dx2 ðRÞ ðRÞ ¼ DDðRÞ þ x31ðRÞ x2 ðRÞ þ BðRÞx1 ðRÞuðRÞ BðRÞx1 ðRÞ dR > > > dx3 ðRÞ : ¼1 dR
Have an associated cost index of
Max Jðt 0 Þ ¼ vðXðT 0 Þ; T 0 Þ þ
ð27Þ
j¼1
ð30Þ
mRf ðRÞ
T0
LðX; U; tÞdt
ð22Þ
t0
where X ¼ ½x1 x2 . . . xn1 T , U ¼ ½u1 u2 un2 T . If the control is unconstrained the condition for optimality is
@H ¼0 @U
ð23Þ
where BðRÞ ¼
keff
Ra ðRÞ
DðRÞ
Rc ðRÞ , u(R) is control vector uðRÞ ¼ mRf ðRÞ keff
the neutron economy point of view, low neutron leakage is expected, which requires low the edge neutron flux u0, so the boundary condition is:
x1 ðR0 Þ ¼ /0
ð31Þ
x2 ð0Þ ¼ 0
where H is Hamiltonian function
. From
Ra ðRÞ
The power density is
H ¼ LðX; U; tÞ þ wT f ðX; U; tÞ
ð24Þ
where w is Lagrange factor, when the control is constrained, Eq. (23) changes into
Pv ðRÞ ¼ F a Ef Rf ðRÞ/ðRÞ
ð32Þ
where Fa is the ratio of the heat of core (mainly fuel element and coolant) to the total heat of the reactor, Ef is the energy released
Table 2 Some typical values of m and Fp. m Fp
0.0 1.088445
0.1 1.085443
0.3 1.073669
0.5 1.048841
0.7 1.068476
0.9 1.093214
1.0 1.107083
129
C. Dai et al. / Annals of Nuclear Energy 50 (2012) 126–132
per fission, and it equal to about 200 MeV. The constraint condition is a. the cross section of control material should not be negative, namely the control domain
qðuÞ ¼ uðRÞ 6 0
ð33Þ
b. the power density should not larger than the maximum power density Pvmax, so that
gðXÞ ¼ F a Ef Rf x1 Pv max 6 0
ð34Þ
The radial PPF which is expected as low as possible could be obtained through dividing the Pvmax by the average power density of the reactor core. As the total area of the core is constant, the minimizing PPF process is then translated to searching the optimum admissible control u(R) to obtain the maximum output power Pt computed by integrating Eq. (32) over the whole core, so
Max Pt ¼ 2pF a Ef ¼ 2pF a Ef
Z Z
R0
from which we can get Rc = 0, that means the BP is needless in this region. ii. l (R) = 0, with the same foundation as case ‘i’, we can get
@H ¼ w2 Bx1 ¼ 0 @u
ð44Þ
Because neither B(R) nor x1(R) is 0, the last variable w2(R) must equal 0. Substituting w2(R) = 0 into Eq. (37), the adjoint equation can be written as
8 dw 1 > ¼ Rf x3 > < dR dw2 ¼ w1 dR > > : dw3 ¼ Rf x1 dR
ð45Þ
As it is impossible for practical engineering situation that dw1 ¼ Rf x3 ¼ 0, case ii does not exist. When g(X) = 0, k – 0, generaldR ized Hamilton is
Rf ðRÞ /ðRÞ R dR 0
H ¼ H þ kðRÞg 00 ðXÞ
R0
Rf ðRÞx1 ðRÞx3 ðRÞdR
ð35Þ
0
4.3. Calculation and analysis
ð46Þ
the optimum condition is 8 < @H ¼ w2 Bx1 þ kRf Bx1 ¼ 0 @u n h io 00 : g ðXÞ ¼ F a Ef 2x2 R0f þ x1 R00f þ Rf 1 þ D0 x2 þ Bx1 uðRÞ Bx1 ¼ 0 x3 D
ð47Þ
Hamilton function is
0 D 1 x2 þ Bx1 uðRÞ B2 x1 þ w3 H ¼ Rf x1 x3 þ w1 x2 þ w2 þ D x3
after calculation, we can obtain
ð36Þ
f 1
dR
f
ð37Þ
x23
The derivatives of g(X) are
g 0 ðXÞ ¼
dgðXÞ ¼ F a Ef ðRf x2 þ x1 R0f Þ dR
ð38Þ
keff
d gðXÞ 2
dR
1 D0 ¼ F a Ef 2x2 R0f þ x1 R00f þ Rf þ x2 þ Bx1 uðRÞ Bx1 x3 D ð39Þ It can be found that g0 (X) does not contains the admissible control vector u(R) explicitly while g 00 ðXÞ does. The augmented Hamilton function is
H ¼ H þ kðRÞg 00 ðXÞ þ lðRÞqðuÞ
ð40Þ
Analysis on the optimum solutions in different conditions is then carried out, with g(X) < 0, kðRÞ ¼ 0, Eq. (40) changes into
H ¼ H þ lðRÞqðuÞ
ð41Þ
Here, two cases are possible for optimum value of u(R), ( i. l (R) – 0, there are the relations in accordance with Eq. (27) @H @u
¼ @H l¼0 @u
qðuÞ ¼ 0
ð42Þ
The compute result is
lðRÞ ¼ w2 Bx1 uðRÞ ¼ 0
02
00
f
f
0
ð48Þ
0
ð43Þ
f
when g(X) = 0, the optimum trajectories of power are in the constraints boundary and the neutron and BP distribution could be obtained from Eq. (48); while if g(X) < 0, the BP is needless and Eq. (33) could provides the optimum power distribution. As keff is unknown, a method so-called ‘‘source iteration’’ is used as follows. The set of ð0Þ initial value keff and (mRfu)(0) makes the fission source mRfu/keff in Eq. (28) known, which could give the solution of (mRfu)(1). Then according to the definition of keff ðnþ1Þ
2
g 00 ðXÞ ¼
0
R 2R R DR > : uðRÞ ¼ 1 BRf R þ BRf2 BRf DBRf
The adjoint equation is
8 dw1 > ¼ Rf x3 w2 ½B þ uðRÞB > dR > < dw 0 2 ¼ w1 þ w2 ðx13 þ DD Þ dR > > > dw3 ¼ R x w2 x2 :
8 P v max > < x1 ðRÞ ¼ F a Ef Rf ðRÞ
R R0 ¼
0 1 ðnÞ keff
mRf ðRÞx1ðnþ1Þ ðRÞ2pRdR mRf ðRÞxðnÞ 1 ðRÞ2pRdR 0
ð49Þ
R R0
ð1Þ
we can update the effective multiplication keff ; so that X can be calculated again. Repeat this process until the convergence criterion satisfy: (a) characteristic value convergence criterion j
e1 ; and (b) neutron flux convergence criterion j
ðnÞ ðn1Þ k eff eff ðnÞ k eff
k
ðnÞ ðn1Þ x1 ðRÞx1 ðRÞ ðnÞ x1 ðRÞ
j < e2 ,
–4
e1 and e2 are small positive value, for example e1 = 10 e2 = 10–3.
j<
and
5. Multi objective optimization of LP and BP The optimum fuel LP and BP are obtained by minimizing Fp and PPF (namely maximizing Pt) respectively in Sections 3 and 4. However LP and BP are coupled in the power distribution optimum calculation, the two objective functions cannot be optimum simultaneously. With Fp reaches to the minimum value, PPF may not be the minimum as the neutron flux at the periphery of core is lower than its limit. On the other hand, the optimum LP required by the objective PPF cannot always satisfy the minimizing of Fp. As a result, a multi objective model is applied to get an eclectic solution in this study.
130
C. Dai et al. / Annals of Nuclear Energy 50 (2012) 126–132
5.1. Multi objective optimization A multi-objective problem (MOP) consists of optimizing (i.e., maximizing or minimizing) several objectives simultaneously, with a number of constraints. MOP can be described as follows (Sanaye and Hajabdollahi, 2010):
Max F i ¼ ðXÞ;
i ¼ 1; 2; . . . ; I
ð50Þ
where X is decision vector
X ¼ ðx1 ; x2 ; . . . ; xn Þ
ð51Þ
The constraint conditions are
g l ðXÞ 6 0 l ¼ 1; 2; . . . ; L
ð52Þ
hj ðXÞ ¼ 0 j ¼ 1; 2; . . . ; s
ð53Þ
Usually, it is impossible to maximize or minimize all the objective function at the same time because of the interaction between them. Therefore the results of optimum designs called ‘‘Pareto optimum solutions’’ are a set of multiple optimum solutions which are compromised by giving different focus to every objective. The genetic algorithm (GA) is a semi-stochastic method based on Darwin’s laws of natural selection. To solve MOP, a non-dominated sorting genetic algorithm which is used in this research was proposed by Srinivas and Deb (1994) (NSGA) and modified by Deb et al. (2002) (NSGA-II).
Fig. 3. The distribution of Pareto optimum points solutions.
Table 3 The optimum values of objective functions for points A–E in Pareto-optimal fronts.
m Fp PPF
O
A
B
C
D
E
0.0 1.088445 1.300000
0.4859 1.050267 1.195738
0.5082 1.047541 1.196277
0.5394 1.044717 1.196919
0.5561 1.043975 1.200540
0.5633 1.043893 1.204250
5.2. Objective functions and constraints Fp and PPF are denoted as the two objectives respectively in the MOP, and as m has an effect on Fp which acts as the basis of the optimization of PPF, it is regarded as the decision variable, subject to
06m61
ð54Þ
Now the MOP is described as
Min s:t:
F 1 ðmÞ ¼ F p F 2 ðmÞ ¼ PPF
m60
ð55Þ
m160 5.3. Multi objective optimum results Using a search population size of M = 200 individuals, crossover probability of pc = 0.9, gene mutation probability of pm = 0.035, the NSGA-II method is performed for 100 generations, the results for Pareto-optimum curve are shown in Fig. 3, which denote the conflict between Fp and PPF clearly. Any change that increases Fp which indicates the fuel consumption leads to a decrease of PPF, and vice versa. Five typical points (A–E), for which the optimum values of Fp and PPF are listed in Table 3, are chosen as object points to study the relationship between optimum Fp and PPF. The minimum Fp exists at design point E (1.043893), where PPF is the highest (1.204250). On the other hand the maximum Fp occurs at point A (1.050267), with a lowest PPF (1.195738). Design points A and E are the optimum situations at which PPF and Fp are single objective function respectively. The distribution of m for Pareto front is shown in Fig. 4, the lower and upper bounds of m are 0 and 1 respectively. m has its optimum values in the range of 0.4859 < m < 0.5633, Fp and PPF in this margin are all lower, compared with point O where m = 0. It means the solutions are optimized successfully by NSGA-II method. The selection of final optimization among the optimum points existing on the Pareto front needs a process of decision-making
Fig. 4. Scattering of m for the Pareto optimum front.
based on the experiences of engineer. In the process, the focus should be given to the importance of each objective function. If Fp is considered to be more important, the points from D to E could be picked up as the final solutions, while points during A and B should be chosen when paying more attentions to PPF. In this study, the importance of the two objective functions is thought to be equivalent, so the points around C are regarded as the final optimum solution. Take typical point C where m = 0.5394 as the discussion object, the non-uniform distribution of the fuel potential K⁄ is shown in Fig. 5. K⁄ firstly increases until P equals Plim, then K⁄ will decrease with the increasing of radial relative position (RRP) in the interval from 0.81 to 1. The minimum K⁄ appears at the edge of core. The relationship between the fuel enrichment and K⁄, which could be utilized to provide distribution of fuel enrichment against the ra-
C. Dai et al. / Annals of Nuclear Energy 50 (2012) 126–132
Fig. 5. The optimum distributions K⁄.
Fig. 7. Distribution of control vector (u(R)).
Fig. 6. The optimum enrichment (U235) distribution.
Fig. 8. Distribution of the optimum power density.
dial relative position, is calculated by Eq. (6). Then, the theoretical relationship between U-235 enrichment and percentage radial distance is shown in Fig. 6. Fig. 6 shows the non-uniform distribution of the fuel optimum enrichment. The enrichment increases when RRP < 0.81, after that it decreases until the edge of core where the enrichment is the lowest. From Figs. 5 and 6, it can be concluded that this LP scheme could not only guarantee that power in the center of core relatively low, but also could make sure the edge neutron flux is not too high as leading to a high neutron leakage. It should be emphasized that the optimal LP obtained here is known as the ‘low leakage loading pattern’. Compared with value which the power peaking is at the center of core (m = 0), Fp and PPF of point C are lower, that means the optimum LP and BP saved the fuel consumption about 4.2% and reduced the radial power peaking factor about 7.93%. As shown in Fig. 7, when RRP < 0.38 or RRP > 0.74, the control vector u(R) = 0 and as a result, the BP is needless; while 0.38 < RRP < 0.74, BP must be arranged to flat the power density peaking caused by the high enrichment fuel. As shown in Fig. 8, power density increases with the increasing of RRP because the enrichment and fission cross section of the fuel rises. When RRP = 0.38, the Pv reaches the maximum value and keeps the maximum value until RRP = 0.74 as the arrangement of BP. When
131
Fig. 9. Distribution of the optimum neutron flux.
0.74 < RRP < 1, Pv decreases and reaches minimum at the edge of core. The normalized neutron distribution in Fig. 9 shows the similar trend as the power distribution in Fig. 6. However, when
132
C. Dai et al. / Annals of Nuclear Energy 50 (2012) 126–132
Table 4 Comparing Fp and PPF with different design for VVER-1000.
VVER-1000 designer Perturbation theory Point C
Fp
PPF
1.04385 1.04391 1.044717
1.22 1.215 1.196919
0.38 < RRP < 0.74, the neutron flux decreases as the optimum trajectory lie in the constraint boundary and the power density is constant while the enrichment increases tardily. There is a relation / (R) = Pvmax/(FaEfRf(R)) in this region, so, obviously, / will decrease with RRP. 5.4. Results and discussion An optimization process for the first operating cycle of a typical VVER-1000 reactor core is performed. The comparison of final optimum results of point C with the references (FSAR, 2003; Hosseini and Vosoughi, 2012) reported load pattern are shown in Table 4. Although the Fp is larger than the two references (about 0.083% and 0.077%), the PPF is reduced by 1.89% and 1.49%, compared with the original VVER-1000 design and perturbation theory optimization respectively. This is result from the coordinating the interaction between Fp and PPF during chosen process of the final optimum design. 6. Conclusions A two-level optimization method is applied to design the initial optimum LP and BP arrangement of VVER-1000. In the lower level optimization, Fp is firstly minimized by searching the optimum LP with the backward diffusion calculation theory, the result shows that when m is in the margin of 0.5 < m < 0.6, Fp is the lowest. On the basis of fuel distribution, BP is optimized to get the minimum radial power peaking factor. Then, considering the interrelationship of each other, an in-depth study which acts as the upper level optimization is carried out to obtain the optimum LP and BP with NSGA-II. A set of Pareto optimum points in the optimum region (0.4859 < m < 0.5633) which clearly revealed the conflict between the objectives is obtained and showed. After the analysis of the relationship of Pareto solutions and weighting of LP and BP, the points around point C are chosen and deemed to be the final optimum solutions. Point C (m = 0.5394) is taken as the typical result for discussion. The absorbing material should be arranged from 0.38 to 0.74 time of core equivalent radius. The calculation shows that not only the fuel consumption is saved about 4.2%, but also the PPF is reduced by 7.93%, compared with a uniform loading pattern core. In the comparison to the
VVER-1000 original design and some other optimization, because of the coordinating the interaction between Fp and PPF during deciding process of the final optimum design, PPF is reduced by 1.89% and 1.49% respectively, however the Fp is larger (about 0.083% and 0.077%). The application of these methods and the optimum solutions will provide theoretical references for improving the initial designs of LP and BP in reactor core. Acknowledgment This research is funded by National Natural Science Foundation of China (10975112). References Boroushaki, M., Ghofrani, M.B., Lucas, C., 2003. Optimum fuel core loading pattern design in PWR nuclear power reactors using genetic algorithm and fuzzy nonlinear programming. J. Intell. Fuzzy Syst. 14, 85–93. Chao, Y.A., Suo, C.A., Penkrot, J.A., 1983. On the theory of interface flux nodal method. In: Proc, CNS/ANS Int. Conf. Numerical Methods in Nuclear Engineering, Montreal, p. 307. Chao, Y.A., Hu, C.W., Suo, C.A., 1986. A theory of fuel management via backward diffusion calculation. Nucl. Sci. Eng. 93, 78–87. Chapot, J.L.C., Silva, F.C., Schirru, R., 1999. A new approach to the use of genetic algorithms to solve the pressurized water reactor’s fuel management optimization problem. Ann. Nucl. Energy 26, 641–655. Deb, K., Pratap, A., Agarwal, S., Meyarivan, T., 2002. A fast and elitist multi objective genetic algorithm: NSGA-II. IEEE Trans. Evol. Comput. 6, 182–197. Final Safety Analysis Report of BNPP-1, 2003. Reactor, 49.BU.1 0. 0. OO.FSAR. RDR001, Rev. 0 (Chapter 4). Haling, R.K., 1964. Operating Strategy for Maintaining Optimum Power Distribution throughout Life. TID 7672. Hosseini, M., Vosoughi, N., 2012. Development of a VVER-1000 core loading pattern optimization program based on perturbation theory. Ann. Nucl. Eng. 39, 35–41. Kim, Y.J., Downar, T.J., Sesonske, A., 1987. Optimization of core reload design for low-leakage fuel management in pressurized water reactor. Nucl. Sci. Eng. 96, 85–101. Kropaczek, D.J., Turinsky, P.J., 1991. In-core nuclear fuel management optimization for pressurized water reactors utilizing simulated annealing. Nucl. Technol. 95, 9–31. Mazrou, H., Hamadouche, M., 2006. Development of a supporting tool for optimum fuel management in research reactors using artificial neural networks. Nucl. Eng. Des. 236, 255–266. Pontryagin, L.S., Boltyanskii, V.G., Gamkrelidze, R.V., Mischenko, E.F., 1962. The Mathematical Theory of Optimum Process. Wiley-Interscience Publishers, New York. Suh, J.S., Levine, S.H., 1990. Optimized automatic reload program for pressurized water reactors using simple direct optimization techniques. Nucl. Sci. Eng. 105, 371–382. Srinivas, N., Deb, K., 1994. Multi objective optimization using non-dominated sorting in genetic algorithms. J. Evol. Comput. 2, 221–248. Sanaye, S., Hajabdollahi, H., 2010. Multi-objective optimization of shell and tube heat exchangers. Appl. Therm. Eng. 30, 1937–1945. Tran, H.N., Yammoto, A., Yamane, Y., 2009. An improved inverse analysis model for fuel loading pattern optimization. J. Nucl. Sci. Technol. 46, 1162–1169. Zhao, F.Y., Ren, X.R., 2009. Optimum performance analysis of a sodium cold trap. Nucl. Eng. Des. 239, 490–494.
Contents lists available at SciVerse ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
The optimum design of power distribution for pressurized water reactor Chunhui Dai ⇑, Xinyu Wei, Yun Tai, Fuyu Zhao School of Energy and Power Engineering, Xi’an Jiaotong University, No. 28, Xianning West Road, Xi’an, Shaanxi 710049, PR China
a r t i c l e
i n f o
Article history: Received 13 September 2011 Received in revised form 25 July 2012 Accepted 29 July 2012 Available online 21 September 2012 Keywords: PWR Power distribution Two-level optimization
a b s t r a c t The aim of this work is to develop a two-level optimization method for designing the optimum initial fuel loading pattern and burnable poison placement in pressurized water reactors. At the lower level, based on the fuel loading pattern (LP) optimized by backward diffusion calculation theory, Pontryagin’s maximum principle is employed to investigate the optimum arrangement of burnable poison (BP) that can generate the lowest radial power peaking factor (PPF). At the upper level a multi-objective problem (MOP), with LP and BP as two objective functions, is proposed by coordinate the interrelationship of LP and BP, and optimized by non-dominated sorting genetic algorithm (NSGA-II). The results of optimum designs called ‘Pareto optimum solutions’ are a set of multiple optimum solutions. After sensitivity analysis is performed, the final optimum solution which is chosen based on a typical VVER-1000 reactor reveals that the method could not only save the fuel consumption but also reduce the PPF in comparison to published data. Ó 2012 Elsevier Ltd. All rights reserved.
1. Introduction The average output power of a reactor with a given volume is inversely proportional to the radial power peaking factor and directly proportional to the maximum power density. To achieve the minimum PPF, the fuel assemblies are usually placed in different areas of the core, with BPs reasonably arranged. The optimization of LP and BP is difficult in reactor design, but several techniques have been developed to facilitate this process. Haling (1964), Kim et al. (1987), and Suh and Levine (1990) decoupled fuel and BP to determine their placement; stochastic algorithms such as genetic algorithm, simulated annealing were used by Chapot et al. (1999) and Kropaczek and Turinsky (1991) to implement the global optimized calculation for searching the optimum arrangement of fuel LP and BP directly. Most of these methods require long computation time. Chao et al. (1986) obtained the theoretically power distributions with the backward diffusion calculation, which could reduce the calculation time. However, it is only applicable to the LP optimization without paying specific attention to BP. In this paper the typical Russian PWR VVER-1000 is taken as the subject to study the one-dimensional optimum initial design parameters of the reactor core. A two-level optimum method is developed and implemented successfully for the multi-objective optimization of power distribution, based on minimizing the total fuel potential required for certain reactivity and obtaining the lowest PPF subject to the constraint condition of the maximum power density and marginal neutron flux of the reactor core. At the lower ⇑ Corresponding author. E-mail address: [email protected] (C. Dai). 0306-4549/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.anucene.2012.07.017
level optimization, backward diffusion calculation theory (Chao et al., 1986) and Pontryagin’s maximum principle (Pontryagin et al., 1962) are individually applied to get the optimum solutions of LP and BP divided from the optimum power distribution problem. On this basis, the MOP with LP and BP treated as two objective functions is proposed and computed by NSGA-II method (Deb et al., 2002) to obtain the optimum designs. Because the results of the optimum designs ‘Pareto optimum solutions are a set of multiple optimum solutions, sensitivity analysis is performed so that the final core design scheme could be chosen and reported. The optimum results can provide theoretical references for improving the initial designs of LP and BP in reactor core. 2. The VVER-1000 reactor core The reactor core of a typical VVER-1000 consists of 163 hexagonal fuel assemblies with equal dimensions. The reactor core includes six different enrichments (FSAR, 2003). Table 1 and Fig. 1 show specifications of fuel assemblies and their typical arrangement in first cycle respectively. 3. Optimization of core fuel LP 3.1. Backward diffusion calculation theory The two group diffusion equations are written as follows:
D1 r2 /1 þ R1 /1 ¼
1 ðm1 Rf 1 /1 þ m2 Rf 2 /2 Þ keff
ð1Þ
127
C. Dai et al. / Annals of Nuclear Energy 50 (2012) 126–132 Table 1 Main description of fuel assemblies. Assembly type
Number of fuel rods (enrichment/%)
Number of BP rods
Density of BP (kg/m3)
16 24 36 24B20 24B36 36B36
311(1.6) 311(2.4) 245(3.7) + 66(3.3) 311(2.4) 311(2.4) 245(3.7) + 66(3.3)
0 0 0 18 18 18
0 0 0 20 36 36
where r is made a dimensionless quantity by dividing actual radial distance (R) by L1. Using Green’s function (Chao et al., 1983), K⁄ could be obtained
PðrÞ
K ðrÞ ¼
/0 II00ðrðrÞ0 Þ þ A1 ðrÞ þ A2 ðrÞ A3 ðrÞ
A1 ðrÞ ¼ K 0 ðrÞ
Z
r
Pðr 0 ÞI0 ðr 0 Þr 0 dr
0
ð8Þ
ð9Þ
0
A2 ðrÞ ¼ I0 ðrÞ
Z
r0
Pðr 0 ÞK 0 ðr 0 Þr 0 dr
0
ð10Þ
r
36
36B36
24B36
24B20
16
16
24
24
16
24
16
24
24B36 16
16
36
24B36
24 16
36
r0 ¼
36B36
v 1 Rf 1 v 2 Rf 2 R12 keff R1
þ
keff R2 R1
k
e ¼ v 2 Rf 21R
ð2Þ
ð3Þ
ð4Þ
12
keff R2 R1
Referring to the previous studies (Tran et al., 2009), the ‘‘effective’’ fast neutron diffusion equation can be written as:
L21 r2 /1 /1 ¼ K /1
ð5Þ
where L1 is the fast diffusion length, which is assumed to be constant as it vary mildly over the entire core; depending only on the group constants, K⁄ is defined as the fuel potential for evaluating the fuel utilization (Chao et al., 1986) L22
K ¼ k1
1 þ eL 2 1
L2
0
ð11Þ
0
Pðr0 ÞI0 ðr 0 Þr 0 dr h i ð1b1 ÞL1 0Þ r 0 I0 ðr 0 Þ II10 ðr þ 2ð1þb ðr 0 Þ ÞD 1 1 0
ð12Þ
R0 L1
ð13Þ
3.2. Objective functions and constrains conditions
where D1, D2 are the fast and thermal diffusion coefficient; u1, u2 are the fast and thermal neutron flux; v1Rf1, v2Rf2 are the fast and thermal fission cross section; R1, R2 are the fast and thermal removal cross section; R1,2 is the transfer cross section from fast to thermal neutron; keff is the effective multiplication factor which is assumed to be unity. In the design of initial fuel LP, keff value should be larger than 1 to guarantee adequate excess reactivity. The infinite multiplication factor of the whole reactor core k1 and the fast fission factor e can be written as:
k1 ¼
Pðr 0 ÞI0 ðr 0 Þr 0 dr
0
where R0 is the equivalent radius of the reactor core.
Fig. 1. Typical fuel assembly arrangement of VVER-1000.
D2 r2 /2 þ R2 /2 ¼ R1;2 /1
r0
where P(r) is one-dimensional radial power distribution; K0 is the zeroth order modified Bessel function of the first kind; I0 is the zeroth order modified Bessel function of the second kind; b1 is the albedo of fast neutron; r0 is
36
24B36
16
Z
R r0
/0 ¼
36
24B36 16
24B36
16 24
16
I0 ðrÞK 0 ðr 0 Þ I0 ðr 0 Þ
36
24
16
A3 ðrÞ ¼
ð6Þ
1 þ eL22 k1
As shown in Eq. (8), the fuel potential K⁄ is a function of the power distribution. Integrating K⁄ over the space and divide it by the area of core, we get the average total fuel potential, denoted as Fp
R r0 Fp ¼
0
2prK ðrÞdr pr20
ð14Þ
Since Fp is the average total potential relevant to keff, the LP optimization can be obtained by seeking proper P(r) to minimize Fp under the power peaking constraint condition 0 6 PðrÞ 6 Plim : For simplicity, the power distribution P(r) is set as a third order polynomial and normalized by Eq. (16)
PðrÞ ¼ a0 þ a1 r þ a2 r 2 þ a3 r 3 R r0 2prpðrÞdr a0 a2 r 20 a3 r 30 0 ¼ 2 þ þ ¼1 2 4 5 pr20
ð15Þ ð16Þ
where a0, a1, a2, a3 are coefficients which indicate the power distribution. In the center of core, the neutron is total reflection, so
P0ðrÞ jr¼0 ¼ a1 ¼ 0
ð17Þ
P(r)should not be negative in the overall core, and suppose P(r) get maximum at rm = m r0,
P0 ðrÞjr¼rm ¼ a1 þ 2a2 þ 3rm a3 ¼ 0
ð18Þ
a0 þ r 2m a1 þ r 2m a2 þ r 3m a3 ¼ Plim
ð19Þ
3.3. Calculation and analysis
1
where L2 is the thermal diffusion length. In cylindrical coordinates, Eq. (5) could be written as 2
d /1 ðrÞ dr
2
þ
1 d / ðrÞ /1 ðrÞ ¼ K ðrÞ/1 ðrÞ r dr 1
ð7Þ
The admissible value of Plim is around 1.3 (Boroushaki et al., 2003; Mazrou and Hamadouche, 2006). Now P(r) could be obtained through calculating Eq. (20) which is combined by Eqs. (16), (18), and (19).
CA¼J
ð20Þ
128
C. Dai et al. / Annals of Nuclear Energy 50 (2012) 126–132
8 _ wðtÞ ¼ @H > @X > < _ XðtÞ ¼ @H @w > > _ _ : _ _ _ HðX ; U ; w; tÞ P HðX ; U; w; tÞ
ð25Þ
where ^ denotes optimum quantities in admissible region. In the g l ðXÞ 6 0 Eq. (23) is replaced case X(t) and U(t) are all constrained qj ðUÞ 6 0 by more general condition (Zhao and Ren, 2009):
@H ¼0 @U
ð26Þ
with
H ¼ H þ
n s X X ki ðtÞpi ðX; UÞ þ lj ðtÞqj ðUÞ i¼1
where ki ðtÞ and lj(t) are Lagrange factors, pi(X, U) is the lowest order derivative of gi(X) with respect to t in the derivatives of containing explicit functions U(t). If gi(X) < 0, ki ðtÞ ¼ 0 satisfy, optimum trajectories do not be in the constrained boundary gi(X) = 0; if qj(U) < 0, lj(t) satisfy, admissible control U(t) do not lie in the con strained boundary qj(U) = 0; if qj(U) = 0 and lj(t) – 0, admissible control U(t) lie in qj(U) = 0.
Fig. 2. The relationship between Fp and m.
where J ¼ ½1 P lim 0T , A ¼ ½a0 a2 a3 T
2
1 r 20 =4 r 30 =5
6 C ¼ 41
r 2m
0
2
3
7 r 2m 5 3r m
4.2. Modeling and calculation
m is the most important parameter as it determines the peaking point of the power distribution, and A is computed by giving m in its allowable domain (0 6 m 6 1). The relationship of Fp and m is shown in Fig. 2. Fp reaches to the minimum value when the value of m is between 0.5 and 0.6, neither m = 0 nor m = 1.0. Some values of Fp are calculated with different m, listed in Table 2. 4. Optimization of BP Burnable poisons should be arranged in the core to compensate the excess reactivity and flatten the power distribution shape on the basis of the optimum fuel LP obtained in Section 3. Using the single group diffusion theory, the optimization of BP is carried out under the constraint condition of the maximum power density. 4.1. Pontryagin’s maximum principle (Pontryagin et al., 1962) Let the plant
X_ ¼ f ðX; U; tÞ
ð21Þ
Single group diffusion equation is
r Dr/ þ Ra / þ Rc / ¼
Z
mRf /
ð28Þ
keff
where Ra is the absorption cross section, vRf is the fission cross section, keff is the effective multiplication factor and D is the diffusion coefficient. A variable denoted as Rc is introduced to represent the burnable poisons. Assuming that there is no leakage in the axial direction of core, we can get
D 0 mRf / D/00 þ D0 þ Ra / Rc / ¼ 0 / þ R keff
ð29Þ
Let x1(R) = u(R), x2(R) = u0 (R), x3(R) = R, Eq. (28) could be rewritten as
8 dx1 ðRÞ > ¼ x2 ðRÞ > > < dR 0 dx2 ðRÞ ðRÞ ¼ DDðRÞ þ x31ðRÞ x2 ðRÞ þ BðRÞx1 ðRÞuðRÞ BðRÞx1 ðRÞ dR > > > dx3 ðRÞ : ¼1 dR
Have an associated cost index of
Max Jðt 0 Þ ¼ vðXðT 0 Þ; T 0 Þ þ
ð27Þ
j¼1
ð30Þ
mRf ðRÞ
T0
LðX; U; tÞdt
ð22Þ
t0
where X ¼ ½x1 x2 . . . xn1 T , U ¼ ½u1 u2 un2 T . If the control is unconstrained the condition for optimality is
@H ¼0 @U
ð23Þ
where BðRÞ ¼
keff
Ra ðRÞ
DðRÞ
Rc ðRÞ , u(R) is control vector uðRÞ ¼ mRf ðRÞ keff
the neutron economy point of view, low neutron leakage is expected, which requires low the edge neutron flux u0, so the boundary condition is:
x1 ðR0 Þ ¼ /0
ð31Þ
x2 ð0Þ ¼ 0
where H is Hamiltonian function
. From
Ra ðRÞ
The power density is
H ¼ LðX; U; tÞ þ wT f ðX; U; tÞ
ð24Þ
where w is Lagrange factor, when the control is constrained, Eq. (23) changes into
Pv ðRÞ ¼ F a Ef Rf ðRÞ/ðRÞ
ð32Þ
where Fa is the ratio of the heat of core (mainly fuel element and coolant) to the total heat of the reactor, Ef is the energy released
Table 2 Some typical values of m and Fp. m Fp
0.0 1.088445
0.1 1.085443
0.3 1.073669
0.5 1.048841
0.7 1.068476
0.9 1.093214
1.0 1.107083
129
C. Dai et al. / Annals of Nuclear Energy 50 (2012) 126–132
per fission, and it equal to about 200 MeV. The constraint condition is a. the cross section of control material should not be negative, namely the control domain
qðuÞ ¼ uðRÞ 6 0
ð33Þ
b. the power density should not larger than the maximum power density Pvmax, so that
gðXÞ ¼ F a Ef Rf x1 Pv max 6 0
ð34Þ
The radial PPF which is expected as low as possible could be obtained through dividing the Pvmax by the average power density of the reactor core. As the total area of the core is constant, the minimizing PPF process is then translated to searching the optimum admissible control u(R) to obtain the maximum output power Pt computed by integrating Eq. (32) over the whole core, so
Max Pt ¼ 2pF a Ef ¼ 2pF a Ef
Z Z
R0
from which we can get Rc = 0, that means the BP is needless in this region. ii. l (R) = 0, with the same foundation as case ‘i’, we can get
@H ¼ w2 Bx1 ¼ 0 @u
ð44Þ
Because neither B(R) nor x1(R) is 0, the last variable w2(R) must equal 0. Substituting w2(R) = 0 into Eq. (37), the adjoint equation can be written as
8 dw 1 > ¼ Rf x3 > < dR dw2 ¼ w1 dR > > : dw3 ¼ Rf x1 dR
ð45Þ
As it is impossible for practical engineering situation that dw1 ¼ Rf x3 ¼ 0, case ii does not exist. When g(X) = 0, k – 0, generaldR ized Hamilton is
Rf ðRÞ /ðRÞ R dR 0
H ¼ H þ kðRÞg 00 ðXÞ
R0
Rf ðRÞx1 ðRÞx3 ðRÞdR
ð35Þ
0
4.3. Calculation and analysis
ð46Þ
the optimum condition is 8 < @H ¼ w2 Bx1 þ kRf Bx1 ¼ 0 @u n h io 00 : g ðXÞ ¼ F a Ef 2x2 R0f þ x1 R00f þ Rf 1 þ D0 x2 þ Bx1 uðRÞ Bx1 ¼ 0 x3 D
ð47Þ
Hamilton function is
0 D 1 x2 þ Bx1 uðRÞ B2 x1 þ w3 H ¼ Rf x1 x3 þ w1 x2 þ w2 þ D x3
after calculation, we can obtain
ð36Þ
f 1
dR
f
ð37Þ
x23
The derivatives of g(X) are
g 0 ðXÞ ¼
dgðXÞ ¼ F a Ef ðRf x2 þ x1 R0f Þ dR
ð38Þ
keff
d gðXÞ 2
dR
1 D0 ¼ F a Ef 2x2 R0f þ x1 R00f þ Rf þ x2 þ Bx1 uðRÞ Bx1 x3 D ð39Þ It can be found that g0 (X) does not contains the admissible control vector u(R) explicitly while g 00 ðXÞ does. The augmented Hamilton function is
H ¼ H þ kðRÞg 00 ðXÞ þ lðRÞqðuÞ
ð40Þ
Analysis on the optimum solutions in different conditions is then carried out, with g(X) < 0, kðRÞ ¼ 0, Eq. (40) changes into
H ¼ H þ lðRÞqðuÞ
ð41Þ
Here, two cases are possible for optimum value of u(R), ( i. l (R) – 0, there are the relations in accordance with Eq. (27) @H @u
¼ @H l¼0 @u
qðuÞ ¼ 0
ð42Þ
The compute result is
lðRÞ ¼ w2 Bx1 uðRÞ ¼ 0
02
00
f
f
0
ð48Þ
0
ð43Þ
f
when g(X) = 0, the optimum trajectories of power are in the constraints boundary and the neutron and BP distribution could be obtained from Eq. (48); while if g(X) < 0, the BP is needless and Eq. (33) could provides the optimum power distribution. As keff is unknown, a method so-called ‘‘source iteration’’ is used as follows. The set of ð0Þ initial value keff and (mRfu)(0) makes the fission source mRfu/keff in Eq. (28) known, which could give the solution of (mRfu)(1). Then according to the definition of keff ðnþ1Þ
2
g 00 ðXÞ ¼
0
R 2R R DR > : uðRÞ ¼ 1 BRf R þ BRf2 BRf DBRf
The adjoint equation is
8 dw1 > ¼ Rf x3 w2 ½B þ uðRÞB > dR > < dw 0 2 ¼ w1 þ w2 ðx13 þ DD Þ dR > > > dw3 ¼ R x w2 x2 :
8 P v max > < x1 ðRÞ ¼ F a Ef Rf ðRÞ
R R0 ¼
0 1 ðnÞ keff
mRf ðRÞx1ðnþ1Þ ðRÞ2pRdR mRf ðRÞxðnÞ 1 ðRÞ2pRdR 0
ð49Þ
R R0
ð1Þ
we can update the effective multiplication keff ; so that X can be calculated again. Repeat this process until the convergence criterion satisfy: (a) characteristic value convergence criterion j
e1 ; and (b) neutron flux convergence criterion j
ðnÞ ðn1Þ k eff eff ðnÞ k eff
k
ðnÞ ðn1Þ x1 ðRÞx1 ðRÞ ðnÞ x1 ðRÞ
j < e2 ,
–4
e1 and e2 are small positive value, for example e1 = 10 e2 = 10–3.
j<
and
5. Multi objective optimization of LP and BP The optimum fuel LP and BP are obtained by minimizing Fp and PPF (namely maximizing Pt) respectively in Sections 3 and 4. However LP and BP are coupled in the power distribution optimum calculation, the two objective functions cannot be optimum simultaneously. With Fp reaches to the minimum value, PPF may not be the minimum as the neutron flux at the periphery of core is lower than its limit. On the other hand, the optimum LP required by the objective PPF cannot always satisfy the minimizing of Fp. As a result, a multi objective model is applied to get an eclectic solution in this study.
130
C. Dai et al. / Annals of Nuclear Energy 50 (2012) 126–132
5.1. Multi objective optimization A multi-objective problem (MOP) consists of optimizing (i.e., maximizing or minimizing) several objectives simultaneously, with a number of constraints. MOP can be described as follows (Sanaye and Hajabdollahi, 2010):
Max F i ¼ ðXÞ;
i ¼ 1; 2; . . . ; I
ð50Þ
where X is decision vector
X ¼ ðx1 ; x2 ; . . . ; xn Þ
ð51Þ
The constraint conditions are
g l ðXÞ 6 0 l ¼ 1; 2; . . . ; L
ð52Þ
hj ðXÞ ¼ 0 j ¼ 1; 2; . . . ; s
ð53Þ
Usually, it is impossible to maximize or minimize all the objective function at the same time because of the interaction between them. Therefore the results of optimum designs called ‘‘Pareto optimum solutions’’ are a set of multiple optimum solutions which are compromised by giving different focus to every objective. The genetic algorithm (GA) is a semi-stochastic method based on Darwin’s laws of natural selection. To solve MOP, a non-dominated sorting genetic algorithm which is used in this research was proposed by Srinivas and Deb (1994) (NSGA) and modified by Deb et al. (2002) (NSGA-II).
Fig. 3. The distribution of Pareto optimum points solutions.
Table 3 The optimum values of objective functions for points A–E in Pareto-optimal fronts.
m Fp PPF
O
A
B
C
D
E
0.0 1.088445 1.300000
0.4859 1.050267 1.195738
0.5082 1.047541 1.196277
0.5394 1.044717 1.196919
0.5561 1.043975 1.200540
0.5633 1.043893 1.204250
5.2. Objective functions and constraints Fp and PPF are denoted as the two objectives respectively in the MOP, and as m has an effect on Fp which acts as the basis of the optimization of PPF, it is regarded as the decision variable, subject to
06m61
ð54Þ
Now the MOP is described as
Min s:t:
F 1 ðmÞ ¼ F p F 2 ðmÞ ¼ PPF
m60
ð55Þ
m160 5.3. Multi objective optimum results Using a search population size of M = 200 individuals, crossover probability of pc = 0.9, gene mutation probability of pm = 0.035, the NSGA-II method is performed for 100 generations, the results for Pareto-optimum curve are shown in Fig. 3, which denote the conflict between Fp and PPF clearly. Any change that increases Fp which indicates the fuel consumption leads to a decrease of PPF, and vice versa. Five typical points (A–E), for which the optimum values of Fp and PPF are listed in Table 3, are chosen as object points to study the relationship between optimum Fp and PPF. The minimum Fp exists at design point E (1.043893), where PPF is the highest (1.204250). On the other hand the maximum Fp occurs at point A (1.050267), with a lowest PPF (1.195738). Design points A and E are the optimum situations at which PPF and Fp are single objective function respectively. The distribution of m for Pareto front is shown in Fig. 4, the lower and upper bounds of m are 0 and 1 respectively. m has its optimum values in the range of 0.4859 < m < 0.5633, Fp and PPF in this margin are all lower, compared with point O where m = 0. It means the solutions are optimized successfully by NSGA-II method. The selection of final optimization among the optimum points existing on the Pareto front needs a process of decision-making
Fig. 4. Scattering of m for the Pareto optimum front.
based on the experiences of engineer. In the process, the focus should be given to the importance of each objective function. If Fp is considered to be more important, the points from D to E could be picked up as the final solutions, while points during A and B should be chosen when paying more attentions to PPF. In this study, the importance of the two objective functions is thought to be equivalent, so the points around C are regarded as the final optimum solution. Take typical point C where m = 0.5394 as the discussion object, the non-uniform distribution of the fuel potential K⁄ is shown in Fig. 5. K⁄ firstly increases until P equals Plim, then K⁄ will decrease with the increasing of radial relative position (RRP) in the interval from 0.81 to 1. The minimum K⁄ appears at the edge of core. The relationship between the fuel enrichment and K⁄, which could be utilized to provide distribution of fuel enrichment against the ra-
C. Dai et al. / Annals of Nuclear Energy 50 (2012) 126–132
Fig. 5. The optimum distributions K⁄.
Fig. 7. Distribution of control vector (u(R)).
Fig. 6. The optimum enrichment (U235) distribution.
Fig. 8. Distribution of the optimum power density.
dial relative position, is calculated by Eq. (6). Then, the theoretical relationship between U-235 enrichment and percentage radial distance is shown in Fig. 6. Fig. 6 shows the non-uniform distribution of the fuel optimum enrichment. The enrichment increases when RRP < 0.81, after that it decreases until the edge of core where the enrichment is the lowest. From Figs. 5 and 6, it can be concluded that this LP scheme could not only guarantee that power in the center of core relatively low, but also could make sure the edge neutron flux is not too high as leading to a high neutron leakage. It should be emphasized that the optimal LP obtained here is known as the ‘low leakage loading pattern’. Compared with value which the power peaking is at the center of core (m = 0), Fp and PPF of point C are lower, that means the optimum LP and BP saved the fuel consumption about 4.2% and reduced the radial power peaking factor about 7.93%. As shown in Fig. 7, when RRP < 0.38 or RRP > 0.74, the control vector u(R) = 0 and as a result, the BP is needless; while 0.38 < RRP < 0.74, BP must be arranged to flat the power density peaking caused by the high enrichment fuel. As shown in Fig. 8, power density increases with the increasing of RRP because the enrichment and fission cross section of the fuel rises. When RRP = 0.38, the Pv reaches the maximum value and keeps the maximum value until RRP = 0.74 as the arrangement of BP. When
131
Fig. 9. Distribution of the optimum neutron flux.
0.74 < RRP < 1, Pv decreases and reaches minimum at the edge of core. The normalized neutron distribution in Fig. 9 shows the similar trend as the power distribution in Fig. 6. However, when
132
C. Dai et al. / Annals of Nuclear Energy 50 (2012) 126–132
Table 4 Comparing Fp and PPF with different design for VVER-1000.
VVER-1000 designer Perturbation theory Point C
Fp
PPF
1.04385 1.04391 1.044717
1.22 1.215 1.196919
0.38 < RRP < 0.74, the neutron flux decreases as the optimum trajectory lie in the constraint boundary and the power density is constant while the enrichment increases tardily. There is a relation / (R) = Pvmax/(FaEfRf(R)) in this region, so, obviously, / will decrease with RRP. 5.4. Results and discussion An optimization process for the first operating cycle of a typical VVER-1000 reactor core is performed. The comparison of final optimum results of point C with the references (FSAR, 2003; Hosseini and Vosoughi, 2012) reported load pattern are shown in Table 4. Although the Fp is larger than the two references (about 0.083% and 0.077%), the PPF is reduced by 1.89% and 1.49%, compared with the original VVER-1000 design and perturbation theory optimization respectively. This is result from the coordinating the interaction between Fp and PPF during chosen process of the final optimum design. 6. Conclusions A two-level optimization method is applied to design the initial optimum LP and BP arrangement of VVER-1000. In the lower level optimization, Fp is firstly minimized by searching the optimum LP with the backward diffusion calculation theory, the result shows that when m is in the margin of 0.5 < m < 0.6, Fp is the lowest. On the basis of fuel distribution, BP is optimized to get the minimum radial power peaking factor. Then, considering the interrelationship of each other, an in-depth study which acts as the upper level optimization is carried out to obtain the optimum LP and BP with NSGA-II. A set of Pareto optimum points in the optimum region (0.4859 < m < 0.5633) which clearly revealed the conflict between the objectives is obtained and showed. After the analysis of the relationship of Pareto solutions and weighting of LP and BP, the points around point C are chosen and deemed to be the final optimum solutions. Point C (m = 0.5394) is taken as the typical result for discussion. The absorbing material should be arranged from 0.38 to 0.74 time of core equivalent radius. The calculation shows that not only the fuel consumption is saved about 4.2%, but also the PPF is reduced by 7.93%, compared with a uniform loading pattern core. In the comparison to the
VVER-1000 original design and some other optimization, because of the coordinating the interaction between Fp and PPF during deciding process of the final optimum design, PPF is reduced by 1.89% and 1.49% respectively, however the Fp is larger (about 0.083% and 0.077%). The application of these methods and the optimum solutions will provide theoretical references for improving the initial designs of LP and BP in reactor core. Acknowledgment This research is funded by National Natural Science Foundation of China (10975112). References Boroushaki, M., Ghofrani, M.B., Lucas, C., 2003. Optimum fuel core loading pattern design in PWR nuclear power reactors using genetic algorithm and fuzzy nonlinear programming. J. Intell. Fuzzy Syst. 14, 85–93. Chao, Y.A., Suo, C.A., Penkrot, J.A., 1983. On the theory of interface flux nodal method. In: Proc, CNS/ANS Int. Conf. Numerical Methods in Nuclear Engineering, Montreal, p. 307. Chao, Y.A., Hu, C.W., Suo, C.A., 1986. A theory of fuel management via backward diffusion calculation. Nucl. Sci. Eng. 93, 78–87. Chapot, J.L.C., Silva, F.C., Schirru, R., 1999. A new approach to the use of genetic algorithms to solve the pressurized water reactor’s fuel management optimization problem. Ann. Nucl. Energy 26, 641–655. Deb, K., Pratap, A., Agarwal, S., Meyarivan, T., 2002. A fast and elitist multi objective genetic algorithm: NSGA-II. IEEE Trans. Evol. Comput. 6, 182–197. Final Safety Analysis Report of BNPP-1, 2003. Reactor, 49.BU.1 0. 0. OO.FSAR. RDR001, Rev. 0 (Chapter 4). Haling, R.K., 1964. Operating Strategy for Maintaining Optimum Power Distribution throughout Life. TID 7672. Hosseini, M., Vosoughi, N., 2012. Development of a VVER-1000 core loading pattern optimization program based on perturbation theory. Ann. Nucl. Eng. 39, 35–41. Kim, Y.J., Downar, T.J., Sesonske, A., 1987. Optimization of core reload design for low-leakage fuel management in pressurized water reactor. Nucl. Sci. Eng. 96, 85–101. Kropaczek, D.J., Turinsky, P.J., 1991. In-core nuclear fuel management optimization for pressurized water reactors utilizing simulated annealing. Nucl. Technol. 95, 9–31. Mazrou, H., Hamadouche, M., 2006. Development of a supporting tool for optimum fuel management in research reactors using artificial neural networks. Nucl. Eng. Des. 236, 255–266. Pontryagin, L.S., Boltyanskii, V.G., Gamkrelidze, R.V., Mischenko, E.F., 1962. The Mathematical Theory of Optimum Process. Wiley-Interscience Publishers, New York. Suh, J.S., Levine, S.H., 1990. Optimized automatic reload program for pressurized water reactors using simple direct optimization techniques. Nucl. Sci. Eng. 105, 371–382. Srinivas, N., Deb, K., 1994. Multi objective optimization using non-dominated sorting in genetic algorithms. J. Evol. Comput. 2, 221–248. Sanaye, S., Hajabdollahi, H., 2010. Multi-objective optimization of shell and tube heat exchangers. Appl. Therm. Eng. 30, 1937–1945. Tran, H.N., Yammoto, A., Yamane, Y., 2009. An improved inverse analysis model for fuel loading pattern optimization. J. Nucl. Sci. Technol. 46, 1162–1169. Zhao, F.Y., Ren, X.R., 2009. Optimum performance analysis of a sodium cold trap. Nucl. Eng. Des. 239, 490–494.