Abstract
In order to study the influence of the circumferential position of the balance hole on the cavitation performance of the semiopen impeller centrifugal pump, a low specific speed semiopen impeller centrifugal pump is taken as the object, and 4 kinds of circumferential positions of balance holes are designed. The SST k-ω turbulence model and the Rayleigh–Plesset cavitation bubble dynamics equation are used to calculate the full flow field of the centrifugal pump. Research shows that, under cavitation conditions, as the circumferential position of the balance hole is farther away from the blade working surface, the cavitation performance of the pump is reduced, and the larger θ (the angle of the balance hole and the leading edge of the blade in the direction of rotation) is, the easier the jet cavitation occurs near the balance hole. On the other hand, with the development of cavitation, the axial force of the impeller has also changed greatly. In contrast, the farther the balance hole is arranged in the circumferential direction (i.e., the greater θ), the more limited is the ability of the balance hole to balance the axial force.
1. Introduction
Due to the simple processing, easy replacement, and cleaning of semiopen impeller centrifugal pumps, they are often used in the petrochemical industry, agricultural production, fire protection, and other fields and are also involved in the aerospace field. The semiopen impeller is different from the closed impeller. Because there is no front cover plate in the semiopen impeller centrifugal pump, in order to avoid friction and collision between the blade and the pump shell, it is necessary to leave a certain blade tip clearance during processing and manufacturing [1–4]. This structure makes the efficiency of the semiopen impeller centrifugal pump generally worse than the same type of closed impeller centrifugal pump, and the rotor axial force is generally larger as well. Therefore, in addition to improving pump efficiency, balancing the axial force of the semiopen impeller and improving the cavitation of the pump have become people's concerns. At present, there have been some studies on the cavitation problem of semiopen impeller centrifugal pumps. Yang [5] and others found that the degree of cavitation in blade tip clearance of the screw centrifugal pump depends on the clearance size. The closer the tip is to 1.2 mm, the more severe the clearance cavitation is, and when the clearance is less than 0.3 mm, the clearance cavitation is less likely to occur. Zhao et al. [6] studied the cavitation phenomenon in tip clearance flow with cavitation number, angle of attack, blade load, and tip clearance as the main factors and proved that the numerical results are in good agreement with the experimental results. Zhang et al. [7] have shown that, under the condition of a small flow rate, leakage vortex cavitation occurs in the tip gap, and then the cavitation is seriously deteriorated. Bubble clouds are broken at the tail of the blade, and gather in the tip area, occupying a large space in the blade passage area. Li et al. [8] used a high-speed digital camera to photograph the tip clearance cavitation flow of an axial flow pump at the flow rate of the best efficiency point. The evolution of the cavitation region and the formation of cavitation bubbles were tracked, which provided an effective reference for the tip clearance cavitation region and the development process of an axial flow pump. Liu [9] et al. conducted a series of experiments in a closed-loop cavitation wind tunnel. The cavity is analyzed by the image analysis method, and the development characteristics of cavitation flow under different gap sizes are compared. Liu et al. [10] proposed that changing the geometry of the microturbine pump impeller can reduce the effect of tip clearance. Luo et al. [11] and others designed an impeller with inclined blades. Compared with conventional impellers, they found that the semiopen centrifugal impeller with inclined blades has better cavitation performance. In the study of the influence of balance hole diameter on the cavitation performance of low specific speed centrifugal pumps, Zhao et al. [12] et al. found that for low specific speed centrifugal pumps, when the specific area (the ratio of the total area of the balance hole to the flow cross section area of the rear ring) of the balance hole is between 2.5 and 4.5, the anticavitation ability of pump and the effect of balancing the axial force are the best. Zhang et al. [13] verified the NPSH curve and cavitation flow field in the tip region of the axial flow pump, and the primary cavitation appeared in the tip region of the axial flow pump. Under different cavitation numbers, the cavitation area coefficient of the suction surface increases from hub to flange, reaches the maximum in the middle of the blade, and then decreases rapidly. Tao [14] found that the balance hole has a significant effect on balancing axial force of the centrifugal pump. Cheng et al. [15] found that the cavitation performance of the high-speed pump can be improved when the section angle of the meridional blade increases within a certain range. Zhou et al. [16] simulations show that the sine wave structure in HCD generates an important flow structure-counter-vortex pair (CVPs), which increases the momentum exchange between the near-wall region and the flow core. Gai et al. [17] found that the collapse of cavitation bubbles in the crevice at the particle surface is structure sensitive.
During the operation of the pump, the impeller drives the liquid to rotate. The suction port of the centrifugal pump impeller is low pressure and the back is high pressure. Due to the different pressures before and after the impeller, a pressure difference will be formed between the front wheel cover and the rear wheel disc of the impeller, which acts on the pressure on the front cover and the rear cover that cannot be balanced with each other which is the main reason for the imbalance of the axial force. As the main means to balance the axial force of the pump, the balance hole is widely used in the impeller pump, but the current research on the balance hole of the semiopen impeller centrifugal pump is very limited. This article takes the semiopen impeller centrifugal pump as the research object. By changing the position of the balance hole, the influence of the circumferential position of the balance hole on the cavitation performance of the semiopen impeller centrifugal pump was studied. On this basis, the variation law of the axial force of the centrifugal pump with different balance hole positions was studied.
2. Geometric Model and Design Scheme
2.1. Basic Parameters of Semiopen Centrifugal Pump
The semiopen impeller centrifugal pump is taken as the research object, and its main parameters are shown in Table 1. The impeller blades are 6 long and 6 short backward curved cylindrical blades, and the tip clearance is uniform and equivalent clearance. The fluid domain of the semiopen centrifugal impeller is modeled by 3D software, as shown in Figure 1.

