Verification of a self-developed CFD-based multi-physics coupled code MPC-LBE for LBE-cooled reactor

Zhi-Xing Gu· Qing-Xian Zhang · Yi Gu · Liang-Quan Ge · Guo-Qiang Zeng ·Mu-Hao Zhang · Bao-Jie Nie

Abstract To perform an integral simulation of a pool-type reactor using CFD code, a multi-physics coupled code MPC-LBE for an LBE-cooled reactor was proposed by integrating a point kinetics model and a fuel pin heat transfer model into self-developed CFD code. For code verification, a code-to-code comparison was employed to validate the CFD code.Furthermore,a typical BT transient benchmark on the LBE-cooled XADS reactor was selected for verification in terms of the integral or system performance. Based on the verification results, it was demonstrated that the MPC-LBE coupled code can perform thermal-hydraulics or safety analyses for analysis for processes involved in LBE-cooled pool-type reactors.

Keywords LBE-cooled pool-type reactor · Computational fluid dynamics · Multi-physics coupling code · Safety analysis code · Verification

1 Introduction

Lead/lead–bismuth eutectic (LBE)-cooled reactors,abbreviated as lead-based reactors,have been announced to be promising advanced nuclear energy systems owing to their superior characteristics in nuclear safety, sustainability, and waste disposal capability [1]. Many countries and organizations around the world are very active in research and development concerning lead-based reactors[2–4].As a type of generation IV advanced nuclear energy system, lead-based reactors tend to employ integral pooltype configurations [2–4]. In the literature, it has been noted that the above configurations may cause a series of complicated multi-dimensional thermal-hydraulics and safety problems concerning liquid–metal-cooled reactors,such as thermal stratifications in the upper/lower plenum,thermal striping at the nuclear core outlet, and cover gas entrainment [5]. To cope with these special problems in pool-type reactors, several system-level safety analysis codes have been developed or secondarily developed,such as SAS4A/SASSYS-1[6],RELAP5[7],and ATHLET[8].However, considering the complicated pool-type thermal–hydraulic and safety characteristics, system analysis codes that always employ 0D lumped parameter models for pool behaviors are relatively insufficiently elaborate in terms of pool-type reactor simulations [5].

Computational fluid dynamics (CFD)-based simulation is considered as an effective method to simulate the multidimensional and complicated behaviors involved in a pooltype reactor because the partial differential equations describing the detailed flow and heat transfer phenomena are solved completely, instead of using lumped parameter models. Owing to their powerful capabilities in geometric modeling and mesh generation, commercial CFD tools have been widely applied in 3D simulations of pool-type reactors. For example, in order to study the thermal stratification and natural circulation characteristics of a pooltype reactor, Fluent was employed by Sakamoto et al. to perform a 3D simulation of the loss of flow transient in the sodium-cooled fast reactor Monju[9].3D mixing and heat transfer process in the upper plenum of MONJU reactor during a turbine trip test was simulated by using CFD techniques [10]. Further simulations of thermal stratification problem in MONJU were carried out using FLUENT[11, 12] and CFX [13] tools. And CFD-based multi-dimensional thermal-hydraulics simulations were also performed in other innovative nuclear energy systems,such as TMSR-SF [14] and pebble-bed reactor [15].

Although CFD tools are very powerful, it is somewhat unrealistic to perform the entire reactor simulation using CFD tools alone owing to the vast computational resources required and the absence of other physical models,such as neutron dynamics models and core heat transfer models.Alternatively, multi-scale coupling simulations have been proven to be a good choice. Some corresponding activities have been carried out by coupling the above CFD tools with system analysis codes. Specifically, CFD tools are used for pool behavior simulation, while system analysis codes are employed to simulate the core and other systems.For example, to perform an elaborate simulation of a loss of flow accident in a sodium-cooled fast reactor, multiscale coupling between SAS4A/SASSYS-1 and STAR-CD was conducted by Fanning T. H. and Thomas J. W. [16].Another two typical multi-scale coupling simulation activities were performed by Pialla D. et al. (system code CATHARE and CFD code TRIO_U; system code ATHLET and CFD code OpenFoam) [17] and by Bandini G.et al. (system code RELAP5/MOD3.3 and CFD code STAR-CCM +) [18]. Multi-scale coupling simulation helps to realize elaborate and sufficient simulations of transient processes concerning pool-type reactors. However, this approach still requires massive computing resources, and numerical errors will inevitably occur during the data mapping process between different dimensionalities and scales. Furthermore, just as a premise, we must have reasonable access to the system analysis program, especially its source code.

