- Email: [email protected]
A nodal modal method for the neutron diffusion equation. Application to BWR instabilities analysis
Annals of Nuclear Energy 29 (2002) 1171–1194 www.elsevier.com/locate/anucene
A nodal modal method for the neutron diffusion equation. Application to BWR instabilities analysis R. Miro´a, D. Ginestarb, G. Verdu´a,*, D. Hennigc a
Departamento de Ingenierı´a Quı´mica y Nuclear, Universidad Polite´cnica de Valencia, Camı´ de Vera, 14, 46022, Valencia, Spain b Departamento de Matema´tica Aplicada, Universidad Polite´cnica de Valencia, Camı´ de Vera, 14, 46022, Valencia, Spain c Laboratory for Reactor Physics and Systems Behaviour (LRS), Paul Scherrer Institut (PSI), CH-5232 Villigen PSI, Switzerland Received 26 July 2001; received in revised form 22 September 2001; accepted 24 September 2001
Abstract Fast codes, capable of dealing with three-dimensional geometries, are needed to be able to simulate spatially complicated transients in a nuclear power reactor. In this paper, we propose a modal method to integrate the neutron diffusion equation in which the spatial part has been previously dicretized using a nodal collocation method. For the time integration of the resulting system of differential equations it is supposed that the solution can be expanded as a linear combination of the dominant Lambda modes associated with a static configuration of the reactor core and, using the eigenfunctions of the adjoint problem, a system of differential equations of lower dimension is obtained. This system is integrated using a variable time step implicit method. Furthermore, for realistic transients, it would be necessary to calculate a large amount of modes. To avoid this, the modal method has been implemented making use of an updating process for the modes at each certain time step. Five transients have been studied: a homogeneous reactor, a non-homogeneous reactor, the 3D Langenbuch reactor and two transients related with in-phase and out-of-phase oscillations of Leibstadt NPP. The obtained results have been compared with the ones provided by a method based on a one-step backward discretization formula. # 2002 Elsevier Science Ltd. All rights reserved.
* Corresponding author. Tel.: +34-9-6387-7635; fax: +34-9-6387-7639. E-mail address: [email protected] (G. Verdu´). 0306-4549/02/$ - see front matter # 2002 Elsevier Science Ltd. All rights reserved. PII: S0306-4549(01)00103-7
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
1172
1. Introduction To simulate the behaviour of a nuclear power reactor it is necessary to be able to integrate the time-dependent neutron diffusion equation inside the reactor core. In this way, it is interesting to develop fast and reliable methods to integrate this equation. In this work, we will develop a modal method based on the assumption that the solution for the neutron flux of the studied transient can be expressed approximately as a linear combination of the dominant eigenfunctions of a static auxiliary problem, known as the Lambda modes, associated with a static reactor configuration (Henry, 1982). To calculate these auxiliary eigenfunctions we have used a nodal serendipity collocation method (He´bert, 1987) to discretize the equations, and an iterative method to obtain the dominant eigenvalues and their corresponding eigenfunctions (Verdu´ et al., 1994, 1999). This modal method is interesting to study transients where the number of necessary eigenmodes is not very large, and it can also be useful to analyze transients with a physical interpretation for the different amplitudes of the neutron flux expansion, as the ones appearing in the study of in-phase and out-of-phase oscillations of boiling water reactors (Miro´ et al., 2000). In the case of one eigenmode, the method could be considered as a version of the quasi-static method to integrate the neutron diffusion equation (Ott, 1969).
2. Nodal modal method Under general assumptions (Stacey, 1969), the neutronic population inside a nuclear power reactor can be modeled by the neutron diffusion equation in the approximation of two energy groups. This model is of the form,
K X : 1 þ L ¼ ð1 ÞM þ lk Ck ;
ð1Þ
k¼1
: C k ¼ k f 1 f 2 lk Ck ; k ¼ 1; . . . K;
ð2Þ
where, K is the number of delayed neutron precursors groups considered, 2 L¼4
~ þ a1 þ 12 ~ D1 r r 12
0
5; ~ þ a2 ~ D2 r r
and M¼
f 1 0
3
f 2 1 1 ; ¼ ; ¼ : 0 2 0
1
21 6 ¼ 4 1 0
0
3
7 1 5; 2
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
1173
The diffusion constants and cross-sections, Dg, 12, ag, fg, g=1,2, appearing in the equations depend on the reactor materials, that is, they are position dependent functions. To solve the problem, we begin with an approximation, consisting of discretizing the reactor in N nodes, where the nuclear properties are supposed to be constant. Next, a nodal collocation method (He´bert, 1987; Verdu´ et al., 1994), that assumes a truncated expansion of the neutron flux in terms of Legendre polynomials is used. This method permits to discretize the spatial part of the equations, approximating Eqs. (1) and (2) by a set of ordinary differential equations of the form
1
:
þ L ¼ ð1 ÞM þ
: XC k ¼ k M lk XCk ;
K X lk XCk ;
ð3Þ
k¼1
ð4Þ
where and Ck are vectors, whose components are the corresponding coefficients of the Legendre polynomials expansions of and Ck in each node, at each time step. L and M have the following block structure,
0 L11 M11 M12 ; L¼ ; M¼ L21 L22 0 0 and we have introduced the matrix X, I X¼ : 0 To integrate Eqs. (3) and (4) using the modal approximation, we assume that can be expressed approximately as Md X nl ðtÞ l ;
ð5Þ
l¼1
where l, l=1, . . ., Md are the eigenvectors associated with the Md dominant eigenvalues of a static problem of the form (Verdu´ et al., 1994) L0
l
¼
1 M0 l ; kl
ð6Þ
where the matrices L0 and M0 are the matrices L and M of Eqs. (3) and (4), for a given configuration of the reactor core. This problem is known as the Lambda modes problem (Henry, 1982; Verdu´ et al., 1994). Since L0 is not a hermitic matrix, we introduce the adjoint problem (Henry, 1982), Ly0
y l
¼
1 y M kl 0
y l;
ð7Þ
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
1174
y m,
to obtain a set of modes, n, and their adjoints, relation, y y m ; M0 n ¼ m ; M0 m n;m ¼ Nm n;m ;
satisfying the biorthogonality
ð8Þ
where n,m is the Kronecker delta function. Multiplying Eqs. (3) and (4) by ym , writing L ¼ L0 þ L;
M ¼ M0 þ M;
and making use of expansion (5), and Eq. (6), we obtain the following set of equations Md X
y m;
1
l¼1
x
Md X
Md X d 1 nl þ l dt k l¼1 l
y m ; M0
l
nl þ ð 1 Þ
l¼1
d dt Md X
y m ; M0
Md X
l
nl þ
y m ; M0
y m ; L l
nl ¼ ð1 Þx
l¼1 y m ; M
l
nl þ
l¼1 y m ; XCk
Md X
K X lk
y m ; XCk
;
k¼1
¼ k l
l¼1
nl þ k
Md X
y m ; M
l
nl lk
y m ; XCk
:
ð9Þ
l¼1
Using the biorthogonality relationship (8) and introducing the notation,
ml ¼ ym ; 1 l ; ALml ¼ ym ; L l ; y Cmk ¼ ym ; XCk ; AM ml ¼ m ; M l ; and the mode m reactivity, defined as m ¼
km 1 ; km
Eqs. (9) can be expressed as, Md Md Md K X X X X d
ml nl ¼ ð m ÞNm nm þ ð1 Þ AM n ALml nl þ lk Cmk ; ml l dt l¼1 l¼1 l¼1 k¼1 Md X d Cmk ¼ k Nm nm þ k AM ml nl lk Cmk ; dt l¼0
k ¼ 1; . . . ; K:
ð10Þ
Varying the index m from 1 to Md, we can rewrite these equations in the following matrix form
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
1175
d½n ¼ ½ 1 ½ I½N½n þ ð1 Þ½ 1 AM ½n ½ 1 AL ½n dt K X þ lk ½ 1 ½Ck ;
ð11Þ
k¼1
d½Ck ¼ k ½N½n þ k AM ½n lk ½Ck ; k ¼ 1; . . . ; K: dt Considering explicitly the K equations for the delayed neutron precursors, Eqs. (11) can be grouped as the following system of differential equations d N ¼ AN ; dt where 2 3 n 6 C1 7 7 N ¼6 4 .. 5; . CK and
2 6 6 A¼6 4
ð12Þ
½ 1 ½ I½N þ ð1 ÞAM AL 1 ½ N þ A M .. . M K N þ A
½ 1 l1 ½ 1 lK l1I
0 . .. . . . ... 0 lK I
3 7 7 7: 5
If the number of dominant modes needed in expansion (5), Md, is not very large, the dimension of the new system of differential Eqs. (12) is quite smaller than the one of the initial system (3) and (4). To solve the system (12) we have used a variable step implicit method (Hindmarsh, 1983) because of the stiffness of the differential equations. To obtain the initial conditions for the time integration, we start from a critical configuration of the reactor. To obtain this critical reactor configuration we solve the Lambda modes problem for a given initial configuration, L0
1
¼
1 M0 k1
1;
for the fundamental mode. Dividing the fission cross-sections by the fundamental eigenvalue, k1, we obtain L0
1
¼ MCrit
1;
ð13Þ
where MCrit is a matrix which components are the components of M divided by k1. In this way we obtain that 1 is a stationary solution of Eqs. (3) and (4) for a critical configuration of the reactor. The subsequent fission cross-sections of Eqs. (12) will be also divided by the eigenvalue k1 to calculate the time evolution of N .
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
1176
Eq. (13), together with the equation 0 ¼ k MCrit
1
lk XCk ;
ð14Þ
constitute the set of equations defining the steady state associated with the critical reactor configuration. Hence, the components of vector [n] for the initial condition are n1 ¼ 1;
nl ¼ 0;
l ¼ 2; . . . ; Md :
Multiplying (14) by y1 we have D E D E 0 ¼ k y1 ; MCrit 1 lk y1 ; XCk ; Thus, the initial values for [Ck] are given by E k D y Crit C1;k ¼ 1 ; 1; M lk Cl;k ¼ 0; l ¼ 2; . . . ; Md : 2.1. Modes updating For realistic transients, the nuclear cross-sections are time dependent functions and to use the modal method proposed above, it would be necessary to calculate a large amount of modes. This is prohibitive from the computational point of view. Thus, instead of calculating a large amount of modes, we calculate only a small number of them and we update these modes each certain time step, ti. In this way, to integrate the neutron diffusion equation in the time interval [ti, ti+1], we use the Lambda modes associated with the problem 1 Li il ¼ i Mi il ; kl where Li and Mi are the matrices associated with the reactor configuration at time ti. Hence, the differential equations to integrate are the following ones, Md X
iy m;
1
l¼1
x
Md X
Md i d i X 1 nl þ l dt ki l¼1 l
iy i m; M
i l
nil þ ð1 Þ
l¼1
d dt Md X l¼1
Md X
iy i m; M
i l
nil þ
iy i m; M
i l
iy i m ; L
i l
nil ¼ ð1 Þx
l¼1 iy i m ; M
i l
nil þ
l¼1 iy m ; XCk
Md X
K X lk
iy m ; XCk
;
k¼1
¼ k Md X nl þ k
iy i m ; M
i l
nil lk
iy m ; XCk
:
l¼1
the operators Li and Mi change along the time interval.
ð15Þ
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
1177
To calculate the initial conditions for nim ðti Þ; we reconstruct the vector ðti Þ ¼
Md X ni1 l ðti Þ
i1 l ;
l¼1
from the variables ni1 l ðti Þ and nim ðti Þ as 1
nim ðti Þ ¼ D
iy i m; M
i m
E
iy i m; M
i1 l
obtained in the last time step, and we calculate
ðti Þ :
ð16Þ
To obtain the variables related to the concentration of the delayed neutron pre i1y cursors, iy ; XC ð Þ; in terms of ; XC ð Þ; which are known quantities t t k i k i m m from the integration in the previous time interval, we suppose that approximately we can write iy m
Md X aml
i1y ; l
l¼1
where aml ¼ D
iy i1 i1 m; M l E; i1y i1 i1 ; M l l
and this yields
iy m ; XCk
ðti Þ
Md D X aml
i1y ; XCk l
E
ðti Þ;
ð17Þ
l¼1
which completes the initial conditions to integrate (15) in [ti, ti+1].
3. Numerical results The exposed nodal modal method has been implemented in a code written in fortran and called modkin, which calculates the dominant Lambda modes associated with a given configuration of the reactor core, making use of the implicit restarted Arnoldi method (Verdu´ et al., 1999). To test the performance of the method, we have studied five transients: one associated with a homogeneous reactor without delayed neutron precursors, another one associated with a non-homogeneous reactor with one group of delayed neutron precursors, the Langenbuch problem (Langenbuch et al., 1977) and two transients corresponding to in-phase and out-of-phase oscillations, whose nuclear data where obtained using a ramona thermal-hydraulic model of Leibstadt NPP. The former three transients have been used as validation of the code. For the first transient we have an exact analytical solution (Miro´, 2001).
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
1178
For the other examples, the results obtained with the nodal modal method have been compared with the ones obtained with a code that integrates Eqs. (3) and (4) using a temporal discretization based on the one-step backward difference formula, called nokin3d (Verdu´ et al., 1995; Ginestar et al., 1997) which is an accurate method to integrate the neutron diffusion equation. In all the numerical examples studied, we have only considered changes in the cross-sections. This kind of transients can be analyzed using only spatial modes. More general perturbations should be addressed with more general models as the multigroup (more than two) neutron diffusion equation or the neutron transport equation. 3.1. Homogeneous reactor This problem consists of a 10060180 cm3 homogeneous reactor, divided in 36 equal nodes (336), as it is shown in Fig. 1. These nodes correspond to only one material, which nuclear properties are given in Table 1. The transient consists of a step perturbation on the thermal fission cross-section f 2 ðtÞ ¼ 0f2 þ f 2 ¼ 0f 2 þ 1 104 HðtÞ; where H(t) is the unit step function, and 0f2 is the initial thermal fission cross-section, which is chosen imposing that the initial configuration of the reactor is critical.
Fig. 1. Homogeneous reactor. Table 1 Group constants for the homogeneous reactor Group
Dg (cm)
ag (cm1)
fg (cm1)
12 (cm1)
1/g (s/cm)
1 2
1.28205 0.666667
0.01 0.1
0.01 0.1097683675455275
0.075 –
1107 1105
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
1179
That is, the initial configuration of the reactor is a steady state configuration described by the first eigenvector of a Lambda modes problem associated with the reactor such that the first eigenvalue of this problem is equal to 1 (Henry, 1982). This problem will be used as a first test for the nodal modal method. The transient has been followed during 1 s using an integration time step of 1103 s, and the results obtained with the nodal modal method using one mode have been compared with the analytical solution for the relative neutronic power evolution. The results are presented in Fig. 2, showing a very good agreement.
Fig. 2. Power evolution for the homogeneous reactor.
Table 2 Nuclear properties of the non homogeneous reactora Mat.
Group
Dg (cm)
ag (cm1)
fg (cm1)
12 (cm1)
1/g (s/cm)
1
1 2
1.695310 0.409718
0.0139530 0.2614097
0.013406855 0.342323710
0.0164444 –
0 4.5455106
2
1 2
1.695310 0.409718
0.0139954 0.2614200
0.013406855 0.342323710
0.0164444 –
0 4.5455106
3
1 2
1.695310 0.409718
0.0139523 0.2614095
0.013406855 0.342323710
0.0164444 –
0 4.5455106
a
1=0.064; l1(s1)=0.08.
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
1180
3.2. Non homogeneous reactor To analyze a problem where a larger number of modes are necessary, we consider a non homogeneous reactor. This reactor is made of three different materials. Table 2 summarizes the nuclear properties of these materials, as well as other data that define the problem. The core is divided in 72 identical nodes (338), with total dimensions 9090240 cm3, with the layout shown in Fig. 3. During the transient analyzed the relative power increases during 2 s and then it decreases during the following 4 s. This is obtained supposing that f1 and f2 of material 1 are the following functions of time 8 0:122 > > t 04t42 1þ < 0:013406855 0:8 f1 ðtÞ ¼ 0:122 > > : 0:013815764 ð t 2Þ 24t46 1 0:8 and
f2 ðtÞ ¼
8 > > < 0:342323710 > > : 0:352764583
0:122 t 1þ 0:8 0:122 ð t 2Þ 1 0:8
04t42 24t46
This transient problem clearly represents an irreal situation, but it was designed to analyze the capabilities of the code to model fast and strong changes in the fission power. The non-homogeneity of this reactor is mainly due to fast changes in the fission cross-sections in region 1.
Fig. 3. Non homogeneous reactor.
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
1181
The integration time step used has been of 1103 s. To study the performance of the method for the non homogeneous reactor, we have integrated this problem using different number of modes with different updating times. In Fig. 4, we compare the power evolution obtained for the transient using different number of modes which have not been updated, with the power evolution obtained with nokin3d, taken as a reference. We observe that the larger is the number of modes used the better solution is achieved. In Fig. 5, we show a detail of the power evolution calculated using 3 modes without updating the modes and updating the modes each 1 s. We observe that the process of updating the modes increases considerably the accuracy of the obtained solution. The updating process is an expensive process from the computational point of view, thus it will be necessary to find an equilibrium between the number of modes and its updating time to optimize the performance of the method. The main differences in the transient appear in the height of the maximum of power achieved in the transient. In Table 3, we show the relative percentage of error
Fig. 4. Power evolution for the non homogeneous reactor without updating the set of modes. Table 3 Relative percentage of error for the height of the maximum of the power for the non homogeneous reactor transient Number of modes
No updating
2s
1s
0.1 s
1 3 5
58.06 19.43 5.28
10.18 4.07 1.54
4.13 1.48 3.89102
2.74 0.16 1.78102
1182
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
Fig. 5. Detail of the power evolution for the non homogeneous reactor obtained using 3 modes without updating and updating each 1 s the set of modes.
for the maximum height for different combinations of number of modes and updating times. As we have already mentioned the relative error decreases as the number of modes increases and the updating time for the modes decreases. In Table 4 we show the CPU times obtained with a CONVEX Exemplar S-Class computer using the quasi-static method (Ott, 1969; Verdu´ et al., 1995) which is a fast method to integrate the neutron diffusion equation. In Table 5 we present the CPU times spent to integrate the transient using the nodal modal method. We observe that concerning the CPU time the quasi-static method performs better than the nodal modal method, but the nodal modal method provides modal information on the transient.
Table 4 CPU time for non homogeneous reactor integrated with the quasi-static method using different updating times for the shape function Updating time
CPU time
2s 1s 0.1 s
3 min 56 s 3 min 58 s 4 min 5 s
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
1183
Table 5 CPU time for non homogeneous reactor for different number of modes and updating times using modkin Number of modes
No updating
2s
1s
0.1 s
1 3 5
5 min 5 s 6 min 50 s 10 min 18 s
5 min 6 s 6 min 51 s 10 min 22 s
5 min 8 s 6 min 52 s 10 min 24 s
5 min 26 s 7 min 22 s 10 min 57 s
Fig. 6. Geometry of Langenbuch reactor.
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
1184
Table 6 Group constants for the Langenbuch problem Region
Group
Dg (cm)
ag (cm1)
fg (cm1)
12 (cm1)
g (cm/s)
Fuel 1
1 2
1.423913 0.3563060
0.01040206 0.08766217
0.006477691 0.1127328
0.01755550 –
1.25107 2.5105
Fuel 2
1 2
1.425611 0.3505740
0.01099263 0.09925634
0.007503284 0.1378004
0.01717768 –
1.25107 2.5105
Absorber
1 2
1.423913 0.3563060
0.01095206 0.09146217
0.006477691 0.11273228
0.01755550 –
1.25107 2.5105
Reflector
1 2
1.634227 0.2640020
0.002660573 0.04936351
0.0 0.0
0.02759693 –
1.25107 2.5105
Del. Prec.
Group 1
Group 2
Group 3
Group 4
Group 5
Group 6
i li (s1)
0.000247 0.0127
0.0013845 0.0317
0.001222 0.115
0.0026455 0.311
0.000832 1.4
0.000169 3.87
3.3. Langenbuch reactor The third problem analyzed was introduced by Langenbuch et al. (1977) and represents an operational transient in a small LWR core with 77 fuel assemblies and two types of fuel. Fig. 6, shows the geometry of the reactor model. The numerical solution has been obtained by many authors but, as far as we know, this problem has not a reference solution. In Table 6, the macroscopic cross-sections and other data defining the transient are shown. All the calculations for the reactor have been performed with a whole core model, using 1190 nodes. The value obtained for the fundamental eigenvalue was 0.9991484. The transient was initiated by the withdrawal of a bank of four partially inserted control rods at a rate of 3 cm/s over 04t426.7 s. A second bank of control rods is inserted at the same rate over 7.54t447.5 s. The transient is followed during 60 s, integrating the equations with a maximum time-step of 0.125 s. In Fig. 7 we show the power evolution for this transient using 5 modes without updating the set of modes compared with the solution obtained by nokin3d (Verdu´ et al., 1995) taken as a reference. We observe that it would be necessary to use a number of modes larger than 5 to increase the accuracy of the obtained solution. But to obtain a large number of modes for this three dimensional problem is very expensive from the computational point of view. In Fig. 8 we show a detail of the power evolution for the transient calculated using 3 modes. We compare the solution obtained without updating the modes and updating the modes each 10 s with the reference solution. We observe that the
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
1185
Fig. 7. Detail of the power evolution for Langenbuch transient using 5 modes without updating the set of modes.
Fig. 8. Detail of the power evolution for Langenbuch transient using 3 modes and updating the set of modes.
1186
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
updating process for the modes also increases the accuracy of the obtained solution for this transient. In Table 7 we show the relative percentage of error calculated for the height of the maximum of the power for Langenbuch transient. We observe that this problem can be solved essentially using the fundamental mode, updating this mode each certain time step. This suggests that the developed nodal modal method generalizes the quasi-static method for the neutron diffusion equation (Ott, 1969; Verdu´ et al., 1995) where it is supposed that the neutronic flux admits a factorization in two functions: the amplitude function which is only a time dependent function, and the shape function which takes account of the spatial dependence, that changes slowly in time. In Table 8 we compare the CPU times necessary to integrate the Langenbuch transient with the nodal modal method using one mode and the quasi-static method. For this case we observe that both methods perform a comparable efficiency from the CPU time point of view. 3.4. Leibstadt reactor Finally, with the aim of a better understanding of the BWR instabilities development process, we have studied the stability behaviour of two operational points of the nuclear power plant (NPP) Leibstadt (cycles 7 and 10). In the first of the operational points studied a limit cycle corresponding to an in-phase oscillation was developed (cycle 10), and for the other the limit cycle obtained corresponds to an out-of-phase oscillation (cycle 7) (Hennig, 1998). With ramona (Wulff et al., 1993) code it is possible to obtain detailed information regarding the state of the reactor
Table 7 Relative percentage of error for the height of the maximum of the power for Langenbuch transient Number of Modes
No updating
10 s
5s
1s
1 3 5
6.538 6.451 6.451
1.259 1.105 1.104
0.313 0.299 0.298
6.18102 5.77102 5.71102
Table 8 CPU time for Lagenbuch reactor for different updating times using modkin and the quasi-static method Updating time
modkin
Quasi-static
No update 10 s 5s 2s 1s
2 2 3 3 4
– 3 3 3 3
min min min min min
34 s 47 s 1s 38 s 42 s
min min min min
8s 35 s 34 s 53 s
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
1187
for each integration time step. In this way, we have obtained, for each time step, the nuclear cross-sections for each node of the 64825 NPP Leibstadt (KKL) core nodes (without reflector) having each node a volume of 15.2415.2415.24 cm3, as well as other data defining the transients. To study the transients we have used the nuclear cross-sections provided by ramona, we have added the nuclear constants for the reflector nodes surrounding the reactor and the whole set is used as an input for the modkin3d code. For more details on these test transients see (Blomstrand, 1992; Wiktor, 1997; Miro´ et al., 2000). For both transients, we have used 5 modes and due to its smooth spatial change it is not necessary to update them. The results have been compared with the ones obtained with the code ramona3-12, making use of a time step of 0.1 s. The power evolution for the in-phase case is represented in Fig. 9, where we observe the good agreement among the results calculated with both codes and clearly shows the maintained limit cycle. Fig. 10 compares the results obtained for the out-of-phase case showing the power evolution obtained with the codes, but in this case the power oscillates at the beginning of the transient, but it tends to a steady state as the time of simulation advances. Figs. 11 and 12 show the 5 mode patterns of the two studied cases in two axial locations, one in the lower part of the reactor (plane 4) and the other in the upper one (plane 24). The first and the second harmonics correspond to azimuthal modes while the third harmonic is related with the first axial mode.
Fig. 9. Power evolution for Leibstadt in-phase case.
1188
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
Fig. 10. Power evolution for Leibstadt out-of-phase case.
To characterize the two studied transients as in-phase and out-of-phase, respectively and also to study the importance of the different modes during the transients, we have also calculated the time dependent amplitudes nl(t) of the vector (t), for the fundamental mode and the first four harmonics, also with the code modkin. Figs. 13 and 14 show the obtained evolution for n1(t), and n2(t), n3(t) and n4(t), respectively, for the in-phase case. We observe that the amplitude n1(t) is clearly the dominant oscillation since n2(t), n3(t) vanish in time, while the amplitude n4(t) is smaller than n1(t), but not rejectable, which makes the contribution of the third subcritical mode very important. In particular, because its spatial distribution along the axial direction (see Fig. 11), there is a phase shift between the oscillation in the lower and the upper part of the reactor, this has been confirmed studying the LPRM response simulated also with ramona and modkin, showing that there is a phase shift of approximately 90 between the LPRM response taken from an axial level below and above the middle plane. The amplitude n5(t) is not represented, but also vanishes in time. Figs. 15 and 16 show the results for the out-of-phase case. The fundamental amplitude n1(t) oscillates at the begining of the transient but later it becomes stable, showing the same behaviour of the relative power (Fig. 10). The amplitudes n2(t) and n3(t) have a phase shift of 180 . They increase and clearly dominate the oscillation as the transient goes on, especially n2(t), which, approximately, doubles n3(t). Amplitudes n4(t) and n5(t) are not shown but their value is closed to zero during all the transient. Furthermore, for this case (Fig. 12), as the second harmonic is rotated 90 respect to the first harmonic, it appears a double oscillation, shifted 180 and
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
Fig. 11. In-phase case: Lambda modes profiles at different reactor levels.
1189
1190
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
Fig. 12. Out-of-phase case: Lambda modes profiles at different reactor levels.
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
Fig. 13. n1(t) Amplitude for in-phase case.
Fig. 14. n2(t), n3(t) and n4(t) amplitudes for in-phase case.
1191
1192
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
Fig. 15. n1(t) Amplitude for out-of-phase case.
Fig. 16. n2(t) and n3(t) amplitudes for out-of-phase case.
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
1193
with double amplitude in the two reactor halves defined by the first harmonic symmetry line, giving a rotation effect in the oscillation with respect to the axial axis.
4. Conclusions In this paper we have proposed a nodal modal method to integrate the neutron diffusion equation. This method starts from the discretized equations result of applying a nodal collocation method to the neutron diffusion equation in the approximation of two groups of energy and assumes that the vector of coefficients, solution of these equations, can be expressed as a finite linear combinations of eigenfunctions, which are solutions of a static neutron diffusion equation associated with a given configuration of the reactor core. Since the calculation of a large number of eigenfunctions for a realistic problem is a very expensive process from the computational point of view, we have proposed an updating process for the eigenfunctions used as a set of basis functions. This process permits to integrate the neutron diffusion equation using a moderate amount of modes. This updating process generalizes the quasi-static method for the neutron diffusion equation and allows to decrease the frequency of the actualization of the shape function without lose of accuracy. This nodal modal method is reasonably accurate and adequate to study transients with a clear physical interpretation for the amplitudes of the neutron flux modal expansion, as the ones appearing in the in-phase or out-of-phase oscillations of boiling water reactors.
References Blomstrand, J., 1992. The KKL Core Stability Test, Conducted in September 1990. (ABB Report BR91245). Ginestar, G., Verdu´, G., Vidal, V., Bru, R., Marı´n, J., Mun˜oz-Cobo, J.L., 1998. High order backward discretization for the neutron diffusion equation. Ann. Nucl. Energy 25 (1), 47–64. He´bert, A., 1987. Development of the nodal collocation method for solving the neutron diffusion equation. Ann. Nucl. Energy 14 (10), 527–541. Hennig, D., 1998. Stability Analysis KKL (Part 2). PSI Technical Report TM-41-98-39 (in German). Henry, A.F., 1982. Nuclear-reactor analysis. The MIT Press, Cambridge, Massachussetts. Hindmarsh, A.C., 1983. ODEPACK, A Systematized Collection of ODE Solvers in Scientific Computing. North-Holland. Langenbuch, S., Maurer, W., Werner, W., 1997. Coarse-mesh flux expansion method for the analysis of space-time effects in large light water reactors cores. Nuclear Science and Engineering 63, 437–456. Miro´, R., Ginestar, D., Hennig, D., Verdu´, G., 2000. On the regional oscillation phenomenon in BWR’s. Progress in Nuclear Energy 36 (2), 189–229. Miro´, R., 2001. Me´todos Modales para el Estudio de Inestabilidades en Reactores Nucleares BWR. PhD thesis, Polytechnic University of Valencia (to be presented). Ott, K.O., 1969. Accuracy of the quasistatic treatment of spatial reactor kinetics. Nuclear Science and Engineering 36, 402–411. Stacey Jr., W.M., 1969. Space-time nuclear reactor kinetics. Academic Press, New York.
1194
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
Verdu´, G., Ginestar, D., Vidal, V., Mun˜oz-Cobo, J.L., 1994. 3D l-modes of the neutron-diffusion equation. Ann. Nucl. Energy 21 (7), 405–421. Verdu´, G., Ginestar, D., Vidal, V., Mun˜oz-Cobo, J.L., 1995. A consistent multidimensional nodal method for transient calculations. Ann. Nucl. Energy 22 (6), 395–411. Verdu´, G., Miro´, R., Ginestar, D., Vidal, V., 1999. The implicit restarted Arnoldi method, an efficient alternative to solve the neutron diffusion equation. Ann. Nucl. Energy 26 (7), 579–593. Wiktor, C.G., 1997. KKL Cycle 7 Stability Test-Core Simulation Input Data. (KKL Report BET/97/ 111). Wulff, W., Cheng, H.S., Diamond, D.J., Khatib-Rhabar, M., 1993. A Description and Assessment of RAMONA-3B Mod.0 Cycle 4: A Computer Code with Three-Dimensinal Neutron Kinetic for BWR System Transients. NUREG/CR-3664, BNL-NUREG-51746 (see also: RAMONA-3 Users’ Manual, Scandpower).
A nodal modal method for the neutron diffusion equation. Application to BWR instabilities analysis R. Miro´a, D. Ginestarb, G. Verdu´a,*, D. Hennigc a
Departamento de Ingenierı´a Quı´mica y Nuclear, Universidad Polite´cnica de Valencia, Camı´ de Vera, 14, 46022, Valencia, Spain b Departamento de Matema´tica Aplicada, Universidad Polite´cnica de Valencia, Camı´ de Vera, 14, 46022, Valencia, Spain c Laboratory for Reactor Physics and Systems Behaviour (LRS), Paul Scherrer Institut (PSI), CH-5232 Villigen PSI, Switzerland Received 26 July 2001; received in revised form 22 September 2001; accepted 24 September 2001
Abstract Fast codes, capable of dealing with three-dimensional geometries, are needed to be able to simulate spatially complicated transients in a nuclear power reactor. In this paper, we propose a modal method to integrate the neutron diffusion equation in which the spatial part has been previously dicretized using a nodal collocation method. For the time integration of the resulting system of differential equations it is supposed that the solution can be expanded as a linear combination of the dominant Lambda modes associated with a static configuration of the reactor core and, using the eigenfunctions of the adjoint problem, a system of differential equations of lower dimension is obtained. This system is integrated using a variable time step implicit method. Furthermore, for realistic transients, it would be necessary to calculate a large amount of modes. To avoid this, the modal method has been implemented making use of an updating process for the modes at each certain time step. Five transients have been studied: a homogeneous reactor, a non-homogeneous reactor, the 3D Langenbuch reactor and two transients related with in-phase and out-of-phase oscillations of Leibstadt NPP. The obtained results have been compared with the ones provided by a method based on a one-step backward discretization formula. # 2002 Elsevier Science Ltd. All rights reserved.
* Corresponding author. Tel.: +34-9-6387-7635; fax: +34-9-6387-7639. E-mail address: [email protected] (G. Verdu´). 0306-4549/02/$ - see front matter # 2002 Elsevier Science Ltd. All rights reserved. PII: S0306-4549(01)00103-7
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
1172
1. Introduction To simulate the behaviour of a nuclear power reactor it is necessary to be able to integrate the time-dependent neutron diffusion equation inside the reactor core. In this way, it is interesting to develop fast and reliable methods to integrate this equation. In this work, we will develop a modal method based on the assumption that the solution for the neutron flux of the studied transient can be expressed approximately as a linear combination of the dominant eigenfunctions of a static auxiliary problem, known as the Lambda modes, associated with a static reactor configuration (Henry, 1982). To calculate these auxiliary eigenfunctions we have used a nodal serendipity collocation method (He´bert, 1987) to discretize the equations, and an iterative method to obtain the dominant eigenvalues and their corresponding eigenfunctions (Verdu´ et al., 1994, 1999). This modal method is interesting to study transients where the number of necessary eigenmodes is not very large, and it can also be useful to analyze transients with a physical interpretation for the different amplitudes of the neutron flux expansion, as the ones appearing in the study of in-phase and out-of-phase oscillations of boiling water reactors (Miro´ et al., 2000). In the case of one eigenmode, the method could be considered as a version of the quasi-static method to integrate the neutron diffusion equation (Ott, 1969).
2. Nodal modal method Under general assumptions (Stacey, 1969), the neutronic population inside a nuclear power reactor can be modeled by the neutron diffusion equation in the approximation of two energy groups. This model is of the form,
K X : 1 þ L ¼ ð1 ÞM þ lk Ck ;
ð1Þ
k¼1
: C k ¼ k f 1 f 2 lk Ck ; k ¼ 1; . . . K;
ð2Þ
where, K is the number of delayed neutron precursors groups considered, 2 L¼4
~ þ a1 þ 12 ~ D1 r r 12
0
5; ~ þ a2 ~ D2 r r
and M¼
f 1 0
3
f 2 1 1 ; ¼ ; ¼ : 0 2 0
1
21 6 ¼ 4 1 0
0
3
7 1 5; 2
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
1173
The diffusion constants and cross-sections, Dg, 12, ag, fg, g=1,2, appearing in the equations depend on the reactor materials, that is, they are position dependent functions. To solve the problem, we begin with an approximation, consisting of discretizing the reactor in N nodes, where the nuclear properties are supposed to be constant. Next, a nodal collocation method (He´bert, 1987; Verdu´ et al., 1994), that assumes a truncated expansion of the neutron flux in terms of Legendre polynomials is used. This method permits to discretize the spatial part of the equations, approximating Eqs. (1) and (2) by a set of ordinary differential equations of the form
1
:
þ L ¼ ð1 ÞM þ
: XC k ¼ k M lk XCk ;
K X lk XCk ;
ð3Þ
k¼1
ð4Þ
where and Ck are vectors, whose components are the corresponding coefficients of the Legendre polynomials expansions of and Ck in each node, at each time step. L and M have the following block structure,
0 L11 M11 M12 ; L¼ ; M¼ L21 L22 0 0 and we have introduced the matrix X, I X¼ : 0 To integrate Eqs. (3) and (4) using the modal approximation, we assume that can be expressed approximately as Md X nl ðtÞ l ;
ð5Þ
l¼1
where l, l=1, . . ., Md are the eigenvectors associated with the Md dominant eigenvalues of a static problem of the form (Verdu´ et al., 1994) L0
l
¼
1 M0 l ; kl
ð6Þ
where the matrices L0 and M0 are the matrices L and M of Eqs. (3) and (4), for a given configuration of the reactor core. This problem is known as the Lambda modes problem (Henry, 1982; Verdu´ et al., 1994). Since L0 is not a hermitic matrix, we introduce the adjoint problem (Henry, 1982), Ly0
y l
¼
1 y M kl 0
y l;
ð7Þ
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
1174
y m,
to obtain a set of modes, n, and their adjoints, relation, y y m ; M0 n ¼ m ; M0 m n;m ¼ Nm n;m ;
satisfying the biorthogonality
ð8Þ
where n,m is the Kronecker delta function. Multiplying Eqs. (3) and (4) by ym , writing L ¼ L0 þ L;
M ¼ M0 þ M;
and making use of expansion (5), and Eq. (6), we obtain the following set of equations Md X
y m;
1
l¼1
x
Md X
Md X d 1 nl þ l dt k l¼1 l
y m ; M0
l
nl þ ð 1 Þ
l¼1
d dt Md X
y m ; M0
Md X
l
nl þ
y m ; M0
y m ; L l
nl ¼ ð1 Þx
l¼1 y m ; M
l
nl þ
l¼1 y m ; XCk
Md X
K X lk
y m ; XCk
;
k¼1
¼ k l
l¼1
nl þ k
Md X
y m ; M
l
nl lk
y m ; XCk
:
ð9Þ
l¼1
Using the biorthogonality relationship (8) and introducing the notation,
ml ¼ ym ; 1 l ; ALml ¼ ym ; L l ; y Cmk ¼ ym ; XCk ; AM ml ¼ m ; M l ; and the mode m reactivity, defined as m ¼
km 1 ; km
Eqs. (9) can be expressed as, Md Md Md K X X X X d
ml nl ¼ ð m ÞNm nm þ ð1 Þ AM n ALml nl þ lk Cmk ; ml l dt l¼1 l¼1 l¼1 k¼1 Md X d Cmk ¼ k Nm nm þ k AM ml nl lk Cmk ; dt l¼0
k ¼ 1; . . . ; K:
ð10Þ
Varying the index m from 1 to Md, we can rewrite these equations in the following matrix form
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
1175
d½n ¼ ½ 1 ½ I½N½n þ ð1 Þ½ 1 AM ½n ½ 1 AL ½n dt K X þ lk ½ 1 ½Ck ;
ð11Þ
k¼1
d½Ck ¼ k ½N½n þ k AM ½n lk ½Ck ; k ¼ 1; . . . ; K: dt Considering explicitly the K equations for the delayed neutron precursors, Eqs. (11) can be grouped as the following system of differential equations d N ¼ AN ; dt where 2 3 n 6 C1 7 7 N ¼6 4 .. 5; . CK and
2 6 6 A¼6 4
ð12Þ
½ 1 ½ I½N þ ð1 ÞAM AL 1 ½ N þ A M .. . M K N þ A
½ 1 l1 ½ 1 lK l1I
0 . .. . . . ... 0 lK I
3 7 7 7: 5
If the number of dominant modes needed in expansion (5), Md, is not very large, the dimension of the new system of differential Eqs. (12) is quite smaller than the one of the initial system (3) and (4). To solve the system (12) we have used a variable step implicit method (Hindmarsh, 1983) because of the stiffness of the differential equations. To obtain the initial conditions for the time integration, we start from a critical configuration of the reactor. To obtain this critical reactor configuration we solve the Lambda modes problem for a given initial configuration, L0
1
¼
1 M0 k1
1;
for the fundamental mode. Dividing the fission cross-sections by the fundamental eigenvalue, k1, we obtain L0
1
¼ MCrit
1;
ð13Þ
where MCrit is a matrix which components are the components of M divided by k1. In this way we obtain that 1 is a stationary solution of Eqs. (3) and (4) for a critical configuration of the reactor. The subsequent fission cross-sections of Eqs. (12) will be also divided by the eigenvalue k1 to calculate the time evolution of N .
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
1176
Eq. (13), together with the equation 0 ¼ k MCrit
1
lk XCk ;
ð14Þ
constitute the set of equations defining the steady state associated with the critical reactor configuration. Hence, the components of vector [n] for the initial condition are n1 ¼ 1;
nl ¼ 0;
l ¼ 2; . . . ; Md :
Multiplying (14) by y1 we have D E D E 0 ¼ k y1 ; MCrit 1 lk y1 ; XCk ; Thus, the initial values for [Ck] are given by E k D y Crit C1;k ¼ 1 ; 1; M lk Cl;k ¼ 0; l ¼ 2; . . . ; Md : 2.1. Modes updating For realistic transients, the nuclear cross-sections are time dependent functions and to use the modal method proposed above, it would be necessary to calculate a large amount of modes. This is prohibitive from the computational point of view. Thus, instead of calculating a large amount of modes, we calculate only a small number of them and we update these modes each certain time step, ti. In this way, to integrate the neutron diffusion equation in the time interval [ti, ti+1], we use the Lambda modes associated with the problem 1 Li il ¼ i Mi il ; kl where Li and Mi are the matrices associated with the reactor configuration at time ti. Hence, the differential equations to integrate are the following ones, Md X
iy m;
1
l¼1
x
Md X
Md i d i X 1 nl þ l dt ki l¼1 l
iy i m; M
i l
nil þ ð1 Þ
l¼1
d dt Md X l¼1
Md X
iy i m; M
i l
nil þ
iy i m; M
i l
iy i m ; L
i l
nil ¼ ð1 Þx
l¼1 iy i m ; M
i l
nil þ
l¼1 iy m ; XCk
Md X
K X lk
iy m ; XCk
;
k¼1
¼ k Md X nl þ k
iy i m ; M
i l
nil lk
iy m ; XCk
:
l¼1
the operators Li and Mi change along the time interval.
ð15Þ
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
1177
To calculate the initial conditions for nim ðti Þ; we reconstruct the vector ðti Þ ¼
Md X ni1 l ðti Þ
i1 l ;
l¼1
from the variables ni1 l ðti Þ and nim ðti Þ as 1
nim ðti Þ ¼ D
iy i m; M
i m
E
iy i m; M
i1 l
obtained in the last time step, and we calculate
ðti Þ :
ð16Þ
To obtain the variables related to the concentration of the delayed neutron pre i1y cursors, iy ; XC ð Þ; in terms of ; XC ð Þ; which are known quantities t t k i k i m m from the integration in the previous time interval, we suppose that approximately we can write iy m
Md X aml
i1y ; l
l¼1
where aml ¼ D
iy i1 i1 m; M l E; i1y i1 i1 ; M l l
and this yields
iy m ; XCk
ðti Þ
Md D X aml
i1y ; XCk l
E
ðti Þ;
ð17Þ
l¼1
which completes the initial conditions to integrate (15) in [ti, ti+1].
3. Numerical results The exposed nodal modal method has been implemented in a code written in fortran and called modkin, which calculates the dominant Lambda modes associated with a given configuration of the reactor core, making use of the implicit restarted Arnoldi method (Verdu´ et al., 1999). To test the performance of the method, we have studied five transients: one associated with a homogeneous reactor without delayed neutron precursors, another one associated with a non-homogeneous reactor with one group of delayed neutron precursors, the Langenbuch problem (Langenbuch et al., 1977) and two transients corresponding to in-phase and out-of-phase oscillations, whose nuclear data where obtained using a ramona thermal-hydraulic model of Leibstadt NPP. The former three transients have been used as validation of the code. For the first transient we have an exact analytical solution (Miro´, 2001).
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
1178
For the other examples, the results obtained with the nodal modal method have been compared with the ones obtained with a code that integrates Eqs. (3) and (4) using a temporal discretization based on the one-step backward difference formula, called nokin3d (Verdu´ et al., 1995; Ginestar et al., 1997) which is an accurate method to integrate the neutron diffusion equation. In all the numerical examples studied, we have only considered changes in the cross-sections. This kind of transients can be analyzed using only spatial modes. More general perturbations should be addressed with more general models as the multigroup (more than two) neutron diffusion equation or the neutron transport equation. 3.1. Homogeneous reactor This problem consists of a 10060180 cm3 homogeneous reactor, divided in 36 equal nodes (336), as it is shown in Fig. 1. These nodes correspond to only one material, which nuclear properties are given in Table 1. The transient consists of a step perturbation on the thermal fission cross-section f 2 ðtÞ ¼ 0f2 þ f 2 ¼ 0f 2 þ 1 104 HðtÞ; where H(t) is the unit step function, and 0f2 is the initial thermal fission cross-section, which is chosen imposing that the initial configuration of the reactor is critical.
Fig. 1. Homogeneous reactor. Table 1 Group constants for the homogeneous reactor Group
Dg (cm)
ag (cm1)
fg (cm1)
12 (cm1)
1/g (s/cm)
1 2
1.28205 0.666667
0.01 0.1
0.01 0.1097683675455275
0.075 –
1107 1105
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
1179
That is, the initial configuration of the reactor is a steady state configuration described by the first eigenvector of a Lambda modes problem associated with the reactor such that the first eigenvalue of this problem is equal to 1 (Henry, 1982). This problem will be used as a first test for the nodal modal method. The transient has been followed during 1 s using an integration time step of 1103 s, and the results obtained with the nodal modal method using one mode have been compared with the analytical solution for the relative neutronic power evolution. The results are presented in Fig. 2, showing a very good agreement.
Fig. 2. Power evolution for the homogeneous reactor.
Table 2 Nuclear properties of the non homogeneous reactora Mat.
Group
Dg (cm)
ag (cm1)
fg (cm1)
12 (cm1)
1/g (s/cm)
1
1 2
1.695310 0.409718
0.0139530 0.2614097
0.013406855 0.342323710
0.0164444 –
0 4.5455106
2
1 2
1.695310 0.409718
0.0139954 0.2614200
0.013406855 0.342323710
0.0164444 –
0 4.5455106
3
1 2
1.695310 0.409718
0.0139523 0.2614095
0.013406855 0.342323710
0.0164444 –
0 4.5455106
a
1=0.064; l1(s1)=0.08.
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
1180
3.2. Non homogeneous reactor To analyze a problem where a larger number of modes are necessary, we consider a non homogeneous reactor. This reactor is made of three different materials. Table 2 summarizes the nuclear properties of these materials, as well as other data that define the problem. The core is divided in 72 identical nodes (338), with total dimensions 9090240 cm3, with the layout shown in Fig. 3. During the transient analyzed the relative power increases during 2 s and then it decreases during the following 4 s. This is obtained supposing that f1 and f2 of material 1 are the following functions of time 8 0:122 > > t 04t42 1þ < 0:013406855 0:8 f1 ðtÞ ¼ 0:122 > > : 0:013815764 ð t 2Þ 24t46 1 0:8 and
f2 ðtÞ ¼
8 > > < 0:342323710 > > : 0:352764583
0:122 t 1þ 0:8 0:122 ð t 2Þ 1 0:8
04t42 24t46
This transient problem clearly represents an irreal situation, but it was designed to analyze the capabilities of the code to model fast and strong changes in the fission power. The non-homogeneity of this reactor is mainly due to fast changes in the fission cross-sections in region 1.
Fig. 3. Non homogeneous reactor.
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
1181
The integration time step used has been of 1103 s. To study the performance of the method for the non homogeneous reactor, we have integrated this problem using different number of modes with different updating times. In Fig. 4, we compare the power evolution obtained for the transient using different number of modes which have not been updated, with the power evolution obtained with nokin3d, taken as a reference. We observe that the larger is the number of modes used the better solution is achieved. In Fig. 5, we show a detail of the power evolution calculated using 3 modes without updating the modes and updating the modes each 1 s. We observe that the process of updating the modes increases considerably the accuracy of the obtained solution. The updating process is an expensive process from the computational point of view, thus it will be necessary to find an equilibrium between the number of modes and its updating time to optimize the performance of the method. The main differences in the transient appear in the height of the maximum of power achieved in the transient. In Table 3, we show the relative percentage of error
Fig. 4. Power evolution for the non homogeneous reactor without updating the set of modes. Table 3 Relative percentage of error for the height of the maximum of the power for the non homogeneous reactor transient Number of modes
No updating
2s
1s
0.1 s
1 3 5
58.06 19.43 5.28
10.18 4.07 1.54
4.13 1.48 3.89102
2.74 0.16 1.78102
1182
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
Fig. 5. Detail of the power evolution for the non homogeneous reactor obtained using 3 modes without updating and updating each 1 s the set of modes.
for the maximum height for different combinations of number of modes and updating times. As we have already mentioned the relative error decreases as the number of modes increases and the updating time for the modes decreases. In Table 4 we show the CPU times obtained with a CONVEX Exemplar S-Class computer using the quasi-static method (Ott, 1969; Verdu´ et al., 1995) which is a fast method to integrate the neutron diffusion equation. In Table 5 we present the CPU times spent to integrate the transient using the nodal modal method. We observe that concerning the CPU time the quasi-static method performs better than the nodal modal method, but the nodal modal method provides modal information on the transient.
Table 4 CPU time for non homogeneous reactor integrated with the quasi-static method using different updating times for the shape function Updating time
CPU time
2s 1s 0.1 s
3 min 56 s 3 min 58 s 4 min 5 s
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
1183
Table 5 CPU time for non homogeneous reactor for different number of modes and updating times using modkin Number of modes
No updating
2s
1s
0.1 s
1 3 5
5 min 5 s 6 min 50 s 10 min 18 s
5 min 6 s 6 min 51 s 10 min 22 s
5 min 8 s 6 min 52 s 10 min 24 s
5 min 26 s 7 min 22 s 10 min 57 s
Fig. 6. Geometry of Langenbuch reactor.
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
1184
Table 6 Group constants for the Langenbuch problem Region
Group
Dg (cm)
ag (cm1)
fg (cm1)
12 (cm1)
g (cm/s)
Fuel 1
1 2
1.423913 0.3563060
0.01040206 0.08766217
0.006477691 0.1127328
0.01755550 –
1.25107 2.5105
Fuel 2
1 2
1.425611 0.3505740
0.01099263 0.09925634
0.007503284 0.1378004
0.01717768 –
1.25107 2.5105
Absorber
1 2
1.423913 0.3563060
0.01095206 0.09146217
0.006477691 0.11273228
0.01755550 –
1.25107 2.5105
Reflector
1 2
1.634227 0.2640020
0.002660573 0.04936351
0.0 0.0
0.02759693 –
1.25107 2.5105
Del. Prec.
Group 1
Group 2
Group 3
Group 4
Group 5
Group 6
i li (s1)
0.000247 0.0127
0.0013845 0.0317
0.001222 0.115
0.0026455 0.311
0.000832 1.4
0.000169 3.87
3.3. Langenbuch reactor The third problem analyzed was introduced by Langenbuch et al. (1977) and represents an operational transient in a small LWR core with 77 fuel assemblies and two types of fuel. Fig. 6, shows the geometry of the reactor model. The numerical solution has been obtained by many authors but, as far as we know, this problem has not a reference solution. In Table 6, the macroscopic cross-sections and other data defining the transient are shown. All the calculations for the reactor have been performed with a whole core model, using 1190 nodes. The value obtained for the fundamental eigenvalue was 0.9991484. The transient was initiated by the withdrawal of a bank of four partially inserted control rods at a rate of 3 cm/s over 04t426.7 s. A second bank of control rods is inserted at the same rate over 7.54t447.5 s. The transient is followed during 60 s, integrating the equations with a maximum time-step of 0.125 s. In Fig. 7 we show the power evolution for this transient using 5 modes without updating the set of modes compared with the solution obtained by nokin3d (Verdu´ et al., 1995) taken as a reference. We observe that it would be necessary to use a number of modes larger than 5 to increase the accuracy of the obtained solution. But to obtain a large number of modes for this three dimensional problem is very expensive from the computational point of view. In Fig. 8 we show a detail of the power evolution for the transient calculated using 3 modes. We compare the solution obtained without updating the modes and updating the modes each 10 s with the reference solution. We observe that the
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
1185
Fig. 7. Detail of the power evolution for Langenbuch transient using 5 modes without updating the set of modes.
Fig. 8. Detail of the power evolution for Langenbuch transient using 3 modes and updating the set of modes.
1186
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
updating process for the modes also increases the accuracy of the obtained solution for this transient. In Table 7 we show the relative percentage of error calculated for the height of the maximum of the power for Langenbuch transient. We observe that this problem can be solved essentially using the fundamental mode, updating this mode each certain time step. This suggests that the developed nodal modal method generalizes the quasi-static method for the neutron diffusion equation (Ott, 1969; Verdu´ et al., 1995) where it is supposed that the neutronic flux admits a factorization in two functions: the amplitude function which is only a time dependent function, and the shape function which takes account of the spatial dependence, that changes slowly in time. In Table 8 we compare the CPU times necessary to integrate the Langenbuch transient with the nodal modal method using one mode and the quasi-static method. For this case we observe that both methods perform a comparable efficiency from the CPU time point of view. 3.4. Leibstadt reactor Finally, with the aim of a better understanding of the BWR instabilities development process, we have studied the stability behaviour of two operational points of the nuclear power plant (NPP) Leibstadt (cycles 7 and 10). In the first of the operational points studied a limit cycle corresponding to an in-phase oscillation was developed (cycle 10), and for the other the limit cycle obtained corresponds to an out-of-phase oscillation (cycle 7) (Hennig, 1998). With ramona (Wulff et al., 1993) code it is possible to obtain detailed information regarding the state of the reactor
Table 7 Relative percentage of error for the height of the maximum of the power for Langenbuch transient Number of Modes
No updating
10 s
5s
1s
1 3 5
6.538 6.451 6.451
1.259 1.105 1.104
0.313 0.299 0.298
6.18102 5.77102 5.71102
Table 8 CPU time for Lagenbuch reactor for different updating times using modkin and the quasi-static method Updating time
modkin
Quasi-static
No update 10 s 5s 2s 1s
2 2 3 3 4
– 3 3 3 3
min min min min min
34 s 47 s 1s 38 s 42 s
min min min min
8s 35 s 34 s 53 s
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
1187
for each integration time step. In this way, we have obtained, for each time step, the nuclear cross-sections for each node of the 64825 NPP Leibstadt (KKL) core nodes (without reflector) having each node a volume of 15.2415.2415.24 cm3, as well as other data defining the transients. To study the transients we have used the nuclear cross-sections provided by ramona, we have added the nuclear constants for the reflector nodes surrounding the reactor and the whole set is used as an input for the modkin3d code. For more details on these test transients see (Blomstrand, 1992; Wiktor, 1997; Miro´ et al., 2000). For both transients, we have used 5 modes and due to its smooth spatial change it is not necessary to update them. The results have been compared with the ones obtained with the code ramona3-12, making use of a time step of 0.1 s. The power evolution for the in-phase case is represented in Fig. 9, where we observe the good agreement among the results calculated with both codes and clearly shows the maintained limit cycle. Fig. 10 compares the results obtained for the out-of-phase case showing the power evolution obtained with the codes, but in this case the power oscillates at the beginning of the transient, but it tends to a steady state as the time of simulation advances. Figs. 11 and 12 show the 5 mode patterns of the two studied cases in two axial locations, one in the lower part of the reactor (plane 4) and the other in the upper one (plane 24). The first and the second harmonics correspond to azimuthal modes while the third harmonic is related with the first axial mode.
Fig. 9. Power evolution for Leibstadt in-phase case.
1188
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
Fig. 10. Power evolution for Leibstadt out-of-phase case.
To characterize the two studied transients as in-phase and out-of-phase, respectively and also to study the importance of the different modes during the transients, we have also calculated the time dependent amplitudes nl(t) of the vector (t), for the fundamental mode and the first four harmonics, also with the code modkin. Figs. 13 and 14 show the obtained evolution for n1(t), and n2(t), n3(t) and n4(t), respectively, for the in-phase case. We observe that the amplitude n1(t) is clearly the dominant oscillation since n2(t), n3(t) vanish in time, while the amplitude n4(t) is smaller than n1(t), but not rejectable, which makes the contribution of the third subcritical mode very important. In particular, because its spatial distribution along the axial direction (see Fig. 11), there is a phase shift between the oscillation in the lower and the upper part of the reactor, this has been confirmed studying the LPRM response simulated also with ramona and modkin, showing that there is a phase shift of approximately 90 between the LPRM response taken from an axial level below and above the middle plane. The amplitude n5(t) is not represented, but also vanishes in time. Figs. 15 and 16 show the results for the out-of-phase case. The fundamental amplitude n1(t) oscillates at the begining of the transient but later it becomes stable, showing the same behaviour of the relative power (Fig. 10). The amplitudes n2(t) and n3(t) have a phase shift of 180 . They increase and clearly dominate the oscillation as the transient goes on, especially n2(t), which, approximately, doubles n3(t). Amplitudes n4(t) and n5(t) are not shown but their value is closed to zero during all the transient. Furthermore, for this case (Fig. 12), as the second harmonic is rotated 90 respect to the first harmonic, it appears a double oscillation, shifted 180 and
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
Fig. 11. In-phase case: Lambda modes profiles at different reactor levels.
1189
1190
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
Fig. 12. Out-of-phase case: Lambda modes profiles at different reactor levels.
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
Fig. 13. n1(t) Amplitude for in-phase case.
Fig. 14. n2(t), n3(t) and n4(t) amplitudes for in-phase case.
1191
1192
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
Fig. 15. n1(t) Amplitude for out-of-phase case.
Fig. 16. n2(t) and n3(t) amplitudes for out-of-phase case.
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
1193
with double amplitude in the two reactor halves defined by the first harmonic symmetry line, giving a rotation effect in the oscillation with respect to the axial axis.
4. Conclusions In this paper we have proposed a nodal modal method to integrate the neutron diffusion equation. This method starts from the discretized equations result of applying a nodal collocation method to the neutron diffusion equation in the approximation of two groups of energy and assumes that the vector of coefficients, solution of these equations, can be expressed as a finite linear combinations of eigenfunctions, which are solutions of a static neutron diffusion equation associated with a given configuration of the reactor core. Since the calculation of a large number of eigenfunctions for a realistic problem is a very expensive process from the computational point of view, we have proposed an updating process for the eigenfunctions used as a set of basis functions. This process permits to integrate the neutron diffusion equation using a moderate amount of modes. This updating process generalizes the quasi-static method for the neutron diffusion equation and allows to decrease the frequency of the actualization of the shape function without lose of accuracy. This nodal modal method is reasonably accurate and adequate to study transients with a clear physical interpretation for the amplitudes of the neutron flux modal expansion, as the ones appearing in the in-phase or out-of-phase oscillations of boiling water reactors.
References Blomstrand, J., 1992. The KKL Core Stability Test, Conducted in September 1990. (ABB Report BR91245). Ginestar, G., Verdu´, G., Vidal, V., Bru, R., Marı´n, J., Mun˜oz-Cobo, J.L., 1998. High order backward discretization for the neutron diffusion equation. Ann. Nucl. Energy 25 (1), 47–64. He´bert, A., 1987. Development of the nodal collocation method for solving the neutron diffusion equation. Ann. Nucl. Energy 14 (10), 527–541. Hennig, D., 1998. Stability Analysis KKL (Part 2). PSI Technical Report TM-41-98-39 (in German). Henry, A.F., 1982. Nuclear-reactor analysis. The MIT Press, Cambridge, Massachussetts. Hindmarsh, A.C., 1983. ODEPACK, A Systematized Collection of ODE Solvers in Scientific Computing. North-Holland. Langenbuch, S., Maurer, W., Werner, W., 1997. Coarse-mesh flux expansion method for the analysis of space-time effects in large light water reactors cores. Nuclear Science and Engineering 63, 437–456. Miro´, R., Ginestar, D., Hennig, D., Verdu´, G., 2000. On the regional oscillation phenomenon in BWR’s. Progress in Nuclear Energy 36 (2), 189–229. Miro´, R., 2001. Me´todos Modales para el Estudio de Inestabilidades en Reactores Nucleares BWR. PhD thesis, Polytechnic University of Valencia (to be presented). Ott, K.O., 1969. Accuracy of the quasistatic treatment of spatial reactor kinetics. Nuclear Science and Engineering 36, 402–411. Stacey Jr., W.M., 1969. Space-time nuclear reactor kinetics. Academic Press, New York.
1194
R. Miro´ et al. / Annals of Nuclear Energy 29 (2002) 1171–1194
Verdu´, G., Ginestar, D., Vidal, V., Mun˜oz-Cobo, J.L., 1994. 3D l-modes of the neutron-diffusion equation. Ann. Nucl. Energy 21 (7), 405–421. Verdu´, G., Ginestar, D., Vidal, V., Mun˜oz-Cobo, J.L., 1995. A consistent multidimensional nodal method for transient calculations. Ann. Nucl. Energy 22 (6), 395–411. Verdu´, G., Miro´, R., Ginestar, D., Vidal, V., 1999. The implicit restarted Arnoldi method, an efficient alternative to solve the neutron diffusion equation. Ann. Nucl. Energy 26 (7), 579–593. Wiktor, C.G., 1997. KKL Cycle 7 Stability Test-Core Simulation Input Data. (KKL Report BET/97/ 111). Wulff, W., Cheng, H.S., Diamond, D.J., Khatib-Rhabar, M., 1993. A Description and Assessment of RAMONA-3B Mod.0 Cycle 4: A Computer Code with Three-Dimensinal Neutron Kinetic for BWR System Transients. NUREG/CR-3664, BNL-NUREG-51746 (see also: RAMONA-3 Users’ Manual, Scandpower).