- Email: [email protected]
Development of a thermal-hydraulic analysis code for CARR
annals of NUCLEAR ENERGY Annals of Nuclear Energy 32 (2005) 261–279 www.elsevier.com/locate/anucene
Development of a thermal-hydraulic analysis code for CARR W.X. Tian a, S.Z. Qiu a
a,b,* ,
Y. Guo a, G.H. Su
a,b
, D.N. Jia
a,b
State key laboratory of Multi Phase Flow in Power Engineering, XiÕan Jiaotong University, XiÕan, Shaanxi 710049, PR China b Department of Nuclear and Thermal Power Engineering, XiÕan Jiaotong University, XiÕan, Shaanxi 710049, PR China Received 3 July 2004; accepted 13 September 2004
Abstract The China advanced research reactor (CARR) being built in Beijing, China, is a multipurpose research reactor for a variety of fields. Theoretical calculation of thermal hydraulic characteristics of CARR is presented in this paper. The theoretical analysis consists of initial steady and transient accidental analyses. Point reactor neutron kinetics model with six groups of delayed neutron is adopted for the solution of reactor power. All possible flow and heat transfer conditions are considered and the corresponding optional models are supplied in the theoretical calculations. A new simple and convenient model is proposed for the resolution of the transient behaviors of main pump instead of the complicated four-quadrant model. Gear method and Adams predictor–corrector method are adopted alternately for a better solution to such ill-conditioned differential equations corresponding to detail process. The initial multi-channel analysis shows that the effects of geometrical size on flow distribution play dominant role and the effects of core power distribution may be neglected. The temperature fields of fuel elements under asymmetrical cooling condition are also obtained, which are the bases for further study on transient-induced stress analysis, etc. Accidental analyses show that the activity of emergency cooling system apparently reduces the peak temperatures of fuel and coolant, peak quality and other operation parameters. Thus it effectively ensures the
*
Corresponding author. E-mail addresses: [email protected] (W.X. Tian), [email protected] (S.Z. Qiu).
0306-4549/$ - see front matter Ó 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.anucene.2004.09.003
262
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
safety in operation of CARR. Because of the adoption of modular programming techniques, this code is expected to be applied to accidental analysis of other types of reactors by easily modifying the corresponding function modules. Also, this code is expected to be validated against experimental data. Ó 2004 Elsevier Ltd. All rights reserved.
Nomenclature Cp De f g Gr h Hp I K Nu Pr q Re DT W
specific heat, J/kg k equivalent diameter, m friction factor, dimensionless gravity acceleration, m/s2 Grashof number, dimensionless heat transfer coefficient, W/m2 k pump head, m inertia of the rotor, kg m2 thermal conductivity, W/m k Nusselt number, dimensionless Prandtl number, dimensionless heat flux, w/m2 Reynolds number, dimensionless temperature difference, K mass flow rate, kg/s
Greek symbols / dimensionless length shown in Fig. 11 du,1 length shown in Fig. 11 du,2 length shown in Fig. 11 e roughness, m b expansion coefficient q density, kg/m3 x rotation speed, rps g running efficiency Superscript/subscripts b main fluid flow g steam/gas o normal working condition vf vapor/steam flow w wall
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
263
1. Introduction The China Advanced Research Reactor (CARR) is being built at the China Institute of Atomic Energy (CIAE, Beijing, China) to satisfy the multipurpose for a variety of fields in neutron scattering measurements, radioisotope production, neutron transmutation doping of silicon, etc. The CARR is a tank-in-pool, inversed neutron-trap-type reactor with a nuclear power of 60 MW. Slightly pressurized light water is used as the primary coolant and moderator, and heavy water as the reflector. The systematic diagram of CARR is shown as Fig. 1. The reactor body is immersed in a pool of 16 m in depth and the core is located 12 m below the pool water surface. The core is about 0.85 m in height and 0.451 m in diameter. The core is surrounded
Fig. 1. Schematic of CARR complex structure.
264
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
by a heavy water reflector tank, in which the under moderated neutrons leaking out through the core surface would be further moderated and the maximum thermal neutron flux is in a peak. Under CARRs normal working condition (full power working condition), the natural circulation crack valve (NCCV) is kept to be closed. The coolant is pumped to flow through the cold leg and the flow guiding tank, downward through the active zone of the core, then through the decay tank, the hot leg, the heat exchanger and recirculated to the main pump. The fluid flow path under normal working condition is denoted in solid arrow is shown in Fig. 2. However, if the CARR is shutdown, the NCCV will be opened. Due to the driven force caused by the density difference between the fluid in the water pool and the coolant in reactor vessel, the natural circulation (NC) will occur, by which the residual heat will be removed. However, the coolant flow direction will inverse during the change from forced circulation (FC) to NC, as shown in Fig. 2. Thus, the coolant flow path is: reactor pool-filter-decay tank-core vessel-flow guiding tank-NCCV-reactor pool. At a certain time in process of establishment of NC after shutdown, the mass flow rate in the core is becoming zero that is very dangerous for CARR. Therefore, safety analysis is very necessary to be carried out to ensure that CARR is in safety. The objective of the present work is to analyze the steady as well as transient thermal-hydraulic behaviors of CARR. For the initial steady analysis, the multi-channel
Hot water layer Reactor pool
Main Pump
Emergency Coolant system
Water
Flow guiding tank
Natural circulation
Nccv
Normal condition
Plate Heat Exchanger
Fluid flow Direction Under Normal Condition Fluid flow Direction Under Natural Circulation Condition
Active core
Filter screen
Decay tank
Fig. 2. Schematic diagram of the CARR.
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
265
analysis is conducted to obtain the detail flow distribution and temperature field of each fuel element. For the transient analysis, the behaviors of both the mean channel (MC) and the hot channel (HC) are obtained and discussed. All analyses are aimed at verifying if the CARR is safety in operation.
2. Theoretical model There are mass, energy and momentum conservation equations in the theoretical models. These equations are mutual coupled and are easy to find in reference (Guo et al., 1997; Su et al., 2001). The characteristics of the present theoretical model are introduced in the followings. 2.1. Physical models 2.1.1. Core power model Point reactor neutron kinetics model with six groups of delayed neutron is adopted for the solution of the power of the core. The reactivity feedback caused by the temperature changes of the moderator and the fuel is specially considered. 2.1.2. Heat transfer and flow friction models The main flow and heat transfer conditions might occur in the transient process are considered and the corresponding optional models are supplied as shown in Tables 1 and 2. The detail expressions of these correlations can be found in reference (Rhosenow and Hartnett, 1973). Linear interpolation has been done to ensure the continuum of heat transfer coefficients in the edge of different flow regimes and thus to prevent the numerical instability caused by the abrupt change of heat transfer coefficients. Chisholm correlation (Chisholm, 1967) is adopted For the calculation of two phase frictional multiplier. 2.1.3. Critical heat flux model Critical heat flux (CHF) is the prominent factor restricting the core power. There are more than 600 empirical or semi-empirical correlations for CHF based on Table 1 Flow regimes and heat transfer correlations Flow/heat transfer regime Single phase water Laminar flow Transition flow Turbulent flow Subcooled boiling Saturated boiling Transition/film boiling Single vapor flow
Heat transfer correlation Nu ¼ 0:17Re0:33 Pr0:43 Linear interpolation 0:11 lb lw 0.25 P/6.2
Pr Nu ¼ fRe 8X
Pr Prw
0:25
Gr0:1 (Collier, 1981)
(Rhosenow and Hartnett, 1973)
(Jens and Lottes, 1951) DT = 25q e Chen correlation, Expression omitted (Chen, 1966) AECL lookup Tablei Groeneveld h (2003) h ¼ 0:13k
q2 gbðT w T g Þ 1=3 l2
Cp l k
(Thomas, 1992)
266
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
Table 2 Flow regimes and friction correlations Flow regimes
Darcy friction coefficient correlations
Laminar flow Transition flow Turbulent flow Cole–Brook
f = C/Re different values for C for various geometrical channels f = 0.048 e 2:51 p1ffiffi ¼ 2 log e=D p ffiffi þ 3:71 f
Re
f
channels with various structures. Each correlation has its own very limited scope of application. It may be very difficult to find a correlation covering all flow conditions might occur during the establishment of NC of CARR. Thus, in the present study, three correlations are adopted and each one aims at different working condition. W-3 equation has been widely used in the design of pressurized power reactor and is adopted in the existing popular code such as RETRAN series, etc. Sudo correlation (Sudo and Kaminaga, 1993) is proposed on the basis of a great deal of CHF experimental data on rectangular channels. Compared with other CHF correlations, Sudo equation gives a relative conservative prediction. It was adopted in the design of the JRR-3M (Kanminaga, 1990) and its applicability under middle and high flow condition has been verified. Another CHF correlation, Kutateladze correlation (Thomas, 1992) is verified to be valid under the condition of pool boiling and low mass flow rate.
2.1.4. Equipment model Four-quadrant characteristic model has been widely adopted in universal commercial software such as RETRAN series for a better calculation of main pump transient performance. However, this model is inconvenient for use because of its strict and rigorous request for the detail pump performance parameters that are very difficult to obtain even for an actual certain pump (Guo, 1995). In this study, a new and convenient model on the computation of pump head is proposed by quoting the similarity law. The detail parameters are not necessary in the new model. Under all working conditions, the following relationship is always available among the pump output power Np, the pump lift Hp and the mass flow rate W N p ¼ W H p g:
ð1Þ
While the main pump being operated under switched-off condition, the effective output pump power equals to the decrease of the inertial kinetic energy of the rotor of the pump, that is N p ¼ Ixg
dx : dt
ð2Þ
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
267
For the same pump, according to the similarity law W x ¼ : ð3Þ W o xo Then, from Eqs. (1)–(3), it can be draw that 2 I g xo dW : ð4Þ Hp ¼ dt g Wo Eq. (4) is finally employed to calculate the transient behaviors of the main pump. It should be pointed out that there are four parallel identical main pumps arranged in the primary loop of CARR. The total mass flow rate equals the sum of the four branches and the equivalent pump head equals to that of any one of the four parallel pumps. The emergency coolant system (ECS) is connected with the CARR main structure at the cold leg as shown in Fig. 2. Under normal working condition, the ECS works as the cooling system of the pool water. When the mass flow rate in the primary loop is less than 83.3 kg/s, the NCCV will be automatically opened, the ECS is brought in service and the pool water is injected into the core. There are four plate heat exchangers with very complex structure placed in the primary loop of CARR. The fluid flow inside the exchanger is also very complex and disordered. Logarithmic temperature difference method (Tao, 1996) is adopted to calculate the performance of the heat exchangers. 2.1.5. Coolant thermophysical property model The correct coolant thermophysical property model is the key of success for a code. The thermophysical property correlations of water and vapor used in this study are chosen from the international standard (Rhosenow and Hartnett, 1973). Comparison between the model prediction and look-up table is conducted to verify and extend the applicability of those models under low pressure. 2.2. Nodalization The main principle of nodalization is that the different geometrical structure and physical boundary conditions of the core and pipes will be divided into different nodes. Two nodalization systems were adopted to meet the fact that the flow direction and flow path under normal working condition are different from those of NC condition. For normal working condition, the system node diagram is shown in Fig. 3. The active zone is modeled by MC and HC including the fuel, the clad and the coolant which consists 17 control volumes, respectively. For NC condition, during the phase of residual heat removal, the system node diagram is shown in Fig. 4. It is similar to that of normal working condition however, the flow direction, the flow path as well as the numeral order of volumes changed. The number of nodes of the core, the exchanger and pipes can be changed according to the calculation requirement.
268
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279 115-119
72 73 74
18 19 20
89 90 91
Plate heat exchanger
Hot channel
Core
92 93
Fuel
Clad
35 36 37
75 76
Coolant
Clad
52 53 54
3
58 59
4 5
Coolant
21 22
Fuel
Mean channel
2 38 39
Cold leg
1 Upper plenum
106 107 108
122-124
119-121
Main pump
55 Lower plenum 56 57 109-114,125-132
Flow direction
Hot leg
Fig. 3. The systematic diagram of nodalization of CARR under normal working condition.
Pipes 114-117
62
28 27 26
11 10 9
7
Coolant
core
65 64 63
Lower plenum
96 95
113 112
82 81 80
99 98 97
NCCV
Hot channel
Coolant
45 44 43
79 78
Fuel
25 24
59 58
61
Clad
42 41
Clad
Upper plenum
Fuel
Mean channel
60
118
Rector pool 119
8
6
Flow direction
Pipes: 1-5
Fig. 4. The systematic diagram of nodalization of CARR under natural circulation condition.
2.3. Calculation procedure The theoretical models are applied to every control volume. Thus, a set of equations with a stiff differential equation will be built. To solve it, a code has been made. The program flow diagram is shown in Fig. 5. The code is composed of several large modules including Data input module, Initialization module, Control system module, Numerical method module, auxiliary correlation module, Coolant thermophysical property module, Output modules, etc. Each module may be easily modified to
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
269
BEGIN
Data input
Initialization
Choice 1. Normal working condition that fluid flow downward through the core just as force convection
Normal/steady condition calculation
Time step setup etc.
working condition
Choice 1
Choice 2
Choice 2. The emergency cooling condition Choice 3
N Transient calculation
Simulation is over or not?
Y
Choice 3. NC working condition that fluid flow upward through the core
Output END
Fig. 5. The Program flow diagram.
be competent for the calculation under different working conditions of CARR or even other kinds of reactors. 2.4. Introduction of numerical method Once the nodalization of a nuclear reactor system was obtained, the transient analysis can be reduced to the solvement of a set of ordinary differential equations with the initial conditions as shown in Eqs. (5) and (6), where y is a vector including mass flow, enthalpy and other parameters of each control volume. Because of the speciality and complexity of nuclear reactor system, the corresponding differential equations are usually stiff differential equations (also called ill-conditioned equations)
270
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279 *
d y * * *0 ¼ f ðt; y ; y ; xÞ; dt *
y0: y ðt0 Þ ¼ ~
ð5Þ ð6Þ
Gear method (Gear, 1971; Hind marah, 1974) equipped with Adams predictor– corrector method was adopted for a better solution of stiff differential equations. Gear method is high efficient for the solution of ill-conditioned problems for its good stability, high precision, etc. It is the first order upwind difference method and the code of Gear method has the merit of self-stability and self-comparatively. Adams method required less computation time for its simpler iterative procedure, however, maybe failed for strong ill-conditioned equations. In this study, the two methods are alternatively used according to the stiffness of the set of equations in the whole solving process. For example, in the original stage after shutdown, the core power sharply decreased because of the large negative reactivity caused by shutdown rods, however, the fuel temperature, the coolant enthalpy, etc., might gently vary because of the relative delay of heat transfer. Gear method is thus used in this case. During the residual heat removal process, all parameters slowly change and Adams method is employed to save the running time. It is found in this study the alternate adoption of Gear method and Adams methods greatly reduce the computation time that is important for conducting a long-term simulation.
3. Calculation results and discussions 3.1. Principal design parameters of CARR The principal design parameters of CARR are shown in Table 3. The layout the reactor core of CARR is represented in Fig. 6. The core is composed of 17 standard fuel assemblies (numbering from 1 to 17), 4 following fuel assemblies (numbering from 18 to 21) and other experimental channels. The fuel element is shaped of plate type. All of those standard fuel assemblies have the same geometrical with symmetrical structure as shown in Fig. 7. Each standard assembly is composed of 20 fuel Table 3 Principal design parameters of CARR Core diameter/height (m) Elevation of reactor pool water surface (m) Core inlet/outlet temperature (°C) Core inlet/outlet pressure (°C) Core nuclear power (MW) Core thermal power (MW) Mass flow rate in primary loop (kg/s) Number of standard fuel element/assembly Number of following fuel element/assembly Type of fuel element
0.399/0.85 13.2 35/56.2 0.89/0.62 60 56.4 600 17 4 Plate
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
271
Fig. 6. Layout o the core of CARR.
Fig. 7. Detail structure of the standarad assembly.
plates that format 21 fuel channels. Several marginal channels have different gap sizes as shown in Fig. 7. Especially, channel No. 6 is only a half of the channel formatted by two extreme outside fuel plates of two adjacent assemblies. In the later steady analysis, the total 6 typical channels in each standard assembly will be analyzed separately. 3.2. Steady calculation result under normal working conditions A multi-channel steady analysis is conducted for CARR. Fig. 8 shows the radial distribution (percent) of the core power of each assembly. By simulating the whole reactor core, the detail flow distribution (percent) in the core was obtained as shown in Fig. 9. Very few difference was found among all mass flow percents of standard assemblies in Fig. 9(a), which indicated the influence of the non-uniform core power
272
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
4.86
5.2
4.75
3.1
5.4
3.11
4.96
5.35 5.37
5.9
5.37
5.37
4.95
3.1
5.41 3.11
4.93
4.8
5.15 4.87
4.92
Fig. 8. The percent of core power for each assembly (%).
0
Percent of mass flow
(a)
5
10
15
20
0.0498
0.0498
0.0495
0.0495
0.0492
0.0492
0.0400
0.0400
0.0395
0.0395
0.0390
0.0390
0.0385
0.0385 0
5
10
15
20
Number of assembly
Percent of mass flow distribution
(b)
0.0030
1
2
3
4
5
6
0.0030
0.0025
0.0025
0.0020
0.0020
0.0015
0.0015
0.0010
0.0010 1
2
3
4
5
6
Number of typical fuel channel of assembly No. 9
Fig. 9. Flow distribution. (a) For 21 assemblies of the core. (b) For 6 channels of No. 9 assembly.
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
273
radial distribution on flow distribution may be neglected. Fig. 9(b) represents the flow distribution (percent) of 6 typical channels for assembly No. 9. The strong influence of the channel gap size on the flow distribution is obvious in Fig. 9(b). Thermal-hydraulic analysis of 6 typical fuel elements (from channels No. 1 to No. 6) under unsymmetrical cooling condition (from channels No. 1 to No. 6) is carried out and the detail results of assembly No. 9 are shown in Figs. 10(a)–(d). Fig. 10(a) represents the coolant temperature along the axial flow direction. Coolant in channel No. 6 has the minimum temperature increase because of its maximum mass velocity. It also can be found that the coolant in channels No. 1 and No. 2 has the similar temperature increase because of the same gap size. Fig. 10(b) represents the pressure drop along the axial flow direction. Pressure dropped about 0.21 MPa through the core under normal working condition. Fig. 10(c) represented the maximal fuel temperature along the axial flow direction. Corresponding to Fig. 10(a), the minimum fuel temperature increase occurs in channel No. 6 and the fuel temperature profiles are similar to those in channels No. 1 and No. 2.
0
2
4
6
8
10
12
14
16
18
0.80
0.80
0.75
0.75
0.70
0.70
60
o
Temperature / C
55
50
Pressure / MPa
Ch 1 Ch 2 Ch 3 Ch 4 Ch 5 Ch 6
45
0.60
40
0
2
4
6
(a)
8
10
12
14
0.60
0.55
0
16
2
/ Location of maximal fuel temperature
o
120
Ch 1 Ch 2 Ch 3 Ch 4 Ch 5 Ch 6
80
60
40 2
4
6
8
10
12
Number of axial volume
14
Ch 1 Ch 2 Ch 3 Ch 4 Ch 5 Ch 6
0.504
(d)
8
10
12
14
16
18
0.5 represent the middle of the fuel plate
0.500
0.496
0.492
0.488
2
16
6
Number of axial volume
0.508
100
4
(b)
Number of axial volume
140
Temperature / C
0.65
0.55
35
(c)
Ch 1 Ch 2 Ch 3 Ch 4 Ch 5 Ch 6
0.65
4
6
8
10
12
14
16
Number of axial volume
Fig. 10. Thermal-hydraulic calculation results of 6 typical channels for No. 9 assembly (a) Axial coolant temperature profiles. (b) Axial pressure profiles (c). Axial maximum fuel temperature profiles. (d) Location of maximum fuel temperature for each control volume.
274
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
o Channel N Coolant
Channel N+1 Clad Fuel block
Clad
Coolant
Clad
Clad Tuo Tf,1
Tf , 2
Tci , 2
Tci ,1 Tco ,1
Tco , 2 u ,1
u ,2
o Fig. 11. Sketch map of fuel elements under asymmetrical cooling condition.
Fig. 10(d) represents the radial location of maximum fuel temperature for each control volume. Under symmetrical cooling condition, the fuel temperature reaches its maximum at the middle of fuel plate. In standard assemblies of CARR, for channels No. 1 to No. 6 are under unsymmetrical cooling condition, the maximum fuel temperature does not occur at the middle of each fuel plate as shown in Fig. 11. For the purpose of further study such as stress analysis of the fuel plate, it is very necessary to deduce the abnormal temperature distribution. In Fig. 10(d), the Y-coordinate variable is defined as d ¼ du;1 =ðdu;1 þ du;2 Þ;
ð7Þ
where du,1 and du,2 was shown in Fig. 11. From Fig. 10(d), it may be found that the location of maximum fuel temperature of channel No. 6 has a relative bigger departure from the middle of the fuel plate. However, for channels 1–3, the fuel temperature reaches maximum just at the middle. Observing the Fig. 10(a)–(d), the thermal hydraulic parameters in channel 1 are just similar to those of channel 2. It is verified reasonable and competent to choose 6 typical channels in one assembly to simulate to the whole core. Thermal hydraulic results of the other assemblies here omitted are similar to those of assembly No. 9. 3.3. Analysis of loss of power accident without emergency cooling The calculation conditions are as followings: From 0 s, the CARR has been operated under full power normal condition for 50 s. At 50 s, the power of main pumps is switched-off, the safety shutdown system
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
275
responds 1s later. The residual heat is removed by the inertial running of the pump and NC process. During the whole process, the emergency cooling system is failed. The transient thermal hydraulic calculation results are shown in Fig. 1. Fig. 12(a) represents the transient behaviors of mass flow rate and core power. In the first 50 s, the mass flow rate in the primary loop and the core power keeps the initial value (660 kg/s in mass flow rate and 56.4 MW in core power) under normal operation. At time of 50 s, when the power of main pump is switchedoff, the mass flow rate decreases with the effect of inertial running of pumps. 1 s Mass flow rate (kg/s) Core power (10-1 MW)
600
Coolant enthalpy in hot channel ( KJ/kg)
Maximal coolant temperature Maximal clad temperature Maximal fuel temperature saturation temperature
200
40 500
Temperature (˚C)
30 400
20 300
10 200
0 150
160
170
180
190
160
120
80
200
100 40
0
(a)
(c)
240
0
100
200
300
0
400
100
200
300
400
Time / s
(b)
Time / s
0.16
900
Top HC enthalpy Bottom HC enthalpy
800
MC HC
0.08
700 600
0.00
quality
-1
Flow rate (kg/s) & Core power (10 MW)
700
500 400 300
-0.08
-0.16
200 -0.24
100 0
50
160
180
200
220
240
Time / s
0
(d)
50
100
150
200
250
300
Time / s
15
MDNBR
12
W-3 equation Sudo correlation Kutateladze correlation
9
6
3
1.0 0
(e)
0
100
200
300
400
Time / s
Fig. 12. Thermal-hydraulic calculation results of loss of power without emergency cooling. (a) Mass flow rate & core power profiles. (b) Maximum temperature profiles. (c) Enthalpy of the top and the bottom of HC profiles. (d) Equilibrium quality in MC & HC profiles. (e) MDNBR profiles.
276
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
later, the shutdown systems respond and the core power sharply decreases to a very low level of less than 10 MW. At about 170 s when the mass flow rate nearly decreases to zero, the NCCV is opened and the effect of NC driven force is prominent, thus the flow direction rapidly inverses and increases to a relative steady value of 20 kg/s. Fig. 12(b) represents the changes of maximum fuel, clad, and coolant temperatures as well as the coolant saturated temperature in HC. It can be found that the system temperature begins to sharply decrease because of the sharp decrease of the core power after 51 s. As the deceasing mass flow rate worsens the heat transfer capability, the coolant temperature quickly increases to the saturation level at time of 170 s when the mass flow rate decreases to near zero. After 170 s, the system temperature first sharply decreases with the increase of NC mass flow rate, then slowly decreases with the core power after the establishment of steady NC. In the whole process, the maximum fuel temperature is 240 °C which is still less than the melting temperature, 400 °C however, the coolant temperature is saturated. Fig. 12(c) represents the behaviors of coolant enthalpy of the inlet and outlet of hot channel. The inlet and outlet of hot channel exchange when the flow direction inverses. In order to avoid the nomenclative confusion, the words ÔtopÕ and the ÔbottomÕ are employed to represent the inlet and outlet of HC under normal operation conditions, respectively. As it can be seen that the coolant enthalpy of the bottom of the HC (outlet under normal condition) reaches to maximum when the mass flow rate equals zero. After the flow direction inverses, the enthalpy of the bottom sharply decreases with the NC mass flow rate increasing. The coolant enthalpy of the top of HC (inlet under normal condition) keeps steady before the flow direction inverses, at 170 s when NC begins, firstly the enthalpy of the top of HC sharply increases and then sharply decreases. Under this special condition, the temperature field is very complex and inversed heat transfer may occur. The transient behaviors of maximum equilibrium quality of coolant in the MC and HC are shown in Fig. 12(d). It is obvious that the boiling phenomena occur in both channels and the maximum quality of the HC reaches about 0.14 at time of about 170 s. Fig. 12(e) represents the transient behaviors of MDNBR in the HC. Results calculated by W-3 equation, Sudo correlation and Kutateladze correlation are illustrated in this figure. Large differences can be found among the three calculation results (even the tendencies are different). However, each correlation gives a safety prediction in its own scope of application. While the mass flow rate is nearly zero at the time of 170 s, the predictive MDNBR of Sudo correlation is less than 1.0, however, observing the predictions of both W-3 equation and Kutateladze correlation, there is still a large enough margin of Departure from Nucleate Boiling (DNB) just in this working condition for CARR. From Fig. 12(a)–(e), the calculation results of this case show that boiling occurs and lasts for about 6 s which are not allowed according to the Ôsafety design regulation of CARRÕ. However, in the actual engineering design, the emergency cooling system is powered by UPS which may last for more than 2 h. The following calcu-
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
277
lations will be carried out to present the important role that the ECS plays to ensure the safety in operation of CARR. 3.4. Analysis of loss of power accident with emergency cooling The calculation conditions are similar to the those in Section 3.3 except that the ECS is placed in service to pump the pool water to the core when the mass flow rate of the primary loop is less than 83.3 kg/s. The emergency cooling pump is switched-off Mass flow rate (kg/s) Core power (10-1MW)
210
600 50
500
Maximal coolant temperature Maximal clad temperature Maximal fuel temperature satuation temperature
180
40
Temperature (°C)
-1
Mass flow (kg/s) & Core power (10 MW)
700
30
400
20
300
10
200
0 7320
7340
7360
7380
7400
150
120
90
100 60 0 0
50
100
(a)
7400
0
7600
50
7200
7400
7600
Maximal quality in mean channel Maximal quality in hot channel
Top hot channel Enthalpy Bottom hot channel Enthalpy
-0.08
350
Equilibruim qualtiy
Ea nthalpy of coolant in hot channel ( kJ/kg)
150
Time / s
400
300
250
200
-0.12
-0.16
-0.20
150
-0.24 100 0
(c)
100
(b)
Time / s
200
400
7200
7300
7400
7500
0
7600
50
100
150
(d)
Time / s 180
7200
7400
7600
Time / s
W-3 equation Sudo correaltion Kutateladze correaltion
160
MDRBN
12
9
6
3
1.0 0 0
(e)
50
100
150
200
250
300
7400
7600
Time / s
Fig. 13. Thermal hydraulic calculation results of loss of power accident with emergency cooling. (a) Mass flow rate & core power profiles. (b) Maximum temperature profiles. (c) Enthalpy of the top and the bottom of HC profiles. (d) Equilibrium quality in MC & HC profiles. (e) MDNBR profiles.
278
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
after 2 h operation, and the removal of residual heat will be continued by NC which will be established after the automatic opening of NCCV. Fig. 13(a) shows the transient behaviors of mass flow rate and the core power. After 50 s normal operation, the power of main pumps is switched-off and the mass flow begins to decrease. 1 s later, the shutdown system responds and the core power sharply decreases to a very low level. At about 125 s when the primary mass flow rate decreases to less than 83.3 kg/s, the emergency cooling pump is being in service to pump the pool water into the primary loop. The mass flow rate keeps constant (83.3 kg/s) for 2 h before the emergency cooling pump is switched-off at 7325 s. Because of the neglected rotor inertia of emergency cooling pump as well as the prominent NC driven force, the mass flow rate sharply decreases to zero in less than 1 s and quickly increases to a steady NC mass flow rate of about 10 kg/s. The transient behaviors of the maximum fuel, clad, and coolant temperatures and the coolant saturation temperature in HC are shown in Fig. 13(b). Compared with Fig. 12(b), the peak system temperature is much lower. The peak fuel temperature is 130 °C (90 °C smaller than that in Fig. 12(b)) and the maximum coolant temperature (80 °C) has a margin of saturation with 40 °C. From Fig. 13(c), it is found the enthalpy of the bottom of the HC decreases slightly and the increase amplitude of the enthalpy of the top of HC is also much less than that in Fig. 12(c). It is obvious in Fig. 13(d) that the maximum quality in HC is less than zero, the coolant keeps to be subcooled and no boiling would occur. Fig. 13(e) represents the transient behaviors of MDNBR of HC. Similar to Fig. 12(e), at the time of 7325 s, while the mass flow rate is nearly zero, the predictive MDNBR of Sudo correlation is less than 1.0, however, predictions of both W-3 equation and Kutateladze correlation are much greater than 1.0 which exhibit a large margin of DNB. Comparisons between Figs. 12 and 13 show that the emergency cooling system plays an important role on ensuring the safety in operation of CARR. The activity of emergency cooling system effectively reduces the peak operational parameters such as the peak fuel temperature, peak coolant temperature and peak quality, etc.
4. Conclusions In this present study, theoretical calculation of thermal hydraulic characteristics of CARR is conducted. The theoretical analysis consists of initial steady and transient accidental analyses. The initial steady calculation result exhibited the important role that structure size played in flow distribution and the influence of uneven core power distribution could be neglected. The temperature field of fuel elements under unsymmetrical cooling condition was also obtained which was the base of further study such as transient-induced stress analysis, etc. Accidental analysis showed that the most dangerous condition occurred when the mass flow rate equaled zero. Result also showed that the service of emergency cooling system apparently reduced the peak thermal hydraulic parameters and thus more effectively ensured the safety of
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
279
CARR. Because of the adoption of modular programming techniques, this code is expected to be applied to accidental analysis of other types of reactors by easily modifying the corresponding function modules. Also, this code is expected to be validated against experimental data.
References Chen, J.C., 1966. Process design development. Int. Eng. Chem. 5, 322. Chisholm D. Naval Electronics Laboratory Report No. 310, 1967. Collier John, G., 1981. Convective Boiling and Condensation, second ed. McGraw-Hill. Gear, C.W., 1971. Numerical Initial Value Problems in Ordinary Differential Equation. Prentice-Hall, Englewood cliffs, NJ. Groeneveld, D.C., 2003. A look-up table for fully developed film-boiling heat transfer. Nucl. Eng. Des. 225 (1), 83–97. Guo, Y.J., 1995. The transient flow model for pumps of reactor system. Nucl. Sci. Eng. 15 (3), 56–61. Guo,Y.J., Zhang, J.L., Qiu, S.Z., Su, G.H., et al., 1997. MITARS-A thermal hydraulic analysis code for nuclear reactor system. In: Proceedings of the 8th International Topical Meeting on Nuclear Reactor Thermal Hydraulics. Sept. 30–Oct. 4, vol. 2, pp. 1212–1219, Kyoto, Japan. Hind marah, A.C. 1974. Gear: ordinary differential equation system solver, Lawrence Livarmore Laboratory, Report UCID-30001, Revision 3, December. Jens W.H., Lottes P.A. 1951. Argonne National Laboratory, Report No. 4627. Kanminaga, M., 1990. COOLOD-N: A computer code for the analysis of steady-state thermal-hydraulics in plate type research reactor. JAERI-M 90 (21). Rhosenow, W.M., Hartnett, J.R., 1973. Hand Book of Heat Transfer. McGraw-Hill, New York. Sudo, Y., Kaminaga, M., 1993. A new CHF correlation sheme proposed for vertical rectangular channels heated from both sides in nuclear research reactors. Trans. ASME. J. Heat Trans., 115–426. Su, G.H., Jia, D.N., Fukuda, K., Guo, Y.J., 2001. Theoretical study on density wave oscillation of twophase natural circulation under low quality conditions. J. Nucl. Sci. Technol. 38 (8), 607–613. Tao, W.Q., 1996. Heat Transfer. High Education Press Hall, Beijing. Thomas, L.C., 1992. Heat Transfer. Prentice-Hall, Englewood Cliffs, NJ.
Development of a thermal-hydraulic analysis code for CARR W.X. Tian a, S.Z. Qiu a
a,b,* ,
Y. Guo a, G.H. Su
a,b
, D.N. Jia
a,b
State key laboratory of Multi Phase Flow in Power Engineering, XiÕan Jiaotong University, XiÕan, Shaanxi 710049, PR China b Department of Nuclear and Thermal Power Engineering, XiÕan Jiaotong University, XiÕan, Shaanxi 710049, PR China Received 3 July 2004; accepted 13 September 2004
Abstract The China advanced research reactor (CARR) being built in Beijing, China, is a multipurpose research reactor for a variety of fields. Theoretical calculation of thermal hydraulic characteristics of CARR is presented in this paper. The theoretical analysis consists of initial steady and transient accidental analyses. Point reactor neutron kinetics model with six groups of delayed neutron is adopted for the solution of reactor power. All possible flow and heat transfer conditions are considered and the corresponding optional models are supplied in the theoretical calculations. A new simple and convenient model is proposed for the resolution of the transient behaviors of main pump instead of the complicated four-quadrant model. Gear method and Adams predictor–corrector method are adopted alternately for a better solution to such ill-conditioned differential equations corresponding to detail process. The initial multi-channel analysis shows that the effects of geometrical size on flow distribution play dominant role and the effects of core power distribution may be neglected. The temperature fields of fuel elements under asymmetrical cooling condition are also obtained, which are the bases for further study on transient-induced stress analysis, etc. Accidental analyses show that the activity of emergency cooling system apparently reduces the peak temperatures of fuel and coolant, peak quality and other operation parameters. Thus it effectively ensures the
*
Corresponding author. E-mail addresses: [email protected] (W.X. Tian), [email protected] (S.Z. Qiu).
0306-4549/$ - see front matter Ó 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.anucene.2004.09.003
262
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
safety in operation of CARR. Because of the adoption of modular programming techniques, this code is expected to be applied to accidental analysis of other types of reactors by easily modifying the corresponding function modules. Also, this code is expected to be validated against experimental data. Ó 2004 Elsevier Ltd. All rights reserved.
Nomenclature Cp De f g Gr h Hp I K Nu Pr q Re DT W
specific heat, J/kg k equivalent diameter, m friction factor, dimensionless gravity acceleration, m/s2 Grashof number, dimensionless heat transfer coefficient, W/m2 k pump head, m inertia of the rotor, kg m2 thermal conductivity, W/m k Nusselt number, dimensionless Prandtl number, dimensionless heat flux, w/m2 Reynolds number, dimensionless temperature difference, K mass flow rate, kg/s
Greek symbols / dimensionless length shown in Fig. 11 du,1 length shown in Fig. 11 du,2 length shown in Fig. 11 e roughness, m b expansion coefficient q density, kg/m3 x rotation speed, rps g running efficiency Superscript/subscripts b main fluid flow g steam/gas o normal working condition vf vapor/steam flow w wall
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
263
1. Introduction The China Advanced Research Reactor (CARR) is being built at the China Institute of Atomic Energy (CIAE, Beijing, China) to satisfy the multipurpose for a variety of fields in neutron scattering measurements, radioisotope production, neutron transmutation doping of silicon, etc. The CARR is a tank-in-pool, inversed neutron-trap-type reactor with a nuclear power of 60 MW. Slightly pressurized light water is used as the primary coolant and moderator, and heavy water as the reflector. The systematic diagram of CARR is shown as Fig. 1. The reactor body is immersed in a pool of 16 m in depth and the core is located 12 m below the pool water surface. The core is about 0.85 m in height and 0.451 m in diameter. The core is surrounded
Fig. 1. Schematic of CARR complex structure.
264
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
by a heavy water reflector tank, in which the under moderated neutrons leaking out through the core surface would be further moderated and the maximum thermal neutron flux is in a peak. Under CARRs normal working condition (full power working condition), the natural circulation crack valve (NCCV) is kept to be closed. The coolant is pumped to flow through the cold leg and the flow guiding tank, downward through the active zone of the core, then through the decay tank, the hot leg, the heat exchanger and recirculated to the main pump. The fluid flow path under normal working condition is denoted in solid arrow is shown in Fig. 2. However, if the CARR is shutdown, the NCCV will be opened. Due to the driven force caused by the density difference between the fluid in the water pool and the coolant in reactor vessel, the natural circulation (NC) will occur, by which the residual heat will be removed. However, the coolant flow direction will inverse during the change from forced circulation (FC) to NC, as shown in Fig. 2. Thus, the coolant flow path is: reactor pool-filter-decay tank-core vessel-flow guiding tank-NCCV-reactor pool. At a certain time in process of establishment of NC after shutdown, the mass flow rate in the core is becoming zero that is very dangerous for CARR. Therefore, safety analysis is very necessary to be carried out to ensure that CARR is in safety. The objective of the present work is to analyze the steady as well as transient thermal-hydraulic behaviors of CARR. For the initial steady analysis, the multi-channel
Hot water layer Reactor pool
Main Pump
Emergency Coolant system
Water
Flow guiding tank
Natural circulation
Nccv
Normal condition
Plate Heat Exchanger
Fluid flow Direction Under Normal Condition Fluid flow Direction Under Natural Circulation Condition
Active core
Filter screen
Decay tank
Fig. 2. Schematic diagram of the CARR.
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
265
analysis is conducted to obtain the detail flow distribution and temperature field of each fuel element. For the transient analysis, the behaviors of both the mean channel (MC) and the hot channel (HC) are obtained and discussed. All analyses are aimed at verifying if the CARR is safety in operation.
2. Theoretical model There are mass, energy and momentum conservation equations in the theoretical models. These equations are mutual coupled and are easy to find in reference (Guo et al., 1997; Su et al., 2001). The characteristics of the present theoretical model are introduced in the followings. 2.1. Physical models 2.1.1. Core power model Point reactor neutron kinetics model with six groups of delayed neutron is adopted for the solution of the power of the core. The reactivity feedback caused by the temperature changes of the moderator and the fuel is specially considered. 2.1.2. Heat transfer and flow friction models The main flow and heat transfer conditions might occur in the transient process are considered and the corresponding optional models are supplied as shown in Tables 1 and 2. The detail expressions of these correlations can be found in reference (Rhosenow and Hartnett, 1973). Linear interpolation has been done to ensure the continuum of heat transfer coefficients in the edge of different flow regimes and thus to prevent the numerical instability caused by the abrupt change of heat transfer coefficients. Chisholm correlation (Chisholm, 1967) is adopted For the calculation of two phase frictional multiplier. 2.1.3. Critical heat flux model Critical heat flux (CHF) is the prominent factor restricting the core power. There are more than 600 empirical or semi-empirical correlations for CHF based on Table 1 Flow regimes and heat transfer correlations Flow/heat transfer regime Single phase water Laminar flow Transition flow Turbulent flow Subcooled boiling Saturated boiling Transition/film boiling Single vapor flow
Heat transfer correlation Nu ¼ 0:17Re0:33 Pr0:43 Linear interpolation 0:11 lb lw 0.25 P/6.2
Pr Nu ¼ fRe 8X
Pr Prw
0:25
Gr0:1 (Collier, 1981)
(Rhosenow and Hartnett, 1973)
(Jens and Lottes, 1951) DT = 25q e Chen correlation, Expression omitted (Chen, 1966) AECL lookup Tablei Groeneveld h (2003) h ¼ 0:13k
q2 gbðT w T g Þ 1=3 l2
Cp l k
(Thomas, 1992)
266
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
Table 2 Flow regimes and friction correlations Flow regimes
Darcy friction coefficient correlations
Laminar flow Transition flow Turbulent flow Cole–Brook
f = C/Re different values for C for various geometrical channels f = 0.048 e 2:51 p1ffiffi ¼ 2 log e=D p ffiffi þ 3:71 f
Re
f
channels with various structures. Each correlation has its own very limited scope of application. It may be very difficult to find a correlation covering all flow conditions might occur during the establishment of NC of CARR. Thus, in the present study, three correlations are adopted and each one aims at different working condition. W-3 equation has been widely used in the design of pressurized power reactor and is adopted in the existing popular code such as RETRAN series, etc. Sudo correlation (Sudo and Kaminaga, 1993) is proposed on the basis of a great deal of CHF experimental data on rectangular channels. Compared with other CHF correlations, Sudo equation gives a relative conservative prediction. It was adopted in the design of the JRR-3M (Kanminaga, 1990) and its applicability under middle and high flow condition has been verified. Another CHF correlation, Kutateladze correlation (Thomas, 1992) is verified to be valid under the condition of pool boiling and low mass flow rate.
2.1.4. Equipment model Four-quadrant characteristic model has been widely adopted in universal commercial software such as RETRAN series for a better calculation of main pump transient performance. However, this model is inconvenient for use because of its strict and rigorous request for the detail pump performance parameters that are very difficult to obtain even for an actual certain pump (Guo, 1995). In this study, a new and convenient model on the computation of pump head is proposed by quoting the similarity law. The detail parameters are not necessary in the new model. Under all working conditions, the following relationship is always available among the pump output power Np, the pump lift Hp and the mass flow rate W N p ¼ W H p g:
ð1Þ
While the main pump being operated under switched-off condition, the effective output pump power equals to the decrease of the inertial kinetic energy of the rotor of the pump, that is N p ¼ Ixg
dx : dt
ð2Þ
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
267
For the same pump, according to the similarity law W x ¼ : ð3Þ W o xo Then, from Eqs. (1)–(3), it can be draw that 2 I g xo dW : ð4Þ Hp ¼ dt g Wo Eq. (4) is finally employed to calculate the transient behaviors of the main pump. It should be pointed out that there are four parallel identical main pumps arranged in the primary loop of CARR. The total mass flow rate equals the sum of the four branches and the equivalent pump head equals to that of any one of the four parallel pumps. The emergency coolant system (ECS) is connected with the CARR main structure at the cold leg as shown in Fig. 2. Under normal working condition, the ECS works as the cooling system of the pool water. When the mass flow rate in the primary loop is less than 83.3 kg/s, the NCCV will be automatically opened, the ECS is brought in service and the pool water is injected into the core. There are four plate heat exchangers with very complex structure placed in the primary loop of CARR. The fluid flow inside the exchanger is also very complex and disordered. Logarithmic temperature difference method (Tao, 1996) is adopted to calculate the performance of the heat exchangers. 2.1.5. Coolant thermophysical property model The correct coolant thermophysical property model is the key of success for a code. The thermophysical property correlations of water and vapor used in this study are chosen from the international standard (Rhosenow and Hartnett, 1973). Comparison between the model prediction and look-up table is conducted to verify and extend the applicability of those models under low pressure. 2.2. Nodalization The main principle of nodalization is that the different geometrical structure and physical boundary conditions of the core and pipes will be divided into different nodes. Two nodalization systems were adopted to meet the fact that the flow direction and flow path under normal working condition are different from those of NC condition. For normal working condition, the system node diagram is shown in Fig. 3. The active zone is modeled by MC and HC including the fuel, the clad and the coolant which consists 17 control volumes, respectively. For NC condition, during the phase of residual heat removal, the system node diagram is shown in Fig. 4. It is similar to that of normal working condition however, the flow direction, the flow path as well as the numeral order of volumes changed. The number of nodes of the core, the exchanger and pipes can be changed according to the calculation requirement.
268
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279 115-119
72 73 74
18 19 20
89 90 91
Plate heat exchanger
Hot channel
Core
92 93
Fuel
Clad
35 36 37
75 76
Coolant
Clad
52 53 54
3
58 59
4 5
Coolant
21 22
Fuel
Mean channel
2 38 39
Cold leg
1 Upper plenum
106 107 108
122-124
119-121
Main pump
55 Lower plenum 56 57 109-114,125-132
Flow direction
Hot leg
Fig. 3. The systematic diagram of nodalization of CARR under normal working condition.
Pipes 114-117
62
28 27 26
11 10 9
7
Coolant
core
65 64 63
Lower plenum
96 95
113 112
82 81 80
99 98 97
NCCV
Hot channel
Coolant
45 44 43
79 78
Fuel
25 24
59 58
61
Clad
42 41
Clad
Upper plenum
Fuel
Mean channel
60
118
Rector pool 119
8
6
Flow direction
Pipes: 1-5
Fig. 4. The systematic diagram of nodalization of CARR under natural circulation condition.
2.3. Calculation procedure The theoretical models are applied to every control volume. Thus, a set of equations with a stiff differential equation will be built. To solve it, a code has been made. The program flow diagram is shown in Fig. 5. The code is composed of several large modules including Data input module, Initialization module, Control system module, Numerical method module, auxiliary correlation module, Coolant thermophysical property module, Output modules, etc. Each module may be easily modified to
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
269
BEGIN
Data input
Initialization
Choice 1. Normal working condition that fluid flow downward through the core just as force convection
Normal/steady condition calculation
Time step setup etc.
working condition
Choice 1
Choice 2
Choice 2. The emergency cooling condition Choice 3
N Transient calculation
Simulation is over or not?
Y
Choice 3. NC working condition that fluid flow upward through the core
Output END
Fig. 5. The Program flow diagram.
be competent for the calculation under different working conditions of CARR or even other kinds of reactors. 2.4. Introduction of numerical method Once the nodalization of a nuclear reactor system was obtained, the transient analysis can be reduced to the solvement of a set of ordinary differential equations with the initial conditions as shown in Eqs. (5) and (6), where y is a vector including mass flow, enthalpy and other parameters of each control volume. Because of the speciality and complexity of nuclear reactor system, the corresponding differential equations are usually stiff differential equations (also called ill-conditioned equations)
270
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279 *
d y * * *0 ¼ f ðt; y ; y ; xÞ; dt *
y0: y ðt0 Þ ¼ ~
ð5Þ ð6Þ
Gear method (Gear, 1971; Hind marah, 1974) equipped with Adams predictor– corrector method was adopted for a better solution of stiff differential equations. Gear method is high efficient for the solution of ill-conditioned problems for its good stability, high precision, etc. It is the first order upwind difference method and the code of Gear method has the merit of self-stability and self-comparatively. Adams method required less computation time for its simpler iterative procedure, however, maybe failed for strong ill-conditioned equations. In this study, the two methods are alternatively used according to the stiffness of the set of equations in the whole solving process. For example, in the original stage after shutdown, the core power sharply decreased because of the large negative reactivity caused by shutdown rods, however, the fuel temperature, the coolant enthalpy, etc., might gently vary because of the relative delay of heat transfer. Gear method is thus used in this case. During the residual heat removal process, all parameters slowly change and Adams method is employed to save the running time. It is found in this study the alternate adoption of Gear method and Adams methods greatly reduce the computation time that is important for conducting a long-term simulation.
3. Calculation results and discussions 3.1. Principal design parameters of CARR The principal design parameters of CARR are shown in Table 3. The layout the reactor core of CARR is represented in Fig. 6. The core is composed of 17 standard fuel assemblies (numbering from 1 to 17), 4 following fuel assemblies (numbering from 18 to 21) and other experimental channels. The fuel element is shaped of plate type. All of those standard fuel assemblies have the same geometrical with symmetrical structure as shown in Fig. 7. Each standard assembly is composed of 20 fuel Table 3 Principal design parameters of CARR Core diameter/height (m) Elevation of reactor pool water surface (m) Core inlet/outlet temperature (°C) Core inlet/outlet pressure (°C) Core nuclear power (MW) Core thermal power (MW) Mass flow rate in primary loop (kg/s) Number of standard fuel element/assembly Number of following fuel element/assembly Type of fuel element
0.399/0.85 13.2 35/56.2 0.89/0.62 60 56.4 600 17 4 Plate
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
271
Fig. 6. Layout o the core of CARR.
Fig. 7. Detail structure of the standarad assembly.
plates that format 21 fuel channels. Several marginal channels have different gap sizes as shown in Fig. 7. Especially, channel No. 6 is only a half of the channel formatted by two extreme outside fuel plates of two adjacent assemblies. In the later steady analysis, the total 6 typical channels in each standard assembly will be analyzed separately. 3.2. Steady calculation result under normal working conditions A multi-channel steady analysis is conducted for CARR. Fig. 8 shows the radial distribution (percent) of the core power of each assembly. By simulating the whole reactor core, the detail flow distribution (percent) in the core was obtained as shown in Fig. 9. Very few difference was found among all mass flow percents of standard assemblies in Fig. 9(a), which indicated the influence of the non-uniform core power
272
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
4.86
5.2
4.75
3.1
5.4
3.11
4.96
5.35 5.37
5.9
5.37
5.37
4.95
3.1
5.41 3.11
4.93
4.8
5.15 4.87
4.92
Fig. 8. The percent of core power for each assembly (%).
0
Percent of mass flow
(a)
5
10
15
20
0.0498
0.0498
0.0495
0.0495
0.0492
0.0492
0.0400
0.0400
0.0395
0.0395
0.0390
0.0390
0.0385
0.0385 0
5
10
15
20
Number of assembly
Percent of mass flow distribution
(b)
0.0030
1
2
3
4
5
6
0.0030
0.0025
0.0025
0.0020
0.0020
0.0015
0.0015
0.0010
0.0010 1
2
3
4
5
6
Number of typical fuel channel of assembly No. 9
Fig. 9. Flow distribution. (a) For 21 assemblies of the core. (b) For 6 channels of No. 9 assembly.
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
273
radial distribution on flow distribution may be neglected. Fig. 9(b) represents the flow distribution (percent) of 6 typical channels for assembly No. 9. The strong influence of the channel gap size on the flow distribution is obvious in Fig. 9(b). Thermal-hydraulic analysis of 6 typical fuel elements (from channels No. 1 to No. 6) under unsymmetrical cooling condition (from channels No. 1 to No. 6) is carried out and the detail results of assembly No. 9 are shown in Figs. 10(a)–(d). Fig. 10(a) represents the coolant temperature along the axial flow direction. Coolant in channel No. 6 has the minimum temperature increase because of its maximum mass velocity. It also can be found that the coolant in channels No. 1 and No. 2 has the similar temperature increase because of the same gap size. Fig. 10(b) represents the pressure drop along the axial flow direction. Pressure dropped about 0.21 MPa through the core under normal working condition. Fig. 10(c) represented the maximal fuel temperature along the axial flow direction. Corresponding to Fig. 10(a), the minimum fuel temperature increase occurs in channel No. 6 and the fuel temperature profiles are similar to those in channels No. 1 and No. 2.
0
2
4
6
8
10
12
14
16
18
0.80
0.80
0.75
0.75
0.70
0.70
60
o
Temperature / C
55
50
Pressure / MPa
Ch 1 Ch 2 Ch 3 Ch 4 Ch 5 Ch 6
45
0.60
40
0
2
4
6
(a)
8
10
12
14
0.60
0.55
0
16
2
/ Location of maximal fuel temperature
o
120
Ch 1 Ch 2 Ch 3 Ch 4 Ch 5 Ch 6
80
60
40 2
4
6
8
10
12
Number of axial volume
14
Ch 1 Ch 2 Ch 3 Ch 4 Ch 5 Ch 6
0.504
(d)
8
10
12
14
16
18
0.5 represent the middle of the fuel plate
0.500
0.496
0.492
0.488
2
16
6
Number of axial volume
0.508
100
4
(b)
Number of axial volume
140
Temperature / C
0.65
0.55
35
(c)
Ch 1 Ch 2 Ch 3 Ch 4 Ch 5 Ch 6
0.65
4
6
8
10
12
14
16
Number of axial volume
Fig. 10. Thermal-hydraulic calculation results of 6 typical channels for No. 9 assembly (a) Axial coolant temperature profiles. (b) Axial pressure profiles (c). Axial maximum fuel temperature profiles. (d) Location of maximum fuel temperature for each control volume.
274
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
o Channel N Coolant
Channel N+1 Clad Fuel block
Clad
Coolant
Clad
Clad Tuo Tf,1
Tf , 2
Tci , 2
Tci ,1 Tco ,1
Tco , 2 u ,1
u ,2
o Fig. 11. Sketch map of fuel elements under asymmetrical cooling condition.
Fig. 10(d) represents the radial location of maximum fuel temperature for each control volume. Under symmetrical cooling condition, the fuel temperature reaches its maximum at the middle of fuel plate. In standard assemblies of CARR, for channels No. 1 to No. 6 are under unsymmetrical cooling condition, the maximum fuel temperature does not occur at the middle of each fuel plate as shown in Fig. 11. For the purpose of further study such as stress analysis of the fuel plate, it is very necessary to deduce the abnormal temperature distribution. In Fig. 10(d), the Y-coordinate variable is defined as d ¼ du;1 =ðdu;1 þ du;2 Þ;
ð7Þ
where du,1 and du,2 was shown in Fig. 11. From Fig. 10(d), it may be found that the location of maximum fuel temperature of channel No. 6 has a relative bigger departure from the middle of the fuel plate. However, for channels 1–3, the fuel temperature reaches maximum just at the middle. Observing the Fig. 10(a)–(d), the thermal hydraulic parameters in channel 1 are just similar to those of channel 2. It is verified reasonable and competent to choose 6 typical channels in one assembly to simulate to the whole core. Thermal hydraulic results of the other assemblies here omitted are similar to those of assembly No. 9. 3.3. Analysis of loss of power accident without emergency cooling The calculation conditions are as followings: From 0 s, the CARR has been operated under full power normal condition for 50 s. At 50 s, the power of main pumps is switched-off, the safety shutdown system
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
275
responds 1s later. The residual heat is removed by the inertial running of the pump and NC process. During the whole process, the emergency cooling system is failed. The transient thermal hydraulic calculation results are shown in Fig. 1. Fig. 12(a) represents the transient behaviors of mass flow rate and core power. In the first 50 s, the mass flow rate in the primary loop and the core power keeps the initial value (660 kg/s in mass flow rate and 56.4 MW in core power) under normal operation. At time of 50 s, when the power of main pump is switchedoff, the mass flow rate decreases with the effect of inertial running of pumps. 1 s Mass flow rate (kg/s) Core power (10-1 MW)
600
Coolant enthalpy in hot channel ( KJ/kg)
Maximal coolant temperature Maximal clad temperature Maximal fuel temperature saturation temperature
200
40 500
Temperature (˚C)
30 400
20 300
10 200
0 150
160
170
180
190
160
120
80
200
100 40
0
(a)
(c)
240
0
100
200
300
0
400
100
200
300
400
Time / s
(b)
Time / s
0.16
900
Top HC enthalpy Bottom HC enthalpy
800
MC HC
0.08
700 600
0.00
quality
-1
Flow rate (kg/s) & Core power (10 MW)
700
500 400 300
-0.08
-0.16
200 -0.24
100 0
50
160
180
200
220
240
Time / s
0
(d)
50
100
150
200
250
300
Time / s
15
MDNBR
12
W-3 equation Sudo correlation Kutateladze correlation
9
6
3
1.0 0
(e)
0
100
200
300
400
Time / s
Fig. 12. Thermal-hydraulic calculation results of loss of power without emergency cooling. (a) Mass flow rate & core power profiles. (b) Maximum temperature profiles. (c) Enthalpy of the top and the bottom of HC profiles. (d) Equilibrium quality in MC & HC profiles. (e) MDNBR profiles.
276
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
later, the shutdown systems respond and the core power sharply decreases to a very low level of less than 10 MW. At about 170 s when the mass flow rate nearly decreases to zero, the NCCV is opened and the effect of NC driven force is prominent, thus the flow direction rapidly inverses and increases to a relative steady value of 20 kg/s. Fig. 12(b) represents the changes of maximum fuel, clad, and coolant temperatures as well as the coolant saturated temperature in HC. It can be found that the system temperature begins to sharply decrease because of the sharp decrease of the core power after 51 s. As the deceasing mass flow rate worsens the heat transfer capability, the coolant temperature quickly increases to the saturation level at time of 170 s when the mass flow rate decreases to near zero. After 170 s, the system temperature first sharply decreases with the increase of NC mass flow rate, then slowly decreases with the core power after the establishment of steady NC. In the whole process, the maximum fuel temperature is 240 °C which is still less than the melting temperature, 400 °C however, the coolant temperature is saturated. Fig. 12(c) represents the behaviors of coolant enthalpy of the inlet and outlet of hot channel. The inlet and outlet of hot channel exchange when the flow direction inverses. In order to avoid the nomenclative confusion, the words ÔtopÕ and the ÔbottomÕ are employed to represent the inlet and outlet of HC under normal operation conditions, respectively. As it can be seen that the coolant enthalpy of the bottom of the HC (outlet under normal condition) reaches to maximum when the mass flow rate equals zero. After the flow direction inverses, the enthalpy of the bottom sharply decreases with the NC mass flow rate increasing. The coolant enthalpy of the top of HC (inlet under normal condition) keeps steady before the flow direction inverses, at 170 s when NC begins, firstly the enthalpy of the top of HC sharply increases and then sharply decreases. Under this special condition, the temperature field is very complex and inversed heat transfer may occur. The transient behaviors of maximum equilibrium quality of coolant in the MC and HC are shown in Fig. 12(d). It is obvious that the boiling phenomena occur in both channels and the maximum quality of the HC reaches about 0.14 at time of about 170 s. Fig. 12(e) represents the transient behaviors of MDNBR in the HC. Results calculated by W-3 equation, Sudo correlation and Kutateladze correlation are illustrated in this figure. Large differences can be found among the three calculation results (even the tendencies are different). However, each correlation gives a safety prediction in its own scope of application. While the mass flow rate is nearly zero at the time of 170 s, the predictive MDNBR of Sudo correlation is less than 1.0, however, observing the predictions of both W-3 equation and Kutateladze correlation, there is still a large enough margin of Departure from Nucleate Boiling (DNB) just in this working condition for CARR. From Fig. 12(a)–(e), the calculation results of this case show that boiling occurs and lasts for about 6 s which are not allowed according to the Ôsafety design regulation of CARRÕ. However, in the actual engineering design, the emergency cooling system is powered by UPS which may last for more than 2 h. The following calcu-
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
277
lations will be carried out to present the important role that the ECS plays to ensure the safety in operation of CARR. 3.4. Analysis of loss of power accident with emergency cooling The calculation conditions are similar to the those in Section 3.3 except that the ECS is placed in service to pump the pool water to the core when the mass flow rate of the primary loop is less than 83.3 kg/s. The emergency cooling pump is switched-off Mass flow rate (kg/s) Core power (10-1MW)
210
600 50
500
Maximal coolant temperature Maximal clad temperature Maximal fuel temperature satuation temperature
180
40
Temperature (°C)
-1
Mass flow (kg/s) & Core power (10 MW)
700
30
400
20
300
10
200
0 7320
7340
7360
7380
7400
150
120
90
100 60 0 0
50
100
(a)
7400
0
7600
50
7200
7400
7600
Maximal quality in mean channel Maximal quality in hot channel
Top hot channel Enthalpy Bottom hot channel Enthalpy
-0.08
350
Equilibruim qualtiy
Ea nthalpy of coolant in hot channel ( kJ/kg)
150
Time / s
400
300
250
200
-0.12
-0.16
-0.20
150
-0.24 100 0
(c)
100
(b)
Time / s
200
400
7200
7300
7400
7500
0
7600
50
100
150
(d)
Time / s 180
7200
7400
7600
Time / s
W-3 equation Sudo correaltion Kutateladze correaltion
160
MDRBN
12
9
6
3
1.0 0 0
(e)
50
100
150
200
250
300
7400
7600
Time / s
Fig. 13. Thermal hydraulic calculation results of loss of power accident with emergency cooling. (a) Mass flow rate & core power profiles. (b) Maximum temperature profiles. (c) Enthalpy of the top and the bottom of HC profiles. (d) Equilibrium quality in MC & HC profiles. (e) MDNBR profiles.
278
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
after 2 h operation, and the removal of residual heat will be continued by NC which will be established after the automatic opening of NCCV. Fig. 13(a) shows the transient behaviors of mass flow rate and the core power. After 50 s normal operation, the power of main pumps is switched-off and the mass flow begins to decrease. 1 s later, the shutdown system responds and the core power sharply decreases to a very low level. At about 125 s when the primary mass flow rate decreases to less than 83.3 kg/s, the emergency cooling pump is being in service to pump the pool water into the primary loop. The mass flow rate keeps constant (83.3 kg/s) for 2 h before the emergency cooling pump is switched-off at 7325 s. Because of the neglected rotor inertia of emergency cooling pump as well as the prominent NC driven force, the mass flow rate sharply decreases to zero in less than 1 s and quickly increases to a steady NC mass flow rate of about 10 kg/s. The transient behaviors of the maximum fuel, clad, and coolant temperatures and the coolant saturation temperature in HC are shown in Fig. 13(b). Compared with Fig. 12(b), the peak system temperature is much lower. The peak fuel temperature is 130 °C (90 °C smaller than that in Fig. 12(b)) and the maximum coolant temperature (80 °C) has a margin of saturation with 40 °C. From Fig. 13(c), it is found the enthalpy of the bottom of the HC decreases slightly and the increase amplitude of the enthalpy of the top of HC is also much less than that in Fig. 12(c). It is obvious in Fig. 13(d) that the maximum quality in HC is less than zero, the coolant keeps to be subcooled and no boiling would occur. Fig. 13(e) represents the transient behaviors of MDNBR of HC. Similar to Fig. 12(e), at the time of 7325 s, while the mass flow rate is nearly zero, the predictive MDNBR of Sudo correlation is less than 1.0, however, predictions of both W-3 equation and Kutateladze correlation are much greater than 1.0 which exhibit a large margin of DNB. Comparisons between Figs. 12 and 13 show that the emergency cooling system plays an important role on ensuring the safety in operation of CARR. The activity of emergency cooling system effectively reduces the peak operational parameters such as the peak fuel temperature, peak coolant temperature and peak quality, etc.
4. Conclusions In this present study, theoretical calculation of thermal hydraulic characteristics of CARR is conducted. The theoretical analysis consists of initial steady and transient accidental analyses. The initial steady calculation result exhibited the important role that structure size played in flow distribution and the influence of uneven core power distribution could be neglected. The temperature field of fuel elements under unsymmetrical cooling condition was also obtained which was the base of further study such as transient-induced stress analysis, etc. Accidental analysis showed that the most dangerous condition occurred when the mass flow rate equaled zero. Result also showed that the service of emergency cooling system apparently reduced the peak thermal hydraulic parameters and thus more effectively ensured the safety of
W.X. Tian et al. / Annals of Nuclear Energy 32 (2005) 261–279
279
CARR. Because of the adoption of modular programming techniques, this code is expected to be applied to accidental analysis of other types of reactors by easily modifying the corresponding function modules. Also, this code is expected to be validated against experimental data.
References Chen, J.C., 1966. Process design development. Int. Eng. Chem. 5, 322. Chisholm D. Naval Electronics Laboratory Report No. 310, 1967. Collier John, G., 1981. Convective Boiling and Condensation, second ed. McGraw-Hill. Gear, C.W., 1971. Numerical Initial Value Problems in Ordinary Differential Equation. Prentice-Hall, Englewood cliffs, NJ. Groeneveld, D.C., 2003. A look-up table for fully developed film-boiling heat transfer. Nucl. Eng. Des. 225 (1), 83–97. Guo, Y.J., 1995. The transient flow model for pumps of reactor system. Nucl. Sci. Eng. 15 (3), 56–61. Guo,Y.J., Zhang, J.L., Qiu, S.Z., Su, G.H., et al., 1997. MITARS-A thermal hydraulic analysis code for nuclear reactor system. In: Proceedings of the 8th International Topical Meeting on Nuclear Reactor Thermal Hydraulics. Sept. 30–Oct. 4, vol. 2, pp. 1212–1219, Kyoto, Japan. Hind marah, A.C. 1974. Gear: ordinary differential equation system solver, Lawrence Livarmore Laboratory, Report UCID-30001, Revision 3, December. Jens W.H., Lottes P.A. 1951. Argonne National Laboratory, Report No. 4627. Kanminaga, M., 1990. COOLOD-N: A computer code for the analysis of steady-state thermal-hydraulics in plate type research reactor. JAERI-M 90 (21). Rhosenow, W.M., Hartnett, J.R., 1973. Hand Book of Heat Transfer. McGraw-Hill, New York. Sudo, Y., Kaminaga, M., 1993. A new CHF correlation sheme proposed for vertical rectangular channels heated from both sides in nuclear research reactors. Trans. ASME. J. Heat Trans., 115–426. Su, G.H., Jia, D.N., Fukuda, K., Guo, Y.J., 2001. Theoretical study on density wave oscillation of twophase natural circulation under low quality conditions. J. Nucl. Sci. Technol. 38 (8), 607–613. Tao, W.Q., 1996. Heat Transfer. High Education Press Hall, Beijing. Thomas, L.C., 1992. Heat Transfer. Prentice-Hall, Englewood Cliffs, NJ.