A 2D finite element based approach to predict the temperature field in polymer spur gear transmissions

This article describes a new numerical approach to determine the temperature field of polymer spur gears during their operation. The approach is based on an uncoupled procedure in which a mechanical problem is solved to determine the amount of heat that is generated by friction during the meshing of the gears, and then this heat is considered as a thermal load to perform a thermal analysis of a finite element model of the transmission. The amount of heat generated by friction is determined from the results of a numerical loaded tooth contact analysis of the transmission, which is based on the finite element method. The generality of the finite element method enables this approach to be applied to any kind of spur gear transmission, regardless of the geometry and the material of the gears and lubrication conditions. The resulting approach is applied to determine the temperature field of a spur gear transmission where polymer and metallic gears are combined, under several different operating conditions. The results obtained from this approach are compared to those obtained from experimental analyses, showing a good degree of similarity between them.


Introduction
In recent times polymer gears have been replacing metallic gears in several industrial applications. According to Davis [1], this is because polymer gears present a number of advantages compared to metallic ones, especially in terms of their low cost and density. However, polymer gears still have some limitations that prevent them from being used in a wider range of applications.
As explained by Tsukamoto [2], one of these limitations is related to the fact that both their mechanical properties and their dimensional stability depend greatly on their operating temperature. Furthermore, the elevation of the temperature of the gears may not only have a negative impact during the operation of the transmission, but can also cause their premature failure [3].
Unlike metallic gears, polymer gears are usually operated without any kind of lubricant. The lack of lubricant not only accentuates the generation of frictional heat during the meshing of the gears, but also implies that less heat is dissipated from the transmission. This fact, combined with the low thermal conductivity of polymers, favors the elevation of the surface temperature of these gears during their operation. For these reasons, the determination of their operating temperature has become an important part of the design process of polymer gears, and several approaches have been developed for such a purpose.
In general, all these approaches consist in two steps. In the first step a mechanical problem is solved to determine the amount of energy that is generated during the meshing of the gears, which is mostly dissipated as heat. In the second step a heat transfer problem is solved to determine the variations in the temperature of the gears produced by the dissipated heat.
It is commonly accepted that the major sources of heat generation in polymer gear transmissions are the dissipation of frictional and hysteretic energy. However, Koffi [4] pointed out that the energy generated by friction is much greater than the energy generated by hysteresis and, in consequence, the latter is neglected in the great majority of cases. Thus, the solution to the mechanical problem is usually carried out to determine the amount of energy that is generated by friction at the contact interface.
In the classical approaches, analytical methods are used to solve both the mechanical and the thermal problems. In almost all of these methods the gears are approached to lumped systems by assuming a uniform temperature across their entire geometry, and the first law of thermodynamics is used to determine the temperature of the gears when thermal equilibrium is reached. Relevant analytical methods are those developed by Hachmann [5], Takanashi [6] and Hooke [7].
The main advantages of the analytical methods are that they are computationally fast and simple to implement. However, although these methods provide insight into the thermal behavior of the gears, their accuracy and applicability are limited by the fulfillment of the hypotheses under which they are formulated.
When more comprehensive descriptions of the temperature field across the gears are sought, it is a common practice among researchers to solve the thermal problem using numerical methods based on the discretization of the domain. For example, in the studies carried out by Taburdagitan [8], Doll [9] and Roda-Casanova [10] the temperature of the gears is determined through a transient thermo-mechanical analysis of a finite element model of the spur gear transmission.
In these transient analyses, the time domain under study is discretized into a set of time increments, and the mechanical and thermal problems are solved for each of them simultaneously. This means that for each time increment a coupled thermo-mechanical problem involving a large number of degrees of freedom and several non-linear contacts must be solved. In consequence, these transient analyses are time consuming and computationally demanding, especially when large time intervals are studied.
An interesting idea to reduce the computational cost of these numerical solutions is to uncouple the solution to the mechanical problem from the solution of the thermal problem. In addition, and taking advantage of the cyclic nature of the gear tooth contact, a further reduction of the computational cost can be achieved by solving the mechanical problem just once to determine the heat generated by the engagement of a pair of teeth, and then use that heat as an input to perform a numerical solution to the thermal problem.
Following these ideas, Patir [11] and Anifantis [12] determined the temperature of metallic gears in running conditions through a thermal analysis of a two-dimensional finite element model of the transmission. Long [13], Luo [14] and Li [15] developed three-dimensional thermal finite element models of metallic gear transmissions and conducted parametric studies of the variation of the temperature of the gears. Shi [16] used a three-dimensional thermal finite element model to analyze the bulk and flash temperature of a locomotive traction gear. This uncoupled strategy has also been applied to determine the temperature of polymer gear transmissions. Mao [17] and Raghuraman [18] determined the temperature field of polymer gears by means of the finite difference method. Evans [19] presented a simplified three-dimensional finite element model to predict the temperature of polymer gears, in which the gear assembly is approached to a rod and disc arrangement. Fernandes [20] used the finite element method to predict the flash and bulk temperature field of plastic gears running in lubricated and unlubricated conditions.
In these studies the solution to the mechanical problem and the determination of the amount of heat generated by friction are performed through a loaded tooth contact analysis (L-TCA) of the transmission. For this purpose, the Hertz contact theory [21] is used to find an analytical solution to the contact problem that arises when two gear teeth are in contact. Hertz contact theory is a fast and computationally efficient method to solve contact problems, but its applicability is limited by the fulfillment of the hypotheses under which this theory is developed. When these hypothesis are fulfilled, Hertz contact theory has proven to be an interesting method to solve gear tooth contact problems [22,23].
But, in general, the hypotheses of Hertz contact theory are not fulfilled in gear transmissions involving polymer gears, because in these cases the contact between the gear teeth can be subjected to large strains, tip contacts and non-negligible friction conditions. Thereby, erroneous results may be obtained when using the Hertz contact theory to conduct a L-TCA of a polymer gear transmission. Furthermore, polymer gears are subjected to other phenomena (such as off line of action contacts, variations in the actual location of the contact area, complex load sharing distributions, etc.) that can have an important influence over the contact solution, but cannot be easily taken into account when performing the L-TCA by means of analytical methods.
These limitations and difficulties evidence that more comprehensive methods should be used to conduct the L-TCA from which the amount of heat generated during the meshing of polymer gears is determined. The finite element method could be a good candidate for this purpose, as this method is capable to overcome the limitations associated to analytical methods, allowing us to study the interaction between two bodies regardless of their geometry, materials or contact conditions. Its adequacy to conduct L-TCA of polymer gear transmissions has been proven by Langlois [24] and Karimpour [25]. Furthermore, the generality of the finite element method makes it suitable to study the meshing of the gears regardless of their geometry, the material in which they are manufactured or the lubrication conditions under which the transmission is operated.
Thus, the objective of this work is to propose a new uncoupled numerical approach to determine the temperature of polymer spur gear transmissions in operating conditions. In this approach the amount of heat generated during the engagement of the gears will be determined from the results of the L-TCA of the transmission, which will be carried out by means of the finite element method. Finally, the calculated heat will be considered as a thermal load to conduct a thermal analysis of a finite element model of the transmission, from where the temperature field of the gears in operating conditions will be determined.

