Analysis of Dynamic Mesh Stiffness and Dynamic Response of Helical Gear Based on Sparse Polynomial Chaos Expansion

: This paper presents an efﬁcient method for obtaining the dynamic mesh stiffness and dynamic response of a helical gear pair. Unlike the traditional dynamic model that utilizes a time-dependent sequence, the mesh stiffness using the presented method is updated according to the gear displacement vector at each sub-step of the numerical calculation. Three-dimensional loaded tooth contact analysis (3D LTCA) is used to determine the mesh stiffness, and a surrogate model based on sparse polynomial chaos expansion (SPCE) is proposed to improve the computational efﬁciency, which is achieved by reducing the number of coefﬁcients in the polynomial chaos expansion (PCE) model though a quantum genetic algorithm. During the calculation, the gear displacement vector at each sub-step is converted into the changes in center distance, misalignment angle, and mesh force, which are then introduced into the SPCE model to update the mesh stiffness for subsequent calculations. The results suggest that the SPCE model exhibits high accuracy and can signiﬁcantly improve the computational efﬁciency of the PCE model, making it suitable for dynamic calculations. Upon updating the mesh stiffness during the dynamic calculation, the mesh stiffness declines, the dynamic transmission error (DTE) increases, and the frequency components of the responses change signiﬁcantly.


Introduction
Gear systems are widely used in engineering, and mesh stiffness is a main source of internal excitation in gear systems. Over the past few decades, many scholars have focused on researching the mesh stiffness of a gear pair. The common methods for calculating time-varying mesh stiffness (TVMS) include the analytical method, finite element (FE) method, and analytical-FE method. The potential energy method [1][2][3][4][5][6] is one of the most widely used analytical methods, and it was expanded from analyzing spur gears to helical gears using the slice method [7][8][9][10][11][12], which divides helical gears into spur gear slices with small widths. To improve the accuracy of the calculation, it is important to consider the coupling effect between the slices in the slice method. Tang et al. [10] considered the deformation of the teeth and the bodies of the slices in the coupling model separately, thus obtaining more accurate calculation results. Yu et al. [11] utilized nonuniformly distributed weighting coefficients in parabolic form to simulate the coupling between the slices. Wang et al. [13] established a slice coupling model based on the Wright formula [14]. The deformation transfer formulas in both contact and non-contact states were presented, and the TVMS of spur gears with misalignment and modification were identified. Jordan et al. [12] combined the exponential decay function and the moment images method to calculate the forces and deformations of the slices and then obtained the mesh stiffness of a helical gear pair. Using the analytical method, scholars have analyzed the effects of Tayari [46] established a torsional dynamic model of a multistage helical gear system and explored the nonlinear behavior of the system by introducing TVMS and backlash. Zhang [47] analyzed the influence of modification on the vibration characteristics of a pure electric vehicle gearbox. The results showed that the modification of the small sun wheel can reduce the vibration and noise generated by gear meshing. Kumar [48] presented the formulation of a carburized spur gear pair for the calculation of TVMS. The results showed that the carburizing makes the TVMS higher and the response amplitude lower. Yuan [49] calculated dynamic TVMS and dynamic load distributions using a combination model of LTCA and Newmark integration considering the influence of tooth flank errors. Wang et al. [50] analyzed the effect of bearing misalignments on the vibration characteristic of a gear system. In addition to substituting TVMS into the dynamic model, the general contact algorithm is also widely used to analyze the vibration generated by mesh excitation. Lundvall [51] established a multibody model by superposing small displacement elasticity obtained using the FE method on the rigid body motions equation and analyzed the influence of friction on the dynamic transfer error. Ooi [52] established an FE model of a portal axle gear system, calculated the modal parameters of the system, and analyzed the stress distribution with different angular positions. Wu [53] calculated the dynamic TVMS and dynamic transmission error (DTE) of a helical gear pair utilizing the three-dimensional contact FE method. Tamarozzi [54] analyzed the importance of the mode selection of residual attachment modes for the accuracy of multi-body dynamics calculation. Liu [55] utilized the low-frequency modal approximation method to reduce the degrees of freedom and simplified the Jacobian matrix by ignoring the effect of inertial forces. Their method has been demonstrated to be highly efficient via correlation with commercial software.
As mentioned above, most studies use the quasi-static method to calculate TVMS, and only the deformation of a gear pair is considered, and the influences of system deformation and gear vibration are often ignored. Although the contact algorithm based on FEM takes these factors into account, it has the problems of low efficiency and limited accuracy with mesh density. In this paper, an efficient method for the analysis of dynamic TVMS and dynamic response is proposed. It establishes a relationship between analytical finite element 3D LTCA and the dynamic model through TVMS. In 3D LTCA, the changes in the center distance, misalignment angle, and mesh force caused by meshing excitation are considered and reflected in the TVMS. Since the efficiency of analytical finite element 3D LTCA cannot meet the requirements of the iteration, SPCE technology is introduced to generate a surrogate model to efficiently obtain the TVMS.

