Abstract

During the working process of the turbine, some types of faults can cause changes in the vibration characteristics of the blades. The real-time vibration monitoring of the blades is of great significance to the stable operation and state-based maintenance. Torsional vibration is a kind of blade vibration modes and results in failures such as cracks easily. Thus, it is important to measure it due to the harmfulness of torsional vibration. Firstly, the principle of blade tip timing (BTT) is introduced, and the models of the blade are built to analyze the characteristics of torsional vibration. Then, the compressed sensing theory is introduced, and its related parameters are determined according to the measurement requirements. Next, based on the condition that the blade rigidity axis is not bent and bent, respectively, the layout method of sensors is proposed and the numerical simulation of the measurement process is performed. The results of the above two types of numerical simulation verify the proposed measurement method. Finally, by analyzing the influencing factors of measurement uncertainty, the optimization method of sensors’ layout is further proposed. This study can provide important theoretical guidance for the measurement of blade torsional vibration.

1. Introduction

The turbine is a type of engines widely used in various fields. As one of the most essential parts, blades compress air step by step to form high pressure gas that can provide power [1]. Under the severe working conditions, blades are prone to high-cycle fatigue and low-cycle fatigue, causing cracks and fractures, and even safety accidents [2]. Some types of faults can be reflected by changes in the vibration characteristics. Thus, the purpose of monitoring the blade state can be achieved by measuring the vibration. The traditional blade vibration measurement method is contact measurement, which is measured by strain gauges attached to the surface of the blade. However, the strain gauges have an impact on the vibration characteristics of the blade. On the contrary, the strain gauges have short life, low reliability, and high maintenance costs. BTT is a technique to measure vibration parameters of the rotating blade. The noncontact measurement based on it overcomes the aforementioned shortcomings. However, the current research in this field is mainly focused on bending vibration, while there is little research on torsional vibration. Therefore, it is essential to propose a torsional vibration measurement method to provide theoretical basis for noncontact measurement of different blade vibration modes.

Blade vibration modes can be divided into three types: bending vibration, torsional vibration, and coupled vibration based on both. Among them, the vibration around the inertia axis of the blade section is called bending vibration, and the vibration around the axis passing through the blade section rigidity along the height direction is called torsional vibration [3]. Nowadays, systematic research findings have been achieved about measuring the bending vibration of blades based on BTT. Rigosi et al. [4] improved the two-parameter-plot method to measure synchronous blade vibration more accurately. Neri [5] used nonharmonic Fourier series to improve the accuracy of frequency estimation based on BTT. Lin et al. [6] proposed a method for reconstructing unknown multifrequency coupling blade vibration signal. Pan et al. [7] proposed a method to reduce the uncertainty in the reconstruction of the blade vibration spectrum. Ou et al. [8] and Wang et al. [9] studied the synchronous blade vibration parameter identification based on BTT. Liu et al. [10] simulated the parameter identification of asynchronous blade vibration. Li et al. [11] simulated the parameter identification of multifrequency coupling blade vibration based on BTT. Liu et al. [12] proposed a data analysis method for a blade vibration measurement system based on BTT. However, there are also faults caused by torsional vibration. Casciar et al. [13] proved that blade torsional vibration has a significant impact on its performance and three-dimensional flow field. Dossena et al. [14] and Kallesøe and Hansen [15] proved that blade torsional vibration may cause problems such as engine performance degradation and even instability.

In recent years, researchers have studied the characteristics of blade torsional vibration. Liu [16, 17] established the mathematical model to study the characteristics of blade torsion based on the finite element method (FEM). Zheng et al. [18] studied the pressure distribution changes of the blade surface caused by blade torsion based on the fluid-solid coupling method. On the contrary, researchers have established different models and conducted corresponding experiments to simulate the measurement process of blade torsional vibration. Pfister et al. [19, 20] used laser Doppler sensors to measure blade deformation and vibration, but the transparent window on the casing required for this method would affect the original structure of the turbine. Ou [3] proposed a sensors’ layout method for blade torsion angle measurement, but did not conduct further research. Su et al. [21] proposed to install several dial gauges to measure the blade torsion angle. Chi et al. [22] built a high-temperature experimental platform to study the creep behavior of torsional blades and proposed a measurement method for the torsion angle. However, the installation of sensors is greatly restricted, which makes the above methods difficult to be used. On the contrary, considering the postprocessing of the measured signal and obtaining the torsional vibration signal efficiently and accurately will be a big challenge.