In order to avoid vast and multi-dimensional data mapping in multi-scale coupling simulation, as another option for CFD application in pool-type reactors,integrated simulations of reactors can be implemented based on CFD alone. As CFD tools do not have additional models for reactor transients, most importantly the fuel pin heat transfer model and neutron dynamics model,it is necessary to conduct multi-physics coupling activities. For example,1D fuel pin heat transfer (HT) and 0D point kinetics (PK)models were integrated into the ANSYS Fluent platform for multi-physics coupling with the use of UDF tools by Chen Z.et al.[19].An advanced 2D HT model and 0D PK model were developed and integrated into ANSYS Fluent by Gu Z. in our previous work [20]. The advantage of CFD-based multi-physics coupling is that powerful preprocessing tools, such as geometric modeling and mesh generation,can be implemented.However,CFD tools used in multi-physics coupling are always commercialized programs whose source code is not publicly available, so it is almost impossible to change their key functionality for advanced coupling algorithms and modeling of key equipment such as pumps and heat exchangers.To address this problem, self-developed CFD-based multi-physics coupling seems to be a good choice.A typical activity was performed during SIMMER-III development, which is based on a 2D axisymmetric cylindrical coordinate system used for transient safety analysis, especially for severe accident analysis of liquid–metal-cooled reactors [21–23].As SIMMER-III is a multi-component,multi-phase,multivelocity CFD model coupled with a neutron transport model and pin fuel heat transfer model, it is very timeconsuming.

In this study, a single-phase CFD code based on a 2D axisymmetric cylindrical coordinate system was developed, and a multi-physics coupled code MPC-LBE was developed for the transient safety analysis of an LBEcooled reactor by integrating a similar 2D axisymmetric fuel pin heat transfer model and an advanced 0D point kinetics model to the CFD code. Verifications were conducted by a code-to-code validation for the CFD code alone and a transient benchmark for the integral performance of the MPC-LBE.

2 Models of MPC-LBE code

2.1 Computational fluid dynamics

2.1.1 Conservation equations and discretization schemes

The LBE coolant is treated as an incompressible fluid whose density is assumed to be determined only by its temperature. As shown below, the physical properties of LBE, namely viscosity, conductivity, and specific heat capacity,are also considered to vary only with temperature[20].

To solve the above conservation equations,first,they are discretized using finite volume methods considering staggered mesh strategies, as shown in Fig. (1). The mass and energy conservation equations are integrated in the same control volume whose centric coordinates are （i，j ）, where variables such as density, temperature, and pressure are assumed to be located. However, the radial and axial momentum conservation equations are integrated in the staggered control volume whose centric coordinates are（i +1／2，j） and （i ，j+1／2）, respectively.

It is assumed that the physical mechanisms of the flow and heat transfer processes are quasi-static,so the mass and momentum conservation equations are considered to be solved separately from the energy conservation equation.The energy conservation equation is solved first. Thus, in this hypothesis, the mass conservation equation is discretized, as shown in Eq. (9). Here, the first-order upwind scheme is adopted, as shown in Eqs. (10) and (11). In the core regions, it is obvious that cell volumes are occupied by both fuel pins and fluid(coolant,namely LBE),the fuel pins are arranged vertically with certain gaps between them, and the areas at the north and south sides of the control volume for the mass conservation equation are modified by using coefficients （αn）i，jand （αs）i，j,respectively.

Fig. 1 Staggered mesh strategies

The energy conservation equation is discretized as shown in Eq. (12), in which we use the coefficient αi，jto account for the existence of the structure in the control volume. In addition, the convective term is treated using the first-order upwind scheme,as demonstrated in Eq. (13).The heat diffusion term is treated using the equilibrium heat transfer principle, as illustrated in Eq. (14).

Finally,we can obtain the algebraic form of Eq. (12),as follows:

The momentum conservation equation in the radial direction is integrated in the control volume centered on（i +1／2，j）, and is discretized as the formation shown in Eq. (16). A semi-implicit scheme is adopted for the discretization of convective terms, where the first-order upwind scheme is also considered, as in Eq. (17). In the core regions, it is obvious that the cell volumes are occupied by both fuel pins and fluid, the former of which are arranged vertically with gaps between them; thus, as seen in Eq. (16), the corresponding control volume for momentum is corrected by the coefficient αi+1／2，j, and the areas at the north and south boundaries of this control volume are also corrected by the coefficients （αn）i+1／2，jand（αs）i+1／2，j, respectively. Using similar methods, we can obtain the discretized formula for the axial momentum conservation equation, as shown in Eqs. (18) and (19).

Finally, we can obtain the algebraic form of the radial and axial momentum conservation equations, as given in Eq. (16) and (18), respectively.

When the Emperor heard from the Prince how he had gained possession of his fair prize, he at once recognized that he had been helped by some magic art, and on the spot gave up all claim to the beautiful mermaid

2.1.2 Boundary conditions

In nuclear reactors, a virtual wall is always used to control the flow and heat transfer.Here,it is defined that if a certain type of virtual wall (north, south, west, or east)exists at an inner boundary of a control volume cell,it will cut off the heat and mass transfer across it completely.For the flow condition in the tangential direction, a no-slip model was used for the virtual wall.To simulate the pumps in the primary vessel of the reactor, a simple pump model was used by providing additional pressure source terms to the cell pressure in certain directions.To simulate the heat exchanger in a nuclear reactor, an ideal heat exchanger model is adopted by assuming the density in certain cells to be infinite.It is obvious that the reactor core is occupied by fluid and structures; therefore, the fraction of fluid in such control volumes should be considered.

2.1.3 Algorithms and solvers

Just as shown in Fig. 2, the typical SIMPLE algorithm was used to solve the above equations. First, the flow and heat transfer processes were treated separately and successively. Thus, by using the flow information in the last time step, such as un，vn, we can obtain the temperature Tn+1for the next time step by iteratively solving Eq. (15).Then, the density distribution ρn+1for the fluid can be updated using formula (6). After solving the energy conservation equation, we begin to solve the momentum equations. First, with the improper pressure field pn,velocities in the two directions are estimated using typical semi-implicit schemes, namely Eq. (20). Then, we obtain u*，v*, which are used to check the mass conservation condition as described below, where ε*ρdenotes the continuity equation residuals. As we know,the current pressure pnis usually not consistent with the mass conservation equation, and the pressure correction equation must be solved to update the cell pressures.

Fig. 2 Algorithms using to solve N-S equations

Substituting them into the algebraic Eq. (9), we obtain Eq. (23).

Just as shown in Eq. (22), the correction values of velocities u′，v′can be approximated using algebraic Eq. (20), where the influence of adjacent cells is not considered. By substituting Eq. (24) into Eq. (23), the final pressure correction equations can be obtained as shown in Eq. (25), which can be transformed into Eq. (26), whose coefficients can be obtained using Eq. (27).

2.2 Fuel pin heat conduction

The core idea of the fuel pin heat conduction model is the same as that of the HC-PK-CFD code,which is a multiphysics coupled code based on the commercial CFD software ANSYS FLUENT platform, created by using the platform’s user-defined function (UFD) tools developed in our previous work [20]. The heat transfer processes in the fuel pellet and cladding are heat conduction. However,only heat conduction by helium with high pressure, whose thermal properties have been published by D’Angelo A.et al.[24],is considered in the fission gas plenum(namely,the fuel-clad gap). As shown in Eq. (28), the nonlinear partial differential equation in an axisymmetric 2D cylindrical coordinate system is employed to describe the heat conduction phenomenon [20].

Another sub-mesh system for the fuel pin heat transfer should be established based on the CFD mesh system.Equation (28) is discretized, as indicated by Eq. (29), in the above sub-mesh system using the finite volume method with a complete implicit scheme. The equilibrium heat transfer principle is employed to calculate the temperatures at the internal cell interfaces,as shown in Eq. (29)[20].As for the boundary conditions, the cladding outer surface is assumed to be thermally coupled with the LBE coolant by the convective heat transfer mechanism described in Eq. (30), while the other boundaries of the fuel pin are considered to be adiabatic [20, 24]. As shown in Eq. (30),hjrepresents the convective heat transfer coefficient at a specific axial position in a certain channel that consists of one sub-assembly or more; TnLBE，jindicates the corresponding temperature of the LBE coolant along the axial direction.

