Abstract
A method to simulate residual stress field of planetary gear is proposed. In this method, the finite element model of planetary gear is established and divided to tooth zone and profile zone, whose different temperature field is set. The gear's residual stress simulation is realized by the thermal compression stress generated by the temperature difference. Based on the simulation, the finite element model of planetary gear train is established, the dynamic meshing process is simulated, and influence of residual stress on equivalent stress of addendum, pitch circle, and dedendum of internal and external meshing planetary gear tooth profile is analyzed, according to non-linear contact theory, thermodynamic theory, and finite element theory. The results show that the equivalent stresses of planetary gear at both meshing and nonmeshing surface are significantly and differently reduced by residual stress. The study benefits fatigue cracking analysis and dynamic optimization design of planetary gear train.
1. Introduction
Residual stresses are included in the advanced design and fatigue failure analysis of components in the gear transmission system, aerospace, and nuclear and automotive industries. Residual stresses have an effect on many aspects of gears but are normally one of the largest unknown stresses, difficult to measure and to estimate theoretically. Residual stresses may play a significant role in analyzing dynamic behavior of gear transmission systems.
A number of studies have been performed on the dynamic behavior of gear transmission systems. Cooley and Parker [1] demonstrated unusual gyroscopic system eigenvalue behavior observed in a lumped-parameter planetary gear model. The focus was on the eigenvalue phenomena that occur at especially high speeds rather than practical planetary gear behavior. Yuksel and Kahraman [2] employed a computational model of a planetary gear set to study the influence of surface wear on the dynamic behavior of a typical planetary gear set. The results for worn gear surfaces indicated that surface wear had an influence on off-resonance speed ranges while its influence diminishes near resonance peaks primarily due to tooth separations. In [3] the mesh phase relations of general compound planetary gears were systematically studied. The results from the study were necessary for analytical or computational multibody studies on the static or dynamic response of general compound planetary gears. In [4] the lumped parameter dynamic model of a gear system actuated by many pinions was presented in this study and strength and had traditionally been obtained by laboratory measurements. The numerical results showed that these parameters could affect the load transmission behavior of every mesh gear pair. Suitable mechanical parameters of gear systems were important for the load sharing performance of multiple mesh gear pairs which could reduce the internal stress in the gear structure and improve their fatigue life. In [5] the discrete-continuous model of a single gear transmission with nonlinear forces occurring between gear teeth was investigated. The range of application for the proposed nonlinear functions was determined. In [6] a general approach to predict the contact fatigue life of the gears in the drive-train system of a wind turbine under dynamic conditions was presented. In [7] a new approach to the analysis of the performance of purse seine fishing gear was presented using a computer simulation. The accuracy of the mathematical model used in the simulation was verified by comparing measured and calculated data. In [8] the dynamic behavior of a single stage spur gear reducer in transient regime was studied. Dynamic responses came to confirm a significant influence of the transient regime on the dynamic behavior of a gear set, particularly in the case of engine acyclic condition. In [9] the behavior of 35CrMo structural steel gears treated through laser surface-melted (LSM) treatment was investigated. The improvement in resistance of the LSM gears was due to grain refinement and solid solution hardening. In [10] the nonlinear tooth wedging behavior and its correlation with planet bearing forces by analyzing the dynamic response of an example planetary gear were studied. The major causes of tooth wedging relate directly to translational vibrations caused by gravity forces and the presence of clearance-type nonlinearities in the form of backlash and bearing clearance. In [11] an analytical model was proposed to investigate the effect of gear tooth crack on the gear mesh stiffness. The results showed that sidebands caused by the tooth crack were more sensitive than the mesh frequency and its harmonics. The developed analytical model could predict the change of gear mesh stiffness with presence of a gear tooth crack and the corresponding dynamic responses could supply some guidance to the gear condition monitoring and fault diagnosis, especially for the gear tooth crack at early stage. In [12] dynamic and three-dimensional finite element (FE) analytical models of cracked gears were established, and the dynamic characteristics of the gear body were investigated when tooth cracks in the gear appeared. In [13] a systematic analysis of the dynamic behavior of a single degree-of-freedom (SDOF) spur gear system with and without nonlinear suspension was performed. The results presented in this study provided an understanding of the operating conditions under which undesirable dynamic motion took place in a spur gear system and therefore served as a useful source of reference for engineers in designing and controlling such systems. In [14] the dynamic responses of a planetary gear were analyzed when component gears had time-varying pressure angles and contact ratios caused by bearing deformations. The effects of bearing stiffness on the pressure angles and contact ratios were also analyzed. In [15] the dynamic modeling of wheeled planetary rovers operating on rough terrain was dealt with. Different implementations, with different degrees of complexity, were presented and compared with each other so that the user could simulate the dynamics of the rover making a tradeoff between simulation accuracy and computer time. In [16] the distinctive wave vibration of a ring gear affected by mesh effect was investigated based on the inherent symmetry of spur planetary gears. The main results were demonstrated with FE examples and comparisons with the existing literature. In [17] an original lumped parameter model of planetary gears, which accounted for planet position errors and simulated their contribution to the quasistatic and dynamic load sharing amongst the planets, was presented. In [18] nonlinear oscillations of spur gear pairs including the backlash nonlinearity were studied. Dynamic system was described through the classical SDOF model in terms of dynamic transmission error. In [19] a generalized dynamic model for herringbone planetary gear train (HPGT) was developed to investigate its modal properties. The new relations between deflections of planets in translational mode were also derived in this research. In [20] a three-dimensional lumped-parameter model for a pair of helical gears was developed and shown to be equivalent to an arbitrary tooth contact force distribution. The twist stiffness fluctuated periodically with gear rotation generating fluctuating moments (shuttling) that could potentially excite vibrations. In [21] the three-dimensional nonlinear vibration of gear pairs was investigated where the nonlinearity was due to portions of gear teeth contact lines. The nonlinear dynamic response was obtained using the discretized stiffness network.
Although references about gear's dynamic mashing property are available, they mainly focused on the influences of time-varying stiffness, meshing error, and rotational speed, but seldom discussed the influence of residual stress. Few reports about simulating residual stress field of planetary gear train by FE division technique are concerned, as well as those about the influence of residual stress on meshing characteristics of planetary gear train.
In this study, FE model of planetary gear was built, realizing the simulation of planetary gear's residual stress and analyzing the influences of temperature on residual stresses of addendum, pitch circle, and dedendum of planetary gear's tooth profile. Then, FE model of planetary gear pair was established to simulate the meshing process, and equivalent stresses of tooth profile were studied with and without residual stress of planetary gear, respectively, comparing the influence of residual stress on equivalent stresses at different meshing positions of planetary gear's tooth profile.
2. Analysis of Thermal Stress
In general, internal thermal stress, which results in internal strain, is generated because of the different thermal expansions caused by temperature gradient. The Hooke relationship between stress and strain in elastic range is expressed as
where {σ} = [σ x , σ y , σ z , σ xy , σ yz , σ xz ] T is column vector of stress, [D] is coefficient of stiffness matrix, and {ε xd } is column vector of corresponding strain, which can be expressed as
where {ε} = [ε x , ε y , ε z , ε xy , ε yz , ε xz ] T and {εth} are column vectors of total strain and thermal strain, respectively. Using (1) and (2) gives the following:
Relationship between thermal strain and temperature is given as
where α x , α y , and α z are thermal expansion coefficient at x-axis, y-axis, and z-axis, respectively, and ΔT is temperature difference.
[D]−1 can be shown as
where E x , E y , and E z are Young's modulus at x-axis, y-axis, and z-axis, respectively. ν xy , ν yx , and ν yz are Poisson ratio, and G xy , G yx , and G yz are shear modulus at xy plane, yz plane and xz plane, respectively.
The relationship among stress, strain, and mechanical characteristic is shown as
where σ x , σ y , and σ z are normal stress at x-axis, y-axis, and z-axis, respectively, and ε x , ε y , and ε z are corresponding strain. σ xy , σ yz , and σ xz are shear stress at xy plane, yz plane, and xz plane, respectively, and ε xy , ε xy , and ε xz are corresponding strain.
Expression of stress can be obtained by substituting (4), (5), and (6) to (3). As is expressed in (7), thermal stress is related to temperature difference, whose variation will result in the change of stress. Accordingly, residual stress after heat treatment can be simulated by setting temperature fields in different zones of material and studying stress's variation with the change of temperature difference:
3. Method to Simulate the Residual Stress
3.1. FE Model of Planetary Gear
To mesh complicated model, FE division technique, which realizes better controlling position and density of single unit, mesh refinement at concerned zones, and endowing different types, and attribute to divided zones, is introduced to divide FE model to several simple zones. Planetary gear in Figure 1 (a) is divided into body zone and tooth profile zone by FE division technique shown in Figures 1 (b) and 1 (c).

