Abstract
The LOCA (loss of coolant accident) is a kind of severe accident in the operation of PHWR (pressurized heavy water reactor) as well as other nuclear facilities, and possible cause of LOCA can be counted on the ballooning of pressure tube (PT) contacted to the outer calandria tube (CT) in the moderator system of CANDU-6 reactors. In the paper, we simulated the 150-kW experimental facility proposed by IAEA/ISCP, modeling the transient creeping behavior of pressurized tube heated with thermal radiation between the gaps of two concentric pipes. The outer boundary is simplified with a switched model that depends on the local temperature. With a multiphysical model supported by a commercial code, COMSOL multiphysics, the unsteady phenomena are simulated with models concerning various kinds of mechanics such as thermodynamics, nonlinear structural dynamics, and two-phase boiling heat transfer models.
1. Introduction
The CANDU (CANadian Deuterium Uranium) reactors have been adopted in Korea since late 1980s [1], operated with four CANDU-6 units generating 2.6 GW overall output in Wolseong nuclear plant. One of the most important safety issues in the CANDU reactors is the moderator subcooling margin [2] during a large loss of coolant accident to ensure the integrity of the fuel channels submerged in the heavy-water moderator.
International Atomic Energy Agency (IAEA) has launched a committee proposing an International Collaborative Standard Problem (ICSP) to provide a new contact boiling experimental data to assess the subcooling requirements for a heated pressure tube plastically deformed to contact with the calandria tube during a postulated large break LOCA condition. Since participating in the first Workshop of IAEA/ICSP on “HWR Moderator Subcooling Requirements to Demonstrate Backup Heat Sink Capabilities of Moderator during Accidents” held in Ottawa, Canada, 2012 [3], the authors have introduced a multiphysical code based on COMSOL [4] to show the solution of blind test, also compared with one-dimensional CATHENA code result [5]. After the workshop, the two-dimensional simulation is evolved to consider the boiling model at the outer calandria tube bounded emerged on the reservoir of coolant in the moderator system.
The plastic deformation of a pressure tube is brought forth with the high temperature from heat transfer of conduction and radiation where the model of Shewfelt and Godin [6] is applied to the structural dynamics equation of COMSOL. In this paper, the strain rate and the contact time of two tubes are computed for various heating-up conditions of the graphite heater located at the center of apparatus to investigate the sensitivity of open experimental conditions.
2. Numerical Method
2.1. IAEA/ICSP Benchmark
Figure 1(a) shows the IAEA/ICSP for the benchmark of the LOCA model that consists of a graphite heater and the surrounding pipes called pressure tube (PT) and calandria tube (CT) [7]. As the gap between PT and CT is filled with carbon dioxide for the insulation, the center axis of graphite heater is located eccentric 9.5 mm lower from the centroid of tubes. However, the convection effect is very small during the simulation, so this error is ignored in the numerical model of this study. The thicknesses of tubes are 4.4 mm (PT) and 1.42 mm (CT), and also the inner diameters mark 103.6 mm (PT) and 129 mm (CT), respectively. The electrical heater exerts the heat power measured for time (red line) in Figure 1(b) where the total time is 140 seconds, so this scenario is modeled similarly for COMSOL (ver. 5.3a) simulation (blue line).
(a)
(b)
2.2. Computational Model
The governing equations are summarized such as (1) heat transfer (conductive for each tube; radiative among heater and tubes) and (2) structural dynamics where the range of thermal stress lies beyond the linear elastic range [5]. Therefore, this problem becomes multiphysics involving the coupling of two kinds of phenomena.In the mechanical equilibrium at each time step, (1): the stress tensor () is conserved for the external volumetric force (), which should be zero for this problem except for boundaries; and the difference stress from the initial condition, denoted as subscript 0 in (2), is linearly proportional to the strain vector () with addition of thermal expansions depending on temperature () with thermal expansion coefficient (), which can be expressed as the function of displacement gradient () in (3). The elastic coefficient for the isotropic materials in (2) is composed of a matrix () dependent on Young’s modulus () and Poisson’s ratio (). This problem becomes multiphysical because (2) contains the temperature (T) that is the solution of (4) where , , and are material density, heat capacity, and thermal conductivity, respectively. The volumetric heat () distribution should be prescribed at the graphite heater located at the central axis in Figure 1(a).
Equation (1) is solved for components of the stress tensor, and the external force is given as boundary conditions. Note that the stress tensor, later reduced to the plane stress, should be three-dimensional though the computational domain is two-dimensional, which is due to the Poisson’s ratio in the linear elastic stiffness matrix of (2). For isotropic materials, the stiffness is expressed as a symmetric square matrix with the six degrees of freedom [7]:where is the Young’s modulus and is the Poisson’s ratio. Plane strain is hypothesized for the two-dimensional approximation as the tubes in Figure 1(a) are slender enough. The computational domain consists of Zr 2.5Nb (PT) mounted concentrically inside a section of Zircaloy 2 (CT): see Figure 2(a), which are given as a material properties in and , functions of local temperature.
(a)
(b)
The boundary conditions at the inner circles of tubes should be prescribed by the displacement () of the radial expansion from the initial radius aswhere only one-dimensional expansion is applied for the radial direction. Note that the thermal expansion term is modified from (2) where is not a constant but a function of . Therefore, four-point Gaussian quadrature is used for the integration in (6) [5].
In Figure 2(a), the boundary is pressurized inside the PT, and the pressure of carbon dioxide gas is almost atmospheric one, between PT and CT filled as a thermal insulator. Since this gas acts as an adiabatic thin layer between tubes, the natural convection is totally ignored in this study, only considering the dominant radiation effect. The volumetric heat () in (4) is modeled as Figure 1(b), and the thermal radiation at each surface is modeled with Stefan-Boltzmann law such aswhere is the outward-normal unit vector and is a constant of . Equation (7) is solved for the local temperature in the boundaries to calculate the surface radiation () in (8), and the incident radiation () is an explicit function of and the view factor determined iteratively by the mutual geometric configuration of surface-to-surface radiation boundaries where is emissivity marked in Figure 2(a), and their values of the core graphite heater, PT, and CT are 1.0, 0.8, and 0.34, respectively. Recall that Figure 2(a) is symmetric in the circumference direction, but is not simply one-dimensional because we should calculate and in (7)-(8) from the two-dimensional configuration of geometry. The outer temperature of water reservoir is fixed as 70°C, or subcooling of about 30°C under the standard atmospheric condition. However, after PT contacts with CT, the heat flux abruptly increases to full up the outer local temperature even to the boiling condition, which will be explained in Section 2.4.
Figure 2(b) is the triangular grids for the computation consisting of 3,424 elements with 952 edge and 20 vertex ones.
2.3. Nonlinear Thermal Deformation Model
As shown in (6), the thermal expansion is not a simple constant at the range of plastic deformation, and this nonlinear behavior is simplified with Shewfelt and Godin creeping model [13]:where the mechanical stress is (: pressure difference; : radius, : thickness of the tube) and and are times at temperature and , respectively. In (9), the time integration has no analytic solution, so the numerical method of four-point Gaussian quadrature is applied to compute it [5].
The strain () is related to the strain rate, in (9), and a function of temperature from (6):The thermal expansion, the unknown in (10), is plotted as a function of temperature in Figure 3, and it is remained constant at the linear elastic range, but increased exponentially at the nonlinear plastic range or . The expansion rate is abruptly decreased at , but increased again as heat is added: definitions of and are given in (9).
As the temperature rate () is computed in the solution of (4), (9) is transformed simply as a function of temperature [5]. Four-point Gaussian quadrature is used for the integration in (9) and (10). The integrated model in (10) is substitute to the right-hand side of (2) to get the stress field data in (1).
Figures 4(a) and 4(b) are the simulation result before and after contact, respectively, and the contact time is 74.6 s. After the expanded PT contacts with CT (see Figure 4(b)), the heat flux at the outer boundary in Figure 2(b) starts to be increased with natural convection in the water reservoir. However, the heat transfer soon transits to nucleate boiling after the outer temperature of CT exceeding the boiling point. It can be hardly observed that the boiling reaches to the critical heat flux. Therefore, the heat resistance between the contacted PT and CT due to the heat conduction at a very thin boundary of thickness , or the heat conduction gradient is initially assumed as , which is estimated from the experimental data of IAEA/ISCP. Because this value was not sufficient to bring about the critical heat flux, we numerically controlled this unknown value within a very restricted time of 0.1 seconds, for example, during the contact. The physical fact relying on this numerical correction is that the contact is initially an incomplete one, and highly amplified as the gap thickness () reaches to a very small value in a local “spots”, emitting very high heat flux to superheat the heavy water.
(a)
(b)
2.4. Boiling Heat Flux Model
When the temperature arise near 100°C, the boiling heat flux model is applied for the wall temperature higher than 100°C, which is the boiling point of water under the standard atmospheric condition [14]. As the CT outer wall temperature () increases in Figure 5, the boundary condition models are switched such as liquid convective, nucleate boiling, transition boiling, and film boiling. In Table 1, the models are listed with references.
2.4.1. Churchill and Chu Correlation
If the wall temperature is less than the boiling temperature, or 100°C, the wall heat coefficient is just that of a cylinder under the single-phase liquid free convection. There are many classical models, and Churchill and Chu Correlation [8] is applied for the mean convection heat transfer coefficient, : where is thermal conductivity; is the outer diameter of CT; and is the mean Nusselt number:In (12), the Nusselt number is correlated as a function of Rayleigh number, , and Prandtl number, :where , gravitational acceleration; , thermal expansion; 70°C, the temperature of reservoir; and , kinematic viscosity. The Prandtl number of heavy water is approximated as
2.4.2. Rohsenow and Zuber Correlation
In the range of , Rohsenow correlation [9] is used for the heat convection coefficient:where the coefficient is calculated from fluid and steam values: , and the subscript NB represents nuclear boiling.
The critical heat flux temperature () is calculated from the solution of a simple equation based on Zuber correlation [10] where the critical heat flux () is specified to from properties of the saturated steam. From the relation in (14) on that is a quadratic function on the differential temperature,Then, solving (15), 125°C is temperature at the critical heat flux.
2.4.3. Bjornard and Griffith Correlation
In the transition region between the critical heat flux and the film boiling, or , Bjornard and Griffith [11] suggest a simple model of interpolation with the temperature ratio, :where the critical heat flux () in (16) is calculated from Zuber correlation in (15).
All the models are plotted as the heat flux and the heat power per unit area versus wall temperature in Figure 6. The heat increases about 40 times at the critical heat flux more than that at the boiling temperature, decreasing in the transitional boiling, and increases slowly in the film boiling.
2.4.4. Bromley Correlation
For the range of , Bromley correlation [12] is used for the expression of heat convection coefficient:where the subscripts and stand for fluid (or saturated liquid water) and gas (or saturated pure steam) phase, respectively: therefore, under a given saturated temperature, . In the vapor mixture, denoted by the subscript , all the values are linearly interpolated with the temperature difference, from the empirical function in the superheated steam table.
The film boiling temperature is specified as 289.1°C, which is obtained from the graphical intersection point with curves generated with (16)-(17) such as Figure 6. This can be computed only after just several iterations.
3. Result and Discussion
3.1. Contact of PT and CT
In Figure 4, the contact occurs at 74.7 s, and the distribution of stress is reversed before (a) and after (b) the event in the radial direction. The maximum von Mises stress lies at the outer surface of PT before the contact: see Figure 4(a), but, however, it is located at the inner surface of PT after the contact: see Figure 4(b). The reason originates from the discontinuity of Shewfelt and Godin creeping model in Figure 3. The thermal expansion rapidly drops near °C, and this shock effects on the reversal of stress in the radial direction such as Figure 4(b).
The contact conductivity gradient elucidated in Section 2.3 is defined here to artificially induce a critical heat flux.where and are the contact time and the time duration of amplified heat conductivity gradient, in this research, respectively.
3.2. Temperature on PT and CT
In the counterpart experiment [3], the temperature is measured at five section rings in spanwise stations of the apparatus shown in Figure 1(a) where the selected data in Figures 7(a) and 7(b) are selected at the midspan. The computational model is a two-dimensional but eccentric one, so the difference in the circumference is not simulated. Instead, we change the amplification value () in (18). However, the slight change in does not show any remarkable change. When , the contact time and the peak temperature on the inner surface of PT are well predicted with small errors in Figure 7(a), but the transient temperature is affected with the high temperature at the outer boundary of CT. The combination of PT and CT is cooled in the computation faster than that of the experiment. As shown in Figure 7(b), the outer coolant water is quickly transit to the film boiling over the critical heat flux, coming back to nuclear boiling within ten seconds.
(a)
(b)
The mode at the outer boundary sensitively depends on the thermal resistivity at the contact surface between PT and CT. During the first stage of contact, for 0.1 seconds, it never reaches to the critical heat flux if we reduce the amplification factor below 15 (): see Figures 8(a) and 8(b). The effect of outer boundary temperature is reduced (Figure 8(a)), and the overshoot of outer temperature is far mitigated for the smaller amplification (Figure 8(b)). As the thermal expansion is not uniform in the three-dimension, the real experimental data lies between these two amplifications at each azimuthal location, and the thermal spots distribute from nucleate to film boiling.
(a)
(b)
4. Concluding Remark
The computational simulation has been performed for the IAEA-ISCP benchmark problem concerning a severe accident such as LOCA with the commercial code COMSOL multiphysics. The pressure tube (PT) expands with the radiative heat transfer from the graphite heater, soon contacting with the Calandria tube (CT), and the boundary condition at the outer CT is modeled with various physical correlations. The experimental measured data are compared with the result of numerical analysis.
The expansion of PT depends on the creeping of Zircaloy at high temperature, and the Shewfelt and Godin model has a genuine discontinuity at the maximum PT temperature. Using four-point Gaussian quadrature, the numerical integration for the thermal expansion shows a sharp discontinuous point to reverse the stress distribution in the radial direction neat the contact. The thermal resistance at the contact point is so important for the change of modes in the boiling process of the coolant water that the amount of change in the amplification of thermal conduction gradient at a thin contact layer, if it exceeds some threshold, brings out sensitively a considerable difference at the outer heat flux of CT.
The heat transfer observed in the parallel experiment only fits with a tuning of an artificial parameter of amplification factor for the heat conductivity gradient at the early stage of contact. The instant discharge of heat exceeding the critical heat flux can only activate the film boiling, or the reservoir water will remain in the nucleate boiling. However, the numerical data of transient surface temperature are compared with the experimental measurement, showing reasonable agreement overall in spite of much approximation of geometry and physics.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by (1) the Human Resources Development Program (Grant no. 20174010201350) of the Korea Institute of Energy Technology Evaluation and Planning (KETEP) grants funded by the Korea government (Ministry of Trade, Industry, and Energy) and (2) the Space Core Technology Development Program (NRF-2017M1A3A3A03015448) of Korea.