Adaptive matrix formation (AMF) method of space–time multigroup reactor kinetics equations in multidimensional model

annals of

NUCLEAR ENERGY Annals of Nuclear Energy 34 (2007) 103–119 www.elsevier.com/locate/anucene

Adaptive matrix formation (AMF) method of space–time multigroup reactor kinetics equations in multidimensional model A.E. Aboanber *, A.A. Nahla Department of Mathematics, Faculty of Science, Tanta University, Tanta 31527, Egypt Received 25 May 2006; accepted 24 July 2006 Available online 10 January 2007

Abstract Adaptive matrix formation (AMF) method has been developed for the numerical solution of the transient multigroup neutron diffusion and delayed precursor equations in two- and three-dimensional geometry. The method is applied to a general class of twoand three-dimensional problems. The results of numerical experiments, as well as comparison with space–time experimental results indicate that the method is accurate and that the two- and three-dimensional calculations can be performed at ‘‘reasonable’’ computer costs. Moreover, the AMF method offers the flexibility of using smaller time steps between flux shape calculations to achieve a specified accuracy and capability, without encountering numerical problems that occur in the other conventional methods. There is a large considerable saving in computer time and costs due to the partitioning of the matrix adopted in the presented AMF method. The two- and three-dimensional problems were analyzed with the present calculations model to illustrate the accuracy and stability of the method. Furthermore, the stability of the investigated method has been tested for sinusoidal, ramp, and step-change reactivity insertions. The results are in a good agreement with those of the other less approximate methods, including the problems in which the reflector zone is perturbed.  2006 Published by Elsevier Ltd.

1. Introduction Reactor dynamics calculations are concerned mainly with calculating the detailed spatial distribution of neutrons, the power production, the temperature, density and pressure distribution at any instant in time. Therefore careful consideration is required in order to develop methods in reactor dynamics which can predict the core behavior with sufficient accuracy and reliability. To be able to predict the dynamic behavior of physical systems, an appropriate mathematical model has to be developed. The model should be a true description of the system or the best estimate of the true description. However, the equations describing the model have to be of a form and level of complexity that can be solved with the *

Corresponding author. Tel.: +20 40 3300407; fax: +20 40 3350804. E-mail addresses: [email protected], aaboanber@yahoo. com (A.E. Aboanbaler), [email protected] (A.A. Nahla). 0306-4549/$ - see front matter  2006 Published by Elsevier Ltd. doi:10.1016/j.anucene.2006.07.012

available resources. These requirements contradict each other, and therefore, the success of any dynamics model in predicting the system behavior is strongly dependent on the assumption made in formulating the model and the problem to which the model or the code is applied. The motivation behind the development of methods for solving the energy-, space-, and time-dependent kinetic equations is not only the challenge of developing a method for solving a large set of coupled partial differential equations, but also a real need to predict the performance and assess the safety of large commercial reactors, both these presently operating and those being designed for the future. Since the early 1960s a host of new methods for spatio-temporal dynamics in two and three dimension were developed. For many of the basics methods a number of improvements were developed to improve solution speed and accuracy. In the following section, the methods for solving the multigroup space–time dependent diffusion equations are

104

A.E. Aboanber, A.A. Nahla / Annals of Nuclear Energy 34 (2007) 103–119

described, which in its one-, two-, or three-dimensional forms have been accepted by most analysts as the best approximation for the prediction of reactor dynamics behavior. The spatial methods are categorized as follows: direct methods, space–time factorization methods, and modal and synthesis methods. Each of the above methods is further subdivided into several subvariants. A large amount of literature and survey papers are available on each of these methods. Some of these methods, such as the modal and synthesis methods (Stacey, 1969), appear no longer to be in widespread use. In addition, the traditional nodal methods have largely been superseded by the more accurate transverse-integrated nodal (Lawrence, 1986; Chao and Shatilla, 1994; Azmy, 1992; DeLorey, 1993) and coarsemesh methods. The flux expansion method is a coarse-mesh method that has been implemented in the QUABOX and CUBBOX transient codes (Langenbuch et al., 1977a,b). However, the older nodal methods are still in use for applications such as simulators where execution speed is paramount. The fine-mesh finite difference method remains widely used when accuracy is important and a fast computer with substantial memory is available. In the finite difference method, the problem space is discretized by superimposing a computational mesh and requiring that material properties be treated as uniform within each mesh box (Adams, 1977; Nakamura, 1977; Diamond, 1996). The advent of advanced nodal homogenization and reconstruction procedures, however, has largely eliminated the limitations of the transverse-integrated nodal methods; and as a result these have replaced fine-mesh calculations for many applications. Coarse-mesh and finite element methods remain viable alternatives for certain types of applications. Another method employed for coarse mesh calculations is the finite element method (FEM) (Kang and Hansen, 1973). The FEM has been employed in several transient codes, including HERMITE (Rohan et al., 1975) and CRONOS (Kavenoky and Lautard, 1982; Lautard et al., 1991). Finally, the use of the Improved Quasistatic Method has been expanded to include its use in conjunction with modern nodal methods (Gehin, 1992; Taiwo and Khalil, 1992; Zimin and Ninokata, 1996; Chen et al., 1997). The methods for numerically advancing the space–time diffusion group diffusion equations, along with their timedependent delayed neutron precursor counterparts, through time have been developed. For each of these considered methods, the system to be solved has been reduced to the matrix form. The h-method is an accurate and efficient finite difference scheme that has been employed in the numerical integration of the time-dependent group diffusion equations since the mid-1960s (Henry and Vota, 1965). First-order finite difference schemes such as the implicit (backward) Euler, explicit (forward) Euler, and Crank–Nicholson (central difference) schemes can be derived easily from the h-method in its general form. A concise discussion of the h-method applied to the time-dependent group diffusion equations can be found in Henry and Vota (1965) and Stacey

(1969). The Alternating Direction Implicit (ADI) scheme is closely related to the h-method. However, the treatment of the resulting linear matrix equations is quite different, imposing an alternating method of obtaining a high degree of implicitness while reducing the computational effort requires inverting the matrix (see e.g., Hageman and Yasinsky (1969); Stacey (1969); Langenbuch et al., 1977a,b). A unique method of numerically integrating the space–time diffusion equations has been developed by Chao and Attard (1985). It is clear from previous paragraphs that the methods for solving space kinetics equations in two and three space dimensions are effective and accurate for few energy and delayed precursor groups, since the number of equations to be solved at each time step is generally quite small. However, a lock of definitive error bounds on final results. On the other hand, direct solution techniques include alternating direction implicit, explicit and Pade´ approximations for exponential of coefficient matrix are developed which solve the system of equations by finite difference approximations. These techniques are characterized by fairly definitive error bounds and are thus valuable as indirect techniques. Most of these methods are taking relatively large computational time requirements when compared with time step calculations. This paper is therefore concerned with presenting the equations that govern the reactor dynamics behavior in various approximations and proposing method and program to solve these equations, for large number of energy and delayed precursor groups, with an excellent accuracy and faster time of calculations using reasonable computer resources. Additionally, this method is characterized by fairly definitive error bounds when compared with other methods. The (AMF) method to be presented here is quite different from those previously discussed and is applicable to a very general class of two- and three-dimensional problems. The amount of computer time required is so small that their primary use will be as simulator for the actual situation in the reactor. Also, it can be used in the safety analysis of final design and thus compared well with the other conventional methods. The use of AMF computational mesh in space and time results in a dramatic decrease in the number of unknowns, and hence computer resources while the high order spatial treatment within a node, coupled with the use of advanced homogenization procedures, maintains the accuracy at the level of the more expensive fine-mesh calculations. In Section 2, the multigroup neutron diffusion kinetics equations with delayed neutrons are derived from neutron transport equation and formulated in their matrix form. A new method, Adaptable Matrix Formation (AMF), for solving this form of multigroup reactor kinetics equations in two- and three-dimensional has been developed. This method is characterized by the expansion of the coefficients matrix for all mesh points of the reactor into two different matrices, specified for each mesh point. Also, the method is

A.E. Aboanber, A.A. Nahla / Annals of Nuclear Energy 34 (2007) 103–119