Blades are usually thin-walled irregular structures. The rigidity axis may bend during torsional vibration, causing blade coupled vibration [23]. In the measurement of bending vibration, only one BTT sensor is arranged along the axial direction of the blade, which is difficult to meet the measurement requirements of the complex blade vibration mode. On the contrary, the blade vibration signal measured based on BTT is undersampled and needs further processing [24, 25]. In other words, it is necessary to study the layout of sensors and signal processing methods for torsional vibration measurement.

The measurement method for bending vibration cannot be directly applied to torsional vibration. Thus, it is necessary to analyze the characteristics of the blade torsional vibration and propose the corresponding measurement method. However, research on this area has rarely been reported so far.

Based on the characteristic analysis of blade torsional vibration, a measurement method is proposed to provide theoretical guidance for the further promotion of BTT. Contents of this paper are organized as follows. In Section 2, the principle of BTT is introduced, and the models of the blade are established to analyze the characteristics of torsional vibration. In Section 3, the compressed sensing theory is introduced, and its related parameters are combined with measurement requirements. In Section 4, based on the condition that the blade rigidity axis is not bent and bent, respectively, the corresponding layout method of sensors is proposed and the numerical simulation of the measurement process is performed. Besides, by analyzing the influencing factors of measurement uncertainty, the optimization method of sensors’ axial layout is further proposed. Finally, conclusions are drawn in Section 5.

2. Introduction to BTT and Analysis of Blade Vibration Characteristics

2.1. Principle of BTT

According to [26], the circumferential layout of sensors when measuring bending vibration is proposed, as shown in Figure 1. The reference sensor is installed directly above the rotating shaft and generates a pulse signal every cycle of the rotating shaft, which acts as a time reference. The number of blades installed on the blisk is , and the installation angle of each blade along the circumferential direction relative to the reference sensor is . The number of BTT sensors installed on the casing is , and the installation angle of each BTT sensor along the circumferential direction relative to the reference sensor is . The rotation direction of the shaft is set to clockwise, and the rotation period is . When the blade does not vibrate, the theoretical time for the blade tip to reach the BTT sensor after rotation periods is calculated as follows:

The continuous theoretical time of arrival is , and the continuous actual time of arrival is . Thus, the difference between the two can be established as equation (2), and the continuous vibration signal of the blade can be established as equation (3):

In equation (3), is the distance from the blade tip to the center of the rotation shaft. For the blade , it is measured by BTT sensors and generates pulse signals in one rotation period. The final measured pulse signal sequence is as follows:where is Dirichlet sampling function and is the measured vibration signal that can be processed to obtain required vibration information. The measurement process is shown in Figure 2.

2.2. Analysis of Blade Vibration Characteristics

The torsional vibration of the blade is different from the bending vibration, which is analyzed to lay groundwork for subsequent research on measurement methods.

2.2.1. Analysis of Vibration Characteristics Based on the Finite Element Model

The structure of blisk area is shown in Figure 3. The model is composed of blades and a blisk. The blades are circumferentially fixed on the blisk by tenon connection and can be driven by it. It can be seen that the blade is thin and long with presetting angle.

Analysis of vibration characteristics is done based on the above model combined with ANSYS. Before the analysis, some parameters of the blade material have to be set, as shown in Table 1.

The blade and blisk are connected by a tenon, and the relative displacement of the two is almost negligible. Thus, fixed supports are set on both sides of the blade root to simulate actual working condition. The element type of the mesh is set to Solid187, and the total number of generated elements and nodes are 60122 and 99525, respectively. After the above steps, a finite element model (FEM) of the blade is generated, as shown in Figure 4.

Modal analysis of the above FEM is performed, and the results are shown in Figure 5 and Table 2. It can be seen that the 3rd and 6th modes are both torsional vibrations. Thus, it is of significance and importance to study them.

The schematic diagram of blade vibration is depicted in Figure 6. It can be found that, during bending vibration, the displacements of all nodes of the blade tip are almost the same, while the displacements of all nodes of the blade tip differ greatly during torsional vibration.

