#### Abstract

Frequent failure of the walking wheel seriously restricts the performance of the whole shearer. To reduce the failure rate and improve the working performance of the walking wheel, the meshing characteristics of the tooth pin are researched. The dynamic state equations of the walking wheel and the contact force model of tooth pin meshing are established. The rigid-flexible coupling simulation model of tooth pin meshing is built. The load distribution characteristics of the walking wheel are analyzed, as well as the effects of impact load amplitude and duration. Results show that the curve of the longitudinal load distribution coefficient (*K*_{β}) of the contact area is W-shaped, with a maximum of 1.325 at the moment of a single tooth contact. The end of the transition curve is the most serious position of the longitudinal load imbalance at the tooth root. In addition, on the impact moment, *K*_{β} tends to decrease and maximum stress obviously increases with the increase in impact load under 40%; the material at the contact position will fail under an extra 39% instantaneous impact load. Furthermore, with the impact load of 30%, the influence of load impact duration under 0.5 s on the meshing characteristics of the walking wheel is relatively faint. The results provide some guidance for the design optimization of the walking wheel and provide a reference for improving the reliability of the shearer.

#### 1. Introduction

Ensuring high production and efficiency of coal mining has important practical significance for world energy security. Shearer is the key equipment used in coal mining [1, 2], and the failure of gear is the main cause of shearer malfunction [3]. The walking wheel meshes with the pin rows on the scraper to drive the shearer to move; compared with other transmission gears, its working environment is worse and its meshing form is special [4]. Therefore, it is necessary to study the dynamic state, meshing mechanism, and contact characteristics of the walking wheel in order to improve the working performance of the shearer and coal mining efficiency.

Scholars from all over the world have researched a lot on the dynamic of the shearer and gear meshing mechanism. Chen et al. [5] established a simulation model of the shearer based on cutting load from an experimental drum. Mechanical characteristics of the walking mechanism were obtained when the pitch angle was −45° to 45° and the roll angle was −15° to 15°, based on considering clearance of the guide shoes. Through structural optimization, the service life of the guiding support plates is increased by 1.5 times. Xuan and Cui [6] obtained the tribological behavior of a sliding boot made of the Fe-Cr-B alloy under 30–60 N and 0.03–0.12 m/s. Reid et al. [7] proposed a real-time estimation method of shearer cutting force, and the result shows that the worst RMS error is 2.59%. Kinga et al. [8] studied the malfunction forms and causes of shearers, employing mathematical statistics, and identified stoppage reasons for 83.94% of the total breaks registered for the longwall shearer, which provided data basis for reliability analysis and fault detection of shearers. Mao et al. tested the walking characteristics of the shearer through the equal proportion fully mechanized mining experiment equipment and found that the longitudinal load of the left guide shoe increases by 884.3% under the oblique cutting condition [9]. They established the simulation model of coupling between the cutting unit and the traction unit and obtained the relationship between cutting force and traction speed [10]. In addition, they obtained the vertical vibration characteristics of shearers at different traction speeds and found that the biggest displacement of the drum is 38.9104 mm, by a six-degree-of-freedom dynamic model [11]. Di et al. [12] established the dynamic transfer model of the traction unit, analyzed its dynamic response under impact load, and proposed a reliability evaluation method based on the high-order stochastic response method. Zeng et al. [13] established the precise parametric models of the spur gear according to the principle of gear shaping, and the error of the driven gear rotation angle between the simulation result and the theoretical value is 8.85%, providing the reference for the precise design and optimization of gear transmission. Wang et al. [14] proposed a new system-structure coupling analysis method based on the consideration of both system excitation and structural parameters, which was used in the structural dynamic analysis of planetary gears, with obvious advantages over the traditional analysis model. Li et al. [15] established the dynamic equilibrium equations of the planetary gear train and found that the dangerous position is the root of the sun gear with the maximum stresses of 517.435 MPa. The result provides the research basis for the fatigue analysis of the gear. Santosh et al. [16] studied the meshing mechanism of helical gears under different friction conditions and analyzed the influence of the helical angle and the friction coefficient on the maximum stress. Results show that the stress increase in the spur gear is 10% and the increase in the helical gear is 8% when the friction coefficient increases from 0 to 0.3.