presented and analyzed from a mathematical point of view; several modifications of the developed method are presented and discussed. Numerical solutions for several two- and three-dimensional sample problems have been obtained using the adopted new methods. In order to assess the validity and stability of the investigated AMF method the results are compared with other methods currently in use. These results are presented and compared in Section 3. Finally, the main conclusions are presented in Section 4. 2. The neutron diffusion kinetics equations 2.1. Derivation from transport equation The principal equation governing the neutron behavior in a nuclear chain reactor is the Boltzmann (Weinberg and Wigner, 1958) or transport equation: 1 oUðr; E; X; tÞ v ot ¼ X  rUðr; E; X; tÞ þ

Z

dE0 dX0 Uðr; E0 ; X0 ; tÞ

 Rs ðr; E0 ! E; X0 ! X; tÞ 

Z

dE0 dX0 Uðr; E; X; tÞ

 Rs ðr; E ! E0 ; X ! X0 ; tÞ þ Uðr; E; X; tÞ  Ra ðr; E; tÞ þ Sðr; E; X; tÞ

ð1Þ

Simplifying this equation by assuming that the differential scattering cross section Rs ðr; E ! E0 ; X ! X0 ; tÞ is independent of the incident angle and by integrating over all directions, the neutron scattering cross section Rs ðr; E; tÞ is obtained at r, E and t. The general source term S(r, E, t) is replaced by three separate and detailed source terms, the fission source, the delayed neutrons source, and the external, i.e. nonfission source, which is of importance when subcritical reactors are calculated. The total cross section Rt ðr; E; tÞ which represents all the processes to remove neutrons from dr to r and from dE to E at time t is defined. So, the adopted balance equation is   o Uðr; E; X; tÞ ot vðEÞ ¼ X  rUðr; E; X; tÞ  Rt ðr; E; tÞUðr; E; X; tÞ Z 1 Z 0 dX dE0 Rs ðr; E0 ! E; X0 ! X; tÞUðr; E0 ; X0 ; tÞ þ X0

Z

0

1 dX þ 0 4p X 

0

Z 0

1

dE0

N X

vn ðEÞð1  bn Þmn ðE0 Þ

n

Z

dX0

X0

Z

1

dE0 bnm mn ðE0 ÞRnf ðr; E; tÞUðr; E; X; tÞ

0

 knm C nm ðr; tÞ

ð3Þ

There is M · N equations for the delayed neutron precursors, where N is the number of fissile isotopes (n = 1, . . . , N) presented in the core and M is the number of delayed neutron groups (m = 1, . . . , M) used in the representation of delayed neutrons, where U(r, E, X, t) is the directional neutron flux at r in the interval dE about energy E, moving in the cone of directions dX about direction of X at time t, Rs ðr; E0 ! E; X0 ! XÞ is the total macroscopic scattering cross section for neutrons, Rt ðr; E; tÞ is the total cross section for neutron interaction with its environment, Rnf ðr; E; tÞ is the macroscopic fission cross section of the nth fissile isotope to neutrons having energy E at r and t, mn(E) is the average number of neutrons emitted by fission from isotope n at energy E, vn(E) is the spectrum of the neutron emitted in the fission process from isotope n, bnm is the fraction of fission neutrons from isotope n born as delayed neutrons into group m, knm is the decay constant of the appropriate delayed neutron precursor, v(E) is the speed of the neutron having energy E and S(r, E, X, t) is the neutron source. The diffusion equation is usually developed from the transport equation by expanding Eq. (2) into spherical harmonics and integrating over all directions X, neglecting several terms (Weinberg and Wigner, 1958; Glasstone and Edlund, 1952; Stacey, 1969). The diffusion equation is rearranged in a form similar to that of the neutron transport equation as 1 oUðr; E; tÞ vðEÞ ot ¼ r  Dðr; E; tÞrUðr; E; tÞ  Rt ðr; E; tÞUðr; E; tÞ Z þ dE0 Rs ðr; E0 ! E; tÞUðr; E0 ; tÞ þ

N X

n

n

v ðEÞð1  b Þ

Z

dE0 mn ðE0 ÞRnf ðr; E0 ; tÞUðr; E0 ; tÞ

n

þ

N X M X n

vnm ðEÞknm C nm ðr; tÞ þ Sðr; E; tÞ

ð4Þ

m

where D(r, E, t) is the diffusion coefficient. For the delayed neutron precursors C nm ðr; tÞ an equation similar to Eq. (3) can be written: oC nm ðr; tÞ ¼ bnm ot

Rnf ðr; E0 ; tÞUðr; E0 ; X0 ; tÞ

N X M 1 X þ vn ðEÞknm C nm ðr; tÞ þ Sðr; E; X; tÞ 4p n m m

oC nm ðr; tÞ ¼ ot

105

Z

1

dE0 mðE0 ÞRnf ðr; E0 ; tÞUðr; E0 ; tÞ  knm C nm ðr; tÞ 0

ð5Þ ð2Þ

The equations for the delayed neutron precursors are given by

The notations were the same as in Eqs. (2) and (3) except that the scalar flux U(r, E, t) is used instead of the directional flux U(r, E, X, t).

106

A.E. Aboanber, A.A. Nahla / Annals of Nuclear Energy 34 (2007) 103–119

Eqs. (4) and (5) usually in its multigroup formulation are the basis for the vast majority of computer programs developed for solving both the reactor steady state and time-dependent equations. By integrating Eqs. (4) and (5) in limited energy intervals from Eg to Eg+1 obtained by subdividing the energy, E, from zero to infinity into G intervals, a set of G coupled differential equations is obtained: 1 owg ðr; tÞ ¼ r  Dg ðr; tÞrwg ðr; tÞ  Rr;g ðr; tÞwg ðr; tÞ vg ot g1 X Rs;g0 !g ðr; tÞwg0 ðr; tÞ þ g0 ¼1

þ ð1  bÞvg

G X

mg0 Rf ;g0 ðr; tÞwg0 ðr; tÞ

g0 ¼1

þ

G X

vm;g km ðrÞC m ðr; tÞ þ S g ðr; tÞ

ð6Þ

tion. Obviously, the most precise representation would be in the three dimensions. This, however, can be extremely time consuming; also, reevaluation of code due to the change of the cross sections in the case of feedback is not practical for several methods which have been proposed by different code developers. To overcome these problems, in the present work the system of multigroup diffusion equations is formulated in the matrix form without approximation which results in a high degree of accuracy, provides significant time saving and their applicability is unlimited. The nuclear reactor kinetics equations (6) and (7) with delayed neutrons can be rewritten in the following general form: G X o Ug ðr; tÞ ¼ vg r  Dg ðr; tÞrUg ðr; tÞ þ T gg0 ðr; tÞUg0 ðr; tÞ ot g0 ¼1

m¼1

þ

and for the delayed neutron precursor:

Id X

F gi C i ðr; tÞ;

g ¼ 1; 2; . . . ; G

ð8Þ

i¼1 G X o C i ðr; tÞ ¼ P ig Ug ðr; tÞ  ki C i ðr; tÞ; ot g¼1

G X oC m ðr; tÞ ¼ bm ðrÞ mg0 Rf ;g0 ðr; tÞwg0 ðr; tÞ  km ðrÞC m ðr; tÞ ot g0 ¼1

i ¼ 1; 2; . . . ; I d ð9Þ

ð7Þ Eqs. (6) and (7) are now written for only one fissionable isotope (N = 1), as opposed to Eqs. (4) and (5) that are written for N fissionable isotopes. The process of obtaining the group energy averaged cross sections for some specific materials composition and geometry is described in many textbooks. 2.2. Adaptive matrix formation (AMF) method The multigroup diffusion equations are used to describe the reactor spatio-temporal behavior, there are some diffi-

r  Dg rUg ffi

1 2

where T g;g ¼ vg vg mg Rfg ð1  bÞ  Rag 

X

! Rsgg0

g0 >g

T g;g0 ¼ vg ðvg mg0 Rfg0 ð1  bÞ þ Rsg0 g Þ; Rsg0 g ¼ 0 F gi ¼ vg fgi ki

and

P ig ¼ bi mg Rfg ;

for g0 > g

where f gi ¼ vgi

The leakage term $ Æ Dg$Ug for three-dimensional geometry is approximated using the seven-point central finite difference approximation, Varga (1962)

h   i Dg;j12;k;l Ug;j1;k;l  Dg;j12;k;l þ Dg;jþ12;k;l Ug;j;k;l þ Dg;jþ12;k;l Ug;jþ1;k;l