2.3 Point reactor kinetics

To perform a transient process simulation, it is necessary to consider the neutron dynamics model. The typical point kinetics model described by Eqs. (31) and (32) was employed to simulate the neutron dynamics behavior. To solve the above equations,a new semi-analytical algorithm with high computational accuracy and efficiency based on the Euler-Maclaurin approximation (SAEMA) was developed in our previous work [25]. The MPC-LBE code was based only on the SAEMA algorithm. However, the corresponding point reactor kinetic module of the HC-PKCFD code is based on another semi-analytical algorithm that uses a one-order Taylor polynomial approximation[20]. For the SAEMA algorithms, first, the value of the neutron density of the next time-step, n（t+h）, can be calculated using Eq. (33), while its derivative with respect to time,n′（t+h）,can be calculated using Eq. (34).Finally,the density of the delayed neutron precursor,Ci（t+h）,can be calculated using n（t+h） and n′（t+h）, as presented in Eq. (35).Detailed algorithms can be found in our previous work [25].

2.4 Multi-physics coupling

At this stage of our work, by employing the same coupling algorithm as the HC-PK-CFD coupled program in our previous work, the above fuel pin heat transfer model and point reactor kinetics model are coupled in the self-developed CFD code [20]. A typical external loose (quasi-static hypothesis)coupling scheme is used in this study,as depicted in Fig. 3.Three significant coupling processes are depicted below. First, with respect to the coupling process between the CFD model and the pin fuel heat transfer model, the former provides the latter with the temperature cladding outer surface, and the latter provides the former with heat flux across the cladding outer surface. Second,the CFD model provides the point reactor kinetics model with coolant temperature and density for feedback calculation, and the point reactor kinetics model provides the CFD code with a heat source. Third, as for the coupling between the pin fuel heat transfer and point reactor kinetics models, the former passes the fuel pellet temperature and density to the latter, while the latter passes the heat source to the former. As for the time-step control strategies, the three modules use the same time-step size (0.01 s) to perform the coupling calculations.

Fig. 3 Coupling scheme of multi-physics module

Fig. 4 (Color online) Comparison of convergence with mesh sizes between MPC_LBE and Fluent

Fig. 5 (Color online) Comparison of velocity and temperature evolutions averagely between MPC_LBE and Fluent

Fig. 6 (Color online) Local comparison of velocity magnitude and temperature evolutions between MPC_LBE and Fluent

3 Verifications on MPC-LBE

3.1 Computational fluid dynamics model

In order to verify the CFD module, a transient natural circulation flow and heat transfer in a cylindrical vessel(radius × height: 0.1 m × 0.1 m) driven by a uniform volumetric heat source equal to 106W/m3(radial:0–0.05 m, axial: 0–0.05 m) are chosen. This was modeled using a 2D cylindrical coordinate system.The mesh for the MPC-LBE code is generated uniformly using an input card in the form of a structured type. The commercial CFD platform Fluent was used to verify the MPC-LBE code. A 2D axisymmetric condition was also used in Fluent. Its computational mesh was also generated uniformly using the ICEM tool with structured mesh partitioning technology. Similar to MPC-LBE, the first upwind scheme was adopted for convective term discretization, the SIMPLE algorithm was used for pressure iteration, and the laminar flow model was considered. To establish the accuracy of MPC-LBE against Fluent, mesh convergence analysis was first carried out,as shown in Fig. 4,where five incremental mesh strategies are investigated. It is obvious that the MPC_LBE code has almost the same trend approaching convergence as the mesh size increases compared with the commercial code of Fluent. As this section focuses on code-to-code verification, further investigation on mesh convergence was not carried out.Second,a comprehensive comparison between the two codes was investigated by choosing an 80 × 80 mesh size. Figure 5 shows a comparison of the evolution of the overall average velocity magnitude and temperature calculated by the two codes. It can be seen from this figure that very good agreement is achieved for the MPC_LBE code benchmarked by Fluent.In addition, detailed local variations of velocity magnitude and temperature computed by the MPC_LBE code are also compared with Fluent, as shown in Fig. 6. An arbitrary position with coordinates of (0.050625, 0.050625) was selected for evaluation. It can be seen from this figure that although there are several slight deviations (especially in velocity around 28 s) among the two codes, overall, the two codes also agree reasonably well with each other.