Existing literature on shearer’s traction unit mainly focuses on the study of the dynamic characteristics and reliability of the transmission system but less on the contact characteristics of the shearer’s walking wheel. In addition, the structure optimization of the walking wheel mostly refers to common gears at present, which has great limitations. Furthermore, the walking wheel-pin row belongs to open meshing transmission, under bad working conditions, and the popularization of large mining height shearers puts forward higher requirements for its reliability. Therefore, it is a great application value to study the meshing characteristics of the tooth pin under impact load.

The rest of this study is organized as follows. Section 2 establishes the dynamic state equations of the walking wheel and presents a theoretical model of the tooth pin contact force. Section 3 builds the rigid-flexible coupling simulation model of tooth pin meshing based on the model of the virtual prototype and gives the supposed conditions. Section 4 analyses stress distribution and longitudinal load characteristics of the walking wheel. Section 5 analyses the effects of load impact duration and impact amplitude on the meshing characteristics. Section 6 concludes the study.

#### 2. Dynamic Analysis of Walking Wheel

##### 2.1. Dynamic State Equation of Walking Wheel

The driving gear provides the torque input for the walking wheel, while the contact force between the walking wheel and the pin rows provides traction power for the shearer.

Based on the equivalent spring-damp model [17], the mechanical state model of the walking wheel is established, as is shown in Figure 1. And, the differential equation of the mechanical state of the walking wheel can be expressed aswhere *J*_{1} and *J*_{2} are the moments of inertia of the driving wheel and the walking wheel, respectively, *R*_{1} and *R*_{2} are the pitch radiuses of the driving wheel and the walking wheel, *T*_{1} and *T*_{2} are torques of the driving wheel and the walking wheel, *k*_{1} and *c*_{1} are the combined meshing stiffness and damping of the driving wheel, *k*_{2} and *c*_{2} are the combined meshing stiffness and damping of the walking wheel, and *θ*_{1} and *θ*_{2} are the rotation angles of the driving wheel and the walking wheel.

During the movement of the shearer, the center of mass of the walking wheel moves in line along pin rows while any one of its nodes rotates around the axis. As is shown in Figure 2, the motion state model of the walking wheel takes pin rows as the reference of the global coordinate system *O*-*xy* and the rotation center of the walking wheel as the reference of the local coordinate system [18] and marks *O*′ (*x*_{1}, *y*_{1}), node (*x*_{2}, *y*_{2}), node *Q* (*x*_{1}, 0), and initial angle *θ*_{o} between and *QO*′. Therefore, the equation of the motion state of the walking wheel can be expressed as

##### 2.2. Contact Force Model between Walking Wheel and Pin Rows

The meshing contact of the tooth pin is a highly nonlinear problem, and the force acting on the tooth surface of the walking wheel is very complicated. During a whole meshing period, contact force on the surface of the tooth includes a normal force *F*_{N} under the stable meshing state; transient impact force *F*_{i} due to the changes of node velocity, meshing-in, and meshing-out; and tangential friction force *F*_{f}. Therefore, the contact force model between the walking wheel and pin rows can be expressed as

On the basis of considering deformation of the gear, the normal contact force of the tooth pin is expressed aswhere *k* is the equivalent contact stiffness coefficient related to the material, *n* is the nonlinear impact force index, *c* is the damping [19], *r*_{1} and *r*_{2} are the curvature radiuses of the wheel and the pin at the contact point, respectively, *μ* is Poisson’s ratio of material, *E* is the elastic modulus of material, and *δ* is the amount of contact deformation [20].

In the process of gear transmission, the impact of nodes and meshing-out has little influence on the final transmission effect [21]. Therefore, this paper mainly considers the impact force of meshing-in of the walking wheel and uses the analysis method of reference [22] to get the impact velocity of the tooth pin at the moment of meshing-in as follows:where is the impact velocity at meshing-in, *ω* is the angular velocity of the walking wheel, *r*_{b} is the basic circle radius of the walking wheel, *α* is the nominal meshing angle, *γ* is the angle between the line from the meshing point to the rotation center and the vertical line [23], and *α*′ is the pressure angle at the meshing point.