ðDxÞ   i 1 h þ Dg;j;k12;l Ug;j;k1;l  Dg;j;k12;l þ Dg;j;kþ12;l Ug;j;k;l þ Dg;j;kþ12;l Ug;j;kþ1;l 2 ðDyÞ   i 1 h þ D D 1 Ug;j;k;l1  1 þ Dg;j;k;lþ1 Ug;j;k;l þ Dg;j;k;lþ1 Ug;j;k;lþ1 g;j;k;l g;j;k;l 2 2 2 2 2 ðDzÞ

culties involved in the use of these equations. The first problem is the choice of the number of energy intervals and their boundaries. It is obvious that the larger the number of energy groups, the more precise the solution. However, the use of many energy groups for time-dependent problems will increase significantly the computer time necessary for a particular solution. The second problem is the choice of the number of dimensions to be used in the solu-

ð10Þ

Rewrite Eqs. (8) and (9) in a matrix form as d W ¼ AW þ B dt

ð11Þ

where W ¼ col½ U1 U2    UG C 1 C 2    C I d  B ¼ col½ W 01 W 02    W 0G 0 0    0 

ð12Þ ð13Þ

A.E. Aboanber, A.A. Nahla / Annals of Nuclear Energy 34 (2007) 103–119

0

W 1 þ T 11 B T 21 B B .. B B . B B T G1 B A¼B B P 11 B B P B 21 B .. B @ . P Id 1

T 12 W 2 þ T 22 .. .

  .. .

T 1G T 2G .. .

F 11 F 21 .. .

F 12 F 22 .. .

  .. .

F 1I d F 2I d .. .

F G1 k1

F G2 0

 

F GI d 0

 .. .

0 .. .

T G2 P 11

   W G þ T GG  P 1G

P 22 .. .

 .. .

P 2G .. .

0 .. .

k2 .. .

P Id 2



P Id G

0

0

 D 1;k;l Ug;j1;k;l þ Dg;jþ1;k;l Ug;jþ1;k;l g;j 2 2 ðDxÞ  2  vg þ D 1;l Ug;j;k1;l þ Dg;j;kþ1;l Ug;j;kþ1;l g;j;k 2 2 ðDyÞ2   vg þ D 1 Ug;j;k;l1 þ Dg;j;k;lþ1 Ug;j;k;lþ1 ; g;j;k;l 2 2 2 ðDzÞ g ¼ 1; 2; . . . ; G  vg  Wg ¼ Dg;j12;k;l þ Dg;jþ12;k;l 2 ðDxÞ  vg  D þ D  1 1 g;j;k2;l g;j;kþ2;l 2 ðDyÞ   vg Dg;j;k;l12 þ Dg;j;k;lþ12 ; g ¼ 1; 2; . . . ; G  2 ðDzÞ The general solution of Eq. (11) will be

W 0g ¼

vg 

WðtÞ ¼ expðtAÞWð0Þ þ ðexpðtAÞ  IÞA1 B

ð15Þ

ð16Þ

ð17Þ

ð18Þ

This equation suggests a form that should be quite suitable for generating an approximate solution of the more general problem (time-dependent problem). It is well known that the exponential of A (eA) is given by the series: 1 X Ak expðAÞ ¼ k! k¼0

1 C C C C C C C C C C C C C C C A

ð14Þ

   kI d

3. Numerical results and discussion

If Wn and Wn+1 denote the solution at times tn and t(n+1) = tn + Dt, respectively, then Wnþ1 ¼ expðDtAÞWn þ ðexpðDtAÞ  IÞA1 B

107

In order to assess the efficiency of the AMF method we developed the FORTRAN computer code that solves the two-dimensional multigroup problems. At present, the code can handle uniform mesh spacing in rectangular geometry. Results for several multigroup trail problems in two- and three dimensions obtained with the AMF code have been compared with the results available in literature. These results were obtained with computer reference codes such as LUMAC (Reed and Hansen, 1970; Wight et al., 1971); MITKIN (Reed and Hansen, 1970); TWIGL (Hageman and Yasinsky, 1969); SADI (Wight et al., 1971) and TSM (Utku and Christenson, 1994). The solutions of the trial problems are presented, discussed and compared in the following part of the work. All calculations presented in the following sections are performed on a personal computer. The exact data for each problem is listed in Appendix. It should be noted that the truncation error discussed herein is the difference between the solution under consideration and the exact solution of different methods. In the case of the homogeneous problem to be presented as a test case, the exact solution can be generated using an eigenvector expansion technique. 3.1. Two-dimension reactor calculations

ð19Þ

so, the solution of Eq. (18) can be obtained, in principle, from the convergent series: " # ðDtAÞ2 ðDtAÞn expðDtAÞ ¼ I þ ðDtAÞ þ þ  þ þ  2! n! ð20Þ A general expression for a special type of functions has been introduced (Aboanber and Nahla, 2006). Such expression permits us to approximate the exponential of a matrix in an economical manner using different types of Pade´ approximations (Pade´11 and Pade´12). The stability of these approximations is also included in reference (Aboanber and Nahla, 2006).

3.1.1. Case 1: Two dimension homogeneous reactor Perturbation: Step change; DRa ¼ 0:369  104 The test case is a two-region, bare, homogeneous system with one precursor group. The initial condition for this case was a critical distribution consistent with the data presented in Appendix. Ten mesh intervals were used in each direction, and the boundary conditions were homogeneous Dirichlet on all sides. At time t = 0, a step change in group two consists of a uniform step decrease in the thermal group capture cross section and had a reactivity worth of about 50 cents. Since the problem is symmetric about the mid-plane in the x-direction, only the right half of the reactor was actually calculated using the two-dimensional AMF program.

108

A.E. Aboanber, A.A. Nahla / Annals of Nuclear Energy 34 (2007) 103–119

Table 1 Thermal flux at the center of the perturbed region for two-dimensional homogenous reactor Time (s)

Pade´ (Dt = 0.004)

MITKIN (Dt = 0.004)

AMF (Dt = 0.0004)

Exact value

0.0 0.08 0.16 0.24 0.32 0.40

0.382 0.609 0.805 0.974 1.122 1.252

0.382 0.526 0.773 1.008 1.205 1.370

0.382 0.609 0.805 0.974 1.122 1.253

0.382 0.613 0.816 0.997 1.158 1.303

The results of Table 1 showed that the rate of convergence of the method to the exact solution is less than the rate of convergence of MITKIN method. The relative percentage error at a time-step size Dt = 0.0004 s is less than 4% which seems perfectly adequate for most purposes. To illustrate the utility of the method for truly space-dependent problems, comparison was made with results obtained from other reference computer codes, Table 1. 3.1.2. Case 2: Two-dimension heterogeneous symmetric reactor This test case is one with spatial dependence on the solutions. Analytic solutions are not available, but approximate results are available from the TWIGL codes. The geometry is given in Fig. 1 and the problem data are listed in Appendix. The reactor consists of square core surrounded by a blanket with blanket in the interior as well. It is completely symmetric. The reactor quadrant modeled in this calculation consists of three regions containing a total of 100 fuel assemblies, each of a width of 8 cm. The transient is induced by changing the thermal cross section in region 1, three different transients were studied: a positive step change in reactivity, a positive ramp change in reactivity and changing the thermal absorption cross section in a sinusoidal form. Solutions were obtained from

TWIGL, LUMAC, MITKIN, SADI, TSM and Pade´ and Cut-product (Aboanber and Nahla, 2006) method and compared with the presented AMF method. 3.1.2.1. A positive step reactivity change. A positive step change in reactivity of about 50 cents is introduced at time zero by reducing the thermal capture cross section in region 1. The perturbation value is DRa = 0.0035 cm1. The thermal flux shape at the center of the reactor is used as basis of comparison because the problem is space dependent and the change in flux shape is quite small. The results of comparison between the developed AMF method and the other conventional methods are listed in Table 2. All the developed AMF, Pade´ and Cut product approximations methods agree quite closely on the basis of their accuracy. On the other side, on the basis of CPU time the superiority is recorded for AMF over the entire other conventional methods. These results of comparison with the other reference computer codes show an excellent advantage for the tested AMF method. Solution to 0.4 s required 0.0011 s per time step of Pentium IV, 2.6 MHz Computer time. 3.1.2.2. A positive ramp reactivity change. A positive ramp change in reactivity of about 2.5 $/s is introduced over the time interval 0 6 t 6 0.2 s by reducing the thermal capture cross section in the region 1 of the reactor. The results are compared in Tables 3 and 4. The agreement with reference computer codes is excellent. The execution time for the AMF code of Pentium IV, 2.6 MHz computing system was 0.0011 s, Tables 3 and 4, which is very small compared with different Pade´ and Cut product approximations and TSM method. 3.1.2.3. A positive sinusoidal ramp reactivity change. The more interesting case of an oscillatory perturbation for this

