Abstract
In order to study the evolution process and hydraulic characteristics of pressurizer insurge in CPR1000 pressurized water reactor (PWR), the computational fluid dynamics (CFD) of three-dimensional unsteady heat transfer was used to capture the temperature and velocity fluctuation intensity of the mixing of the hot and cold. The results show that Realizable k–ε Turbulence Model combined with VOF multiphase heat transfer model can effectively predict the development trend of pressurizer insurge process. The small diameter pressurizer surge line of CPR1000 enhances the intensity of velocity fluctuation. From the essence of flow and heat transfer, it is concluded that buoyancy force can increase the degree of fluctuation and make an accelerated effect on the influx cold fluid. The electric heater inside the pressurizer should be arranged as far as possible in z<-0.45m and z>0.45m; it is beneficial to improve its harsh operating environment. This research can provide reference for the structural design of pressurizer and the layout optimization of the electric heater.
1. Introduction
There are several types of nuclear power plants under construction or operation in China, such as CPR1000, EPR1000, AP1000, Hualong-1, and so on, most of which are based on the technology of pressurized water reactor (PWR) power plants [1]. Among all these PWR units, CPR1000 units are currently dominant in quantity in China, which are generation II+ pressurized water reactor, based on the design basic accident analyses under full range of operation conditions [2]. To further enhance CPR1000’s safety, after the accident of Fukushima, the Chinese government started a national project which targeted to promote CPR1000’s safety system. As a part of this work, the safety analysis of pressurizer is emphasized.
The pressurizer is the main equipment of PWR power plant pressure safety system [3] and the electric heaters are its central part. During actual operation, the electric heaters are subjected to ultra-high water pressure, transient extreme temperature changes, high radiation doses, and other severe tests. Once the electric heaters fail, it will affect the ability to maintain and control the operating pressure of reactor coolant system and even lead to the primary circuit overpressure shutdown [4, 5]. The power deviation of CPR1000 electric heater is plus or minus 5-6%. Such deviation will increase alternating heat and cold stress of the electric heater cladding during operation, making it easier for the electric heater cladding to expand and rupture, resulting in the casing to mechanically deform, increasing the probability of coolant leakage caused by a crevasse [6, 7]. And these cases of electric heater failure that have occurred in nuclear power plants in recent years have attracted more and more attention of scholars [8]. Therefore, on the basis of ensuring the quality of electric heating elements, it is also necessary to avoid the direct impact of inflow and outflow cold and hot fluids on the electric heating rods and to reduce the influence of velocity and temperature oscillations on the electric heater which generate during the mixing of cold and hot fluids. To prevent the repeated vibration, thermal fatigue failure and penetrating crack of the electric heater cladding and casing surface from damaging the welding seam between the electric heater and the casing, electric heater casing, and pressurizer. So it is very important to reveal the fluctuation mechanism of the mixing process in the pressurizer, which is of great significance to the specific arrangement of the electric heating element and the safe operation of CPR1000 unit.
In this paper, a full-scale three-dimensional model of CPR1000 pressurizer and its vertical surge line is established. Without considering the electric heating elements, the mixed fluctuation of cold and hot fluid is analyzed from the source. The hydraulic characteristics and heat transfer mechanism of the mixing of cold and hot fluids of pressurizer insurge are studied by 3D full-scale transient numerical simulation, comparing the fluctuations of three different surge line diameters, and the spatial distribution of time-averaged and RMS (root-mean-square) value of the velocity and temperature were obtained to analyze the strength of different position velocity and temperature fluctuation in fluid, in order to provide reference for the optimization and safe operation of electric heating elements in the actual project.
2. Physical Geometry and Numerical Modeling
2.1. Geometric Feature
From Table 1, it can be seen that CPR1000 has the smallest pressurizer volume, only half the volume of AP1000 pressurizer. The ratio of its volume to the reactor power is 1.5×10−2 m3/M; while the third generation nuclear power unit AP1000 is 1.75×10−2m3/M, EPR1000 is 1.67×10−2m3/M [9]. The operating environment of CPR1000 pressurizer is more severe than other types of pressurizers, so it is particularly important to analyze the transient fluctuation characteristics in the pressurizer.
Figure 1 shows the three-dimensional layout of CPR1000 pressurizer surge line. Based on the full demonstration and detailed assessment of the existing data, six locations (Section 1~Section 6) which are prone to fatigue damage are selected; the simulation condition of this paper is derived from the measured data of Section 1. The surge line connects the hot fluid in the pressurizer and the cold fluid in the heat pipe section to form a complete coolant circulation loop. When the coolant water temperature rises, the volume of coolant swells, leading to the phenomenon that coolant water fluctuates into the pressurizer through the surge line, which is called pressurizer insurge.
2.2. The Flow Domain for Calculation
In the present study, under the premise of satisfying the flow field analysis and calculation requirements, the actual pressurizer and its surge line were partially simplified. According to the design parameters of CPR1000, EPR1000, AP1000 pressurizer, and their surge line, three models which have the different pressurizer size and surge line diameter were established. The calculated area includes the entire pressurizer (upper head, barrel, and bottom head), the vertical part of surge line (above Section 1), and its length is . The geometric modeling and physical properties of medium are shown in Figure 2 and Table 2, respectively; it uses the center of the central cross-section of pressurizer barrel as the coordinate origin to establish a coordinate system. The flowing coolant inside surge line along y-axis and the direction of gravity are y-axis negative. is the total length of the pressurizer, is the inner diameter of barrel, and is the inner diameter of surge line; the specific parameters are shown in Table 1.
2.3. Simulation Conditions
Loss of offsite power design transient process is the expected transient condition of medium frequency during reactor operation. CPR1000 accident analysis assumes the loss of offsite power for all applicable events [10, 11]. Relevant tests and researches have shown that, in the initial stage of loss of offsite power, the peak velocity and flow rate of instantaneous coolant that flow into pressurizer are the largest in all transient conditions [12]. Meanwhile, the transient thermal fluctuation and velocity fluctuation of pressurizer are the most intense, which can easily lead to the instability of electric heater, thermal fatigue, and vibration of the support plate. Therefore, it is necessary to study the mixing process of cold and hot fluids in pressurizer when the loss of offsite power transient occurs. The full-scale three-dimensional transient simulation condition is derived from the thermal functional test data of the nuclear power plant, when pressurizer insurge occurs, the vertical surge line and half pressurizer are filled with thermal fluid, and the other 50% is steam space. The liquid phase and vapor phase are in saturated equilibrium, and the temperature of the steam and coolant is equal to the saturated temperature of the rated operating pressure (corresponding to 617.95K at 15.5 MPa). Suddenly, 595K low-temperature water from the surge line flows into pressurizer through the inlet of vertical surge line at a rate of 5.5m/s.
2.4. Experimental Setup
In order to obtain verification data for numerical simulation, experiment of mixing of hot and cold fluids in the pressurizer was carried out. The test section is made of transparent acrylic resin and has a water jacket around the pressurizer. The test section consists of a vertical pressurizer, thermal sleeve connecting to a vertical surge line with diameter of 0.366m, buffer tanks, high power electric heaters, high-pressure pump, liquidometer, pressure gage, and five thermocouples mainly. Schematic of the test section is shown in Figure 3 [13]. During the test, the pressure and temperature in the pressurizer rise to the rated parameters of the nuclear power plant by high-pressure pump and electric heaters (corresponding to 617.95K at 15.5 MPa).
The flow velocity distribution and the fluid temperature in the pressurizer are, respectively, measured by PIV and thermocouples in measurement planes. These thermocouples are installed at the central axis of y= -17DR, -16., -16DR, -15., and -15DR plane, respectively; detailed temperature field in the pressurizer can be measured. The diameter of thermocouple is 0.25 mm. The time response of these thermocouples is around 0.03 s and accuracy is less than 0.2°C by relative calibration. The temperature sampling frequency was 100 Hz and the measurement time period was 10 s. The PIV system consists of a double pulse YAG laser, a CCD camera, and a timing controller. The nylon powder of around 30 μm in diameter was used as tracer particle. The cross-correlation and subpixel methods [14] were used for the PIV data analysis. The spatial error of the correlation was about 0.2 pixels by using the subpixel method, and the velocity measurement error was less than 0.04m/s. The velocity vector fields of 256 sets were obtained with interval of 0.066 s (15 Hz). Influence of the refraction at the round pipe wall was reduced by the rounded water jacket.
3. Mathematical Models and Calculation Methods
3.1. Mesh Generation Methodologies
In order to improve the accuracy of calculation, to avoid the effect of pseudodiffusion and to reduce the amount of computation, the model adopts a hexahedron structured, block-based grid division. The standard wall function method was used to process near-wall flow. The height of the first mesh cell was determined by the Frank formula [15]. Boundary layer grids are set on the inner wall of bottom head and surge line, which are divided by o-block. The distance from the first layer to the wall is 0.05 mm, and the thickness growth coefficient of the boundary layer is 1.15. In consideration of relatively rapid change of flow field and temperature field at the joint of the bottom head and the surge line, finer meshes were distributed in the area by second cut and topological operation. The grid cells are getting smaller gradually and analyzing the result of each calculation. Compared with the previous calculation result, as two calculation results were within 5%, the grids will be adopted to ensure both accuracy and convergence [16–18], and 511250 hexahedral grid cells were generated for CPR 1000, 546596 for EPR1000, and 601942 for AP1000. The overall grid quality is between 0.6 and 1, and the aspect ratio is 1.2~18, which meets the requirement that maximum aspect ratio of three-dimensional heat transfer calculation is less than 35:1. The overall grid of pressurizer modeling and cross-section grid are presented in Figure 4.
(a) Overall grid
(b) Intersection grid
(c) Grid of surge line cross-section
3.2. Numerical Calculation Methods
Since there is an obvious vapor-liquid interface in the pressurizer, the simulation is based on RANS equation of VOF multiphase flow model. When doing VOF multiphase flow transient calculations, the time step selection is based on Global Courant Number. Finally, two time steps of 1ms and 5ms are selected, and Global Courant Number of each time step is below 2, which can completely satisfy the requirements of the simulation calculation [19]. Realizable k–ε Turbulence Model can eliminate the negative positive stress of Standard k–ε which may appear when the average strain rate is especially large by correlating the coefficient of turbulent dynamic viscosity formula and strain rate. Compared to the Standard k–ε, Realizable k–ε can accurately predict the diffusion of planar and circular jets, and it has a good simulation effect for the flow in the pipeline, flow separation, free flow of mixed flow, and the secondary flow [20, 21], so it was selected to deal with the turbulent flow of liquid phase. The finite volume method (FVM) is used to establish the solving control equation and the commercial CFD code ANSYS FLUENT is employed to solve the Reynolds-averaged Navier–Stokes equations with the Realizable k–ε turbulence model.
The focus of this study is not on the fluctuation of the vapor-liquid interface. Therefore, in order to save computing cost, the free liquid surface reconstruction adopts modified HRIC format. In addition, the PISO algorithm deals with the coupling of pressure and velocity, and the QUICK scheme was used to discrete each control equation.
In the mixed convection process of hot and cold fluids, buoyancy is one of the important analytical indicators. The flow strength induced by buoyancy can be determined by the ratio of buoyancy to inertia force, namely, the dimensionless Richardson number (Ri), and as shown in (1)-(3). It can be seen that dimensionless Gr and Re are measures of buoyancy lift and fluid inertia forces, respectively. When Ri approaches or exceeds 1.0, buoyancy has a great influence on flow [22].
where is the acceleration of gravity, 9.8 m/s2; β is the coefficient of thermal expansion of water, K−1; is the fluid viscosity, m2/s; is the average velocity of the fluid, m/s; is the characteristic length, m; is the hot fluid temperature, K; is the cold fluid temperature, K.
Considering the great temperature difference between the hot and cold fluids in this simulation condition, the variation of fluid density with temperature is quite large. If the Boussinesq hypothesis is used to approximate the buoyancy force, a large error will occur [23–25]. Therefore, this present simulation abandons the Boussinesq model and treats the density as a polynomial function of temperature to estimate the impact of the buoyancy term, as shown in Table 2.
3.3. The Dimensionless Time-Averaged and RMS Values
For the convenience and accuracy of analysis, the dimensionless time-averaged value was used to describe the average velocity and temperature, and the dimensionless RMS (root-mean-square) value was used to describe the fluctuation intensity in the mixing process. The dimensionless time-averaged and RMS temperature are as follows:
where is the dimensionless temperature.
In the same way, the dimensionless time-averaged and RMS velocities are defined as follows:
where is the velocity of insurge fluid.
In the preceding equation, N is the total number of collection points in the calculation time and is equal to 4400.
In this simulation, data were sampled along z-axis on the intersection lines of pressurizer vertical center section (x=0m) and y/ = -17, -16.5, -16, -15.5, and -15 sections, the sampling interval is 0.001s, and the total sampling time is from 0.6s to 5s.
3.4. Numerical Model Validation
In the present study, experimental research was adopted to understand the dynamic process of pressurizer insurge under the simulated condition parameters, and the instant temperature and velocity of the bottom head are obtained. Also the VBA function was used to extract and integrate data through writing the macrofile. In order to verify the validity and accuracy of the numerical model, five turbulence models were used to predict the flow and heat transfer process of the mixing of cold and hot fluids in AP1000 pressurizer, and the results were compared with the experimental data. In the unsteady state simulation, half pressurizer was filled with hot coolant at 617.95 K and the other 50% was steam space at the beginning of calculation. The surge line inlet was set as velocity inlet with 5.5 m/s at 595 K; meanwhile the pressurizer reference pressure was set to 15.4MPa. The ambient temperature outside pressurizer was set at 40°C and the heat transfer coefficient of the pressurizer wall was 0.01 W/(m2·K).
As shown in Figures 5 and 6, the vertical RMS velocity and temperature at the intersection line of the central section x= 0 and section y/ = -16 are plotted as samples, all five turbulence models can reflect the general trend of velocity and temperature fluctuation. The S-A (Spalart-Allmaras) model can simulate the fluctuation on both sides better, but the simulation result in the middle of pressurizer differ greatly from the test. The results of RNG k–ε and Realizable k–ε turbulence models are very similar, with differences only at the highest point of fluctuation. In general, the result with Realizable k–ε model is more close to the experimental data, so the model has been used in the subsequent calculations.
Figures 7 and 8 show that the distribution of RMS velocity and temperature, respectively, in the simulation and experiment. It can be seen from the figures that RMS velocities and temperature from the simulation agree well with the experiment, although there is a certain deviation between the two. In Figure 7, the average error of velocity fluctuation intensity between three simulation curves and measured value is 4.6% and the maximum error is 14.8%. In Figure 8, the average error of temperature fluctuation intensity between calculated and measured value is 2.3% and the maximum error is 6%. The error is mainly concentrated in the geometric center of the pressurizer. Compared with the experimental data, the numerical results are within the acceptable error range of 15%, which indicates that Realizable k–ε combined with VOF multiphase flow heat transfer model is reliable and economical for predicting the complex turbulent mixing process of cold and hot fluids in the pressurizer.
Figures 9 and 10 show the comparisons of dimensionless RMS velocity and temperature of y-direction at the intersection line of the central section x= 0 and section y/= -16 at different time steps between experimental results and simulated results. The error between numerical simulation results and experimental values will decrease as the time step decreases. Therefore, in the following analyses, 1ms is selected as the calculation time step.
4. Results and Analyses
According to the velocity comparison between three different surge line diameters of typical section lines at the vertical center cross-section (x=0) of the pressurizer (Figure 11), it can be seen that the trend of different diameters is basically consistent, the middle is high, and two sides are low. As the diameter of surge line increases, the peak value of dimensionless RMS velocity also decreases, indicating that when = 0.366, the fluctuation degree of velocity is the smallest. In other words, the velocity fluctuation in CPR1000 pressurizer is the most intense, which increases the possibility of pulsating water directly rushing toward the vapor-liquid interface. Next, the mixed fluctuation characteristics of hot and cold fluids of pressurizer insurge in CPR1000 pressurizer ( = 0.284) are analyzed in detail.
(a) y/= -16.5
(b) /= -16
4.1. Mesh Sensitivity Study
To verify that the calculation results are independent of mesh number in this paper, a mesh sensitivity study has been performed by refining the meshes in specific portions of the domain. Mesh 4 with 511250 cells was chosen as the reference mesh. Mesh 5 with 531940 cells was obtained by refining Mesh 4 in regions of surge line and bottom head. Mesh 6 with 639835 cells was generated by further refining Mesh 5 in the same regions by decreasing the minimum and maximum cell size parameters by half. The major differences between finer and coarser meshes are boundary layer number increased from 10 to 25 layers; the base cell size. The number of different types of meshes is shown in Table 3. The mesh sensitivity study has been conducted by comparing the distributions of dimensionless RMS velocity and temperature along y/= -16 intersection line using different meshes, as shown in Figures 12 and 13. It can be seen that very little improvement was obtained using finer meshes (Meshes 5 and 6). Thus, Mesh 4 was chosen as the final mesh and used for the prediction of mixing phenomena in CPR1000 pressurizer in this study.
4.2. Velocity Analysis
The y-velocity flow field contours of the plane x=0m at different times are presented in Figures 14(a)–14(f). As can be seen from the figures, the high-speed cold fluid in the surge line is vertically rushed into the pressurizer and complete longitudinal circulations appear on the left and right sides of the water tongue which are formed under the combination of inertial force and buoyancy lift in the mixing process. At the time of 0.75 s, the center vortex of the circulation close to the wall of bottom head, the junction of the bottom head, and surge line will have frequent thermal shock, so thermal sleeve often is provided at the nozzle to reduce thermal stress [26]. As the time series progresses, the water flow spirals upward, and the central vortex of the longitudinal circulation is offset upwards from the bottom head wall, and the circulation gradually occupies the bottom head.
It can be seen from Figure 14 that the velocity of inflow cold fluid has a tendency to decrease after a short period of acceleration at the exit of surge line, but a stable accelerating “elastic zone” would be formed gradually at the end of main axis region. This is because as the cold fluid and the surrounding water are sucked and mixed, the energy of insurge is transmitted quickly to the surrounding medium, so the velocity gradually diffuses and decays toward the surroundings. At the end of each moment, the Ri number increases significantly due to the velocity difference of interface between hot and cold fluids at the end of the water tongue decreases sharply, which leads to buoyancy lift has a greater impact on the fluid than inertial forces. In addition, due to the low density of thermal fluid in the pressurizer, the up-down convection of hot and cold fluids occurs under the effect of buoyancy lift, and the cold fluid tends to fluctuate upward, resulting in the formation of the accelerated zone with maximum velocity of up to 6 m/s.
Figure 15 shows dimensionless time-averaged velocity of y-direction at the x=0m plane. It can be seen that the time-averaged velocity of y-direction increases sharply and the gradient of change is large in the range of -0.15m<z<0.15m. As a whole, the time-averaged velocity in the y-direction of different section lines is consistent along z-axis and the velocity decreases from the middle to both sides. But the distribution trend is more gradual with the increase of y, indicating that, with the advancement of the inflow process, the kinetic energy of insurge is consumed gradually, so the velocity becomes more uniform and the fluctuation is mainly concentrated in the bottom head. In the Z>0.2m and Z<-0.2m areas, the y-direction time-averaged velocity is negative, indicating that there is recirculation in the areas; as the y increases, the starting point of the recirculation zone gradually expands to both sides, inflow fluid in the upper part of the vortex is accelerated to a higher velocity than the main flow velocity, and this trend corresponds to the distribution characteristics of the velocity flow field in Figure 14, further illustrating that the cold and hot fluids form a significant vortex during the mixing process.
The dimensionless time-averaged velocity can describe the average velocity of pressurizer insurge process, but it cannot reflect the intensity of velocity fluctuations over a period of time. Figure 16 shows dimensionless RMS velocity of y-direction in plane x= 0m. It can be seen that dimensionless RMS velocity of -0.25m<z<0.25m is relatively large, indicating that y-direction velocity fluctuation in this region is more severe, so the arrangement of the electric heating elements within the pressurizer should avoid this area, thereby reducing the impact of fluid confluence on them. On the other hand, the RMS velocity of z<-0.45m and z>0.45m is small, indicating that the velocity fluctuation intensity in these areas is weak. The RMS velocity in the y-direction of different section lines is consistent along z-axis, but the distribution trend becomes steeper with the increase of y gradually, indicating that, with the advancement of the inflow process, the velocity fluctuation is getting bigger and bigger at the bottom head. This is because, as the jet flows away from the surge line, the kinetic energy of insurge can be gradually absorbed by the surrounding fluid. In addition, the oppression of original thermal fluid can make the velocity of hot fluid away from the surge line not change rapidly; instead, it fluctuates repeatedly under the influence of gravity and buoyancy, so the degree of fluctuation will increase.
4.3. Temperature Analysis
Figure 17 shows the distribution of dimensionless time-averaged temperature at different positions of the center section of the pressurizer along z-axis. In the range of z<-0.4m and z>0.4m, the dimensionless time-averaged temperature at all positions remains the same, which is the initial dimensionless temperature, indicating that there is no mixing effect of cold and hot fluid in these areas. In general, the distribution trend of each section line’s dimensionless temperature along z-axis is approximately the same, and there are obvious inflection points at the center of geometry, indicating that hot and cold fluids in the area are mixed violently. At y/ = -15.5 and -15 far from the surge line, the dimensionless time-averaged temperature along z-direction gradually flattens out with the increase of y, and the trough of dimensionless temperature rises gradually. The range in which dimensionless temperature is 1 is also shrinking, further confirming that the energy of insurge gradually diffuses into the surrounding medium, resulting in a gradual expansion of the mixing range of cold and hot fluids.
The dimensionless RMS temperature can describes the temperature fluctuation intensity, which is similar to dimensionless RMS velocity which reflects velocity fluctuation intensity. As can be seen in Figure 18, there are differences in the fluctuation range and peaks of the z-direction temperature fluctuation curve of each section line. In the range of z<-0.45m and z>0.45m, the dimensionless RMS temperature is zero, indicating that the area is completely filled with hot fluid; in other words, there is no temperature fluctuation in the area. With y position goes up, the hot fluctuation area shows a tendency of expansion obviously; the distribution of dimensionless RMS temperature along z-axis becomes more and more convex, and the degree of overall thermal fluctuation increases; the mixing process occurred in the entire area of -0.45m < z < 0.45m. In the range of -0.2m<z<0.2m, the mixing of hot and cold fluids is most intense; the hot fluctuations at y/ = -15 and -15.5 become more severe, further confirming that, as the advancement of inflow process, the energy of insurge is gradually absorbed by the surrounding medium; the mixing away from the surge line mainly depends on the effects of gravity and buoyancy. From the above analyses, it should be avoided that the electric heaters are installed in a relatively intense temperature fluctuation range of −0.45m<z<0.45m to reduce the degree of thermal stress occurring on the casing wall and the probability of fatigue failure of pipe material.
5. Conclusions
This paper selects the maximum inflow velocity of pressurizer insurge during the initial period of loss of offsite power as the experimental and simulation conditions. The CFD method of multiphase flow transient heat transfer was adopted to perform three-dimensional numerical simulation on the mixing process of cold and hot fluids to obtain the instantaneous velocity field and the distributions and variations of temperature. By defining the dimensionless time-averaged and RMS parameters to reflect the fluctuations of velocity and temperature over time at different cross-sections lines, the whole transient process was analyzed and the following conclusions are drawn:
There is certain difference between the numerical simulation results and experimental results, but the overall trend is in good agreement, indicating that Realizable k–ε turbulence model combined with VOF multiphase flow heat transfer model can effectively simulate the flow and heat transfer process of cold and hot fluids mixing in the pressurizer.
The symmetrically distributed longitudinal circulations are formed on both sides of the wave inflow axis, which can enhance the mixing and entrainment effect. As the transient process proceeds, the circulation will occupy the entire bottom head. Ri has great influence on the mixing of cold and hot fluids; its size determines the magnitude of the effect of buoyancy lift on the fluid. Under the action of buoyancy lift, a stable accelerating “elastic zone” will be formed gradually at the end of water tongue.
The mixing of hot and cold fluids occurs mainly in the area of -0.45m<z<0.45m, in which there is a clear backflow. With the upward movement of y, the starting point of backflow area gradually extends to both sides, so the backflow is mainly concentrated in the bottom head. The most intense fluctuations of temperature and velocity are concentrated in -0.2m<z<0.2m, and the peak value of fluctuation gradually increases with the advancement of inflow process.
In order to ensure the reliability of the electric heating element and reduce the influence of velocity and thermal oscillations on the electric heater generated during the mixing of the hot and cold fluids, the electric heater within CPR1000 pressurizer should be arranged in the z<-0.45m and z>0.45m as much as possible.
Data Availability
The data used to support the findings of this study are included within the article. All data included in this study are available upon request by contact with the corresponding author.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by China National Nuclear Corporation (CNNC) Thermal Hydraulic Key Laboratory Fund Project (20130901).