The impact dynamic model of the walking wheel at meshing-in is shown in Figure 3. The mass of the shearer translation system can be converted into the equivalent rotational inertia around the revolving center of the walking wheel. Therefore, the kinetic energy generated by the conversion of the sum of rotational inertia into the equivalent mass of the unit tooth width on the meshing line is equivalent to the kinetic energy of meshing impact, as follows:where *E*_{k} is the kinetic energy of the meshing impact, *m*_{red} is the equivalent mass of the unit tooth width, *J* is the rotational inertia of the walking wheel, *J*_{t} is the equivalent rotational inertia of the shearer translation system, *b* is the tooth width, *m*_{t} is the comprehensive mass of the shearer translation system, is the radius of the actual base circle, and *r*_{b} is the radius of the theoretical base circle.

According to the theory of impact mechanics, maximum impact force, maximum deformation, and impact kinetic energy have the following relationships:where *δ*_{m} is the maximum deformation of the walking wheel at the moment of meshing-in and *q*_{s} is the comprehensive flexibility of the walking wheel at the meshing point.

According to formulas (8)–(10), the meshing impact force of the walking wheel can be expressed as follows:

In order to fully reflect the sliding friction between the teeth, the tangential friction on the tooth surface can be expressed by the Coulomb friction model as follows:where *μ* is the sliding friction coefficient, *F*_{n} is the normal force of the contact point at any meshing time, and is the tangential velocity of the meshing point.

#### 3. Modeling and Simulation

##### 3.1. Model Simplification

As is shown in Figure 4, a simplified virtual prototype model of tooth pin meshing transmission is established. The tooth profile of the walking wheel is involute. The number of the teeth is 11, and the modulus is 46.8 mm. In order to prevent undercutting, appropriate modification has been made. The pin tooth is the involute profile with a pitch of 147 mm. Only considering the main functional structural features, chamfers and other auxiliary structural features are idealized. Finally, the meshing system of the tooth pin with the positive transmission is constructed according to the requirements of no backlash center distance and standard top clearance.

##### 3.2. System Simulation

The three-dimensional model is divided into the hexahedron grid, with the walking wheel set as flexible and the inner surface of the axle hole set as rigid for the convenience of adding the simulation drive. The pin row is set as rigid because of actual constraints and relatively small deformation. In addition, the density, elastic modulus, and Poisson’s ratio of the wheel and the pin row are defined, respectively, according to the material properties of 18Cr_{2}Ni_{4}MA and 40CrMNMo. Furthermore, the relative motion relationship between the walking wheel and the pin row is equivalently treated in order to facilitate simulation setting. Therefore, a rotating pair around the *Z*-axis is added to the walking wheel, and a moving pair along the *X*-axis is added to the pin row to simulate the relative motion between the pin row and the walking wheel. Moreover, because the influences of the swing angle of shearer’s rocker arm and coal seam inclination are not considered, traction load is equivalent to the resistance opposite to the direction of movement, and the rated load is set as 380 kN. The simulation system of tooth pin meshing transmission is driven by angular speed of the walking wheel. The angular speed is 0.647 rad/s, obtained by shearer’s traction speed of 10 m/min. Finally, the contact type is set as automatic surface to surface and solved by the LS-DYNA solver in ANSYS.

#### 4. Analysis of Load Distribution

##### 4.1. Stress State Analysis of Walking Wheel

The equivalent stress state of the walking wheel at the moment of the single tooth contact is shown in Figure 5. This simulation results show that there is obvious concentrated stress on both sides of the tooth root and the contact area. It is consistent with the actual easy wear and fracture position of the walking wheel and with the results of existing modeling [24]. This result reflects the feasibility of this simulation method to a certain extent. The stress distribution of the tensile side and the compressed side of the tooth root is asymmetrical at the moment of meshing. The phenomenon is caused by the obvious difference of equivalent moment between two sides of the tooth root, according to the calculation model of gear bending stress based on the cantilever beam hypothesis [25].

The stress distribution in the contact area is band-shaped along the *Z* direction and narrower in the middle. It indicates that the tooth pin meshing is not an ideal line contact, which is caused by the elastic deformation of the contact position. It provides a reference for further study of the contact state of the walking wheel. In addition, stress on the contact end surface of the walking wheel is larger and the concentrated stress is more obvious, similar to the result of reference [26]. It is due to the “edge effect” in the theory of finite length linear contact mechanics [27]. It shows that there is a complex material deformation on the end face according to the classical Hertzian contact theory. Conversely, the characteristic of bending stress on the tensile side of the tooth root is large in the middle and small on the end face, which is related to the alternating tension and compression at two sides of the tooth root. Furthermore, the stress of the walking wheel shows good symmetry in the direction of the tooth width, which is different from the result of reference [28]. This is because the influence of the rotating shaft is neglected in this modeling in order to emphasize the difference of longitudinal load caused by the structural deformation of the walking wheel; however, in order to obtain the mechanical state of different stages, the contact force in reference [28] is added in different regions.