y 3 20 19 18 17 2 1 1 16 15 14 13 12 3 2 2 3 3 11 10 9 8 7 2 1 1 6 5 4 3 3 2 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 x Fig. 1. Two-dimensional geometry for TWIGL problem (Dx = 8.0 cm, Dy = 8.0): 1 – driver region; 2 – core region; 3 – blanket region.

A.E. Aboanber, A.A. Nahla / Annals of Nuclear Energy 34 (2007) 103–119

109

Table 2 Response of the thermal flux at the center of two-dimensional TWIGL reactor, point (80, 80), for a positive step reactivity Method

Dt

Time (s) 0.00

0.01

0.02

0.03

0.04

0.05

0.1

0.2

AMF Method

0.0001

16.75 (0.0%)

27.17 (1.76%)

31.35 (1.85%)

32.91 (1.57%)

33.51 (1.09%)

33.75 (0.62%)

34.02 (0.03%)

34.31 (0.0%)

Pade´11

0.001

16.75 (0.0%)

27.17 (1.76%)

31.32 (1.75%)

32.89 (1.51%)

33.49 (1.03%)

33.74 (0.60%)

34.01 (0.0%)

34.30 (0.03%)

Pade´112

0.001

16.75 (0.0%)

27.21 (1.91%)

31.35 (1.85%)

32.90 (1.54%)

33.50 (1.06%)

33.74 (0.60%)

34.01 (0.0%)

34.30 (0.03%)

Cutp22

0.001

16.75 (0.0%)

27.16 (1.72%)

31.31 (1.72%)

32.88 (1.48%)

33.48 (1.0%)

33.73 (0.57%)

34.00 (0.03%)

34.29 (0.06%)

Cutp33

0.001

16.75 (0.0%)

27.00 (1.12%)

31.04 (0.85%)

32.56 (0.50%)

33.14 (0.02%)

33.38 (0.48%)

33.64 (1.09%)

33.92 (1.13%)

33.15

33.54

34.01

34.31

34.03 (0.06%) 33.64 (1.09%)

34.33 (0.06%) 33.93 (1.11%)

TWIGL

0.001

16.75

26.70

30.78

32.40

LUMAC

0.00008

16.75 (0.0%)

27.29 (2.21%)

31.48 (2.27%)

33.06 (2.04%)

MITKIN

0.0002

16.75 (0.0%)

27.32 (2.32%)

31.50 (2.34%)

33.13 (2.25%)

33.97 (2.47%)

34.63 (3.25%)

SADI

0.00025

27.15 (1.69%)

TSM

0.0025

16.75 (0.0%) 16.75 (0.0%)

33.21 (7.89%) 31.93 (3.74%)

33.94 (4.75%) 33.00 (1.85%)

33.90 (2.26%) 33.34 (0.57%)

33.88 (1.01%) 33.46 (0.24%)

Pade´112 refers to different type of Pade´ approximations (Aboanber and Nahla, 2006).

Table 3 The thermal flux at the center of two-dimensional TWIGL reactor, point (80, 80), for a ramp reactivity insertion Method

Dt

CPU (s)

0.0011

Time (s) 0.00

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

16.75 (0.0%)

18.76 (0.0%)

21.73 (0.05%)

25.94 (0.09%)

32.31 (0.19%)

34.10 (0.15%)

34.25 (0.03%)

34.40 (0.02%)

34.54 (0.0%)

AMF method

0.00025

Pade´11

0.001

55.14

16.75 (0.0%)

18.77 (0.05%)

21.75 (0.05%)

25.97 (0.05%)

32.37 (0.0%)

34.08 (0.09%)

34.24 (0.0%)

34.38 (0.03%)

34.53 (0.03%)

Pade´112

0.001

79.89

16.75 (0.0%)

18.78 (0.11%)

21.76 (0.09%)

25.98 (0.08%)

32.38 (0.03%)

34.08 (0.09%)

34.24 (0.0%)

34.38 (0.03%)

34.53 (0.03%)

Cutp22

0.001

96.94

16.75 (0.0%)

18.77 (0.05%)

21.75 (0.05%)

25.97 (0.05%)

32.36 (0.03%)

34.07 (0.06%)

34.23 (0.03%)

34.37 (0.06%)

34.52 (0.06%)

Cutp33

0.001

113.4

16.75 (0.0%)

18.70 (0.30%)

21.93 (0.88%)

26.42 (1.77%)

32.44 (0.22%)

33.72 (0.98%)

33.87 (1.1%)

34.01 (1.1%)

34.15 (1.1%)

TWIGL

0.01

–

16.75

18.76

21.74

25.96

32.37

34.05

34.24

34.39

34.54

LUMAC

0.00008

–

16.75 (0.0%)

MITKIN

0.0001

–

16.75 (0.0%)

18.79 (0.16%)

21.75 (0.05%)

25.95 (0.02%)

32.31 (0.16%)

34.12 (0.22%)

33.57 (1.9%)

33.23 (3.37%)

33.19 (3.9%)

SADI

0.00025

–

16.75 (0.0%)

18.70 (0.33%)

21.66 (0.37%)

25.84 (0.47%)

32.15 (0.66%)

34.58 (1.58%)

34.29 (0.14%)

34.44 (0.15%)

34.60 (0.2%)

MADI

0.001

–

18.75 (0.05%)

18.77 (0.05%)

21.79 (0.23%)

TSM

0.01

106.3a

16.75 (0.0%)

18.75 (0.05%)

21.71 (0.14%)

21.73 (0.03%)

32.49 (0.38%)

34.83 (1.72%)

32.43 (0.19%) 25.89 (0.27%)

32.22 (0.46%)

34.30 (0.18%)

CPU (s) = Execution time per time step, all calculations were done under the same conditions. Pade´112 refers to different type of Pade´ approximations (Aboanber and Nahla, 2006). a The computation was carried out on a VAX 3100 computer system.

110

A.E. Aboanber, A.A. Nahla / Annals of Nuclear Energy 34 (2007) 103–119

Table 4 Response of the thermal flux at a point (40, 40), two-dimensional TWIGL reactor, for the ramp reactivity insertion Method

Dt

AMF

0.00025

Pade´11

0.001

Pade´112

CPU (s)

Time (s) 0.00

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

5.39 (%)

6.15 (0.0%)

7.25 (0.0%)

8.81 (0.11%)

11.18 (0.18%)

11.79 (0.17%)

11.84 (0.0%)

11.89 (0.0%)

11.94 (0.0%)

55.14

5.39 (0.0%)

6.15 (0.0%)

7.26 (0.14%)

8.82 (0.0%)

11.20 (0.0%)

11.78 (0.08%)1

11.84 (0.0%)

11.89 (0.0%)

11.94 (0.0%)

0.001

79.89

5.39 (0.0%)

6.15 (0.0%)

7.26 (0.14%)

8.83 (0.11%)

11.21 (0.09%)11

1.79 (0.17%)1

11.84 (0.0%)

11.89 (0.0%)

11.94 (0.0%)

Cutp22

0.001

96.94

5.39 (0.0%)

6.15 (0.0%)

7.25 (0.0%)

8.82 (0.0%)

0.20 (0.0%)

1.78 (0.08%)1

11.83 (0.08%)1

11.89 (0.0%)

11.94 (0.0%)

Cutp33

0.001

113.4

5.39 (0.0%)

6.13 (0.32%)

7.31 (0.82%)

9.00 (2.0%)

11.24 (0.35%)

1.68 (0.76%)

1.73 (0.92%)

11.78 (0.92%)

11.83 (0.92%)

TWIGL ADI

0.01 0.0005

– –

5.39 5.39 (0.0%)

6.15 6.16 (0.16%)

7.25 7.25 (0.0%)

