Abstract
The transient behaviors of natural circulation loop (NCL) are important for the system reliability under postulated accidents. The heat loss and structure thermal inertia may influence the transient behaviors of NCL greatly, so a transient analysis model with consideration of heat loss was developed based on the MATLAB/Simulink to predict the thermal-hydraulic characteristic of liquid metal NCL. The transient processes including the start-up, the loss of pump, and the shutdown of thermal-hydraulic ADS lead bismuth loop (TALL) experimental facility were simulated by using the model. A good agreement is obtained to validate the transient model. The appended structure would provide significant thermal inertia and flatten the temperature distribution in the transients. The oscillations of temperature and flow rate are also weakened. The temperature difference between hot leg and cold leg would increase with the decrease of heat loss, so the flow rate increases as well. However, a significant increase of hot section temperature may cause a failure of facility integrity due to the decrease of heat loss. Hence, the full power of the core tank may also be limited.
1. Introduction
Natural circulation is physical phenomena that occur in a gravity environment when a geometrically distinct heat sink and heat source are connected through a flow path. In NCL the fluid motion is generated by density differences in the fluid due to the temperature gradients or the phase change, which also means that no external sources of mechanical energy for the fluid motion are involved. Therefore, the natural circulation principles are very attractive to the designers of nuclear power plant systems and components; not only could the system cost be reduced but also the reliability may be enhanced.
Liquid metals such as lead (Pb), lead bismuth eutectics (LBE), sodium (Na), and sodium-potassium alloy (NaK) are candidate coolants for the next generation liquid metal reactors [1–5]. In general, a high mass flow rate could be provided in liquid metals natural circulation loops due to their significant thermal expansion coefficient which make it possible to use natural circulation based primary coolant systems and natural circulation based residual heat removal systems in the design of nuclear power plants. We also proposed a conceptual design of a small modular natural circulation liquid metal fast reactor with alkali-metal thermal-to-electric conversion (AMTEC) units and provided a straightforward comparison of steady state natural circulation characteristics for different liquid metals in previous studies [6]. On the other hand, the transient behaviors of NCL would be more important because it can influence the system performance under postulated accidents directly. Unlike the forced circulation, the operation of NCL may not be easily controlled. The transient behaviors of NCL strongly depend on the interaction between buoyancy and frictional forces, so the performance of NCL may be influenced more significantly by other factors such as flow resistance, thermal inertia, and environmental heat loss.
Most of the experimental investigations and theoretical analysis methods have been carried out for symmetrical and asymmetric (determined by the relative position of heater and cooler) water NCL. Vijayan et al. investigated the steady state and stability characteristics of single-phase natural circulation in a rectangular loop with different heater and cooler orientations experimentally and theoretically [7]. Basu et al. theoretically investigated the dynamic response of a single-phase rectangular natural circulation loop to different excitations of input power [8]. Misale and Frogheri studied the influence of pressure drop on the transient behavior and stability of water experimental NCL by changing the diameter of orifices in the loop [9]. Besides, some computer codes, such as RELAP5 system code, could be applied to simulate the transient behavior of the water NCL, and the calculations could provide results which agree well with the experimental data [10–12].
Some investigations are also carried out for the liquid metal NCL. Ma et al. experimentally investigated the capability and stability of TALL experimental facility and used RELAP5 to estimate the natural circulation flow rate [13]. Yiqiong used TRAC/AAA code to build a more detailed model to simulate the transient behavior of the loop. A good agreement between simulation results and experimental data had been obtained [14]. Abánades and Peña used computational fluid dynamic codes to analyze the free convection cooling mode of the reference liquid metal-cooled accelerator driven system [15]. Wu et al. employed the linear analysis method to study the stability of Argonne LBE NCL and give the steady state fluid velocity and temperature difference between hot leg and cold leg under different heating power [16, 17]. Although there are many researches of natural circulation of liquid metals, the study of transient behavior for liquid metal NCL is still limited. Furthermore, the influence of structure thermal inertia and heat loss to ambient on the performance of liquid metal NCL has not been paid much attention to in some of the previous studies. Therefore, more research is needed to analyze factors that influence the transient behaviors and steady state conditions of the liquid metal NCL.
In order to analyze the loop transient behavior and evaluate the natural circulation performance during various transient processes, a liquid metal NCL model was developed. The loop model including core tank model, pipe model, and intermediate heat exchangers (IHXs) model is built based on the MATLAB/Simulink platform. In order to evaluate the effectiveness of the loop model, the transient processes including the start-up, the loss of pump, and the shutdown of TALL experimental facility were simulated. Moreover, the influence of appended structure thermal inertia and heat loss to ambient on the performance of NCL is evaluated in present work.
2. Calculation Model and Method
A general liquid metal NCL (TALL test facility) is shown in Figure 1. The primary loop of the natural circulation system mainly consists of core tank (heat source), IHX (heat sink), pipes (fluid flow path), electromagnetic (EM) pump (optional), and so forth. During operation, the coolant initially at rest may be accelerated by the EM pump. At the same time, it would be heated by heaters to a specific temperature when it flows through the core tank. Then, the coolant flows from the outlet of core tank to the expansion tank through a long vertical pipe. The temperature of coolant would be decreased to a specific value after it crosses the IHX. Finally, the coolant flows down through a vertical pipe and returns to the inlet of core tank to complete the circuit.
It can be seen that once the temperature difference is established, the coolant could be driven by the buoyancy force without the operation of EM pump. So the natural circulation transient behaviors and steady state conditions are mainly determined by the interaction between buoyancy and frictional forces. The temperature and flow rate oscillations would always exist before the loop reaches steady state. Moreover, a fluctuation of temperature may cause a flow rate oscillation again even when the loop is initially at a steady state and vice versa. Therefore, the transient process and steady state could be predicted very well only when the temperature distribution and friction losses are exactly evaluated.
2.1. Governing Equations and Correlations for Different Models
As can be seen in Figure 1, the core tank holds the immersion heater which simulates the fuel rod and provides specific heating power to heat the coolant directly. The heater consists of NiCr, filler of MgO, and stainless steel 316 cladding. The transient heat conduction equation of structure neglecting axial heat conduction is expressed as follows:
The coolant in the annular channel would be heated by the heater. It should be noted that in the TALL facility some guide heaters and insulation are placed on different sections of the primary loop in order to reduce the environmental heat loss. It means that the pipe wall, structure components, and insulation are also heated by the coolant. The energy equation with consideration of the boundary conditions for the incompressible coolant in the core tank would be where the heat transfer coefficient could be obtained by estimating the Nusselt number [18]
As for the pipe wall, structure components, and insulation, the energy equation is the same as (1) except that for wall and insulation. The outer surface of the pipe insulation is assumed to be cooled by the air in room temperature (~300 K). The Nusselt number for natural convection heat transfer from a horizontal and vertical cylinder can be expressed in the following general form: where the coefficients and depending on Grashof number could be found in [19].
In the pipe model, the energy equation for the coolant is
The heat transfer through the pipe, structure components, and insulation could also be described by using (1).
In the TALL facility, the IHX is a single tube LBE-Glycerol heat exchanger with the primary coolant in the inner tube and the secondary coolant in the annulus. In the IHX model the energy equation for primary side and secondary side would be where the heat transfer coefficient in the secondary side of IHX could be estimated by using the following correlations [20]: where . As for the tube wall the governing equation would be the same as (1) with .
In the simulation the coolant is assumed to be incompressible so the mass flow rate is only a function of time and independent of space coordinate; namely,
On the other hand, the momentum equation for one-dimensional incompressible flow is When Boussinesq approximation is employed, the density is taken to be constant except in the body force term. Then integrating on both sides of (9) along the loop we could obtain where . For laminar flow and for turbulent flow . is the difference in pressure between the outlet and inlet of EM pump. It would be equal to zero when the pump is out of work and the coolant is completely driven by the buoyancy.
In addition, the liquid metal NCL model also includes some other submodels. The junction model is mainly used to estimate the coefficient of local resistance. For sudden expansion and contraction of a stream with uniform velocity distribution, the coefficient would be
The resistance coefficient of the bend is
As for the multihole plates the resistance coefficients could be obtained by the following correlations: where the cross-section coefficient is equal to the ratio of flow area of the obstruction cross-section to the area of the obstruction front. The coefficients and are a function of , respectively, and could be determined on the graph of [21].
2.2. Discretization of the Energy Equation
The liquid metal NCL model is mainly composed of core tank model, IHX model, and pipe model. Figure 2 shows the calculation cells for these models.
(a)
(b)
(c)
Figure 2(a) shows the radial nodding in an axial slice of the core tank. The NiCr part of the heater is divided into three (user input) nodes; the filler, clad, coolant, pipe wall, and appended structure have one node each. The insulation is divided into three (user input) nodes. Because the SIMULINK is a time-based simulation platform, the discretization of the partial derivative with respect to time that appears in the governing equations is not necessary [22]. By using the finite volume method the general form of the discretization of the energy equation for structure could be expressed as where for heat conduction and for convective heat transfer. It should be noted that the thermal conductivity of the adjacent structures may be very different. So the numerical treatment of coefficient at interface could be expressed as follows [23].
According to Fourier’s law the heat flux density across the interface (Figure 3) would be So the right side coefficient would be . It can also be seen that the temperature at the interface would be Obviously, the interface temperature is determined by the conductivity of materials (also a function of temperature) and the node positions. As for the coolant in the channel (Figure 2(a)) the discretization of (2) would be where the average coolant temperature in the control volume . In addition, the boundary conditions need to be specified to solve the discretized equations, so the convective heat transfer coefficient of the exterior of insulation and the air temperature are specified. As can be seen in Figure 4 that five nodes are used to simulate the heaters of core tank, three nodes are used to simulate the upper part of core tank and junction 001 would be used to estimate the pressure drop coefficients when the coolant flows through the multihole plates.
In the radial direction of IHX model (Figure 2(b)) two nodes are used to represent the primary coolant and secondary coolant, respectively, and one node is used for the tube wall. The same discretization method is employed to obtain discretized equations in (6) which are similar to (17) except that there is only one convection heat transfer surface. It can be noted that the whole secondary loop is not simulated in this loop model because the main function of the secondary loop is to provide the steady coolant for the LBE-Glycerol heat exchanger. On the other hand, the heat loss through the insulation of the secondary side of IHX is neglected due to a relatively lower secondary coolant temperature and smaller overall heat transfer coefficient compared to primary loop.
In Figure 2(c) the pipe model is represented by six mesh intervals: coolant, pipe wall, and appended heat structure (represent the guider heaters, valve, thermocouples, etc.) have one node each and another three are put for modeling insulation. The discretized equations for pipe model could also be obtained by using the same method that is employed in core tank model. In Figure 4 the pipes are connected by many junctions which are used to estimate the coefficient of local resistance.
2.3. Solution Procedures
The transient temperature distribution at the next time step can be numerically obtained by solving the discretized energy equations for given initial conditions and boundary conditions. Then, the buoyancy term in (10) can be solved by integrating the temperature distribution along the entire length of the loop. At the same time, the mass flow rate at the next time step through the loop can be obtained by solving (10). Finally, the mass flow rate is substituted into the energy equations for the coolant in the different sections of the loop and the next iteration is started.
It can be seen from (14) and (17) that the discretized energy equations for structure and coolant could be rewritten in the form of matrix multiplication expressed as follows: where ; is a matrix and the general source term . Figure 5(a) shows the simplified solution diagram for the heated part in core tank model by using MATLAB/Simulink.
(a)
(b)
In core tank model the coefficient matrix and general source term are solved, respectively. The heat source and the boundary conditions including the inlet temperature of core tank and the environment temperature are used to calculate the general source term. The temperature distribution would be used to determine the material properties and the buoyancy term in the momentum equation. Meanwhile, the outlet temperature of the model would be considered as the inlet temperature of the next adjacent model. As a result, the NCL is composed of these connected models which represent the transient temperature distribution in different sections of the loop.
Figure 5(b) shows the Simulink model for momentum equation solver. The driving force is mainly provided by the pump and the buoyancy caused by the temperature difference between the hot leg and cold leg. The buoyancy would be determined by the temperature distribution in the core tank, IHX, and pipes. At the same time, the flow losses take place either as friction loss due to wall friction or as local flow loss depending on channel geometry. The friction loss in each section of the loop has been calculated in the corresponding model. The local flow losses are calculated in the junction models. Then the transient mass flow rate would be determined by the interaction between driving forces and resistance. Finally, the mass flow rate at the next time step is transferred to other models to obtain the temperature distribution in each section of the loop. Consequently, the transient behaviors of the liquid metal NCL can be predicted by using the Simulink model for the loop. On the other hand, because the problem is supposed to be stiff, an implicit continuous variable-step solver (ode15s) is chosen for the simulation. Some other parameters such as loop geometry and material properties (the appended structure is assumed to be composed of stainless steel 316) could be found in [14, 24–26].
3. Results and Discussion
Some transient experiments had been carried out on the TALL facility to study the steady state and transient thermal-hydraulics performance of LBE-cooled reactors [14]. Meanwhile, the TRACE code was used to simulate the transient behavior of TALL facility without consideration of thermal inertia caused by appended structure such as rope heaters, valves, and sensors. In the present study the start-up, loss of primary pump, and shutdown transient experiments are simulated with appended structure to test the liquid metal NCL model and study the influence of thermal inertia on the performance of NCL.
3.1. Start-Up Experiment
In the start-up experiment the loop would be started by forced circulation. The initial LBE temperature in the primary loop is about 200°C. The initial mass flow rate for primary and secondary loop coolant is 0 kg/s. Meanwhile, the heater in the core tank was out of work at the initial state. At about 140 s the transient experiment was started. The primary and secondary loop pumps were switched on to work, whereas the power was supplied to the core tank heater. The mass flow rate of LBE and Glycerol would immediately reach their final value of 0.91 kg/s and 0.69 kg/s, respectively. A constant power of 8.5 kW is supplied to the heater at the same time. Figure 6 illustrates the experiment and simulation results of the start-up forced circulation.
(a)
(b)
It can be seen from Figure 6 that at the beginning of the experiment the outlet temperature of IHX would decrease due to the cooling by the secondary coolant flow. After a few seconds, the cold coolant reaches the inlet of core tank. It should be noted that the inlet temperature of core tank does not decrease sharply as the outlet temperature of IHX. It seems that the structure of the cold leg plays an important part in flattening the temperature distribution due to the thermal inertia. On the other hand, the minimum inlet temperature of core tank is slightly higher than that of IHX. This is because of the heat transfer from the cold leg structure to the coolant due to a higher temperature of appended structure. The outlet temperature of core tank would increase due to the heat supply. However, because of the decrease of the inlet temperature, the core tank outlet temperature would stop increasing in a short period of time and then keep increasing. The inlet temperature of IHX also increases with the outlet temperature of core tank. Finally, a steady state loop temperature distribution would be established after about 4000 s. By comparing Figures 6(a) and 6(b) we can see that the outlet temperature of core tank/IHX is higher than the inlet temperature of IHX/core tank at steady state. The heat loss to ambient would decrease the temperature of coolant when it flows through the pipes.
It can be seen from Figure 6 that both of the simulation results show good agreement with the experimental data. Moreover, both results indicate that the outlet temperature of core tank/IHX is higher than the inlet temperature of IHX/core tank at steady state because the heat loss to ambient is considered in present simulations and [14]. However, the steady state temperatures are still somewhat different which may be caused by the different calculation methods and assumptions in the simulations. Furthermore, the influence of appended structure thermal inertia was also considered in present study which, however, was not included in [14]. Therefore, the operating temperature of the loop in the transient would rise more slowly with consideration of appended structure which is closer to the experimental process. It is clear that in the start-up transient not only the coolant in the loop but also all the structure surrounding the loop would be heated. Therefore, the thermal inertia of the structure will extend the time required to reach steady state. It should be noted that a completely stable initial condition is not easy to be established because of a small natural circulation residual flow and the use of rope heaters in the experiment. On the other hand, it is difficult to predict the structure surrounding the loop and the environmental heat loss very exactly; thus, a longer time is taken to reach the steady state in the experiment than in both of the simulation results. In addition, it seems that there is more heat loss in the experiment than in the prediction. As a result, a good simulation of thermal inertia and heat loss is needed to obtain an accurate prediction of the transient behavior of the NCL.
3.2. Loss of Primary Pump Experiment
When the primary pump is out of work, the coolant would be completely driven by the buoyancy. It means that the power of the heater would be removed by natural circulation and the loss of flow accidents can be well prevented. Therefore, a natural circulation based primary circuit is often preferred in some designs of advanced nuclear power plants.
In the loss of primary pump experiment the loop would initially work at steady state forced circulation with a heater power of 12.4 kW. At about 550 s the primary pump was switched off, and then the rotational speed of primary pump would reduce to 0 in about 2 s. Figure 7 illustrates the variation of temperature and velocity.
(a)
(b)
(c)
It can be seen that the velocity of coolant would decrease dramatically when the pump stops working. Therefore, the outlet temperature of core tank increases when the heater power is kept constant, while the outlet temperature of IHX decreases suddenly because the primary coolant is further cooled by the secondary coolant when the velocity of primary coolant decreases. It should be noted that Glycerol was chosen as the intermediate fluid in the secondary loop. Its boiling point (~290°C) is higher than the melting point of LBE (~125°C). It means that the minimum temperature in the loop could be much greater than the melting point of LBE so that the solidification of LBE in the heat exchanger can be avoided. It is clear that the structure again causes a thermal inertia in the transient experiment. The heat supplying and absorbing of the structure would influence the transient temperature distribution in cold leg and hot leg. Meanwhile, the temperature difference between the inlet and outlet of core tank/IHX at steady state would increase greatly when the natural circulation flow rate is only about half of the flow rate provided by forced circulation. It can also be noted that the inlet temperature of IHX is obviously lower than the outlet temperature of core tank under natural circulation condition.
A good agreement between the simulation results and the experimental data can be obtained. In addition, the simulation results show that the oscillations of temperature and flow rate would be more obvious (amplitude and frequency) when lower thermal inertia is employed. As we know that the structure including the pipe wall was not generally considered in most of previous studies on the stability of NCL [7, 8, 10, 16, 17]. It seems that the thermal inertia would play an important part in the stability of NCL. This may also explain why the instability of NCL is sometimes overpredicted in the simulation compared to the experiment. From Figure 7 we can also see that the simulation results reach the stable state earlier than the experimental data. This may be because the thermal and coolant inertia caused by the auxiliary equipment (e.g., expansion tank) are not accurately calculated in the simulation. On the other hand, the appended structure is assumed to be uniformly distributed on the outer surface of pipe in the simulation which is different from the actual conditions.
3.3. Shutdown Experiment
The shutdown procedure of the TALL facility is also simulated in present work. The initial power is 17 kW when the mass flow rate is about 0.91 kg/s. At the beginning of the transient experiment, the power would be decreased to zero at about 280 seconds while the forced natural circulation working condition was kept. Figure 8 presents the inlet and outlet temperature variation of core tank and IHX.
(a)
(b)
It can be seen from Figure 8 that the outlet temperature of core tank would decrease dramatically due to the shutdown of heating power while the rate of decrease of IHX inlet temperature is a little slower due to the thermal inertia. On the other hand, it can be observed that the cold coolant at the outlet of IHX would be somewhat reheated by the pipe and structure when it flows through the cold leg. Hence, its temperature rises when the coolant reaches the inlet of the core tank. An opposite feature is observed in the start-up experiment.
Generally, both of the simulation results show good agreement with the experimental data in Figure 8. However, the operating temperature of the loop would be decreased more quickly with lower thermal inertia. The simulation results would be closer to the experimental data when the appended structure is considered.
4. Influence of Heat Loss on the Performance of NCL
In most of previous studies the heat loss to ambient in the NCL was not considered. Generally, the hot leg and cold leg are assumed to be perfectly adiabatic in the theoretical analysis. It had been proved to be reasonable as the calculation results show good agreement with the experimental data. It seems that the heat loss to ambient is limited when the loop operating temperature is not much higher than ambient temperature or the experimental loop size is not very large. Hence, the influence of heat loss on the performance of NCL would depend on the specific case.
In the loss of primary pump transient of TALL facility it can be seen that the inlet temperature is significantly lower than the outlet temperature of core tank especially when a stable natural circulation is established in primary loop. In order to evaluate the effect of heat loss on the performance of NCL, the thermal conductivity of the insulation is set to be a minimum value (e.g., 10−4 W/m/K) to simulate the loop performance without heat loss. Table 1 shows the thermal conductivity of insulation material (asbestos) at different temperatures. Figure 9 illustrates the calculation results with and without consideration of heat loss.
(a)
(b)
(c)
(d)
It is clear that the heat loss plays an important part in the performance of the loop, no matter what the loop is at the condition of forced circulation or natural circulation. Obviously, the average operating temperature of the NCL increases a lot (~30°C) when the hot leg and cold leg are assumed to be adiabatic. Under this condition the total heat is completely removed by the IHX, so the primary coolant temperature would increase in order to extract more heat to the secondary loop. On the other hand, the average hot leg temperature would increase ~40°C while the cold leg temperature increases ~25°C under natural circulation. It can be concluded that more heat is dissipated from the hot section due to its higher temperature. In addition, from Table 1 we can see that the thermal conductivity of insulation material will increase with temperature which may also cause an increase of heat loss through the hot leg.
The capability of the natural circulation of the loop will also be influenced by the heat loss to ambient. It can be seen from Figure 9(c) that the temperature difference between hot leg and cold leg would be increased by decreasing the heat loss through pipes. It means that the flow rate or the capability of the natural circulation could also be increased due to the increase of buoyancy (Figure 9(d)). It is preferred to have as high natural circulation flow as possible in a reactor design. However, reducing the heat loss may also bring some undesirable results. The corrosion rate may also be increased with the operating temperature which undermines the integrity of the facilities. For example, the outlet temperature of the core tank is restricted to ~460°C in the TALL facility. It also means that the total core tank power will therefore be limited.
5. Conclusions
A liquid metal NCL computing approach with the heat-loss model, including core tank model, IHX model, and pipe model, is developed based on the MATLAB/Simulink platform to investigate the transient behaviors and steady state characteristics of the loop. Three types of transient experiments (the start-up, the loss of pump, and the shutdown) on the TALL facility are simulated to validate the model and evaluate the influence of appended structure thermal inertia and heat loss on the performance of NCL.
A good agreement between the simulation results and the experimental data of the transient behaviors for the primary LBE coolant is obtained, which provides good validation for the MATLAB/Simulink models. It is found that the thermal inertia greatly influence the temperature variation during the transient processes. The coolant temperature distribution in the pipe would be flattened due to heat exchange with structure. It seems that the simulation results reach the steady state earlier than the experimental data. This may be because the thermal and coolant inertia caused by the auxiliary equipment are not accurately considered in the simulation. On the other hand, the structure thermal inertia often neglected in the previous analysis of the stability of NCL may play an important role in the stability behavior of NCL. The oscillations of temperature and flow rate in the experiment may be not obvious which is unlike the simulation results. This may explain why the NCL is supposed to be unstable but the instability behavior is not observed in the experiment.
The influence of heat loss through the pipes on the performance of NCL is also evaluated. Obviously, the average operating temperature of the NCL would increase with the decrease of heat loss. Furthermore, the temperature of hot section would increase more significantly compared to the cold section temperature. It may even exceed the upper temperature limit and so that undermine the integrity of the facilities, so the heat power is also limited at the same time. However, the capability of the NCL could be increased due to the decrease of heat loss.
Nomenclature
| Area (m2) | |
| : | Specific heat capacity (J·kg−1·K−1) |
| : | Loop diameter (m) |
| : | Friction coefficient |
| : | Gravity (m·s−2) |
| : | Grashof number |
| : | Heat transfer coefficient (W·m−2·K−1) |
| : | Thermal conductivity (W·m−1·K−1) |
| : | Local resistance coefficient |
| : | Length (m) |
| : | Nusselt number |
| : | Pressure (Pa) |
| : | Prandtl number |
| : | Heat flux density (W·m−2) |
| : | Volumetric rate of heat generation (W·m−3) |
| : | Radius (m) |
| : | Mean radius of curvature of the bend (m) |
| : | Reynolds number |
| : | Space coordinate (m) |
| : | Time (s) |
| : | Temperature (°C) |
| : | Velocity (m·s−1) |
| : | Volume (m3) |
| : | Mass flow rate (kg·s−1). |
| : | Structure |
| : | Fluid |
| : | Heat transfer |
| : | Primary loop |
| : | Secondary loop |
| 0: | Reference value. |
| : | Thermal expansion coefficient (K−1) |
| : | Angle of bend of the curved channel (°) |
| : | Inclination angle (°) |
| : | Fluid density (kg·m−3). |
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.