Fig. 7 Comparison of velocity distribution of fluid at 122 s between MPC_LBE and Fluent

Fig. 8 (Color online)Comparison of temperature distribution of fluid at 122 s between MPC_LBE and Fluent

Fig. 9 (Color online) Computational model for MPC-LBE code

Finally,comparisons in terms of the overall distribution of velocity and temperature were investigated,as shown in Figs. 7 and 8. For the two figures, the results at 122 s computed by the two codes are illustrated with the same coordinate scale and variable magnitude. Obviously, it is difficult to find obvious differences between the two codes.As demonstrated in Fig. 8, even better agreement is achieved for comparisons of temperature distribution between the two codes.In this section,it can be concluded that the CFD module of the MPC_LBE code can effectively simulate the flow and heat transfer phenomenon.

3.2 Integral coupling level

3.2.1 Geometric model

To verify the capability and correctness of the MPCLBE code at an integral coupling level, a typical international transient benchmark on the beam trip (BT) accident for the LBE-cooled XADS reactor was employed[24].The geometric model of this benchmark is described as a single averaged fuel pin channel of the XADS reactor [24]. We adopted it for the HC-PK-CFD coupled program in our previous work,in which it was assumed to be only a single flow channel with inlet and outlet boundaries [20]. However, in order to validate the MPC-LBE code based on the entire reactor primary vessel level, a new enclosed geometric model was established, as shown in Fig. 9. The radius and height of this reactor primary vessel are 1.96 m and 5.72 m, respectively. As shown in Fig. 9, this primary vessel is mainly composed of a reactor core,hot pool,cold pool, lower plenum, pump, pump pipes, and a heat exchanger.Specifically,for the reactor core,neutron source sub-assemblies are located at channels 1–3, fuel subassemblies are located at channels 4–7, 9–11, and 13–15,reflector sub-assemblies are located at channels 16–21,shielding sub-assemblies are located at channels 22–24,and the control rod sub-assemblies are located at channels 8 and 12.It is assumed that all the sub-assemblies have the same axial profiles for volumetric power density values,but the radial profiles for these values vary, and the normalized amplitude factors for these values in each channel from No. 1 to No. 24 are 0.01, 0.008, 0.008, 1.0, 1.0, 0.9,0.9, 0.001, 0.78, 0.78, 0.78, 0.001, 0.52, 0.52, 0.52, 0.001,0.001, 0.001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, and 0.0001, respectively. It should be noted that the object channel is located at channel 4,including a certain number of single fuel pin channels,similar to the geometric model in the benchmark report mentioned above [24]. To ensure the same flow boundary as that of the benchmark model,the resistance coefficient (modeled by setting the orifice coefficients, which are Cr，orifice/Cz，orificein formula (5)), in the objective channel must be adjusted constantly to keep the inlet velocity unchanged.

3.2.2 Material thermal properties and heat transfer coefficient

To maintain consistency with the benchmark, the heat source distribution of the pin fuel in this new model was established based on the data published in the benchmark report mentioned above [20].The thermal property data of the LBE coolant,fuel pellet,fission gas,and cladding were also derived from the benchmark report [20, 24]. The relationship of the convective heat transfer coefficient between the fuel pin and LBE coolant was established based on that in the benchmark report [20, 24], and is demonstrated in Eq. (36), where h（z ）, Nu（z ）, Pe（z ）,ρCoolant（z ）, kCoolant（z ）, vCoolantz（ ）, and DCellare the convective heat transfer coefficient, Nusselt number, Peclet number,coolant density,coolant heat conductivity,coolant velocity, and hydraulic diameter of each CFD cell along the axial direction.

Fig. 10 (Color online) Material axial temperature distribution (steady state)*. *Results given by MPC-LBE are shown in the corresponding figure of the benchmark report [24] using the Origin tool

Fig. 11 (Color online) Steady state simulation results:a velocity distribution;b temperature distribution

Fig. 12 (Color online) Fuel pin temperature distribution (Steady state)

3.2.3 Neutron dynamic parameters