Overview of the proposed method for the determination of the working temperature of spur gears
In a spur gear transmission, the teeth of the drive gear enter in contact with the teeth of the driven gear in a combination of rolling and sliding motions. As a consequence of this sliding motion, frictional energy is generated that is mostly dissipated as heat, which is released over the contact surfaces of the gears and produces an increase in their working temperature.
This article describes a numerical approach to determine the temperature field of spur gears during their operation, based on a thermal analysis of a two-dimensional finite element model of the gear of interest. By using a two-dimensional approximation to the thermal problem it is assumed that neither heat flow nor temperature change is produced in the width direction of the gears. However, it is known that in reality there exists a variation in the temperature along the width direction of the gear, which is caused by the convective heat transfer across its end faces.
Some authors [11,15,26] have noted that this temperature variation along the width direction of the gear can be approached to a parabolic function, which can be characterized by the temperature at the mid-section of the gear (that defines the maximum value of the parabola) and a parabola coefficient (that defines the shape of the parabola). For gears having moderate to large face widths, Patir [11] showed that the convective heat transfer across their end faces does not have a severe impact on the temperature level at their mid-section, although it may have a certain influence over the shape of the parabola. In these cases, a two-dimensional thermal analysis can successfully predict the temperature field at the mid-section of the gear. In contrast, a loss of accuracy can be expected when applying this two-dimensional approximation to narrow gears, and this is one limitation of this approach.
In this thermal analysis, the heat supplied to each tooth of the gear by the dissipation of frictional energy is simulated through a moving heat source that acts over its contact profiles. For the definition of the moving heat source, the contact profile of the gear teeth is parametrized through an intrinsic coordinate ξ, as shown in Fig. 1a, in such a way that any point on the contact profile has an associated value of ξ. This intrinsic coordinate has length units, and it is defined for an interval ξ = [ξ 0 , ξ f ], where ξ 0 represents the start point of the contact profile and ξ f represents its end point.
Taking this parametrization of the contact profile into account, the moving heat source is represented by a heat flux function q(ξ, t) that describes the amount of energy that is released as heat at each point of the contact profile of a gear tooth (defined by the intrinsic coordinate ξ) in a given instant of time (defined by variable t). Because of the cyclic nature of the gear tooth contact, the heat flux function is a periodic function that is repeated for each gear cycle. A gear cycle is defined as a complete revolution of the gear, and its duration is denoted by ∆t cycle . Each gear cycle starts at time t 0 and ends at time t f . Figure 1b shows a sample heat flux function during a gear cycle. The thermal power Q(t) that is delivered to the contact profile at a given instant of time by the moving heat source can be calculated using Eq. 1a. In a similar way, the thermal energy intensity E(ξ) that is released over each point of the contact profile by the moving heat source during a gear cycle can be calculated using Eq. 1b.
The determination of the heat flux function is the first step of the proposed procedure, and is discussed in sections 3 and 4 of this work. Once the heat flux function has been determined, it can be considered as a thermal load to perform a thermal analysis of the gear using the finite element method. As described in section 5, the thermal analysis can be performed under either steady-state or transient conditions, depending on whether thermal equilibrium is assumed to be reached or not. As a result of the thermal analysis, the temperature field across the gear geometry is obtained.

