Time-dependenT reliabiliTy of spur gear sysTem based on gradually wear process

Friction and wear are unavoidable for transmitting power in gear systems. Severe wear can cause mechanical component damage. The dynamic characteristic of gear system is affected not only by friction and wear but also by variation of geometrical shape and dimension parameter caused by mild wear, which lead to more serious wear. Therefore, the coupled relation between wear and dynamic load should be taken into consideration to study gear wear. Gear parameters might not be accurate after machining in practice, which is one of reasons why product failure still occurs. That is to say, parameters have uncertainty and randomness. Besides, wear depth accumulates over time gradually, so gear wear is a random process, which has random statistical regularity. As a result, more practical surface wear model and dynamic reliability model with gradually parameters can be established from the definition of reliability, the reliability of mechanical component can be observed all over its life cycle, and the tolerance of designed parameters can be determined to reduce failure probability. The wear model should be established based on failure mechanism to study gear wear. In recent 30 years, wear and failure mechanism of gears have been studied based on lots of laboratory tests, Zurowski[28] et al. employed TT-3 tester to research wear resistance of C45(norm)/145Cr6 and C45(600)/145Cr6 matchings and indicated that wear resistance of the two matchings had a significant dependence on friction area temperature and material hardness. However, research on establishing models and calculation methods for gear wear have little been touched on. In general, these models break down into three categories: The first model was established on undetermined coefficient or regression analysis method according to test data. The conclusion was convincing but only obtained laws for intermediate variables. For example, Põdra[15] studied wear had a linear correlation with normal load and had incomplete correlation with sliding velocity. The second was developed to discuss the influence of parameters from energy loss aspect. For example, Onishchenko[14] studied the effect of machine operation, corresponding tribological theories, the eccentricity of pitch circle and the instant temperature in the contact ZHU L, ZHANG Y, ZHANG R, ZHANG P. Time-dependent reliability of spur gear system based on gradually wear process. Eksploatacja i Niezawodnosc – Maintenance and Reliability 2018; 20 (2): 207–218, http://dx.doi.org/10.17531/ein.2018.2.05.


