# Prediction Mathematic Model and Influencing Factors of Contact Stress of Cylindrical Gear with Arc Tooth

Authors: Qi Zhang; Guang Wen; Shuang Liang; Qin Tong; Li Hou; Guangchun Yang
Journal: Complexity
 (2020)

---

## Abstract

Given the absence of a theoretical formula to analyze the influence of parameters on the contact stress of cylindrical gear with arc tooth, an explicit mathematical model of cylindrical gear with arc tooth between the design parameters and the contact stress is established based on Kriging surrogate model. The parameters of the variation function of Kriging model are optimized by using the whale optimization algorithm (WOA), and the explicit mathematical model accuracy between the design parameters and the contact stress of the gear is in turn optimized by the improved Kriging surrogate model. The influence of design parameters on the contact stress of cylindrical gear with arc tooth is analyzed based on the established mathematical model. The proposed algorithm was realized via the programming platform MATLAB; the simulation results indicate that the precision evaluation indexes (the correlation coefficient (R2), root mean square error (RMSE), and the relative maximum absolute error (RMAE)) of the proposed Kriging model are improved, in addition to the error range which is narrowed from (−2, 4) to (0, 3). As the tooth width, modulus, pressure angle, and tooth line radius increased, the contact stress of the cylindrical gear with arc tooth gear declined, which was negatively correlated with the design parameter. The amplitude of contact stress of the cylindrical gear with arc tooth was the largest due to the change of tooth radius, followed by the change of modulus, while the influence of tooth width was less. Finally, the influence of modulus-tooth line radius interaction and pressure angle-tooth line radius interaction on contact stress of cylindrical gear with arc tooth was significant.


---

## Body

## 1. Introduction

Gears are widely used in various industries [1], and the involute cylindrical gear is widely used at present. However, in the involute cylindrical gear, the bearing capacity of spur gear is limited, the axial force of helical gear subsists, and the processing technology of herringbone gear is complex. In order to solve the problem of involute gear, Japanese scholar Yoshiro Hasegawa and others put forward a new type of circular arc gear transmission device. This new type of gear transmission has the advantages of good meshing performance, high coincidence, no axial force, stable transmission, and so on [2]. For this kind of gear, scholars have carried out in-depth research. Tseng et al. have studied its mathematical model, undercutting conditions, contact, and other aspecst [3–6]. Lei et al. have studied the accurate three-dimensional modeling method of circular tooth profile gear [7]. Dengqiu et al. have studied its contact distribution [8]. The contact performance of gears was studied by Chen Shuai and Wang [9, 10, 11, 12]. Surrogate model is an approximate modeling method that can generate mathematical models with a small number of samples to reduce the number of physical experiments and improve efficiency, and the commonly used surrogate model mainly includes response surface, Kriging, Support Vector Machine, Artificial Neural Network, etc. [13, 14]. In order to study the influence of design parameters on the contact stress of circular arc tooth line gear, this paper proposes using Kriging surrogate model to establish the mathematical model of contact stress between design parameters (tooth width, modulus, pressure angle, and tooth radius) of cylindrical gear with arc tooth. The Kriging model with the best effect is selected by comparing the prediction effect of these common surrogate model methods on contact stress of cylindrical gear with arc tooth. With a view to improving the modeling accuracy of Kriging model, intelligent algorithm is used to optimize the variation function of Kriging; the commonly used intelligent algorithm mainly includes genetic algorithm, immune algorithm, simulated annealing algorithm, particle swarm algorithm, fish swarm algorithm, ant colony algorithm, cuckoo algorithm, and so forth [15]. In this paper, whale optimization algorithm (WOA) was used for the parameters of the variation function of the traditional Kriging model. The Kriging surrogate model with optimized parameters of the variation function is used to establish the explicit mathematical model between the design parameters and the contact stress of the gear. The influence of design parameters (tooth width, modulus, pressure angle, and tooth radius) on the contact stress of circular arc tooth line cylindrical gear is discussed based on the established mathematical model.


## 2. Mathematical Model of Circular Arc Gear

According to its forming principle [16, 17], the coordinate system of forming principle of circular arc tooth line cylindrical gear is shown in Figure 1. In Figure 1, SO−XYZ is the static coordinate system, S1O1−X1Y1Z1 is the solidification coordinate of the gear blank, and STOT−XTYTZT is the tool coordinate, which moves in relation to SO−XYZ coordinate at the speed of VT = R × ω.Figure 1
Coordinate system of forming principle of cylindrical gear with arc tooth.In the coordinate systemSTOT−XTYTZT, the vector expression of the parametric equation of the surface formed in the cutting process is(1)rT⟶=−±qsinα+RT±πm4cosθiT⟶+±qsinα+RT±πm4sinθjT⟶+qcosαkT⟶,where a is the pressure angle of the cutting edge of the cutting tool; θ is the angle from the current enveloping point of the gear blank to the middle section; and q is the length of the enveloping reference point from the point on the cutting surface of the cutting tool along the direction of the generatrix of the cutting cone in the coordinate system. R is dividing circle radius of gear, RT is cutter head radius of machining tool, Rn is inner edge radius in the tangent direction of dividing circle, Rw is outer edge radius in the tangent direction of dividing circle, Rw = RT + πm/4, m is gear module, ω is cutter head rotation angle speed of machining tool, and VT is cutter head moving speed of machining tool.TransformSTOT−XTYTZT to SO−XYZ, and its coordinate conversion relationship is as follows:(2)T0T=100Rφ+RT001R01000001.TransformSO−XYZ to S1O1−X1Y1Z1, and its coordinate conversion relationship is(3)T10=cosφ−sinφ00sinφcosφ0000100001.IfSTOT−XTYTZT is converted to SO−XYZ, then the coordinate transformation matrix is T1T=T10T0T.


### 2.1. Unit Normal Vector of Tool Surface