Theoretical background for the determination of the heat flux function
In order to establish a theoretical background for the determination of the heat flux function, let us consider a spur gear transmission such as the one shown in Fig. 2, composed of a drive gear and a driven gear. The drive gear is rotating around point O 1 with a constant angular velocity ω 1 and, consequently, the driven gear is rotating around point O 2 with an angular velocity ω 2 . The position of point O 1 is coincident with the origin of the global coordinate frame. A driving torque C 1 is applied over the drive gear, and a resistive torque C 2 is applied over the driven gear.
The heat flux function is determined by taking one of the gears of the transmission as a reference, which is defined as the master gear, whereas the other gear is defined as the slave gear. The contact profile of the master gear is denoted by Σ m , and the contact profile of the slave gear is denoted by Σ s . Without loss of generality and for illustrative purposes, the driven gear is selected as the master gear in Fig. 2. The contact profile of the master gear is parametrized through the intrinsic coordinate ξ m , as illustrated in Fig. 1a. Considering this coordinate ξ m , the position of any point M ∈ Σ m at a given instant of time (defined by the variable t) is defined by the vector r m (ξ m , t). In a similar way, the contact profile of the slave gear is parametrized through the intrinsic coordinate ξ s , so the position of any point S ∈ Σ s at instant of time t is defined by the vector r s (ξ s , t).
Taking advantage of the parametrization of the contact profile of the master gear, the normal contact pressure distribution that arises when the gear teeth are in contact can be described by the function p(ξ m , t). Additionally, if the effect of friction is taken into account, a tangential contact stress distribution arises that can be described by the function τ (ξ m , t).
Consider a pointŜ in the contact profile of the slave gear, which is defined as the closest point in Σ s to point M . The position of pointŜ is defined by the position vector r s (ξ s , t), whereξ s is the value of the profile intrinsic coordinate ξ s that satisfies the minimum distance condition: Assume now that, at a given instant of time, point M lies within the contact area. In that case point M and pointŜ are coincident, and the sliding velocity v S (ξ m , t) of the contact profile Σ s with respect to point M can be calculated as: where v m (ξ m , t) and v s (ξ s , t) are the velocities of points M andŜ, respectively. The modulus of the sliding velocity defines the sliding speed v S (ξ m , t), and its direction vector defines the sliding direction t S (ξ m , t).
As the frictional energy generated at the contact interface depends on the tangential contact stress and the sliding speed, and only a portion of this energy is released as heat at the contact surfaces of the master gear, the heat flux function can be defined as: where λ represents the fraction of the dissipated energy that is converted to heat and ϕ(ξ m , t) is a heat partition factor for the distribution of this heat between the interacting surfaces. Assuming the gear tooth contact to be a particular type of band contact, heat can be partitioned following Blok's equation [27]: where e m and e s are the thermal effusivities of the master and slave gears, respectively. On the other hand, v m T (ξ m , t) and v s T (ξ s , t) are defined as:

Determination of the heat flux function through a numerical loaded tooth contact analysis of the transmission
As has been shown in the previous section, the determination of the heat flux function (according to Eq. 4) requires the calculation of the tangential contact stress distribution that arises when two gear teeth are in contact, the sliding speed at the contact interface, and the heat partition factor. The determination of these parameters requires the solution of the contact problem that arises during the engagement of a pair of gear teeth.
Several approaches have been developed to solve the contact problem in gear drives, under the name of loaded tooth contact analysis (L-TCA). Some of them are analytical methods that are based on Hertz contact theory [21], whereas others are numerical methods such as the finite element method.
Many authors have relied on analytical methods to solve the contact problem when determining the heat flux function [13,14,17,20]. For a given contact position, these methods determine the theoretical contact point through an unloaded tooth contact analysis (U-TCA) of the transmission. Then, Hertz contact theory [21] is applied to determine the contact pressure distribution that is formed around the theoretical contact point, and the tangential contact stress distribution is determined from the contact pressure distribution following Coulomb's law. Finally, the sliding velocity and the heat partition factors are determined using rigid body kinematics, and the heat flux function is then determined following Eq. 4.
The main advantage of the analytical methods is that they are able to provide a fast and computationally efficient solution to the contact problems. However, their applicability is limited by the hypotheses of the underlying theory. Although there are some cases of gear tooth contact that can be analyzed under these hypotheses, in general they are not fulfilled in the case of polymer gears because: (i) Hertz contact theory was developed to solve frictionless contact problems, but polymer gears are usually operated without any kind of lubricant, and in consequence the effect of friction can be significant. (ii) The reduced stiffness of polymer gears (compared to the stiffness of metallic gears) can imply large deformations that are outside the scope of the linear small strain theory of elasticity under which the Hertz contact theory is formulated. (iii) Hertz contact theory can only be applied when the contact surfaces are continuous and non-conforming.
However, the contact between gear teeth can be produced in areas of the contact profile where these requirements are not fulfilled and, consequently, the contact problem cannot be solved by means of Hertz contact theory. This issue is especially relevant when the contact is produced in the vicinity of the tip of the tooth.
Furthermore, there are other issues that can make it difficult to apply analytical methods to solve contact problems involving polymer gears. One of these difficulties consists in the determination of the acual location of the contact area. On the one hand, the large deflections of the polymer gears under load cause the theoretical contact area determined from the U-TCA to differ from the actual contact area. On the other hand, the low stiffness of the polymer gears causes the gears to mesh prematurely and go out of mesh beyond the theoretical point of last contact [4,24], and this phenomenon cannot be easily taken into account when performing a U-TCA.
Moreover, solving the tooth contact problem through Hertz contact theory requires an estimation of the load that is borne in each of the simultaneous contacts that are produced between the gear teeth. As explained by Walton [28], in metallic gears the theoretical load sharing distribution does not differ substantially from the actual load sharing distribution, and in consequence the former can be successfully used to perform the loaded tooth contact analysis of the transmission. However, in polymer gears the differences between the theoretical and the real load sharing distribution might be significant, leading to erroneous solutions to the contact problem.
All these difficulties lead us to think that more comprehensive methods should be used to solve the contact problem in polymer gears in order to determine the heat flux function. For these reasons, in this work the solution to the contact problem and the determination of all the parameters that are involved in the calculation of the heat flux function are performed through a set of static analyses of a two-dimensional finite element model of the transmission, like the one shown in Fig. 3. As illustrated by Langlois [24] and Karimpour [25], all the limitations and difficulties discussed above are overcome by using the finite element method to perform the loaded tooth contact analysis, although the computational cost of the approach may be globally increased. For the definition of this finite element model, the geometries of both the drive and the driven gears are generated and then discretized into linear quadrilateral finite elements using the method described by Argyris [29]. Each gear geometry in the finite element model is composed of at least three teeth, which makes it possible to consider a realistic load sharing distribution between teeth when simultaneous contacts are produced.
Rigid edges are defined in the bore of the drive and driven gears, where their respective shafts should be placed, as shown in Fig. 3. The movements of these rigid edges are coupled to the movement of reference nodes that are coincident with the center of rotation of the gears.
The analysis is performed under plane strain hypotheses. There are several material models that can be considered for the gears, depending on the conditions under which they are studied. If the deformation of the gears is assumed to be so small that they can be studied under the assumptions of the small strain theory, a linear elastic material model can be specified for the gears. If large strains are expected, hyperelastic material models should be adopted in the way illustrated in Ref. [10,30].
The boundary conditions of the finite element model are specified at the gear reference nodes. These reference nodes have all their translational degrees of freedom restricted. A constant torque C 2 is applied over the reference node of the driven gear, and the rotation of the system is restricted by prescribing a fixed value θ i for the rotational degree of freedom of the drive gear reference node.
Surface-to-surface finite sliding contact pairs are specified between the teeth of both gears. A penalty method is established to take into account the tangential behavior of the contact, as recommended by Trobentar [30]. In the simplest of cases, a constant value is defined for the coefficient of friction. However, provided that enough information is available, the coefficient of friction can also be defined as a function of the slip ratio, the contact pressure and the operating temperature.
The static analysis is repeated for several contact positions distributed along the engagement of a pair of teeth. Each contact position i is characterized by a value of the angle θ i (Fig. 3), which varies between θ ini (which defines the initial contact position) and θ end (which defines the final contact position). To ensure that the whole engagement of the pair of teeth is analyzed, the initial contact position is defined as the position in which the contact between the previous pair of teeth is produced at the pitch circle, and the final contact position is defined as the position in which the contact between the following pair of teeth is produced at the pitch circle.
The analyzed contact positions form an increasing sequence in which θ i < θ i+1 . Each analyzed contact position has an associated time increment t i . Assuming that the transmission under study is operating at a constant angular velocity ω 1 , the time increment t i associated with contact position i may be calculated as: where t ini is the time instant associated to the initial contact position, which by default is set to t ini = 0 s. After completing the set of static analyses of the described finite element model, one of the gears is selected as the master gear, and then the following results are retrieved from the analyses: • A tangential contact stress value τ j,i for each node j of the contact profile of the master gear at contact position i.
• A vector r m j,i that describes the deformed position of each node j of the contact profile of the master gear for each contact position i.
• A vector r s j,i that describes the deformed position of each node j of the contact profile of the slave gear for each contact position i.
The velocity v α j,i of any node j in the contact profile α (α = m for the contact profile of the master gear and α = s for the contact profile of the slave gear) at contact position i can be determined as the time derivative of the position vector r α j,i , which can be calculated using the central difference approximation: The deformed position r α (ξ α , t i ) and the velocity v α (ξ α , t i ) of any point along the contact profiles with associate intrinsic coordinate ξ α can be determined as a linear interpolation of the deformed position and the velocity of the surrounding nodes: where ξ α j and ξ α j+1 are the profile intrinsic coordinates of nodes j and j + 1, respectively. It is important to note that equations 9a and 9b are only valid as long as ξ α j ≤ ξ α ≤ ξ α j+1 . From these results it is possible to determine a heat flux value q j,i for each node j of the contact profile of the master gear (with associated contact profile intrinsic coordinate ξ m j ) and contact position i (with associated time increment t i ). For this purpose, Eq. 4 is rewritten as: where τ (ξ m j , t i ) is the tangential contact stress, v S (ξ m j , t i ) is the sliding speed, and ϕ(ξ m j , t i ) is the heat partition factor at node j of the contact profile of the master gear and contact position i.
As has been mentioned earlier, the tangential contact stress at node j and contact position i is obtained directly from the static analyses of the finite element model, in such a way that τ (ξ m j , t i ) = τ j,i . Note that in those cases where τ j,i = 0, the heat flux q j,i is null and no further calculations will be required. In other cases where τ j,i > 0, the sliding speed v S (ξ m j , t i ) and the heat partition factor ϕ(ξ m j , t i ) must be calculated in order to determine q j,i .
Before the sliding speed can be calculated, it is necessary to find the pointŜ ∈ Σ s that satisfies the minimum distance condition (given by Eq. 2) with respect to each node j at contact position i. Under a discretized domain, the position of pointŜ can be found as the orthogonal projection of node j over Σ s , as illustrated in Fig. 4a. The magnitude of the intrinsic coordinateξ s that defines the position ofŜ along Σ s can be retrieved using Eq. 9a. The sliding speed v S (ξ m j , t i ) is determined as the modulus of the sliding velocity v S (ξ m j , t i ), which is calculated using Eq. 3. For this purpose, the velocity of node j of the contact profile of the master gear is calculated using Eq. 8, and the velocity of the corresponding pointŜ is calculated using Eq. 9b. Finally, the heat partition coefficient ϕ(ξ m j , t i ) is determined using Eq. 5. When the heat flux q j,i is determined for all the nodes in the contact profile of the master gear and all the analyzed contact positions, the heat flux function is interpolated from the collection of values of q j,i using bilinear interpolation, as shown in Fig. 4b.