Schematics of FE division of planetary gear. (a) Planetary gear. (b) Tooth profile of planetary gear. (c) Body of planetary gear.
3.2. Simulation Method
In this study, the compression stress is used to simulate the residual stress. The compression stress at high-temperature zone results from the thermal stress difference. The thermal stress is produced by the temperature difference between two adjacent zones. The FE model is established and planetary gear is divided into body zone and tooth profile zone (shown in Figure 1). The temperature of each zone is set using temperature field predefined function in Abaqus to introduce the temperature difference, by which the thermal stress is produced. The thermal stresses at two adjacent zones are different because of the temperature difference, and compression stress is brought at high-temperature zone.
Parameters of planetary gear for simulating residual stress are listed in Table 1. The residual stress of planetary gear is simulated by Abaqus. It is shown in Figure 2 that residual stress is produced at addendum, pitch circle and dedendum of tooth profile if temperature field at body zone is set to be different from that at tooth zone.
Parameters of planetary gear for simulating residual stress.

Residual stress simulation of planetary gear. (a) Stress of planetary gear. (b) Residual stress at dedendum of planetary gear.
3.3. Influence of Temperature Difference on Residual Stress
In order to study the influence of temperature difference on residual stress, the entire simulation process is divided into multiple analysis steps. The method with multistep analysis is adopted, setting different temperature in the two zones in each step. The temperature of each zone is set using temperature field predefined function in Abaqus to introduce the temperature difference. Temperature difference for multiple analysis steps is shown in Table 2.
Temperature difference for multiple analysis steps.
The FE simulation is conducted to study the effect of temperature differences on the residual stress around addendum of planetary gear. Five nodes are chosen around addendum of planetary gear, numbered by 1006, 1005, 1004, 1003, and 1002, respectively, as is shown in Figure 3 (a). The relationship between residual stress of nodes and temperature differences are presented in Figure 3 (b), where it can be found that the residual stress value gains an obvious rise as temperature difference increases from 20°C to 120°C.