The target area of the vibration measurement based on BTT is the blade tip, so the mode shapes of the blade tip during torsional vibration are extracted, as shown in Figure 7. Sort the displacement of all the nodes to find the node with the smallest displacement, as shown in Table 3. It can be seen that the node with the smallest displacement is the same, which is marked in Figure 7. Furthermore, for the 6th mode, the displacement of this node is only 1.21% of the maximum displacement of the blade tip, which can be ignored. Thus, the 3rd and 6th modes are coupled vibration and torsional vibration, respectively. On the contrary, it is easy to be found that the axial distance between the node and both sides of the blade tip is 13.5 mm and 14.5 mm, respectively, when the blade is stationary.

2.2.2. Analysis of Vibration Characteristics Based on the Geometric Model

The blade vibrates torsionally around the rigidity axis, while the rigidity axis and the gravity axis may not coincide. The rigidity axis may bend due to inertial force, resulting in the blade coupled vibration [23]. During torsional vibration, the end point on the side of the blade tip of the rigidity axis is the center of twist, and it can be seen that its displacement is the minimum of all nodes. Thus, the node with the smallest displacement at the blade tip in Table 3 is the end point. In order to ignore the influence of the bending of rigidity axis, the following assumptions can be proposed:(1)It is assumed that the rigidity axis does not deform(2)It is assumed that the rigidity axis is straight and fixed vertically on the plane of the fixed end(3)It is assumed that the blade is rigid

Based on the above assumptions, the geometric model of the blade is built, as shown in Figure 8. The blade root is set as the origin of the coordinate system, the rigidity axis is set as Y-axis, and the distance between a certain blade cross section and the origin is set as . Researchers usually describe the torsional vibration based on the torsion angle. Considering the movement of , the vibration amplitude of the torsion angle is different due to the difference in the shape of each cross section. Before the further analysis, some parameters of blade torsional vibration have to be set, as shown in Table 4.

According to the geometric model of blade shown in Figure 8, the motion equation of blade cross section can be established as

The vibration amplitude of the torsion angle of the blade tip is set as ; then, the torsional vibration equation of the blade tip can be established as

The measurement target is vibration amplitude and vibration frequency .

3.1. Principle of Compressed Sensing Theory

For the measured blade, when BTT sensors are installed along the circumferential direction of the casing, only vibration displacements can be obtained per rotation period. Since the rotation frequency of the rotation shaft is much lower than the vibration frequency of the blade tip, the measured signal is an undersampled signal whose vibration parameters cannot be directly extracted. And, it is difficult to install multiple BTT sensors to meet Shannon’s sampling theory. Lin et al. [6] proposed a method to reconstruct the original blade vibration signal based on compressed sensing theory. It breaks through the limitation of Shannon’s sampling theory and is a promising method [27]. A brief description is given below.

Sparsity means that most of the elements of a signal are zero [28]. Most natural signals are not sparsity, but they can be sparsity based on a certain transform domain. The specific method is to disperse the original signal and multiply it by a transformation matrix to obtain a sparse vector. is set as a nonsparse vector, which can be expressed as follows after being decomposed by a sparse matrix :where is the sparse vector and the number of its nonzero elements is the sparsity. If the sparsity satisfies the property that , then matrix is the sparse matrix of vector .

The process of compressed sensing is described as follows. Firstly, the sparse representation of the original signal is obtained by measurement matrix . Next, the sparse vector of the original signal is obtained by sparse reconstruction. Finally, the original signal is obtained by inverting the sparse matrix.

The model of sparse reconstruction is as follows:where is the sensing matrix and can be calculated as and is a measurement vector whose number of elements is far less than that of . Equation (8) is described intuitively, as shown in Figure 9.

The difficulty of solving equation (8) is high (it is an NP-hard problem), and it cannot be solved in the certain polynomial time. However, is sparse; the only optimal solution can be obtained through the optimization models. Researchers have proposed some algorithms including greedy algorithms, convex optimization algorithms, and probabilistic algorithms to solve the problem. It is necessary to study the layout method of BTT sensors (determining the measurement matrix), the construction of the sparse matrix, and the selection of reconstruction algorithm.

3.2. Determination of the Related Parameters
3.2.1. Determination of the Sparse Matrix

