Abstract
The graphite oxidation of fuel element has obtained high attention in air ingress accident analysis of high temperature gas-cooled reactor (HTR). The shape function, defined as the relationship between the maximum and the average of the oxidation, is an important factor to estimate the consequence of the accident. There are no detailed studies on the shape function currently except two experiments several decades ago. With the development of computer technology, CFD method is used in the numerical experiment about graphite oxidation in pebble bed of HTR in this paper. Structured packed beds are used in the calculation instead of random packed beds. The result shows the nonuniform distribution of oxidation on the sphere surface and the shape function in the condition of air ingress accident. Furthermore, the sensitive factors of shape function, such as temperature and Re number, are discussed in detail and the relationship between the shape function and sensitive factors is explained. According to the results in this paper, the shape function ranges from 1.05 to 4.7 under the condition of temperature varying from 600°C to 1200°C and Re varying from 16 to 1600.
1. Introduction
Pebble-bed reactor is a typical design of high temperature gas-cooled reactor (HTR), including HTR-MODULE, PBMR, and HTR-PM. Spherical fuels are used in these HTRs and were packed in the core randomly. In case of HTR-PM, there are 420,000 fuel elements with a diameter of 6 cm in the equilibrium core. The average filling rate of the core is about 0.61. Each fuel element has two zones, where the inside one is the fuel zone of mixing TRISO and matrix graphite and the outside is the fuel-free zone with pure graphite.
During the air ingress accident, which is a specific beyond design basis accident for HTRs, air will enter the reactor core and react with the graphite of fuel elements. System codes are utilized for the oxidation calculation in the core, such as THERMIX/REACT, TINTE, and GAMMA. Porous model is built in these codes to simplify the complex phenomenon in the pebble bed and shows a good agreement with the experiments in terms of flow and heat transfer. The oxidation result in porous model represents an average value in a mesh or on one fuel pellet. However, for most cases, the oxidation does not proceed homogeneously among the surface, because the flow field is not homogeneous around the fuel ball. So the result may be nonconservative considering this phenomenon. The so-called shape function is defined as the relation between the maximum and the average burn-off of a pebble. According to Moormann and Hinssen H’s research based on SUPERNOVA experiment (Juelich Research center), the shape function is measured to be about 3 [1], as shown in Figure 1. The experiment is carried out with high flow rate (Re ≈ 1600) and high temperature ( K) and the total burn-off is 15%~25%. Unlike Moormann and Hinssen H’s result, another experiment is held in the facility VELUNA (University of Duisburg). The shape function measured in this experiment is about 9 under lower flow rate [2]. Due to the lack of more experimental research, the factor between 3 and 9 is recommended to estimate the heterogeneous oxidation on the pebble in paper [1]. There is no further research for the nonuniform oxidation after that.
Due to the rapid development of computer technology, computational fluid dynamics (CFD) become quite popular and effective for solving complex problem including pebble bed. With the method of CFD, detailed studies were performed for the flow and heat transport characteristics in the pebbles during the last decade. For example, Dixon and Nijemeisland performed a series of numerical researches on the flow and heat transfer processes in the fixed-bed reactors [3–5]. The velocity and temperature details inside the pores were obtained and the mass and heat transfer performance were carefully analyzed. Yang and his team conducted experiments and CFD simulation for the pressure drop and heat transfer in structured packed beds [6–8]. Typical structural models were built and computed with the code CFX. Similar to the pervious study for flow and heat transport, oxidation on the surface of fuel elements is simulated with CFD method in this paper. Instead of randomly packed bed, structured packed bed is used to get the representative result and simplify the calculation. Three typical stacking modes are in consideration, which are FCC (face-centered cubic), BCC (body center cubic), and SC (simple cubic). Their porosity is 0.2595, 0.3202, and 0.4764, respectively, covering the range of randomly packed bed.
2. Physical and Computing Model
2.1. Oxidation Model
During the air ingress accident, the graphite will react with the oxygen in the core. Ignoring other oxidizing gases in the air, two reactions are considered in the air ingress accident:
The reaction rate of reactions (1)~(2) is usually written in the form of (3)~(4), which shows the reaction rate of graphite and the ratio of the two reactions, respectively,
The matrix graphite of fuel element is A3-3. Based on the experiment of Juelich Research Center, reaction rate of graphite is shown as (5) [9]. This formula is widely used in system codes of HTR, such as TINTE and THERMIX/REACT. For the lack of suitable experiment, the ratio fitted by experiment of IG-110, another kind of nuclear grade graphite, is used for the calculation of A3-3 [10]
Reaction rate is highly dependent on the graphite temperature and partial pressure of oxygen in the gas flow. As discussed by some researchers, the reaction mechanism can be divided into three zones according to the temperature. At low temperature (zone I), which is less than 650°C, the reaction rate is so small that oxygen consumption does not change the oxygen distribution in the gas stream. In this zone, temperature is the mean factor affecting the reaction rate. At high temperature (zone III), which is higher than 750°C, the reaction rate is too large that the concentration of oxygen near the wall is much lower than the mainstream because of the consumption. In this zone, reaction rate is restricted by the speed of mass transfer near the pellets surface according to the theory of boundary layer.
In the air ingress accident of HTR-PM, the temperature of most fuel elements is higher than 750°C, and the region of interest, where the maximum oxidation occurs in the core, is oxidized in zone III [11].
2.2. Computing Model
Structured packed bed is used to simulate the oxidation of fuel elements in the core. Three typical arrangements are selected, including simple cubic (SC), body center cubic (BCC), and face center cubic (FCC), as shown in Figure 2. Referring to the SUPERNOVA experiment, only one sphere reacts with oxygen (one graphite sphere only and aluminum for other spheres) in calculation zone. In order to eliminate the effects of entrance and outlet, three cells are connected in series along the flow direction. According to the symmetry of the structure, 1/4 model of the cells are selected for the numerical calculation. Three-dimensional geometry model and the meshing are established by code ANSYS 16.0.
It is difficult to generate high-quality computational grids near the contact point between the spheres, because it is highly skewed. Several methods are listed to modify the contact point in Bu et al.’s paper [8]. As the first typical modification, the diameter of spheres reduces a little and there will be a small gap between the two spheres instead of the contact point. Another method is modified as follows: enlarge the diameter and there will be an overlap replacing the contact point. These two methods may change the porosity of the packed bed and affect the flow field in the pore. An effective modification is named bridge. A cylindrical bridge connects the two spheres, making the space near the contact point less skewed. The accuracy can be adjusted by changing the diameter of the bridge with this method. Opposite to bridge, if the cylinder is cut off at the contact point, there will be gap with constant distance between the two spheres. This modification is called cap. The latter two methods do not change the diameter of spheres, so the porosity of CFD model is more similar to the real one. The bridge method is recommended in Bu et al.’s paper since its result is more close to the experimental data in local flow and temperature distributions. Bridge method is used in this paper and the mesh near the contact point is shown in Figure 3.
The grid-independence is checked before the calculation. The fluid zone is filled with unstructured meshes and inflation meshes are placed on the surface of spheres and bridges. In the case of BCC, three sets of grids with total cell number of 817235, 1309163, and 1810201 are used for the grid-independence test. The total amount of oxidation on the graphite sphere, as well as the maximum of reaction rate on the surface, is compared among the cases, as shown in Table 1. The results of oxidation distribution of different grids are nearly the same and the difference between fine and coarse grids is less than 0.5%. Richardson extrapolation is used to calculate the grid convergence value, estimating the real result. It is clearly indicated that the results of fine grid are accurate enough with the error less than 0.2%. The fine grid of BCC is used in the further calculation and the same grid setting is used for FCC and SC model in this paper.
With the modification of bridge, reaction area on graphite sphere in computing model is smaller than that in real geometry. The difference is determined by the diameter of cylinder. The diameter of graphite sphere is 6 cm in HTR-PM. The cases with different bridge diameters, varying from 6 mm to 12 mm, are discussed in this section, whose results are listed in Table 2. The area loss with different bridges varies from 2% to 8.15%. However, the oxidation results, the integral and the max value, have little difference with each other. So, the diameter of bridge has only small influence on the distribution of oxidation. The bridge diameter of 8 mm is used in this paper.
The setting of CFD is listed below, which is verified in paper [11]: CFD code: ANSYS FLUENT 16.0. Solver: pressure-based. Viscous: SST -omega. Coupling scheme: SIMPLEC. Discretization: second-order upwind.
3. Result and Discussion
3.1. Oxidation Distribution on the Sphere Surface
BCC packing form is used as a basic case, because its porosity is similar to that of pebble-bed reactor core. According to the research of air ingress accident, the mass flow rate of stable air ingress is about 0.1 kg/s and the temperature of pebbles remains high relatively [12]. Therefore, low mass flow rate (Re ≈ 16) and high temperature (°C) are suitable for the numerical experiments.
Ignoring the change of shape, the distribution of reaction rate is the same as the distribution of oxidation with low burn-off. Figure 4 shows the oxidation distribution on the sphere surface. It is obvious that oxidation on the surface is uneven seriously. In the case of BCC arrangement, oxidation mainly locates in the flow regions far away from the contact point in the front of the sphere. It has lower reaction rate near the contact point especially in the tail region behind the bridge. This result is similar to the phenomenon that no oxidation is found at the contact point up to a total burn-off of 20% in SUPERNOVA experiment.
(a) Geometry of flow path
(b) Oxidation rate
Based on the theory of graphite oxidation, reaction rate is affected by temperature and partial pressure of oxygen on the surface. With the constant temperature in this case, the distribution of reaction rate is determined by oxygen partial pressure which is affected by two factors. The first factor is the oxygen concentration in the main flow. Oxygen consumes rapidly when flowing though the graphite surface, causing the decreasing of oxygen concentration along the flow direction, as shown in Figure 5. Therefore, the oxidation on the windward side is larger than that on leeward side. According to the structure of BCC, mass flow rate near the edge of the model is much large than other places (Figure 6). Since the adequate supply of oxygen, oxidation concentrates in this area. On the contrary, the oxidation in the zone behind the bridges is weaker since it has lower flow rate.
Another important factor is mass transfer coefficient near the graphite surface. Concentration boundary layer is shown clearly in Figure 5 and it is affected by local Re significantly. Therefore, the thickness of boundary layer on leeward side is much thicker than that in windward, leading to lower reaction rate there. Stagnation zone existing in the forefront of sphere has lower Re since the blocking by the front sphere. The oxidation is lower compared to the zone nearby, although the oxygen concentration of main flow there is the highest.
The maximum and the average oxidation on the sphere surface will be calculated in this section. For the modification of contact point, the area of computation model is slightly less than the real one. The average value will be modified by the sphere area. In this case, the max and average oxidation are, respectively, kg/(m2·s) and kg/(m2·s). The shape function is 1.95 according to its definition.
With the analysis of the result, the distribution of oxidation is affected by some factors, such as Re, temperature, and arrangement of spheres. So the relationships between them will be discussed in the next sections.
3.2. The Influence of Re
Re number in pebble bed is determined by mass flow rate. Figure 7 shows the distribution of reaction rate with different Re numbers. The result in former section shows the oxygen consumption in low Re. With Re increasing, more oxygen is supplied in the air flow. Besides, mass transfer is enhanced near the surface of graphite with the strengthening of turbulence. So the oxidation rate, both the maximum and the average value, gets stronger with the increase of Re number (Figure 8(a)). Furthermore, turbulence is strengthened with the increase of Re. With the influence of turbulence, the stagnation zone in the forefront of the sphere becomes smaller. The vortex behind the bridge is also strengthened, enhancing the mixing process in the tail region. The region with low oxidation gets smaller and the region of irregular oxidation gets bigger due to the vortex with the Re increasing.
(a) Average and max oxidation
(b) Shape function
Overall, maximum and average oxidation will increase with the increase of Re. As a result of the quantitative calculation, the relationship between shape function and Re number is shown in Figure 8(b). The shape function decreases with the Re number increasing and shows a linear relationship with logarithmic Re. The data of SUPERNOVA and VELUNA experiments shows a similar trend to this simulation, although the result of VELUNA experiment is larger than that in the simulation.
3.3. The Influence of Solid Temperature
Temperature is another factor that influences oxidation rate significantly. Oxidation rate is shown in Figure 9 with different temperatures from 600°C to 1200°C. The reaction rate is particularly sensitive to temperature. At low temperature, the reaction rete is relatively low as well. The oxygen consumption is too small to change the concentration in the main flow. Thus, the oxidation rate on the sphere surface is uniform relatively. As the temperature increases, the supply of oxygen becomes inadequate, leading to the inhomogeneity of oxidation rate. After it reaches the transition temperature, reaction rate is limited by oxygen concentration on the surface seriously and the growth rate of oxidation becomes much lower with the temperature increasing. This trend is explained by Arrhenius’s theory. In the case of Re = 16, the transition temperature may be some value between 800°C and 1000°C.
Shape function is calculated at different temperatures and different Re numbers. The relationship between shape function and temperature is shown in Figure 10. The shape function will increase monotonically as the temperature increases with all Re values. However, the transition temperature is different with different Re numbers. For higher Re, the limitation of temperature will be higher relatively. In the case of high temperature, the relationship with Re number shows different discipline above the transition point.
3.4. Results of Different Arrangements
Three typical arrangements are involved in this research. The result of BCC structure is explained above in detail. In this section, the results and trends of FCC and BC are introduced. The distribution of reaction rate is shown in Figure 11, which shows similar discipline to BCC arrangement. The oxidation on windward side is much larger than that on leeward side. And the oxidation concentrates in some certain regions according to the contacts of spheres. The locations of those regions are different between FCC and SC. The oxidation appears to be largest in the forefront of graphite sphere in FCC, while that region in SC is on the windward side of the sphere. The regions with lowest oxidation locate behind the bridge for both cases, the same as BCC.
(a) FCC
(b) SC
There are many similar results between different models. Firstly, the relationship between shape function and temperature is the same as the model of BCC. And at low temperature, the discipline between Re and shape function is very similar (Figure 12).
(a) FCC
(b) SC
Also, there are some differences between the results of the three models. The sequence of the three models by complexity of flow path is FCC > BCC > SC. Concentration layer is affected more seriously by Re in complex channel. So the Re value has more influence in FCC model. The porosity of the models is 0.2595, 0.3202, and 0.4764, respectively. According to the calculation results of FCC, BCC, and SC, the shape function is increasing as the porosity increases, especially at high temperature. However, the amount of samples in this research is not large enough to get an accurate conclusion. More calculations or experiments with random-pebbles should be carried out to verify this prediction
4. Conclusion
The distribution of oxidation on the surface of fuel element in air ingress accident is studied in this paper. CFD method is used to simulate the reaction rate by utilizing the code ANSYS 16.0. Structured packed beds (BCC, FCC, and SC) are adopted to model the basic cell of pebble beds, instead of random packed beds. Oxidation rate and shape function (the relationship between maximum and average of oxidation on the sphere surface) are calculated with the model of BCC. Furthermore, the sensitive factors of shape function, such as Re, temperature, and arrangement, are discussed in detail and the result shows agreement with SUPERNOVA experiment. The main conclusions in this paper are as below.
The distribution of oxidation rate is nonuniform on the surface of element in air ingress accident. The region with highest oxidation always locates on the windward side and far away from the contact point. On the other side, the region with lowest oxidation locates near or behind the contact point.
Re and temperature are sensitive factors of shape function. The shape function always increases with the solid temperature increasing. The growth rate decreases obviously when temperature reaches a transition point. At low temperature, the shape function decreases with Re increasing. However, the discipline will change when the temperature is higher than the transition point.
Different arrangements are involved in the calculation. The discipline of oxidation distribution is similar in the model of FCC, BCC, and SC. The shape function ranges from 1.05 to 4.7 in the condition of temperature varying from 600°C to 1200°C and Re varying from 16 to 1600.
Nomenclature
| : | Ratio of two reactions |
| : | Partial pressure of oxygen (Pa) |
| : | Reaction rate (kg/(m2 s)) |
| : | 8.314 (J/(K mol)) |
| : | Reynolds number |
| : | Solid temperature (K). |
Additional Points
Two aspects need further research in this paper. Firstly, the calculations with structured packed beds are not insufficient for the research of the relationship between shape function and arrangement. So model of random packed beds will be built with DEM method and be used in the oxidation calculation in the next step. Moreover, the computation method in this research is fit for oxidation on fuel element surface with low burn-off, because the shape of sphere remains the same in the oxidation process. To simulate the oxidation process with high burn-off accurately, transient calculation is required. Multiphase flow and dynamic mesh may be used to describe the changes of shape.
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
This paper is supported by the National S&T Major Project (Grant no. ZX06901) and the National Natural Science Foundation of China (Grant no. 51306100).