8.82 8.82 (0.0%)

11.20 11.18 (0.18%)

11.77 11.80 (0.25%)

11.84

11.89

11.94

MADI

0.001

–

5.39 (0.05%)

6.12 (0.49%)

7.20 (0.69%)

TSM

0.01

106.3

5.39 (0.0%)

6.16 (0.16%)

7.25 (0.0%)

0.0011

11.01 (1.70%) 8.80 (0.23%)

11.15 (0.45%)

11.71 (1.1%)

Pade´112 refers to different type of Pade´ approximations (Aboanber and Nahla, 2006).

reactor configuration has been also investigated. The numerical experiments results for this case of sinusoidal change were performed according to the form:    2pt Ra2 ðtÞ ¼ Ra2 ð0Þ 1  0:2 sin T or,

   2pt Ra2 ðtÞ ¼ Ra2 ð0Þ 1 þ 0:2 sin ; T

T ¼ 102

The behavior of the fast and thermal neutron flux shapes obtained by AMF method is studied and the cross section perturbation is introduced in the four cross-hatched seed regions in an oscillatory (sinusoidal) fashion. Fig. 2, illustrates the time history of the thermal neutron flux over the first full period for the case where the amplification factor of the sinusoidal term was chosen to be large enough to create a complete distribution of the flux shape and thus produce a stringent test of the ability of the AMF to handle situations in which there are simultaneously large change in both shape and amplitude of the group fluxes. The first peak of the thermal flux for the positive change of reactivity appears in the blanket region at a short distance from the blanket core interface; while the maximum peak arises at the center of the core region followed by a shorter peak at the core–blanket interface. Also, the time history of the thermal flux at the center of the perturbed region over three cycles appears in the same figure. 3.1.3. Case 3: OBLONG – nonhomogeneous, nonsymmetrical reactor The TWIGL problem is completely symmetric and shows very little change in flux shape over transient

(<5%). A severe test of a space–time method requires a sample problem with no symmetries whatsoever, and with a significant change in flux shape. The OBLONG problem was designed for that purpose. The reactor is a four-energy-group, one-delayed-group system. This is an OBLONG problem with 160 · 80 cm. As shown in Fig. 3, there are two fuel materials and a reflector zone. Region 1 is driver region, region 2 is a blanket region, and regions 3 and 4 are a water reflector. The cross sections of the materials are quite different and there is a drastic variation of the flux in all groups. The perturbation is caused by changing the capture cross section in group 4 of a part of the reflector region, region 4, according to the following manner: DRc ðtÞ ¼ 

0:003 t cm1 0:2

0 6 t 6 0:2

and D Rc ðtÞ ¼ 0:003 cm1

t > 0:2

Results are shown in Table 5 four group 1 and group 4 at space points (88, 16) cm and (16, 64) cm, respectively. The tabulated data illustrate the magnitude of the shape and spectrum change over the transient. Furthermore, a sinusoidal change in the group 1 and 4 are introduced in Fig. 4, where the thermal absorption cross section of the thermal group is perturbed according to the shape mentioned in the previous section. These results showed that rather large changes occurred in the spatial shape and spectrum over the transient and verified the data obtained in Table 5. The proximity of the solutions at 0.0001 s time step sizes recorded for the proposed method implies that

A.E. Aboanber, A.A. Nahla / Annals of Nuclear Energy 34 (2007) 103–119

450

111

Point(80,80)

400 350 Thermal Flux

300 Point(40,40)

250 200 150 100 50 0 0

0.01 0.02 Time (Second)

0.03

Fig. 2. Thermal flux in two-dimensional for sinusoidal change for TWIGL.

y 10 9 1

8

2

7 6 5

Unperturbed Region

4 2

3

m Perturbed Region

2 1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

3 17 18

( m = 3, 0 ≤ z ≤ 56 cm; m = 4 , 56 ≤ z ≤ 120 cm ) Fig. 3. Geometry for OBLONG reactor problem.

19

20

x

112

A.E. Aboanber, A.A. Nahla / Annals of Nuclear Energy 34 (2007) 103–119

Table 5 Flux of group 1 and group 4 at different perturbed points for two-dimensional OBLONG reactor, for the ramp reactivity insertion Methods

Time (s)

Point (88, 16)

Point (16, 64)

Group 1

Group 4

Group 1

Group 4

LUMAC (EP1 = 0.00008) MITKIN (Dt = 0.0005) AMF (Dt = 0.0001)

0.0

0.1314 0.1314 0.1341

0.9684 0.9684 0.9684

0.4463 0.4463 0.4463

0.03594 0.03594 0.03594

LUMAC (EP1 = 0.00008) MITKIN (Dt = 0.0005) AMF (Dt = 0.0001)

0.05

0.1385 0.1385 0.1382

1.056 1.054 1.052

0.4569 0.4566 0.4555

0.03680 0.03677 0.03667

LUMAC (EP1 = 0.00008) MITKIN (Dt = 0.0005) AMF (Dt = 0.0001)

0.10

0.1453 0.1433 0.1429

1.166 1.153 1.150

0.4730 0.4679 0.4664

0.03810 0.03767 0.03755

LUMAC (EP1 = 0.00008) MITKIN (Dt = 0.0005) AMF (Dt = 0.0001)

0.15

0.1499 0.1488 0.1484

1.278 1.270 1.266

0.4830 0.4807 0.4791

0.03890 0.03869 0.03857

LUMAC (EP1 = 0.00008) MITKIN (Dt = 0.0005) AMF (Dt = 0.0001)

0.20

0.1551 0.1553 0.1548

1.410 1.408 1.404

0.4943 0.4955 0.4941

0.03980 0.03988 0.03976

LUMAC (EP1 = 0.00008) MITKIN (Dt = 0.0005) AMF (Dt = 0.0001)

0.30

0.1605 0.1553 0.1587

1.451 1.409 1.400

0.5123 0.4957 0.4947

0.04120 0.03989 0.03980

Fig. 4. Flux in two-dimensional for sinusoidal change for OBLONG.

the results are quite accurate, and the close agreement with the results obtained with the LUMAC code is a strong evidence that a good approximation to the exact solution has been obtained. Regarding that the MITKIN results are among the most accurate codes available. The obtained AMF results showed a considerable decrease of five folds in the time step than that obtained using the MITKIN code, Table 5. However, the recorded maximum error is less than 1% and the generation of the above solutions to 0.3 s required 0.0016 s execution time per time step of computer time. These computational results showed that this

method can reach the same accuracy as the standard reference methods in sensationally shorter computational times. 3.2. Three-dimensional reactor calculations The above three two-dimensional problems were extended to three-dimensional geometry and analyzed by AMF code as following. 3.2.1. Case 1: Three-dimensional homogeneous reactor Perturbations: Step change, DRa2 ¼ 0:369  104

A.E. Aboanber, A.A. Nahla / Annals of Nuclear Energy 34 (2007) 103–119

3.2.2. Case 2: Three-dimensional heterogeneous symmetric reactor Geometry: Fig. 5, Dx = Dy = Dz = 8.0 cm, height = 160 cm. Perturbation: Ramp change

Table 6 Thermal flux at the center of the perturbed region for three-dimensional homogenous reactor Time (s)

3DKIN (Dt = 0.001)

AMF (Dt = 0.0001)

Exact value

0.0 0.05 0.10 0.15 0.20 0.30 0.40

0.816 1.124 1.406 1.660 1.890 2.289 2.622

0.816 1.127 1.407 1.660 1.890 2.287 2.620

0.816 1.127 1.407 1.660 1.890 2.288 2.620

0:0045 t for 0 6 t 6 0:2 s 0:2 DRa ðmaterial 1; group 2Þ ¼ 0:0045 for t > 0:2 s DRa ðmaterial 1; group 2Þ ¼ 

