- Email: [email protected]
Reactivity insertion limits in a typical pool-type research reactor cooled by natural circulation
annals of
NUCLEAR ENERGY Annals of Nuclear Energy 33 (2006) 252–261 www.elsevier.com/locate/anucene
Reactivity insertion limits in a typical pool-type research reactor cooled by natural circulation H. Kazeminejad
*
Gamma Irradiation Center, Atomic Energy Organization of Iran, P.O. Box 11365-3486, Tehran, I.R. Iran Received 29 August 2005; received in revised form 10 October 2005; accepted 14 November 2005 Available online 27 December 2005
Abstract The course of reactivity insertion in a pool type research reactor, with scram disabled under natural circulation condition is numerically investigated. The analyses were performed by a coupled kinetic–thermal–hydraulic computer code developed specifically for this task. The 10-MW IAEA MTR research reactor was subjected to unprotected reactivity insertion (step and ramp) for both low and highenriched fuel with continuous reactivity feedback due to coolant and fuel temperature effects. In general, it was found that the power, core mass flow rate and clad temperature under fully established natural circulation are higher for high-enriched fuel than for low enriched fuel. This is unlike the case of decay heat removal, where equal clad temperatures are reported for both fuels. The analysis of reactivity represented by the maximum insertion of positive reactivity ($0.73) demonstrated the high inherent safety features of MTR-type research reactor. Even in the case of total excess reactivity without scram, the high reactivity feedbacks of fuel and moderator temperatures limit the power excursion and avoid consequently escalation of clad temperature to the level of onset of nucleate boiling and sub-cooled void formation. The code can also be modified to provide an accurate capability for the analyses of research reactor transients under forced convection. Ó 2005 Elsevier Ltd. All rights reserved.
1. Introduction The complex set of physical phenomena that occur in a gravity environment when geometrically distinct heat sink and heat source are connected by a fluid flow path can be identified as natural circulation (NC). In all the adopted configuration of existing water-cooled reactors, NC allows the removal of the decay heat produced by the core. Furthermore, NC is the working mode for steam generators in existing pressurized heavy and light water reactors. It is necessary also for the core cooling system in an unlikely loss of primary coolant accident to occur in all the reactor types. The new frontier of nuclear technology represented by advanced reactor designs, with the aim of simplifying the systems, of reducing costs and of increasing the safety level, *
Tel.: +98 21 88004065; fax: +98 21 88009054. E-mail address: [email protected].
0306-4549/$ - see front matter Ó 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.anucene.2005.11.002
foresees an extended use of the NC principles. This includes, in some cases, the nominal operation of the reactors and, in all cases, the reactor cooling by passive systems (Jafari et al., 2003). Although the application of passive safety systems is a potential means of achieving simplification and competitive economics in new nuclear power plant designs, the use of passive systems is not entirely new, and is not unique to any particular line of new reactor designs. But an increased reliance on this approach, making safety functions less dependent on active components like pumps and diesel generators, is potentially an important means to achieve reduced costs for future nuclear power plants. Many new reactor designs incorporate passive systems based on NC. The study of NC in research reactors has up to now been mainly focused on decay heat removal from the core (Langerman, 1994; Woodruff et al., 1996; Housiadas, 2000; Hamidouche et al., 2004) and no consideration is given to the issue of reactivity insertion under NC. The main objec-
H. Kazeminejad / Annals of Nuclear Energy 33 (2006) 252–261
253
Nomenclature a cp Ci De f fa fe fr F(t) g G hs hg k L m N p P q00i q00l q000 q000 q000 c r r0 R
constant for planar or cylindrical geometry specific heat ith group precursor density equivalent diameter friction factor axial peaking factor engineering peaking factor radial peaking factor external neutron source gravitational acceleration defined in Eq. (3) surface heat transfer coefficient defined in Eq. (6) thermal conductivity fuel length time step number of collocation points reactor power pressure incipient boiling heat flux coolant heat flux volumetric heat generation average volumetric heat generation volumetric heat generation at the fuel center radial coordinate dimensionless radial coordinate (r/R) fuel radius, fuel half thickness
tive of this study is to establish the state of the art in the area of NC by considering the reactor start-up and operation utilizing NC mode. Attempts were also made to establish operating limits of the single-phase NC mechanism as applied to nuclear research reactors during a transient step and ramp positive reactivity insertions. Results were obtained for 10-MW IAEA benchmark cores with both high enriched uranium (HEU) and low enriched uranium (LEU) fuel with a failure to scram. An estimate of the reactivity insertion threshold without exceeding the onset of nucleate boiling (ONB) is obtained. Also, the code developed for this purpose can easily be modified and adapted for the testing of methods and for subsequent use in the detail analysis of transient behavior (loss of flow or reactivity insertion under forced flow) in research reactors. 2. Mathematical models 2.1. Heat transfer model Many problems in a nuclear reactor analyses require the knowledge of the temperature distribution in the fuel, cladding and the coolant during normal reactor transients and
Re t T V w Yc z
Reynolds number time temperature velocity water channel thickness clad thickness axial coordinate
Greek letters b total delayed neutron fraction bi ith group of delayed neutron fraction bl coolant volumetric thermal expansion coefficient h defined in Eq. (3) ki ith group decay constant K prompt neutron life time l dynamic viscosity q reactivity or fluid density s reactor period Subscripts b bulk condition f fuel c clad l coolant o initial or reference conditions
in accident situations such as reactivity insertion accidents (RIA). Temperature distribution is an important phenomenon, which influences the fuel behavior and integrity during all stages of the facility lifetime. This information is required to assess the safety margins in fuel plates both during normal reactor transients and when exposed to hypothetical reactor accident condition. The behavior of the cladding is also very important, because it is the first barrier to release of fission products. The exact calculation of temperature variation in fuel and cladding is rather complex. The complexity arises due to the fact that the heat generation in fuel assembly, the temperature variation and thermal–hydraulic behavior of the core are interdependent. The problem is very simplified by considering a single coolant flow channel surrounded by two fuel plates and a heat transfer coefficient between the cladding and the coolant. The coolant temperature in the reactor core is taken as the boundary condition, which supplies the fuel assembly with the sink temperature. The heat generation in the fuel is assumed to have a cosine distribution or to be a known function of time. The thermal behavior of the fuel is described by heat conduction equation in fuel and clad-
254
H. Kazeminejad / Annals of Nuclear Energy 33 (2006) 252–261
ding. One-dimensional equations for mass, momentum and energy are solved in the coolant. Because of the complexity of equations, boundary conditions and temperature dependence of physical parameters, an analytical solution cannot be obtained for the differential equations. The problem must be approximated by numerical solutions. The heat transfer model and numerical procedure used in developing a general computer code suitable for research reactor are similar to thermal–hydraulic codes for power reactors (COBRA-IV, 1976) in that the transient thermal response of the fuel is coupled with the hydraulics through the heat transfer coefficients. The fuel plates heat conduction model includes axial conduction and variable thermal conductivity model. The accuracy of the present results is improved by using an orthogonal collocation numerical procedure. The physical properties of the coolant are calculated by interpolation of the liquid property table. As the fluid is heated in each computational grid, this process is repeated, employing linear interpolation where necessary, to determine the fluid properties in terms of the increasing temperature. The fuel heat transfer model allows the computation of the fuel plate heat fluxes and internal temperature distributions based on the heat generation in the fuel and the surface heat transfer coefficient. The inclusion of axial conduction and the temperature dependence of fuel thermal conductivity provide greater detail and allow consideration of wider class of problems. The method of solution of the governing differential equations is based on the method of weighted residuals (MWRs) as a replacement to the finite difference scheme in the radial direction. Finite difference scheme is still used for time derivatives and axial, z, space derivatives. In MWR, orthogonal collocation is employed to determine the form of the approximated polynomial solution. MWR affords a higher order of accuracy by using the roots of orthogonal polynomial as the nodal position where the solution is evaluated, which eventually provides accuracy comparable to finite element techniques. In what follows a brief description of the governing equations, method of solution and the results obtained from its main application will be presented and discussed.
8.05
0.038 0.223
0.051
6.30
6.65
7.60
0.475
Standard Fuel Element: 23 fuel plates Control Element: 17 Fuel Plates 4 Al Plates Fig. 1. View of the unit cell representing the 10-MW benchmark MTR fuel element (dimensions are in cm).
w
Clad L
Fuel Coolant
Z
R
V, To, P
2.1.1. Fuel interior The view of a unit cell representing the 10-MW benchmark MTR fuel element is shown in Fig. 1. A single channel geometry as shown in Fig. 2 is chosen to represent the behavior of any region of the core. The z direction is the direction of flow along the channel and the R direction is the direction normal to the plate. The heat transfer model can be derived from the fundamental heat conduction equation: oT qcp ¼ r ðkrT Þ þ q000 ; ot
ð1Þ
where k is the local thermal conductivity and q000 is the volumetric rate of heat generation.
Fig. 2. Flow model.
Eq. (1) is expanded in two dimensions for both planner and cylindrical geometry with the radial coordinate nondimensionalized using the relationship r = r 0 /R, where r 0 is the radial coordinate and R is the fuel half thickness or radius: oT 1 o oT o oT ¼ 2 a1 ra1 k f ðT Þ þ k f ðT Þ þ q000 ; ð2Þ ot or ot oz oz Rr where a = 1 for planar; a = 2 for cylindrical. Eq. (2) is reduced to a linear partial differential equation using KirchoffÕs transformation:
ðqcp Þf
H. Kazeminejad / Annals of Nuclear Energy 33 (2006) 252–261
1 h¼ ko
Z
T
kðT Þ dT ¼ GðT Þ;
ð3Þ
To
where ko is the conductivity at reference temperature To. Differentiating Eq. (3) with respect to r and t and substituting into the differential equation (2), yields the transformed equation: ðqcp Þf
k o oh k o o a1 oh o oT ¼ r þ k f ðT Þ þ q000 . or oz oz k f ðT Þ ot R2 ra1 or
ð4Þ
Eq. (4) can be solved numerically by a combination of techniques. In the orthogonal collocation method, the radial positions are taken to be the roots to orthogonal polynomials, which affords accuracy comparable to finite element techniques. The radial positions or collocation points are defined by Finlayson (1972) for a second-order approximation (N = 3). The axial conduction term is calculated by central finite difference and the time derivative is approximated by a forward finite difference. The boundary condition at the fuel surface is handled by a lumped resistance technique. The governing equation is k o oh ¼ hg ðT N T N þ1 Þ; ð5Þ R or where N and N + 1 denote the fuel and clad exterior surfaces, respectively. The term hg is defined by 1 1 Yc ¼ þ ; hg hc k c
ð6Þ
where Yc and kc are cladding thickness and conductivity and hc is the fuel-clad gap conductance if any. 2.1.2. Cladding A transient energy balance for the lumped clad is: a1 oT N þ1 hg rN hs ¼ ðqcp Þc ðT N T N þ1 Þ ðT N þ1 T l Þ ot Y c rN þ1 Yc 2
þ k c ðT Þ
o T N þ1 ; oz2
ð7Þ
where hs is the clad surface heat transfer coefficient and Tl is the local fluid temperature. To evaluate the implicit temperatures TN and TN + 1 in Eqs. (5) and (7), a truncated Taylor series is used: T ¼ Tm
GðT m Þ hðT m Þ ; G0 ðT m Þ
ð8Þ
where the function G is defined in Eq. (3) and G 0 is the derivative of G with respect to T. After applying the orthogonal collocation method to boundary equations (5) and (7), they are combined with the differential equation (2) to yield the matrix equation of the form [A]{h} = {Q} for the transformed temperatures at one axial level. This equation is solved by an iterative Gauss-Seidel procedure. Once the {h} solution vector is determined, the temperature solution is evaluated using Eq. (8) with h(Tm) specified. If any elements of [A] or {Q} are temperature dependent, the procedure is repeated
255
until the resulting temperature change remains within a set limit and the solution is considered converged. The volumetric rate of heat generation q000 and fluid temperature, Tl at each axial node are taken explicitly from the previous time step solutions of the kinetics and thermal–hydraulics, respectively. The sequence is repeated at each axial level to the ends of the fuel with additional axial sweeps if axial conduction is being considered. Normally, for plate type fuel elements the axial conduction has little effect and only two or three additional iterations are necessary. 2.2. Thermal–hydraulics model Thermal–hydraulic modeling is based on the equations of conservation of mass, momentum and energy, formulated for a one-dimensional axial homogeneous flow through a vertical channel (Fig. 2) proposed by Nakamura (1977): Continuity equation: oq oðqV Þ þ ¼ 0. ot oz
ð9Þ
Using Boussinesq approximation, the differentials in density in the continuity equation (9) are of the order of b and hence can be neglected, to give oV/oz = 0.0 as for an incompressible fluid. Langerman (1994) show that, for decay heat removal by natural circulation in reactor fuel plates, the difference between the temperature ratios calculated with the Boussinessq and the non-Boussinesq models is rather small for water flow. z-momentum equation (vertical direction): oV oV 1 oP fV 2 þV ¼ þ gbl ðT l T o Þ. ot oz qo oz 2De
ð10Þ
Energy equation: oT l oT l 2hs ðT c T l Þ þV ¼ . qcp w ot oz
ð11Þ
The Boussinesq approximation is incorporated into the momentum equation (10). It assumes that the liquid density q varies with temperature Tl according to q = qo[1 bl(Tl To)], but all other liquid properties are invariant within the specified temperature range. qo is the density at reference temperature To and bl is the coolant volumetric thermal expansion coefficient. The coupled differential equations (9)–(11) with initial conditions V(0) = 0 and Tl(0) = To, were solved by the Lax-Wendroff explicit method (Anderson et al., 1984), which is a second-order finite difference method, for time evolution of liquid velocity, V, and liquid temperature Tl. 2.3. Reactor kinetics model The reactor power dynamics under transient conditions is derived through a numerical solution of the conventional point kinetics formula. The point kinetics equations for six
256
H. Kazeminejad / Annals of Nuclear Energy 33 (2006) 252–261
delayed groups and light water moderator are (Lamarch, 1982): 6 X dp qðtÞ b ¼ pþ ki C i þ F ðtÞ; ð12Þ dt K i¼1 dC i bi ¼ p ki C i i ¼ 1; 2; . . . ; 6; ð13Þ dt K where p is the reactor power, Ci is the ith precursor density, q(t) is the time dependent P reactivity function, bi is the ith 6 delayed fraction, and b ¼ i¼1 bi is the total delayed fraction. In addition, K is the neutron generation time, ki is the ith group decay constant, and F(t) is the time dependent neutron source function. The appropriate initial conditions that accompany Eqs. (12) and (13) are: b pð0Þ ¼ po ; C i ð0Þ ¼ i po . ð14Þ ki K Eqs. (12) and (13) are a system of coupled linear ordinary differential equations. The time dependent parameters in this system are the reactivity and the neutron source terms. An important property of the point kinetics equations is the stiffness of the system. It is well known that stiffness is a severe problem in numerical solution of the point kinetics equations and results in the need for small time steps in a computational scheme. There has been a great deal of research done in the field of reactor kinetics. Much of the work has been focused on eliminating the aforementioned problem of stiffness. A modified Runge–Kutta method to avoid the problem of stiffness devised and implemented by Sanchez (1989) is used to derive the solution of the set of differential equations. This method constitutes an easy to implement algorithm that provides results with sufficient accuracy for most applications. The main advantage of this method is that it allows step size control to be performed in a systematic manner and estimation of the truncation error can be readily obtained at each time step. The reactivity is estimated by taking into account the different feedback effects as follows (Housiadas, 2002): qðtÞ ¼ qex ðtÞ þ aM dql aT dT l aF dT f .
ð15Þ
In Eq. (15), the first term, qex(t), is the externally introduced reactivity (for example, the step and ramp reactivity insertions or shutdown reactivity if scram is enabled) and the remaining terms are the various feedback contributions. As can be seen three feedback mechanisms are taken into consideration, change of coolant density dql, change of coolant temperature dTl (spectrum effects only) and changes of fuel temperature dTF (Doppler effects). The corresponding constant coefficients of reactivity are, respectively, aM, aT and aF. In applying Eq. (15) the changes in feedback reactivity have to be determined by considering the variations of the mean (whole volume) coolant and fuel temperatures. 2.4. Heat generation in fuel plates Axial variation of the volumetric heat generation is assumed to have a cosine distribution or to be given as any function of time, using the exponential representation:
q000 ¼ q000 c cos
pz . L
ð16Þ
In Eq. (16), q000 and q000 c are the axial variation of volumetric heat generation and at the fuel center, respectively. The local average heat generation depends on the axial location z along the fuel element and on the location of the fuel element in the core. The code uses a two-channel model, one channel to represent the hottest plate and flow channel and the other average channel to represent the remaining fuel plates in the core. The axial source distribution was represented by 21 axial regions and a chopped cosine shape which has an axial power peaking of, fa, for both the average channel and the hot channel. For the hot channel, the axial distribution was multiplied by the other specified hot channel factors. Taking into account these variations, the hot channel heat generation at the fuel center is given by q000 q000 fa fr fe ; c ¼
ð17Þ
where q000 is the average volumetric heat generation of the fuel in the core. To compute the hot channel peak temperatures, first an average channel power trace was obtained, and then a second power driven transient run with hot channel factors was carried out. 2.5. Heat transfer and friction factor correlations Surface heat transfer coefficient is required to couple the hydraulic solutions to the fuel and clad model. Heat transfer from the clad surface to the coolant is described by NewtonÕs Law of Cooling: q00l ¼ hs ðT c T l Þ;
ð18Þ
where ql is the convection heat flux, hs is the convective heat transfer coefficient, Tc and Tl are the local clad and liquid temperatures, respectively. hs will depend upon the properties and flow conditions of the coolant. The coolant flow between the coolant passages of a pool type research reactor core in natural convection is laminar (Reynolds number, Re 6 2000) but the existing heat transfer correlations for natural convection are inadequate and predict the heat transfer coefficients that are too low. Therefore, in line with the single-phase heat transfer model in PARET, an approximation of the heat transfer coefficient proposed by Woodruff (1984) was used 1=2 kqcP hs ¼ . s
ð19Þ
The heat transfer coefficient predicted by Eq. (19) is higher than the existing heat transfer correlations for natural convection. The use of Eq. (19) requires the selection of an appropriate period, s. This period is chosen as the minimum sm such that |sm + 1 sm|/sm is less than 0.01 where sm is the instantaneous reactor period at time step m. This has proven to be a reasonable criterion under most conditions. At certain pressure and temperature conditions, an incipient heat flux, qi, exists at which isolated vapor nucle-
H. Kazeminejad / Annals of Nuclear Energy 33 (2006) 252–261
ation occurs on the cladding surface. The Bergles–Rohsenow correlation is used to predict qi as follows (Collier, 1980): q00i ¼ 15:60P 1:156 DT as
a ¼ 2:30P 0:0234 ;
ð20Þ
where P is the absolute pressure (kPa), DTs is the degree of superheat (K) and qi is the incipient heat flux. Onset of nucleate boiling (ONB) is not a limiting criterion in the design of a fuel element. However, it is a heat transfer regime, which should be identified for proper hydraulic and heat transfer considerations, i.e., single phase flow versus two-phase flow. The actual axial location at which ONB will occur depends upon the axial heat flux distribution, the coolant velocity and the pressure drop along the channel. For simplicity, the heat flux for ONB can be calculated conservatively by using the worst combination of parameters, i.e., ONB occurs at the channel exit with peak heat flux, lowest pressure and saturation temperature, and highest coolant temperature rise. The friction factor for laminar flow between rectangular cross-sections is determined as f = 96/Re(lc/lb), where Re is the Reynolds number and subscripts b and c refer to coolant properties evaluated at the bulk and clad temperatures, respectively. 3. Benchmark problem description The benchmark problem considered in this investigation consists of a 10-MW light water, pool type research reactor specified by the IAEA (1992). Two similar core configurations with different fuel enrichments HEU (93% of U235) and LEU (20% of U235) and fuel material type but with the same geometric dimensions (Fig. 1) and operating conditions are considered. These cores have been used as benchmark for the transient analyses of research reactors. The detailed specifications are given in Table 1. The benchmark problem is performed for transients governed by the core kinetics and thermal–hydraulics. 4. Results and discussions As stated earlier, the main aim of the present work is to study the safety aspects of 10-MW benchmark MTR research reactor under NC condition. The safety of the reactor demands that the stagnant liquid inside the core does not reach its boiling point. Therefore, it is important to find the level of reactivity insertion where the liquid starts boiling in a research reactor under NC condition. The transients performed for the benchmark cores are self-limiting and the limits imposed on reactivity insertion by the ONB condition. Figs. 3a–c provide a comparison of power, reactivity, core mass flow rate and peak clad temperature for both HEU and LEU for a step reactivity insertion of $0.5. Fig. 3a shows the development of core average power and reactivity during the initial start up phase. The situa-
257
Table 1 IAEA research reactor benchmark problem operating conditions Core material Nuclear fuel Fuel element Coolant Moderator
MTR Plate type clad in AL Light water Graphite-light water
Core thermal–hydraulics Fuel thermal conductivity (W/cm K) Cladding thermal conductivity (W/cm K) Fuel specific heat (J/g K) Cladding specific heat (J/g K) Fuel density (g/cm3) Cladding density (g/cm3) Radial peaking factor Axial peaking factor Engineering peaking factor Inlet coolant temperature (°C) Operating pressure (bar) Fuel element dimensions Length (cm) Width (cm) Height (cm) Number of plates/fuel ECS/ECN Plate meat (mm) Width (cm) active/total Height (cm) Water channel thickness (mm) Plate clad thickness (mm)
8.00 7.60 60.0 23/17 0.51 6.30/6.65 60.0 2.23 0.38
Core kinetics Effective delayed neutron fraction Prompt neutron generation time (ls) Density feedback coefficient, aM ($/% void) Doppler feedback coefficient, af ($/°C) Coolant feedback coefficient, aT ($/°C)
HEU 1.58
LEU 0.5
1.80
1.80
0.728 0.892 0.68 2.7 1.4 1.5 1.2 38.0 1.7
0.340 0.892 4.45 2.7
HEU 7.607E 3
LEU 7.275E 3
55.96
43.74
0.3257
0.4047
3.6E 5
3.31E 3
1.537E 2
1.082E 2
tion simulates the reactor start up when withdrawing the control rods to increase the reactor power. After a short coolant heat up time of about 50 s, the core mass flow rate (Fig. 3b) and the peak clad temperature (Fig. 3c) start to increase due to onset of NC. At the same time, the reactivity starts to decrease (Fig. 3a) due to high negative reactivity feedback of fuel and moderator temperatures. The negative feedback effects are able to consume the entire available excess reactivity forcing the power to stabilize and demonstrating the inherent safety of MTR research reactors. For LEU, the power increases to a maximum of 188 kW at 75 s after the start of transient and then attains a nearly constant value of 155 kW after about 140 s. The variation of core mass flow rate is similar to that of power and after about 140 s in transient, it reaches a semi-constant state characterizing a fully developed NC. It is also interesting to note that in a fully developed NC, the net reactivity for both LEU and HEU cores is zero, which
258
H. Kazeminejad / Annals of Nuclear Energy 33 (2006) 252–261
means that the effect of external reactivity is almost balanced by the effects of reactivity feedbacks. In the fully developed NC regime all reactor variables reach a steady state value after about 150 s. For HEU, the peak clad temperature (Fig. 3c) is higher for LEU due to higher power produced by the HEU core despite the fact that, the mass flow rate induced by the HEU core is much higher than that of the LEU core. The peak clad temperatures for both
cores are well below the ONB condition for a reactivity insertion of $0.5. Figs. 4a–c show the time evolution for power, reactivity, core mass flow rate and peak clad temperature through a period of 300 s for a step reactivity insertion of $0.73. This reactivity insertion limit is imposed by the ONB. It should be noted that any slight increase in the value of 0.73 would cause the peak clad temperature for both cores to exceed
Fig. 3. (a) Power and reactivity; (b) mass flux; and (c) peak clad temperature for step reactivity insertion of $0.5.
H. Kazeminejad / Annals of Nuclear Energy 33 (2006) 252–261
the ONB condition and hence making the single-phase fluid flow assumption invalid. As expected, the maximum values attained by the power, core mass flow rate and clad temperature are higher as compared with $0.5 reactivity insertion. For LEU and HEU cores, the maximum powers are 518 kW and 593.5 kW occurring at about 21 s and 21.9 s, with steady-state values of 211 kW and 417 kW, respectively. The peak clad temperatures for LEU and HEU cores are 115.2 °C and 113.3 °C occurring at bout
259
51 s and 44 s, respectively. For both LEU and HEU cores, the clad temperature first decreases to a minimum value of 98 °C and 105.5 °C at about 87 s and 60.5 s before reaching a steady-state value of 108 °C and 115.6 °C, when NC is fully established. The peak clad temperature for both cores remain below the ONB condition of 116.5 °C during the entire transient, which is also shown in Fig. 4c. Figs. 5a–c provide a comparison of HEU and LEU benchmark cores for a ramp positive reactivity insertion
Fig. 4. (a) Power and reactivity; (b) mass flux; and (c) peak clad temperature for step reactivity insertion of $0.73.
260
H. Kazeminejad / Annals of Nuclear Energy 33 (2006) 252–261
of $0.73/0.5 s. The variation of core power, reactivity, core mass flow rate and peak clad temperature are almost similar to those of step reactivity insertion. This is due to the fact that the ramp insertions of short duration are equivalent to a step insertion. The entire ramp is inserted before the power, feedback and temperatures have increased substantially and the limiting reactivity insertion remains constant. Although the results for reactivity insertion of $0.73 show differences between the HEU and LEU cores at the
early stages of the transient, the peak clad temperatures reached in the fully developed NC for both cores remain below the ONB condition, with the HEU core attaining a higher peak clad temperature than that of the equivalent LEU core. This is unlike the case for decay heat removal under NC, where the reported peak clad temperatures for both cores are the same. The results for the analyses of loss of flow transients for an MTR research reactor benchmark problem using PARET and RELAP5/3.2 mode are
Fig. 5. (a) Power and reactivity; (b) mass flux; and (c) peak clad temperature for ramp reactivity insertion of $0.73/0.5 s.
H. Kazeminejad / Annals of Nuclear Energy 33 (2006) 252–261
reported in the literature by Woodruff et al. (1996) and Hamidouche et al. (2004). In the flow decay, due to the combined effect of constant decay heat and continuous reduction of the core flow rate, the peak-clad temperature exhibits a second rise after the scram of the reactor power. The increase is further sustained as the flow regime passes to laminar regime and to mix flow when the natural convection valve opens. The peak-clad temperature begins to decrease only when the NC flow is fully established. Indeed, this phenomenon is governed by the balance between the pressure and the buoyancy forces involved during the flow decay process. For both HEU and LEU, the peak-clad temperatures for fully established NC predicted by PARET and RELAP codes in slow and fast flow decay processes are the same. 5. Conclusions The results for a simple physical model developed to simulate thermal–hydraulic and neutron dynamic features for a 10-MW IAEA research reactor under NC mode are presented. Two cores with similar configuration and operation but different fuel enrichments (HEU and LEU) are considered. The analyses were carried out for self limiting cases and the limits imposed on the reactivity insertion by the ONB clad surface temperature. Based on model predictions, it can be concluded that a research reactor with either HEU or LEU core can be safely operated under NC mode for a reactivity insertion as high as $0.73 without exceeding the ONB condition. The results for step reactivity insertion were found to be similar to those of ramp reactivity insertion with short duration. For the estimated reactivity insertion, the peak clad temperatures for both HEU and LEU cores under fully established NC attained a constant value with the HEU core predicting a higher peak clad temperature (115.6 °C) than that of the equivalent LEU (108 °C) core. However, for decay heat removal under NC condition, the peak-clad temperatures reported in the literature for both cores are almost equal. The high negative reactivity feedbacks of fuel and moderator temperature were able to consume the entire available excess reactivity, limiting the power excursion and forcing the reactor to stabilize during the reactor start-up in NC mode. The results presented clearly demonstrated the inherent safety features of MTR research reactors utilizing NC for the recommended reactivity insertion. Also, the use of relatively simple physical models were shown to provide a
261
useful capability for the analyses of most research reactors encountered in practice under transient conditions. Acknowledgements This work is supported by the research grant from Atomic Energy Organization of I.R. Iran (AEOI) whose assistance is hereby gratefully acknowledged. References Anderson, D.A., Tannehill, J.C., Pletcher, R.H., 1984. Computational fluid mechanics and heat transferSeries in Computational Methods in Mechanics and Thermal Sciences. Hemisphere, New York. COBRA-IV, 1976. An interim version of COBRA for thermal–hydraulic analysis of rod bundle nuclear fuel elements and cores, BNWL-1962. Collier, J.G., 1980. Convective Boiling and Condensation. McGraw-Hill, New York. Finlayson, B., 1972. The method of weighted residuals and variational principlesMathematics in Science and Engineering Series, vol. 87. Academic Press, New York. Hamidouche, T., Bousbia-Salah, A., Adorni, M., DÕAuria, F., 2004. Dynamic calculations of the IAEA safety MTR research reactor benchmark problem using RELAP5/3.2 code. Annals of Nuclear Energy 31, 1385–1402. Housiadas, C., 2000. Simulation of loss of flow transients in research reactors. Annals of Nuclear Energy 27, 1683–1693. Housiadas, C., 2002. Lumped parameters analysis of coupled kinetics and thermal–hydraulics for small reactors. Annals of Nuclear Energy 29, 1315–1325. IAEA-TECDOC-643, 1992. Research reactor core conversion guide book, vol. 3, Analytical Verification. Jafari, J., DÕAuria, F., Kazeminejad, H., Davilu, H., 2003. Reliability evaluation of a natural circulation system. Nuclear Engineering and Design 224, 79–104. Lamarch, J.R., 1982. Introduction to Nuclear Engineering, second ed. Addison-Wiley, Reading, MA. Langerman, M.A., 1994. Calculating decay heat removal rates via natural circulation along variable heat flux reactor fuel plates. Nuclear Engineering and Design 150, 61–68. Nakamura, S., 1977. Computational Methods in Engineering and Science with Applications to Fluid Dynamics and Nuclear Systems. Wiley– Interscience, New York, pp. 171–175. Sanchez, J., 1989. On the numerical solution of the point reactor kinetics equations by generalized Runge–Kutta method. Nuclear Science and Engineering 103, 94–99. Woodruff, W.L., 1984. A kinetics and thermal–hydraulics capability for the analysis of research reactors. Nuclear Technology 64, 196–206. Woodruff, W.L., Hanan, N.A., Smith, R.S., Matos, J.E., 1996. A comparison of the PARET/ANL and RELAP/MOD3 codes for the analysis of IAEA benchmark transients. In: International Meeting on Reduced Enrichment for Research Reactors, Seoul, Republic of Korea, October 7–10.
NUCLEAR ENERGY Annals of Nuclear Energy 33 (2006) 252–261 www.elsevier.com/locate/anucene
Reactivity insertion limits in a typical pool-type research reactor cooled by natural circulation H. Kazeminejad
*
Gamma Irradiation Center, Atomic Energy Organization of Iran, P.O. Box 11365-3486, Tehran, I.R. Iran Received 29 August 2005; received in revised form 10 October 2005; accepted 14 November 2005 Available online 27 December 2005
Abstract The course of reactivity insertion in a pool type research reactor, with scram disabled under natural circulation condition is numerically investigated. The analyses were performed by a coupled kinetic–thermal–hydraulic computer code developed specifically for this task. The 10-MW IAEA MTR research reactor was subjected to unprotected reactivity insertion (step and ramp) for both low and highenriched fuel with continuous reactivity feedback due to coolant and fuel temperature effects. In general, it was found that the power, core mass flow rate and clad temperature under fully established natural circulation are higher for high-enriched fuel than for low enriched fuel. This is unlike the case of decay heat removal, where equal clad temperatures are reported for both fuels. The analysis of reactivity represented by the maximum insertion of positive reactivity ($0.73) demonstrated the high inherent safety features of MTR-type research reactor. Even in the case of total excess reactivity without scram, the high reactivity feedbacks of fuel and moderator temperatures limit the power excursion and avoid consequently escalation of clad temperature to the level of onset of nucleate boiling and sub-cooled void formation. The code can also be modified to provide an accurate capability for the analyses of research reactor transients under forced convection. Ó 2005 Elsevier Ltd. All rights reserved.
1. Introduction The complex set of physical phenomena that occur in a gravity environment when geometrically distinct heat sink and heat source are connected by a fluid flow path can be identified as natural circulation (NC). In all the adopted configuration of existing water-cooled reactors, NC allows the removal of the decay heat produced by the core. Furthermore, NC is the working mode for steam generators in existing pressurized heavy and light water reactors. It is necessary also for the core cooling system in an unlikely loss of primary coolant accident to occur in all the reactor types. The new frontier of nuclear technology represented by advanced reactor designs, with the aim of simplifying the systems, of reducing costs and of increasing the safety level, *
Tel.: +98 21 88004065; fax: +98 21 88009054. E-mail address: [email protected].
0306-4549/$ - see front matter Ó 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.anucene.2005.11.002
foresees an extended use of the NC principles. This includes, in some cases, the nominal operation of the reactors and, in all cases, the reactor cooling by passive systems (Jafari et al., 2003). Although the application of passive safety systems is a potential means of achieving simplification and competitive economics in new nuclear power plant designs, the use of passive systems is not entirely new, and is not unique to any particular line of new reactor designs. But an increased reliance on this approach, making safety functions less dependent on active components like pumps and diesel generators, is potentially an important means to achieve reduced costs for future nuclear power plants. Many new reactor designs incorporate passive systems based on NC. The study of NC in research reactors has up to now been mainly focused on decay heat removal from the core (Langerman, 1994; Woodruff et al., 1996; Housiadas, 2000; Hamidouche et al., 2004) and no consideration is given to the issue of reactivity insertion under NC. The main objec-
H. Kazeminejad / Annals of Nuclear Energy 33 (2006) 252–261
253
Nomenclature a cp Ci De f fa fe fr F(t) g G hs hg k L m N p P q00i q00l q000 q000 q000 c r r0 R
constant for planar or cylindrical geometry specific heat ith group precursor density equivalent diameter friction factor axial peaking factor engineering peaking factor radial peaking factor external neutron source gravitational acceleration defined in Eq. (3) surface heat transfer coefficient defined in Eq. (6) thermal conductivity fuel length time step number of collocation points reactor power pressure incipient boiling heat flux coolant heat flux volumetric heat generation average volumetric heat generation volumetric heat generation at the fuel center radial coordinate dimensionless radial coordinate (r/R) fuel radius, fuel half thickness
tive of this study is to establish the state of the art in the area of NC by considering the reactor start-up and operation utilizing NC mode. Attempts were also made to establish operating limits of the single-phase NC mechanism as applied to nuclear research reactors during a transient step and ramp positive reactivity insertions. Results were obtained for 10-MW IAEA benchmark cores with both high enriched uranium (HEU) and low enriched uranium (LEU) fuel with a failure to scram. An estimate of the reactivity insertion threshold without exceeding the onset of nucleate boiling (ONB) is obtained. Also, the code developed for this purpose can easily be modified and adapted for the testing of methods and for subsequent use in the detail analysis of transient behavior (loss of flow or reactivity insertion under forced flow) in research reactors. 2. Mathematical models 2.1. Heat transfer model Many problems in a nuclear reactor analyses require the knowledge of the temperature distribution in the fuel, cladding and the coolant during normal reactor transients and
Re t T V w Yc z
Reynolds number time temperature velocity water channel thickness clad thickness axial coordinate
Greek letters b total delayed neutron fraction bi ith group of delayed neutron fraction bl coolant volumetric thermal expansion coefficient h defined in Eq. (3) ki ith group decay constant K prompt neutron life time l dynamic viscosity q reactivity or fluid density s reactor period Subscripts b bulk condition f fuel c clad l coolant o initial or reference conditions
in accident situations such as reactivity insertion accidents (RIA). Temperature distribution is an important phenomenon, which influences the fuel behavior and integrity during all stages of the facility lifetime. This information is required to assess the safety margins in fuel plates both during normal reactor transients and when exposed to hypothetical reactor accident condition. The behavior of the cladding is also very important, because it is the first barrier to release of fission products. The exact calculation of temperature variation in fuel and cladding is rather complex. The complexity arises due to the fact that the heat generation in fuel assembly, the temperature variation and thermal–hydraulic behavior of the core are interdependent. The problem is very simplified by considering a single coolant flow channel surrounded by two fuel plates and a heat transfer coefficient between the cladding and the coolant. The coolant temperature in the reactor core is taken as the boundary condition, which supplies the fuel assembly with the sink temperature. The heat generation in the fuel is assumed to have a cosine distribution or to be a known function of time. The thermal behavior of the fuel is described by heat conduction equation in fuel and clad-
254
H. Kazeminejad / Annals of Nuclear Energy 33 (2006) 252–261
ding. One-dimensional equations for mass, momentum and energy are solved in the coolant. Because of the complexity of equations, boundary conditions and temperature dependence of physical parameters, an analytical solution cannot be obtained for the differential equations. The problem must be approximated by numerical solutions. The heat transfer model and numerical procedure used in developing a general computer code suitable for research reactor are similar to thermal–hydraulic codes for power reactors (COBRA-IV, 1976) in that the transient thermal response of the fuel is coupled with the hydraulics through the heat transfer coefficients. The fuel plates heat conduction model includes axial conduction and variable thermal conductivity model. The accuracy of the present results is improved by using an orthogonal collocation numerical procedure. The physical properties of the coolant are calculated by interpolation of the liquid property table. As the fluid is heated in each computational grid, this process is repeated, employing linear interpolation where necessary, to determine the fluid properties in terms of the increasing temperature. The fuel heat transfer model allows the computation of the fuel plate heat fluxes and internal temperature distributions based on the heat generation in the fuel and the surface heat transfer coefficient. The inclusion of axial conduction and the temperature dependence of fuel thermal conductivity provide greater detail and allow consideration of wider class of problems. The method of solution of the governing differential equations is based on the method of weighted residuals (MWRs) as a replacement to the finite difference scheme in the radial direction. Finite difference scheme is still used for time derivatives and axial, z, space derivatives. In MWR, orthogonal collocation is employed to determine the form of the approximated polynomial solution. MWR affords a higher order of accuracy by using the roots of orthogonal polynomial as the nodal position where the solution is evaluated, which eventually provides accuracy comparable to finite element techniques. In what follows a brief description of the governing equations, method of solution and the results obtained from its main application will be presented and discussed.
8.05
0.038 0.223
0.051
6.30
6.65
7.60
0.475
Standard Fuel Element: 23 fuel plates Control Element: 17 Fuel Plates 4 Al Plates Fig. 1. View of the unit cell representing the 10-MW benchmark MTR fuel element (dimensions are in cm).
w
Clad L
Fuel Coolant
Z
R
V, To, P
2.1.1. Fuel interior The view of a unit cell representing the 10-MW benchmark MTR fuel element is shown in Fig. 1. A single channel geometry as shown in Fig. 2 is chosen to represent the behavior of any region of the core. The z direction is the direction of flow along the channel and the R direction is the direction normal to the plate. The heat transfer model can be derived from the fundamental heat conduction equation: oT qcp ¼ r ðkrT Þ þ q000 ; ot
ð1Þ
where k is the local thermal conductivity and q000 is the volumetric rate of heat generation.
Fig. 2. Flow model.
Eq. (1) is expanded in two dimensions for both planner and cylindrical geometry with the radial coordinate nondimensionalized using the relationship r = r 0 /R, where r 0 is the radial coordinate and R is the fuel half thickness or radius: oT 1 o oT o oT ¼ 2 a1 ra1 k f ðT Þ þ k f ðT Þ þ q000 ; ð2Þ ot or ot oz oz Rr where a = 1 for planar; a = 2 for cylindrical. Eq. (2) is reduced to a linear partial differential equation using KirchoffÕs transformation:
ðqcp Þf
H. Kazeminejad / Annals of Nuclear Energy 33 (2006) 252–261
1 h¼ ko
Z
T
kðT Þ dT ¼ GðT Þ;
ð3Þ
To
where ko is the conductivity at reference temperature To. Differentiating Eq. (3) with respect to r and t and substituting into the differential equation (2), yields the transformed equation: ðqcp Þf
k o oh k o o a1 oh o oT ¼ r þ k f ðT Þ þ q000 . or oz oz k f ðT Þ ot R2 ra1 or
ð4Þ
Eq. (4) can be solved numerically by a combination of techniques. In the orthogonal collocation method, the radial positions are taken to be the roots to orthogonal polynomials, which affords accuracy comparable to finite element techniques. The radial positions or collocation points are defined by Finlayson (1972) for a second-order approximation (N = 3). The axial conduction term is calculated by central finite difference and the time derivative is approximated by a forward finite difference. The boundary condition at the fuel surface is handled by a lumped resistance technique. The governing equation is k o oh ¼ hg ðT N T N þ1 Þ; ð5Þ R or where N and N + 1 denote the fuel and clad exterior surfaces, respectively. The term hg is defined by 1 1 Yc ¼ þ ; hg hc k c
ð6Þ
where Yc and kc are cladding thickness and conductivity and hc is the fuel-clad gap conductance if any. 2.1.2. Cladding A transient energy balance for the lumped clad is: a1 oT N þ1 hg rN hs ¼ ðqcp Þc ðT N T N þ1 Þ ðT N þ1 T l Þ ot Y c rN þ1 Yc 2
þ k c ðT Þ
o T N þ1 ; oz2
ð7Þ
where hs is the clad surface heat transfer coefficient and Tl is the local fluid temperature. To evaluate the implicit temperatures TN and TN + 1 in Eqs. (5) and (7), a truncated Taylor series is used: T ¼ Tm
GðT m Þ hðT m Þ ; G0 ðT m Þ
ð8Þ
where the function G is defined in Eq. (3) and G 0 is the derivative of G with respect to T. After applying the orthogonal collocation method to boundary equations (5) and (7), they are combined with the differential equation (2) to yield the matrix equation of the form [A]{h} = {Q} for the transformed temperatures at one axial level. This equation is solved by an iterative Gauss-Seidel procedure. Once the {h} solution vector is determined, the temperature solution is evaluated using Eq. (8) with h(Tm) specified. If any elements of [A] or {Q} are temperature dependent, the procedure is repeated
255
until the resulting temperature change remains within a set limit and the solution is considered converged. The volumetric rate of heat generation q000 and fluid temperature, Tl at each axial node are taken explicitly from the previous time step solutions of the kinetics and thermal–hydraulics, respectively. The sequence is repeated at each axial level to the ends of the fuel with additional axial sweeps if axial conduction is being considered. Normally, for plate type fuel elements the axial conduction has little effect and only two or three additional iterations are necessary. 2.2. Thermal–hydraulics model Thermal–hydraulic modeling is based on the equations of conservation of mass, momentum and energy, formulated for a one-dimensional axial homogeneous flow through a vertical channel (Fig. 2) proposed by Nakamura (1977): Continuity equation: oq oðqV Þ þ ¼ 0. ot oz
ð9Þ
Using Boussinesq approximation, the differentials in density in the continuity equation (9) are of the order of b and hence can be neglected, to give oV/oz = 0.0 as for an incompressible fluid. Langerman (1994) show that, for decay heat removal by natural circulation in reactor fuel plates, the difference between the temperature ratios calculated with the Boussinessq and the non-Boussinesq models is rather small for water flow. z-momentum equation (vertical direction): oV oV 1 oP fV 2 þV ¼ þ gbl ðT l T o Þ. ot oz qo oz 2De
ð10Þ
Energy equation: oT l oT l 2hs ðT c T l Þ þV ¼ . qcp w ot oz
ð11Þ
The Boussinesq approximation is incorporated into the momentum equation (10). It assumes that the liquid density q varies with temperature Tl according to q = qo[1 bl(Tl To)], but all other liquid properties are invariant within the specified temperature range. qo is the density at reference temperature To and bl is the coolant volumetric thermal expansion coefficient. The coupled differential equations (9)–(11) with initial conditions V(0) = 0 and Tl(0) = To, were solved by the Lax-Wendroff explicit method (Anderson et al., 1984), which is a second-order finite difference method, for time evolution of liquid velocity, V, and liquid temperature Tl. 2.3. Reactor kinetics model The reactor power dynamics under transient conditions is derived through a numerical solution of the conventional point kinetics formula. The point kinetics equations for six
256
H. Kazeminejad / Annals of Nuclear Energy 33 (2006) 252–261
delayed groups and light water moderator are (Lamarch, 1982): 6 X dp qðtÞ b ¼ pþ ki C i þ F ðtÞ; ð12Þ dt K i¼1 dC i bi ¼ p ki C i i ¼ 1; 2; . . . ; 6; ð13Þ dt K where p is the reactor power, Ci is the ith precursor density, q(t) is the time dependent P reactivity function, bi is the ith 6 delayed fraction, and b ¼ i¼1 bi is the total delayed fraction. In addition, K is the neutron generation time, ki is the ith group decay constant, and F(t) is the time dependent neutron source function. The appropriate initial conditions that accompany Eqs. (12) and (13) are: b pð0Þ ¼ po ; C i ð0Þ ¼ i po . ð14Þ ki K Eqs. (12) and (13) are a system of coupled linear ordinary differential equations. The time dependent parameters in this system are the reactivity and the neutron source terms. An important property of the point kinetics equations is the stiffness of the system. It is well known that stiffness is a severe problem in numerical solution of the point kinetics equations and results in the need for small time steps in a computational scheme. There has been a great deal of research done in the field of reactor kinetics. Much of the work has been focused on eliminating the aforementioned problem of stiffness. A modified Runge–Kutta method to avoid the problem of stiffness devised and implemented by Sanchez (1989) is used to derive the solution of the set of differential equations. This method constitutes an easy to implement algorithm that provides results with sufficient accuracy for most applications. The main advantage of this method is that it allows step size control to be performed in a systematic manner and estimation of the truncation error can be readily obtained at each time step. The reactivity is estimated by taking into account the different feedback effects as follows (Housiadas, 2002): qðtÞ ¼ qex ðtÞ þ aM dql aT dT l aF dT f .
ð15Þ
In Eq. (15), the first term, qex(t), is the externally introduced reactivity (for example, the step and ramp reactivity insertions or shutdown reactivity if scram is enabled) and the remaining terms are the various feedback contributions. As can be seen three feedback mechanisms are taken into consideration, change of coolant density dql, change of coolant temperature dTl (spectrum effects only) and changes of fuel temperature dTF (Doppler effects). The corresponding constant coefficients of reactivity are, respectively, aM, aT and aF. In applying Eq. (15) the changes in feedback reactivity have to be determined by considering the variations of the mean (whole volume) coolant and fuel temperatures. 2.4. Heat generation in fuel plates Axial variation of the volumetric heat generation is assumed to have a cosine distribution or to be given as any function of time, using the exponential representation:
q000 ¼ q000 c cos
pz . L
ð16Þ
In Eq. (16), q000 and q000 c are the axial variation of volumetric heat generation and at the fuel center, respectively. The local average heat generation depends on the axial location z along the fuel element and on the location of the fuel element in the core. The code uses a two-channel model, one channel to represent the hottest plate and flow channel and the other average channel to represent the remaining fuel plates in the core. The axial source distribution was represented by 21 axial regions and a chopped cosine shape which has an axial power peaking of, fa, for both the average channel and the hot channel. For the hot channel, the axial distribution was multiplied by the other specified hot channel factors. Taking into account these variations, the hot channel heat generation at the fuel center is given by q000 q000 fa fr fe ; c ¼
ð17Þ
where q000 is the average volumetric heat generation of the fuel in the core. To compute the hot channel peak temperatures, first an average channel power trace was obtained, and then a second power driven transient run with hot channel factors was carried out. 2.5. Heat transfer and friction factor correlations Surface heat transfer coefficient is required to couple the hydraulic solutions to the fuel and clad model. Heat transfer from the clad surface to the coolant is described by NewtonÕs Law of Cooling: q00l ¼ hs ðT c T l Þ;
ð18Þ
where ql is the convection heat flux, hs is the convective heat transfer coefficient, Tc and Tl are the local clad and liquid temperatures, respectively. hs will depend upon the properties and flow conditions of the coolant. The coolant flow between the coolant passages of a pool type research reactor core in natural convection is laminar (Reynolds number, Re 6 2000) but the existing heat transfer correlations for natural convection are inadequate and predict the heat transfer coefficients that are too low. Therefore, in line with the single-phase heat transfer model in PARET, an approximation of the heat transfer coefficient proposed by Woodruff (1984) was used 1=2 kqcP hs ¼ . s
ð19Þ
The heat transfer coefficient predicted by Eq. (19) is higher than the existing heat transfer correlations for natural convection. The use of Eq. (19) requires the selection of an appropriate period, s. This period is chosen as the minimum sm such that |sm + 1 sm|/sm is less than 0.01 where sm is the instantaneous reactor period at time step m. This has proven to be a reasonable criterion under most conditions. At certain pressure and temperature conditions, an incipient heat flux, qi, exists at which isolated vapor nucle-
H. Kazeminejad / Annals of Nuclear Energy 33 (2006) 252–261
ation occurs on the cladding surface. The Bergles–Rohsenow correlation is used to predict qi as follows (Collier, 1980): q00i ¼ 15:60P 1:156 DT as
a ¼ 2:30P 0:0234 ;
ð20Þ
where P is the absolute pressure (kPa), DTs is the degree of superheat (K) and qi is the incipient heat flux. Onset of nucleate boiling (ONB) is not a limiting criterion in the design of a fuel element. However, it is a heat transfer regime, which should be identified for proper hydraulic and heat transfer considerations, i.e., single phase flow versus two-phase flow. The actual axial location at which ONB will occur depends upon the axial heat flux distribution, the coolant velocity and the pressure drop along the channel. For simplicity, the heat flux for ONB can be calculated conservatively by using the worst combination of parameters, i.e., ONB occurs at the channel exit with peak heat flux, lowest pressure and saturation temperature, and highest coolant temperature rise. The friction factor for laminar flow between rectangular cross-sections is determined as f = 96/Re(lc/lb), where Re is the Reynolds number and subscripts b and c refer to coolant properties evaluated at the bulk and clad temperatures, respectively. 3. Benchmark problem description The benchmark problem considered in this investigation consists of a 10-MW light water, pool type research reactor specified by the IAEA (1992). Two similar core configurations with different fuel enrichments HEU (93% of U235) and LEU (20% of U235) and fuel material type but with the same geometric dimensions (Fig. 1) and operating conditions are considered. These cores have been used as benchmark for the transient analyses of research reactors. The detailed specifications are given in Table 1. The benchmark problem is performed for transients governed by the core kinetics and thermal–hydraulics. 4. Results and discussions As stated earlier, the main aim of the present work is to study the safety aspects of 10-MW benchmark MTR research reactor under NC condition. The safety of the reactor demands that the stagnant liquid inside the core does not reach its boiling point. Therefore, it is important to find the level of reactivity insertion where the liquid starts boiling in a research reactor under NC condition. The transients performed for the benchmark cores are self-limiting and the limits imposed on reactivity insertion by the ONB condition. Figs. 3a–c provide a comparison of power, reactivity, core mass flow rate and peak clad temperature for both HEU and LEU for a step reactivity insertion of $0.5. Fig. 3a shows the development of core average power and reactivity during the initial start up phase. The situa-
257
Table 1 IAEA research reactor benchmark problem operating conditions Core material Nuclear fuel Fuel element Coolant Moderator
MTR Plate type clad in AL Light water Graphite-light water
Core thermal–hydraulics Fuel thermal conductivity (W/cm K) Cladding thermal conductivity (W/cm K) Fuel specific heat (J/g K) Cladding specific heat (J/g K) Fuel density (g/cm3) Cladding density (g/cm3) Radial peaking factor Axial peaking factor Engineering peaking factor Inlet coolant temperature (°C) Operating pressure (bar) Fuel element dimensions Length (cm) Width (cm) Height (cm) Number of plates/fuel ECS/ECN Plate meat (mm) Width (cm) active/total Height (cm) Water channel thickness (mm) Plate clad thickness (mm)
8.00 7.60 60.0 23/17 0.51 6.30/6.65 60.0 2.23 0.38
Core kinetics Effective delayed neutron fraction Prompt neutron generation time (ls) Density feedback coefficient, aM ($/% void) Doppler feedback coefficient, af ($/°C) Coolant feedback coefficient, aT ($/°C)
HEU 1.58
LEU 0.5
1.80
1.80
0.728 0.892 0.68 2.7 1.4 1.5 1.2 38.0 1.7
0.340 0.892 4.45 2.7
HEU 7.607E 3
LEU 7.275E 3
55.96
43.74
0.3257
0.4047
3.6E 5
3.31E 3
1.537E 2
1.082E 2
tion simulates the reactor start up when withdrawing the control rods to increase the reactor power. After a short coolant heat up time of about 50 s, the core mass flow rate (Fig. 3b) and the peak clad temperature (Fig. 3c) start to increase due to onset of NC. At the same time, the reactivity starts to decrease (Fig. 3a) due to high negative reactivity feedback of fuel and moderator temperatures. The negative feedback effects are able to consume the entire available excess reactivity forcing the power to stabilize and demonstrating the inherent safety of MTR research reactors. For LEU, the power increases to a maximum of 188 kW at 75 s after the start of transient and then attains a nearly constant value of 155 kW after about 140 s. The variation of core mass flow rate is similar to that of power and after about 140 s in transient, it reaches a semi-constant state characterizing a fully developed NC. It is also interesting to note that in a fully developed NC, the net reactivity for both LEU and HEU cores is zero, which
258
H. Kazeminejad / Annals of Nuclear Energy 33 (2006) 252–261
means that the effect of external reactivity is almost balanced by the effects of reactivity feedbacks. In the fully developed NC regime all reactor variables reach a steady state value after about 150 s. For HEU, the peak clad temperature (Fig. 3c) is higher for LEU due to higher power produced by the HEU core despite the fact that, the mass flow rate induced by the HEU core is much higher than that of the LEU core. The peak clad temperatures for both
cores are well below the ONB condition for a reactivity insertion of $0.5. Figs. 4a–c show the time evolution for power, reactivity, core mass flow rate and peak clad temperature through a period of 300 s for a step reactivity insertion of $0.73. This reactivity insertion limit is imposed by the ONB. It should be noted that any slight increase in the value of 0.73 would cause the peak clad temperature for both cores to exceed
Fig. 3. (a) Power and reactivity; (b) mass flux; and (c) peak clad temperature for step reactivity insertion of $0.5.
H. Kazeminejad / Annals of Nuclear Energy 33 (2006) 252–261
the ONB condition and hence making the single-phase fluid flow assumption invalid. As expected, the maximum values attained by the power, core mass flow rate and clad temperature are higher as compared with $0.5 reactivity insertion. For LEU and HEU cores, the maximum powers are 518 kW and 593.5 kW occurring at about 21 s and 21.9 s, with steady-state values of 211 kW and 417 kW, respectively. The peak clad temperatures for LEU and HEU cores are 115.2 °C and 113.3 °C occurring at bout
259
51 s and 44 s, respectively. For both LEU and HEU cores, the clad temperature first decreases to a minimum value of 98 °C and 105.5 °C at about 87 s and 60.5 s before reaching a steady-state value of 108 °C and 115.6 °C, when NC is fully established. The peak clad temperature for both cores remain below the ONB condition of 116.5 °C during the entire transient, which is also shown in Fig. 4c. Figs. 5a–c provide a comparison of HEU and LEU benchmark cores for a ramp positive reactivity insertion
Fig. 4. (a) Power and reactivity; (b) mass flux; and (c) peak clad temperature for step reactivity insertion of $0.73.
260
H. Kazeminejad / Annals of Nuclear Energy 33 (2006) 252–261
of $0.73/0.5 s. The variation of core power, reactivity, core mass flow rate and peak clad temperature are almost similar to those of step reactivity insertion. This is due to the fact that the ramp insertions of short duration are equivalent to a step insertion. The entire ramp is inserted before the power, feedback and temperatures have increased substantially and the limiting reactivity insertion remains constant. Although the results for reactivity insertion of $0.73 show differences between the HEU and LEU cores at the
early stages of the transient, the peak clad temperatures reached in the fully developed NC for both cores remain below the ONB condition, with the HEU core attaining a higher peak clad temperature than that of the equivalent LEU core. This is unlike the case for decay heat removal under NC, where the reported peak clad temperatures for both cores are the same. The results for the analyses of loss of flow transients for an MTR research reactor benchmark problem using PARET and RELAP5/3.2 mode are
Fig. 5. (a) Power and reactivity; (b) mass flux; and (c) peak clad temperature for ramp reactivity insertion of $0.73/0.5 s.
H. Kazeminejad / Annals of Nuclear Energy 33 (2006) 252–261
reported in the literature by Woodruff et al. (1996) and Hamidouche et al. (2004). In the flow decay, due to the combined effect of constant decay heat and continuous reduction of the core flow rate, the peak-clad temperature exhibits a second rise after the scram of the reactor power. The increase is further sustained as the flow regime passes to laminar regime and to mix flow when the natural convection valve opens. The peak-clad temperature begins to decrease only when the NC flow is fully established. Indeed, this phenomenon is governed by the balance between the pressure and the buoyancy forces involved during the flow decay process. For both HEU and LEU, the peak-clad temperatures for fully established NC predicted by PARET and RELAP codes in slow and fast flow decay processes are the same. 5. Conclusions The results for a simple physical model developed to simulate thermal–hydraulic and neutron dynamic features for a 10-MW IAEA research reactor under NC mode are presented. Two cores with similar configuration and operation but different fuel enrichments (HEU and LEU) are considered. The analyses were carried out for self limiting cases and the limits imposed on the reactivity insertion by the ONB clad surface temperature. Based on model predictions, it can be concluded that a research reactor with either HEU or LEU core can be safely operated under NC mode for a reactivity insertion as high as $0.73 without exceeding the ONB condition. The results for step reactivity insertion were found to be similar to those of ramp reactivity insertion with short duration. For the estimated reactivity insertion, the peak clad temperatures for both HEU and LEU cores under fully established NC attained a constant value with the HEU core predicting a higher peak clad temperature (115.6 °C) than that of the equivalent LEU (108 °C) core. However, for decay heat removal under NC condition, the peak-clad temperatures reported in the literature for both cores are almost equal. The high negative reactivity feedbacks of fuel and moderator temperature were able to consume the entire available excess reactivity, limiting the power excursion and forcing the reactor to stabilize during the reactor start-up in NC mode. The results presented clearly demonstrated the inherent safety features of MTR research reactors utilizing NC for the recommended reactivity insertion. Also, the use of relatively simple physical models were shown to provide a
261
useful capability for the analyses of most research reactors encountered in practice under transient conditions. Acknowledgements This work is supported by the research grant from Atomic Energy Organization of I.R. Iran (AEOI) whose assistance is hereby gratefully acknowledged. References Anderson, D.A., Tannehill, J.C., Pletcher, R.H., 1984. Computational fluid mechanics and heat transferSeries in Computational Methods in Mechanics and Thermal Sciences. Hemisphere, New York. COBRA-IV, 1976. An interim version of COBRA for thermal–hydraulic analysis of rod bundle nuclear fuel elements and cores, BNWL-1962. Collier, J.G., 1980. Convective Boiling and Condensation. McGraw-Hill, New York. Finlayson, B., 1972. The method of weighted residuals and variational principlesMathematics in Science and Engineering Series, vol. 87. Academic Press, New York. Hamidouche, T., Bousbia-Salah, A., Adorni, M., DÕAuria, F., 2004. Dynamic calculations of the IAEA safety MTR research reactor benchmark problem using RELAP5/3.2 code. Annals of Nuclear Energy 31, 1385–1402. Housiadas, C., 2000. Simulation of loss of flow transients in research reactors. Annals of Nuclear Energy 27, 1683–1693. Housiadas, C., 2002. Lumped parameters analysis of coupled kinetics and thermal–hydraulics for small reactors. Annals of Nuclear Energy 29, 1315–1325. IAEA-TECDOC-643, 1992. Research reactor core conversion guide book, vol. 3, Analytical Verification. Jafari, J., DÕAuria, F., Kazeminejad, H., Davilu, H., 2003. Reliability evaluation of a natural circulation system. Nuclear Engineering and Design 224, 79–104. Lamarch, J.R., 1982. Introduction to Nuclear Engineering, second ed. Addison-Wiley, Reading, MA. Langerman, M.A., 1994. Calculating decay heat removal rates via natural circulation along variable heat flux reactor fuel plates. Nuclear Engineering and Design 150, 61–68. Nakamura, S., 1977. Computational Methods in Engineering and Science with Applications to Fluid Dynamics and Nuclear Systems. Wiley– Interscience, New York, pp. 171–175. Sanchez, J., 1989. On the numerical solution of the point reactor kinetics equations by generalized Runge–Kutta method. Nuclear Science and Engineering 103, 94–99. Woodruff, W.L., 1984. A kinetics and thermal–hydraulics capability for the analysis of research reactors. Nuclear Technology 64, 196–206. Woodruff, W.L., Hanan, N.A., Smith, R.S., Matos, J.E., 1996. A comparison of the PARET/ANL and RELAP/MOD3 codes for the analysis of IAEA benchmark transients. In: International Meeting on Reduced Enrichment for Research Reactors, Seoul, Republic of Korea, October 7–10.