Multiphysics Investigation on Coolant Thermohydraulic Conditions and Fuel Rod Behavior During a Loss-of-Coolant Accident

This article simulates the multiphysics coolant thermohydraulic conditions and fuel performance of a pressurized water reactor (PWR) during a loss-of-coolant accident (LOCA). In the coolant channel of a PWR, the coolant undergoes a series of different boiling regimes along the axial direction. At the inlet of the coolant channel, heat exchange between the cladding wall and coolant is based on single-phase forced convection. As the coolant flow distance increases, the boiling regime gradually converts to nucleate boiling. When a LOCA occurs, on the one hand, the coolant flux and coolant pressure decrease sharply; on the other hand, the heat flux at the cladding wall decreases relatively slowly. They both contribute to a swift increase in coolant temperature. As a consequence, a boiling crisis may occur as critical heat flux (CHF) decreases. In this article, the void fraction along the length of coolant channel in a reactor and mechanical performance of Zr cladding enwrapping UO2 fuel are investigated by establishing a fully coupled multiphysics model based on the CAMPUS code. Physical models of coolant boiling regimes are implemented into the CAMPUS code by adopting different heat transfer models and void fraction models. Physical properties of the coolant are implemented into the CAMPUS code using curve-fitting results. All physical models and parameters related to solid heat transfer are implemented into the CAMPUS code with a 2D axisymmetric geometry. The modeling results help enhance our understanding of void fraction along the length of the coolant channel and mechanical performance of Zr cladding enwrapping UO2 fuel under different coolant pressure and mass flux conditions during a LOCA.


INTRODUCTION
Developing a computational code fully coupling multi-physical fields in the pressurized water reactor (PWR) is of great interest. It is essential to the simulation of fuel performance under a LOCA, which makes a significant difference to both the safety and economic effect of a reactor. In the past, several codes have been developed to simulate separately normal operation fuel behaviors (e.g., FRAPCON, ENIGMA, COMETHE, and FEMAXI) and accident operation behaviors in a reactor (e.g., FRAPTRAN, TRANSURANUS, SCANAIR, TESPA-ROD, and SFPR). In recent years, some of these codes in both categories tend to develop toward a unified code that could simulate both normal operation and accident conditions (Van Uffelen et al., 2019). However, most of these codes require the implementation of specific models and have an increased complexity level, which requires reduction (Van Uffelen et al., 2008;Van Uffelen et al., 2019;Pastore et al., 2021). Therefore, a simplified thermo-fluid-coupled code involving key fuel performance phenomena and continuous flow and phase change of coolant is highly desired. This is beneficial for the more realistic simulation under a LOCA.
A LOCA results from the break on the pipes of the loop in a reactor. Once a LOCA happens, the core scrams, and a drop of pressure and flux in the circulation loop occurs without appropriate coolant supplements. In the condition of depressurization, the flow of the coolant might suffer from a transition from single-phase forced convection to boiling regime and even encounter critical heat flux (CHF). The appearance of CHF results in the increase of temperature at the outer surface of cladding, which leads to fuel cracking, pellet-cladding interaction, and the release of fission product gases (Belle, 1961;Holden, 1966;Frost, 1982;Bailly et al., 1999;Liu et al., 2016a).
We investigated the void fraction of the coolant and the mechanical properties of fuel cladding before the flow boiling reaches CHF in this article, fully coupling key fuel performance phenomena, and cladding-coolant heat transfer. A set of heat transfer correlations were adopted to describe the heat transfer condition prior to CHF. In the regime of single-phase forced convection, the most widely used correlation is the Dittus-Belter correlation (Dittus and Belter, 1985); it is applicable when the coolant is still in the liquid phase. In the regime of subcooled boiling, many studies were also carried out to predict the heat transfer coefficient for different geometries, flow, and heating conditions; among which, the most widely used are those proposed by Bergles and Rohsenow (Bergles and Rohsenow, 1964), Jens and Lottes (Jens and Lottes, 1951), Thom et al. (Thom et al., 1965), and Shah (Shah, 1977;Shah, 2017). In the saturated boiling region, Schrock-Grossman correlation (Schrock and Grossman, 1962) and Chen's correlation (Chen, 1966) have been applied in a large range in the engineering calculation. In terms of void fraction, Levy (Chen, 1967), Kroeger and Zuber (Kroeger and Zuber, 1968), Saha and Zuber (Saha and Zuber, 1974), and Lahey and Moody (Lahey and Moody, 1977) proposed different models to estimate the void fraction in subcooled boiling regime based on the Zuber and Findlay (Zuber and Findlay, 1965) drift flux model. Later, Manon (Manon, 2000) modified Lahey and Moody's model (Lahey and Moody, 1977) by combining the Griffith et al. (1958) model. In this work, the model was established based on the CityU Advanced Multiphysics Nuclear Fuels Performance with Userdefined Simulation (CAMPUS) code. The CAMPUS code is based on the framework of COMSOL Multiphysics, which is, in general, a finite element analysis solver and simulation software. Previous results (Liu et al., 2016a) calculated by the CAMPUS code using the constant heat transfer coefficient has been proven reliable by comparing with the results calculated using BISON, ABAQUS, and FRAPCON. The modeling results of a previous work were adopted and combined with different heat transfer models and void fraction models to predict the void fraction along the length of the coolant channel in a reactor and mechanical performance of Zr cladding enwrapping UO 2 fuel under normal operating conditions and conditions with a drop of pressure and coolant mass flux. Our understanding is toward void fraction along the length of coolant channel and mechanical performance of Zr cladding enwrapping UO 2 fuel cladding under different coolant pressures and mass fluxes through fully coupled modeling including a LOCA. This is useful for the prediction of CHF and the safety of cladding material.