Sparse polynomial Chaos Expansion (SPCE)
Polynomial chaos expansion (PCE) is a modeling method widely used in the study of sensitivity and uncertainty analysis. It uses a set of orthogonal polynomials to fit the output of random variables, which can be expressed as: where a denotes the polynomial coefficients, ω is the random variable, and I n is the orthogonal polynomial of order n. The most commonly used polynomials are Hermite polynomials, which are associated with Gaussian distribution and are suitable for most cases. However, their accuracy may not be adequate for complex conditions. As a solution, mixed polynomials are often used to enhance the model's expressive capability. In this paper, two types of orthogonal polynomials are used, namely, Hermite and Legendre polynomials. Figure 1 shows the first five orders of the polynomials.
By substituting random variables and their response Y into Equation (1), the surrogate model can be transformed into a linear regression problem, and its solution calculated using the least-squares method can be expressed as: where â is the coefficient vector, and Ψ + is the generalized inverse of the polynomial matrix Ψ. In practical applications, the accuracy of the PCE model is greatly influenced by the maximum polynomial order. The relationship between the number of coefficients P and the maximum polynomial order p and the number of factors M can be expressed as: It can be seen that as M and p increase, P significantly increases. However, most of the coefficients are related to high-order polynomials and cross-terms that have a negligible effect on the accuracy of the surrogate model but hinder computational efficiency. Therefore, a sparsification method is proposed based on a quantum genetic algorithm to reduce the number of coefficients. Compared with traditional genetic algorithms, quantum algorithms use qubit codes instead of binary gene codes and quantum rotation gates instead of crossovers. This approach offers fast convergence and high stability.
In a quantum genetic algorithm, information is stored in the form of qubits and is in the superposition of two quantum states at the same time: where 0 and 1 are spin-down and spin-up states, respectively, and (α, β) are two complex constants, which represent the probability amplitudes of the two states and satisfy 22 1  +=