The tool surface equation is(4)rdq⟶=±sinαcosθiT⟶±sinαsinθjT⟶+cosαkT,⟶rdθ⟶=−±qsinα+RT+πm4sinθiT⟶+±qsinα+RT+πm4cosθjT⟶.Then, the unit normal vector of the tool surface is(5)n⟶=rdq⟶×rdθ⟶rdq⟶×rdθ⟶=−cosαcosθiT⟶−cosαsinθjT⟶±sinαkT⟶.


### 2.2. Relative Speed of Tool and Gear at Meshing Point

The direction vector can be expressed as(6)λ⟶=O1OT⟶=Rφ+RTiT⟶+RkT⟶.Therefore, the relative speed of the tool and the gear at the meshing point is(7)v⟶12=ω⟶12×r⟶1+dλdt−ω⟶2×λ⟶,=−ω1qcosαiT⟶−ω1±qsinα+RT+πm4cosθkT⟶+ω1RiT⟶−ω1RiT⟶−ω1Rφ+RTkT⟶,=−qcosαω1iT⟶−±qsinα+RT+πm4cosθ+Rφ+RTω1kT.⟶


### 2.3. Meshing Function

Based on the meshing principle, the meshing function is expressed as follows:(8)Γ=n⟶⋅v⟶12=0.Becauseω1 ≠ 0, it can be concluded from the above formula that(9)q=∓sinαcosθRT+πm/4+Rφ+RTcosθ.


### 2.4. Conjugate Surfaces

Combining equations (1) and (9), the contact line equation between the tool and the gear to be machined is as follows:(10)xT=−±qsinα+RT±πm4cosθ,yT=±qsinα+RT±πm4sinθ,zT=qcosα,q=∓sinαcosθRT+πm/4+Rφ+RTcosθ.By converting the coordinates inSTOT−XTYTZT to S1O1−X1Y1Z1, the tooth surface equation of the cut gear can be obtained as follows:(11)x1=−±qsinα+RT±πm4cosθ+Rφ+RTcosφ−qcosα+Rsinφ,y1=−±qsinα+RT±πm4cosθ+Rφ+RTsinφ+qcosα+Rcosφ,z1=±qsinα+RT±πm4sinθ,q=∓sinαcosθRT+πm/4+Rφ+RTcosθ.


### 2.5. Instantaneous Contact Line

Meshing functionГ=Гq,θ,φ is a function of variables q,θ,φ. As φ1=ω1t, φ1 is regarded as a function of entering Г in a certain instant, thus the function expression of the instantaneous contact line between the rack cutter and the gear blank can be obtained. If the value φ1 at any time is substituted into equation (7), the contact line equation of the whole surface can be obtained as follows:(12)xT=−1+sin2αRT±πm4+sin2αRφ+RTcosθcosθ,yT=1+sin2αRT±πm4+sin2αRφ+RTcosθsinθ,zT=∓sinαcosαRT±πm4+Rφ+RTcosθ.


### 2.6. Tooth Profile Equation

In the axial middle section of the gear, from its extended coordinate system, it can be seen thatb = 0; then, θ is also equal to zero. Substituting these two parameters into equation (11), the tooth profile equation of the middle section is obtained:(13)x1=−±qsinα+RT±πm4+Rφ+RTcosφ−qcosα+Rsinφ,y1=−±qsinα+RT±πm4+Rφ+RTsinφ+qcosα+Rcosφ,q=∓sinαRT+πm4+Rφ+RT.According to equation (13), it can be seen that the tooth profile of the axial symmetry surface of the gear is involute.Similarly, on the axial asymmetric plane, letz1 = b; according to the expressions of Z1 and q in equation (11), we can get(14)φ=−b/tanθ+RT±πm/4cosθ/sin2α−RT±πm/4cosθ−RTR.Thus, the expression of the tooth profile of the non-intermediate section can be obtained as follows:(15)xT=−±qsinα+RT±πm4cosθ+Rφ+RTcosφ−qcosα+Rsinφ,yT=−±qsinα+RT±πm4cosθ+Rφ+RTsinφ+qcosα+Rcosφ,q=∓sinαcosθRT+πm/4+Rφ+RTcosθ,φ=−b/tanθ+RT±πm/4cosθ/sin2α−RT±πm/4cosθ−RTR.


## 2.1. Unit Normal Vector of Tool Surface

The tool surface equation is(4)rdq⟶=±sinαcosθiT⟶±sinαsinθjT⟶+cosαkT,⟶rdθ⟶=−±qsinα+RT+πm4sinθiT⟶+±qsinα+RT+πm4cosθjT⟶.Then, the unit normal vector of the tool surface is(5)n⟶=rdq⟶×rdθ⟶rdq⟶×rdθ⟶=−cosαcosθiT⟶−cosαsinθjT⟶±sinαkT⟶.


## 2.2. Relative Speed of Tool and Gear at Meshing Point

The direction vector can be expressed as(6)λ⟶=O1OT⟶=Rφ+RTiT⟶+RkT⟶.Therefore, the relative speed of the tool and the gear at the meshing point is(7)v⟶12=ω⟶12×r⟶1+dλdt−ω⟶2×λ⟶,=−ω1qcosαiT⟶−ω1±qsinα+RT+πm4cosθkT⟶+ω1RiT⟶−ω1RiT⟶−ω1Rφ+RTkT⟶,=−qcosαω1iT⟶−±qsinα+RT+πm4cosθ+Rφ+RTω1kT.⟶


## 2.3. Meshing Function

Based on the meshing principle, the meshing function is expressed as follows:(8)Γ=n⟶⋅v⟶12=0.Becauseω1 ≠ 0, it can be concluded from the above formula that(9)q=∓sinαcosθRT+πm/4+Rφ+RTcosθ.


## 2.4. Conjugate Surfaces

Combining equations (1) and (9), the contact line equation between the tool and the gear to be machined is as follows:(10)xT=−±qsinα+RT±πm4cosθ,yT=±qsinα+RT±πm4sinθ,zT=qcosα,q=∓sinαcosθRT+πm/4+Rφ+RTcosθ.By converting the coordinates inSTOT−XTYTZT to S1O1−X1Y1Z1, the tooth surface equation of the cut gear can be obtained as follows:(11)x1=−±qsinα+RT±πm4cosθ+Rφ+RTcosφ−qcosα+Rsinφ,y1=−±qsinα+RT±πm4cosθ+Rφ+RTsinφ+qcosα+Rcosφ,z1=±qsinα+RT±πm4sinθ,q=∓sinαcosθRT+πm/4+Rφ+RTcosθ.