Determination of the temperature of the gears through a thermal finite element analysis
The temperature field of the gear studied in operating conditions is determined through the thermal analysis of a two-dimensional finite element model of the gear. As explained in section 2, the thermal analysis can be performed in either transient or steady-state conditions. When the analysis is performed in steady-state conditions, it is assumed that the gear has reached its thermal equilibrium, and in consequence the temperature field across the gear geometry remains constant in time. Hence, the time dependency of the problem can be removed, and the temperature field across the gear is obtained by solving the steady-state heat conduction equation.
When the analysis is performed in transient conditions it is assumed that the gear has not reached its thermal equilibrium, and in consequence its temperature field across the gear geometry can vary with time. The time domain of the problem must be discretized into a set of time increments, and for each of them the temperature field across the gear is obtained by solving the transient heat conduction equation. For the first time increment, an initial temperature must be defined for each node of the finite element model, which enables the transient analysis to be performed from two different points of view: (i) A uniform arbitrary temperature can be defined for all the nodes of the finite element model to observe how the gear is heated during a certain period of time. (ii) The temperature field obtained from the steady-state thermal analysis can be defined as an initial temperature to observe the cyclical variations in the temperature field once the thermal equilibrium is reached.
For the construction of the finite element model, the geometry of the gear is generated and then discretized into linear quadrilateral finite elements using the method described by Argyris [29]. In this case, only one tooth needs to be considered in the finite element model, and it is not mandatory that the finite element mesh of this tooth coincides with the mesh used for the loaded tooth contact analysis. The rest of the gear geometry is simulated through a cyclic symmetry boundary condition that is established at the sides of the tooth, as shown in Fig. 5a. Of the three heat transfer mechanisms, only conduction and convection are considered in the proposed finite element model. Radiation is not taken into account because it is assumed that the amount of heat dissipated by radiation is small compared to the heat dissipated by the other heat transfer mechanisms.
Heat conduction is characterized by the thermal conductivity, the specific heat and the density of the gear material. Heat convection is modeled through three surface film conditions (denoted by Γ A , Γ B and Γ C ) that are applied to the external edges of the gear, as shown in Fig. 5a. Surface film condition Γ A is specified for the leading side of the gear tooth, Γ B is specified for the top land and Γ C is specified for the trailing side of the tooth. Each surface film condition is characterized by the environment temperature T env and a heat transfer coefficient (h A , h B and h C , respectively). Several methods for the calculation of the heat transfer coefficients can be found in the literature, and choosing one or another depends mostly on whether the transmission is operating in lubricated [13,14,20] or non-lubricated [6,9,18] conditions. As explained in section 3, the heat supplied to the gear as a consequence of frictional dissipation is simulated through a moving heat source that acts over its contact profile, and is represented by a heat flux function (the calculation of which has been described in section 4). Under a discretized domain, this  heat flux function must be converted into a set of nodal heat fluxes Q j , which are applied over each node j of the contact profile of the gear, as illustrated in Fig. 5b. In the transient analysis, the nodal fluxes are time-dependent and can be calculated using the following equation: where ξ j is the intrinsic coordinate of the node under consideration, ξ j+1 is the intrinsic coordinate of the following node, and ξ j−1 is the intrinsic coordinate of the previous node. On the other hand, b represents the face width of the gears. In contrast, in the steady-state analysis, all the variables that define the finite element model must be time-independent, including the nodal fluxes. For this reason, the heat flux function that is used in Eq. 11 to determine the nodal fluxes is replaced by an averaged heat flux function q(ξ), which, being constant in time, releases at each point of the contact profile the same thermal energy intensity as the original heat flux function. This average heat flux function is calculated as: After the analysis of this finite element model has been performed, the temperature field across the geometry of the gear tooth is obtained.