Residual stress of nodes around addendum of planetary gear. (a) The number of nodes around addendum of planetary gear (b) Relationship between residual stress and temperature differences.
The influence of temperature difference on the residual stress around pitch circle of planetary gear is analyzed by using FE simulation. Five nodes are chosen around pitch circle of planetary gear, numbered by 4650, 4575, 7697, 3916, and 4741, respectively, as is shown in Figure 4 (a). The relationship between residual stress of nodes and temperature differences are shown in Figure 4 (b), where it can be found that the residual stress value gains an obvious rise as temperature difference increases from 20°C to 120°C, showing similar variation.

Residual stress of nodes around pitch circle of planetary gear. (a) The number of nodes around pitch circle of planetary gear. (b) Relationship between residual stress and temperature differences.
The effect of temperature differences on the residual stress around dedendum of planetary gear is studied by the FE simulation. Five nodes are chosen around dedendum of planetary gear. The nodes are numbered by 978, 977, 80, 976, and 79, as is shown in Figure 5 (a). The relationship between residual stress of nodes and temperature differences are presented in Figure 5 (b). The residual stress of each node around dedendum similar to the law around pitch circle, increases as temperature difference increases from 20°C to 120°C.

Residual stress of nodes around dedendum of planetary gear. (a) The number of nodes around dedendum of planetary gear (b) Relationship between residual stress and temperature differences.
4. The FE Model of Planetary Gear Pair
4.1. The Establishment of the FE Model
The FE model of the planetary gear pair is established by Abaqus, as shown in Figure 6, which includes sun gear, planetary gear, and ring gear.

FE model of planetary gear pair.
The parameters of FE model used in the study are shown in Table 3. For each CAD model, a high quality mesh is generated by Abaqus, which fully interfaces with Proe program. An appropriate element type and a degree of refinement characterized by an adequately small element size are sought in different contact areas to ensure convergence of the results. In accordance with the principle that mesh density of sun gear should be greater than that of the planetary gear and mesh density of the planetary gear should be bigger than the ring gear, mesh generation is made on the sun gear, planetary gear, and ring gear, respectively, by means of hexahedron and tetrahedron.
Parameters of planetary gear pair.
For the current stage of the analysis, due to the lack of information concerning the precise material properties for some of the model components, all the materials were assumed to be linear, elastic, homogeneous, and isotropic. The sun gear model has a total of 181,184 elements with 207,675 nodes, the planetary gear model 126,776 elements with 126,776 nodes, and the ring gear model 455,727 elements with 115,010 nodes.
Table 4 is used to specify the inputs, boundary conditions, material properties, and load conditions. The U i (i = 1, 2, 3) represent the translational degrees of freedom in the x-axis orientation, in the y-axis orientation, and in the z-axis orientation; the UR i (i = 1, 2, 3) represent the rotational degrees of freedom in the x-axis orientation, in the y-axis orientation, and in the z-axis orientation, respectively. The model of planetary gear pair with constraints is shown in Figure 7.
Parameters of the FE model and boundary conditions.

