Abstract
High efficiency and low fuel consumption make the contrarotating open rotor (CROR) system a viable economic and environmentally friendly powerplant for future aircraft. While the potential benefits are well accepted, concerns still exist with respect to the vibrations and noise caused by the aerodynamic interactions of CROR systems. In this paper, emphasis is placed on the detailed analysis of the aerodynamic interactions between the front and aft propellers of a puller CROR configuration. For the first step, unsteady Reynolds-averaged Navier-Stokes (URANS) simulations coupled with dynamic patched grid technology are implemented on the isolated single-rotating propeller (SRP) configuration in various operating conditions in order to test the accuracy and feasibility of the numerical approach. The numerical results are verified by a wind tunnel test, showing good agreements with the experimental data. Subsequently, the URANS approach is applied to the CROR configuration. The numerical results obtained through the URANS approach help to improve the understanding of the complex flow field generated by the CROR configuration, and the comparison of SRP flow field and CROR flow field allows for a detailed analysis of the aerodynamic interactions of the front propeller blade wakes and tip vortices with the aft propeller. The main reason of the aerodynamic interactions is the mutual effects of the blade tip vortices, and the aft propeller reduces the strength of the blade tip vortices of the front propeller. Aerodynamic interactions will lead to the periodic oscillations of the aerodynamic forces, and the frequency of the oscillations is linked to the blade numbers. In addition, a CROR has a larger thrust and power coefficient than that of the SRP configuration in the same operating conditions. The URANS approach coupled with a dynamic patched grid method is tested to be an efficient and accurate tool in the analysis of propeller aerodynamic interactions.
1. Introduction
In modern aircrafts, besides supersonic and transonic mainline airliners, regional airliners have the properties of smaller weight and size, lower flight speed and altitude, and good maneuverability; as a result, propeller-driven aircrafts can meet these requirements very well. Due to their shorter take-off distance, larger climbing speed, longer voyage, less fuel consumption, and lower requirements for runways, propeller-driven aircrafts play an important role in both civil and military applications [1–3]. Along with people’s attentions to environment issues in recent decades, countries around the world are calling for reducing emissions from aircrafts; therefore, the demands of novel and efficient propeller systems have greatly increased. It is on this premise that CROR systems arouse people’s interest as a possible fuel-efficient and environmentally friendly propulsion system for future aircrafts [4, 5].
The efficiency of the CROR system mainly comes from the reduction of the front propeller slipstream swirl by the impingement of the aft propeller slipstream, which will also lead to the extra axial acceleration of the flow downstream when it passes through the aft propeller; therefore, CROR systems can generate greater thrust [6–9]. However, the aerodynamic interactions caused by the front propeller blade tip vortices impinging on the blade tip vortices generated by the aft propeller will result in the noise emissions as well as vibrations [10, 11], making it a main problem for the development of the CROR system in the application of modern aircraft [12–15]. What is more, the mutual interactions between airframe components and the CROR have a remarkable impact on the aerodynamic performance of a CROR. Such installation effects will considerably influence the inflow of the propeller compared to an isolated CROR [16–18]. Some key issues need to be addressed and solved, such as figuring out the flow field around the propeller system and analyzing the slipstream flow interactions between the two propellers in order to design a good CROR system [19, 20].
In recent years, the rapid developments of computer hardware and numerical computation approach have enabled the computational fluid dynamics (CFD) method to become one of the most effective methods to simulate and analyze propeller slipstream flow field [21]. When compared to the traditional wind tunnel test, a CFD method is more economic and time-saving and more importantly, the CFD method can capture the details of flow field while the wind tunnel test cannot [22]. Propeller slipstream flow field is rather complicated due to the strongly perturbed inflow by the blades’ periodical affection. Propeller slipstream has unsteady properties; therefore, unsteady numerical approaches should be used in the investigation of propeller slipstream flow field. The URANS calculations have been implemented in the numerical simulations of propeller slipstream, but the mechanism analysis of the complex aerodynamic interactions of CROR systems needs to be further enriched.
In this paper, the URANS code coupled with a dynamic patched grid method is first used on an isolated SRP configuration with experiment verification. The same numerical method is applied to an isolated CROR configuration in different operating conditions with the specific focus on analyzing the mechanism of the complex flow field and the aerodynamic interactions between the front propeller and the aft propeller [23, 24]. Then, some conclusions are drawn with respect to the CROR flow field and aerodynamic interactions.
2. SRP Numerical Simulations with Experiment Verification
2.1. SRP Wind Tunnel Test
An isolated SRP configuration wind tunnel test was conducted in the DNW-LST wind tunnel in Netherland in 2014. The main purpose of the wind tunnel test was to gain access to the thrust and power coefficient as well as the efficiency of the isolated propeller in different rotor speeds under certain blade pitch angles in various operating conditions. The DNW-LST wind tunnel is a low-speed return-flow wind tunnel, which has a 3.0 m (width) × 2.25 m (height) × 8.75 m (length) test section, as shown in Figure 1. The wind tunnel test model was using a 1/6th scale model that consisted of a nonrotating center shaft and a rotor equipped with 6 blades of the same blade geometry at a fixed diameter of m with a faring in front of the center shaft. The wind tunnel test model was powered by an air turbine that could drive the rotor speeds to at most 6400 rpm, as shown in Figure 2. The wind tunnel test of the 1/6th scale SRP configuration provides a broad scope of database to evaluate the numerical approaches’ ability to predict the aerodynamic performance of the SRP configuration as well as the related CROR configuration. Specific emphasis was placed on the flow conditions of the wind tunnel test in order to accommodate the numerical simulations; the flow speed of the wind tunnel test could be set ranging from 1.5 m/s to 80 m/s, and the flow precision can be controlled in the range of 0.2%. The yaw angle of the horizontal and vertical flow is less than 0.1°, while the turbulence level is controlled between 0.02% and 0.03%. The center shaft was connected to a strut, which was connected to an external balance. In this particular wind tunnel test, attentions were paid to aerodynamic performance by installing a rotated shaft balance inside the faring. Thrust coefficients of different advance ratios are obtained via the rotated shaft balance. The results obtained by the rotated shaft balance were verified through the comparison of the results obtained by the external balance.
The experiment Mach number was set to at an angle of attack °. Several wind tunnel test conditions were conducted due to various rotor speeds, listed in Table 1. In this wind tunnel test, blade pitch angle was defined as 70% radius blade pitch angle, and in these series of experiments, the blade pitch angle was set as °.
The propeller thrust coefficient, power coefficient, and efficiency have the following forms: where is the thrust generated by the propeller and , along with denote the air density, propeller rotational speed, and propeller diameter, respectively. In addition, denotes the shaft power and denotes the torque, whereas denotes free stream velocity and is the efficiency.
2.2. SRP Numerical Simulation Configuration
The SRP numerical simulations mainly focus on the investigation of the propeller flow field and blade aerodynamics, and the most important aspect is the validation of the URANS approach with the wind tunnel test. For the numerical studies, a configuration for numerical simulations was approximated as closely to the wind tunnel test model as possible. Modifications were made to the numerical simulation configuration to ensure that there was no adverse impact on the performance of the propeller as well as the flow field in their direct vicinity. The length of the center shaft was chosen to be (2 times of the propeller diameter) whereas the strut to which the SRP was connected in the wind tunnel test was removed. The numerical SRP configuration is shown in Figure 3.
(a) Side view
(b) Front view
2.3. Grid Generation
The accuracy and reliability of solutions obtained by means of CFD have close ties to the density and quality of the computational grids used. A dynamic patched grid method has the primary advantage that the grids of each subdomain can be generated separately without concerning the topology relations among adjacent blocks. One-to-one correspondence of grid points is not required on the two sides of the patched surfaces, and fluxes are interpolated near the patched surfaces in the process of computation; therefore, the flow information of adjacent subdomains is coupled and exchanged. For the numerical simulation described here, the baseline SRP configuration grids have one computational domain and one rotating domain. The computational domain, which is also called far-field, is a cylindrical region that has ( is the propeller diameters) long upstream and long downstream of the propeller, and the diameter of the computational domain is approximately , as shown in Figure 4.
The entire grids are included in the computational domain. Except the grids in the rotating domain, the other grids in the computational domain remain static. The rotating domain is a cylindrical region that contains the propeller, and the grids are generated in the cylindrical region, while fluxes are interpolated through the cylinder surfaces, which are patched surfaces, between the grids in the computational domains and rotating domain. The sketch map of the patched grid approach is shown in Figure 5. The rotating cylindrical domain is marked in purple.
The block-structured grids of this SRP configuration have been created with the software ANSYS ICEM CFD. The whole computational grids consist of two parts: the blocks around the rotating propeller and the outer blocks around the static center shaft. Attentions should be paid to the grid generation of the blade passages to make sure that each of the blade passages is created as equal as possible to allow for a reliable spatial resolution for the blades of the propeller. Refined grids are applied in the boundary layer and propeller slipstream region to improve blade wake and tip vortex resolution. The surface grid of the SRP configuration is shown in Figure 6(a). The grid shown in grey in Figure 6(b) is an enlarged image of a far-field section grid. The grid patched in the rotating domain is the propeller section grid, and the grid shown in pink is the propeller section grid. The complete block-structured grids of the isolated SRP configuration consist of 6,664,206 nodes.
(a) Surface gird
(b) Far-field section grid
2.4. URANS Approach
For the simulations described here, URANS calculations are applied to the isolated SRP configuration in order to testify the feasibility and accuracy of this numerical method, spatial discretization of the convective fluxes is done using a third-order upstream MUSCL scheme, and the viscous fluxes are discretized with central differences. A Spalart–Allmaras (SA) model is chosen as the turbulence model. A fully implicit LU-SGS method with subiteration in pseudotime is implemented as a time-marching method, and the well-established dual-time approach is used to enhance the precision of time discretization. In the course of the URANS simulations, the time step size is set equivalent to a propeller rotation of deg per physical time step and the numerical results are obtained after 8 full revolutions of the SRP propeller. The numerical approach is implemented on 4 cases, whose numerical parameters are in accordance with the wind tunnel test. The SRP numerical simulations are considered converged if the propeller thrust coefficients show no notable differences. The SRP simulations were run on 64 processors of Inspur TS10000 Cluster in Center for High Performance Computing, Northwestern Polytechnical University. For each case, the approximate total wall-clock runtime for the simulation of the 8 full revolutions was 1 day.
2.5. Result Comparison
The comparison of the SRP wind tunnel test and URANS approach in terms of thrust coefficient and power coefficient is shown in Figure 7, and the results obtained through numerical simulations show a good agreement with the wind tunnel test.
(a) Thrust coefficient
(b) Power coefficient
(c) Efficiency
The detailed results of the wind tunnel test and numerical simulations are listed in Table 2. The relative errors of propeller thrust and power coefficient of the 4 operating conditions are quite small, which means that the grid generation and URANS numerical approaches are validated to be accurate in the investigation of the SRP configuration.
3. CROR Numerical Simulations
3.1. CROR Numerical Simulation Configuration
In order to demonstrate the capability of previously used URANS simulations in the analysis of aerodynamic phenomena and interactions of a CROR, an isolated 6 × 6 puller CROR configuration is put forward built on the isolated SRP configuration mentioned above by installing another propeller of the same blade geometry and blade numbers behind the propeller of the SRP configuration, as shown in Figure 8. This CROR configuration is only fictitious and has never actually been built. The geometric parameters of the CROR are consistent with the SRP configuration. The length of the center shaft is , while the distance between the two propellers is . The blade pitch angles of the front and the aft propellers are both set as °. The rotational speed is set to be identical for both propellers, and the two propellers rotate in the opposite directions, counterclockwise for the front propeller and clockwise for the aft one. The main reason for using the same blade numbers for a CROR with the identical diameters and blade pitch angles is that the CROR configuration is built on the SRP configuration. The mechanism of aerodynamic interactions of a CROR should be compared with the SRP. The impacts of blade numbers, blade diameters, and blade pitch angles as well as rotor speeds on the aerodynamic interactions have been researched, and the research papers will be released later [25].
(a) Side view
(b) Front view
3.2. Grid Generation
In the investigation of a CROR, the grid generation approach placed primary emphasis on the adequate resolution of the flow phenomena that played a significant role in the interactions between the front and aft propellers. Thus, the main effort goes into obtaining high-quality and high-resolution propeller grids to address requirements for good aerodynamic analysis results, which are a significant aspect for the design of a high efficient CROR system. The surface grid and section grid of the CROR configuration are shown in Figure 9. The rotational domains of the front and aft propellers are created of the same size to ensure identical grid densities and qualities. Blade passages are created as equal as possible in order to capture the periodic nature of the CROR configuration. The CROR numerical simulations mainly focus on the investigation of the propeller flow field and aerodynamic interactions between the two propellers. Refined grids are applied in the propeller slipstream region of the CROR configuration and designed to capture the flow features especially the interactions between the front and aft propellers in terms of blade wake and tip vortices during the process of numerical simulations. The grids used in this investigation are block-structured and generated through ANSYS ICEM CFD and consist of 17,723,511 nodes.
(a) Surface gird
(b) Far-field section grid
3.3. URANS Approach
In the investigation of the CROR configuration, all the numerical settings are aligned with the settings implemented in the SRP simulations. The CROR flow field is obtained after 12 full propeller revolutions, and the simulations are considered converged if the periodic fluctuations of the blade loadings for the 4 test cases show no obvious changes. The settings of CROR numerical simulation are listed in Table 3. The CROR simulations were run on 64 processors of Inspur TS10000 Cluster in Center for High Performance Computing, Northwestern Polytechnical University. For each case, the approximate total wall-clock runtime for the simulation of the 12 full revolutions was 5 days.
4. CROR Aerodynamic Interactions
4.1. Flow Field Analysis
The isosurfaces of the -criterion of the isolated SRP configuration and CROR configuration are shown in Figure 10 at the advance ratio , which demonstrate numerical simulations of the configuration flow field characterized by Mach number. The vorticity distributions on both sides of the patched surfaces did not break or dislocate; hence, the blade tip vortices and blade wakes smoothly got through the patched surfaces, which indicated the proper applications of the dynamic patched technology in the grid generation of SRP and CROR configurations and validated the numerical approach in terms of URANS used in the investigation. For the SRP configuration, as shown in Figure 10(a), the free stream passed through the propeller, then twisted by the rotating blades, and formed a helical vortex-tube downstream. While for the CROR configuration, the two propellers rotate in the opposite directions, the blade wakes and tip vortices of the front propeller are “chopped” by the blade wakes, and tip vortices generated by the aft propeller and an interwoven net-structured vortex system were formed as shown in Figure 10(b), which lead to periodic unsteady oscillations. Iso- surfaces reveal the interactions of the front propeller blade tip vortices and blade wakes with the aft propeller. The front blade wakes played a dominant role in the interaction with the aft blades, so the unsteady flow structures are mainly characterized by the wake impingement of the front propeller on the aft propeller. What should be emphasized is that the blade tip vortices of both propellers are almost perpendicular which happened in each of the test cases regardless of the propeller rotational speeds; the blade tip vortices of the front propeller progressively merged with the tip vortices produced by the aft propeller yielding the net-structured pattern observed downstream.
(a) SRP flow field
(b) CROR flow field
4.2. Blade Tip Vortex Development
In order to investigate and analyze the vortex systems in depth, instantaneous value of the vorticity magnitude of the cross section is plotted, as shown in Figure 11. The position of the blade tip vortices can be seen clearly for both SRP and CROR configurations by the red spots of vorticity contours. Clearly the blade wakes of each propeller are visible in the slipstream as well as their deformation with respect to the flow acceleration. Visible in Figure 11(a), the blade tip vortices generated on the SRP shows no obvious dissipation as they develop downstream and the diameter of slipstream remains almost the same. However, contraction of slipstream can be seen from Figure 11(b), as the tip vortex emanating from the front rotor interacts with the vortex from the second rotor with the flow develops downstream. According to the law of conservation of mass, the mass flow passing through the slipstream tube remained unchanged and the flow velocity increased, which lead to the decrease of the radius of the slipstream tube and hence the contraction of slipstream. Visible in the plot of vorticity magnitude contour in Figure 11(b), compared to the aft propeller, the front propeller shows a much more rapid decay in the blade tip vortices. Generally, the aerodynamic interactions between the two propellers result in the obvious contraction of the slipstream and the rapid decay of the blade tip vortices of the front propeller. The vorticity magnitude greatly increased in the hub downstream region. As the flow moves downstream, with an increasing axial distance from the propellers, blade wakes and tip vortices are gradually dissipating.
(a) Vorticity magnitude of the SRP at longitudinal position
(b) Vorticity magnitude of the CROR at longitudinal position
The blade tip vortex development of SRP and CROR configurations is plotted in terms of vorticity magnitude at axial positions of , , and as shown in Figure 12. As shown in Figure 12(a), the vorticity level is gradually decreasing at axial positions of , , and (marked as black, red, and blue solid lines) with the peak vorticity locates approximately in the radial position of . While in Figure 12(b), the black solid line shows the interception of the front blade tip vortex before it impinges on the blade tip vortex generated by the aft propeller at axial position of , while the red solid line demonstrates the interception of the front and aft blade tip vortices at axial position of . At , the peak vorticity level of the front blade tip vortex locates in the radial position of , and when the front blade tip vortex interacts with the aft blade tip vortex at axial position of , the peak vorticity of front blade tip vortex locates in the radial position of . The peak vorticity of front blade tip vortex at is lower than that of the front blade tip vortex at axial position of (this is due to the reduction in the front blade tip vortex strength by the aft blade tip vortex) and validates the contraction of slipstream. Also, visible in Figure 12(b), the peak vorticity level of the aft blade tip vortex locates in the radial position of and is lower than that of the front blade tip vortex at ; this means that the aft propeller has a function of tip vortex dissipation.
(a) Blade tip vortex development of the SRP
(b) Blade tip vortex development of the CROR
4.3. Aerodynamic Characteristics
Unsteady thrust coefficients are plotted in order to analyze the aerodynamic interactions between the two propellers, as shown in Figure 13. These unsteady blade loadings can be easily obtained from numerical simulations. Figure 13 shows the front and aft propeller thrust coefficient development during one full revolution with the rotational speed of 3290/3290 rpm marked as solid purple and dash-dotted blue lines, respectively. Visible in Figure 13, the thrust coefficients of the two propellers have unsteady but periodic oscillations with respect to aerodynamic interactions between the front and aft propellers. During one full revolution, the thrust coefficient of each propeller will fluctuate twelve times and the thrust coefficient of the two propellers reached maximum or minimum almost at the same time. Due to the 6 × 6 puller CROR configuration studied here, the oscillations are periodic, recurring every 30 deg of propeller rotation. Strong oscillations of the propeller performance are observed for the aft propeller. For the front propeller, the unsteady loadings are a result of the aft propeller flow field; as for the aft propeller, the much stronger oscillation is caused by the impingement of the front blade wakes and tip vortices. The front propeller delivers a mean thrust coefficient of , while the aft propeller delivers a larger mean thrust coefficient of . Compared to the CROR configuration, the propeller of the SRP configuration has a relatively constant blade loading in terms of thrust coefficient, which delivers a value of marked as red a solid line in Figure 13.
Figure 14 shows the time-averaged plots during one full rotation of the unsteady numerical simulations in terms of Mach number and axial velocity magnitude with the particular focus on the investigation of flow acceleration. As shown in Figure 14(a), the value of Mach number barely increased when the flow passes through the propeller of the SRP configuration. The plot of the time-averaged axial velocity of one full rotation, normalized with free stream Mach number at certain axial positions of , , and , indicates that slipstream is accelerated slightly as it develops downstream, as shown in Figure 14(b). As shown in Figure 14(c), the value of Mach number increased when the flow passes through the front propeller, and Mach number increases a second time when the flow passes through the aft propeller of the CROR configuration. The acceleration of the flow by the front propeller as well as the further subsequent increase in the Mach number of the aft propeller can be seen clearly from Figure 14(c). Figure 14(d) shows the plot of time-averaged axial velocity magnitude for one complete revolution, normalized with free stream Mach number, at axial positions of , , and downstream of the front propeller. The black solid line indicates that the flow is first accelerated when it passes through the front propeller in terms of the increase in the axial velocity, while the red solid line shows that the flow is further accelerated when it passes through the aft propeller. The slipstream will be accelerated downstream of the front propeller which is demonstrated by the blue solid line. In addition, with the increasing distance downstream of the front propeller, the slipstream acceleration will gradually weaken. The plot of time-averaged axial velocity magnitude clearly shows that slipstream acceleration is most notable between radial positions of and .
(a) SRP cross section Mach number distribution
(b) SRP time-averaged slipstream velocity
(c) CROR cross section Mach number distribution
(d) CROR time-averaged slipstream velocity
In the numerical investigation of SRP and CROR configurations, the figures of velocity profiles are plotted in terms of time-averaged slipstream-normalized velocity components of one complete revolution. Figures 15(a) and 15(b) show the tangential velocity at the same axial positions of , , and . For the SRP, the propeller rotates counterclockwise in a right-handed coordinate system aligned with the -axis, the values of tangential velocity at the three axial positions are positive between radial positions of and , and the three curves are almost coincided, which means that the tangential velocities remain nearly unchanged at each axial position. For the CROR, the front propeller rotates counterclockwise and the aft propeller rotates clockwise in a right-handed coordinate system aligned with the -axis. This brings about the fact that the tangential velocity is positive at the axial position of between radial positions of and , while the tangential velocity is negative at and between radial positions of and . The aerodynamic interactions lead to the decrease of the values between radial positions of and at and compared to the absolute value at ; this is because the front blade tip vortex encounters the aft blade tip vortex and vortex strength is reduced. The tangential velocity becomes positive near the blade tips between and at axial positions of and , because the strength of the mutual interactions of blade tip vortices has weakened and this indicates that the front propeller of the CROR plays a dominant role in the aerodynamic interactions between the two propellers.
(a) SRP velocity component
(b) CROR velocity component
(c) SRP velocity component
(d) CROR velocity component
Figure 15(c) shows the plot of time-averaged slipstream-nomalized radial velocity of the SRP, slipstream contraction results in the negative values of the radial velocity component outboard of the radial position of at the axial position of . The radial velocity component inboard of shows a positive value because of the higher pressure near the center shaft. At further downstream positions of and , the radial velocity component decreases close to zero as the flow tends to be stable, which represents the completion of slipstream contraction.
Figure 15(d) is the plot of time-averaged slipstream-nomalized radial velocity of a CROR, the maximum radial velocity component outboard of the radial position of at the axial position of is nearly two times as close as that of the SRP due to the increase of Mach number in the region between the two blades. When it comes to the downstream positions of and , the radial velocity component shows an obvious fluctuation along the radial direction because of the strong interference and the slipstream shows a more intense contraction.
Figure 16 shows the comparison of CROR and SRP thrust and power coefficient in regard to advance ratios. As shown in Figure 16, the SRP thrust coefficient and power coefficient are larger than that of the CROR front propeller but smaller than that of the CROR aft propeller in the same advance ratio; this is due to the recovery of the swirl flow behind the front propeller into axial momentum. Naturally, in the given operating conditions, the CROR can generate much larger thrust than SRP does. The efficiency of a CROR is higher than that of the SRP due to the reduction of blade tip vortices by the aft propeller.
(a) Thrust coefficient
(b) Power coefficient
(c) Efficiency
5. Conclusions
An analysis of aerodynamic interactions was conducted for a 6 × 6 puller CROR configuration by using a dynamic patched grid method-based URANS numerical simulation approach. The URANS approach has been applied to an isolated SRP configuration with wind tunnel test verification, and this numerical approach is validated to be an efficient and accurate tool in the investigation of propeller slipstream and aerodynamic characteristics. The results obtained through the URANS approach help to improve the understanding of the CROR slipstream flow field and aerodynamic interactions between the two propellers.
The isocontours of vorticity magnitude characterize the vortex system of CROR slipstream flow field, which shows the complex interaction of the front propeller blade wakes and tip vortices with the aft propeller. The mutual interactions between the two propellers result in unsteady periodic blade loading oscillations during one full rotation. Due to the equal diameter and rotational speed of the two propellers, the frequency of blade loading oscillations depends on the blade numbers. The instantaneous contours of vorticity magnitude of the CROR blade tip vortex development reveal the vortex dissipation caused by the aerodynamic interactions. Time-averaged three-component velocity profiles at certain axial positions have been analyzed, and the results show that the freestream is accelerated twice when it passes through the two propellers. With the slipstream developing downstream, the diameter of the slipstream will decrease until the completion of the contraction. The upstream front propeller blade wakes play a dominant role in the interaction with the aft propeller, and these wakes are then ingested by the wakes of the aft propeller and have a strong impact on the vortex system. Generally, a CROR will produce much larger thrust than that of a SRP under the same circumstances, and the higher the rotational speed is, the greater the power consumption will be.
Data Availability
The data of the SRP propeller wind tunnel test used to support the findings of this investigation has not been made available because it involves commercial secrets.
Conflicts of Interest
The authors declare that no conflict of interest regarding the publication of this paper.
Acknowledgments
This study is supported by the State Key Development Program for Basic Research of China (no. 2015CB755800).