The blade vibration signal is sparse in the frequency domain, which can be realized through the Fourier basis matrix . The method of obtaining this matrix is to perform Fourier transform on an identity matrix with the same dimension as the original signal. The sparse representation of the original signal can be calculated as follows:where is the discretized original signal and is the corresponding sparse signal.

3.2.2. Determination of the Measurement Matrix

Set the number of virtual BTT sensors evenly arranged along the circumferential direction on the casing as , and set th one of them as a total of real BTT sensors. Then, the corresponding measurement matrix when sampling the blade vibration signals of rotation periods is as follows:

The number of rows and columns of this measurement matrix are and , which can be calculated as and , respectively. There is exactly one element in each row of this matrix that is 1, and the column coordinate can be calculated as . Then, is the time when the target blade reaches the BTT sensor in each rotation period.

The uncorrelation or weak correlation between and can keep the original structure of the low-dimensional linear projection of the signal, which is an important condition for the application of compressed sensing theory. Candes and Tao [29] introduce restricted isometry property (RIP) to describe it. If there is a minimum constant that can make the following equation hold for all K-order sparse vectors , the matrix satisfies the K-order constraint isometry property as follows:where is the constrained isometry constant of matrix , and the matrix satisfies 2K-order constraint isometry property, satisfying the equation as follows:

The mathematical explanation of the above equation is as follows. The Euclidean distance between any pair of K-order sparse vectors is almost unchanged after linear transformation based on measurement matrix , which is of great significance to noise reduction. When , a unique K-order sparse vector can be reconstructed. However, the application of RIP is difficult, and the following equivalent method is proposed.

When the RIP is satisfied, each column vector of the sensor matrix is close to orthogonality [30]. In order to make sensing matrix satisfy the property, it is normalized as follows [31]:

The coherence coefficient of Q is used as the basis for the sensors’ layout, which can be calculated as follows [32]:where and are the vectors of column and column . The coherence coefficient is the normalized inner product between any two column vectors of the matrix , which reflects the correlation between them. The smaller the coherence coefficient, the stronger the noncorrelation of the sensing matrix and the greater the probability of accurately reconstructing the original signal.

The number of sensors proposed by Pan [33] is applied, the number of virtual BTT sensors is set to 25 and the 1st, 3rd, 9th, and 19th one of them are set as real BTT sensors.

3.2.3. Determination of Reconstruction Algorithm

The matching pursuit (MP) algorithm is used as the reconstruction algorithm in this study, which is a greedy algorithm based on iteration proposed by Mallet and Zhang [34, 35]. It regards a signal as linear combination of the elements of a dictionary, and the dictionary used for sparse reconstruction is the sensing matrix, in which each column vector is an element. After the measurement signal is obtained through compressed sensing, the signal can be expressed through sparse reconstruction.

In each iteration of the MP algorithm, a column vector with the greatest correlation with the residual is selected from the dictionary. The residual is the unexplained part of the original signal, which can be calculated as follows:where is the column of matrix . When the column that satisfies the above equation is determined, a new coefficient is added to the index set for approximating the original signal, and the residual and reconstruction result are updated accordingly:

The residual becomes smaller after iteration, and the algorithm stops when the norm of the residual is less than the preset threshold or the number of iterations is equal to the set value [36].

4. Analysis of Sensors Layout Based on Compressed Sensing Theory and Numerical Simulation

4.1. Analysis of Sensors Layout When Rigidity Axis Is Not Bent

The vibration amplitude of torsion angle is so small under real conditions that can be simplified as a straight line [37]. Furthermore, the axial layout of the BTT sensor without bending of the rigidity axis can be proposed, as shown in Figure 10, where point D is the end point on the side of the blade tip of the rigidity axis, point A is any point that does not coincide with point D, point C is the midpoint of blade tip, and point E and point F are the two end points of the blade tip. The BTT sensor is installed in the casing directly above point A, and the blade tip torsionally vibrates around point D. The circumferential layout of the BTT sensors refers to the layout used for bending vibration measurement.

Instantaneous torsional angle of blade tip can be calculated as follows:where is the presetting angle, is the axial distance between point A and D, is the displacement signal measured by the BTT sensor, and is the axial distance between the two end points of blade tip and point D when the blade is stationary, where and correspond to point E and point F, respectively. is the absolute value of the axial displacement of the two end points of blade tip when the blade is vibrating, where and correspond to point E and point F, respectively. is the axial distance between the two end points of blade tip and point D when the blade is vibrating, where and correspond to point E and point F, respectively. Substitute the displacement measured by the BTT sensor into equation (17), and solving it, the instantaneous torsional angle of blade tip can be obtained. It should be noted that the blade tip torsionally vibrates in two directions, so can be negative.