## 2.5. Instantaneous Contact Line

Meshing functionГ=Гq,θ,φ is a function of variables q,θ,φ. As φ1=ω1t, φ1 is regarded as a function of entering Г in a certain instant, thus the function expression of the instantaneous contact line between the rack cutter and the gear blank can be obtained. If the value φ1 at any time is substituted into equation (7), the contact line equation of the whole surface can be obtained as follows:(12)xT=−1+sin2αRT±πm4+sin2αRφ+RTcosθcosθ,yT=1+sin2αRT±πm4+sin2αRφ+RTcosθsinθ,zT=∓sinαcosαRT±πm4+Rφ+RTcosθ.


## 2.6. Tooth Profile Equation

In the axial middle section of the gear, from its extended coordinate system, it can be seen thatb = 0; then, θ is also equal to zero. Substituting these two parameters into equation (11), the tooth profile equation of the middle section is obtained:(13)x1=−±qsinα+RT±πm4+Rφ+RTcosφ−qcosα+Rsinφ,y1=−±qsinα+RT±πm4+Rφ+RTsinφ+qcosα+Rcosφ,q=∓sinαRT+πm4+Rφ+RT.According to equation (13), it can be seen that the tooth profile of the axial symmetry surface of the gear is involute.Similarly, on the axial asymmetric plane, letz1 = b; according to the expressions of Z1 and q in equation (11), we can get(14)φ=−b/tanθ+RT±πm/4cosθ/sin2α−RT±πm/4cosθ−RTR.Thus, the expression of the tooth profile of the non-intermediate section can be obtained as follows:(15)xT=−±qsinα+RT±πm4cosθ+Rφ+RTcosφ−qcosα+Rsinφ,yT=−±qsinα+RT±πm4cosθ+Rφ+RTsinφ+qcosα+Rcosφ,q=∓sinαcosθRT+πm/4+Rφ+RTcosθ,φ=−b/tanθ+RT±πm/4cosθ/sin2α−RT±πm/4cosθ−RTR.


## 3. Contact Stress Analysis of Cylindrical Gear with Arc Tooth Based on Finite Element Method

### 3.1. Model Material Attribute Definition

Establish material information in ABAQUS, such as modulus of elasticity (E = 2.08 MPa), Poisson’s ratio (0.298), and other parameters.


### 3.2. Establishing Analytical Steps and Interactions

Establish analysis steps, which mainly include defining analysis type (static), defining analysis increment step, determining iteration method, creating field variable, and historical variable and determining output parameters. Turn on nonlinearity, and define interaction as contact; then set contact type of finite element analysis of circular arc tooth line gear pair as “friction free.” The MPC (multipoint constraint) is established by the rotating centers of the active and the driven.


### 3.3. Applying Restraint and Moment

The torque of magnitude 6.08 × 104 N·mm is applied to the driving wheel. At the same time, MPC is added to the driving wheel: the direction of its rotation axis is set as free, and other rotation and translation are set as fixed. In the driven wheel, all directions of MPC constraints are fixed, and no motion is allowed in any direction.


### 3.4. Gridding

The method of sweeping is used to divide the mesh of gear pair, with C3D8I as the mesh type. The overall unit size of the gear is set to 2 mm, while the local subdivision size of the contact area is 0.02 mm. If the requirements are not met, continue to adjust the parameters. In the subsequent calculation, the mesh can be further refined to improve the quality of division. Through multiple trial calculation and comparative analysis results, the grid can be used as the final analysis grid if the calculation result is nearly stable. In the analysis process, due to the local subdivision of the contact area, the number of grids for each pair of gears is about 1.2 × 106. The result of gear grid division is shown in Figure 2.Figure 2
Meshing results of gears. (a) Global mesh results. (b) Local mesh results.
(a)(b)


### 3.5. Solution and Visualization

Considering that the mesh number of each pair of gears is about 1.2 × 106, a parallel calculation is used in the solution. The calculated stress nephogram of driving wheel and that of driven wheel are shown in Figure 3. Figure 3(a) shows the dynamic contact stress diagram of driving wheel, and Figure 3(b) shows the contact stress diagram of driving wheel.Figure 3
Stress neutrogram analysis of driving wheel and driven wheel. (a) Driving wheel. (b) Driven wheel.
(a)(b)According to the simulation analysis, contact area is near dividing circle. As shown in Figure3(a), the driving wheel of the contact area of the above standard pitch circle has a maximum stress value of 503.3 MPa. As shown in Figure 3(b), driven wheel contact area of the below standard pitch circle has a maximum stress value of 501 MPa. Therefore, the maximum contact stress is 503.3 MPa, and the difference in value between active and driven wheel is 2.3 MPa, only equal to 0.4570% of the maximum contact stress, which almost can be ignored. Similarly, it can be seen from Figures 3(a) and 3(b) that the contact area of the gear studied in this paper is an ellipse under the action of load, which confirms the fact that the contact area of point contact gear under the action of load is ellipse.


## 3.1. Model Material Attribute Definition

Establish material information in ABAQUS, such as modulus of elasticity (E = 2.08 MPa), Poisson’s ratio (0.298), and other parameters.


## 3.2. Establishing Analytical Steps and Interactions

Establish analysis steps, which mainly include defining analysis type (static), defining analysis increment step, determining iteration method, creating field variable, and historical variable and determining output parameters. Turn on nonlinearity, and define interaction as contact; then set contact type of finite element analysis of circular arc tooth line gear pair as “friction free.” The MPC (multipoint constraint) is established by the rotating centers of the active and the driven.


## 3.3. Applying Restraint and Moment

The torque of magnitude 6.08 × 104 N·mm is applied to the driving wheel. At the same time, MPC is added to the driving wheel: the direction of its rotation axis is set as free, and other rotation and translation are set as fixed. In the driven wheel, all directions of MPC constraints are fixed, and no motion is allowed in any direction.


