Abstract
With the rapid development of nuclear power industry, simulating and analyzing the reactor transient are of great significance for the nuclear safety. The traditional diffusion theory is not suitable for small volume or strong absorption problem. In this paper, we have studied the application of discrete ordinates method in the numerical solution of space-time kinetics equation. The fully implicit time integration was applied and the precursor equations were solved by analytical method. In order to improve efficiency of the transport theory, we also adopted some advanced acceleration methods. Numerical results of the TWIGL benchmark problem presented demonstrate the accuracy and efficiency of this methodology.
1. Introduction
In reactor transient, neutron’s movement can be described by the time-dependent Boltzmann transport equation or the simplified diffusion equation in many cases. For example, when control rod ejects from the core, the flux and power will rise to a dangerous degree during a very short time. Diffusion theory is widely used in reactor transient analysis because it needs less computing resources. However, diffusion theory has its own shortage under special circumstances, for example, the small and high-leakage reactor [1]. Except for the long computing time, transport theory performs well in these problems.
Numerical calculation methods and computer technology make calculation faster. Transport theory has greater and greater potential in solving the space-time reactor kinetics equation. Some famous international codes based on transport theory have been developed in recent years, for instance, TDTORT [2] in 2001 and DORT-TD [3] in 2003.
There are mainly three methods [4] in solving the kinetics equation. The simplest is the direct numerical processing method, which transforms transport equation into algebraic equations by applying finite difference scheme. The second is the node method. It expands the angular neutron flux by introducing a series of trial functions. The quasistatic method or the improved quasistatic method decomposes neutron flux into amplitude function and shape function. Actually, the latter two methods are more efficient because they introduce several approximations, which in turn make them not as precise as the direct numerical processing. But recently, the new generation of nuclear power stations put forward higher requirements to numerical methods. Take AP1000 as an example. There are many complicated cells in AP1000 assembly, for example, the IFBA bars and the WABA bars. Besides, by using various acceleration techniques, the calculation time of direct numerical processing is completely acceptable. This is also the reason why we choose the fully implicit method on the basis of transport theory.
The fully implicit method for time variable has not been widely applied in mature transport code. DORT-TD and TORT-TD [5] developed by GRS in Germany are the most influential ones. In addition, Wu et al. [6], Liu et al. [7], and Liao [1] in China have also done some research and developed corresponding codes. We did some related analysis about this methodology in detail. To test accuracy and efficiency, a transient code was also developed. The influence of acceleration methods and important parameters was presented.
2. Solution Method
Considering the delayed neutron, the time-dependent Boltzmann transport equation has the following form: The physical meaning of each parameter shows the following: = angular neutron flux, = scalar neutron flux, = neutron velocity, = delayed neutron group number, = time, = macroscopic total cross section, = macroscopic scattering cross section, = energy spectrum for prompt neutron, = energy spectrum for precursor group , = decay constant for precursor group , = concentration for precursor group , = fraction for delayed neutron group , = external source term.
As for difference scheme, explicit, semi-implicit, and implicit schemes have been adopted. In order to achieve stability, the explicit scheme needs extremely small time step, wasting plenty of time. And the semi-implicit scheme is relatively complicated. So the unconditionally stable implicit scheme is chosen in this paper.
Let denote time point, let denote the coefficient matrix, and let denote neutron flux or concentration for precursor group , and implicit method can be described as follows: where After applying (2) to (1), the following equations can be obtained: Then we redefined , , and : Equation (4) can be written as
It can be seen from (5) and (7) that if and are known from time point , can be calculated by (7). Then it is quite convenient to get by solving (5) through analytical method. So if the original and are known, solution for the whole time domain is available.
As for energy, the multigroup approximation is adopted. The space variable is handled by finite difference method. With regard to direction, we use the discrete ordinates method. It has been applied in many famous codes, for example, ANISN [8], DORT [9], and TORT [10]. ARES [11] is a multigroup particle transport code system with arbitrary order anisotropic scattering. ARES2D is the 2D transport module of ARES code system. It is used to solve the first-order form of the steady-state Boltzmann transport equation in the two-dimensional geometry with uniform rectangular mesh. ARES2D can do both fixed source calculation and critical calculation. The methodology in this paper is based on the ARES2D.
We also want to say some words about the solution of (7). It is obviously a fixed source problem in each time step because the external source is only determined by precursor concentration and neutron flux of previous time step. However, the subroutine for fixed source calculation did not consider the fission source. An additional power iteration was added outside the primary inner iteration. The convergence criterion is where denotes mesh, denotes the number of outer iteration, denotes convergence criteria for power iteration, and denotes the fission density: where is determined by the user and it is set to in this paper.
3. Measures for Accuracy
The flux changes fast and there are lots of time steps in transient calculation. For example, if the total time is 5 s and the fixed time step is 1 ms, 5000 steps are needed to complete this calculation. Without rigorous convergence criteria, the error will accumulate rapidly. The convergence criteria for power iteration are and for eigenvalue are .
Besides, under the actual situation, the reactivity introduced is quite small, and the value changing at each step is about -. Such small change will enlarge the truncation error easily. So we applied the double-precision data structure.
4. Acceleration Technique
As is shown in Section 3, the rigorous convergence criteria and thousands of time steps will result in considerable calculation time. This is the most serious weakness of fully implicit method. In order to decrease the calculation time, it is necessary to use these methods below.
In general, the flux of step is passed to step directly. This is enough when the flux changes smoothly. However, in most transient analyses, that is invalid. We may estimate an approximate flux for step via a simple relationship: where is a special parameter that measures the rate of change of the angular flux. Then where is the approximate flux needed. If fixed time step is adopted, we can get the following simpler equation:
In addition, when is relatively small, that is to say, the gradient is small, we may adopt a larger time step. But the fully implicit method demands a very fine time step size, as already mentioned above. Considering this contradiction, the adaptive time step technique works very well. Another parameter is the error limit . It is determined by user. A reasonable value can guarantee a satisfactory result.
For outer iteration, many acceleration methods can be introduced. One of them is the fission source extrapolation method [1]. The “fission source” here is actually the fission density mentioned in Section 2: where is the fission density used in time step , is the new fission density calculated after time step , and is the extrapolation coefficient: where is the so-called dominance ratio, which is estimated by the following relationship:
After accomplishing these steps, fission source for calculation is calculated: where is fission spectrum.
5. Results and Analysis
A concise calculation flow chart for the transient code is shown in Figure 1.
As the results can be compared with other codes, the famous TWIGL benchmark problem has been calculated. The geometry of the reactor is shown in Figure 2.
TWIGL problem is a 2D transient benchmark based on diffusion theory, which was originally put forward by Hageman and Yasinsky in 1969 [12]. It models a 160.0 cm square seed-blanket reactor which consists of three materials. The primary fissile material is filled in seeds 1 and 2. The reactivity is introduced in region 1. The central and outer regions are blanket, which also contains fissile material. The geometry of the reactor is shown in Figure 2 and the information of cross section as well as precursors is given in Tables 1 and 2. They are all taken from [2].
Two kinds of perturbations are modeled. The step perturbation is introduced by reducing the of material 1 by 0.0035 at time zero. In the ramp perturbation, the same is done for but the duration is from time 0.0 s to time 0.2 s. The total time for both perturbations is 0.5 s. The reference values are taken from [2]. The normalized reactor power versus time or the flux distribution will be presented in following analysis.
The computer hardware configuration is i3-2120, 3.3 GHz, and 4G. It is worth mentioning that, before the transient calculation, we always perform a steady calculation first and then divide fission cross section by to make system critical. Figures 3 and 4 are the normalized power in the step and ramp perturbation.
- approximation is introduced: corresponds to the order of the expansion in Legendre polynomials of the scattering cross-section matrix and represents the order of the flux angular discretization. Fully symmetrical quadrature sets were introduced.
The computing time is 27 seconds and 24 seconds, respectively. We can see from the two figures that the results are in good agreement with reference. However, there is a slight error accumulation during the later stage of transient. This can be caused by the collocation of several parameters and acceleration methods, which will be presented later. Besides, the power in both cases in 0.5 seconds is closed to each other. In other words, the same reactivity change, no matter what the speed is, will bring in the same effect.
Tables 3 and 4 show normalized power with different discretization order. Discretization order has little influence on the results but the approximation will save much time.
We have applied the uniform mesh size. The results for different mesh size are given in Tables 5 and 6. The time step is kept in 0.1 ms.
It is obvious that the smaller the mesh size is, the longer the calculation time will be. Moreover, the time is inversely proportional to the square of the mesh size approximately. In addition, the power deviation in small mesh size gets bigger. This is because the error will accumulate with superabundant mesh in steady calculation. So considering efficiency, we do not recommend the excessively fine mesh except for some special cases, for example, the control rod area.
Tables 7 and 8 show the power distribution with different time step. The mesh size is 8 cm. It can be seen from the two tables that when the time step size is extremely small, the step number is enormous, which will also result in error accumulation. And when the time step is 0.5 ms, the calculation time is only a few seconds. For discrete ordinates method, this is completely acceptable. Unfortunately, the results become unreasonable when the time step gets bigger.
We have also adopted several acceleration methods. When time step is 0.1 ms and the mesh size is 8 cm, we get the results of comparison.
Tables 9 and 10 show the difference between with flux estimation and without flux estimation. We can see clearly that in step perturbation the speedup ratio is 2.6 and in ramp perturbation it is 7.09. This is because the fundamental of flux estimation is , which is similar to the inverse reactor period. By constructing the reactor period in step , we estimate the flux of step . For step perturbation, the power changes too fast. The estimation is not as good as ramp perturbation.
The influence of fission source extrapolation is shown in Tables 11 and 12. Results show that the extrapolation method does not work in TWIGL problem. That is related to the system dominance ratio. When the dominance ratio is quite small, the extrapolated fission source will not be calculated. But this method will be effective for big dominance ratio system.
The precursor concentration and flux distribution in transient process are shown in Figures 5, 6, and 7. Because the precursor concentration does not change much in the transient, we only give the steady distribution here. As the ramp perturbation is similar to the step perturbation, only the power distribution in step perturbation is illustrated.
(a) s
(b) s
(c) s
(d) s
(a) s
(b) s
(c) s
(d) s
The thermal neutron flux in seed area is obviously lower than that in blanket area. Fast neutron flux in area 3 is sunken, so does the precursor concentration. Because of the delayed neutron, the distribution and amplitude of the flux are violently changed from 0.0 seconds to 0.1 seconds and keep smooth in the next 0.4 seconds. But the first 0.1 seconds takes plenty of calculation time.
6. Conclusion and Outlook
In this paper, we have done some research on the numerical methods for transient based on discrete ordinates method. The results of the code are in good agreement with the reference in TWIGL problem. As a method of solving the transport equation, compared with the diffusion theory, the discrete ordinates method with fully implicit method spends relatively more time generally. But by applying various acceleration methods, the calculation time can be a few seconds, which is totally acceptable. The dual accuracy of the fully implicit method and transport theory will make this method play a major role in many fields in the future.
Among these methods, the fission source extrapolation only works in system with big dominance ratio. The flux estimation method performs well and its form is very simple when uniform mesh is used. The excessively fine mesh will result in slight error accumulation and large amount of time. The time is inversely proportional to the square of the mesh size approximately. So the time step cannot be too large, which is subjected to the feature of fully implicit method.
As for future work, more advanced acceleration methods will be adopted to improve efficiency, for example, the multigrid and the diffusion synthetic acceleration. There is still great potential in decreasing the calculation time. To simulate the real situation, this method should be coupled with thermal hydraulic. The neutron velocity and effective delayed neutron fraction will be calculated accurately. We will also try to apply this transient method in core to fully exert its advantage.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
The authors are grateful for the support of this research by the Mayor National S&T Specific Program of Large Advanced Pressurized Water Reactor Nuclear Power Plant (2011ZX06004-007).