Model Geometry
The model used in this work is established based on a 2D axisymmetric geometry. The whole fuel rod is represented with a single pellet by applying periodic boundary conditions, as is shown in Figure 1 (Liu et al., 2016a;Liu et al., 2016b;Liu and Zhou, 2017;Liu et al., 2018). The previous study focused on the fuel performance of thermal conductivity-enhanced material in light water reactor has been done in Ref. (Liu et al., 2016a;Liu et al., 2016b;Liu and Zhou, 2017;Liu et al., 2018). They used a self-defined multiple physics model based on COMSOL to simulate the fuel performance of UO 2 -BeO and UO 2 -SiC. Almost all physical phenomena related to actual fuel rods were taken into consideration, which included heat generation and conduction, mechanical properties, fission gas release, cladding irradiation creep, and oxidation. However, the coolant properties in the previous work are not sophisticated enough. As a result, in this work, we tried to build a thermalfluid-coupled multiphysics model to estimate the fuel performance by adding to the previous model, more complex fluid properties, and heat conduction on boundaries. The assemblage geometry was assumed to be rectangular, as shown in Figure 2. Therefore, the hydraulic diameter was calculated by the equation: In our model, a linear average power of 200 W/cm was assumed to be reached within 1,000 s, and lasted over 3 years. Figure 3 shows the evolution of wall temperature under 15.5 MPa. Obviously, the wall temperature increases sharply after 10 7 s because of the rapid increase of zirconium oxide layer thickness. This leads to a rapid diminution of the coolant heat transfer coefficient. However, the increase of zirconium oxide layer thickness is mild enough between 10 4 and 10 7 s. As a result, we could consider that the evolution of wall temperature during this time is stationary so that different sections of the pellet on the fuel rod could be represented by changing the inlet coolant temperature on a single pellet. Some of the calculated results were therefore averaged between 10 4 and 10 7 s. Besides, the position of the pellet on the fuel rod was associated with coolant temperature by an energy balance equation followed assuming that q ″ was maintained constant along the height of the fuel rod.
where the inlet temperature was taken as T in 566 K and z in 0 m.

Properties of the Coolant
A coolant in this work is assumed as water. Several physical properties of water were needed for the modeling of heat transfer between fuel cladding and coolant. Properties of both liquid and vapor of water as a function of temperature and pressure were implemented in the model using the curve fitting method.

Density
The density of the liquid phase was assumed to be in the form of ρ l a ρ l T 3 + b ρ l T 2 + c ρ l T + d ρ l . The coefficients a, b, c, and d are different under different coolant pressures. The fitting results are valid for a pressure range between 15.5 and 10 MPa and temperature between 293.15 K and T sat . The density of saturated vapor took the form ρ g a ρ g exp(b ρ g T sat ). Then ρ g was applied to an unsaturated coolant temperature to obtain the corresponding density of unsaturated vapor. The fitting results are valid between 293.15 K and critical temperature. Values of constants are given in Table 1. Figure 4 gives the figure representation of density for both liquid phase and saturated vapor.