## 3.4. Gridding

The method of sweeping is used to divide the mesh of gear pair, with C3D8I as the mesh type. The overall unit size of the gear is set to 2 mm, while the local subdivision size of the contact area is 0.02 mm. If the requirements are not met, continue to adjust the parameters. In the subsequent calculation, the mesh can be further refined to improve the quality of division. Through multiple trial calculation and comparative analysis results, the grid can be used as the final analysis grid if the calculation result is nearly stable. In the analysis process, due to the local subdivision of the contact area, the number of grids for each pair of gears is about 1.2 × 106. The result of gear grid division is shown in Figure 2.Figure 2
Meshing results of gears. (a) Global mesh results. (b) Local mesh results.
(a)(b)


## 3.5. Solution and Visualization

Considering that the mesh number of each pair of gears is about 1.2 × 106, a parallel calculation is used in the solution. The calculated stress nephogram of driving wheel and that of driven wheel are shown in Figure 3. Figure 3(a) shows the dynamic contact stress diagram of driving wheel, and Figure 3(b) shows the contact stress diagram of driving wheel.Figure 3
Stress neutrogram analysis of driving wheel and driven wheel. (a) Driving wheel. (b) Driven wheel.
(a)(b)According to the simulation analysis, contact area is near dividing circle. As shown in Figure3(a), the driving wheel of the contact area of the above standard pitch circle has a maximum stress value of 503.3 MPa. As shown in Figure 3(b), driven wheel contact area of the below standard pitch circle has a maximum stress value of 501 MPa. Therefore, the maximum contact stress is 503.3 MPa, and the difference in value between active and driven wheel is 2.3 MPa, only equal to 0.4570% of the maximum contact stress, which almost can be ignored. Similarly, it can be seen from Figures 3(a) and 3(b) that the contact area of the gear studied in this paper is an ellipse under the action of load, which confirms the fact that the contact area of point contact gear under the action of load is ellipse.


## 4. The Prediction Model of Contact Stress of Cylindrical Gear with Arc Tooth Based on Surrogate Model

### 4.1. Test Design and Response Results

The common test methods mainly include uniform test design, orthogonal test design, Latin square test design, and other methods. This paper studies the relationship between tooth width, modulus, pressure angle, tooth radius, and contact force. Four factors and three levels are selected for the simulation scheme by using the principle of regression orthogonal test method. The level table of gear design factors is shown in Table1.Table 1
Gear design factor level.
Level orderTooth width (mm)ModulusPressure angle (°)Tooth radius (mm)Level 138.836852.9673718.83685191.8425Level 2403.220250Level 341.163153.4326321.16315308.1575In this paper, the orthogonal test is used to sample, and then the finite element method is used to get the gear joint stress under different sample data. According to the orthogonal table, 9 samples were processed, and the contact stress simulation results ofL9 (34) orthogonal test arc gear are shown in Table 2.Table 2
Contact stress simulation results ofL9 (34) orthogonal test.
OrderTooth width (mm)ModulusPressure angle (°)Tooth radius (mm)Contact force (MPa)138.836852.9673718.83685191.8425542.7238.836853.220250457.2338.836853.4326321.16315308.1575396.84402.9673720308.1575452.95403.221.16315191.8425494.16403.4326318.83685250449.1741.163152.9673721.16315250469.6841.163153.218.83685308.1575440.2941.163153.4326320191.8425476.7


### 4.2. Improved Kriging Model

#### 4.2.1. Kriging Surrogate Model

Kriging model can obtain ideal fitting results when solving the problem of high nonlinearity. Its interpolation results are defined as the linear weighting of the response value of the known sample function, which is(16)yx=∑j=1kβjfjx+Zx.In the formula,fjx is function, βj is corresponding coefficient, and Zx is static random process, whose mean value is 0, and its variance is σ2. The covariance between the random variables corresponding to two different points in the design space is(17)CovZxi,Zxj=σ2Rxi,xj,(18)Rxi,xj=∏l=1NsRlθl,xli,xlj,where Rxi,xj is correlation function, which represents the correlation between random variables at different positions. The commonly used correlation function is Gaussian-type function.In equation (17), the parameter θ of the variation function of Kriging model is determined by the maximum likelihood estimation method to solve the optimization problem:(19)max−Ns2lnσ˜2−12lnR−Ns2ln2π,θ>0,σ˜2=1Nsys−Fβ˜TR−1ys−Fβ˜.In order to minimize the root mean square error (RMSE) between Kriging prediction value and actual function value, the approximate expression of Kriging model can be obtained as(20)y˜x=fxTβ˜+Rx,x1,Rx,x2,…,Rx,xNsTR−1ys−Fβ˜,where F is matrix composed of basis function vector fxT; β˜ is least square estimate of KRG model coefficient; R is correlation matrix; and Ns is sample number.


#### 4.2.2. Whale Algorithm

The principle of whale algorithm is derived from the“bubble net,” foraging behavior of humpback whales. The mathematical expression based on this special strategy is as follows:(21)D=C⋅X∗t−Xt,Xt+1=X∗t−A⋅D,where t is the number of current iterations; Xt is the coordinate vector of the current humpback whale; Xt+1 is the target coordinate vector after the next iteration; X∗t is the best position vector obtained so far, which is updated with time; D is the distance between the current humpback whale and the best position; and A and C are the number of systems, which are, respectively, expressed as(22)A=2a⋅r−a,C=2r,where a is parameter in the value range (0, 2) and linearly decreasing with the iteration time and r is random vector in the interval (0, 1).WhenA>1, whales have the behavior of wandering and foraging. Using the random individual coordinates of the population to locate and navigate food, the mathematical expression is as follows:(23)Xt+1=Xrandt−A⋅D.WhenA<1, there were two behaviors of encircling and attacking prey, and the mathematical model is described as follows:(24)Xt+1=X∗t−D⋅ebl⋅cos2πl.In the formula, b is the constant of the spiral shape andl is the random number on the interval (−1, 1).Since the contraction and spiral position updating of whales are carried out synchronously, mathematically selecting the same way of probability to update their position, then we can get the following expression:(25)Xt+1=X∗t−A⋅D,p<0.5,X∗t−D⋅ebl⋅cos2πl,p≥0.5.According to the characteristics of Kriging and WOA algorithm, the Kriging surrogate model based on WOA algorithm is used for parameter impact analysis, and the flow chart is shown in Figure4.Figure 4
Improvement process of Kriging surrogate model based on WOA algorithm.