The time-varying stress in the middle and the end of the tooth width of the walking wheel is extracted as shown in Figure 6, by choosing the position of two sides of the tooth root, the pitch line, and the front and the back of the pitch line. The stress curves of adjacent teeth show good periodicity according to Figure 6(a). The meshing frequency of the single tooth of the walking wheel is 1.06 Hz, basically consistent with the theoretical value of 1.14 Hz. The results show that there are phenomena of premature meshing-in or delayed meshing-out, in the meshing process, due to elastic deformation. The result is consistent with the assumption that the meshing position shifts to the tooth root in the impact dynamics model in Figure 3. In addition, in the position of the pitch line of the walking wheel, the stress of the end face increases sharply at the moment of contact, while the stress in the middle of the tooth width decreases slightly. It indicates that the impact of nodes on the end face has a greater influence than that in the middle of the tooth width. Figure 6(b) is the maximum stress curve in the contact area and the tooth root of the wheel. The maximum stress on both sides of the tooth root changes little, with a maximum of 334 MPa at the moment of meshing-out. In the contact area, the maximum is 814 MPa at the beginning of meshing, indicating that the impact of meshing-in has the most obvious effect on the final contact characteristics.

**(a)**

**(b)**

**(c)**

**(d)**

The variation curve of equivalent stress and its component force of the walking wheel near the pitch line is shown in Figure 6(c). The results show that the stress duration is more than 0.38 s, and the maximum stress is 810 MPa, in the front of the pitch line; and in the position of the theoretical pitch line, the maximum is reduced to 670 MPa. This indicates that the pitch line has the tendency of directional tooth root deviation in the actual meshing process, which is aligned with the analysis in reference [23]. In addition, the variation of the *Y*-direction component of contact stress lags behind the *X*-direction component by 0.8 s, and the amplitude of the *Y*-direction component is significantly smaller than that of the *X*-direction component. Moreover, the variation trend of the *X*-direction component near the pitch line is basically the same as that of the total stress, and the peak value remains basically stable; the *Y*-direction component begins to decay rapidly when the theoretical pitch line enters into meshing. The reason is that the angle between the contact surface and the vertical plane is close to zero when the theoretical pitch line participates in meshing. The small change rate of the cosine function of the angle makes the change of the *X*-direction component smaller, while the large change rate of the sinusoidal function makes the change of the *Y*-direction component more obvious.

Figure 6(d) is the variation curve of equivalent stress and its component on both sides of the tooth root of the walking wheel. The results show that the stress on both sides of the root tends to increase in a single cycle, with a maximum of 323 MPa on the compressed side. This is because the contact position moves from the root to the top of the teeth during the meshing process of the walking wheel, which leads to the increase in torque at the root position generated by the contact force. In addition, stress of the tensile side sharply increases by more than 80 MPa when the pitch line enters meshing, which indicates that the tensile side is more affected by the impact of nodes than the compressed side.

##### 4.2. Longitudinal Load Analysis of Walking Wheel

Longitudinal load distribution is an important reference factor in the gear design. Thus, it is necessary to analyze the degree of longitudinal load imbalance in different areas of the walking wheel. According to the principle of the gear design, the ratio of the maximum stress to the average stress in the same tooth width direction is defined as the longitudinal load distribution coefficient (*K*_{β}), and the simulation results are shown in Figure 7.

**(a)**

**(b)**

**(c)**

The *K*_{β} fitting curve of the contact line position shows W-shape, as is shown in Figure 7(a). The peak time is basically consistent with the time of meshing-in impact, nodes impact, and meshing-out impact, respectively. In addition, the increase in *K*_{β} is more obvious in the midmeshing period with the maximum of 1.325 because it is single-tooth meshing in this period; the longitudinal load distribution tends to be uniform after the pitch line enters into meshing because the next tooth starts to engage in meshing at this stage, the walking wheel is double-teeth meshing, and its mechanical state is relatively stable.

