Numerical Thermal 3D Modeling of Plastic Gearing

This paper will primarily focus on the prediction of gear temperature of plastic gears using a numerical heat transfer model based on 3D Finite Difference (FDE) method. It is quite common that in most of the applications, the plastic gears are self-lubricated. Tooth surface wear is an important failure mode in plastic gears, and this primarily occurs because, at significantly higher loads, the surface temperature might increase to a value close to the melting point of the material, thereby changing the surface behavior. Thus, it is critical to compute the temperature of the gear pair in an accurate fashion. The heat source is the frictional heat dissipation due to sliding of the gears. The model is capable of solving for both helical and spur gears.

Since the gear tooth experiences a repetitive heating and cooling cycle for every rotation, the heat input is averaged over one rotation of tooth and is independent of time. The program computes both the surface and body temperature of the gears as a function of space and time. Because of the inherent nature of the implicit FDE method, there is no restriction on the discretization both in the time and space domains. This reduces the simulation time to a great extent without much compromise in accuracy. The results are correlated with experimental data, and the good agreement is achieved between the test and simulation results for different cases of load and speed. This simulation was developed for plastic gears and can be extrapolated to metal gears, with the greatest challenge being obtaining the accurate heat transfer coefficients for lubricated gears.

1: Introduction

The popularity of plastic gears is slowly rising in the gearing community, and they are constantly in the realm of replacing metallic gears in many of the low-load applications. This is primarily because of their low cost and weight, their quietness in operation, and ability to function in the absence of any lubrication. But to its downside, their performance under higher loading conditions are questionable. Adding to this, their complex thermo-mechanical and visco-elastic behavior pose a great challenge when it comes to the design of plastic gears.

Much documented research was published from the latter half of the 20th century on plastic gear design. Gauvin et al. [1] worked on the prediction of peak surface temperature of plastic gear teeth. Park et al. [2] made an experimental verification of the surface temperature and wear of plastic gears after a number of million cycles and compared them with other computer models developed for polymer gears. Tooth surface temperature of gear pairs has been experimentally verified by various authors [3-5]. Koffi et al. [6] developed a rudimentary model to predict gear temperature by considering both frictional and hysteresis heating. He concluded that the contribution of hysteresis heating in plastics is much less than that of frictional heating. Hooke et al. [7-8] developed a model to predict the surface temperature in polymer gears and its relationship to gear wear — along with comparison with experiments. Takanashi et al. [9-10] studied the heat generated in plastic gears due to friction and predicted the equilibrium temperature of a plastic gear pair.

The numerical models to predict gear temperature, specified in the above mentioned literature are mainly empirical relationships based on experimental data and are not mathematical models that involve solving physics-based heat equations. The primary disadvantage of having an empirical model is to justify its use over a wide range of operating conditions. The accuracy of the model is highly influenced by the sample of operating conditions considered in the experiments. Doll et al. [11] has done significant work on developing numerical FEA models to compute temperature of gears in contact using ANSYS, but the downside of that method is the simulation time. This paper talks about an efficient Finite Difference-based approach that provides a numerical solution to the heat equations that govern the heat transfer mechanism in a gear pair.

Previous work done by the authors [12-13] involves the development of a calculation procedure to account for transient temperature and humidity in the load distribution analysis of plastic and metal gears. This work acts as a pre-processor to that study by providing the temperature distribution of the gears in mesh.

The temperature distribution predicted by this model is compared with test results published by Hooke et al. [7-8]. The steady state results predicted by the model is in good agreement with the experimental measurements obtained using infrared camera and thermocouples. Case studies that show these comparisons are elaborated in sections below.

Methodology

This section explains the methodology involved in this model. It is important to identify the modes of heat transfer for a gear pair in mesh in order to predict the transient thermal behavior of the system. In this model, friction is assumed to be the only source of heat generation. Since most plastic gear applications are self-lubricated, the heat is lost due to convection to the atmosphere and conduction to the supporting shaft. Another assumption in this model is that a temperature distribution is imposed on the shaft that can either be time-dependent or time-independent.