(7)
Using two qubit states to indicate whether coefficients are considered, the sparsification can be expressed as a qubit code: The Hermite polynomial can be expressed as: The Legendre polynomial can be expressed as: By substituting random variables and their response Y into Equation (1), the surrogate model can be transformed into a linear regression problem, and its solution calculated using the least-squares method can be expressed as: whereâ is the coefficient vector, and Ψ + is the generalized inverse of the polynomial matrix Ψ.
In practical applications, the accuracy of the PCE model is greatly influenced by the maximum polynomial order. The relationship between the number of coefficients P and the maximum polynomial order p and the number of factors M can be expressed as: It can be seen that as M and p increase, P significantly increases. However, most of the coefficients are related to high-order polynomials and cross-terms that have a negligible effect on the accuracy of the surrogate model but hinder computational efficiency. Therefore, a sparsification method is proposed based on a quantum genetic algorithm to reduce the number of coefficients. Compared with traditional genetic algorithms, quantum algorithms use qubit codes instead of binary gene codes and quantum rotation gates instead of crossovers. This approach offers fast convergence and high stability.
In a quantum genetic algorithm, information is stored in the form of qubits and is in the superposition of two quantum states at the same time: where | 0 and | 1 are spin-down and spin-up states, respectively, and (α, β) are two complex constants, which represent the probability amplitudes of the two states and satisfy Machines 2023, 11, 736 5 of 20 Using two qubit states to indicate whether coefficients are considered, the sparsification can be expressed as a qubit code: The quantum rotation gate can be expressed as: where θ qi is the rotation angle in generation i, which indicates the direction of evolution. Figure 2 illustrates the flow of the proposed sparsification method. Firstly, a set of qubit codes are initialized. In order to make the probability of the two states equal, the initial qubit codes of all coefficients are: Then, by observing each qubit code, a series of sparse coefficient vectors is obtained, and PCE models are developed accordingly. Additionally, R 2 (the determination coefficient, R-squared) and RMSE (root-mean-square error) are calculated from the test set to evaluate the accuracy of the models, which can be expressed as:  (12) where ŷ is the estimated value, y is the observed value, and y is the average of the observed value. R 2 is within the interval [0, 1], and the closer it is to 1, the more accurate the model, and the smaller the RMSE is, the more accurate the model. The best result is sought from these observations, which can be represented as an optimization problem: where Ri 2 and RMSEi are the R 2 and RMSE of observation i, and R0 2 and RMSE0 are the initial R 2 and RMSE before sparsification. Finally, all qubit codes pass through the quantum rotation gate and evolve toward the best result. Then, new quantum bits are generated, forming a cycle to support the iterative solution.  Figure 3 illustrates the entire modeling process of the SPCE model. Initially, the training set is generated using the uniform design method to ensure the representativeness of Then, by observing each qubit code, a series of sparse coefficient vectors is obtained, and PCE models are developed accordingly. Additionally, R 2 (the determination coefficient, R-squared) and RMSE (root-mean-square error) are calculated from the test set to evaluate the accuracy of the models, which can be expressed as: whereŷ is the estimated value, y is the observed value, and y is the average of the observed value. R 2 is within the interval [0, 1], and the closer it is to 1, the more accurate the model, and the smaller the RMSE is, the more accurate the model. The best result is sought from these observations, which can be represented as an optimization problem: where R i 2 and RMSE i are the R 2 and RMSE of observation i, and R 0 2 and RMSE 0 are the initial R 2 and RMSE before sparsification. Finally, all qubit codes pass through the quantum Machines 2023, 11, 736 6 of 20 rotation gate and evolve toward the best result. Then, new quantum bits are generated, forming a cycle to support the iterative solution. Figure 3 illustrates the entire modeling process of the SPCE model. Initially, the training set is generated using the uniform design method to ensure the representativeness of the samples. Next, a polynomial matrix is constructed based on Hermite and Legendre polynomials, and the coefficient vectors are calculated to develop the PCE surrogate model. Finally, a sparsification method is introduced to reduce the number of coefficients while retaining the accuracy of the model, which is verified with the R 2 and RMSE of the testing set. In contrast with the training set, the sample points of the testing set are selected randomly.

Loaded Tooth Contact Analysis
In this paper, the mesh stiffness of a helical gear pair is calculated using 3D LTCA. The general deformation compatibility principle leads to the following equations: where λb is the global compliance matrix, λc is the contact compliance matrix, F is the load distribution vector, ste is the static transmission error, E is the clearance distribution vector, f is the mesh force, km is the mesh stiffness, and nlste is the non-loaded static transmission error, which can be calculated via nlste = min(E). λc is calculated using the nonlinear Hertzian contact formula: where E is the elastic modulus, Fi is the elements of F, Li is the face width of the i-th segment, and n is the number of contact points. λb can be expressed as: where subscripts "p" and "g" denote the pinion and gear, and i and j denote two possible contact nodes on the tooth surface. In conformity with the definition of compliance, λij can be seen as the displacement of point i when a unit force is applied to point j. Herein, the unit force method based on the finite element model is employed to obtain λb, and the iterative method is employed to solve Equation (14). (Interested readers may refer to Refs. [38,39]).

Loaded Tooth Contact Analysis
In this paper, the mesh stiffness of a helical gear pair is calculated using 3D LTCA. The general deformation compatibility principle leads to the following equations: k m = f ste − nlste (15) where λ b is the global compliance matrix, λ c is the contact compliance matrix, F is the load distribution vector, ste is the static transmission error, E is the clearance distribution vector, f is the mesh force, k m is the mesh stiffness, and nlste is the non-loaded static transmission error, which can be calculated via nlste = min(E). λ c is calculated using the nonlinear Hertzian contact formula: where E is the elastic modulus, F i is the elements of F, L i is the face width of the i-th segment, and n is the number of contact points. λ b can be expressed as: where subscripts "p" and "g" denote the pinion and gear, and i and j denote two possible contact nodes on the tooth surface. In conformity with the definition of compliance, λ ij can be seen as the displacement of point i when a unit force is applied to point j. Herein, the unit force method based on the finite element model is employed to obtain λ b , and the iterative method is employed to solve Equation (14). (Interested readers may refer to Refs. [38,39]).