### 4.3. Accuracy Evaluation Indexes of Kriging

Generally, the evaluation indexes of correlation coefficient (R2), root mean square error (RMSE), and relative maximum absolute error (RMAE) are used to evaluate the accuracy of Kriging.


#### 4.3.1. The Correlation Coefficient (R2)

(26)R2=1−∑i=1nfi−f^i2∑i=1nfi−f¯2.


#### 4.3.2. The Root Mean Square Error (RMSE)

(27)RMSE=∑i=1nfi−f^i2n.


#### 4.3.3. The Relative Maximum Absolute Error (RMAE)

(28)RMAE=maxf1−f^1,f2−f^2,…,fn−f^n1/n∑i=1nfi−f¯2.In equations (26) to (28), n is number of test sample points, fi is the truth value of the actual model at the ith test sample point, and f^i is the predicted value of the surrogate model at the ith test sample point. The closer R2 is to 1, the better the global approximation effect of the surrogate model will be. The smaller the RMSE is, the higher the precision of the surrogate model will be, and the smaller the RMAE is, indicating that the higher the precision of the surrogate model will be.


### 4.4. Numerical Simulation

Using the MATLAB digital simulation platform, Kriging and the improved Kriging algorithm based on WOA algorithm are, respectively, used to establish the contact stress prediction model of circular arc gear. The evaluation indexes of the correlation coefficient (R2), root mean square error (RMSE), and the relative maximum absolute error (RMAE) before and after optimization are shown in Table 3.Table 3
Evaluation indexes before and after optimization.
IndexesKrigingWOA-KrigingOptimization effect (%)R20.99220.99740.52RMSE2.85691.654042.11RMAE0.13220.075442.97It can be seen from Table3 that the correlation coefficient (R2) has increased from 0.9922 to 0.9974, which is 0.52% higher and closer to 1. The correlation coefficient is closer to 1 implying that the global approximation ability of the optimized Kriging model is better. RMSE is used to indicate the impairment between the estimated value and the actual value. The higher the value of RMSE, the greater the impairment. After optimization, RMSE went from 2.8569 to 1.6540, reducing by 42.11%, which shows that the optimized Kriging model can better estimate the samples. After optimization, RMAE decreased from 0.1322 to 0.0754, which is 42.97% lower and closer to 0. This shows that the local error is less after optimization. Therefore, Kriging algorithm is improved based on WOA algorithm, which enhances the fitting ability and accuracy of Kriging algorithm.The residual graph histogram of the test set for the contact stress prediction model is established by using Kriging and the improved Kriging algorithm based on WOA, as shown in Figures5(a) and 5(b).Figure 5
Residual diagram of contact stress prediction model. (a) Residual diagram of contact stress prediction model based on Kriging ((A) residual histogram, (B) residual with sequence, (C) residual, (D) sequence of observations). (b) Residual diagram of contact stress prediction model based on Kriging based on WOA ((A) residual histogram, (B) residual with sequence, (C) residual, (D) sequence of observations).
(a)(b)It can be seen from Figure5 that the error dispersion of the optimized Kriging based on the improved WOA to establish the contact stress prediction model of the circular arc gear is more centralized, and the optimized error range is reduced from (−2, 4) to (0, 3), which significantly improved accuracy.


### 4.5. Finite Element Verification

In order to further verify the validity of the established mathematical model, the finite element method is used to analyze the gear. The definition of material attribute of the model, the establishment of analysis steps and interactions, the application of constraints and moments, and mesh division are the same as those in Section3. The gear parameters are shown in Table 4, and the results of the finite element analysis are shown in Figure 6.Table 4
Gear parameters of finite element analysis.
OrderTooth width (mm)ModulusPressure angle (°)Tooth radius (mm)138.8368503.220308.157500238.8368502.96737018.836850308.1575003403.43263021.163150308.1575004403.43263020.000000191.8425005402.96737021.163150250Figure 6
Stress nephogram of driving wheel and driven wheel. (a) Driving wheel. (b) Driven wheel.
(a)(b)As shown in Figure6, the maximum contact stress of the gear is 472.2 MPa. The mathematical model established in Section 4.3 is used to predict the maximum contact stress of the gear in Table 4.As shown in Table5, the mathematical model established in Section 4.3 is used to predict the maximum contact stress of gears in Table 4. Compared with the results of finite element analysis, the maximum error is 3.2013 MPa, which shows that the mathematical model can effectively predict the contact stress of gears.Table 5
Finite element prediction results.
OrderTooth width (mm)ModulusPressure angle (°)Tooth radius (mm)Finite element stress (MPa)Predicted stress (MPa)Error (MPa)138.836853.220308.157500427.2429.79352.5935238.836852.96737018.836850308.157500462.3465.50133.20133403.43263021.163150308.157500396.4394.7222−1.67784403.43263020.000000191.842500479.2478.3456−0.85445402.96737021.163150250470471.80391.8039


## 4.1. Test Design and Response Results

