Progress of Particle Flow, Fluid/Solid Mechanics, and Heat Transfer in Advanced Gas/Water Nuclear ReactorsView this Special Issue
Research Article | Open Access
Shiming Wang, Cheng Ren, Yangfei Sun, Xingtuan Yang, Jiyuan Tu, "A Study on the Instantaneous Turbulent Flow Field in a 90-Degree Elbow Pipe with Circular Section", Science and Technology of Nuclear Installations, vol. 2016, Article ID 5265748, 8 pages, 2016. https://doi.org/10.1155/2016/5265748
A Study on the Instantaneous Turbulent Flow Field in a 90-Degree Elbow Pipe with Circular Section
Based on the special application of 90-degree elbow pipe in the HTR-PM, the large eddy simulation was selected to calculate the instantaneous flow field in the 90-degree elbow pipe combining with the experimental results. The characteristics of the instantaneous turbulent flow field under the influence of flow separation and secondary flow were studied by analyzing the instantaneous pressure information at specific monitoring points and the instantaneous velocity field on the cross section of the elbow. The pattern and the intensity of the Dean vortex and the small scale eddies change over time and induce the asymmetry of the flow field. The turbulent disturbance upstream and the flow separation near the intrados couple with the vortexes of various scales. Energy is transferred from large scale eddies to small scale eddies and dissipated by the viscous stress in the end.
With the increasing acute global energy crisis, more and more attention is paid to seek a green environment-friendly and low-carbon development pattern. Nuclear energy is identified as a competitive alternative in the energy markets. The uranium-based, graphite-moderated, helium-cooled pebble bed reactor high temperature reactor (HTR) has been proved to be one of the most promising nuclear power projects in China due to its high efficiency, attractive economics, high levels of passive safety, and potential usage for hydrogen production . Also, some test reactors, such as AVR in Germany , MPBR in USA , and PBMR in South Africa , have been demonstrated around the world. The Chinese high-temperature gas-cooled reactor-pebble bed module (HTR-PM) demonstration project, oriented by Institute of Nuclear and New Energy Technology of Tsinghua University (INET), is very attractive and competitive in the nuclear power market and is going well .
The HTR-PM has one steam turbine-generator and two independent primary loops, which have the same construction design and are arranged in the so-called “side-by-side” way . Both of the loops have their own reactor modules and steam generators, as shown in Figure 1(a). The helium in the loop is heated to 750°C by the high temperature fuel spheres in the core of the reactor and flows into the chamber of the spiral tube exchanger in the steam generator through the inner pipe of the coaxial hot gas duct. Then the hot helium releases its thermal power to the water in the secondary loop and is cooled down to around 250°C. The cooled helium from 14 riser pipes collects into the header box and is pumped back to the reactor by the blower. Each of the riser pipes is connected to the header box by a 90° elbow, as shown in Figure 1(b), which is used as the sensor to measure the helium flow rate by monitoring the pressure difference between its intrados and extrados [7–9].
The elbow sensor does not introduce any additional pressure loss in the pipe system other than that already caused by the elbow. As an enormous amount of heat must be extracted out of the core, enough helium is needed to meet the requirement. Hence, the accurate measurement of helium gas flow rate in the primary loop is closely related to the safe operation and power control of the reactor system. Considering that the elbows arranged behind the riser pipes is used as the sensor to measure the helium flow rate, the characteristics of the turbulent flow in a 90° elbow must be studied to guarantee the accuracy and reliability.
In a 90° elbow, the central parts of the flow which go more rapidly are forced outwards by the centrifugal actions, while the slower parts near the extrados are forced inwards to the intrados along the circumference wall by the pressure difference between extrados and intrados. Then the so-called secondary flow takes place in either a clockwise or an anticlockwise sense about the axis [10, 11]. These behaviors induce a squeeze effect to the extrados and a suction effect to the intrados. Along with the flow separation occurring on the intrados, the instantaneous velocity field and pressure field in the elbow turned out to be very complicated.
The development of secondary flow and flow separation is classified as transient process. The interaction between them and the turbulent disturbance from upstream are factors of significant influence on the flow field. There have been a number of representative studies on the flow field in elbow pipes due to the progresses in both computational fluid dynamics and experimental methods. Pigott , Kirchbach , and Schubart  measured the pressure drop of flowing through a 90° elbow and provide the detailed results including the methods used to calculate the elbow pipe pressure drop. Enayet et al.  and Sudo et al.  investigated the turbulent flow in a circular-sectioned 90° elbow by using laser Doppler velocimetry and single inclined hot wire, respectively. Eguchi et al.  applied the LES method to simulate the high-Re flow in a short-elbow pipe and the mechanism of the pressure fluctuation and the origin of the fluid force were discussed. And the same method was utilized by Rütten et al. [18, 19] to investigate the unsteady flow separation, unstable shear layers, and an oscillation of the Dean vortices in 90° elbow pipes.
However, there are few works on the asymmetry of the instantaneous pressure field and velocity field in 90° elbow pipes. Tunstall and Harvey  studied the fully developed turbulent flow through a mitred sharp bend experimentally and declared that the secondary flow was dominated by a single circulation about the axis in either a clockwise or an anticlockwise sense, between which it switches abruptly at a low, random frequency. Ono et al.  investigated the influence of elbow curvature on flow structure by using PIV method and concluded that the separation region was generated constantly in short elbows but was formed intermittently in long elbows. All of the characteristics mentioned above have significant influence on the flow field in the elbow, which is the critical factor to decide the helium flow rate in the HTR-PM. Unfortunately, these studies are limited only to global flow field structure. Little attention has been paid to the local instantaneous information of turbulent flow in 90° elbow pipes.
Under general operation conditions (250°C, 7.0 MPa), the velocity of helium in the riser pipes of HTR-PM reaches 24.69 m/s, to which , the corresponding Reynolds number, is about 1.3 × 106. Considering the limitation of the drive capability of the pump arranged in the experimental loop, too large Reynolds number cannot be achieved in our experiments. So the focus of this paper is concentrated on the detailed information on instantaneous pressure field and velocity field through a circular-sectioned 90° elbow but not on simulating the working condition in the HTR-PM. The purpose is achieved by applying the large eddy simulation approach coupled with experimental data and RANS results.
In the following, the experimental apparatus and numerical description are explained in detail. The accuracy of the numerical method is demonstrated by comparing numerical and experimental findings. Detailed information on the instantaneous flow field is analyzed by a detailed discussion of the instantaneous pressure and velocity distribution on the cross section of the elbow.
2. Experimental Apparatus
The compressibility of helium can be left out of account due to the relatively low velocity in the practical situation of HTR-PM, to which the corresponding Mach number is less than 0.3. Hence, the incompressible helium flow can be imitated by the water flow in specific test loop.
The experiments were conducted on a self-built water circulating loop. The experimental section is illustrated in Figure 2. The configuration of the test pipe is circular and its inner diameter, , is 80 mm. It is made up of a 90° elbow of curvature radius mm and two straight pipes in the upstream tangent and downstream tangent with the length of . The upstream and downstream tangents are arranged long enough to eliminate the disturbance generated by chokes or bends upstream and to guarantee the full development of turbulent flow downstream. The ratio of bending radius to diameter is 1.5 to mimic the geometric construction of elbows utilized in HTR-PM. The test section was arranged in vertical plane. The water flows along the horizontal and then turns to vertical due to the constraints of the elbow. The circulation is sustained by the constant pressure difference between the elevated water tank and the storage tank at lower place to guarantee the stable and continual flow. The flow rate is adjusted by the electric control valve located downstream of the test section and is measured by an electromagnetic flowmeter installed on the horizontal pipe. The water in the storage tank is pumped to the elevated water tank by a variable frequency pump.
There are 8 pressure tappings scattered along the intrados and extrados of the elbow at different bending angle and 4 tappings scattered on the inner and outer side of the downstream tangent, as illustrated in Figure 2. The pressure difference, , between each of the 12 tappings and the identical pressure reference point was measured by identical Rosemount 3051S pressure difference transmitter, whose accuracy was within 0.04%. The pressure reference point is located 4 upstream of the inlet of the elbow. Also, the system pressure and temperature were measured to decide the physical property parameter of water.
Measurements were performed for a water flow with the bulk mean velocity m/s, to which the corresponding Reynolds number was 44805. Special attention has been paid to decreasing the measuring error, like averaging the observed values of different measurements and maintaining the stability of the power supply.
3. Numerical Description
3.1. Geometrical Model and Mesh Division
The computational domain extracted from the experimental apparatus has the same geometric parameters with the test section except the length of upstream and downstream tangents, which were shortened to and , respectively, to reduce the computational amount. This reduction can be acceptable owing to the undisturbed flow upstream in the numerical case, and it was proved to be applicable by comparing the results to the experimental data.
Great care must be taken when setting the size, , of the wall-normal mesh near the wall because the flow characteristics in the boundary layer near the wall may interact with the flow separation at an intense level. Eguchi et al.  claimed that, under the condition of a fully developed velocity profile, the relation between Re and the pipe diameter in terms of the wall unit, , was described as below:
Here, the constants and are the same as those in the widely used wall law, that is, typically 0.4 and 5.5. is the pipe diameter in terms of the wall unit and is obtained by where is the friction velocity and is the kinematic viscosity. The required mesh size normalized by the pipe diameter can be roughly estimated asaccording to Piomelli and Balaras . Then the required mesh size is obtained as illustrated in (1) and (3), which suggests that should be approximately 8 × 10−4 at Re = 44805 for a resolved LES computation. Therefore, the mesh size normal to the wall was set to 0.0001 in the present mesh design to satisfy the requirements of LES via ICEM CFD software. The entirety of the domain was divided into about 3.13 million hexahedral elements finally, as illustrated in Figure 3.
In this paper, the main attention is paid to the influence of the large scale vortexes induced by the elbow on the instantaneous flow field. Hence, the mesh scales are designed elaborately to catch the main information of large scale vortexes. check was made after the calculation. The vast majority of is less 1 and the maximum is 2.5 only distributed on minority of the intrados of the elbow. By comparing with the experimental results, the mesh strategy was proved to be fine enough to capture the pressure characteristics adjacent to the wall.
3.2. Numerical Model and Boundary Conditions
The large eddy simulation method embedded in ANSYS FLUENT 14.5 software was used to solve the incompressible Navier-Stokes equations for the instantaneous flow field. In LES, large eddies are resolved directly, while small eddies are modeled. The governing equations for LES employ the Boussinesq hypothesis as in the RANS models. The dynamic Smagorinsky-Lilly model was selected to model eddy viscosity because of the fact that the Smagorinsky model constant, , can be dynamically computed based on the information provided by the resolved scales of motion. The dynamic procedure therefore obviates the need for users to specify the model constant in advance. Also, the RNG model based on the second-order center-difference scheme was utilized to calculate the RANS results which was compared with the LES results.
We have assumed steady inlet velocity without any temporal fluctuation superimposed and the spatial variations of velocity were taken into consideration. This is mainly because the inlet disturbance is expected to be low due to the long straight pipe arranged before the elbow. For purpose of comparing the numerical results with the experimental data, the boundary conditions of the simulation cases were set identically with the experimental cases. A steady state flow in a 12.5 long straight pipe with the same diameter and inlet velocity was also solved by using the RNG model. The velocity profile on the cross section located 2 upstream of the outlet of the straight pipe was set as the velocity inlet boundary conditions for the inlet of the elbow. A nonslip boundary condition was employed for the wall boundary and the pressure outlet was inducted for the outlet. There is no heat transfer happening in the computational domain.
Water at ambient conditions (25°C, 0.1 MPa) was used as the working fluid ( kg/m3 and kg/(m·s)) in the simulation cases. Turbulent intensity () based on the elbow pipe diameter was specified for an initial assumption of turbulent quantities.
4. Results and Discussion
4.1. Longitudinal Distribution of Wall Static Pressure
Nondimensional pressure difference normalized by 0.5, namely, the pressure coefficient, , was introduced to describe the longitudinal distribution of wall static pressure. The wall static pressure at various values gained from LES method, experiments, and RANS method is illustrated against the locations along the flow direction in Figure 4. The results of LES method were time-averaged to meet the requirement for comparison. Figure 4 shows that the results of LES method are in good agreement with that of experiments and RANS method between and 90°. The accuracy of the LES results is verified and more effective information of the flow field can be extracted from the LES results.
Upstream of and downstream of , the LES results are relatively smaller than that of RANS method but more close to the experimental data. The deviation between the three methods is induced by their own characteristics. For experiments, the tappings drilled in the wall change the surface structure to some degree. This exerts an influence on the boundary layer flow. RANS method solves the time-averaged equations of motion for fluid flow and loses some of the transient information. However, the LES method takes the interactions among the various scale vortexes into consideration and accords with the actual conditions.
Under the impact of the centrifugal force produced throughout the swerve, the extrados is squeezed and the intrados is sucked, both by the fluid flowing outwards to the outer wall. This induces a pressure swell on the extrados and leads to a pressure sink on the intrados along the flow direction, just as shown in Figure 4.
Before = −1.5, the flow is not yet disturbed by the elbow and the wall static pressures on the inner wall and the outer wall remain equal. However, the static pressure on extrados exceeds that on intrados rapidly after this location. Temporary pressure increase and pressure decrease occur, respectively, on the extrados and intrados. At the location of on extrados, the static pressure is up to its maximum. Hence, the fluid near the extrados between the location of = −1.5 and that of is decelerated slightly due to the adverse pressure gradient. Simultaneously, at the location of , the static pressure on the intrados reaches its minimum and the fluid between the location of = −1.5 and that of is accelerated in accordance with the positive pressure gradient.
Downstream of the location of , the distorted wall static pressure begins to recover and ends up with a uniform pressure field. Then between the location of and that of , the fluid near the extrados is accelerated slightly because of the positive pressure gradient, while, the fluid near the intrados is decelerated according to the adverse pressure gradient between the location of and that of .
On the intrados, between and 90°, the static pressure has a little subsidence, as marked in blue in Figure 4. The abnormality is caused by the flow separation occurring near the intrados. The recirculation zone has two stagnation points, of which the velocities are both zero. And at the middle of the recirculation zone, the back-flow velocity is relatively high. Then the pressure here is lower than other locations in the zone according to the Bernoulli equation.
4.2. Circumferential Distribution of Wall Static Pressure
Wall static pressure on specific locations is extracted from the LES results. The schematic diagram of the distribution of the pressure monitors on the circumference of the cross section at different bending angle is illustrated in Figure 5. There are 8 monitors distributed on each circumference. Every two of them, which are located at the same level, form a symmetric point pair and are connected to each other by a horizontal center line, just as shown in Figure 5.
The fluctuation of pressure on the wall makes it especially hard to compare the instantaneous pressures on the monitors at a specific time, while, a global time-averaged value of LES results will lose all of the instantaneous information. So a short-time-averaged method was introduced to make the comparison realizable. In this paper, we randomly intercepted a continuous data flow with a length of 2 seconds from the integral LES results. The intercepted data was averaged over time and this conduct made it feasible to compare the instantaneous pressures without losing instantaneous information. The comparison of the static pressure on the symmetric point pair is shown in Table 1, where and represent the short-time-averaged results of the pressure coefficient on the monitors distributed on left and right semicircle, respectively, and is the mean of and . The relative deviation is described as
For the cross section sliced at , 48°, 60°, and 78°, the relative deviation related to the symmetric point pairs on the upper semicircle is quite small. This indicates that the instantaneous pressure distribution on the upper semicircle of the cross section is symmetrical. However, related to the symmetric point pairs on the lower semicircle is relatively large, especially to the pairs whose is 210° and 330°, as marked by bold font in Table 1. related to this pair is about 3~10 times larger than that related to other pairs and decreases gradually along with the increasing of . The asymmetry of the instantaneous pressure distribution on the lower semicircle of the cross section is caused by the counterrotating vortexes with different swirling intensity and the flow separation occurring near the intrados. The separation grows stronger gradually along with the flow developing downstream and interacts with the counterrotating vortexes. Large scale vortexes are broken into small scale vortexes. The momentum and energy are transferred among them and the fluctuation energy is dissipated rapidly, which is conductive to the recovery of the distorted pressure field.
The pressure contours on the cross section sliced at , 75°, and 90° are illustrated in Figure 6. The data was extracted at the time s. The static pressure value demonstrated here refers to the pressure of the outlet of the domain. Corresponding to the former conclusion, the pressure field of the cross section is symmetrical in the upper semicircle and asymmetrical in the lower semicircle.
Figure 6 indicates that the zone of disturbed pressure distribution spreads towards the upper part along with the increasing of . And the asymmetrical area of the instantaneous pressure field extends gradually. Simultaneously, the contours thin and the difference between the maximum and minimum of the contours diminishes gradually. In conclusion, the instantaneous pressure distribution on the cross section is no longer globally symmetrical due to the secondary flow and the flow separation. The influence of flow separation on the flow field focuses on the lower semicircle. Along with the development of flowing, the instantaneous pressure distribution has a tendency to homogenize.
4.3. The Velocity Vector Fields on the Cross Section
The velocity vector fields projected onto the cross sections sliced at , 75°, and 90° are illustrated in Figure 7. The fluid near the axis of symmetry is forced to flow towards the extrados under the impact of the centrifugal force, while the fluid near the wall of left and right semicircles is transported to the intrados along the circumferential wall under the differential pressure. Then, the large scale vortexes rotating in contrary direction come up in the elbow, known as the Dean vortex.
On the lower semicircle of the cross sections sliced at and 75°, two counterrotating vortexes with small scale are observed. These vortexes are induced by the relatively high intensity secondary flow near the intrados. Figure 6 shows that a low pressure area presents at the intrados. Then the fluid adjacent to the left and right semicircle wall is accelerated by the pressure difference near the intrados. When it comes to the appropriate location at the lower semicircle, the fluid deflects towards extrados rapidly and flows away. In the process, the fluid that is deflecting with relatively high velocity rotates the neighbouring fluid through the shearing action and promotes the generation of the two vortexes in the end. The cores of the two vortexes are of axial symmetry but the intensity of them has some difference. On the cross section at , the two vortexes interact with each other to some degree. Their velocity fields begin to merge and the location of the cores moves to an upper level comparing to that of the section at .
The flow field of the lower semicircle of the cross section at is more complicated than that upstream. Secondary vortexes with smaller scale are observed in this section. These vortexes are the outcome of the interaction between the secondary flow and the separation vortexes. The separation vortexes are broken into small eddies by the shearing action. Then small eddies merge while the others break with certain probability. Actually, the probability of eddy breakdown is much higher than that of eddy merging. Kinetic energy is transferred from eddies with relatively large scale to eddies with small scale and this transfer process continues until the rotational energy is dissipated by the viscous stress.
Both of the pressure field and velocity field are of homogeneous distribution on the upper semicircle of the cross section. Near the extrados, the pressure is relatively high due to the squeeze effect. Then the velocity of this region is relatively low. In addition, the pressure gradient adjacent to the extrados is lower than that adjacent to the intrados. Hence, the flow near the extrados is much more stable and this is the primary cause of the homogeneous flow field on the upper semicircle of the cross section.
The LES method combined with RANS method and experiments was applied to research the complicated flow field in a 90° elbow and the instantaneous turbulence information was obtained. The results were clarified to be reliable by comparing with that of RANS method and experiments. The conclusions obtained in this study are summarized as follows:(1)The instantaneous flow field in a 90° elbow is globally asymmetrical. The pattern and intensity of the Dean vortex and the small scale eddies change over time and induce the asymmetry of the flow field.(2)The counterrotating vortexes and various scale eddies induced by shearing action exist in the region adjacent to the intrados.(3)The turbulent disturbance upstream and the interaction between the secondary flow and the flow separation are the dominant factors distorting the flow field.(4)The flow field on the upper semicircle of the cross sections is homogeneous due to the relatively lower velocity and pressure gradients of the fluid adjacent to the extrados.
The authors declare that they have no competing interests.
This work is carried out in the frame of a current research project funded by the National Natural Science Foundation of China under Grant no. 51506114, the Specialized Research Fund for the Doctoral Program of Higher Education under Grant no. 20130002120015, and the Science Fund for Creative Research Groups of National Natural Science Foundation of China (Grant no. 51321002).
- S. Y. Jiang, X. T. Yang, Z. W. Tang et al., “Experimental and numerical validation of a two-region-designed pebble bed reactor with dynamic core,” Nuclear Engineering and Design, vol. 246, pp. 277–285, 2012.
- R. Schulten, “Pebble bed HTRs,” Annals of Nuclear Energy, vol. 5, no. 8–10, pp. 357–374, 1978.
- A. C. Kadak and M. V. Berte, “Modularity in design of the MIT pebble bed reactor,” in Proceedings of the American Nuclear Society Meeting, Reno, Nev, USA, November 2001.
- A. Koster, H. D. Matzner, and D. R. Nicholsi, “PBMR design for the future,” Nuclear Engineering and Design, vol. 222, no. 2-3, pp. 231–245, 2003.
- Z. Y. Zhang, Z. X. Wu, Y. L. Sun, and F. Li, “Design aspects of the Chinese modular high-temperature gas-cooled reactor HTR-PM,” Nuclear Engineering and Design, vol. 236, no. 5-6, pp. 485–490, 2006.
- Z. Zhang, Z. Wu, D. Wang et al., “Current status and technical description of Chinese 2 × 250 MWth HTR-PM demonstration plant,” Nuclear Engineering and Design, vol. 239, no. 7, pp. 1212–1219, 2009.
- D. C. Taylor, M. B. McPherson, and E. B. Meier, “Elbow meter performance [with Discussion],” Journal of the American Water Works Association, vol. 46, no. 11, pp. 1087–1095, 1954.
- K. Rup and L. Malinowski, “Fluid flow identification on base of the pressure difference measured on the secant of a pipe elbow,” Forschung im Ingenieurwesen, vol. 70, no. 4, pp. 199–206, 2005.
- L. Malinowski and K. Rup, “Measurement of the fluid flow rate with use of an elbow with oval cross section,” Flow Measurement and Instrumentation, vol. 19, no. 6, pp. 358–363, 2008.
- S. A. Berger, L. Talbot, and L. S. Yao, “Flow in curved pipes,” Annual Review of Fluid Mechanics, vol. 15, pp. 461–512, 1983.
- H. Ito, “Flow in curved pipes,” JSME International Journal, vol. 30, no. 262, pp. 543–552, 1987.
- R. J. S. Pigott, “Pressure losses in tubing, pipe and fittings,” Transactions of the ASME, vol. 72, pp. 679–688, 1950.
- H. Kirchbach, Loss of Energy in Miter Bends, Transactions of the Munich Hydraulic Institute, American Society of Mechanical Engineers, 1935.
- W. Schubart, “Energy loss in smooth and rough surfaced bends and curves in pipe-lines,” Transactions of the Hydraulic Institute of the Munich Technical University, vol. 3, pp. 81–99, 1935.
- M. M. Enayet, M. M. Gibson, A. M. K. P. Taylor, and M. Yianneskis, “Laser-Doppler measurements of laminar and turbulent flow in a pipe bend,” International Journal of Heat and Fluid Flow, vol. 3, no. 4, pp. 213–219, 1982.
- K. Sudo, M. Sumida, and H. Hibara, “Experimental investigation on turbulent flow in a circular-sectioned 90-degree bend,” Experiments in Fluids, vol. 25, no. 1, pp. 42–49, 1998.
- Y. Eguchi, T. Murakami, M. Tanaka, and H. Yamano, “A finite element les for high-Re flow in a short-elbow pipe with undisturbed inlet velocity,” Nuclear Engineering and Design, vol. 241, no. 11, pp. 4368–4378, 2011.
- F. Rütten, M. Meinke, and W. Schröder, “Large-eddy simulation of 90° pipe bend flows,” Journal of Turbulence, vol. 2, 2001.
- F. Rütten, W. Schröder, and M. Meinke, “Large-eddy simulation of low frequency oscillations of the Dean vortices in turbulent pipe bend flows,” Physics of Fluids, vol. 17, no. 3, Article ID 035107, 2005.
- M. J. Tunstall and J. K. Harvey, “On the effect of a sharp bend in a fully developed turbulent pipe-flow,” Journal of Fluid Mechanics, vol. 34, no. 3, pp. 595–608, 1968.
- A. Ono, N. Kimura, H. Kamide, and A. Tobita, “Influence of elbow curvature on flow structure at elbow outlet under high Reynolds number condition,” Nuclear Engineering and Design, vol. 241, no. 11, pp. 4409–4419, 2011.
- Y. Eguchi, T. Murakami, M. Tanaka, and H. Yamano, “A finite element les for high-Re flow in a short-elbow pipe with undisturbed inlet velocity,” Nuclear Engineering & Design, vol. 241, no. 11, pp. 4368–4378, 2011.
- U. Piomelli and E. Balaras, “Wall-layer models for large-eddy simulations,” Annual Review of Fluid Mechanics, vol. 34, pp. 349–374, 2002.
Copyright © 2016 Shiming Wang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.