2.1: Boundary Specifications

Figure 1 shows a schematic of the tooth section. Each tooth is divided into eight sections, which has its own set of governing equations to define the heat transfer mechanism. Based on Figure 1, the following sections are elaborated below.

2.2: Grid Layout

The gear tooth is discretized in space in r, Ø, and z coordinates which are radial, angular, and face-width directions, respectively. As far as this model is concerned, the domain is mapped by cylindrical grids as shown in Figure 2.

The mesh is generated using the given tooth geometry, and the nodes form the center of each finite control volume. This process is repeated in all the three directions. The mesh can be non-uniform, which means that the nodes need not be equi-spaced. Figure 3 shows an example generated mesh diagram of the tooth section. The number of nodes along the active profile, top and bottom, land and face width, are specified by the user. The involute shape of the active and inactive profile is discretized based on the number of points specified, and a straight line connection is assumed between the nodes.

2.3 Governing Equations and Discretization

As stated previously, different regions of the tooth section exhibit different boundary conditions. Each of these boundary conditions are discretized with respect to space and time, and the resulting equations are solved simultaneously to obtain the temperature distribution. As far as the discretization is concerned, implicit finite difference method [14] is used in this model. Some of the main advantages of this method are:

The method is unconditionally stable. Unlike explicit method, implicit method is stable for all values of time and space discretization steps.

This method is not restricted by the discretization lengths in both space and time domains. This increases the computational efficiency, as the time steps can be higher — to quickly attain steady state condition.

Figure 4 shows the nodal arrangement for an interior node (i, j, k). Let i, j, and k denote the nodes along the angular, radial, and face-width direction.

In the above figure, “L” represents the distance between any two nodes in a specific direction. “A” represents the cross-sectional area (2D) of the heat transfer face. i.e.:

All the interior nodes lose/gain heat only through conduction. The heat transfer mechanism can be described by the following equation:

Where, a is the thermal diffusivity of the gear material. Superscript “n” represents the temperature at n^th time instant. Equation 1 can be discretized as follows:

 

Equations 3-6 are combined to obtain Tⁿ_i,j,k as a function of the temperature in the next time step. Combining the above equations, we get:

Where

“V” is the control volume expressed in SI units. Equation 7 gives the temperature of every interior node as a function of mesh geometry and the temperature in the next time step. Similarly, the temperature of every node on each of the boundaries are obtained as a function of mesh geometry and temperature of that node in the next time step. So, if the mesh contains 500 nodes, we will be having 500 equations of the form:

Where C_i is the coefficient of neighboring nodal temperature.

For a node along active profile, the governing equation can be found using the energy balance method. Figure 5 shows the nodal arrangement of a node along the active profile. It can be seen that there is external heat that is being added to the node (i, j, k).

This external heat is because of the friction, due to sliding between the contacting surfaces. It is a function of friction coefficient, sliding velocity, and the contact pressure. The frictional heat can be computed using the relation:

It is important to note that the value of “Q” varies for every node along the profile. Since sliding is zero at pitch point, Q will be zero at a node located in the operating pitch point of the gear. Consecutively, Q is maximum either at EAP or SAP. The contact pressure distribution as a function of roll angle and face width can be obtained from WindowsLDP, which is a gear contact solver developed by The Gear and Power Transmission Laboratory at The Ohio State University [15]. Figure 6 shows a sample heat rate distribution along the active profile (no micro-geometry modification).

