Abstract
Finger seals are a new type of seal with good sealing performance and long service life. The noncontacting feature relies on the gas film force. However, when the seal works in an unsuitable environment or its design parameters are not reasonable, the lifting pad may not be able to generate sufficient air film force. This causes contact between the fingers and the rotor, resulting in a reduced service life of the seal. In view of this situation, this paper proposes a method that can quickly determine whether there is enough gas film force to lift the sealing finger at the design stage. The aeroelastic coupling characteristics of the noncontacting finger seal are studied in conditions where contact exists between the fingers and the rotor. The influences of various environmental and key structural parameters on the number of contact fingers, leakage, bearing force, and friction moment are studied. The results show that the pressure difference, eccentricity, and key design parameters have important effects on the number of contact fingers. The effect of rotation speed is relatively small. This paper provides a time-efficient tool for the design of noncontacting finger seals, which can quickly predict the performance of the sealing system.
1. Introduction
Gas turbines are widely used for high-power and high-efficiency applications. To ensure the proper operation of gas turbines, stringent requirements are made on their internal components. As essential components, advanced and well-designed seal equipment can effectively reduce the leakage of gas turbines. In recent years, many researchers have widely studied a new adaptive seal, namely, finger seal [1], due to its low cost and good sealing performance. A spring-mass-damper equivalent dynamic model is used to obtain the seal performance. Braun et al. [2] studied the effects of the fluid stiffness, fluid damping, and finger stiffness on the phase shift and displacement transmissibility. Chen et al. [3] studied the displacement responses, leakage clearances, and contact pressures of finger elements through a distributed mass model. The calculation results agree well with the experimental data. Experimental tests [4] show that the sealing performance of the finger seal is acceptable for engines, but the seal exhibits some hysteresis, and wear occurs at the beginning of the test. Guoqing et al. [5] experimentally studied the leakage and wear characteristics of finger seals in hot/cold states. One interesting finding is that double-laminate finger seals can achieve lower leakage in cold conditions, while triple-laminate finger seals perform better in hot conditions. A pressure-balanced finger seal [6] is proposed, which can reduce the hysteresis. To reduce wear and obtain a long-life application, a noncontacting finger seal was presented by Proctor MP and Steinetz BM [7], which has hydrodynamic lifting capabilities and a noncontacting nature.
Dynamic models have been proposed to study noncontacting finger seals. Braun et al. [8] studied the thermofluid and dynamic behavior of a two-laminate finger seal. Marie [9] provided theoretical supports for finger seal design by parametrically studying a noncontacting finger seal assembly. Du et al. [10] developed a semianalytical model that can save computation time. The dynamic properties of noncontacting finger seals were found through the model. In order to have an excellent sealing performance, the ratio of friction force to finger stiffness and the rotor excitation amplitude should be as small as possible; the finger natural frequency should be as large as possible [11]. The fluid-solid interaction method is another commonly used way to study noncontacting finger seals. Yue et al. [12] used commercial software to analyze the leakage, the pressure distribution under lifting pads, the fluid flow through the finger seal, and the displacement of lifting pads. Braun et al. [13] developed solid models to study the deformations, stress of fingers, and lifting capabilities for a better parametric design of a finger seal. Zhang et al.[14] used a two-way fluid-structure interaction method to study the deformations of fingers and the flow field through seal assembly. The effects of different groove structure lifting pads on leakage flows and the lifting force on the finger were also studied. Studies show the possibility of fingers contacting the rotor [15]. Jia et al. [16] proposed a three-layer finger seal with grooved lifting pads. A numerical method was used to investigate the effects of the grooved structures on lifting and leakage performance. The grooved structure can significantly increase the lifting capacity of the finger seal and does not increase leakage too much. Temis et al. [17] proposed a two-way fluid-solid interaction method to calculate the seal characteristics of a noncontacting finger seal. The results of the method are very close to those obtained by finite element software. Zhao et al. [18] considered the side leakage flow into the seepage flow through a porous medium to analyze the total leakage of finger seals. Leakage is increased by 3 to 6% when considering side leakage. Fleming [19] presented an approximate solution to account for fluid inertia in gas flow through rectangular seal pads. The results show that fluid inertia has beneficial influences on noncontacting finger seals. Proctor et al. [20] tested a noncontacting finger seal with a grooved rotor, and the noncontacting operation was achieved. Smith and Braun [21] tested four-finger seal configurations whose lifting properties and structure compliance were investigated through parametric studies.
Above all, we can see that the leakage performance, lifting ability, and effects of the structure and environmental parameters of noncontacting finger seals are widely studied. However, the noncontacting feature requires specific operating conditions to be achieved because the lifting forces generated by gas films are limited. When the rotor is at high eccentricity or the seal is operating in an unsuitable working condition, the lifting forces may not be enough to keep the fingers noncontact with the rotor. It results in contact during an operating rotor system. The contact is inconsistent with the original purpose and should be avoided at the design stage. Thus, how to determine if there is contact and the effect of this contact should be studied.
In this paper, a method to determine whether there is contact between the fingers and the rotor is proposed. The effects of contact on the static characteristics of the noncontacting finger seal are studied. First, a basic analysis unit is divided to calculate the finger deformation. The formula of gas films between the lifting pads and the rotor is deduced by considering the deflection, rotation angle, and torsion of the fingers. The film forces and pressure distributions of the lifting pads are obtained through gas lubrication theory. Then, the proposed method is used to determine whether the fingers are in contact with the rotor. Based on the judgment results, the static characteristics of each finger are calculated. For fingers not in contact with the rotor, an aeroelastic model is used. For fingers in contact with the rotor, a Coulomb friction model is used. Finally, the static characteristics of the entire seal can be obtained by summing the previous results of each finger. The effects of structural parameters, working conditions, and the number of contact fingers are parametrically studied. Compared with the three-dimensional finite element method, the method presented in this paper can judge whether the designed noncontacting finger seal can achieve the noncontacting feature under specified working conditions at the design stage. Static characteristics under various operating conditions when the contact exists can be quickly calculated, thus providing a time-efficient approach for the design of noncontacting finger seals.
2. Method
2.1. Basic Analysis Units
As shown in Figure 1, the noncontacting finger seal studied in this paper consists of one high-pressure laminate, one low-pressure laminate, and one backplate with a pressure-balanced design. Uniform cuts of the laminates form flexible elements, which are called fingers. The gaps of the laminate are covered by the fingers of the other laminate. Each low-pressure laminate finger has a lifting pad to generate film force by the gas film between the lifting pad and the rotor. The film forces lift the fingers, resulting in no contact between the fingers and the rotor. However, film forces are limited, and it may not be possible to keep the finger out of contact with the rotor in some cases. For example, the fingers need a much larger film force to generate enough deformation when the rotor is under large eccentricity conditions. To better investigate the static performance of the noncontacting finger seal, the finger model and the fluid model need to be built first.
The structure of both laminates is periodically symmetrical, and the basic unit is divided according to the structural characteristic. As shown in Figure 2, each basic unit consists of two split high-pressure fingers and one low-pressure finger. During the calculation, it is assumed that there is no relative displacement between the fingers of the same basic unit. The backplate provides support for the fingers and greatly reduces their axial deformation. Thus, the axial deformations of fingers are neglected during the analysis.
2.2. Forces Calculation under Fluid Lubrication Condition
Based on the assumptions above, the deformation of the basic units can be obtained by the low-pressure finger. Figure 3 shows the schematic of a low-pressure finger under fluid lubrication. The finger is subject to three forces in the y-z plane. Fk represents the force acting on the low-pressure finger by the high-pressure finger. The direction of Fk is shown in Figure 3. The magnitude of Fk iswhere k is the stiffness of the high-pressure finger and u is the deformation of the high-pressure finger.
The stiffness of the finger can be obtained by applying a radial unit force at the end of the high-pressure finger and calculating the corresponding deformation. The deformation is calculated in the same way as the deformation of the low-pressure finger is obtained in equation (9).
Ff represents the friction force between the low-pressure finger and the backplate. It is greatly reduced by the pressure-balanced dam and is much smaller than Fk and the film force (Fh) which means it has little effect on the deformation. Ff always hinders the finger deformation. Thus, the direction of Ff is considered downward during the calculation. The magnitudes of Ff are calculated by the following equation:where μlb is the friction coefficient between the low-pressure fingers and backplate, Plb is the pressure that pushes the finger laminate against the backplate, and Alb is the contact area between the low-pressure laminate and the backplate.
Fh represents the film force generated by the film between the lifting pad and the rotor surface. The film consists of four parts: the rotor eccentricity, finger deformations, wedge shape of the lifting pad, and initial installation gap. As shown in Figure 4, the film thickness of a basic unit can be deduced:where uif is the finger radial deformation, e is the rotor eccentricity, φ is the attitude angle, hrif is the film change due to rotation angle of fingers, is the wedge under lifting pads, T is the film change due to the torsion of fingers, hinit is the initial installation gap, αr is the finger repeat angle, and nf is the number of fingers.
(a)
(b)
The wedge under each lifting pad is shown in Figure 5. The dashed line represents the original shape of the lifting pad, and and are the minimum and maximum values of the wedge, respectively.
The film force and pressure distribution under lifting pads are determined by the nonlinear compressible Reynolds equation:
The static compressible dimensionless form is as follows:where , , , and .
Equation (6) is linearized by using the Newton–Raphson method [22]:where , and ξ is the number of iterations.
The finite difference method is used to discretize the linearized equation. The iterative formula can be obtained:where are indices for the grid in the circumferential and radial directions, respectively, and are coefficients shown in the appendix.
The boundary condition of the calculation area is defined in this way: among the four sides of the lifting pad, the pressure of the high-pressure side is related to the axial thickness of the laminate and the axial length of the lifting pad. The pressures of the other three sides are equal to the downstream pressure. The boundary conditions can be obtained as follows:where b is the axial width of the finger laminate, L is the axial length of the lifting pad, and n and m are the numbers of grids.
Then, the deformation of the low-pressure finger is calculated by using Mohr’s integrals [23], and the torsional deformation is obtained using elasticity theory:where uh, , , and are finger deformation in the radial direction, the axial force, bending moment, and shear force under the fluid lubrication condition, respectively; , , and are the internal forces due to the unit load, and A and I are the area and moment inertia of a cross section; E and G are Young’s modulus and shear modulus of elasticity, for rectangular cross section η = 1.2, αul is the curvature of the fingers, τ is the twist angle, MT is the torque, bs and hs are the cross-sectional parameters, and βT is the factor determined by bs and hs.
2.3. Forces Calculation under Contact Condition
When the finger is in contact with the rotor, the film force is replaced by the contact force applied by the rotor. Due to the small circumferential length of the lifting pad, the surface of the rotor can be considered a plane with respect to the lifting pad during the contact judgment. Since the axial deformation is neglected, the torsion does not affect the judgment of the contact position, which means that T in equation (3) is taken as a unit matrix in the process of judgment. Thus, it can be seen from Figure 5 that the rotor first contacts the finger at position 1 or 2.
When the rotation angle is small, it has little effect on the film thickness, which means that the first contact is at position 1. From (3), the finger deformation under this condition can be obtained:where θ1 is the angle of position 1.
Similarly, the first contact position is at position 2 when the rotation angle is large. The finger deformation iswhere θ2 is the angle of position 2.
The process for determining the contact position is shown in Figure 6. Take the direction of the green arrow in Figure 4 as the positive direction. For each finger, the possible deformation (u1 or u2) when contact occurs can be calculated by equations (10) and (11). Due to the pressure difference, the finger deformations should be positive. If u1 and u2 are negative, it means that the finger needs to produce a negative displacement to contact the rotor. This is inconsistent with the above conclusion. In other words, there is no contact between the finger and the rotor. If u1 and u2 have positive and negative values, the contact is at the position corresponding to the positive value. If both u1 and u2 are positive, the contact is assumed to be at position 1, and then the film thickness is calculated. If contact does occur at position 1, the film thickness at position 1 (h1) should be less than that at position 2 (h2). Then, contact is considered to occur at position 1; otherwise, contact occurs at position 2. The deformation of the contact finger (uc) is obtained.
Using uc, Fk is determined by equation (1). When there is contact between the rotor and the finger, the force acting on the finger is denoted by Fc. Fc satisfies the condition that the deformation of the finger under the forces of Fk, Ff, and Fc is equal to uc. Figure 7 shows the process of finding Fc. After entering the initial values, the finger deformation when the finger contacts the rotor (ucnew) is calculated bywhere Nc, Mc, and Qc are the axial force, bending moment, and shear force of the finger under the contact condition, respectively.
Compare uc and ucnew. If the absolute value of the relative error is greater than 10−6, update the value of Fc according to the formula shown in Figure 7. Calculate ucnew again with new Fc (Fcnew). This process is repeated until the relative error is less than 10−6. Then, the contact force Fc is obtained.
To judge whether a finger contacts the rotor, the relationship between the film force and the film thickness should be obtained first. The film thickness ratio describes the relationship between the minimum film thickness and surface roughness [24]:where γ is the film thickness ratio and Rar and Rap are the surface roughness of the rotor and the lifting pad, respectively.
When the lifting pad is not in fluid lubrication, a gas film between the lifting pad and the rotor cannot be formed, which results in contact. This paper does not consider the mixed lubrication state. The value of γ is taken as 3 to determine the minimum film thickness. As shown in Figure 8(a), the film force is larger than the contact force to the left of the dotted line, which means that the film can generate enough film force to lift this finger. As shown in Figure 8(b), the film between the lifting pad and the rotor is divergent due to the large rotation angle of the fingers. The film force increases with increasing film thickness. The contact force is always larger than the film force over the entire range of film thicknesses. The film force is not great enough to lift this finger. Thus, the finger is judged to contact the rotor.
(a)
(b)
The comparison between the previous simulation results and the results calculated using this paper’s data and methods is shown in Figure 9. The figure shows the force acting on the rotor in the z-direction (shown in Figure 4) at different attitude angles of the rotor. It can be seen that the variation trend of the forces is in good agreement. When contact exists between the lifting pads and the rotor, the force acting on the rotor increases significantly.
2.4. Static Performance
This paper calculates four static characteristics of noncontacting finger seals: leakage, lift force, friction moment, and bearing force. Leakage under fluid lubrication is defined as follows:
When the average axial flow velocity is larger than the local sound velocity, the gap between the lifting pad and rotor is considered chocked. The local sound velocity is used to simplify and replace the average axial velocity of gas flow. The leakage when chocked flow occurs is defined aswhere ρ is the density of the upstream gas, V is the average axial flow velocity, and As is the sectional film area.
After the finger contacts the rotor, its deformation is complex to predict. On the other hand, the thickness under the contacting finger is much smaller than that under the noncontacting finger. Thus, the leakage of the fingers under contact conditions is ignored.
The lift force of each gas film is defined aswhere Ap is the area of the lifting pads and pL is the downstream pressure.
For the fingers under fluid lubrication, the bearing force and the friction moment acting on the rotor by each finger are defined as follows:
For the fingers under the contact condition,where μfr is the friction coefficient between the fingers and rotor.
Thus, the bearing force and friction moment acting on the rotor arewhere nh and nc are the number of fingers under the fluid lubrication and contact conditions, respectively.
The aeroelastic coupling calculation flow is shown in Figure 10. The convergence conditions are as follows:where λ is the number of iterations.
3. Results
The parameters of the noncontacting finger seal in this paper are given in Table 1.
Figure 11 shows the effects of eccentricity at different pressure differences. As shown in Figure 11(a), the leakage decreases first and then increases at 50 kPa pressure difference. As the eccentricity increases, some of the clearances are reduced, resulting in a decrease in leakage. The clearances under the other lifting pads are increased, resulting in an increase in leakage. Figure 11(d) shows the leakage of each finger at different eccentricity and pressure differences (fingers 1 to 80 are numbered according to the circumferential direction). The reduced leakages through fingers 1 to 20 and 61 to 80 are less than the increased leakages through fingers 40 to 60 at 0.16 eccentricity and 50 kPa pressure difference. Thus, the total leakage is reduced. The leakages through clearances are reduced because of the contact (fingers 1 to 14 and 70 to 80) and the eccentricity (fingers 15 to 20 and 61 to 69) at 0.64 eccentricity and 150 kPa pressure difference. However, the total leakage is still greater than that at 0 eccentricity due to the increased film thickness caused by the increased eccentricity. When the pressure difference is equal to 50 kPa, a slight decrease in leakage with increasing eccentricity appears in the eccentricity range of 0.08 to 0.24. When the pressure difference is equal to 150 kPa, this change occurs in the eccentricity range of 0.32 to 0.56. This is due to the larger finger deformations and the greater clearance between the lifting pads and the rotor under the high-pressure difference.
(a)
(b)
(c)
(d)
(e)
Figures 11(b) and 11(c) show the friction moment and bearing force for various eccentricities at different pressure differences. They increase with increasing eccentricity and are greater at 50 kPa pressure difference than at 150 kPa pressure difference. This is because the film forces are smaller at low-pressure differences, resulting in more fingers in contact with the rotor. As shown in Figure 11(e), when there is no contact between the rotor and the fingers, the pressure distributions are similar for most gas films, making the bearing force, that is, the combined force acting on the rotor, small. When contact occurs between the fingers and the rotor, the contacting force is greater than the film force, resulting in a significant increase in the bearing force. The higher the eccentricity is, the higher the number of sealing fingers in contact is, and the higher the bearing force and friction moment are.
Figure 12 shows the effects of rotation speed at 100 kPa pressure difference and different eccentricities. As shown in Figure 12(a), the leakage decreases with increasing rotation speed. The rotation angles of the fingers are large at 100 kPa pressure difference, which means that the gas films are divergent. The higher the rotation speed is, the lower the film force is and the smaller the clearance is, resulting in decreased leakage. When the rotor is at high eccentricity, the decrease in the film force leads to contact between the finger and the rotor. From Figure 12(b), the friction moment increases with increasing rotation speed. At 0.2 eccentricity, the increase in friction moment is continuous because there is no contact between fingers and the rotor. Negative friction moments occur at a low rotation speed, which means that the direction of the friction moment is the same as the rotation speed. For a divergent gas film at 100 kPa pressure difference, the friction moment generated by the pressure difference is in the same direction as the rotating speed (expressed as the negative direction in this paper). When the rotating speed is low, the friction moment is mainly generated by the pressure difference. With increasing rotation speed, the friction moment generated by rotation increases, and the friction moment becomes positive (opposite to the direction of rotation). At 0.6 eccentricity, the friction moment consists of two types: the friction moment caused by contact and the friction moment caused by the gas film. The higher the rotation speed is, the bigger the number of contact fingers is and the greater the Coulomb friction of all contact fingers is. As shown in Figure 12(c), the rotation speed has no significant effect on the bearing force at 0.2 eccentricity, which means that the film force is mainly generated by the pressure difference. There is no contact between the fingers and the rotor. At 0.6 eccentricity, the bearing force increases and then does not change significantly with increasing rotation speed. This is because the increase in the number of contact fingers at 0.6 eccentricity leads to an increase in bearing force.
(a)
(b)
(c)
Figure 13 shows the effect of pressure difference. As shown in Figure 13(a), the leakage increases, and the number of contact fingers decreases as the pressure difference increases. The increase in pressure difference results in increased film forces, which reduces the number of contact fingers. Figure 13(d) shows the leakage of each finger at different pressure differences and eccentricities. Leakage at high eccentricity is not always high, which corresponds to the discussion of eccentricity above. As shown in Figure 13(b), the friction moment decreases with increasing pressure difference. However, the reasons for the decreases are different. At 0.2 eccentricity, all the fingers are under fluid lubrication. With increasing pressure difference, the film thickness increases, which means that the friction moment caused by the rotation speed decreases. At 0.6 eccentricity, the decrease is due to the reduction in the number of contact fingers. As shown in Figure 13(c), the bearing force decreases with increasing pressure difference. This is because the forces acting on the rotor by each finger are gradually the same. This phenomenon is further illustrated in Figure 13(e), which shows each finger’s force acting on the rotor at different pressure differences and eccentricities. When the pressure difference is 13 kPa, the bearing force is the combined force of the film forces at 0.2 eccentricity; the bearing force is the combined force of the film forces and the contact forces at 0.6 eccentricity. When the pressure difference is 200 kPa, the forces acting on the rotor by each finger are almost the same.
(a)
(b)
(c)
(d)
(e)
Figure 14 shows the effect of the number of fingers at a 100 kPa pressure difference, eccentricity of 0.6, and rotation speed of 8000 rpm. The average film force decreases with an increasing number of fingers because the circumferential length of the lifting pad (B) decreases. L/B in Figure 14 represents the ratio of the axial length to the circumferential length of the lifting pads. As shown in Figure 14, the leakage increases with the number of fingers; the friction moment and bearing force decrease with the number of fingers. The leakages of the square pads (L/B = 1) are greater than those of the wider pads (L/B = 0.5) because the square pads can generate more film force to lift the fingers, resulting in a higher increase in the clearances. The finger stiffness decreases with the increase in the number of fingers, and the finger deformations increase with the number of fingers under the same working environment. Thus, the clearances are increased, and the pressure distributions of each finger are almost the same, which results in the increase in the leakage and the decrease in the bearing force and friction moment.
(a)
(b)
(c)
4. Conclusion
In this paper, the static characteristics of noncontacting finger seals are studied when the seal works in an inappropriate environment or its design is unreasonable, which means contact may exist between the fingers and the rotor. A method to judge whether the fingers contact the rotor is proposed. The effects of various structural and environmental parameters are studied. Detailed conclusions are summarized as follows:(1)The contact between the fingers and the rotor reduces the leakage and significantly increases the bearing force and friction moment. This reduction or increase is proportional to the number of contact fingers.(2)The higher the pressure difference is and the lower the eccentricity is, the less likely the finger is to contact the rotor. The effect of rotation speed on the number of contact fingers is significantly smaller than pressure difference and eccentricity. Therefore, the pressure difference and rotor vibration amplitude should be considered the main environmental parameters in the design phase.(3)As the number of fingers increases, the film force, friction moment, and bearing force decrease, and the leakage increases. Wider pads can generate more film force than square pads. The square pads have smaller bearing forces and friction moments but have higher leakages than the wider pads. The greater the number of fingers is or the longer the axial length of the lifting pads is, the less likely the fingers come into contact with the rotor. The noncontacting feature is easier to achieve.
Appendix
The coefficients of equation (7) are as follows:
Abbreviations
| A: | Cross-sectional area |
| b: | Axial thickness of seal laminate |
| e: | Eccentricity distance |
| Ff: | Friction force between low-pressure finger laminate and backplate |
| Fhc, Fvc: | Horizontal and vertical components of contact force |
| Fh: | Film force |
| h: | Local film thickness |
| hif: | Film thickness of the i-th finger |
| hrif: | Film change due to rotation angle |
| : | Minimum pad wedge |
| I: | Moment of inertia |
| L: | Axial length of lifting pad |
| : | Unit bending moment |
| Mh: | Bending moment of noncontact finger |
| N: | Axial force |
| nf: | Number of fingers |
| O: | Diameter center of seal outer diameter |
| k: | Stiffness of high-pressure finger |
| p: | Pressure |
| ph: | Upstream pressure |
| Rs: | Stick arc radius |
| Rf: | Foot upper radius |
| Rcc: | Radius of circle of centers |
| uc: | Finger deformation of contact finger |
| uh: | Finger deformation under fluid lubrication |
| αul: | Curvature of the finger |
| μfr: | Friction coefficient between finger and rotor |
| θ1, θ2: | Circumferential angles of the contact position |
| θ: | Circumferential angle |
| φ: | Attitude angle |
| As: | Sectional film area |
| Alb: | Contact area between the low-pressure laminate and the backplate |
| E: | Young’s module |
| W: | Bearing force |
| Fk: | Force between two laminates |
| Fhh, Fvh: | Horizontal and vertical components of film force |
| Fc: | Contact force |
| : | Dimensionless film thickness |
| hinit: | Initial installation gap |
| : | Wedge under lifting pad |
| : | Maximum pad wedge |
| if: | Finger number |
| M: | Bending moment |
| m: | Number of circumferential grids |
| Mc: | Bending moment of contact finger |
| : | Unit axial force |
| n: | Number of axial grids |
| O2: | Diameter center of rotor |
| Plb: | Pressure that pushes the finger laminate against the backplate |
| : | Dimensionless pressure |
| pl: | Downstream pressure |
| Rr: | Rotor radius |
| Rb: | Finger base radius |
| u: | Deformations of fingers |
| u1, u2: | Finger deformations of positions 1 and 2 |
| αr: | Finger repeat angle |
| μ: | Gas viscosity |
| μlb: | Friction coefficient between low-pressure laminate and backplate |
| ε: | Eccentricity |
| θp: | Arc degree of lifting pad |
| Λx: | Compressibility number. |
Data Availability
The data calculation results 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 National Natural Science Foundation of China (Grants nos. 11972131 and 12072089) and National Science and Technology Major Project (Grant no. 2017-IV-0010-0047).