Introduction
Friction and wear are unavoidable for transmitting power in gear systems.Severe wear can cause mechanical component damage.The dynamic characteristic of gear system is affected not only by friction and wear but also by variation of geometrical shape and dimension parameter caused by mild wear, which lead to more serious wear.Therefore, the coupled relation between wear and dynamic load should be taken into consideration to study gear wear.
Gear parameters might not be accurate after machining in practice, which is one of reasons why product failure still occurs.That is to say, parameters have uncertainty and randomness.Besides, wear depth accumulates over time gradually, so gear wear is a random process, which has random statistical regularity.As a result, more practical surface wear model and dynamic reliability model with gradually parameters can be established from the definition of reliability, the reliability of mechanical component can be observed all over its life cycle, and the tolerance of designed parameters can be determined to reduce failure probability.
The wear model should be established based on failure mechanism to study gear wear.In recent 30 years, wear and failure mechanism of gears have been studied based on lots of laboratory tests, Zurowski [28] et al. employed TT-3 tester to research wear resistance of C45(norm)/145Cr6 and C45(600)/145Cr6 matchings and indicated that wear resistance of the two matchings had a significant dependence on friction area temperature and material hardness.However, research on establishing models and calculation methods for gear wear have little been touched on.In general, these models break down into three categories: The first model was established on undetermined coefficient or regression analysis method according to test data.The conclusion was convincing but only obtained laws for intermediate variables.For example, Põdra [15] studied wear had a linear correlation with normal load and had incomplete correlation with sliding velocity.
The second was developed to discuss the influence of parameters from energy loss aspect.For example, Onishchenko [14] studied the effect of machine operation, corresponding tribological theories, the eccentricity of pitch circle and the instant temperature in the contact ZHU L, ZHANG Y, ZHANG R, ZHANG P. Time-dependent reliability of spur gear system based on gradually wear process.Eksploatacja i Niezawodnosc -Maintenance and Reliability 2018; 20 (2) on wear.The model stated that wear is simply proportional to specific power, which needed further study.
The third was the Archard wear model [3] which had been widely adopted.The formula had evolved from severe sliding of intermeshing tooth.Andersson [2] obtained analytical formula of sliding distance which varied with different mesh position.Flodin and Andersson [9] developed a numerical model for wear prediction of spur gear and the contact were modelled by Winkler's elastic foundation model.Park [16] combined Archard's wear model with a finite-element based hypoid gear contact model for simulation of surface wear of hypoid gear pairs.Most of the above adopted equal or linear distributed load [2,4,6], without coupled relations between tooth load and wear.Flodin [8,9] studied the effect of spur gears on contact condition, whereas, the sequence of contact, clearance size between two adjacent teeth and time-varying stiffness were not considered.
The parameters were deterministic in the wear models above, randomness of excited load and gradual wear were also not considered.As a matter of fact, wear-based reliability of mechanical product is varying gradually with time [22] .There are many methods to build time-dependent reliability model, such as Markov theory [19], dynamic fault-tree model [7,13], stochastic MCS model [5], Go-flow model [11], stochastic petri nets [20].As mechanical structure is so complex that the methods above are hardly applied to dynamic reliability analysis for mechanical systems and components.One is because the relationship between response and parameters are nonlinear and there is no analytic expression for the state limited function in general.The other is lack of statistical data and distribution law for parameters.
Reliability analysis is based on probability analysis and mathematical statistics, so it needs enough samples to determine probability characteristic of variables and response.Therefore, it would be time-consuming to obtain enough samples for complex structure.To establish reliability model for complex structure system, the accurate numerical model is replaced with one of surrogate models based on response surface, neutral network, radial basis function network, support vector machine and kriging method, and reliability are obtained by surrogate-model-based monte carlo simulation method.Gomes and Awruch [10] presented response surface and artificial neural network techniques to solve complex and more elaborated problems and carried out comparison using FORM, direct MCS and MCS with adaptive importance sampling.Tan [18] proposed radial basis function networks and support vector machines for reliability analysis.Zhu and Du [27] applied the kriging-based MCS method which was suitable for highly nonlinear limit-state functions to reduce the computational cost.Because the nonlinear response with complex structure can be rapidly predicted in any design point using a kriging interpolator and the iteration process has higher stabilization and faster convergence speed than other methods, the Kriging model is employed in the present paper.
To solve the problem of unknown distribution, on the basis of the surrogate model, Edgeworth series and Four moment techniques are suitable for reliability analysis.The method combines analysis method and moment method and Zhang [26] developed a statistical fourth moment method to examine reliability of the rotor-stator systems with rubbing.All the models and methods above have not been used in analysing wear-based reliability of gear system.
In the present paper, a numerical wear model for a pair of spur gears is established considering coupled relation of wear and dynamic load.A time-dependent reliability model with gradual parameters is also build under the condition that original parameters are stochastic and wear is gradually changing.The time-varying reliability curves for spur gears are obtained during the life cycle.