Thermal Conductivity
Thermal conductivity of the liquid phase was assumed to be in the form of k l a k l T 2 + b k l T + c k l . The values of constants a, b, and c are different under different coolant pressures. The fitting results are valid for pressure between 15.5 and 7 MPa and temperature between 293.15 K and T sat . Values of constants are given in Table 1. Figure 5A gives the figure representation of curve fitting results under 15.5 MPa.

Specific Heat Capacity
We have carried out curve fitting of specific heat capacity for the liquid phase between 15.5 and 10 MPa, 293.15 K and T sat , which took the following form: c p,l a c p,l + b c p,l T/1 + c c p,l T + d c p,l T 2 . Values of coefficients are given in Table 1, and the fitting results are presented in Figure 5B.

Viscosity
Curve fitting of dynamic viscosity for supersaturated vapor under different coolant pressures was carried out. Then it was applied to T sat to obtain viscosity of saturated vapor. It was presented in the form of μ g a μ g T sat + b μ g . Values of constants are shown in Table 1. Besides, dynamic viscosity of the liquid phase recommended in the study by Fox et al. (1998) was adopted in this work, which is shown as follows:

Surface Tension
The surface tension of water in this work was needed to determine the void fraction. We estimated the surface tension of water at coolant temperature and corresponding saturation pressure. The surface tension of water under these conditions was adopted as follows: The curve fitting result is shown in Figure 6.

Specific Enthalpy
Specific enthalpy of liquid and saturated gas was needed to calculate void fraction and other relevant quantities. Specific enthalpy of liquid was adopted as h l a h l T 3 + b h l T 2 + c h l T + d h l . The curve fitting was carried out between 473.15 K and T sat in order to make results 1.02 ± 0.21 1.02 ± 0.21  more precise. In addition, specific enthalpy of saturated vapor is only a function of pressure, which is h g a hg P 2 + b hg P + c hg . It is valid for the pressure between 15.5 and 10 MPa. Values of constants and figure representation are shown in Table 1 and Figure 7, respectively.

Saturation Temperature
Saturation temperature was also obtained by fitting the following equation: T sat 5.587 × ln(P P 0 ) 2 + 19.553 × ln(P P 0 ) + 377.448, (5) where P 0 1 bar. The curve fitting result is shown in Figure 8.

Heat Convection With Coolant
The single-channel model was used mathematically in this work to calculate the boundary condition. Several heat transfer coefficients on the outer surface of cladding were adopted depending on the heat transfer regime. We consider only the heat transfer conditions prior to the point of critical heat flux (CHF), and the corresponding coolant heat transfer coefficient is described as follows.

Single-Phase Forced Convection
Dittus-Belter correlation was adopted to calculate the heat transfer coefficient under the single-phase forced flow   condition (Dittus and Belter, 1985), as shown in Eq. 6. This equation is applicable for 0.7 < Pr < 100 and Re > 10000. Properties of coolant are evaluated at the film temperature T f (T W + T b )/2.

Subcooled Boiling Regime
Several studies were carried out to estimate the heat transfer coefficient under the subcooled boiling regime depending on coolant pressure and heat flux at the boundary. The most widely used correlations were proposed by Jens-Lottes and Thom et al. (Jens and Lottes, 1951;Thom et al., 1965). The expressions are described as follows: Jens-Lottes correlation: Thom correlation: Both Jens-Lottes correlation and Thom correlation cover typical PWR-type conditions. However, since Thom correlation is capable of calculating the coolant heat transfer coefficient in both subcooled boiling regime and saturated boiling regime (Aounallah, 2010), it could serve as a transitional correlation between two regimes. Therefore, it was adopted in this work.

Saturated Boiling Regime
Both Schrock-Grossman correlation and Chen correlation are suitable for the simulation in this region. However, since Schrock-Grossman correlation has a more concise expression and is well performed at low and intermediate pressure ranges (Aounallah, 2010), it was adopted in this work to calculate the coolant heat transfer coefficient in this regime (Schrock and Grossman, 1962), which is given by In this correlation, h coef , DB is the heat transfer coefficient in the liquid phase under the same heat flux, namely, the heat transfer coefficient calculated using Dittus-Belter correlation. The values of constants a 1 , a 2 , and b are a 1 7400, a 2 1.11, and b 0.66.