The typical point reactor kinetics model was adopted in the MPC-LBE code to simulate the neutron dynamics behaviors.The shape of the neutron density profile remains unchanged throughout. The effective delayed neutron fraction was set to 350 pcm based on the benchmark report mentioned above [20, 24]. For the initial state, the XADS reactor operates at a sub-criticality of approximately 8$ ($indicates the measurement standard of reactivity which equals to the total fraction of delayed neutrons in reactor)[24]. The detailed operating data, such as the neutron generation time and delayed neutron parameters, were derived from that benchmark report [24]. The reactivity feedback model for the fuel pellet and coolant was established as shown in Eqs. (37) and (38) based on that of the benchmark report [24].

Fig. 13 (Color online)Objective channel normalized power evolution in different cases of BT*.*Results given by MPC-LBE are shown in the corresponding figure of the benchmark report[24]using the Origin tool

Fig. 14 (Color online) Fig. 13 enlargement*. *Results given by MPC-LBE are shown in the corresponding figure of the benchmark report[24]using the Origin tool

Fig. 15 (Color online)Objective channel outlet coolant temperature evolution in different cases of BT*.*Results given by MPC-LBE are added to the figure of the benchmark report[24]using the Origin tool

3.2.4 Computation results and code verification

First, the steady state was calculated using the MPCLBE code,and the material temperature axial distributions in the objective channel of the new geometric model were compared with those calculated using ten codes from the benchmark report[24].In Fig. 10,as in our previous work,the steady-state temperature distributions are added to the figure in the benchmark report using the Origin tool[20, 24]. From this figure, it is demonstrated that satisfactory agreement is achieved between the results calculated by the MPC-LBE code and the codes in the benchmark report [24]. As an integral simulation of the primary coolant system was conducted using our self-developed CFD code, MPC-LBE, the corresponding temperature and velocity distributions are shown in Fig. 11. It can be seen from this figure that distinct thermal stratification occurred in the hot pool region, and the flow distribution was realized by the MPC-LBE code. The fuel pin temperature distribution in 2D axisymmetric cylindrical coordinates is shown in Fig. 12.

Furthermore, as in our previous work [24], a typical accident, namely beam trip, in the XADS reactor wassimulated using MPC-LBE based on the above steady state.Four cases were considered in terms of different beam recovery times: 1 s, 3 s, 6 s, and 12 s. In Figs. 13 and 14,the evolutionary process of the normalized power in the objective channel is shown. As depicted in these two figures, the normalized powers for different cases calculated by the MPC-LBE code agree well with those calculated by the other ten codes in the benchmark reports [24].

Fig. 17 (Color online)Objective channel mid-plane fuel internal temperature evolution in different cases of BT. *Results given by MPCLBE are added to the figure of the benchmark report[24]using the Origin tool

In the objective channel, the outlet coolant temperature evolutionary processes for different BT cases are shown in Fig. 15,and the mid-plane fuel pellet external and internal temperatures in different BT cases are depicted in Figs. 16 and 17, respectively. A detailed description of the evolutionary processes can be found in our previous work[24]. From these figures, it can be observed that the MPCLBE agrees well with the other ten codes in the benchmark report [24].

4 Conclusion

In order to avoid vast and multi-dimensional data mapping in multi-scale coupling simulation, as well as to provide another option for CFD application in pool-type reactors, it is important to implement an integrated reactor simulation based on CFD platforms. In this study, a selfdeveloped CFD-based multi-physics coupled code MPCLBE was developed for the CFD-based integral simulation of thermal-hydraulics and safety analysis of LBE-cooled pool-type reactors. Detailed models and algorithms for the MPC-LBE code are presented. The first-stage verification activity for the MPC-LBE code was conducted successfully.First,a code-to-code verification work was employed to verify the capability and correctness of the self-developed CFD code alone.Then,for verification in terms of the integral or system level, a transient benchmark on the typical BT accident with different cases in the LBE-cooled XADS reactor was selected. In this verification, a new geometric model was established to simulate the primary coolant system based on the single-channel geometric model in the benchmark.The results of the proposed MPCLBE code agree very well with those of the other ten codes in the above benchmark report. In conclusion, it has been demonstrated that the MPC-LBE coupled code is capable of performing sufficiently accurate simulations of thermalhydraulics and safety analysis problems involved in pooltype LBE-cooled reactors. In our future work, studies on implicit coupling strategies will be conducted.

Author contributionsAll authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Zhi-Xing Gu, Qing-Xian Zhang and Yi Gu. The first draft of the manuscript was written by Zhi-Xin Gu and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.