Load distribution theory
As shown in Figure 1, in one meshing cycle, single teethmeshing area and double teeth-meshing area appear alternately.In single teeth-meshing area, the transmission load is carried only by one pair of mesh gear.In double teeth-meshing area, the transmission load is carried by two pairs of mesh gear.If only considering elastic deformation without wear, the two pairs of mesh gear can be regarded as parallel springs.The distribution load is calculated according to time-varying stiffness of every pair.However, clearance is formed among teeth pair caused by wear in gears, which are run in low speed and heavy load, thus meshing point will deviate from its normal involute location.So, it is necessary to consider time-varying meshing stiffness, the sizes of the two clearances and the contact sequence of the two teeth pairs.Take the contact ratio 1< ε < 2, for instance, as shown in Figure 1.The total transmission load equals the sum of load carried by two adjacent teeth pairs, that is: Where W is total normal load per unit width, W 1 is normal load shared by the first gear pair, W 2 is normal load shared by the second gear pair.As shown in Figure 2, In case of e t1 <e t2 , the first gear pair contacts firstly, e t1 and e t2 are composite error of gear pair 1 and 2, where the direction along gear surface dent is positive.δ 1 and δ 2 are elastic deformation of the two gear pairs.k 1 (t) and k 2 (t) are meshing stiffness per unit width.
According to the principle that the load is proportional to the deformation, W 1 and W 2 can be established as: ( ) ( ) Load distribution coefficient β can be obtained by equations ( 1) to (3), that is: Where Δ=e t1 -e t2 is geometric clearance while meshing.Because of Δ<0, the clearance is occurred in the second gear pair.
The distribution load between the two gear pairs are deduced respectively as: As shown in Figure 3, the second gear pair mesh firstly, the clearance is occurred in the first gear pair because Δ is more than 0. Load distribution coefficient β is: The distribution load between the two gear pairs are deduced respectively as: Considering the clearance caused by gear wear, mathematical expression of clearance Δ can be deduced as follows: When the meshing point is at section B 1 C, clearance Δ is: When the meshing point is at section DB 2 , clearance Δ is: Where y is the distance between contact point i and pitch point i' along the direction of meshing line.h(y) is wear depth per unit width.

A numerical prediction model of wear in spur gear
Under high contact pressure, relative rolling and sliding motion occurs at the surfaces of meshing teeth.Most gear sets have oil and grease lubrication; however, lubrication condition is boundary or mixed lubrication.It indicates that meshing surface can't be fully separated by lubricant and occurs directly metal rubbing.Thus, material transfer and peeling metal will be generated at teeth surface under the low speed and heavy-duty condition, which is called as gear wear.
A generalized wear equation, which is called Archard model, is usually used to predict gear wear: Where V is wear volume loss of material, s is sliding distance of meshing gear teeth, W is the applied normal load, H is material hardness of observed surface.K is dimensionless wear coefficient, which is related to lubrication condition and wear mechanism.
Flodin [9] proposed a mild wear model in spur gear.Supposing k is constant, and surface pressure p i and sliding velocity v i remain unchanged in a very short time, the wear depth can be expressed as: Where h i,n-1 is wear depth of mesh point i after n-1 wear cycles, n is present wear cycles, Δt is time increment, N is running revolutions at every interval, during which the profile of gear teeth will not be updated.j is the point within Hertzian contact radius, p i,j is surface pressure of point j, v i,j is sliding velocity of point j.
Surface pressures of contact point i and j within contact radius can be calculated by the Winkler model [12,17].(0.59 E ) ( ) ( ) Where, j is the meshing point, a is contact radius, x i,j is distance from point i in contact radius to meshing point j, E* is equivalent elastic modulus, R is equivalent radius of two rotating cylinders.The computation expressions of indirect variable in equation ( 12) are: Where W represents W 1 or W 2 , W 1 is shared load by the first teeth pair, W 2 is shared load by the second teeth pair, which can be obtained by equations ( 5) and ( 7), E 1 , E 2 , μ 1 , μ 2 are respectively elastic modulus and poisson ratio of pinion and gear, R b1 , R b2 are radiuses of base circle, β 1i , β 2i are pressure angle of the mesh point i, Φ is engagement angle.
The sliding velocity v i of meshing point i can be determined by gear mesh theory: (18) The flow chart for wear calculation is provide in Figure 4.

Artificial neural network
Because of the capabilities in pattern classification and function approximation, artificial neural network models are being one of most mature and widely used network models.A simple three-layer BP network can approach any complicated nonlinear function relations.
where δ is S-type logarithmic function, f is linear function.