4.2. Numerical Simulation When Rigidity Axis Is Not Bent

Before the simulation, some parameters have to be set, as shown in Table 5. It can be found by Figure 10 that and are the minimum when the torsion angle is the maximum, which can be obtained such that and based on the FEA in Section 2. BTT sensor cannot measure the vibration signal when is too large, so there is a layout constraint as follows:

On the contrary, the installation position of the BTT sensor needs to correspond to the direction of selected as the constraint.

The sampling rate of virtual BTT sensors and real BTT sensors can be calculated, respectively, as follows:

On the basis of satisfying Shannon’s sampling theory, the upper limits of the measurable vibration frequency can be calculated, respectively, as follows:

Substitute the parameters in Table 5 into the above equation and solve it, the upper limits of the measurable vibration frequency with and without compressed sensing technology are calculated to be 781.25 Hz and 125 Hz, respectively, which shows that compressed sensing technology is of great significance for measurement under the condition of the limited sampling rate.

The simulation process is shown in Figure 11. Original vibration signal of the blade tip torsion angle is generated based on the parameters in Table 5, and it is measured to obtain the undersampled signal based on compressed sensing technology and equation (17). The above simulation results are both intercepted from the initial 0.05 s in order to better observe; the subsequent simulation results are still processed in the same way.

The undersampled signal of torsion angle is reconstructed to obtain the original signal, and the result is shown in Figure 12.

It is easy to find that the original and reconstructed signals in time domain and frequency domain are almost identical. The residual obtained by calculating the cumulative deviation of all sampling points in time domain is 2.8302e − 13, which meets the measurement requirements.

Next, the axial distance between point A and D is changed to verify the reliability of the proposed measurement method. The set value of d1 is shown in Table 6.

The reconstruction residual is calculated, as shown in Figure 13. It can be found that the results meet the measurement requirements.

The real working conditions are noisy; enhancing the signal is a way to improve the signal-to-noise ratio. Based on Figure 10, it is easy to find that, with the increase of d1, the displacement signal measured by the BTT sensor also increases so that the signal-to-noise ratio can be increased.

A white noise signal is added to the original signal of torsion angle, as shown in Figure 14. Afterwards, numerical simulation is performed based on d1 in Table 6 and the reconstruction residual trend is plot, as shown in Figure 15. It can be found that, as d1 increases, the reconstruction residual is in a decreasing trend.

The positions of point E and point F when is the minimum can be determined based on FEA. For the blade model in Section 2, the 6th-order mode shape of the blade tip under a certain working condition is shown in Figure 16. The absolute values and of the maximum axial displacement of point E and F are 0.37 mm and 0.79 mm, respectively. Substitute them with the values of and when the blade is stationary into equation (21) and solve it; the values of and can be obtained to be 13.13 mm and 13.71 mm, respectively. When the installation position of the BTT sensor satisfies equation (22), the signal-to-noise ratio can be maximized under the premise of meeting the measurement requirements:

4.3. Analysis of Sensors Layout When Rigidity Axis Is Bent

When the rigidity axis is bent, it is difficult to measure torsion angle of the blade tip with only one BTT sensor. Furthermore, the axial layout of the BTT sensors with bending of the rigidity axis is proposed, as shown in Figure 17. Point D is the end point on the side of the blade tip of the rigidity axis, where point D′ and D″ correspond to the stationary and vibrating blade tip; point X is the intersection point of the radial projection when the blade tip is stationary and vibrating. Point C is the midpoint of blade tip. Point A and point B are arbitrary points, where point A does not coincide with point D. The two BTT sensors are installed in the casing directly above point A and point B, respectively. The circumferential layout of the BTT sensors refers to the layout used for bending vibration measurement.