The heat rate constantly traverses across the teeth depending on the speed of rotation of the gear. A single tooth undergoes cyclic heating and cooling for one mesh cycle. In this model, the heat rate that is applied on the active profile is assumed not to change with time. The benefit of this approach is that the model is insensitive to the time step that is chosen for transient temperature calculation. But for this approach, the time step involved in the analysis will be in the order of the mesh cycle, and this would drastically increase the computational time. Essentially, by this idea, the boundary conditions on the active profile are made to not vary with time, whereas in reality, it does, as the gear makes one full rotation. In order to make the heat flux independent of time, a net heat rate is calculated, which is a time-based average of the heating and cooling that is seen by the nodes on the active profile. The calculation method is shown in Equation 11:

Where t₁ is the time taken for the node to heat up and t₂ is the time taken for the node to cool down. t_meshcycle is the time taken for one complete rotation of the gear. It can be deduced that t_meshcycle is essentially the sum of t₁ and t₂. The cooling rate can be calculated based on the convection relation as shown in Equation 12:

The shape of the active profile is assumed to be a perfect involute, ignoring any micro-geometry modifications that might be incorporated in the gear pair. However, this modification is ignored only during mesh generation, and its effects are considered in the temperature prediction calculations. The contact stress profile generated by WindowsLDP accounts for the micro-geometry modification, and this directly affects the temperature distribution on the surface of the gear pair.

By solving the energy balance for every node on the active profile, the governing equation can be found out and it is shown in Equation 13:

Where

k_cond being the thermal conductivity of the gear material. For plastics, the thermal conductivity is in the order of 0.5 W/mk. This low thermal conductivity causes most of the heat generated to be retained at the surface of the tooth. This is even more critical when it comes to predicting the gear failure modes like tooth wear, which is caused by the tooth surface melting at high operating temperatures.

For the nodes along the inactive profile, the methodology is essentially the same, except for the fact that the nodes are always cooled. So, the term Q_avg in Equation 13 is replaced by Q_cool from Equation 12. Figure 7 shows the nodal arrangement of the nodes in the inactive profile.

It can be seen that the heat is constantly removed from the surface through convection. Using energy balance method, the governing equation can be written as shown in Equation 15:

Where the terms R_1..5 are explained in Equation 15. For the nodes along the top land and bottom land, the governing equations resemble the ones for the nodes on inactive profile, as the heat transfer mechanism is the same in both cases — convective cooling. Figure 8 shows the nodal arrangement of the nodes along top land and root fillet boundaries.

The discretized form of the PDEs for the nodes in this boundary is given in Equation 16:

Where

It is assumed that all the nodes that are in contact with the shaft are at a specified temperature at any given instant of time. This specified temperature can either be a constant number or a time-dependent distribution. For the nodes along the shaft:

For the nodes along the front face of the tooth (“B” from section 2.1), the nodal arrangement is shown in Figure 9. The nodal arrangement for the nodes along the rear face is exactly the same but for the node (i, j, k+1). Instead, the node (i, j, k−1) will take its place.

Again, the heat transfer mechanism is exactly the same as the ones on the top and bottom land. The discretized form of the equations goes like this:

Where the terms R_1-6 has their original meaning, it is important to keep in mind that the Q_cool is computed using Equation 12. The heat transfer coefficients for each of the boundaries can either be specified by the user, or the values that are computed by the program for a given gear can be used. The methodology to compute heat transfer coefficients is elaborated in the next section.

Finally, the nodes along the boundaries C, F (from Section 2.1) exhibit symmetry in the angular direction, i.e. the temperature profile for every tooth section is identical. The discretized form of the governing PDEs are shown in Equation 20:

Where the suffix (n_ang-1) represents the node before the last node in the angular direction as shown in Figure 10.

After obtaining the equations for all the nodes at current time steps, they are solved simultaneously to compute the temperature of a given node at every time step. The connectivity matrix is generated by filling out the right rows and columns that match the spatial dependence of every node according to each of its governing equations. The connectivity matrix looks like the one shown in Equation 21:

The above equation is of the form KX = F, where X is the unknown. In the above case, Tⁿ⁺¹_i,j,k is the unknown, and the K matrix can be inversed to get the necessary solution at every time step. This process is repeated for every time step, and the transient temperature is obtained.