The common test methods mainly include uniform test design, orthogonal test design, Latin square test design, and other methods. This paper studies the relationship between tooth width, modulus, pressure angle, tooth radius, and contact force. Four factors and three levels are selected for the simulation scheme by using the principle of regression orthogonal test method. The level table of gear design factors is shown in Table1.Table 1
Gear design factor level.
Level orderTooth width (mm)ModulusPressure angle (°)Tooth radius (mm)Level 138.836852.9673718.83685191.8425Level 2403.220250Level 341.163153.4326321.16315308.1575In this paper, the orthogonal test is used to sample, and then the finite element method is used to get the gear joint stress under different sample data. According to the orthogonal table, 9 samples were processed, and the contact stress simulation results ofL9 (34) orthogonal test arc gear are shown in Table 2.Table 2
Contact stress simulation results ofL9 (34) orthogonal test.
OrderTooth width (mm)ModulusPressure angle (°)Tooth radius (mm)Contact force (MPa)138.836852.9673718.83685191.8425542.7238.836853.220250457.2338.836853.4326321.16315308.1575396.84402.9673720308.1575452.95403.221.16315191.8425494.16403.4326318.83685250449.1741.163152.9673721.16315250469.6841.163153.218.83685308.1575440.2941.163153.4326320191.8425476.7


## 4.2. Improved Kriging Model

### 4.2.1. Kriging Surrogate Model

Kriging model can obtain ideal fitting results when solving the problem of high nonlinearity. Its interpolation results are defined as the linear weighting of the response value of the known sample function, which is(16)yx=∑j=1kβjfjx+Zx.In the formula,fjx is function, βj is corresponding coefficient, and Zx is static random process, whose mean value is 0, and its variance is σ2. The covariance between the random variables corresponding to two different points in the design space is(17)CovZxi,Zxj=σ2Rxi,xj,(18)Rxi,xj=∏l=1NsRlθl,xli,xlj,where Rxi,xj is correlation function, which represents the correlation between random variables at different positions. The commonly used correlation function is Gaussian-type function.In equation (17), the parameter θ of the variation function of Kriging model is determined by the maximum likelihood estimation method to solve the optimization problem:(19)max−Ns2lnσ˜2−12lnR−Ns2ln2π,θ>0,σ˜2=1Nsys−Fβ˜TR−1ys−Fβ˜.In order to minimize the root mean square error (RMSE) between Kriging prediction value and actual function value, the approximate expression of Kriging model can be obtained as(20)y˜x=fxTβ˜+Rx,x1,Rx,x2,…,Rx,xNsTR−1ys−Fβ˜,where F is matrix composed of basis function vector fxT; β˜ is least square estimate of KRG model coefficient; R is correlation matrix; and Ns is sample number.


### 4.2.2. Whale Algorithm

The principle of whale algorithm is derived from the“bubble net,” foraging behavior of humpback whales. The mathematical expression based on this special strategy is as follows:(21)D=C⋅X∗t−Xt,Xt+1=X∗t−A⋅D,where t is the number of current iterations; Xt is the coordinate vector of the current humpback whale; Xt+1 is the target coordinate vector after the next iteration; X∗t is the best position vector obtained so far, which is updated with time; D is the distance between the current humpback whale and the best position; and A and C are the number of systems, which are, respectively, expressed as(22)A=2a⋅r−a,C=2r,where a is parameter in the value range (0, 2) and linearly decreasing with the iteration time and r is random vector in the interval (0, 1).WhenA>1, whales have the behavior of wandering and foraging. Using the random individual coordinates of the population to locate and navigate food, the mathematical expression is as follows:(23)Xt+1=Xrandt−A⋅D.WhenA<1, there were two behaviors of encircling and attacking prey, and the mathematical model is described as follows:(24)Xt+1=X∗t−D⋅ebl⋅cos2πl.In the formula, b is the constant of the spiral shape andl is the random number on the interval (−1, 1).Since the contraction and spiral position updating of whales are carried out synchronously, mathematically selecting the same way of probability to update their position, then we can get the following expression:(25)Xt+1=X∗t−A⋅D,p<0.5,X∗t−D⋅ebl⋅cos2πl,p≥0.5.According to the characteristics of Kriging and WOA algorithm, the Kriging surrogate model based on WOA algorithm is used for parameter impact analysis, and the flow chart is shown in Figure4.Figure 4
Improvement process of Kriging surrogate model based on WOA algorithm.


## 4.2.1. Kriging Surrogate Model

Kriging model can obtain ideal fitting results when solving the problem of high nonlinearity. Its interpolation results are defined as the linear weighting of the response value of the known sample function, which is(16)yx=∑j=1kβjfjx+Zx.In the formula,fjx is function, βj is corresponding coefficient, and Zx is static random process, whose mean value is 0, and its variance is σ2. The covariance between the random variables corresponding to two different points in the design space is(17)CovZxi,Zxj=σ2Rxi,xj,(18)Rxi,xj=∏l=1NsRlθl,xli,xlj,where Rxi,xj is correlation function, which represents the correlation between random variables at different positions. The commonly used correlation function is Gaussian-type function.In equation (17), the parameter θ of the variation function of Kriging model is determined by the maximum likelihood estimation method to solve the optimization problem:(19)max−Ns2lnσ˜2−12lnR−Ns2ln2π,θ>0,σ˜2=1Nsys−Fβ˜TR−1ys−Fβ˜.In order to minimize the root mean square error (RMSE) between Kriging prediction value and actual function value, the approximate expression of Kriging model can be obtained as(20)y˜x=fxTβ˜+Rx,x1,Rx,x2,…,Rx,xNsTR−1ys−Fβ˜,where F is matrix composed of basis function vector fxT; β˜ is least square estimate of KRG model coefficient; R is correlation matrix; and Ns is sample number.


## 4.2.2. Whale Algorithm