Numerical examples and discussion of results
The approach developed in this study has been applied to determine the working temperature of a spur gear transmission that is defined by the set of parameters shown in Tab. 1. Both drive and driven gears have standard involute profiles, and it is assumed that the transmission is free from manufacturing and assembly errors. The drive gear is made of AISI 1040 steel, whereas the driven gear can be made of either high-density polyethylene (HDPE) or polyoxymethylene (POM). The mechanical and thermal properties of these materials are shown in Tab. 2.
The transmission has been analyzed under several operating conditions, which are defined by a combination of input velocity ω 1 , output torque C 2 , and environment temperature T env . The operating conditions The combination of the driven gear material and one of the selected operating conditions constitutes a case of study. Considering the materials selected for the driven gear (HDPE or POM) and the operating conditions shown in Tab. 2, 16 different cases of study have been analyzed in this work. For each of them, the temperature field of the driven gear has been determined using steady-state and transient thermal analyses.
The temperature results are compared to the temperature results obtained experimentally by Singh [31] for the same transmission and operating conditions. As a result of his experiments, Singh provides the range of temperatures measured for each case of study once thermal equilibrium is reached, characterized by the maximum and minimum temperature.
In the following sections, the results obtained from the loaded tooth contact analysis, the steady-state thermal analysis, and the transient thermal analysis of the selected cases of study are discussed.

Results obtained from the loaded tooth contact analysis
The procedure described in section 4 has been applied to determine the heat flux function for each of the described cases of study. For this purpose, a finite element model of the transmission has been developed and then analyzed in 150 different contact positions that are uniformly distributed between the initial and the final contact positions. Each gear geometry has been discretized into 2208 elements and 2445 nodes, with 45 nodes in each of their contact profiles. As the temperature field will be determined for the driven gear, this gear is selected as the master gear in these numerical examples.
A constant coefficient of friction has been considered for the analysis of the whole mesh cycle. In reality, the magnitude of the coefficient of friction varies depending on the contact conditions (that include temperature, contact pressure and relative velocity between the contacting bodies, among others). However, it is a common practice among the researchers [15,17,18,19] to use an averaged coefficient of friction for the whole mesh cycle. The values considered in this work for the averaged coefficient of friction in POM/steel and HDPE/steel interactions are shown in Tab. 2, and they have been selected according to the values proposed by Singh [31]. Figure 6 summarizes the results of the loaded tooth contact analysis for a representative case of study. One of the first results that is obtained from this analysis is the tangential contact stress τ j,i associated to each node j of the contact profile of the master gear and analyzed contact position i. In Fig. 6a, each black dot indicates the position of the nodes where contact is detected (τ j,i > 0), and the collection of black dots constitutes the contact path along the contact profile.
In addition, the theoretical contact path has been calculated through an unloaded tooth contact analysis (U-TCA) of the transmission [32]. Theoretically, the unloaded contact between the profiles of spur gears is produced at a single point, whose evolution over time traces the red line shown in Fig. 6a. The start and the end of the theoretical contact path are marked with dashed vertical lines.
On the other hand, Fig. 6b shows a contour plot representing the calculated heat flux function for the same case of study. According to previous works [14,20], the heat flux function is calculated considering that 95% of the dissipated frictional energy is converted into heat (represented by parameter λ in Eq. 10). The heat flux function reaches its maximum value at the tip of the tooth (ξ = 3.473 mm), and its minimum value in the vicinity of the pitch point (ξ = 1.238 mm), where a pure rolling contact is produced and no frictional energy is generated.
The results of the U-TCA displayed in Fig. 6a show that the engagement of the teeth starts at t = 0.55 ms and ends at t = 4.44 ms, having a duration of 3.89 ms. In contrast, according to the results of the L-TCA the engagement of the teeth starts at t = 0.23 ms and ends at t = 4.83 ms. This means that the duration of the engagement predicted by the U-TCA is 15% lower than the duration of the engagement predicted by the L-TCA. Moreover, it can be appreciated in Fig. 6b that the parts of the engagement that are not considered by the U-TCA coincide with those where the heat flux function reaches its maximum values.
Moreover, it is a common assumption in analytical methods used to determine the heat flux function that the contact area is formed symmetrically around the theoretical contact point predicted by the U-TCA. However, it can be observed in Fig. 6a that the contact area predicted by the L-TCA is not symmetrical with respect to the theoretical contact point predicted by the U-TCA.
From the heat flux functions obtained, the thermal power and the thermal energy intensity functions have been calculated for each case of study considering Eq. 1a and Eq. 1b, respectively. Figure 7a shows the evolution of the thermal power delivered to the driven gear for some representative cases of study. It can be observed that for those cases of study having the same output torque, the thermal power delivered to the driven gear is greater when this gear is manufactured from POM than when it is made of HDPE. This is mainly because the coefficient of friction is larger for POM-steel contacts than for HDPE-steel contacts, and consequently larger tangential contact stresses are produced under the same applied torque, which increase the dissipated frictional energy. Furthermore, the difference in the modulus of elasticity of the driven gear also has an impact on the instantaneous power functions. On the one hand, the lower modulus of elasticity of HDPE favors a greater tooth deflection that tends to increase the duration of the engagement of the teeth, thereby expanding the time in which thermal powered is delivered. On the other hand, a smoother load sharing distribution between teeth is achieved as the rigidity of the gears is reduced. This explains why some peaks that are observed in the thermal power functions corresponding to the POM gears are not observed in the functions corresponding to the HDPE gears. These effects, which were also observed by Diez-Ibarbia [33] for metallic gear transmissions, evidence the importance of taking into account the deflection of the gears when determining the heat generated by friction. And this is especially true for polymer gear transmissions, because polymer gears may be subjected to large deflections as a consequence of their reduced rigidity.
Finally, Fig. 7b shows the thermal energy intensity distribution along the contact profile for some representative cases of study. In this figure it can be observed that the highest peaks of thermal energy intensity are produced at the tip of the tooth (ξ = 3.473 mm) and in the vicinity of the tooth root (ξ 0.234 mm). It is interesting to note that these peaks would not be reflected if the loaded tooth contact analysis were carried out using analytical methods: (i) On the one hand, as explained in section 4, Hertz contact theory cannot be applied when the tooth contact is produced at the tip of the tooth, because its hypotheses are not fulfilled. (ii) On the other hand, Fig. 6b shows that the high values of the heat flux function that produce the peak of thermal energy intensity in the vicinity of the root of the tooth are left out of the analytical study conducted using the U-TCA results.
The differences between the results obtained from both methods evidence some of the limitations of the analytical methods based on a U-TCA to determine the heat flux function in polymer gears.