Influencing Factors
The influences of three factors on mesh stiffness are now considered, which are mesh force, angular misalignment, and center distance. Equation (14) directly introduces the influence of mesh force into 3D LTCA, so it is not discussed further.
A change in center distance can lead to a variation in the working pressure angle, as shown in Figure 4, which can be expressed as: where α is the working pressure angle, and α and a 0 are the initial working pressure angle and initial center distance. Due to the change in the working pressure angle, the transverse contact ratio also changes, which can be expressed as: where z is the tooth number, α a is the addendum pressure angle, and the subscripts "p" and "g" represent the pinion and gear, respectively.

Influencing Factors
The influences of three factors on mesh stiffness are now considered, which are mesh force, angular misalignment, and center distance. Equation (14) directly introduces the influence of mesh force into 3D LTCA, so it is not discussed further.
A change in center distance can lead to a variation in the working pressure angle, as shown in Figure 4, which can be expressed as: where α′ is the working pressure angle, and α and a0 are the initial working pressure angle and initial center distance. Due to the change in the working pressure angle, the transverse contact ratio also changes, which can be expressed as: where z is the tooth number, αa is the addendum pressure angle, and the subscripts "p" and "g" represent the pinion and gear, respectively. The misalignment angle can be divided into two parts, on the plane of action and on the off plane of action. According to Refs. [23,56], the misalignment angle on the plane of action has a great influence on the mesh stiffness, while the other can be ignored. The effect of the misalignment angle is introduced into LTCA via the mesh clearance, which can be expressed as: where ei is the mesh clearance, bi is the coordinates along the tooth width, θ is the misalignment angle on the plane of action, and β is the helix angle. According to the projection of the space angle, θ can be calculated from the angular displacement θx and θy (see Figure  5b where ϕ is the angle between the action panel and the x-axis: The misalignment angle can be divided into two parts, on the plane of action and on the off plane of action. According to Refs. [23,56], the misalignment angle on the plane of action has a great influence on the mesh stiffness, while the other can be ignored. The effect of the misalignment angle is introduced into LTCA via the mesh clearance, which can be expressed as: where e i is the mesh clearance, b i is the coordinates along the tooth width, θ is the misalignment angle on the plane of action, and β is the helix angle. According to the projection of the space angle, θ can be calculated from the angular displacement θ x and θ y (see Figure 5b,c): where φ is the angle between the action panel and the x-axis:

Dynamic Model
In this paper, the dynamic model consisting of gear-shaft elements is utilized, as depicted in Figure 6, wherein both the helical gear and shaft are simulated with the Timoshenko beam element. The helical gear is composed of two beam elements with the radius of a reference circle, and the two intermediate nodes of the gears are connected with a spring-damping system, which represents the meshing of the gear pair.

Dynamic Model
In this paper, the dynamic model consisting of gear-shaft elements is utilized, as depicted in Figure 6, wherein both the helical gear and shaft are simulated with the Timoshenko beam element. The helical gear is composed of two beam elements with the radius of a reference circle, and the two intermediate nodes of the gears are connected with a spring-damping system, which represents the meshing of the gear pair.

Dynamic Model
In this paper, the dynamic model consisting of gear-shaft elements is utilized, as depicted in Figure 6, wherein both the helical gear and shaft are simulated with the Timoshenko beam element. The helical gear is composed of two beam elements with the radius of a reference circle, and the two intermediate nodes of the gears are connected with a spring-damping system, which represents the meshing of the gear pair. The mass matrix of the Timoshenko beam element can be expressed as: Machines 2023, 11, 736 9 of 20 a = 13 35 + 7 10 ϕ + 1 3 ϕ 2 + 6 5 (r g /l) 2 (1 + ϕ) 2 , b = 9 70 + 3 10 ϕ + 1 6 ϕ 2 − 6 5 (r g /l) 2 (1 + ϕ) 2 c = 11 210 + 11 120 ϕ+ 1 24 where A s is the sectional area, l is the element length, I is the cross-sectional moment of inertia, J is the torsional moment of inertia, and ϕ is the shear influence coefficient. The stiffness matrix of the Timoshenko beam element can be expressed as: The coordinate of the helical gear pair q can be written as: q = x p , y p , z p , θ xp , θ yp , θ zp , x g , y g , z g , θ xg , θ yg , θ zg T (26) Subsequently, the stiffness matrix of the mesh model can be expressed as: where K m is the 12 × 12 stiffness matrix of the mesh element, and v m is the projection vector along the mesh line, which can be expressed as: v m = [sin Ψ cos β, cos Ψ cos β, sin β, -r p sin Ψ sin β, -r p cos Ψ sin β, r p cos β, − sin Ψ cos β, − cos Ψ cos β, − sin β, -r g sin Ψ sin β, -r g cos Ψ sin β, r g cos β] where r is the base radius. The dynamic transfer error can be written as: The dynamic mesh force can be written as: where c m is the mesh damping. Rayleigh damping is adopted here, which can be expressed as: where α c and β c are the damping coefficients, and c m can be calculated using: The dynamic mesh stiffness can be written as: The proposed method for updating mesh stiffness is depicted in Figure 7. The dynamic response of the system is computed via Newmark iteration. The mesh force, angular misalignment, and center distance are extracted from each sub-step of the Newmark iteration and then input into the SPCE model to obtain the mesh stiffness. The DTE, dynamic mesh force, and dynamic mesh stiffness of the gear pair are calculated using Equations (29), (30), and (33).

Numerical Examples
The dynamic responses of the proposed method and traditional method were compared through numerical examples, as shown in Figure 8. The parameters of the gear pair and the shaft are shown in Tables 1 and 2. In addition, the elastic modulus is 210 GPa, the density is 7850 kg/m 3 , Poisson's ratio is 0.3, the speed is 1200 rpm, and the torque is 500 Nm.
l Bearing Figure 7. Proposed method for updating mesh stiffness.

Numerical Examples
The dynamic responses of the proposed method and traditional method were compared through numerical examples, as shown in Figure 8. The parameters of the gear pair and the shaft are shown in Tables 1 and 2. In addition, the elastic modulus is 210 GPa, the density is 7850 kg/m 3 , Poisson's ratio is 0.3, the speed is 1200 rpm, and the torque is 500 Nm.

Numerical Examples
The dynamic responses of the proposed method and traditional method were compared through numerical examples, as shown in Figure 8. The parameters of the gear pair and the shaft are shown in Tables 1 and 2. In addition, the elastic modulus is 210 GPa, the density is 7850 kg/m 3 , Poisson's ratio is 0.3, the speed is 1200 rpm, and the torque is 500 Nm.   The SPCE model for mesh stiffness requires four inputs: mesh cycle, mesh force, angular misalignment, and center distance. Among them, the relationship between the mesh cycle and mesh stiffness is the most complex. In order to determine the maximal polynomial order for the mesh cycle, PCE models with varying maximal polynomial orders were generated and compared, in which the mesh cycle was the only input parameter. RMSE and R 2 were then employed as accuracy indices to evaluate the accuracy of the generated models (see Figure 9). It should be noted that the PCE model contained two sub-models because of the change in the number of mesh teeth, and the polynomial order of the two sub-models is the same. The results indicate that the accuracy of the PCE model improved with an increase in the polynomial order. In order to ensure both computational efficiency and accuracy, the maximal polynomial order of the mesh cycle was set to 8. At this order, the value of R 2 was 0.9994, and the RMSE was 5.75 × 10 5 , indicating the high level of accuracy of the PCE model.
If the maximal polynomial order of each input parameter is 8, the coefficient vector will be too numerous to be calculated. To enhance efficiency, the maximal polynomial order of the rest input parameters was set to 3. The samples of the training set were selected based on the uniform design, and the parameter intervals are listed in Table 3. As a result, 21 groups of mesh stiffness with different input parameters were calculated, and each group contained the mesh stiffness in an entire mesh cycle.
If the maximal polynomial order of each input parameter is 8, the coefficient vector will be too numerous to be calculated. To enhance efficiency, the maximal polynomial order of the rest input parameters was set to 3. The samples of the training set were selected based on the uniform design, and the parameter intervals are listed in Table 3. As a result, 21 groups of mesh stiffness with different input parameters were calculated, and each group contained the mesh stiffness in an entire mesh cycle.  The PCE model for mesh stiffness was generated using the training set, and the influence of each factor is presented in Figure 10a-c. To assess the accuracy of the model, test points were randomly selected on the response surface, and the accuracy indexes were calculated as shown in Figure 10d. It can be seen that mesh stiffness increased with increasing load but decreased with an increasing misalignment angle and center distance. Furthermore, the accuracy indexes indicate that the model accurately reflects the factors that influence mesh stiffness. The proposed sparsification method was then applied to the PCE model, and 21 groups of mesh stiffness with random input parameters were calculated and employed as the testing set. As shown in Figure 11, the coefficients of the PCE model decreased gradually with the increase in iterations. During the sparsification process, R 2 remained almost unchanged, and the RMSE decreased slightly, which suggests that the sparse model maintains the accuracy of the original model.
Utilizing the proposed method (depicted in Figure 7), the dynamic mesh stiffness, DTE, and dynamic mesh force were extracted from the results of different cases, as illustrated in Figure 12. The average computation time using the SPCE model was 4.28 s, while the average computation time using the PCE model was 19.27 s. This indicates that the SPCE model significantly improves computational efficiency. Additionally, it is evident that the dynamic mesh stiffness obtained through the proposed method was lower than the static mesh stiffness, and the DTE obtained through the proposed method was higher than that using static mesh stiffness. Moreover, the system became more susceptible as it became more flexible. Case 1 was the least affected, and Case 4 was the most affected.
The acceleration signal at the bearing of the driving shaft in Figure 8 was extracted and presented in Figure 13. In Figure 13a, the highest vibration acceleration response is obtained in Case 1, whereas the lowest response is exhibited in Case 4. This discrepancy can be attributed to the fact that the shaft in Case 1 was the shortest and the bearing was positioned closest to the gear. Furthermore, the significant impact of variations in mesh stiffness on the dynamic response calculations can be observed in Figure 13b-e. In the spectrum, the components corresponding to the first four mesh frequencies display relatively minor changes, while the fifth and sixth mesh frequencies experience more pronounced alterations. The reason for this discrepancy is that changes in mesh stiffness affect both the dynamic characteristics and internal excitation of the system. that influence mesh stiffness. The proposed sparsification method was then applied to the PCE model, and 21 groups of mesh stiffness with random input parameters were calculated and employed as the testing set. As shown in Figure 11, the coefficients of the PCE model decreased gradually with the increase in iterations. During the sparsification process, R 2 remained almost unchanged, and the RMSE decreased slightly, which suggests that the sparse model maintains the accuracy of the original model.  Utilizing the proposed method (depicted in Figure 7), the dynamic mesh stiffness, DTE, and dynamic mesh force were extracted from the results of different cases, as illustrated in Figure 12. The average computation time using the SPCE model was 4.28 s, while the average computation time using the PCE model was 19.27 s. This indicates that the SPCE model significantly improves computational efficiency. Additionally, it is evident that the dynamic mesh stiffness obtained through the proposed method was lower than the static mesh stiffness, and the DTE obtained through the proposed method was higher than that using static mesh stiffness. Moreover, the system became more susceptible as it Figure 11. The process of sparsification. trated in Figure 12. The average computation time using the SPCE model was 4.28 s, while the average computation time using the PCE model was 19.27 s. This indicates that the SPCE model significantly improves computational efficiency. Additionally, it is evident that the dynamic mesh stiffness obtained through the proposed method was lower than the static mesh stiffness, and the DTE obtained through the proposed method was higher than that using static mesh stiffness. Moreover, the system became more susceptible as it became more flexible. Case 1 was the least affected, and Case 4 was the most affected. The acceleration signal at the bearing of the driving shaft in Figure 8 was extracted and presented in Figure 13. In Figure 13a, the highest vibration acceleration response is In order to explore the influence of the loads on the dynamic mesh stiffness, different loads of 200 Nm, 300 Nm, 400 Nm, and 500 Nm were substituted into Case 1-4, and the dynamic mesh stiffness values were obtained and are shown in Figure 14. It can be seen that Case 4 was the most sensitive to the change in loads. In Cases 1-3, the mesh stiffness was relatively smaller in the area with more mesh teeth when subjected to a greater load. Conversely, in the area with fewer mesh teeth, the mesh stiffness was relatively larger under a higher load. In Case 4, the mesh stiffness decreased throughout the whole mesh cycle with the increase in load. Figure 10 illustrates that the increase in the mesh force resulted in an increase in mesh stiffness, whereas the increases in the misalignment angle and center distance led to a decrease in stiffness. The phenomenon depicted in Figure 14 is clearly related to the combined effect of several factors. In Case 4, greater flexibility resulted in higher system deformation when subjected to the increase in load. This deformation subsequently increased the misalignment angle and center distance, leading to a lower mesh stiffness with greater load. Conversely, in Cases 1-3 with relatively small flexibility, the increase in the load in the area with less mesh teeth led to an increase in the mesh force. The influence of the mesh force was greater than that of the misalignment angle and center distance, resulting in an increase in mesh stiffness.
can be attributed to the fact that the shaft in Case 1 was the shortest and the bearing was positioned closest to the gear. Furthermore, the significant impact of variations in mesh stiffness on the dynamic response calculations can be observed in Figure 13b-e. In the spectrum, the components corresponding to the first four mesh frequencies display relatively minor changes, while the fifth and sixth mesh frequencies experience more pronounced alterations. The reason for this discrepancy is that changes in mesh stiffness affect both the dynamic characteristics and internal excitation of the system. In order to explore the influence of the loads on the dynamic mesh stiffness, different loads of 200 Nm, 300 Nm, 400 Nm, and 500 Nm were substituted into Case 1-4, and the dynamic mesh stiffness values were obtained and are shown in Figure 14. It can be seen that Case 4 was the most sensitive to the change in loads. In Cases 1-3, the mesh stiffness was relatively smaller in the area with more mesh teeth when subjected to a greater load. Conversely, in the area with fewer mesh teeth, the mesh stiffness was relatively larger under a higher load. In Case 4, the mesh stiffness decreased throughout the whole mesh cycle with the increase in load. Figure 10 illustrates that the increase in the mesh force resulted in an increase in mesh stiffness, whereas the increases in the misalignment angle is clearly related to the combined effect of several factors. In Case 4, greater flexibility resulted in higher system deformation when subjected to the increase in load. This deformation subsequently increased the misalignment angle and center distance, leading to a lower mesh stiffness with greater load. Conversely, in Cases 1-3 with relatively small flexibility, the increase in the load in the area with less mesh teeth led to an increase in the mesh force. The influence of the mesh force was greater than that of the misalignment angle and center distance, resulting in an increase in mesh stiffness.

Conclusions
In this study, a dynamic mesh stiffness and dynamic response analysis method was presented, in which the mesh stiffness at each sub-step is updated according to the gear displacement vector. The mesh stiffness is determined through 3D LTCA, and the computational efficiency is improved through an SPCE surrogate model. The main conclusions are summarized as follows: (1) The SPCE model proposed in this paper has high accuracy and efficiency. It can accurately reflect the influence of mesh force, angular misalignment, and center distance on mesh stiffness. Moreover, the proposed sparsification method can effectively reduce the number of coefficients without affecting the accuracy of the model. (2) Compared with the dynamic model with static mesh stiffness, the proposed method resulted in lower mesh stiffness and a larger transmission error, and the system became more susceptible as it became more flexible. Additionally, the proposed method also led to a significant difference in the system response. In the cases analyzed in this paper, the 5th and 6th mesh frequency components were significantly reduced when considering the changes in the engagement force, center distance, and angular misalignment during the calculation, while the 1st, 2nd, 3rd, and 4th mesh frequency components also exhibited significant changes, albeit to a much lesser degree.

Conclusions
In this study, a dynamic mesh stiffness and dynamic response analysis method was presented, in which the mesh stiffness at each sub-step is updated according to the gear displacement vector. The mesh stiffness is determined through 3D LTCA, and the computational efficiency is improved through an SPCE surrogate model. The main conclusions are summarized as follows: (1) The SPCE model proposed in this paper has high accuracy and efficiency. It can accurately reflect the influence of mesh force, angular misalignment, and center distance on mesh stiffness. Moreover, the proposed sparsification method can effectively reduce the number of coefficients without affecting the accuracy of the model. (2) Compared with the dynamic model with static mesh stiffness, the proposed method resulted in lower mesh stiffness and a larger transmission error, and the system became more susceptible as it became more flexible. Additionally, the proposed method also led to a significant difference in the system response. In the cases analyzed in this paper, the 5th and 6th mesh frequency components were significantly reduced when considering the changes in the engagement force, center distance, and angular misalignment during the calculation, while the 1st, 2nd, 3rd, and 4th mesh frequency components also exhibited significant changes, albeit to a much lesser degree. (3) An increase in mesh force leads to an increase in mesh stiffness and greater system deformation, which leads to a decrease in mesh stiffness. The interaction between these factors makes the dynamic mesh stiffness not always increase or decrease with the increase in torque. For a gear system with greater rigidity, the mesh stiffness may increase slightly with the increase in torque at some moments, while for a gear system with greater flexibility, the mesh stiffness will decrease with the increase in torque.