The principle of whale algorithm is derived from the“bubble net,” foraging behavior of humpback whales. The mathematical expression based on this special strategy is as follows:(21)D=C⋅X∗t−Xt,Xt+1=X∗t−A⋅D,where t is the number of current iterations; Xt is the coordinate vector of the current humpback whale; Xt+1 is the target coordinate vector after the next iteration; X∗t is the best position vector obtained so far, which is updated with time; D is the distance between the current humpback whale and the best position; and A and C are the number of systems, which are, respectively, expressed as(22)A=2a⋅r−a,C=2r,where a is parameter in the value range (0, 2) and linearly decreasing with the iteration time and r is random vector in the interval (0, 1).WhenA>1, whales have the behavior of wandering and foraging. Using the random individual coordinates of the population to locate and navigate food, the mathematical expression is as follows:(23)Xt+1=Xrandt−A⋅D.WhenA<1, there were two behaviors of encircling and attacking prey, and the mathematical model is described as follows:(24)Xt+1=X∗t−D⋅ebl⋅cos2πl.In the formula, b is the constant of the spiral shape andl is the random number on the interval (−1, 1).Since the contraction and spiral position updating of whales are carried out synchronously, mathematically selecting the same way of probability to update their position, then we can get the following expression:(25)Xt+1=X∗t−A⋅D,p<0.5,X∗t−D⋅ebl⋅cos2πl,p≥0.5.According to the characteristics of Kriging and WOA algorithm, the Kriging surrogate model based on WOA algorithm is used for parameter impact analysis, and the flow chart is shown in Figure4.Figure 4
Improvement process of Kriging surrogate model based on WOA algorithm.


## 4.3. Accuracy Evaluation Indexes of Kriging

Generally, the evaluation indexes of correlation coefficient (R2), root mean square error (RMSE), and relative maximum absolute error (RMAE) are used to evaluate the accuracy of Kriging.


### 4.3.1. The Correlation Coefficient (R2)

(26)R2=1−∑i=1nfi−f^i2∑i=1nfi−f¯2.


### 4.3.2. The Root Mean Square Error (RMSE)

(27)RMSE=∑i=1nfi−f^i2n.


### 4.3.3. The Relative Maximum Absolute Error (RMAE)

(28)RMAE=maxf1−f^1,f2−f^2,…,fn−f^n1/n∑i=1nfi−f¯2.In equations (26) to (28), n is number of test sample points, fi is the truth value of the actual model at the ith test sample point, and f^i is the predicted value of the surrogate model at the ith test sample point. The closer R2 is to 1, the better the global approximation effect of the surrogate model will be. The smaller the RMSE is, the higher the precision of the surrogate model will be, and the smaller the RMAE is, indicating that the higher the precision of the surrogate model will be.


## 4.3.1. The Correlation Coefficient (R2)

(26)R2=1−∑i=1nfi−f^i2∑i=1nfi−f¯2.


## 4.3.2. The Root Mean Square Error (RMSE)

(27)RMSE=∑i=1nfi−f^i2n.


## 4.3.3. The Relative Maximum Absolute Error (RMAE)

(28)RMAE=maxf1−f^1,f2−f^2,…,fn−f^n1/n∑i=1nfi−f¯2.In equations (26) to (28), n is number of test sample points, fi is the truth value of the actual model at the ith test sample point, and f^i is the predicted value of the surrogate model at the ith test sample point. The closer R2 is to 1, the better the global approximation effect of the surrogate model will be. The smaller the RMSE is, the higher the precision of the surrogate model will be, and the smaller the RMAE is, indicating that the higher the precision of the surrogate model will be.


## 4.4. Numerical Simulation

Using the MATLAB digital simulation platform, Kriging and the improved Kriging algorithm based on WOA algorithm are, respectively, used to establish the contact stress prediction model of circular arc gear. The evaluation indexes of the correlation coefficient (R2), root mean square error (RMSE), and the relative maximum absolute error (RMAE) before and after optimization are shown in Table 3.Table 3
Evaluation indexes before and after optimization.
IndexesKrigingWOA-KrigingOptimization effect (%)R20.99220.99740.52RMSE2.85691.654042.11RMAE0.13220.075442.97It can be seen from Table3 that the correlation coefficient (R2) has increased from 0.9922 to 0.9974, which is 0.52% higher and closer to 1. The correlation coefficient is closer to 1 implying that the global approximation ability of the optimized Kriging model is better. RMSE is used to indicate the impairment between the estimated value and the actual value. The higher the value of RMSE, the greater the impairment. After optimization, RMSE went from 2.8569 to 1.6540, reducing by 42.11%, which shows that the optimized Kriging model can better estimate the samples. After optimization, RMAE decreased from 0.1322 to 0.0754, which is 42.97% lower and closer to 0. This shows that the local error is less after optimization. Therefore, Kriging algorithm is improved based on WOA algorithm, which enhances the fitting ability and accuracy of Kriging algorithm.The residual graph histogram of the test set for the contact stress prediction model is established by using Kriging and the improved Kriging algorithm based on WOA, as shown in Figures5(a) and 5(b).Figure 5
Residual diagram of contact stress prediction model. (a) Residual diagram of contact stress prediction model based on Kriging ((A) residual histogram, (B) residual with sequence, (C) residual, (D) sequence of observations). (b) Residual diagram of contact stress prediction model based on Kriging based on WOA ((A) residual histogram, (B) residual with sequence, (C) residual, (D) sequence of observations).
(a)(b)It can be seen from Figure5 that the error dispersion of the optimized Kriging based on the improved WOA to establish the contact stress prediction model of the circular arc gear is more centralized, and the optimized error range is reduced from (−2, 4) to (0, 3), which significantly improved accuracy.


## 4.5. Finite Element Verification

In order to further verify the validity of the established mathematical model, the finite element method is used to analyze the gear. The definition of material attribute of the model, the establishment of analysis steps and interactions, the application of constraints and moments, and mesh division are the same as those in Section3. The gear parameters are shown in Table 4, and the results of the finite element analysis are shown in Figure 6.Table 4
Gear parameters of finite element analysis.
OrderTooth width (mm)ModulusPressure angle (°)Tooth radius (mm)138.8368503.220308.157500238.8368502.96737018.836850308.1575003403.43263021.163150308.1575004403.43263020.000000191.8425005402.96737021.163150250Figure 6
Stress nephogram of driving wheel and driven wheel. (a) Driving wheel. (b) Driven wheel.
(a)(b)As shown in Figure6, the maximum contact stress of the gear is 472.2 MPa. The mathematical model established in Section 4.3 is used to predict the maximum contact stress of the gear in Table 4.As shown in Table5, the mathematical model established in Section 4.3 is used to predict the maximum contact stress of gears in Table 4. Compared with the results of finite element analysis, the maximum error is 3.2013 MPa, which shows that the mathematical model can effectively predict the contact stress of gears.Table 5
Finite element prediction results.
OrderTooth width (mm)ModulusPressure angle (°)Tooth radius (mm)Finite element stress (MPa)Predicted stress (MPa)Error (MPa)138.836853.220308.157500427.2429.79352.5935238.836852.96737018.836850308.157500462.3465.50133.20133403.43263021.163150308.157500396.4394.7222−1.67784403.43263020.000000191.842500479.2478.3456−0.85445402.96737021.163150250470471.80391.8039