Results from the steady-state thermal analysis
The heat flux functions obtained from the loaded tooth contact analysis have been used to conduct a steady-state thermal analysis of the driven gear. For this purpose, a finite element model of the gear has been developed and analyzed following the ideas outlined in section 5. In this case, the gear tooth has been discretized into 736 elements and 822 nodes, with 45 nodes in the contact profile.
Following the ideas presented in Ref. [18], the heat transfer coefficients that define the surface film conditions Γ A , Γ B and Γ C (Fig. 5a) have been calculated using Eq. 13. This empirical equation was proposed by Takanashi [6] and relates the heat transfer coefficient h to: • the gear geometry, characterized by the normal modulus (m n ) and the face width (b) of the gear, • the thermal properties of the surrounding air, which are defined by its thermal conductivity (k air ), kinematic viscosity (ν air ) and thermal diffusivity (α air ), and • the linear velocity of the gear at the pitch point (v).
As a result of these analyses, the temperature field of the driven gear when thermal equilibrium is reached is obtained for each case of study. Figure 8 shows the temperature field determined for some representative cases of study. The points of the gears where the maximum temperature (T max ) and minimum temperature (T min ) are found are marked in this figure. In addition, Fig. 8 shows a reference point within the gear body that has been selected to measure the average bulk temperature (T bulk ). It can be observed that in all the cases of study the points of maximum and minimum temperature correspond to the tip and the root of the gear teeth, respectively. Other researchers [15,16,34] have also found that these are the locations where the maximum and minimum temperatures usually take place in standard spur gears. The maximum, minimum and average bulk temperature results obtained from the steady-state thermal analysis for all the cases of study are represented in Fig. 9. In these figures, the temperature ranges obtained by Singh [31] from the experimental analyses are also displayed. In general, a good degree of similarity can be observed between the experimental results and the results obtained using the proposed method. This similarity tends to be better when the input speed is ω 1 = 1200 rpm than when the input speed is ω 1 = 600 rpm. The maximum deviation observed between the average bulk temperature and the temperature interval predicted by Singh is 3K, and it is produced at case of study POM-OC4. Figure 9 shows that the obtained maximum, minimum and average bulk temperature values increase linearly with the torque applied. This linear relation between the applied torque and the temperature of the gears has been previously reported by Hooke [7], Thyla [26] and Singh [31], among others. In addition, it can be observed that under the same operating conditions, the driven gears made of POM reach higher temperatures than the gears made of HDPE. This fact is not only favored by the larger values of the thermal energy intensity that is released over the POM gears (as observed in Fig. 7b), but also by the lower thermal conductivity of POM compared to the thermal conductivity of HDPE.

