# Stress Analysis of Gear Meshing Impact Based on SPH Method

Authors: Rahmatjan Imin; Mamtimin Geni
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Based on the kinetic equations of the gear mesh impact, SPH discrete equations were established. Numerical simulation was carried out on the meshing impact process of the gear, and stress and strain of each discrete point were obtained. After data processing, stress propagation was calculated, which shows stress distribution on tooth-profile surface. It is concluded that the stress concentrate mainly occurs in the pitch circle. The paper provides an effective new numerical simulation algorithm to gear mechanical properties analysis.


---

## Body

## 1. Introduction

In researches studying gear transmissions which bear heavy load and are responsible for dynamic transfer, analytical methods or simple numerical methods usually do a great deal of simplifications and could not properly reflect the actual situation [1]. With the development of computer technology and the emergence of new numerical algorithms, computer simulation has shown a lot of advantages in gear transmission analysis. For example, with elastic contact finite element method, it is possible to analyze displacement and stress change in gear contact zone [2, 3]. Analytical methods and finite element method could be used to solve surface deformation and load distribution problem on tooth surface of meshing gear and also could be used to analyze effects of tooth surface deformation on the gear transmission [4] ANSYS/LSDYNA software can simulate gears meshing impact [5]. However, with methods based on the meshing (finite difference method, finite element method, boundary element method, etc.), it is hard to track gear transmission process, especially when calculating large deformation and crack propagation because of mesh tangling and distortion. In this case, remeshing is needed in calculation process which is computationally expensive and affects calculation accuracy.Smoothed particle hydrodynamics (known as SPH) is one of the meshless methods developed in recent years. The basic idea of SPH method is to discretize the gear model into a set of particles, each of them having properties like mass, position, velocity, and density [6, 7]. During computation, particles are linked by a kernel function. The process does not need meshing, is easy to program, can track gear transmission process every time, and guarantees calculation accuracy.In this paper, meshing impact of a pair of gear teeth is analyzed using SPH method. Dynamic numerical simulation was conducted regarding meshing impact process in different rotation speeds, and the stress distribution on the tooth profile and changes was analyzed. The simulation results meet actual meshing rules. The results of this study show that SPH method has very good potential in analyzing impact of gear meshing.


## 2. SPH Description of Correct Gear Meshing and Contact

### 2.1. Fundamentals of SPH Method