Instantaneous torsional angle of blade tip can be calculated as follows:where is the presetting angle, is the axial distance between point A and B, is the axial distance between point B and X, is the axial distance between point D′ and D″, which is easy to find that the two are positively correlated, and are the displacement signals measured by BTT sensor A and B, respectively, is the axial distance between the two end points of blade tip and point B when the blade is stationary, where and correspond to point E and point F, respectively, is the absolute value of the axial displacement of the two end points of the blade tip when the blade is vibrating, where and correspond to point E and point F, respectively, and is the axial distance between the two end points of blade tip and point B when the blade is vibrating, where and correspond to point E and point F, respectively. Substitute the displacements and measured by the BTT sensors into equation (23) and solve it; the instantaneous torsional angle of blade tip can be obtained.

4.4. Numerical Simulation When Rigidity Axis Is Bent

Before the simulation, some parameters have to be set, as shown in Table 7. After that, and are calculated by the proposed parameters and the FEM in Section 2. It can be found by Figure 17 that and are the minimum when the torsion angle is the maximum, which can be obtained as and . On the contrary, BTT sensors cannot measure the vibration signal when is too large, and there is a layout constraint as follows:

The simulation process is shown in Figure 18. Original vibration signal of the blade tip torsion angle is generated based on the parameters in Table 7, and it is measured to obtain the undersampled signal based on compressed sensing technology and equation (23). The above simulation results are both intercepted from the initial 0.1 s in order to better observe; the subsequent simulation results are still processed in the same way.

The undersampled signal of torsion angle is reconstructed to obtain the original signal, and the result is shown in Figure 19.

It is easy to find that the original and reconstructed signals in time domain and frequency domain are almost identical. The residual obtained by calculating the cumulative deviation of all sampling points in time domain is 7.1804e − 04, which meets the measurement requirements.

The positions of point E and point F when is the minimum can be determined based on FEA. For the blade model in Section 2, the 3rd-order mode shape of the blade tip under a certain working condition is shown in Figure 20. The absolute values and of the maximum axial displacement of point E and F are 0.06 mm and 0.18 mm, respectively. Substitute them with the values of and when the blade is stationary into equation (25) and solve it; the values of and can be obtained to be 4.44 mm and 23.32 mm, respectively. When the installation position of the BTT sensor satisfies equation (26), the signal-to-noise ratio can be maximized under the premise of meeting the measurement requirements.

4.5. Optimization of Sensors’ Layout for Measurement Uncertainty

The bending of the rigidity axis causes actual sampling time to deviate from theoretical sampling time, which makes measured values of torsion angle deviate from the true values, resulting in measurement uncertainty. Next, change maximum axial distance between point B and X to further study the influence of uncertainty. The set value of d3max is shown in Table 8.

Calculate the reconstruction residual and plot the residual trend, as shown in Figure 21. It can be found that reconstruction residual increases with d3max.

In order to eliminate the influence of measurement uncertainty, the position of point X must be obtained, which is difficult to achieve based on BTT. However, from the previous analysis, it can be found that if the BTT sensor B is installed in the casing directly above point D′ to make as small as possible, the uncertainty can be minimized within the controllable range.

5. Conclusions

The vibration of blades is a threat to the operation of turbines, and some faults can also be found through it. Therefore, measuring blade vibration efficiently and accurately is significant. Thus, a measurement method for torsional vibration based on BTT is proposed, which provides theoretical guidance for the promotion of BTT. Firstly, the blade models are built to analyze the characteristics of torsional vibration. Then, based on the conditions of no bending and bending of rigidity axis, the layout method of sensors is proposed, respectively. Next, numerical simulation is performed to verify the proposed measurement method, and the sensor axial layout optimization method that can improve the signal-to-noise ratio is proposed. Finally, an optimization method for sensors’ axial layout is proposed based on the analysis of measurement uncertainty. Conclusions can be drawn as follows:(1)Based on the measurement method proposed in this paper, torsional vibration of the blade can be measured to obtain vibration parameters including frequency and amplitude(2)The installation position of BTT sensors can affect the signal-to-noise ratio, and an optimization method that can improve it is proposed(3)Aiming at measurement uncertainty, an optimization method of BTT sensors’ axial layout that can reduce it is proposed

The findings in this research provide theoretical guidance for the measurement of blade torsional vibration. In future work, based on the more complex blade structure and working conditions, the proposed blade model can be improved to further study the measurement method, and the corresponding experimental verification can be carried out, too.

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This work was supported by the Natural Science Foundation of Hunan Province, China (Grant no. 2019JJ50721).