2.2. Scheme Design
The diameter and the number of the balance hole are determined to be 6 mm and 6, considering the impeller structure and the effect on the external performance of the pump. Calculate the diameter of the balance hole according to the following empirical formula. According to the calculation result, the diameter of the balance hole is 6 mm:
Among them, dk is the diameter of the balance hole and s is the cross-sectional area of the seal ring gap.
The location of the balance holes is optimized in this study. On this basis, the influence of the circumferential position of the balance hole on the cavitation performance of the semiopen impeller centrifugal pump is studied. As shown in Figure 2, the angle θ of the balance hole and the leading edge of the blade in the direction of rotation is defined as the circumferential position angle of the balance hole. In this study, four different circumferential positions were designed to arrange the balance holes. The layout position is shown in Figure 2. The angle θ between the origin and the center of one of the balance holes OA and the leading edge of the blade in the rotation direction is defined as the circumferential position angle of the balance holes. With 0° as the starting position, the design balance holes are rotated anticlockwise from the starting position, and the other five balance holes are uniformly distributed in array mode. The opening positions of balance holes in the four schemes are shown in Table 2, where D represents the distance of OA.

3. Numerical Calculation Method
3.1. Computational Domain Meshing
The fluid domain of the semiopen impeller centrifugal pump model is mainly composed of six parts: inlet section, outlet section, tip clearance, impeller, volute, and back cavity. ICEM software is used to mesh the fluid domain. The calculation domain of the tip clearance and the semiopen impeller adopts a hexahedral structured grid, and the volute, back cavity, inlet, and outlet sections adopt unstructured grids. The grid in the tip clearance area is locally encrypted. The number of encryption layers is 20. The height of the first layer near the wall of the tip clearance calculation domain is 6 × 10−5m, and the y + value is between 1 and 10. Finally, the total number of grids generated in the tip clearance calculation domain is 1.69 million. The grid diagram of the key water-passing parts is shown in Figure 3.

Considering the influence of grid density on the calculation results, four kinds of grids of 5.4 million, 6.0 million, 6.8 million, and 8.6 million are used to verify the independence. Figure 4 shows that, with the increase of grid number, the head of the model pump shows a downward trend, but the change rate is gradually decreasing. When the grid number reaches about 6.8 million, the head tends to be stable. Therefore, in order to ensure the accuracy of numerical calculation and choose the grid scheme with lower calculation cost, the number of grids in the calculation domain is about 6.8 million.