Kriging model
On the basis of the kriging method, the stochastic response model is denoted as two parts.The first part is the parametric regression model F(β), x.The second part is the non-parametric stochastic process model ( ) z x .The response is defined as: Where x is random variable and h is stochastic response, both of them are assumed to satisfy the normalization conditions.β is regression parameter.( ) f x is polynomial equation of x. ( ) z x is random process and assumed to have mean zero and covariance is: cov ( ), ( ) ( , , ) Where R is correlation function with 7 common forms, θ k is the key parameter of the function.
To keep the predictor unbiased, The Kriging approximations are defined as: Where is the vector of correlation between untried points and design points.R m×m is correlation matrix between design points.H is response vector of design points.F is design matrix of design points.
β * is deduced according to the equation below: Where R k i j ( , , ) θ x x is spatial correlation function of any two design points (

, m i j s s ∈
).And the generalized correlation function is Gaussian function.

R
x R k i j ( , , ) θ x x is assumed to be assembled into the correlation matrix R: The optimal coefficients θ * of the correlation function is derived by maximum likelihood estimation method: Where | R | is the determinant of R.
If the form of correlation function, the order of polynomial and design samples are given, β * can be calculated by equation (26).And h(x) ˆ can be obtained by equation (25) when untried points x are given.

Reliability analysis method based on the Edgeworth series and four moment techniques
Gear wear accumulates gradually along with time.System failure is said to have occurred due to wear over the threshold value.As the wear speed and threshold of driving wheel are different from driven wheel, the reliability model of the driving and driven wheel need to be built respectively.Taking driving wheel as an example, according to the safety criterion that maximum wear depth on the tooth profile are not allowed to exceed specified clearance, the limited state function of driving wheel is defined as: Where (,) i h xt is wear depth of meshing point i at any time t, sni E is lower deviation of tooth thickness and maximum clearance allowed.
The failure probability of driving wheel is formulated as follows: p p p =P ( , ) 0 =F (0) The reliability of driving wheel can be expressed by equation (32): According to equation (32), the expression of limited state function and probability distribution of original variables are necessary precondition to calculate reliability.However, both of them can't be obtained due to lack of the statistical data and complex structure.To solve the problem, the arbitrary distribution function of the standard random variable is approximately expressed by the standard normal distribution function using the Edgeworth series [25] as follows: where ϕ(•) is the standard normal probability density function, and H i (y) is the Hermite polynomial, σ g θ g η g are the variance, the third moment, and the fourth moment of the limited state function, which can be expressed by the corresponding moments of the random variables as follows based on the four moment method: Where Var(X) is the variance matrix that includes all variances and covariances of the random parameters.C 3 (X),C 4 (X) are the third and fourth moment matrices comprised by all the third and fourth central moments of each random parameter.Var(X), C 3 (X), C 4 (X) can be figured out by parameter estimation method.The Kronecker power is (•

sciENcE aNd tEchNology
As the time-varying stiffness and the distribution load are both considered in present paper, wear depths at the conversion position of single teeth-meshing area and double teeth-meshing area are fluctuating slightly and don't appear sudden change while this phenomenon has appeared in reference [2] where the load is distributed equally.The wears are varying from pitch to addendum but drastically from root to pitch.
Based on the gear parameters in reference [9], wear depth with the method in the present paper are compared with Foldin's.Considering distributed load between two gear pairs, wear depths of the pinion with the two methods are shown in figure 8, and the results are very close.

., , )
i n h x x x t is treated as the output parameter of network.Three-layer BP neutral network models with a 6-13-1 form are structured when sample sizes are 50,100 and 500 respectively (see Appendix A).Log-sigmoid function is selected as the transfer function of hidden layer.Purelin function is selected as the transfer function of output layer.Trainlm func- To determine the suitable sample size, both of the accuracy and goodness of fit needed to be tested for the three network models with different samples.The accuracy means comparison of gradual wear process between exact and surrogate models when original stochastic variables { } T 12 , ,..., n x x x are deterministic and equal to the mean values.The goodness of fit means comparison of wear depths between the two models with another set of samples when time is deterministic.
The accuracy of driving and driven wheels with ANN model are figured as follows: As shown in figures 9 and 10, the accuracy is low with 50 samples whereas high with 100 and 500 samples.And the latter two make accurate prediction on the developing trend of wear depth.
The accuracy of driving and driven wheels with Kriging model are figured as follows: As shown in figures 11 and 12, all of three models have high accuracy and make accurate prediction on the developing trend.It shows that Kriging method needs less samples than ANN method under the same accuracy.
Figures 13 and 14 shows goodness of fit based on ANN and Kriging surrogate models with 100 samples.Results of the other two are shown in Table 2.
The goodness of fit for ANN models with 100 and 500 samples can reach more than 0.9999.But the 50-samples model has a low goodness of fit so that can't be used.The goodness of fit for Kriging models with all samples can reach more than 0.96.In a comprehensive view, the suitable sample size for ANN model is more than 100, for Kriging more than 50.