It is important to have knowledge on the heat transfer coefficients and the friction coefficients that are involved in the process. The next section elaborates on how the convective coefficients are obtained.

2.4 Heat Transfer Coefficient

For plastic gear applications, Takanashi et al. [9-10] formulated an empirical relationship for the heat transfer coefficients that essentially depend on the speed and the properties of the cooling medium. Since most of the plastic gearing applications are self-lubricated, the cooling medium is assumed to be air. The convective heat transfer coefficient depends on the following:

• Flow Geometry: It is the shape of the surface over which the fluid flow occurs.
• Material Property: It includes the property of the material such as surface roughness, etc., and also the property of fluid — the Prandtl number. It is defined as the ratio of momentum and thermal diffusivity of the fluid.
• Reynolds Number: The boundary layer conditions are strongly influenced by this parameter. It is important to predict whether the flow is laminar or turbulent.

It can be intuitively stated that the gears cool faster when they run at greater speeds. However, to ascertain this statement, works from the past have been studied, which predict the relation between heat transfer coefficient and the rotational speed of the cooling. Takanashi et al. [9-10] have come up with an equation that relates the heat transfer coefficient, the rotational speed and other fluid properties. The overall heat transfer coefficient is calculated as shown in Equation 22:

Where RP is the pitch radius of the gear, w is the rotational speed, v_a is the kinematic viscosity of air and a_air is the thermal diffusivity of air. m and FW represent the module and face width of the gear respectively. The next section elaborates the functionalities of the program by analyzing a gear pair and comparing it with the experimental data found in the literature search.

3: Case Study – Resultsand Discussion

Mao et al. [7–8] conducted a series of experiments to predict the gear temperature of plastic gear pairs in mesh. These results have been used as a metric for validating the temperature profiles predicted by the model presented in this paper. Some of the other capabilities of the program are also presented in this paper which include:

• Prediction of transient gear bulk temperature.
• Effect of micro-geometry modification on the surface temperature distribution of the gear.
• Prediction of micro-geometry due to temperature change based on previous work by the author [12-13] and Kashyap [16].

Table 1 summarizes the gear geometry that was used in this paper. As mentioned before, it is the same gear pair that was used by Mao et. al [7-8].

The ambient temperature conditions are summarized in Table 2.

The ambient and the shaft temperature values are taken from the experiments conducted by Mao et al. and are imposed in the model. The gears were tested at various levels of speed and torque and measured the non-contact surface temperatures (inactive flank) using infrared video cameras. Figure 11 shows the comparison between their experimental results and the temperature values generated by the program. The values shown are the maximum surface temperature of the non-contacting flank of the gear. The predicted surface temperatures of the inactive flank are in good agreement with the experimental data.

For the 500, 1,000, and 1,500 rpm speed cases, the experimental data shows a change in slope at about 7 Nm torque. This might be due to change in boundary condition with respect to ambient temperature during the testing. The airflow in the gap between the gear teeth as the gears come out of mesh will be turbulent in nature. It is fair to assume that the temperature of the air within each pocket is a constant. Theoretically, the heat transfer coefficients will vary across the gear flank as the gears act like an air pump. But, capturing this physics is beyond the scope of this work, and one might need to incorporate detailed CFD models to accurately predict the heat transfer coefficients. One other important assumption of this work is that the coefficient of friction is assumed to be 0.11 for Acetal-Acetal surface pair based on the Delrin data sheet [17]. In reality, the coefficient of friction depends on load and speed conditions, but it is assumed to be constant in this model. The model relies on the user to provide reasonable coefficient of friction values.

The steady state temperature values of the inactive profile as predicted by the model is shown in Figure 11. It can be observed that they are close to being linear with respect to torque for a given speed. This can be attributed to the fact that the change in torque level is reflected only in the contact stress distribution that is fed to the model to predict the heat flux generated due to friction. The model takes approximately six seconds to compute the temperature distribution for a single load case and is therefore a computationally efficient method. This comparison checks for the following:

• Reasonable temperature prediction within the realm of experimental data.
• Dependence of temperature on load and speed of the gear pair.

This comparison is to verify the model behavior and not to correlate with this set of experimental data because of the unknowns that exist with respect to boundary conditions, surface parameters, etc. Figure 12 shows the heat transfer coefficients that were used in this model as a function of gear speed. This model assumes the heat transfer coefficients are the same for all boundaries of the tooth section that experience convective cooling. However, the user can override the program to specify different heat transfer coefficients for different boundaries. The relationship with the speed is non-linear, as can be seen from Equation 22.

To illustrate the capabilities of the model, a speed of 1,000 rpm and a torque of 7 Nm was chosen. Figure 13 shows the steady state temperature distribution of the active flank of the gear pair.

The maximum temperature occurs at the EAP of the teeth due to the peak sliding velocity occurring at either EAP or SAP. At the SAP, heat is conducted to the shaft and convectively cooled by the ambient air. This causes the EAP to have higher temperatures than the SAP. The temperature at the edges of the teeth are approximately 10 degrees less than the temperature at the center. This can be attributed to the fact that nodes at the edges cool faster because of the convective heat loss to the atmosphere. The minimum temperature on the active profile is observed at the pitch point. This is congruent to the fact that sliding is zero at pitch point. This can be seen from the heat flux distribution at the center of the teeth, as shown in Figure 14.

In order to study the effect of micro-geometry on the surface temperature, a profile crown of 25 microns and a lead crown of 40 microns is used on the active flank to predict the contact stress distribution for input to the heat transfer model. Figure 15 shows the heat flux distribution as a result of the micro- geometry modification.

The effect of lead crown can be seen in the flux distribution. However, the effect of the profile crown is not that evident because of the multiplication of sliding velocity with contact stress. Since sliding is maximum at SAP, EAP, and zero at pitch, this type of flux distribution is seen. Figure 16 shows the temperature distribution after applying the micro-geometry modification.

It is clearly observed that temperature at the edges of face width are drastically lower than the center, and the overall temperature of the surface is lower when compared with the case of no micro-geometry. This can be attributed to the reduction of overall contact zone in the gear mesh. These results show the sensitivity of the model to micro-geometry and its effect on the temperature distribution.

Due to increase in surface temperature of plastic gears and subsequent thermal expansion, the tooth form will be modified from its nominal involute form at ambient temperature. The resulting modification to the new micro-geometry can be predicted because of the thermal expansion. This algorithm was initially structured by Kashyap [16] and was further developed as a part of the author’s master’s thesis. The new modification, along the surface normal direction, is calculated for every point on the involute profile based on its temperature. Figure 17 shows the micro-geometry modification at steady state because of thermal expansion.

The negative sign indicates the growth in material from pure involute form. However, the main assumption in this calculation is that the tooth surface is intact and did not undergo any form of wear or any other surface abnormality during the operation. It can be seen that the zone with the highest temperature has expanded at a higher rate than the other regions. Ideally, one approach is to iteratively compute contact stress distribution due to this micro-geometry and re-calculate temperatures at regular intervals. This is one of the future goals of this work.

4: Summary

The plastic gear surface temperature distribution can be predicted using the developed model. The finite difference method was used to perform thermal analysis on a gear pair for predicting surface and bulk temperatures. This is a physics-based model and does not use any empirical relationship to compute temperature.

The model predicts that the steady state temperature increases with torque, which confirms what is expected. The relationship between rotational speed and temperature is non-linear, due to the speed dependency of heat transfer coefficients. The model also predicts the effect of micro-geometry modification on temperature distribution, and vice versa. Profile or lead crown on the gear teeth can cause the maximum temperature to increase on the active flank, but at the same time, reducing the temperatures at the edges. The scope of this work can be expanded by developing simple analytical models to predict the heat transfer and friction coefficients for various speed and load conditions that can be fed as an input to this model.