3.2. Turbulence Model
The turbulence model used in this paper is the SST k-ω turbulence model. Compared with the general k-ε model, the SST k-ω model has a higher simulation accuracy for the free flow near the wall and can capture the vorticity distribution inside the tip clearance more effectively. It is more widely used in the research on the clearance problem.
The SST k-ω model has a similar form to the standard k-ω model, including k equation and ω equation [18], and its basic form is as follows:
In the formula, ω is defined as a specific dissipation rate. Pk represents the generation term of turbulent pulsation kinetic energy k. Pω is the generation term of turbulent pulsation frequency. Γk and Γω represent the effective diffusion coefficients of k and ω, respectively. Yk and Yω represent the dissipation term of k and ω, respectively. Dω represents the orthogonal diffusion term.
3.3. Cavitation Model
Based on the homogeneous flow hypothesis, the cavitation flow is generally regarded as a unidirectional flow of homogeneous fluid, which is used to calculate the mass transfer process between gas and liquid phases. In order to describe the process of cavitation bubble development and collapse more accurately, the cavitation model proposed by Rayleigh–Plesset is adopted. The following formula is an idealized basic equation of cavitation bubble dynamics without empirical coefficients [18]:where .
In the formula, R is the cavity radius; ρ is the liquid density; μ is the kinematic viscosity of the liquid; PB and P∞ represent the cavity internal pressure and hydrostatic pressure, respectively; σ is the surface tension coefficient.
3.4. Boundary Conditions
This study uses related software to perform numerical simulation calculations on the entire fluid domain. Set the boundary conditions of the model pump. The inlet is set as the total pressure inlet, the outlet is set as the mass flow outlet, the wall conditions are selected as the solid wall without slippage, and the fluid medium is 25°C (constant temperature) clear water. The saturated vapor pressure is 3169 Pa, and the system reference pressure is set to 0. The result based on the noncavitation steady calculation is set as the initial value of the cavitation calculation, and the total inlet pressure of the calculation model is reduced gradually to cause cavitation in the pump, by observing that the outlet pressure tends to be stable or the convergence residual value is less than 10−5, it is determined that the solution has converged.
3.5. Validation of Numerical Simulation
The initial model impeller with the balance hole circumferential position angle θ = 30° is selected as the test object to carry out the external characteristic test of the centrifugal pump. The test system is a visual test bench for centrifugal pumps, as shown in Figure 5. In the test, the test data of the flow rate of 0–9 m3/h is recorded, and the corresponding pump head and efficiency inlet and outlet pressure and flow rate are calculated through actually measured pump. The cavitation performance test system consists of a visual test bench. First, carry out the external characteristic test of the pump on the test bench. After the pump is started normally, keep the inlet valve unchanged, adjust the outlet valve to change the flow, and record the performance parameters of the pump at different flow rates. The cavitation performance test system consists of a visual test bed and a high-speed camera system. First, fully open the inlet valve, and adjust the outlet valve to control the flow at 6 m3/h. Then, gradually close the inlet valve under the condition of constant flow, observe the change in vacuum, and stop when the calculated equivalence is reached. Recording data such as flow rate, inlet and outlet pressure, and torque, and measure 15 times in this cycle until the vacuum level no longer rises. The casing of the model pump of the fully transparent test bed is made of plexiglass. The high-speed camera shooting system consists of a high-speed digital camera, a light source, and a computer. The high-speed camera adopts the Phantom v1210 high-speed camera launched by Vision Research, which can shoot up to 1000000 frames per minute. The sampling frequency set in this experiment is 1000 Hz.

(a)

(b)
Select the flow rate from 4 m3/h to 9 m3/h for numerical simulation. Figure 6 shows the external characteristic curve obtained by numerical calculation of different schemes of a semiopen impeller centrifugal pump at a speed of 1500 rpm, as well as the flow-head curve and flow-efficiency curve measured by the test. From the change trend of the external characteristic curve in the figure, the test results of head and efficiency are in good agreement with the simulation results, especially when the design flow rate is around 6 m3/h. The difference between the head test value and the head calculated value is within 1%, and the difference between the efficiency test value and the efficiency calculated value is within 5%, indicating that the current calculation result has a certain degree of accuracy and reliability. The selected numerical method can be used in this research. The formula for calculating the head of a centrifugal pump iswhere H is head, m; Z is position head, m; p is pressure, Pa; is velocity, m/s; ρ is fluid density, kg/m3, subscript 1 indicates pump inlet, and subscript 2 indicates pump outlet (the reference plane is selected on the horizontal plane where the centerline of the pump shaft is located).