The stress of transition position from the root to the abdomen of the tooth is complex. Thus, *K*_{β} of the starting point, the midpoint, and the terminal point of the transition curve from the root to the abdomen are extracted in turn. As is shown in Figure 7(b), at the starting and the middle of the tensile side transition curve, the *K*_{β} curve fluctuates gently, with the value close to 1. This shows that the longitudinal load distribution is relatively uniform, and the impact of meshing is relatively low. In addition, at the end of the tensile side transition curve, *K*_{β} shows a downward trend as a whole; and the overall distribution of the longitudinal load is uniform after the pitch line engages in contact. This indicates that in the early stage of meshing, longitudinal unbalance load is more serious due to the meshing-in impact.

As is shown in Figure 7(c), in the compressed side, the overall distribution of the longitudinal load is more uniform. At the starting point of the transition curve, the *K*_{β} curve tends to increase after the pitch line engages in contact; at the terminal point of the transition curve, *K*_{β} tends to decrease. However, the changes at both locations are smaller relative to the tension side. Therefore, the analysis shows that compared with the compressed side, the stress on the transition position between the tooth root and the abdomen of the tensile side is more complex, which should be put into consideration in the design process.

The stress distribution and unbalance load of the walking wheel are closely related to its longitudinal deformation. Therefore, the displacement of each element in the same tooth width position is extracted along the *Z*-direction, and the maximum displacement difference (∆*s*) in the *Z*-direction in different areas of the walking wheel is obtained, as is shown in Figure 8. ∆*s* in the contact area of the gear abdomen is obviously larger than that in the bending area of the root. The peak of ∆*s* is 27.38 *μ*m near the pitch line and deviates to the tooth root, and its time is ahead of the peak time of unbalance load by 0.2 s. Along the tooth profile from the root to the top, the peak of ∆*s* tends to decrease. This indicates that the serious unbalance load which is near the pitch line and deviates to the tooth root is related to the *Z*-direction displacement difference of this position, which provides the reference for the further study of the actual meshing state of the walking wheel.

#### 5. Analysis of Meshing Characteristics under Impact Load

##### 5.1. The Influence of Impact Load Amplitude on Meshing Characteristics

In this simulation, 10%, 20%, 30%, and 40% extra impact loads are added at 0.5 s based on the rated load of 380 kN, and the duration of impact load is 0.01 s. The maximum equivalent stress of the walking wheel in the contact area of the tooth abdomen and in the bending area of the tooth root is extracted and analyzed on nonlinear regression, as is shown in Figure 9. The results show that the instantaneous impact amplitude has a significant effect on the maximum equivalent stress in the contact area of the tooth abdomen, but has a little effect on the tooth root. According to the regression analysis, the maximum equivalent stress under different impact magnitudes can be predicted. The regression equation shows that the material in the contact area of the tooth abdomen will fail when the walking wheel bears 39% extra instantaneous impact load.

In order to study the failure mechanism of the walking wheel overload fracture, *K*_{β} of the tooth root tension side under different amplitudes of instantaneous impact loads is extracted, as is shown in Figure 10. The results show that *K*_{β} tends to increase under instantaneous load impact, relative to the ideal condition; however, within 30% of the overload range, the larger the overload amplitude, the smaller the increase in *K*_{β}. This is because the overall increase in contact stress results in the weakening of the longitudinal load imbalance, within 30%. And, when the overload is 40%, some materials have exceeded the allowable limit and thus failure occurs. In addition, the effect of impact load on *K*_{β} disappears at 0.54 s, and the disappearance time lags behind the impact time by 0.3 s. Based on the Kelvin–Voigt viscoelastic model, the phenomenon is related to the lag of material deformation recovery to stress disappearance. The results provide the reference for optimizing the load coefficient in checking the bending fatigue strength of the walking wheel under impact load and for preventing overload breakage of running wheels after impact load occurs.