Regime Boundaries in Subcooled Boiling
The onset of nucleate boiling (ONB) is defined as the position where the first bubble of vapor appears on the heating wall. According to the study by Al-Yahia and Jo (2017), the evolution of heat flux and wall temperature around the ONB both present a turning point. T ONB was estimated as the intersection point of coolant temperature between single-phase-forced convection regime and subcooled-boiling regime. Furthermore, an energy balance provided the axial location of the ONB on the cladding wall: The onset of significant void (OSV) is defined as the point where the void fraction begins to increase significantly. The specific enthalpy of the coolant at the OSV was evaluated applying the Saha and Zuber model (Saha and Zuber, 1974): where Pe GD h c p,l /k l . Then we deduced T OSV from h OSV . z OSV could also be obtained using energy balance: The relation between G and G m is The flow area of the core and a single channel were estimated as follows: Typical parameters of assemblages in the PWR are shown in Table 2.

Void Fraction
Plenty of studies have been carried out to calculate the void fraction in the reactor core. In this work, we adopted the model proposed by Zuber and Findlay (Zuber and Findlay, 1965) and Manon (Manon, 2000). The void fraction was calculated as follows: FIGURE 9 | (A) Experimental data of wall temperature in subcooled boiling regime from the study by Cubizolles et al. (2009); the experimental data are compared against Jens-Lottes correlation. This test is subsequent to the boiling test shown in (C). (B) comparison of wall temperature predicted using Jens-Lottes correlation and Thom correlation under the normal PWR-type condition in this work. The temperature is compared at 10 7 s. (C) Experimental data of wall temperature in saturated boiling regime in the study by Cubizolles et al. (2009). (D) Simulation results of wall temperature under several PWR-type conditions in this work. In our model, the area of the entire core was estimated as 4.7 m 2 . Therefore, 17,000, 15,000, 12,000, 10,000, and 8,000 kg/s correspond, respectively, to 3,617, 3,191, 2,553, 2,128, and 1,702 kg/m 2 /s. Frontiers in Energy Research | www.frontiersin.org June 2021 | Volume 9 | Article 691813 where C 0 was given by Dix (Dix, 1971) β and b were calculated as The formula of drift velocity V gj was given by The actual steam quality proposed by Levy was calculated based on equilibrium steam quality: where x e is the equilibrium steam quality defined as However, this expression assumes that the actual steam quality is equal to zero at the OSV, which does not conform to experimental results. In order to obtain a nonzero value at the OSV, Manon (Manon, 2000) proposed the following expression for actual steam quality based on the Levy model (Chen, 1967):  where x(z OSV ) is calculated by where α(z OSV ) is the void fraction at the OSV. The expression proposed by Manon was adopted in this work. Since the actual steam quality and the void fraction increase slightly from the ONB to OSV, they were calculated linearly between the ONB and OSV, assuming that they were both equal to zero at the ONB. In order to determine α(z OSV ), we suggested to use the Griffith et al. (1958) model. This model estimates the void fraction at the OSV as follows: a q ″ k l Pr l 1.07 h 2 coef ,DB (T sat − T OSV ) .
The single-phase heat transfer coefficient h coef ,DB was calculated using the Dittus-Boelter correlation. All liquid properties were evaluated at OSV temperature. The empirical constant of 1.07 was used to adjust the curve to fit into experimental data.