The efficiency calculation formula of the pump iswhere QV is flow rate, H is head, m, P is input power, kw, and ρ is fluid density.
According to the curve changes of the four schemes, it can be found that when θ = 10°, the head and efficiency of the centrifugal pump are the highest. As the value of θ increases, the head and efficiency decrease. When θ = 30°, the head and efficiency are the lowest among the four schemes. When the value of θ is further increased to 45°, the head and efficiency begin to increase. It shows that the closer the opening position of the balance hole is to the working surface of the long blade, the better its head and efficiency is. But when the position is far closer to the back of the next blade, the head and efficiency will be improved. Therefore, the farther the balance hole is placed from the working surface of the long blade, the changes in efficiency are uncertain. Next, we will analyze the cavitation performance of centrifugal pumps equipped with balance holes in different circumferential positions under design conditions ( = 6 m3/h).
4. Numerical Calculation Results and Analysis
4.1. Static Pressure Distribution of the Impeller Back Shroud
It can be seen from the static pressure distribution on the axial surface of the impeller back shroud that when the balance holes are installed at different circumferential positions, the pressure distribution on the impeller will be affected to a certain extent, and the low-pressure distribution at the impeller inlet has obvious changes (see Figure 7). Among them, when θ = 30°, the area of the low-pressure area near the impeller inlet is the largest. Taking θ = 30° scheme as the boundary, when the value of θ is reduced or increased, the low-pressure area is significantly reduced, and the pressure distribution at the inlet also becomes more uniform, while the cavitation conditions are improved. Comparing the solutions of θ = 15° and θ = 45°, it can be found that although the low-pressure area in the flow passage near the impeller inlet is obviously reduced, the degree of reduction is different. When θ = 45°, the low-pressure area occupies more flow passages, and the local low-pressure areas are concentrated near the balance hole. Because when the fluid enters the suddenly reduced balance hole passage from the large space passage of the back cavity, the fluid will form a jet at the balance orifice position to generate a vortex, causing a local low-pressure area. The existence of such a low-pressure area can easily cause jet cavitation near the balance hole. This shows that selecting an appropriate circumferential position of the balance hole can provide a more stable inflow condition for the impeller, thereby improving the cavitation performance.

4.2. Cavitation Performance Analysis
The cavitation performance of the pump is usually described by the relationship curve between the effective cavitation margin and the head of the pump. This paper draws the cavitation margin-head curve and the cavitation margin-flow curve through research, as shown in Figure 8. Among them, when the head drops by 3%, the corresponding effective cavitation margin is defined as the critical cavitation margin, and the effective cavitation margin is calculated aswhere represents the saturated vapor pressure of the fluid medium at 25°C and Vin represents the inlet velocity.

Ignoring the second-order term and surface tension term, the equation can be simplified to
Total mass transfer rate between phases:where
For the collapse and condensation process of bubbles,
In summary, the final expression of the cavitation model iswhere Re, Rc are steam generation rate and condensation rate, respectively. αruc is nucleation site volume fraction. Fvap is the empirical coefficient of the steaming process, taking 50. Fcond is the empirical coefficient of the coagulation process, taking 0.01. The cavitation radius R is 0.001 mm.
The change trend of the cavitation characteristic curve under the four schemes is basically the same. In the early stage, the head change is less than 0.1%; as the NPSHa decreases, the head begins to drop sharply until it drops to 0 m. The whole process can be understood as four stages: the noncavitation stage where the head is no change, the initial stage of cavitation (a-b) where the head remains stable within a certain range, the critical cavitation stage (b-c) where the head drops by 3%, and the complete cavitation stage (after point c) where the head drops steeply. The dotted line in the figure represents the experimental results. Compared with the numerical calculation results of the same model, it is found that the head variation law is consistent, and the maximum error between the data is about 1.5%, indicating that the numerical calculation method used in this article is feasible. As the effective cavitation margin decreases, the θ = 30° scheme first enters the critical cavitation stage, then the θ = 10° scheme, and finally the θ = 15° scheme. It can be seen that, compared with other schemes, the pump head and efficiency are the highest but the anticavitation ability is not the best when θ = 10°. In contrast, when θ = 15°, the centrifugal pump has the best cavitation performance.
Figure 9 shows the distribution of cavitation bubble in the tip clearance and the impeller under cavitation conditions. Choose the situation when the NPSHa = 2.71 m, 0.66 m, and 0.46 m for analysis. It can be seen from the figure that, at the initial stage of cavitation when NPSHa = 2.71 m, there is no obvious cavitation area in the tip clearance and the impeller flow passage, but the head begins to decrease. According to the results in Figure 8, when NPSHa = 0.66 m, the θ = 30° scheme has reached the critical cavitation state, and the other three schemes have not yet entered the critical cavitation state. However, it can be seen from the cavitation bubble distribution diagram that there is an obvious cavitation area in the impeller flow passage under the four solutions at this time, and the bubble is distributed asymmetrically on the back of the blade. There is also an obvious leakage vortex cavitation in the tip clearance, and its distribution position is consistent with that in the impeller flow passage, which shows that the cavitation of the tip clearance leakage vortex and the cavitation in the impeller flow passage affect each other. The rotation of the impeller causes dynamic and static interference between the impeller and the volute, which makes the flow state in the centrifugal pump vary greatly along the circumferential position. In particular, the complex structure of the volute tongue makes the fluid flow near the tongue of the volute the most affected by the volute. Therefore, the cavitation in the flow channel is distributed asymmetrically. There are obviously fewer bubbles in the flow channel near the distance from the tongue. On the contrary, the symmetrical flow channel has a larger number of bubbles. Comparing the bubble distribution under different schemes, the cavitation area occupies the most impeller flow passage when θ = 30°; the cavitation area of the impeller is the largest but the tip clearance cavitation is better than other schemes when θ = 45°. In contrast, the cavitation performance is better when θ = 15° and θ = 10°. This is consistent with the changing law of the cavitation characteristic curve shown in Figure 8.