References

1. Gauvin, R., Girard, P., and Yelle, H. “Prediction of the peak temperature on the surface of thermoplastic gear teeth”, AGMA, 1983.
2. Park J.Y., Plastic gear surface temperature and wear behavior. Master’s thesis, The Ohio State University, 2007.
3. Yousef, S. S., Burns, D.J., and McKinlay, W., “Techniques for assessing the Running Temperature and Fatigue Strength of Thermoplastic gears.” Mechanism and Machine Theory, 1973, pp. 175–185.
4. Hackman, H., and Strickle, E., “Nylon Gears”, Konstruktion, Vol. 3, No. 18, 1966, pp. 81–94.
5. Akozan, M., Yelle, H., and Gauvin, R., “Heat transfer between a Plastic on a Metal Gear and Ambient Air”, ASME Int. Power and Transmission gearing conf., Cambridge, MA, Oct. 1984.
6. Koffi, D., R. Gauvin, and H. Yelle, “Heat generation in thermoplastic spur gears”, Journal of Mechanisms Transmissions and Automation in Design – Transactions of the ASME 107, 1985, pp. 31–37.
7. Mao, K., Li, W., Hooke, C.J., Walton, D., “Polymer gear surface thermal wear and its performance prediction”, Tribology International 43(2010) 433–439.
8. Mao, K., Langlois, P., Hu, Z., Alharbi, K., Xu, X., Milson, M., Hooke, C.J., “The wear and thermal mechanical contact behaviour of machine cut polymer gears”
9. Takanashi, S. and A. Shoji, “Temperature Increase of Teeth of Plastic Gears.1. Heat Generated by Friction on Plastic Gear Tooth Surface”, Science Reports of The Research Institutes of Tohuku University Series a- Physics Chemistry and Metallurgy 28, no.1, 1979, pp. 93–102.
10. Takanashi, S. and A. Shoji, “Temperature Increase of Teeth of Plastic Gears.2. Equilibrium Temperature of Plastic Gear Tooth Surface”, Science Reports of The Research Institutes of Tohuku University Series a- Physics Chemistry and Metallurgy 28, no.1, 1979, pp. 103–115.
11. [Doll, P., Neil, “Modeling Thermomechanical behavior of polymer gears”. Master’s thesis, University of Wisconsin, Madison, 2015.
12. Raghuraman, N. “A Calculation Procedure to Account for Tolerances, Transient Temperature and Humidity in the Load Distribution Analysis of Plastic and Metal Gears”, Master’s thesis, The Ohio State University, 2013.
13. Houser, D., Raghuraman, N., Smith Z., “A Procedure For Evaluating The Effect Of Microgeometry Modifications On Mixed Material Gears Operating Under Extreme Conditions Of Tolerance Stack-Up And Temperature”, ANTEC-SPE Conference 2013.
14. Frank P. Incropera, David D. DeWitt, “Fundamentals of Heat and Mass Transfer”.
15. Houser, D.R., “LDP User’s Guide.” Gear and Power Transmission Research Laboratory, The Ohio State University, Columbus Ohio, 2012.
16. Kashyap, S. “Development of a procedure to describe plastic gear geometry after a temperature change with application to the prediction of gear load distribution”, Master’s thesis, The Ohio State University, 2011.
17. Product Information, DuPont Delrin 100 NC010 Acetal Resin

*Printed with permission of the copyright holder, the American Gear Manufacturers Association, 1001 N. Fairfax Street, Suite 500, Alexandria, Virginia 22314. Statements presented in this paper are those of the authors and may not represent the position or opinion of the American Gear Manufacturers Association. (AGMA) This paper was presented October 2016 at the AGMA Fall Technical Meeting in Pittsburgh, Pennsylvania. 16FTM19