Model Validation
In order to verify and validate the capability and reliability of the CAMPUS code, some comparisons of simulation results are needed. Validity of the CAMPUS code in predicting fuel rod performance has been proved in the study by Liu et al. (2015); Liu et al. (2016a); Liu et al. (2016b) by comparing the simulation results of CAMPUS with that of BISON, ABAQUS, and FRAPCON. In the study by Cubizolles et al. (2009), a series of fuel rod bundle heat transfer tests in the OMEGA-2 test facility were carried out. The experiments were carried out for a 5 × 5 fuel rod bundles, and they covered typical PWR-type conditions: coolant pressure ranges from 100 to 155 bar, mass flux ranges from 3,000 kg/m 2 /s to 4,600 kg/m 2 /s, and heat flux varies from 570 kW/m 2 to 1,400 kW/m 2 . The tests also covered singlephase-forced convection, subcooled boiling, and saturated boiling. Cubizolles et al. (2009) compared the experimental single-phase heat transfer data against Dittus-Boelter correlation. The comparison shows that under the typical PWR-type condition, all experimental data of heat transfer coefficient were distributed within 10 % of the upper and lower values predicted by the Dittus-Boelter correlation. In the subcooled boiling regime, according to the study by Cubizolles et al. (2009), wall temperature predicted using Jens-Lottes correlation was a little overestimated. However, the agreement was still within 1 K (see Figure 9A). The comparison of wall temperature predicted in this work using Jens-Lottes correlation and Thom correlation is shown in Figure 9B. We can see that under the normal PWR-type condition, the difference is within 2 K, which leads to the validity of Thom correlation. Experimental data of boiling tests in the study by Cubizolles et al. (2009) are shown in Figure 9C. Under the condition of 15.5 MPa and 3,440 kg/m 2 / s, the wall temperature remains almost constant with an extremely slight decrease. Figure 9D shows one of the simulation results of wall temperature in this work. Since the coolant heat transfer coefficient predicted by saturated boiling correlation predominates and increases only at high-steam quality, the wall temperature in the height range studied is maintained almost as a constant. Therefore, it is believed that simulation results in the saturated boiling regime are verified.

Void Fraction
In this section, the stationary modeling results of void fraction for UO 2 fuel under a drop of coolant pressure and mass flux are presented to provide a reference for the LOCA. Steam quality and void fraction are averaged between 10 4 and 10 7 s because calculated results show that their variation in this period is negligible. Averaged void fraction as a function of x e under different coolant pressures and mass fluxes is shown in Figure 10. The figure presents only the calculated void fraction before reaching the saturated boiling regime because the distinction mainly occurs in the subcooled boiling region. From the figure, we observe that as coolant pressure decreases, x e at both the ONB and OSV increases. However, the variation of x e at the OSV is much smoother than that of x e at the ONB, which leads to a shrink of the length of the low subcooled boiling region. The calculated results also show that the void fraction presents less dependency on mass flux with the increase of mass flux. This phenomenon is observed over all coolant pressures investigated.
The simulation results of void fraction are shown in Figure 11, as a function of axial location on the fuel rod. We took inlet coolant temperature as 566 K. Under 15.5 MPa, with mass flux investigated, the flow pattern of coolant passes from singlephase-forced convection to two-phase flow. As the pressure and mass flux decrease, the region of single-phase flow is increasingly replaced by two-phase flow. When the pressure decreases to 10 MPa, almost the entire fuel rod cladding is covered by the coolant in subcooled boiling and saturated boiling. This might provide a reference for the prediction of the starting point of critical heat flux under a LOCA.