In the SPH method, the first step is to approximate a function and its gradient using integration of kernel function based on interpolation theories. Continuous partial differential equation is transformed into integral equations. At the second step continuous forms of integral equations are discretized into discrete equations using particle approximation method [6].For any functionf
(
x
), the value at a point x can be approximated as
(1)
〈
f
(
x
)
〉
=
∫
Ω
f
(
x
′
)
W
(
x
-
x
′
,
h
)
d
x
′
,
where 〈
〉 denotes kernel approximation of the function, x
′ is position vector, and W
(
x
-
x
′
,
h
) is kernel function which depends on distance |
x
-
x
′
| and smoothed length h.The most frequently used kernel function in SPH method is the cubic B-spline kernel [8]. It has the form as shown in (2); Figure 1 shows that choosing different smoothed length (
h
=
1 and h
=
0.5, resp., and x
′
=
0
) affects influence radius and shape of cubic B-spline kernel. Consider
(2)
W
(
x
-
x
′
,
h
)
=
α
{
2
3
-
R
2
+
1
2
R
3
0
≤
R
<
1
1
6
(
2
-
R
)
3
1
≤
R
<
2
0
R
≥
2
,
where α is equal to 1
/
h
,
15
/
7
π
h
2, and 3
/
2
π
h
3 in one-, two-, and three-dimensional cases, respectively; R is the relative distance between two points (particles) at points x and x
′; R
=
r
/
h
=
|
x
-
x
′
|
/
h where r is the distance between the two points.Figure 1
Figure of cubic B-spline kernel function in 1D case.Using integration by parts, Gauss theorem, and property of kernel function, gradient∇
f
(
x
) can be approximated as
(3)
〈
∇
f
(
x
)
〉
=
∫
Ω
f
(
x
′
)
∇
W
(
x
-
x
′
,
h
)
d
x
′
.
Function f
(
x
) and its gradient ∇
f
(
x
), continuous equations (1) and (3) can be discretized using SPH method as
(4)
〈
f
(
x
i
)
〉
=
∑
j
=
1
N
m
j
ρ
j
f
(
x
j
)
W
i
j
〈
∇
f
(
x
i
)
〉
=
∑
j
=
1
N
m
j
ρ
j
f
(
x
j
)
∇
i
W
i
j
,
where i and j are the particle indices, m
j and ρ
j are mass and density of the particles j
,
N is the number of particles in the influence domain of particle i
,
W
i
j
=
W
(
x
i
-
x
j
,
h
)
=
W
(
r
i
j
,
h
)
,
∇
i
W
i
j
=
(
x
i
j
/
r
i
j
)
(
∂
W
i
j
/
∂
r
i
j
)
,
r
i
j
=
x
i
-
x
j. In the above equation, the sum is calculated only in influence area of W and is not in the entire computational area; the influence area is determined by the radius of kernel function [9], as shown in Figure 2 in 2D case.Figure 2
SPH particle kernel approximation.


### 2.2. SPH Discretization of Impact Gear Meshing Equations Based on Elastic Dynamics Theory

Based on continuum mechanics and linear elastic dynamics theory, the kinetic equation and constitutive equation of gear meshing impact can be written as follows [5].Conservation of mass:(5)
∂
v
β
∂
x
β
=
-
1
ρ
d
ρ
d
t
.
Conservation of momentum:
(6)
d
v
α
d
t
=
1
ρ
∂
σ
α
β
∂
x
β
+
f
α
.
Strain-displacement equation:
(7)
ε
α
β
=
1
2
(
∂
v
α
∂
x
β
+
∂
v
β
∂
x
α
)
.
Constitutive equation:
(8)
σ
α
β
=
E
1
+
υ
(
υ
1
-
2
υ
ε
α
β
δ
α
β
+
ε
α
β
)
.
where v
,
ρ
,
σ
,
f
,
E, and υ are the velocity, density, stress, external force, elastic modulus, and Poisson ratio, respectively; δ
α
β
=
1 if α
=
β and 0 if α
≠
β.Considering the artificial viscosity effectΠ
i
j, using (4), a set of SPH particle discreet equations such as conservation of mass, conservation of momentum, and strain-displacement equation could be obtained:
(9)
d
ρ
i
d
t
=
∑
j
=
1
N
m
j
(
v
i
β
-
v
j
β
)
∂
W
i
j
∂
x
i
β
d
v
i
α
d
t
=
∑
j
=
1
N
m
j
(
σ
i
α
β
ρ
i
2
+
σ
j
α
β
ρ
j
2
+
Π
i
j
)
i
∂
W
i
j
∂
x
i
β
+
f
i
α
ε
i
α
β
=
1
2
∑
j
=
1
N
m
j
ρ
j
(
v
j
i
α
∂
W
i
j
∂
x
i
β
+
v
j
i
β
∂
W
i
j
∂
x
i
α
)
,
where
(10)
Π
i
j
=
{
-
α
c
-
i
j
μ
i
j
+
β
μ
i
j
ρ
-
i
j
μ
i
j
<
0
,
0
μ
i
j
≥
0
,
μ
i
j
=
h
(
v
i
-
v
j
)
·
(
x
i
-
x
j
)
(
x
i
-
x
j
)
2
+
ε
h
2
,
c
-
i
j
=
1
2
(
c
i
+
c
j
)
,
ρ
-
i
j
=
1
2
(
ρ
i
+
ρ
j
)
,
where c
i and c
j denotes the sound speeds at particles i and j respectively, x
i and x
j are coordinates of particles, v
i and v
j are velocities of particles, h is smoothed length of kernel function, and α
,
β
,
ε are adjustable parameters [10].


## 2.1. Fundamentals of SPH Method

In the SPH method, the first step is to approximate a function and its gradient using integration of kernel function based on interpolation theories. Continuous partial differential equation is transformed into integral equations. At the second step continuous forms of integral equations are discretized into discrete equations using particle approximation method [6].For any functionf
(
x
), the value at a point x can be approximated as
(1)
〈
f
(
x
)
〉
=
∫
Ω
f
(
x
′
)
W
(
x
-
x
′
,
h
)
d
x
′
,
where 〈
〉 denotes kernel approximation of the function, x
′ is position vector, and W
(
x
-
x
′
,
h
) is kernel function which depends on distance |
x
-
x
′
| and smoothed length h.The most frequently used kernel function in SPH method is the cubic B-spline kernel [8]. It has the form as shown in (2); Figure 1 shows that choosing different smoothed length (
h
=
1 and h
=
0.5, resp., and x
′
=
0
) affects influence radius and shape of cubic B-spline kernel. Consider
(2)
W
(
x
-
x
′
,
h
)
=
α
{
2
3
-
R
2
+
1
2
R
3
0
≤
R
<
1
1
6
(
2
-
R
)
3
1
≤
R
<
2
0
R
≥
2
,
where α is equal to 1
/
h
,
15
/
7
π
h
2, and 3
/
2
π
h
3 in one-, two-, and three-dimensional cases, respectively; R is the relative distance between two points (particles) at points x and x
′; R
=
r
/
h
=
|
x
-
x
′
|
/
h where r is the distance between the two points.Figure 1
Figure of cubic B-spline kernel function in 1D case.Using integration by parts, Gauss theorem, and property of kernel function, gradient∇
f
(
x
) can be approximated as
(3)
〈
∇
f
(
x
)
〉
=
∫
Ω
f
(
x
′
)
∇
W
(
x
-
x
′
,
h
)
d
x
′
.
Function f
(
x
) and its gradient ∇
f
(
x
), continuous equations (1) and (3) can be discretized using SPH method as
(4)
〈
f
(
x
i
)
〉
=
∑
j
=
1
N
m
j
ρ
j
f
(
x
j
)
W
i
j
〈
∇
f
(
x
i
)
〉
=
∑
j
=
1
N
m
j
ρ
j
f
(
x
j
)
∇
i
W
i
j
,
where i and j are the particle indices, m
j and ρ
j are mass and density of the particles j
,
N is the number of particles in the influence domain of particle i
,
W
i
j
=
W
(
x
i
-
x
j
,
h
)
=
W
(
r
i
j
,
h
)
,
∇
i
W
i
j
=
(
x
i
j
/
r
i
j
)
(
∂
W
i
j
/
∂
r
i
j
)
,
r
i
j
=
x
i
-
x
j. In the above equation, the sum is calculated only in influence area of W and is not in the entire computational area; the influence area is determined by the radius of kernel function [9], as shown in Figure 2 in 2D case.Figure 2
SPH particle kernel approximation.


## 2.2. SPH Discretization of Impact Gear Meshing Equations Based on Elastic Dynamics Theory

Based on continuum mechanics and linear elastic dynamics theory, the kinetic equation and constitutive equation of gear meshing impact can be written as follows [5].Conservation of mass:(5)
∂
v
β
∂
x
β
=
-
1
ρ
d
ρ
d
t
.
Conservation of momentum:
(6)
d
v
α
d
t
=
1
ρ
∂
σ
α
β
∂
x
β
+
f
α
.
Strain-displacement equation:
(7)
ε
α
β
=
1
2
(
∂
v
α
∂
x
β
+
∂
v
β
∂
x
α
)
.
Constitutive equation:
(8)
σ
α
β
=
E
1
+
υ
(
υ
1
-
2
υ
ε
α
β
δ
α
β
+
ε
α
β
)
.
where v
,
ρ
,
σ
,
f
,
E, and υ are the velocity, density, stress, external force, elastic modulus, and Poisson ratio, respectively; δ
α
β
=
1 if α
=
β and 0 if α
≠
β.Considering the artificial viscosity effectΠ
i
j, using (4), a set of SPH particle discreet equations such as conservation of mass, conservation of momentum, and strain-displacement equation could be obtained:
(9)
d
ρ
i
d
t
=
∑
j
=
1
N
m
j
(
v
i
β
-
v
j
β
)
∂
W
i
j
∂
x
i
β
d
v
i
α
d
t
=
∑
j
=
1
N
m
j
(
σ
i
α
β
ρ
i
2
+
σ
j
α
β
ρ
j
2
+
Π
i
j
)
i
∂
W
i
j
∂
x
i
β
+
f
i
α
ε
i
α
β
=
1
2
∑
j
=
1
N
m
j
ρ
j
(
v
j
i
α
∂
W
i
j
∂
x
i
β
+
v
j
i
β
∂
W
i
j
∂
x
i
α
)
,
where
(10)
Π
i
j
=
{
-
α
c
-
i
j
μ
i
j
+
β
μ
i
j
ρ
-
i
j
μ
i
j
<
0
,
0
μ
i
j
≥
0
,
μ
i
j
=
h
(
v
i
-
v
j
)
·
(
x
i
-
x
j
)
(
x
i
-
x
j
)
2
+
ε
h
2
,
c
-
i
j
=
1
2
(
c
i
+
c
j
)
,
ρ
-
i
j
=
1
2
(
ρ
i
+
ρ
j
)
,
where c
i and c
j denotes the sound speeds at particles i and j respectively, x
i and x
j are coordinates of particles, v
i and v
j are velocities of particles, h is smoothed length of kernel function, and α
,
β
,
ε are adjustable parameters [10].


## 3. Establishing SPH Discrete Particle Model of Two Meshing Gears and Setting Initial, Boundary Conditions

### 3.1. SPH Model of Two Correct Meshing Gears

Establishing SPH discrete particle model is to discretize the gear model into a set of particles. In this paper, by analyzing every curve of involute gear, tooth profile curve equation and tooth root transition curve equation were determined; gear 3D mathematical discrete particle model was established and programmed [11]. Model data was visualized using open source software RASMOL, as shown in Figure 3.Discrete model of gear contact.
(a)
Two correct meshing gears
(b)
Enlarged partial viewAfter discretizing the gear model into discrete particles, it is difficult to introduce into whole body parameters like torque load (including resistance torque) and related conditions of fixed axis rotation; at the same time it will add additional computation work too. To save the modeling and calculation work, our numerical simulation considers only a pair of meshing impact gears as shown in Figure4.Figure 4
Initial and boundary conditions of the gear meshing impact.


### 3.2. Initial Conditions

At the beginning of the calculation, the driving gear was rotating with certain initial angler velocity and impact to the driven gear which in static state(11)
v
i
(
x
,
0
)
=
ω
×
r
i
x
∈
Driving
gear
v
i
(
x
,
0
)
=
0
x
∈
Driven
gear
,
where ω is angler velocity of the driving gear which rotates around gear center C; r
i is position vector of particle i.


### 3.3. Boundary Conditions

During calculation process particles on the driven gear root were fixed and kept motionless in the whole calculation process:(12)
v
i
(
x
,
t
)
=
0
x
∈
S
,
0
≤
t
<
+
∞
.


### 3.4. Interface Coupling

In particular, when a pair of gears contacts each other’s, interaction of particles which belongs to different gears will happen near the contact surface and the particles will participate in calculation. If the influence radiuses of all particles are the same as the original one, the interaction will occur between particles of two gears at contact area earlier, which will cause the two gears not directly contacting each other on the surface and still remains the gap during the meshing process as shown in Figure5(a), so this will lead to errors on gear meshing and contact.Comparison of selecting different influence radius.
(a)
Original method
(b)
New methodIn general, the two problems occur when using the general SPH influence radius as mentioned above, one is correct gear meshing and the other one is correct contact with the tooth surface during the dynamic contact and meshing process. To describe dynamic contact phenomenon of gear meshing and impact process by using SPH method, defining different kernel radius method is used for correct coupling of the interaction and contact areas of gears. In this method, it not needed to define any special contact pair or contact area due to SPH kernel function characteristics, so the choice of the influence radius of kernel is very important. The influence radius affects the accuracy and efficiency of the calculation directly.From the other point of view, different behaviors can be easily obtained by tuning the kernel effective radiusk
h. A larger value of k
h will create averaged interactions of each particle in the body, so the body will be bonded into more strong continuous body such as evaluation of particles in the inside body of independent gear. And a small value of k will create very local interactions of each particle, so the body will be divided more easily into discrete pieces such as evaluation of particles near the two gears meshing at contact area. In this paper due to the above unique characteristics of SPH kernel, different kernel radiuses are defined by choosing different smoothed length as shown in Figure 6 and expressed as (13) to evaluate two gears correctly meshing and impact problems.Figure 6
Selection of influence radius.Consider(13)
k
h
i
j
=
k
[
h
11
h
12
h
21
h
22
]
,
where k
h
i
i
,
i
=
1,2 is larger kernel radius which is used to calculate and evaluate the interaction of particles inside gear and represent strong continuous bodies itself, the coefficient k is determined by corresponding kernel function which is used in calculation, k
=
2 in cubic B-spline, smoothed length equals to 1.2 times of the distance between two particles [12]. Then the k
h
i
j
,
i
,
j
=
1,2
(
i
≠
j
) is smaller kernel radius which is used to calculate and evaluate the interaction between two gears at contact surface where two gears are not continuous and separated each other; generally choose an appropriate smoothed length which keeps kernel radius equal to one time of the distance between two particles. The errors caused by selection of influence radius were complemented by corrective SPH method [13].The results were shown in Figure5(b); the two gears correctly meshing and directly contacting each other’s after different kernel radiuses are used which will effectively avoid interaction of noncontacting particles participating in calculation.


### 3.5. Material Properties and Gear Parameters

Material properties and gear parameters used in numerical simulation are listed as follows: elastic modulusE
=
206 GPa, Poisson ratio: υ
=
0.3, density: ρ
=
7870 kg/m3, gear module: M
=
0.04 mm, numbers of tooth: z
1
=
21
,
z
2
=
42, and pressure angle: α
=
2
0
°.


## 3.1. SPH Model of Two Correct Meshing Gears

Establishing SPH discrete particle model is to discretize the gear model into a set of particles. In this paper, by analyzing every curve of involute gear, tooth profile curve equation and tooth root transition curve equation were determined; gear 3D mathematical discrete particle model was established and programmed [11]. Model data was visualized using open source software RASMOL, as shown in Figure 3.Discrete model of gear contact.
(a)
Two correct meshing gears
(b)
Enlarged partial viewAfter discretizing the gear model into discrete particles, it is difficult to introduce into whole body parameters like torque load (including resistance torque) and related conditions of fixed axis rotation; at the same time it will add additional computation work too. To save the modeling and calculation work, our numerical simulation considers only a pair of meshing impact gears as shown in Figure4.Figure 4
Initial and boundary conditions of the gear meshing impact.


## 3.2. Initial Conditions

At the beginning of the calculation, the driving gear was rotating with certain initial angler velocity and impact to the driven gear which in static state(11)
v
i
(
x
,
0
)
=
ω
×
r
i
x
∈
Driving
gear
v
i
(
x
,
0
)
=
0
x
∈
Driven
gear
,
where ω is angler velocity of the driving gear which rotates around gear center C; r
i is position vector of particle i.


## 3.3. Boundary Conditions

During calculation process particles on the driven gear root were fixed and kept motionless in the whole calculation process:(12)
v
i
(
x
,
t
)
=
0
x
∈
S
,
0
≤
t
<
+
∞
.


## 3.4. Interface Coupling

In particular, when a pair of gears contacts each other’s, interaction of particles which belongs to different gears will happen near the contact surface and the particles will participate in calculation. If the influence radiuses of all particles are the same as the original one, the interaction will occur between particles of two gears at contact area earlier, which will cause the two gears not directly contacting each other on the surface and still remains the gap during the meshing process as shown in Figure5(a), so this will lead to errors on gear meshing and contact.Comparison of selecting different influence radius.
(a)
Original method
(b)
New methodIn general, the two problems occur when using the general SPH influence radius as mentioned above, one is correct gear meshing and the other one is correct contact with the tooth surface during the dynamic contact and meshing process. To describe dynamic contact phenomenon of gear meshing and impact process by using SPH method, defining different kernel radius method is used for correct coupling of the interaction and contact areas of gears. In this method, it not needed to define any special contact pair or contact area due to SPH kernel function characteristics, so the choice of the influence radius of kernel is very important. The influence radius affects the accuracy and efficiency of the calculation directly.From the other point of view, different behaviors can be easily obtained by tuning the kernel effective radiusk
h. A larger value of k
h will create averaged interactions of each particle in the body, so the body will be bonded into more strong continuous body such as evaluation of particles in the inside body of independent gear. And a small value of k will create very local interactions of each particle, so the body will be divided more easily into discrete pieces such as evaluation of particles near the two gears meshing at contact area. In this paper due to the above unique characteristics of SPH kernel, different kernel radiuses are defined by choosing different smoothed length as shown in Figure 6 and expressed as (13) to evaluate two gears correctly meshing and impact problems.Figure 6
Selection of influence radius.Consider(13)
k
h
i
j
=
k
[
h
11
h
12
h
21
h
22
]
,
where k
h
i
i
,
i
=
1,2 is larger kernel radius which is used to calculate and evaluate the interaction of particles inside gear and represent strong continuous bodies itself, the coefficient k is determined by corresponding kernel function which is used in calculation, k
=
2 in cubic B-spline, smoothed length equals to 1.2 times of the distance between two particles [12]. Then the k
h
i
j
,
i
,
j
=
1,2
(
i
≠
j
) is smaller kernel radius which is used to calculate and evaluate the interaction between two gears at contact surface where two gears are not continuous and separated each other; generally choose an appropriate smoothed length which keeps kernel radius equal to one time of the distance between two particles. The errors caused by selection of influence radius were complemented by corrective SPH method [13].The results were shown in Figure5(b); the two gears correctly meshing and directly contacting each other’s after different kernel radiuses are used which will effectively avoid interaction of noncontacting particles participating in calculation.


## 3.5. Material Properties and Gear Parameters

Material properties and gear parameters used in numerical simulation are listed as follows: elastic modulusE
=
206 GPa, Poisson ratio: υ
=
0.3, density: ρ
=
7870 kg/m3, gear module: M
=
0.04 mm, numbers of tooth: z
1
=
21
,
z
2
=
42, and pressure angle: α
=
2
0
°.


## 4. Numerical Simulation and Analyzing Results

In numerical simulation three different driving gear angular velocities,ω
=
60 r/min, ω
=
120 r/min, and ω
=
240 r/min, were used in calculation and analysis.Numerical calculation results were shown in Figure7. Maximum stress and maximum equivalent stress change with time after impact of two meshing gears. We can know from the figure that at the beginning of meshing, the biggest impact stress increases with the extension of engagement time. Meanwhile, enlarging impact speed will increase impact stress gradually.Maximum stress change with time fromt
=
0 μsec. to t
=
20 μsec.
(a)
Stressσ
x
(b)
Equivalent stressFigure8 shows selected consecutive figure of the equivalent stress contours. These stress variations in different time show the propagation of elastic waves in the meshing gears.Propagation of elastic waves represented by equivalent stress distribution.
(a)
t
=
2 μs
(b)
t
=
4 μs
(c)
t
=
8 μsFigure9 shows comparison of stress σ
x and equivalent stress σ
e at time t
=
1
0 μs. It could be seen that impact stress increases with the increase of impact velocity; increase of driving gear initial impact velocity causes propagation of tooth impact stress.Comparison of stress att
=
10 μsec.
(a)
ω
=
60 r/min
(b)
ω
=
12
0 r/min
(c)
ω
=
240
0 r/min
(d)
ω
=
60 r/min equivalent stress
(e)
ω
=
12
0 r/min equivalent stress
(f)
ω
=
24
0 r/min equivalent stressFigure10 shows comparison and calculation result of stress on the middle line along the contact faces of the driven gear and the driven gear’s tooth profile at t
=
1
0 μs in different rotation speed. The results indicate that maximum stresses appear near the pitch circle; when the impact velocity becomes bigger, the stress also becomes bigger.Stress on the middle line along the contact faces of gear’s tooth att
=
10 μsec.
(a)
Stressσ
x in the driving gear
(b)
Stressσ
x in the driven gear
(c)
Equivalent stressσ
e in the driving gear
(d)
Equivalent stressσ
e in the driven gear


## 5. Conclusion and Discussion

In this research, a SPH discrete model of impact of gear meshing was established, and automatic modeling program was developed. A simulation and analysis software of gear meshing which is capable of simulating gear mesh impact problem of different gear parameters, material property, and rotation velocity was programmed.The results of numerical simulation comply with real gear behavior, which verified the correctness of the model and method. A new numerical simulation algorithm based on SPH method is provided for design and optimization of gear transmission.



---

*Source: 328216-2014-04-17.xml*