In order to analyze the cavitation in the balance hole more intuitively, the isosurface with the cavitation bubble volume fraction VOF = 10% is taken to observe the cavitation bubble distribution in the balance hole when NPSHa = 0.66 m. It can be seen from Figure 10 that, under the same effective cavitation margin, the cavitation bubble distribution in the balance hole flow passage in different schemes is different. As the value of θ increases, the cavitation bubble in the balance hole gradually increases, and the flow in the balance hole is gradually blocked, which weakens the pressure relief ability of the balance hole. At the same time, it can be found that, as the value of θ increases, the cavitation near the impeller inlet extends toward the impeller outlet on the one hand and also develops near the balance hole on the other hand, indicating that the position of the balance hole affects the development of cavitation in the centrifugal pump flow passage.

4.3. Axial Flow of Centrifugal Pump Flow Passage
Figure 11 shows the streamline distribution at different axial sections of the centrifugal pump tip clearance, impeller, and back cavity flow passage when the cavitation margin is 0.66 m. The schematic diagram of the axial section position is shown in Figure 12. It can be seen from Figure 11(d) that, at section D-D, the streamline distribution at the impeller inlet position is different under different schemes. When θ = 10°, there is no obvious vortex area at the impeller inlet; as the value of θ increases, vortices appear gradually. When θ = 30°, the vortex area at the entrance of the D-D section is the most serious. When θ = 45°, affected by the balance hole, the flow field flow is disturbed, the vortex area decreases, and the position gradually moves upward. In the remaining three sections, the flow conditions in the same axial flow channel under different schemes are basically the same, and the flow field at the impeller inlet is not greatly interfered by the balance hole. From the flow conditions in the balance hole under different schemes, it can be found that the larger the value of θ, the higher the flow velocity in the balance hole, which will lead to a local low-pressure area at the balance hole position, which is prone to jet cavitation. Therefore, when θ = 30° and 45°, cavitation is more likely to occur near the balance hole, which is consistent with the gas volume distribution shown in Figure 9.


4.4. Analysis of the Axial Force of the Impeller
In solving the problem of axial force imbalance, opening a balance hole on the back cover has always been an effective way to balance the axial force, but for the semiopen impeller centrifugal pump without front cover, the axial force generation position is different from that of the closed impeller. In order to study the influence of the balance hole on the axial force characteristics of the pump rotor, it is necessary to analyze the force of the pump. The force distribution of the model pump is shown in Figure 13, the axial force of the semiopen impeller is composed of the forces F1 acting on the rear cover and F2 acting on the front side of the rear cover, the force F3 exerted on the impeller by the fluid in the tip clearance and formed by the pressure difference between the working surface and the back of the blade, and the dynamic reaction force F4. The calculation formula of the cover force (T1) is as follows [18]:

Among them, R2 is the impeller outlet radius, Rm is the impeller inlet radius, and Rh is the shaft radius.
The above formula is the calculation from theoretical experience. In order to verify the reliability of the numerical calculation, the numerical calculation results of the cover force T1 under different schemes are compared with the formula calculation results (as shown in Table 3). According to the data in Table 3, it is found that the error between the numerical calculation result and the theoretical formula is within 5%. It shows that the axial force obtained by numerical calculation is reliable. Table 4 shows the changes of the rotor axial force under different schemes. It can be seen from the table that, compared with the model without the balance hole, installing the balance hole can effectively balance the axial force, and the axial force is reduced by 60%. At the same time, changing the circumferential position of the balance hole also has a certain effect on the balance of the axial force, but the degree of change is small. In contrast, the θ = 45° scheme is slightly better than other schemes in the effect of balancing the axial force. Figure 14 shows the change curve of the impeller axial force and the amount of balance hole leakage in cavitation. The axial force data in the figure is recorded from the beginning of cavitation (NPSHa = 2.71 m). It can be seen that, with the NPSHa decreasing, the axial force on the impeller shows a downward trend. At the same time, the amount of balance hole leakage is also decreasing. The reason is, with the development of cavitation, the bubbles at the inlet of the impeller increase, which block the flow channel and weakens the function of the impeller, and both the back cavity pressure and the inlet pressure of the impeller decrease. This leads to a reduction in the static head of the impeller. According to the theory of the axial force of the cover, the cover force (T1) is reduced (see formula (15)), so the total axial force of the impeller is reduced. On the other hand, when cavitation occurs, bubbles will gradually block the flow in the balance hole (see Figure 10), resulting in a reduction in the amount of the balance hole leakage, and its pressure relief effect is gradually weakened. When the pressure in the back cavity drops to a certain value, the pressure will no longer decrease, while the impeller inlet pressure is still in a downward trend. Therefore, after the cavitation develops to a certain extent, the change amplitude in the axial force of the impeller will begin to shrink. It can also be found from Figure 14 that, with the development of cavitation, the reduction amplitude of axial force in θ = 45° scheme first trends to be stable, which shows that the back cavity pressure is first reduced to the lowest value under this scheme. It is further explained that when θ = 45°, the cavitation performance of the impeller is worse than other schemes. Under cavitation conditions, the larger θ is, the greater the pressure relief capacity of the balance hole will be affected by cavitation.

5. Conclusion
(1)In the process of changing the circumferential position of the balance hole, the efficiency and head of the centrifugal pump change regularly. With the decrease of θ, the efficiency and head gradually increase, but the change is not great, so the circumferential position of the balance hole has little effect on the external characteristics of the semiopen impeller centrifugal pump.(2)As the circumferential position of the balance hole changes, the anticavitation ability of the semiopen impeller centrifugal pump has changed. Compared with other schemes, the anticavitation ability of the centrifugal pump is weak when θ = 30°. With the decrease of θ, the internal cavitation of the centrifugal pump gradually improves. When θ = 15°, the anticavitation ability of the centrifugal pump is the best. The larger the θ of the balance hole is, the more easily the jet cavitation occurs near the balance hole.(3)When changing the circumferential position of the balance hole, with the increase of θ, the axial force on the impeller decreases slightly, indicating that changing the circumferential position of the balance hole can improve the impeller axial force condition. However, the development of cavitation has a greater impact on the axial force of the impeller. Under cavitation conditions, the farther the balance hole is located in the circumferential direction (i.e., the greater θ), the more limited the ability to balance the axial force.(4)For low specific speed semiopen impeller centrifugal pumps, considering the influence of the balance hole on its external characteristics, cavitation performance, and axial force, the circumferential position of the balance hole closed to the blade working surface (the circumferential position angle is 15°) is the best location.Nomenclature
: | Design flow rate (m3/h) |
H: | Design head (m) |
n: | Rotational speed (rpm) |
ns: | Specific speed |
θ: | Circumferential position angle of the balance hole (°) |
Φ: | Diameter of balance hole (mm) |
ρ: | Fluid density. |
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the Chinese National Natural Science Foundation (no. 51469013).