As mentioned in the previous section, the original twodimensional version of this problem has been used to test several two-dimensional solution methods and codes. This is a heterogeneous problem with octant radial symmetry. The overall reactor geometry is cubic with a side length of 160 cm. An axial blanket material of 24 cm in thickness is presented at the top and the bottom. Four regions containing material 1 are perturbed. They are located symmetrically along the diagonal of the four quadrants of the core. In the three dimensions geometry, this configuration was made of 112 cm in thickness to the z-direction, and a blanket of 24 cm in thickness was added to the top and to the bottom. Thus, the overall reactor is cubic, 160 cm on a side and it has two neutron groups and one delayed precursor group. The four perturbed regions containing material 1, each of 32 · 32 · 112 cm in size, are located symmetrically with respect to the central z-plane. The perturbation is a linear ramp with the above cross section changed effected in 0.2 s and the perturbed cross section is maintained thereafter. Due to the symmetry, only the right half of the cube was considered with a homogenous Neumann boundary condition imposed at the exposed mid-plane. The initial flux distribution and eigenvalue were computed with the steady-state option of AMF test code.

Geometry: This case is a bare homogenous cube of side length 200 cm on each side with two neutron groups and one precursor group. Ten mesh intervals were used in each direction and the boundary conditions were homogenous Dirichlet on all six sides. The results of AMF runs using two different time-step size at various times in the transient are compared with the 3DKIN (Ferguson and Hansen, 1973; Bucknel and Stewart, 1976) in addition to the exact method, Table 6. The tabulated values of the AMF code presented the thermal group fluxes at the center point of the reactor and were found to be in a good agreement with the results of the other conventional codes. Even when the AMF time step was small, the CPU computer time required to run the code was very small. The next three investigated cases are all spatially dependent problems. Test cases 2 and 3 are three-dimensional versions of problems already used to test some or all of the methods discussed in the previous section of the twodimensional geometries.

x cm

21

3

2

8

3

3

3

3

3

3

14

D∇φ = 0

y mesh points

D∇φ = 0

y mesh points

14

3

3

18

1

2

x cm

21

3

18

113

3

8

1

2 4

4

3

1 1

5

3

1 8

11

x-mesh points

1

5

8

11

x-mesh points

Fig. 5. Three-dimensional geometry for TWIGL reactor problem (a) x-y plane for 24 6 z 6 136 cm, (b) x-y plane for 0 6 z 6 24 cm and 136 6 z 6 160 cm.

114

A.E. Aboanber, A.A. Nahla / Annals of Nuclear Energy 34 (2007) 103–119

The calculations of this test case were carried out up to 0.3 s with steps of 0.001 s and 0.0001 s and the result of the AMF runs for these two time-step sizes are presented in Table 7. The tabulated values are corresponding to the thermal flux values and the recorded average execution time per time step is 0.0121 s. Comparison with the results of the ADI (3DKIN) code is quite good for this case with

a recorded deviations (relative percentage error) of the order of approximately 1% are observed for AMF code, Table 7. This could have been caused by the oscillations in the 3DKIN flux values, which are due to the transformation of fluxes used there in. The agreement with 3DKIN appears to be very good at different times during the transient. The z-planes 16 and 144 are at 16 cm below

Table 7 Thermal flux at different points for TWIGL reactor in the three-dimensional Methods

Time

Point

(s)

(0, 80, 16)

(40, 40, 16)

(0, 80, 80)

(40, 40, 80)

(0, 80, 144)

(40, 40, 144)

3DKIN (Dt = 0.001) AMF (Dt = 0.001) AMF (Dt = 0.0001)

0.0

0.347 0.347 0.347

0.245 0.244 0.244

1.279 1.279 1.279

0.422 0.421 0.421

0.347 0.347 0.347

0.245 0.244 0.244

3DKIN (Dt = 0.001) AMF (Dt = 0.001) AMF (Dt = 0.0001)

0.05

0.398 0.393 0.400

0.284 0.280 0.284

1.467 1.473 1.478

0.496 0.497 0.498

0.395 0.405 0.401

0.283 0.289 0.285

3DKIN (Dt = 0.001) AMF (Dt = 0.001) AMF (Dt = 0.0001)

0.10

0.483 0.472 0.487

0.349 0.340 0.350

1.780 1.774 1.800

0.616 0.612 0.621

0.479 0.489 0.488

0.348 0.352 0.352

3DKIN (Dt = 0.001) AMF (Dt = 0.001) AMF (Dt = 0.0001)

0.15

0.619 0.602 0.626

0.454 0.440 0.457

2.284 2.274 2.319

0.809 0.804 0.820

0.615 0.629 0.629

0.452 0.459 0.459

3DKIN (Dt = 0.001) AMF (Dt = 0.001) AMF (Dt = 0.0001)

0.20

0.867 0.827 0.882

0.643 0.611 0.651

3.197 3.149 3.269

1.162 1.142 1.184

0.860 0.878 0.887

0.640 0.650 0.656

3DKIN (Dt = 0.001) AMF (Dt = 0.001) AMF (Dt = 0.0001)

0.25

0.994 1.046 0.997

0.737 0.773 0.736

3.666 3.862 3.686

1.330 1.397 1.333

0.987 1.042 0.998

0.734 0.769 0.737

3DKIN (Dt = 0.001) AMF (Dt = 0.001) AMF (Dt = 0.0001)

0.30

0.991 0.999 1.008

0.735 0.737 0.745

3.655 3.689 3.727

1.326 1.334 1.348

0.984 0.998 1.009

0.732 0.737 0.745

Fig. 6. Thermal flux in three-dimensional for sinusoidal change for TWIGL.

A.E. Aboanber, A.A. Nahla / Annals of Nuclear Energy 34 (2007) 103–119

and above the core respectively, while z-plane 80 is the central z-plane. Points (40, 40, z) are on the central axis of the one of perturbed regions (region 1). Also, the points (0, 80, z) are on the central of z-axis. The variation of the thermal flux in three directions for sinusoidal change at different z levels and times t = (T/2) is shown in Fig. 6. There is a gradual increase in the thermal flux at the reactor center z = 80 cm and in the center of the perturbed region. On the other hand, there is a gradual decrease in the thermal flux at the reflector region; the short peak at the interface in this reflector region is due to the thermalization process. 3.2.3. Case 3: Three-dimensional heterogeneous nonsymmetric reactor This problem, with four neutron groups and one precursor group, is a three-dimensional version of a problem used to compare several methods in the two dimension calculations. The dimension is 160 · 80 · 120 cm. The problem is asymmetric in all three directions, so a full reactor with homogenous Dirichlet boundary conditions had to be considered. Geometry: Fig. 3, Dx = Dy = Dz = 8.0 cm, height = 120 cm. Perturbation: Ramp change, DRa ¼ ðmaterial 4; group 4Þ ¼ 0:0035  ðt=0:2Þ for 0 6 t 6 0:2 s DRa ¼ ðmaterial 4; group 4Þ ¼ 0:0035

for t P 0:2 s

115

>Material 1 and 2 are fuel. Material 1 is a highly enriched material so that group 1 fluxes are initially more than five times higher than the group 4 fluxes in it. On the other hand, material 3 and its perturbed zone 4 are strong moderators so that group 4 flux are peaks in them. The perturbation is a linear ramp in the absorption cross section of material zone 4 with dRa being 0.0035 in 0.2 s in the axial direction zone 4 is restricted to only 64 cm and hence the bottom 56 cm is not perturbed. Given these spectral variations in the initial condition, which is computed with AMF, and the unsymmetrical perturbation, large spatial and energy spectrum changes are results as it is expected. The average execution time per time step is 0.042 s of computer time. The results of the AMF runs up to 0.3 s with time-step size of 0.001 s and 0.0001 s are shown in Table 8. Points (16, 64, z) are near the center line of the material zone 1; z-plane at 32 cm height is the mid-plane of the unperturbed lower portion, while z-plane at 88 cm high is near the center of the upper 64-cm region. As expected, this transient resulted in rather severe spectral changes. At point (16, 64, 32) the group 1 and group 4 fluxes grew by only @6%. Meanwhile, the group 1 and group 4 fluxes at point (88, 16, 88) grew by 11% and 45%, respectively. Fig. 7 shows the three-dimensional thermal and fast fluxes for sinusoidal change perturbation at different planes in the z-direction and different energy groups and times. Data showed that the results of AMF code at the two time-step sizes are in good agreement with the other tested codes and are thought to represent a good approximation to the exact solution. It should be noticed here that in this

Table 8 Flux of group 1 and group 4 at different points for OBLONG reactor in three-dimensional Methods

Time (s)

Point (16, 64, 32)

Point (88, 16, 32)

Point (16, 64, 88)

Point (88,16,88)

Group 1

Group 4

Group 1

Group 4

Group 1

Group 4