Reliability curve
The permitted maximum clearances of driving and driven wheels are -142μm and -143μm respectively.On the basis of ANN surrogate models, combining with Edgeworth series and four moment method, the time-dependent reliability curve can be obtained under the condition that the four moments of original variables are known but the probability distribution unknown.Reliability curve can be also obtained with Kriging-based MCS method.Supposing probability distribution are known is convenient to make comparison of the ANN-Edgeworth Series-Four moment method and Kriging-MCS method, which are shown in Figures 15 and 16.
Reliability remains to be 1 in the beginning for a long time.After 20 hours, reliability decreases gradually and becomes 0 at 50 hours.It is important to note that the applied load in the present paper is much larger than that in practice to save time.And product life under normal load are longer than test load.However, the  present paper is focus on studying gradually law of wear reliability and choosing the suitable method to solve complex problem.According to the results with two methods, ANN-Edgeworth Series-Four moment method is a little bigger than Kriging-MCS method.That is to say, the difference is only a little.Supposing the time which is spent running the exact wear model once is 20 minutes, consuming time of every method for reliability is listed in Table 3.
The time for building up original samples account for absolute proportion, about more than 99% of the total.If more original samples the method needs, more time it will consume.As the direct-MCS method needs a large number of samples, the cost of time is as high as 3.8 years, which is impossible to complete.The former two only need tens of hours, which has a high efficiency.

7.conclusions
The present paper studied dynamic evolution law of mechanical reliability caused by wear.A numerical wear model has been established considering dynamic distribution load of gear tooth.Wear depths and gradual law has been obtained.Wear random process has been introduced.The problem of time-dependent reliability analysis with complex structural has been solved.The main conclusions are: (1) Load and wear are coupled with each other.The dynamic characteristic of gear system is affected by variation of geometrical shape and dimension parameter caused by mild wear and dynamic load will lead to more serious wear.Thus, time-varying stiffness and dynamic distribution load should be both considered into wear model.As a result, sudden change hasn't appeared at the conversion point of single and double meshing area.
(2) Wear calculation spends a lot of time so that direct MCS method is not suitable for reliability analysis.It necessary to develop a surrogate model to replace exact model.Both of ANN and Kriging models have a high accuracy and goodness of fit.Kriging model needs less samples than ANN.
(3) Under the condition that the probability functions of random variables are unknown, Edgeworth series and four moment methods can be used to calculate reliability.
(4) Both of ANN-Edgeworth series-Four moment method and Kriging-MCS method are suitable for reliability with complex structure and have little difference.

Fig. 2 .Fig. 3 .
Fig. 2. Force diagram of meshing gears (condition 1) ..., n x xx and time t is treated as the input parameter of network, and wear depth response 1 2 ( , ,..., , ) i n h x x x t is treated as the output parameter of network.The expression of nonlinear function is obtained after training:

Fig. 6 .
Fig. 6.Wear depth of the pinion after n wear cycles

2 (
t is treated as the input parameter of network, and wear depth response 1 , ,..

Fig. 13 .Fig. 14 .
Fig. 13.Approximate values with Kriging method and Actual values of driving wheel

Table 2 .
Goodness of fit for ANN and Kriging models with different sample sizes

Table 3 .
Consuming time of the three methods

Table 1 .
Original samples for construct ANN and Kriging models (50samples)

Table 2 .
ANN network parameters of driving wheel

Table 3 .
ANN network parameters of driving wheel