In order to study the influence of the impact amplitude on the longitudinal deformations of the walking wheel, the change of maximum displacement difference in the *Z*-direction (∆*s*) near the pitch line with load impact is extracted to compare with that under ideal conditions, as is shown in Figure 11. The results show that ∆*s* on the pitch line increases at the moment of impact, and the phenomenon becomes more obvious with the increase in the impact amplitude. This is because the position is not actually involved in contact at this time, and the *Z*-direction deformation fluctuation occurs when this position is extruded by the front material element. However, at the location near the pitch line and inclining to the root side, ∆*s* decreases obviously when the impact occurs, and the influence of the impact on this position lasts longer than that on the pitch line. This is because the meshing contact has occurred at this position, and the overall increase in the *Z*-direction deformation leads to the decrease in ∆*s*. And, the periodic fluctuation between the deformation and recovery of elastic materials results in the impact load having a more lasting influence on this position.

##### 5.2. The Influence of Impact Load Duration on Meshing Characteristics

The rated load is set at 380 kN, 30% impact load is added at 0.5 second, and the duration is 0.2 s, 0.3 s, 0.4 s, and 0.5 s, respectively. The maximum equivalent stress fitting curves of different areas of the walking wheel are obtained, as is shown in Figure 12. The result shows that the maximum equivalent stress decreases slightly with the increase in the impact duration under the same impact amplitude. The reason is that the step function is used to simulate the impact load, and the load change rate decreases with the increase in impact duration. It is revealed that under the same impact load amplitude, the appropriate extension of impact duration has little effect on the maximum stress of the walking wheel and will not significantly aggravate the pitting and wear failure of the local position.

Under the same impact amplitude and diffident impact duration, *K*_{β} of the tensile side of the root is obtained, as is shown in Figure 13. The result shows that *K*_{β} increases slightly at the beginning of the impact load, with no obvious difference among the four curves at other times. The maximum displacement difference in the *Z*-direction (∆*s*) near the pitch line of the walking wheel with different impact load duration is shown in Figure 14. The result shows that *K*_{β} increases with the increase in load impact duration, and the increment on the pitch line is significantly higher than that before the pitch line. It is indicated that the duration of the impact load has a great influence on the impact of nodes, resulting in the larger longitudinal deformation on the pitch line of the walking wheel.

#### 6. Conclusions

To reduce the failure rate and improve the working performance of the walking wheel of the shearer, the dynamic state equations and the contact force theoretical model of the walking wheel are established; the rigid-flexible coupling simulation model of tooth pin meshing is established; the load distribution characteristics of the walking wheel and its response characteristics to impact load are analyzed. The following conclusions are drawn:(1)The stress distribution in the contact area of the walking wheel is small in the middle of the tooth width and large on the edge, with the maximum of 814 MPa at the moment of meshing-in; the end part of the contact face is more seriously affected by the meshing impact, and the actual meshing points tend to move towards the tooth root. Furthermore, in the meshing period of one tooth, the equivalent stress of the tooth root tends to increase as a whole and sharply increases by more than 80 MPa at the moment of the single tooth contact, especially at the tension side. It provides guidance for the design optimization of the walking wheel.(2)In the contact area, the longitudinal unbalance load is the most serious at the time of meshing-in, single-tooth meshing, and meshing-out, and the curve of *K*_{β} is W-shaped, with the maximum of 1.325 at the moment of the single tooth contact. In addition, on the tensile side, at the end of the transition curve of the tooth root, the imbalance of the longitudinal load is more serious before the pitch line enters the meshing. Furthermore, the distribution of longitudinal loads on the compressed side is relatively uniform.(3)With the amplitude of the impact load increasing, the maximum stress increases, *K*_{β} on the tensile side of the root tends to decrease, and the maximum displacement difference in the *Z*-direction on the pitch line increases. In addition, the duration of short-term load impact has little effect on the meshing characteristics of the walking wheel.

In the study, rigid-flexible coupling simulation analysis of the meshing characteristics between the shearer’s walking wheel and the pin row is carried out. It provides guidance for studying the failure mechanism of the walking wheel and improving its reliability. And, it provides some reference for the study of meshing contact characteristics with a special structure. In future research, it will be necessary to consider the influence of the cutting part on the traction load and optimize the load addition mode. And, it will be necessary to further explore the coupling effect of pin deformation on meshing characteristics.

#### Data Availability

The data used to support the findings of this study 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

This work was supported by the Innovative Team Development Project of the Ministry of Education (Grant no. IRT_16R45), Key Research and Development of Shandong Province (Grant no. 2018GGX103002), and Postgraduate Science and Technology Innovation Project of Shandong University of Science and Technology (Grant no. SDKDYC190328).