## 5. Influence of Design Parameters on Contact Stress

In this section, the influence of design parameters on contact stress is discussed based on the explicitly mathematical model between the gear design parameters (tooth width, modulus, pressure angle, and tooth line radius) and the contact stress of cylindrical gear with arc tooth by using the proposed WOA-Kriging model. When the parameters (tooth width, modulus, pressure angle, and tooth line radius) are changed separately, the influence of tooth width, modulus, pressure angle, and tooth line radius on the contact stress of cylindrical gear with arc tooth is shown in Figure7.Figure 7
Influence of design parameters on contact stress. (a) Influence of tooth width on contact stress contact stress of cylindrical gear with arc tooth. (b) Influence of modulus on contact stress contact stress of cylindrical gear with arc tooth. (c) Influence of pressure angle on contact stress contact stress of cylindrical gear with arc tooth. (d) Influence of tooth line radius on contact stress contact stress of cylindrical gear with arc tooth.
(a)(b)(c)(d)As seen in Figure7, with the increase of tooth width, modulus, pressure angle, and tooth line radius, the contact stress of the cylindrical gear with arc tooth gear was in decline, which was negatively correlated with the design parameter, the amplitude of contact stress of the cylindrical gear with arc tooth is the largest due to the change of tooth radius, followed by the change of modulus, and the influence of tooth width is less.The influence of change in tooth width, modulus, pressure angle, and tooth line radius on contact stress of the cylindrical gear with arc tooth is shown in Figure8.Figure 8
Influence of design parameters interaction on contact stress of cylindrical gear with arc tooth. (a) Influence of tooth width-modulus interaction on contact stress of cylindrical gear with arc tooth. (b) Influence of tooth width-pressure angle interaction on contact stress of cylindrical gear with arc tooth. (c) Influence of tooth width-tooth line radius interaction on contact stress of cylindrical gear with arc tooth. (d) Influence of modulus-pressure angle interaction on contact stress of cylindrical gear with arc tooth. (e) Influence of modulus-tooth line radius interaction on contact stress of cylindrical gear with arc tooth. (f) Influence of pressure angle-tooth line radius interaction on contact stress of cylindrical gear with arc tooth.
(a)(b)(c)(d)(e)(f)As seen in Figure8, when tooth width interacts with modulus, pressure angle, and tooth line radius, respectively, the variation range of contact stress for the cylindrical gear with arc tooth gear is in decline. The response surface is not large and the complexity of surface is not high, indicating that the influence of tooth width interacting with modulus, pressure angle, and tooth line radius, respectively, on the contact stress is not significant. When modulus interacts with modulus, pressure angle, and tooth line radius, respectively, the complexity of the stress response surface is relatively high, which indicates that the influence on the stress is significant, compared with tooth width interacting with modulus, pressure angle, and tooth line radius, respectively. When pressure angle interacts with tooth line radius, the complexity of the stress response surface is relatively high, which also indicates that the influence on the stress is significant. Influence of modulus-tooth line radius interaction on contact stress of cylindrical gear with arc tooth and influence of pressure angle-tooth line radius interaction on contact stress of cylindrical gear with arc tooth are significant. In the contour map above, the smaller the distance between the contour lines is, the greater the influence of their interaction on the contact stress will be.


## 6. Conclusion

In order to study the influence of design parameters of cylindrical gear with arc tooth on its contact stress, a mathematical model of tooth width, modulus, pressure angle, radius of tooth line, and contact stress of gear is established by using Kriging surrogate model. Meanwhile, a Kriging model modeling method based on the whale optimization algorithm (WOA) is proposed. The parameters of variation function of traditional Kriging model are optimized by whale optimization calculation to improve the modeling accuracy of Kriging model. Based on the established mathematical model, the influence of design parameters (tooth width, module, pressure angle, and tooth radius) of cylindrical gear with arc tooth on the contact stress of circular arc tooth line cylindrical gear is discussed.(1)
Based on the improved Kriging model of WOA algorithm, the correlation coefficient (R2) is increased from 0.9922 to 0.9974, MSE is reduced from 2.8569 to 1.6540, and RMAE is also reduced from 0.1322 to 0.0754. The correlation coefficient (R2) of the improved algorithm is improved, and root mean square error (RMSE) and relative maximum absolute error (RMAE) are improved to some extent, which improves the global approximation ability of the traditional Kriging algorithm, reduces the local error, and improves the fitting accuracy.(2)
Both the traditional Kriging algorithm and the improved Kriging algorithm based on WOA algorithm can establish the contact stress prediction model of gear, and the accuracy is within the applicable range. The error range of Kriging is (−2, 4), while the error range of Kriging algorithm based on WOA algorithm is (0, 3). The accuracy of the improved Kriging has been optimized obviously, and the accuracy of the model is verified by the finite element method, where the maximum error is only 2.5935 MPa.(3)
Based on the established mathematical model, the influence of design parameters (tooth width, modulus, pressure angle, and tooth radius) on the contact stress of circular arc tooth line cylindrical gear is discussed. With the increase of tooth width, modulus, pressure angle, and tooth line radius, the contact stress of the cylindrical gear with arc tooth gear was in decline, which was negatively correlated with the design parameter, the amplitude of contact stress of the cylindrical gear with arc tooth is the largest due to the change of tooth radius, followed by the change of modulus, and the influence of tooth width is the least. And influence of modulus-tooth line radius interaction and pressure angle-tooth line radius interaction on contact stress of cylindrical gear with arc tooth is significant.



---

*Source: 8888407-2020-11-03.xml*