Group 1

Group 4

3DKIN (Dt = 0.001) AMF (Dt = 0.001) AMF (Dt = 0.0001)

0.0

1.402 1.402 1.402

0.112 0.112 0.112

0.384 0.384 0.384

2.742 2.742 2.742

1.772 1.772 1.772

0.142 0.142 0.142

0.486 0.486 0.486

3.467 3.467 3.467

3DKIN (Dt = 0.001) AMF (Dt = 0.001) AMF (Dt = 0.0001)

0.05

1.419 1.410 1.417

0.114 0.113 0.114

0.390 0.387 0.389

2.781 2.765 2.778

1.795 1.786 1.793

0.144 0.143 0.144

0.497 0.494 0.496

3.764 3.749 3.762

3DKIN (Dt = 0.001) AMF (Dt = 0.001) AMF (Dt = 0.0001)

0.10

1.438 1.428 1.436

0.115 0.114 0.115

0.396 0.393 0.395

2.824 2.806 2.821

1.820 1.810 1.818

0.146 0.145 0.146

0.509 0.506 0.508

4.112 4.095 4.110

3DKIN (Dt = 0.001) AMF (Dt = 0.001) AMF (Dt = 0.0001)

0.15

1.459 1.449 1.458

0.117 0.116 0.117

0.403 0.400 0.402

2.872 2.853 2.870

1.848 1.837 1.847

0.148 0.147 0.148

0.523 0.519 0.522

4.521 4.500 4.520

3DKIN (Dt = 0.001) AMF (Dt = 0.001) AMF (Dt = 0.0001)

0.20

1.484 1.472 1.483

0.119 0.118 0.119

0.410 0.407 0.410

2.928 2.906 2.926

1.881 1.869 1.880

0.151 0.150 0.151

0.538 0.534 0.537

5.008 4.982 5.007

3DKIN (Dt = 0.001) AMF (Dt = 0.001) AMF (Dt = 0.0001)

0.25

1.484 1.485 1.486

0.119 0.119 0.119

0.411 0.411 0.411

2.930 2.931 2.932

1.881 1.883 1.883

0.151 0.151 0.151

0.538 0.538 0.538

5.012 5.019 5.020

3DKIN (Dt = 0.001) AMF (Dt = 0.001) AMF (Dt = 0.0001)

0.30

1.486 1.486 1.486

0.119 0.119 0.119

0.410 0.411 0.411

2.934 2.932 2.933

1.883 1.883 1.884

0.151 0.151 0.151

0.539 0.538 0.538

5.019 5.021 5.022

116

A.E. Aboanber, A.A. Nahla / Annals of Nuclear Energy 34 (2007) 103–119

Fig. 7. Flux in three-dimensional for sinusoidal change for OBLONG.

four-group problem the cross sections are deliberately chosen to be drastically different in different materials. 4. Conclusions In this work, a computer neutronic code with a new scheme based on adaptive matrix formation (AMF) is developed for the numerical solution of the transient multigroup neutron diffusion and delayed precursor equations in two- and three-dimensional geometry. This code is characterized by the expansion of the coefficients matrix for all mesh points of the reactor into two different matrices, specified for each mesh point of the reactor. The algorithm is tested with three different benchmark problems in twoand three-dimensional geometry and the results are compared with those provided by other codes, with an overall agreement and good performance. The stability of the method has been tested for sinusoidal, ramp, and step-change reactivity insertions. These

include subprompt critical transients with highly localized perturbations. During any part of the transient when the rate of reactivity change is suddenly altered, the time steps must be decreased in size. This is necessary if accuracy is to be retained while the frequencies again seek the new rates of flux change. This must be done to control the damped oscillations that arise if a time step too large is used through this part of the transient. Although, the AMF code has a large time-step range between transients in some applications, only a very small time is required to run on the computer system. The CPU time required by AMF code is 0.0011 s per time step (Dt = 0.00025 s) compared with the same method without partition of the matrix which required 55.14 s per time step (Dt = 0.001 s). We can conclude that AMF is a versatile fast code which can be coupled as neutronic module to a standard thermalhydraulic code involving feedback effect. This task is now being carried out.

Appendix Data for homogenous reactor (two-energy group and one-precursor group) Dimensions Two-dimension

Three dimension

Group 1

D (cm) Ra ðcm1 Þ m Rf ðcm1 Þ Rj!jþ1 ðcm1 Þ s Velocity (cm/s) Fission spectrum v Precursor constant

1.35 0.001382 2.41 0.000242 0.0023 3.0 · 107 1.0 f11 = 1.0 k1 = 0.08 U1 = 0 and U2 = 0 at x = 0, Dx = 20 and Dy = 20

Geometry (cm)

Group 2

Group 1

1.08 0.0056869 2.41 0.00428 0.0 2.2 · 105 0.0 f21 = 0.0 b1 = 0.0064 200 and y = 0, 200

Group 2

1.35 1.08 0.001382 0.0054869 2.41 2.41 0.000242 0.00408 0.0023 0.0 3.0 · 107 2.2 · 105 1.0 0.0 f11 = 1.0 f21 = 0.0 k1 = 0.08 b1 = 0.0064 U1 = 0 and U2 = 0 at x = 0, 200; y = 0, 200 and z = 0, 200 Dx = 20; Dy = 20 and Dz = 20

Data for TWIGL reactor (two-energy group and one-precursor group) Dimensions Two-dimension Material properties

Material 1, 2 Group 1

D (cm) Ra ðcm1 Þ m Rf ðcm1 Þ Rj!jþ1 ðcm1 Þ s Velocity (cm/s) Fission spectrum Precursor constant Geometry (cm)

Three-dimension Material 3 Group 1

Group 1

0.4 1.3 0.5 0.15 0.008 0.05 2.1877 2.1877 2.1877 0.1 0.0015 0.03 0.0 0.01 0.0 2.0 · 105 1.0 · 107 2.0 · 105 0.0 1.0 0.0 f11 = 1.0 f21 = 0.0 k1 = 0.08 b1 = 0.0075 U1 = 0 and U2 = 0 at x = 0, 160 and y = 0, 160 Dx = 8 and Dy = 8

1.4 0.01 2.40 0.0035 0.01 1.0 · 107 1.0

Material 3 Group 2

Group 1

Group 2

0.4 1.3 0.5 0.15 0.008 0.05 2.40 2.40 2.40 0.1 0.0015 0.03 0.0 0.01 0.0 2.0 · 105 1.0 · 107 2.0 · 105 0.0 1.0 0.0 f11 = 1.0 f21 = 0.0 k1 = 0.08 b1 = 0.0075 U1 = 0 and U2 = 0 at x = 0, 160; y = 0, 160 and z = 0, 160 Dx = 8; Dy = 8 and Dz = 20

117

Group 2

1.4 0.01 2.1877 0.0035 0.01 1.0 · 107 1.0

Group 2

Material 1, 2

A.E. Aboanber, A.A. Nahla / Annals of Nuclear Energy 34 (2007) 103–119

Material properties

118

Data for LUMAC reactor (four-energy group and one-precursor group) Dimensional material properties Two-dimension Material 1

Material 2

Group 1 Group 2 Group 3 Group 4 Group 1 Group 2 2.7778 0.00266 1.4507 0.00136 0.0586 1.0 · 109 0.755 f11 = 0.0;

1.0753 0.00297 1.4507 0.00197 0.0828 1.0 · 108 0.245 f21 = 1.0;

Group 3 Group 4 Group 1 Group 2 Group 3 Group 4

0.64103 0.1626 3.3333 1.3889 0.83333 0.20833 4.1667 2.0833 0.0359 0.655 0.00135 0.0014 0.0176 0.332 0.00077 0.00072 1.4507 1.4507 1.4507 1.4507 1.4507 1.4507 0.0 0.0 0.0262 0.540 0.0007 0.0009 0.0131 0.274 0.0 0.0 0.0850 0.0 0.0586 0.0828 0.0850 0.0 0.0570 0.0822 5.0 · 106 2.0 · 105 1.0 · 109 1.0 · 108 5.0 · 106 2.0 · 105 1.0 · 109 1.0 · 108 0.0 0.0 0.755 0.245 0.0 0.0 0.755 0.245 f31 = 0.0 and f41 = 0.0; k1 = 0.08; b1 = 0.0064 U1 = 0; U2 = 0; U3 = 0 and U4 = 0 at x = 0, 160 and y = 0, 80; Dx = 8 and Dy = 8