Material Performance
In this section, UO 2 fuel performances such as hoop stress, axial stress, and radial stress under a drop of coolant pressure and mass flux were investigated. In our model, the simulation time lasts for 10 7 s, eventually reaching the burnup of nearly 120 WM h/kg(U). Figure 12 presents the evolution of hoop stress with fuel burnup for two operating conditions: one at 15.5 MPa and 17000 kg/s, and the other at 15.5 MPa and 12000 kg/s. Hoop stress was calculated as follows: Calculated results are compared at some similar positions in both cases. First, we can see that once the fuel power reaches the expected value, hoop stress in both cases is maintained almost as a constant within 120 WM h/kg(U), and relative variations of hoop stress remain within 5 %. Therefore, it is reasonable to study the average value of hoop stress in this range of burnup. Besides, it seems that mass flux has little influence on the evolution of hoop stress. Dependency on height is not apparent either.
Averaged hoop stress under different coolant pressures and mass fluxes is shown in Figure 13. The average value of hoop stress between 10 4 and 10 7 s is investigated since wall temperature is believed to reach the stationary regime. We can see that under normal coolant pressure (15.5 MPa), hoop stress increases almost linearly with height at the beginning part of the fuel rod. Nevertheless, after a certain point, it reaches saturation and then keeps fluctuating around the saturation line, whatever the mass flux is. This is similar to the distribution of wall temperature along axial location, as shown in Figure 9D. Under normal coolant pressure, within the range of height investigated, wall temperature increases linearly with height in the singlephase-forced convection region and reaches saturation eventually. Therefore, it could be inferred that the distribution of averaged hoop stress indicates exactly the distribution of wall temperature. The point where the hoop stress takes a turn is the starting point of subcooled boiling. We can also observe that hoop stress presents strong dependence on coolant pressure. The   Figure 14 presents the evolution of axial stress with fuel burnup under different coolant pressures: 15.5, 13, and 10 MPa. Under each pressure, two operating conditions are compared: one is 17000 kg/s and the other is 8,000 kg/s. The figure shows simulation results at some similar axial locations of the fuel rod. It can be remarked that axial stress has a relatively large variation within the fuel burnup investigated, especially in the cases of high pressure, different from the evolution of hoop stress. Under the condition of 15.5 MPa and 17000 kg/s, at 4.85 m, axial stress varies from −190 MPa to −130 MPa eventually. The evolution of axial stress with burnup may depend on axial location. At a low axial location of the fuel rod, axial stress decreases more smoothly with burnup, different from the case at a high axial location, as shown in Figure 14A. At higher location such as 3.80 and 4.85 m, the evolution curves even coincide with each other. Combining with Figure 9D, we conclude that the entering of subcooled boiling regime in these two locations leads to the similarity of axial stress evolution. As shown in Figures 14B, D-F 2), 4), 5), and 6), axial stress presents similar trend of evolution at different axial locations since single-phase convection no longer occurs at these positions. Once the coolant flow enters the two-phase region, the evolution of axial stress with burnup at different axial locations follows the same pattern, also independent of mass flux.
Comparing the modeling results in Figures 14B, D, F 2), 4), and 6), we can see that the drop of coolant pressure seems to retard the decrease rate of axial stress with burnup in the twophase flow region. In Figure 15, axial stress at 117 WM h/kg(U) is compared. We can see that axial stress distribution almost coincides with wall temperature distribution, similar to hoop stress. However, at a high mass flux, axial stress distribution in single-phase regions is not as linear as wall temperature.
The evolution of radial stress with fuel burnup is shown in Figure 16, and the simulated results under three coolant pressures are presented: 15.5, 13, and 10 MPa. Under each pressure, modeling results at a similar height with two different mass fluxes are compared: 17,000 and 8,000 kg/s. We can see that radial stress does not show a substantial variation with burnup, similar to hoop stress. Compared to hoop stress and axial stress, the average value of radial stress is less remarkable. Therefore, it is less influential to the cladding material.

CONCLUSION
In summary, based on a previous work by Liu et al. (2016a), we have built up a self-defined multiple physics model based on COMSOL to simulate fuel performance of UO 2 under normal conditions and a drop of both coolant pressure and mass flux. By developing the model and analyzing the modeling results, a more comprehensive understanding of both flow regimes before CHF and fuel performance under a loss of coolant pressure and mass flux is presented in this work, providing a reference for the LOCA.
Prediction of void fraction shows that x e at the ONB and OSV increases with the diminishing of coolant pressure at the same mass flux. However, x e at the ONB changes more rapidly so that the length of the low subcooled boiling region shrinks. We also observe that void fraction becomes less dependent on mass flux as mass flux increases.
The study of material performances, hoop stress, radial stress, and axial stress is investigated. First of all, hoop stress and axial stress play obviously more important roles in the safety of the cladding material. In addition, the evolution of stress on Zr cladding may also serve as an indicator of flow pattern, since both hoop stress and axial stress have the same distribution with wall temperature along the axial direction of the fuel rod, determined by the flow pattern. The turning point signifies the start of the subcooled boiling regime. At last, under the condition of the LOCA, on the one hand, the hoop stress on Zr cladding is influenced the most by the drop of coolant pressure. Both hoop stress and axial stress are limited by a value dependent on the saturated boiling regime. On the other hand, the drop in coolant mass flux makes nearly no difference on radial stress, different from hoop stress and axial stress.
In our future work, more research including thermal fluid coupling after CHF, transient change of coolant pressure and mass flux, and other fuel materials with high thermal conductivity like UO 2 -BeO will be carried out to simulate the LOCA more realistically.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
The work is mainly done by MY. LQ contributed to the accomplishment of this work, and WZ is the advisor of MY.