Results from the transient thermal analysis
Finally, the heat flux functions obtained from the loaded tooth contact analysis have also been used to conduct a transient thermal analysis of the driven gear. For this purpose, a finite element model of the driven gear has been developed following the guidelines presented in section 5. As explained in section 5, the transient thermal analysis can be performed from two points of view, depending on the temperature that is specified at each node at the beginning of the analysis.
First, a constant temperature has been defined for all the nodes of the finite element model, equal to the environment temperature T env . Then, a transient analysis of the finite element model has been conducted during 2000 gear cycles, by dividing each gear cycle into 200 time increments. Figure 10 shows the evolution of the temperature field of the driven gear at four different instants of time for a representative case of study. It can be observed that the evolution of the temperature field as the number of gear cycles is increased follows a similar pattern as the one obtained by Li [15]. Figure 10 shows a progressive increase of the gear temperature as the number of cycles is increased. Theoretically, this increase of the gear temperature continues until the thermal equilibrium is achieved. If the transient thermal analysis is conducted for a large number of gear cycles, the thermal equilibrium will be finally reached and the obtained temperature results will converge towards those obtained from a steady state thermal analysis. However, the number of cycles required to achieve the thermal equilibrium is usually large, and conducting a transient analysis for such a number of cycles would imply an elevated computational cost. For example, Singh [31] determined that 3 · 10 4 cycles are required to achieve the thermal equilibrium in case of study POM-OC8. For this reason it is recommended to use a steady-state thermal analysis (instead of a transient analysis) to determine the gear temperature at thermal equilibrium.
Second, the temperature field obtained from a steady-state thermal analysis has been used as a predefined field to conduct the transient thermal analysis of the finite element model. Figure 11 shows the evolution of the temperature at the maximum and minimum temperature points for two representative cases of study during two gear cycles. As these transient analyses are performed assuming that the gears have reached their thermal equilibrium, cyclical temperature variations are obtained, which are repeated for each gear cycle. The peaks that are observed in this figure represent the flash temperature, caused by the instantaneous rise in temperature produced when the contact between the gear teeth is produced. It can be appreciated that this rise in temperature is only observed at the maximum temperature point (and in general, at the nodes in the contact profile), indicating the localized effect of the flash temperature.
The temperature range obtained experimentally by Singh [31] for these cases of study is also displayed in Fig. 11. It can be observed that in both cases there is a good level of correlation between the experimental results and the results obtained from the presented approach.

Conclusions
In this work, a new numerical approach has been developed to determine the temperature field of spur polymer gears during their operation. This approach is based on a thermal analysis of a finite element model of a gear tooth, in which the frictional heat supplied to the gear is simulated through a moving heat source that is described by a heat flux function. The thermal analyses can be performed either in steady-state or transient conditions, depending on whether it is assumed that the gears have reached their thermal equilibrium or not.
Unlike in other approaches, the determination of the heat flux function is performed through a L-TCA of the transmission, which is carried out by a set of static finite element analyses. The generality of the finite element method enables the developed approach to be applied to any kind of spur gear transmission, independently of the geometry of the gears, the material in which they are manufactured or the lubrication conditions under which the transmission is operated. The accuracy of the results will be conditioned by the accuracy in which the parameters that define the response of the finite element model of the transmission are determined.
The developed approach has been applied to determine the temperature field of a spur gear transmission where polymer and metallic gears are combined, under several different operating conditions. The results obtained from these analyses enable us to draw the following conclusions: (i) The results of the numerical L-TCA reveal that the large deflections that affect polymer gears produce variations in the location of the contact area, change the load sharing distribution between teeth and increase the duration of their engagement. These effects have impact over the amount of heat that is generated during the meshing of the gears, and thereby over their temperature during the operation. However, these effects cannot be considered when the frictional heat is calculated by means of analytical methods based on U-TCA. (ii) The developed approach is suitable to determine the temperature field of the gears after their thermal equilibrium is reached. The results obtained from the application of this approach have been compared to experimental results, and a maximum deviation of 3K is observed in the calculated average bulk temperature. (iii) The developed approach allows to determine the temperature field of the gears in unsteady thermal conditions. The analyses can be conducted either to study the temperature rise of the gears when the operation is started, or to observe the cyclical temperature variations that are produced once the thermal equilibrium is reached. The results obtained under these conditions demonstrate a good level of correlation with experimental results.
This approach can help polymer gear designers to study how the material selected for the gears affects to the thermal behavior of the transmission, to propose topological modifications of the gears in order to minimize the generation of frictional heat and power loss, or to estimate the temperature of the gears under certain lubrication conditions. Further work on this topic could be directed towards adapting the proposed approach for other types of gears that cannot be analyzed by means of a two-dimensional model. Furthermore, including other phenomena (such as hysteretic heating or thermal radiation) in the finite element models could be an interesting extension of this work.