1.0753 0.00051 0.0 0.0 0.0847 5.0 · 106 0.0

0.26247 0.012 0.0 0.0 0.0 2.0 · 105 0.0

Data for LUMAC reactor (four-energy group and one-precursor group) Dimensional material properties Three-dimension Material 1

Material 2

Group 1 Group 2 Group 3 Group 4 Group 1 Group 2 D (cm) Ra ðcm1 Þ m Rf ðcm1 Þ Rj!jþ1 ðcm1 Þ s Velocity (cm/s) Fission spectrum Precursor constant Geometry (cm)

2.7778 0.00266 1.6 0.00136 0.0586 1.0 · 109 0.755 f11 = 0.0;

Material 3 (m) Group 3 Group 4 Group 1 Group 2 Group 3 Group 4

1.0753 0.64103 0.1626 3.3333 1.3889 0.83333 0.20833 4.1667 0.00297 0.0359 0.655 0.00135 0.0014 0.0176 0.332 0.00077 1.6 1.6 1.6 1.6 1.6 1.6 1.6 0.0 0.00197 0.0262 0.540 0.0007 0.0009 0.0131 0.274 0.0 0.0828 0.0850 0.0 0.0586 0.0828 0.0850 0.0 0.0570 1.0 · 108 5.0 · 106 2.0 · 105 1.0 · 109 1.0 · 108 5.0 · 106 2.0 · 105 1.0 · 109 0.245 0.0 0.0 0.755 0.245 0.0 0.0 0.755 f21 = 1.0; f31 = 0.0 and f41 = 0.0; k1 = 0.08; b1 = 0.0064 U1 = 0; U2 = 0; U3 = 0 and U4 = 0 at x = 0, 160; y = 0, 80 and z = 0, 120; Dx = 8;

2.0833 0.00072 0.0 0.0 0.0822 1.0 · 108 0.245

1.0753 0.00051 0.0 0.0 0.0847 5.0 · 106 0.0

Dy = 8 and Dz = 8

0.26247 0.012 0.0 0.0 0.0 2.0 · 105 0.0

A.E. Aboanber, A.A. Nahla / Annals of Nuclear Energy 34 (2007) 103–119

D (cm) Ra ðcm1 Þ m Rf ðcm1 Þ Rj!jþ1 ðcm1 Þ s Velocity (cm/s) Fission spectrum Precursor constant Geometry (cm)

Material 3 (m)

A.E. Aboanber, A.A. Nahla / Annals of Nuclear Energy 34 (2007) 103–119

References Aboanber, A.E., Nahla, A.A., 2006. Solution of two-dimensional space– time multigroup reactor kinetics equations by generalized Pade´ and cut-product approximations. Ann. Nucl. Energy 33, 209–222. Adams, C.H., 1977. Current trends in methods for neutron diffusion calculations. Nucl. Sci. Eng. 64, 552–562. Azmy, Y.Y., 1992. A nodal integral method for neutron diffusion in hexagonal geometry. In: Proc. Topical Meeting on Advances in Reactor Physics, 1.L, Charleston, SC, March 8–11, pp. 509–519. Bucknel, M.R., Stewart, J.W., 1976. Multidimensional space–time nuclear-reactor kinetics studies; Part I: Theoretical. Nucl. Sci. Eng. 59, 289–297. Chao, Y.-A., Attard, A., 1985. A resolution of the stiffness problem of reactor kinetics. Nucl. Sci. Eng. 90, 40–46. Chao, Y.-A., Shatilla, Y.A., 1994. The theory of ANC-H: a hexagonal nodal diffusion code using conformal mapping. Topical Meeting on Advances in Reactor Physics, IlL 324-336, Knoxville, TN, April. Chen, G.S., Christenson, J.M., Yang, D.Y., 1997. Application of preconditioned transpose-free quasi-minimal residual method for 2group reactor kinetics. Ann. Nucl. Energy 24 (5), 339–359. DeLorey, T.F., 1993. A transient, quadratic nodal method for triangularZ geometry. Ph.D.Thesis, Nuclear Engineering, Massachusetts Institute of Technology. Diamond, D.J., 1996. Multidimensional reactor kinetics modeling. In: Proc. OECD/CSNI Workshop on transient thermal-hydraulic and neutronic codes requirements, Annapolis, USA. Ferguson, D.R., Hansen, K.F., 1973. Solution of the space-dependent reactor kinetics equations in three dimension. Nucl. Sci. Eng. 51, 189–205. Gehin, J.C. 1992. A quasi-static polynomial nodal method for nuclear reactor analysis. Ph.D. Thesis, Nuclear Engineering, Massachusetts Institute of Technology. Glasstone, S., Edlund, M.C., 1952. The Elements of Nuclear Reactor Theory. D. Van Nostrand, New York. Hageman, L.A., Yasinsky, J.B., 1969. Comparison of alternating-direction time-differencing methods with other implicit methods for the solution of the neutron group-diffusion equations. Nucl. Sci. Eng. 38, 8–32. Henry, A.F., Vota, A.V., 1965. WlGL2 – a program for the solution of the one-dimensional, two-group, space–time diffusion equations accounting for temperature, xenon, and control feedback. Bettis Atomic Power Laboratory Report WAPD-TM-532. Kang, C.M., Hansen, K.E., 1973. Finite element methods for reactor analysis. Nucl. Sci. Eng. 51, 456–495.

119

Kavenoky, A., Lautard, J.J., 1982. The neutron kinetics and thermalhydraulic transient computational module of the neptune system: CRONOS. In: Proc. Topical Meeting on Advances in Reactor Physics and Core Thermal Hydraulics, Kiamesha Lake, NY, 22–24 September, pp. 781–792. Langenbuch, S., Werner, W., Maurer, W., 1977a. Coarse-mesh fluxexpansion method for the analysis of space–time effects in large light water reactor cores. Nucl. Sci. Eng. 63, 437–456. Langenbuch, S., Werner, W., Maurer, W., 1977b. High-order schemes for neutron kinetics calculations, based on a local polynomial approximation. Nucl. Sci. Eng. 64, 508–516. Lautard, J.J., Loubiere, S., Fedon-Magnaud, C., 1991. Three dimensional pin by pin core diffusion calculation. In: Proc. International Topical Meeting on Advances in Mathematics, Computations, and Reactor Physics, 2, Pittsburgh, 28 April–2 May, pp. 6.11, 1–12. Lawrence, R.D., 1986. Progress in nodal methods for the solution of the neutron diffusion and transport equations. Prog. Nucl. Energy 17, 271–301. Nakamura, S., 1977. Computational Methods in Engineering and Science. Wiley and Sons Inc., New York. Reed, Wm.H., Hansen, K.F., 1970. Alternating direction methods for the reactor kinetics equations. Nucl. Sci. Eng. 41, 431–442. Rohan, P.E., Kang, C.M., Shober, R.A., Wagner, S.G., 1975. A multidimensional space–time kinetics code for PWR transients. In: Proc. Conf. on Computational Methods in Nuclear Engineering, VI, pp. 69– 82. Stacey Jr., W.M., 1969. Space–Time Nuclear Reactor Kinetics. Academic Press, New York. Taiwo, T.A., Khalil H.S., 1992. An improved quasistatic option for the DIF3D Nodal Kinetics Code. In: Proc. Topical Meeting on Advances in Reactor Physics, vol. 2, pp. 469–481. Utku, H., Christenson, J.M., 1994. An investigation of the temporal subdomain method for solving the finite element formulation of the space–time reactor kinetics equations. Nucl. Sci. Eng. 116, 55–66. Varga, R.S., 1962. Matrix Iterative Analysis. Prentice-Hall, Englewood Cliffs, NJ. Weinberg, A.M., Wigner, E.P., 1958. The Physical Theory of Neutron Chain Reactors. University of Chicago Press, Chicago. Wight, A.L. et al., 1971. Application of alternating-direction implicit methods to the space-dependent kinetics equations. Nucl. Sci. Eng. 44, 239–251. Zimin, V.G., Ninokata, H., 1996. Acceleration of the outer iterations of the space dependent neutron kinetics equations solution. Ann. Nucl. Energy 23 (17), 1407–1420.