Abstract
Blade tip timing (BTT) technology is the most effective means for real-time monitoring of blade vibration. Accurately extracting the time of blade tip reaching the sensors is the key to ensure the accuracy of the BTT system. The tip clearance changes due to various complex forces during high-speed rotation. The traditional BTT signal extraction method does not consider the influence of tip clearance change on timing accuracy and introduces large timing errors. To solve this problem, a quadratic curve fitting timing method was proposed. In addition, based on the measurement principle of the eddy current sensors, the relationship among the output voltage of the eddy current sensor, tip clearance, and the blade cutting magnetic line angle was calibrated. A multisensor vibration parameter identification algorithm based on arbitrary angular distribution was introduced. Finally, the experiments were conducted to prove the effectiveness of the proposed method. The results show that in the range of 0.4 to 1.05 mm tip clearance change, the maximum absolute error of the timing values calculated by the proposed method is 26.0359 us, which is much lower than the calculated error of 203.7459 us when using the traditional timing method. When the tip clearance changed, the constant speed synchronous vibration parameters of No. 0 blade were identified. The average value of the vibration amplitude is 1.0881 mm. Compared with the identification results without changing tip clearance, the average value error of the vibration amplitude is 0.0017 mm. It is proved that within the blade tip clearance variation of 0.4 to 0.9 mm, the timing values obtained by the proposed timing method can accurately identify the vibration parameters of the blade.
1. Introduction
Blades are considered to be the core component of aeroengine compressor. Real-time monitoring of blade vibration state can effectively prevent blade failures [1]. At present, the BTT technology is the most promising real-time monitoring method for blade vibration [2, 3]. The BTT technology installs several sensors on the stationary casing, records the time series of blade tips reaching the sensors, and then determines the vibration parameters of the blade [4]. Therefore, accurate measurement of the blade tip timing time is the key of the BTT technology. At present, the optical method [5], the capacitance method [6], the microwave method [7], and the eddy current method [8] can realize the measurements of blade vibration and tip clearance. In the heavily polluted environment, due to the surface pollution of the probe, the measurement accuracy of the optical method decreases. The measurement accuracy of the capacitance method is affected by insulation, stray capacitance, and edge effects. The measurement accuracy of the microwave method is affected by the spatial filtering effect, and the system is expensive. Compared with the above three methods, the eddy current method is more suitable for long-term monitoring of blades. In the working process of rotating machinery, high-speed rotating blades not only are subject to the unbalanced centrifugal force of the rotor but also interact with the airflow to form an uneven and unstable airflow force [9]. These forces lead to changes in the tip clearance between the sensor probe and the blade. The changes in the tip clearance will change the amplitude of the pulse signals measured by the sensor and lead to the introduction of timing errors when extracting the blade timing time.
In addition to the traditional timing identification, there are also constant ratio timing identification and double-threshold leading-edge timing discrimination. The constant ratio timing identification method [10] uses half of the height of the rising edge of the pulse signal as the blade arrival time. However, this method has high signal processing difficulty and is still affected by the change of the blade tip clearance. The double-threshold leading-edge timing discrimination method [11] belongs to the error compensation identification method, which uses multiple thresholds to calibrate the error compensation relationship of the same signal and eliminates the time drift error caused by the single threshold. However, this method uses twice single-threshold leading-edge timing discrimination, which will also introduce new errors while compensating for errors. In recent years, several new time identification methods have been proposed. Duan et al. [12] proposed a method to use the intersection of the linear partial extension line of the leading and trailing edges of the pulse signal as the arrival time of the blade based on the fiber bundle sensor. However, the difference between the blades and the difference in the vibration state of the blade make it difficult to determine the appropriate slope ratio k of the linear portion of the leading and trailing edges. Zhang et al. [13] proposed a timing method for the middle time of the pulse signal based on the microwave sensor. This method compares the leading and trailing edges of the pulse signal with a fixed trigger voltage, and the time average value obtained is approximately the peak time of the pulse signal. However, this method does not consider the influence of the blade vibration phase on the measurement results. Zhang et al. [14] proposed an error compensation method based on the eddy current sensor. This method uses tip clearance information to compensate for tip timing information. However, new errors will inevitably be introduced in the process of error compensation.
In this paper, the influence of the tip clearance change on the accuracy of the BTT measurement system is studied. According to the collected pulse signal sampling points, based on the eddy current sensor, a method for accurately extracting the blade timing time when the blade tip clearance changes is proposed. The timing time obtained by this method can be used to identify the constant speed synchronous vibration parameters of the blade when the tip clearance changes. The paper is organized as follows. Section 2 briefly describes the principle of BTT technology, analyzes the timing errors introduced by the traditional rising edge timing method when the tip clearance changes, and proposes a new timing method. In Section 3, the static calibration experiment based on the eddy current sensor is conducted. Section 4 introduces the multisensor vibration parameter identification algorithm for arbitrary angular distribution. In Section 5, the experiments are conducted to demonstrate the effectiveness of the proposed method. Finally, the conclusions are drawn in Section 6.
2. Blade Tip Timing Technology
2.1. Principle of Blade Tip Timing Measurement
The principle of blade tip timing measurement is shown in Figure 1. A speed synchronization sensor is installed at the position corresponding to the rotating shaft to record the speed synchronization signals. Several BTT sensors (TIP0∼TIP3) are installed along with the stationary casing to record the BTT signals generated when the tips of the rotating blades sweep through the sensor probes [15]. If the blade does not vibrate, the time for the blade tip to reach the BTT sensor can be calculated based on the mounting position of the blade on the rotor and the rotor speed. If the blade vibrates, the tip of the blade will reach the BTT sensor (leading or lagging) [16, 17]. At this time, the actual arrival time of the blade is not equal to the theoretical arrival time when it is not vibrating and forms a time difference . The time differences recorded by multiple BTT sensors are collated to form a time difference sequence . Through the calculation of the time difference sequence , we can get the vibration displacements of the blade, and then with the help of the corresponding blade vibration parameter identification algorithm, we can identify the blade vibration amplitude and other parameters [18].
2.2. Error Caused by Variation of Tip Clearance
The blade tip clearance of the high-pressure rotor of the engine is changed due to the influence of centrifugal deformation and thermal deformation during actual rotation. The traditional timing method of processing the BTT signal takes the moment of the intersection of a certain fixed threshold cutting level and the pulse signal rising edge as the timing time of the blade. This method does not take into account the effect of the tip clearance change on the timing accuracy of the measurement system. It leads to the introduction of large timing errors when extracting the time difference sequence which eventually reduces the accuracy of blade vibration parameter identification.
The traditional BTT signal processing method is shown in Figure 2. When the tip clearance does not change, the BTT sensor measures the pulse signals which are shown in I, and the time interval between the two pulse signals is . When the tip clearance changes, the BTT sensor measures the pulse signals which are shown in II, and the time interval between the two pulse signals is . The reason is that when the tip clearance changes, the circumferential position of the blade trigger level is different from that when the tip clearance does not change. It leads to the introduction of a time error in the calculation of the blade vibration displacement.
The vibration displacement errors due to the variation of the tip clearance can be calculated using the following equation:where refers to the timing errors introduced by the traditional BTT signal processing method when the tip clearance changes and is the rotational angular velocity of the rotating blade.
2.3. Quadratic Curve Fitting Timing Method
The peak voltage and peak time of the pulse signal can be obtained without determining polynomial order and solving the polynomial by using the quadratic curve fitting timing method to process the pulse signal. This method solves the problem of difficulty in determining the peak time because the peak part of the pulse signal is relatively flat. This method discards the part of the pulse signal that differs significantly from the parabola, and quadratic curve fitting is performed on the remaining sampling points to obtain the peak time and the peak voltage of the pulse signal. This method has been used to solve the half-sine shock waveform peak calculation problem [19].
The sampling sequence of the output voltage of the BTT sensor is , , …, , and the sampling sequence of the corresponding time is , , …, . The function expression of the quadratic curve fitting can be expressed using the following equation:
From the coefficients of each term in equation (2), it can be determined that is the peak time of the fitted waveform and is the peak voltage of the fitted waveform. When using this method, the range of interception being too small or too large is not conducive to the accuracy of the quadratic curve fitting. Therefore, intercepting 20% to 50% of the peak amplitude of the pulse signal can be expected to yield good fitting results.
3. Experimental Measurement of Static Calibration Based on Eddy Current Method
3.1. Measurement Principle of Eddy Current Method
The principle of the eddy current method is the use of metal cutting magnetic induction lines to produce magnetic field changes. When the metal conductor is placed in the changing magnetic field, the induced current of the closed circuit is formed in the conductor. This phenomenon is called the eddy effect. The eddy current sensors are established based on the eddy effect. The eddy current sensor probe is placed near the blade. When the probe coil is applied with alternating current, the space around the coil produces a sinusoidal alternating magnetic field. The blade continuously cuts the magnetic induction lines through rotation and radial displacement in this alternating magnetic field and generates the induced current. This induced current generates an alternating magnetic field in the opposite direction of the original magnetic field and counteracts part of the original field [20]. It leads to the change in the equivalent impedance of the coil.
The functional expression of the equivalent impedance Z can be expressed using the following equation:where is the electrical conductivity of the blade, is the magnetic permeability of the blade, is the passing frequency, is the size factor of the coil, and is the tip clearance.
For the known eddy current sensor and the blade, the coil size factor , the blade electrical conductivity , the blade magnetic permeability , and the passing frequency are fixed. Therefore, the equivalent impedance becomes a single-valued function of the tip clearance [21]. Through the eddy current sensor conversion circuit, the change in the equivalent impedance is finally converted into a change in the output voltage . However, the sensor probe diameter selected in this paper is 8 mm, and the thickness of the blade to be measured is 2 mm. The blade thickness is much smaller than the sensor probe diameter. Therefore, the sensor output voltage is related to the tip clearance , the blade thickness , and the angle of the blade cutting the magnetic induction line [22]. Therefore, the output voltage of the eddy current sensor can be expressed using the following equation:
For a defined blade to be measured, the blade thickness is constant. Therefore, the output voltage of the eddy current sensor can be further expressed using the following equation:
As can be seen from equation (5), when the blade thickness is constant, the output voltage of the eddy current type sensor is only related to the tip clearance and the angle of the blade cutting the magnetic induction lines. Therefore, as long as one of the variables is guaranteed to be constant, the output voltage of the eddy current sensor remains a single-valued function of another variable.
The centerline of the blade and the centerline of the eddy current sensor probe are kept in the same line, and the tip of the blade is gradually moved away from the eddy current sensor probe. A voltage value is obtained each time the blade tip clearance is changed. The tip clearance is the horizontal coordinate, and the eddy current sensor output voltage is the vertical coordinate. The relationship between the output voltage of the eddy current sensor and the tip clearance can be expressed using the following equation:
Keeping the tip clearance constant, the blade starts with one end of the eddy current sensor probe and rotates at the same angle each time. Record the voltage value of each angle until the blade is rotated to the other end of the eddy current sensor probe. The angle of the blade cutting magnetic induction lines is the horizontal coordinate, and the output voltage of the eddy current sensor is the vertical coordinate. The relationship between the output voltage of the eddy current sensor and the angle of the blade cutting magnetic induction lines can be expressed using the following equation:
3.2. Building Static Calibration System Experimental Device
The experimental device of the static calibration system is shown in Figure 3. It mainly contains two parts: the radial moving system and the circumferential rotation system. The XYZ three-axis displacement table of the radial movement system and the R-axis rotation table of the circumferential rotation system are fixed on both sides of the base plate. The disk is fixed on the R-axis rotation bench, and the blade is fixed on the disk. The eddy current sensor is fixed on the XYZ three-axis displacement bench through the L-shaped bracket. The static calibration of the blade is completed by the regulation of the XYZ three-axis displacement bench and the R-axis rotation bench. To facilitate the measurement and reading, only one blade is installed on the disk each time during the static calibration experiment. All eight blades are installed on the disk when the BTT measurement experiment is conducted.
3.3. Static Calibration Experiments
The macroscopic process of blade static radial calibration is shown in Figure 4. The static radial calibration experiment takes No. 0 blade as an example. In the measurement range of the eddy current sensor, according to the order of the tip clearance from small to large, the voltage values corresponding to 10 tip clearances in the range of 0.1 mm to 0.55 mm were read with 0.05 mm as the step length. To ensure the accuracy of the experimental data, the measurement was repeated five times for each tip clearance, and the average value of the five measurements was taken as the voltage value corresponding to the blade under the tip clearance. After completing the static calibration of No. 0 to No. 7 blades in turn and obtaining all data, the data were fitted with the help of MATLAB. The relationship curves between the output voltage and the tip clearance of No. 0 to No. 7 blades were measured, respectively. The static radial calibration curve of No. 0 blade is shown in Figure 5. The fitting equation corresponding to the curve is presented using the following equation:
The square of the correlation coefficient is 0.9994. It can be seen from Figure 5 that the output voltage of the eddy current sensor is about −1000 mV when the tip clearance is 0.1 mm, and the output voltage of the eddy current sensor decreases by about 450 mV for every 0.05 mm increase in the tip clearance. Also, the output voltage varies nearly linearly with the change of the tip clearance .
The macroscopic process of blade static circumferential calibration is shown in Figure 6. The static circumferential calibration experiment also takes No. 0 blade as an example. Turning the R-axis rotation bench would gradually turn the blade out of the measurement range of the eddy current sensor. The rotation stopped when the voltage value displayed in the software interface was about −5000 mV. The angle at this position was recorded as 0°. Each scale of the R-axis rotation bench is 0.46, and the blade was recorded an output voltage value for every five scales of rotation. When the output voltage value displayed on the software interface was again close to −5000 mV, the blade was at the edge of the eddy current sensor measurement range, the rotation was stopped, and all circumferential data under this blade tip clearance were collected. Experimental measurements were performed every 0.05 mm to complete the static circumferential calibration of No. 0 blade at the tip clearances of 0.1 to 0.55 mm. To ensure the accuracy of the experimental data, the calibration was repeated five times for each group of experiments and the average value was taken as the static circumferential calibration data of No. 0 blade at this tip clearance. After finishing the calibration of No. 0 to No. 7 blades in turn and obtaining all data, the data were fitted with the help of MATLAB. The relationship curves between the output voltage and the angle of the blade cutting magnetic inductance of No. 0 to No. 7 blades were measured, respectively. The static circumferential calibration curve of No. 0 blade is shown in Figure 7. The corresponding fitting equations of the curves are shown in equations (9)–(18).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
The square of the correlation coefficient of the fitting equations is not less than 0.9992 for any tip clearance from 0.1 to 0.55 mm. It can be seen from Figure 7 that in the blade static axial calibration experiment, since the range of magnetic field cut by the blade gradually decreases with the increase of the blade tip clearance , the number of voltage values collected by the eddy current sensor probe also gradually decreases. Also, the output voltage varies parabolically with the blade cutting magnetic induction line angle .
4. Constant Speed Blade Synchronous Vibration Parameter Identification Method
Ou et al. [23] proposed a multisensor vibration parameter identification algorithm based on arbitrary angular distribution, which does not restrict the installation angle and layout method of the BTT sensors and can identify vibration parameters for different states and forms of blades. In this paper, this parameter identification algorithm is used to identify the constant speed blade synchronous vibration parameters when the tip clearance changes.
The speed synchronization sensor is installed at the position corresponding to the rotating shaft. The installation angle between the speed synchronization sensor and the BTT sensor TIP0 is . The vibration displacement can be calculated using the following equation:where is the distance from the tip of the blade to the center of the rotating shaft, is the rotational angular velocity of the rotating blade, is the time interval between the blade tip passing through the speed synchronization sensor and TIP0, is the angle between the speed synchronization sensor and TIP0, is the vibration amplitude, is the vibration frequency doubling, is the initial phase, and is the constant offset.
The angle between the speed synchronization sensor and the BTT sensor TIPi is, and the time interval between the blade tip passing through the speed synchronization sensor and the BTT sensor TIPi is . Therefore, equation (19) is expressed as follows:
Expanding equation (20), we get
Substituting, , , …, and into equation (21), we get
Write equation (22) in the form of matrix equation Y = BX, where the matrices of Y, B, and X are expressed as equations (23)–(25).
From the matrix equation , it is known that matrix is the vibration displacement of the blade, the parameter in matrix is unknown, and the parameters , , and in matrix are unknown. There are four unknown parameters. Therefore, the vibration parameters of the blade can be identified by using four BTT sensors. Among them, the vibration frequency doubling is a positive integer, which will be traversed in a certain range of frequency doubling and bring all the values of in the range into the matrix equation . The solution vector can be calculated using the least squares method. The calculation equation is
Substitute the solution vector into the matrix equation and subtract the actual vibration displacement of the blades to obtain their residuals. The residuals can be calculated using the following equation:where .
The variance of the residual is the error between the traversal estimate value and the actual measured value. The variance can be calculated using the following equation:
Let the traversal frequency doubling ; in the range of traversal frequency doubling , the frequency doubling when the variance of residual obtains the minimum value is the actual frequency doubling of blade vibration. At this time, the traversal frequency doubling , and use it to find out the solution vector . Finally, the vibration frequency, vibration amplitude , initial phase, and constant offset are shown in equations (29)–(32).
Vibration frequency :
Vibration amplitude :
Initial phase :
Constant offset :
5. BTT Measurement Experiment
5.1. Building BTT Measurement Experiment Bench
The BTT measurement experimental bench is divided into the rotating blade noncontact measurement device shown in Figure 8(a) and the signal acquisition and data processing system shown in Figure 8(b). The rotating blade noncontact measurement device consists of a bladed disk, the blade tip timing sensors, the noncontact electromagnetic exciters, a rotating shaft, a speed synchronization sensor, a V-belt, and a servo motor. The servo controller controls the servo motor speed. The signal acquisition and data processing system consists of a power amplifier, a data collector, an eddy current sensor power supply box, and the Adcras data processing software from the CRAS software library. The experiments of BTT measurement are divided into the experiment of improving the accuracy of the BTT system and the experiment of constant speed blade synchronous vibration parameter identification.
(a)
(b)
5.2. Experiment of Improving Accuracy of BTT System
This experiment is based on the principle of noncontact measurement of rotating blade vibration. The speed synchronization sensor was installed at the position corresponding to the rotating shaft, and the installation angle was 0°. The reflective marker was affixed to the same horizontal position with No. 0 blade, and No. 0 blade was taken as the experimental research object. TIP0 was installed at the horizontal position on the stationary casing, and the installation angle was also 0°. The sampling frequency of the experiment was 25.6 kHz, and the revolution speed was 400 r/min.
References [24, 25] show that an increase in blade tip clearance reduces aeroengine efficiency and increases fuel consumption. The initial blade tip clearance of aeroengine compressors is generally about 1 mm. With the increase in rotational speed, the blade tip clearance decreases due to the increase in centrifugal force. Through the actual measurements of low rotational speed and the premise of ensuring the safe operation of the experimental device, the dynamic tip clearance of No. 0 blade is measured in the range of 0.4 to 1.05 mm. The initial tip clearance was 0.4 mm, and the tip clearance was increased by 0.05 mm each time. A total of 14 groups of experiments were conducted. The time of No. 0 blade arrived at the speed synchronization sensor was recorded as , and the pulse signals of No. 0 blade at different tip clearance were measured by TIP0. The highest point voltage collected was taken as the peak voltage, and the sampling points in the range of 20%, 30%, 40%, and 50% down of the peak voltage were intercepted. The quadratic curve fitting of the intercepted sampling points was performed to calculate the timing value of No. 0 blade. When the tip clearance was 0.4 mm, the fitting curves of the experimental data and the quadratic fitting curves of different interception ranges are shown in Figure 9.
(a)
(b)
(c)
(d)
As can be seen from Figure 9, the curves obtained by quadratic fitting of the sampling points in the range of 20%, 30%, 40%, and 50% down of the peak voltage basically overlap with the experimental data fitting curves; the smaller the intercepted range, the higher the degree of overlap. Therefore, the peak value coordinate of the quadratic fitting curve is approximated as the peak value coordinate of the pulse signal. In theory, the peak time of the pulse signal is used as the timing time of the blade without the influence of the tip clearance change. Therefore, this method makes the timing easier while greatly reducing the influence of the tip clearance change on the timing accuracy. When the tip clearance was 0.4 mm, the quadratic fitting results of different interception ranges of the pulse signal are shown in Table 1.
As can be seen from Table 1, the peak time obtained by quadratic fitting of the sampling points in the range of 20%, 30%, 40%, and 50% down of the peak voltage is very close to each other, and the maximum value of the absolute error is 1.1911 us. The peak voltage decreases gradually as the intercept range increases. If a good interception range is chosen, the dynamic tip clearance of the blade is expected to be measured based on the fitting peak voltage.
The timing value error comparison curves for different interception ranges are shown in Figure 10. As can be seen from Figure 10, the quadratic curve fitting timing method is very little influenced by the change of the tip clearance, and in the range of 0.4 to 1.05 mm of the tip clearance change, the maximum absolute error of the timing value is 26.0359 us.
These four sets of curves were averaged and compared with the traditional rising edge timing method. The timing value error comparison curves between the traditional rising edge timing method and the quadratic curve fitting timing method are shown in Figure 11.
As can be seen from Figure 11, when the traditional rising edge timing method is used, the larger the change of the tip clearance, the more obvious the error. When the change value of the tip clearance reaches 0.65 mm, the timing value error is 203.7459 us. When the quadratic curve fitting timing method is used, the maximum absolute error of the timing value is 23.5606 us, which reduces the timing error introduced by the change of the tip clearance.
5.3. Experiment of Constant Speed Blade Synchronous Vibration Parameter Identification
TIP0 to TIP3 were installed along with the stationary casing in order of 0°, 14°, 29°, and 43°, The tip clearances between TIP0 to TIP3 and No. 0 blade were adjusted to 0.4 mm, 0.55 mm, 0.75 mm, and 0.9 mm in turn. The rotational speed of the experiment was 456 r/min. To make the blade vibrate during rotation, the exciter is connected to a power amplifier and placed on one side of the bladed disk to excite the blade as it rotates.
The speed synchronization signals and the BTT signals were collected in the rotating state, and the time of No. 0 blade arriving at the speed synchronization sensor was recorded as . The tip timing data of No. 0 blade were obtained by TIP0 to TIP3. The experimental data were processed by the quadratic curve fitting timing method, and the obtained timing values were used for the identification of constant speed synchronous vibration parameters. The parameter identification results of No. 0 blade are shown in Table 2.
As can be seen from Table 2, in the range of 20%, 30%, 40%, and 50% down of the peak voltage, the vibration frequency doubling of No. 0 blade is obtained as 8 and the vibration amplitudes are 1.0806 mm, 1.0957 mm, 1.0873 mm, and 1.0889 mm. The average value of the vibration amplitude is 1.0881 mm. However, the timing value obtained by the traditional method is greatly affected by the change of the blade tip clearance, and it has a large identification error. When TIP0 to TIP3 and No. 0 blade have the same tip clearance, the average value of the vibration amplitude of No. 0 blade identified is 1.0864 mm. The error of the average value of vibration amplitude is 0.0017 mm. The identification results are very close. It is proved that within a certain range of the tip clearance change, the timing values obtained by the quadratic curve fitting timing method can accurately identify the vibration parameters of the blade.
6. Conclusions
The tip clearance between the blade tip and the casing changes during the actual rotation of the engine high-pressure rotor. It causes the traditional timing method to introduce a large timing error, and the result is that the vibration parameters of the blade are no longer accurately identified. To solve this problem, a new timing method was proposed and validated by experiments. The main contents and conclusions of this study can be summarized as follows:(1)The influence of tip clearance variation on the timing accuracy of the measurement system was analyzed, and a quadratic curve fitting timing method was proposed. This method can directly derive the peak voltage and peak time from the collected waveform data.(2)The eddy current method of the tip clearance measurement was analyzed. The output voltage of the eddy current sensor was calibrated to have a linear relationship with the tip clearance and a parabolic relationship with the angle of the blade cutting magnetic induction lines.(3)In the experiment of improving the accuracy of the BTT system, the timing errors of the two methods were compared when the tip clearance changed. The experimental results show that the quadratic curve fitting timing method effectively improves the accuracy of the BTT system when the tip clearance changes.(4)A multisensor vibration parameter identification algorithm based on arbitrary angular distribution was introduced, and the parameter identification experiments were conducted. The experimental results show that within the blade tip clearance variation of 0.4 to 0.9 mm, the timing values obtained by the quadratic curve fitting timing method can accurately identify the vibration parameters of the blade.
Data Availability
The data used to support the findings of this paper are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The authors are grateful for the support provided by the National Natural Science Foundation of China (no. 51505206), Liaoning Provincial Natural Science Foundation Guidance Program Project (no. 2019-ZD-0694), Liaoning Provincial Department of Education Scientific Research General Project (no. L2014246), and Liaoning University of Technology Teacher Research Start Fund Project (no. X201202).