Model of planetary gear pair with constraints.
4.2. Validation of the FE Model
During dynamic FE analysis, zero-energy modes without stiffness may appear, that is, the so-called “hourglass model” that can extend out through the grid, so that the calculated result is not significant [22]. When the pseudostrain energy accounts for about 1% of the total internal energy in Abaqus, the hourglass mode can be considered to have little effect on the calculated result. When the pseudostrain energy occupies 5% to 10% of the total internal energy, the analytical result is deemed effective. When the pseudostrain energy exceeds 10% of a total internal energy, the analytical result is considered ineffective, and appropriate measures must be taken to solve it. Through computational analysis of the model shown in Figure 2, it can be known that a pseudo strain ratio to the total internal energy is less than 6%; therefore, the analysis result was valid.
5. Influence of Residual Stress on Tooth Contact Characteristic of Planetary Gear
Contact pair consists of the master surface and slave surface in Abaqus. In the process of simulation, nodes on the slave surface could not pass through the master surface but nodes on the master surface could go across the slave surface. When both the master surface and slave surface were defined, it should be met that the selected surface with greater stiffness serves as master surface; the surface with coarser mesh should be selected to act as master surface for the contact surface with similar stiffness. The master surface of contact pairs must be continuous and the normal directions should be opposite.
The parameters of the FE model and boundary conditions are shown in Table 4. Temperature difference between body and tooth zone of planetary gear is set at 120°C and 0°C, respectively. Contact characteristics of planetary gear train at both conditions are analyzed by FE method. As is shown in Figures 8, 9, 10, and 11, no matter at meshing surface or at nonmeshing surface, equivalent stress at tooth profile of both external and internal meshing planetary gear all significantly decreases when residual stress is available, compared to when it is unavailable.

Influence of residual stress on contact stress of tooth profile at meshing surface of external-meshing planetary gear. (a) Addendum of meshing surface of planetary gear. (b) Pitch circle of meshing surface of planetary gear. (c) Dedendum of meshing surface of planetary gear.

Influence of residual stress on contact stress of tooth profile at nonmeshing of external-meshing planetary gear. (a) Addendum of nonmeshing surface of planetary gear. (b) Pitch circle of nonmeshing surface of planetary gear. (c) Dedendum of nonmeshing surface of planetary gear.

Influence of residual stress on contact stress of tooth profile at meshing of internal-meshing planetary gear. (a) Addendum of meshing surface of planetary gear. (b) Pitch circle of meshing surface of planetary gear. (c) Dedendum of meshing surface of planetary gear.

Influence of residual stress on contact stress of tooth profile at nonmeshing surface of internal-meshing planetary gear. (a) Addendum of nonmeshing surface of planetary gear. (b) Pitch circle of nonmeshing surface of planetary gear. (c) Dedendum of nonmeshing surface of planetary gear.
6. Conclusion
The residual stress field of planetary gear is simulated using the compressive stress at the high temperature zone. The method is with a series of advantages, such as universality, operability, no necessity to analyze complicated heat conduction, and ability to simulate residual stress field at each condition.
Equivalent stresses at tooth profile of both external and internal meshing planetary gear all decrease with residual stress compared to those without residual stress, no matter at meshing surface or at nonmeshing surface. For external-meshing planetary gear, residual stress influences maximum equivalent stress at addendum of meshing surface most, followed by that at pitch circle of nonmeshing and least effect is gained to that at pitch circle of meshing surface. While for internal-meshing planetary gear, the largest decline of maximum equivalent stress is obtained at addendum of nonmeshing surface, the second largest one is at addendum of meshing surface, and the smallest decline is at pitch circle of meshing surface.
It is necessary to control residual stress of planetary gear reasonably, benefiting improving dynamic characteristics and life of planetary gear train.
Footnotes
Acknowledgments
The authors gratefully acknowledge the support of the Chinese National Science Foundation (no. 51175299), the Shandong Provincial Natural Science Foundation, China (no. ZR2010EM012), the Independent Innovation Foundation of Shandong University (IIFSDU2012TS044), and the Graduate Independent Innovation Foundation of Shandong University, GIIFSDU (no. yzc10117).
