Advanced model for the calculation of meshing forces in spur gear planetary transmissions

This paper presents a planar spur gear planetary transmission model, describing in great detail aspects such as the geometric definition of geometric overlaps and the contact forces calculation, thus facilitating the reproducibility of results by fellow researchers. The planetary model is based on a mesh model already used by the authors in the study of external gear ordinary transmissions. The model has been improved and extended to allow for the internal meshing simulation, taking into consideration three possible contact scenarios: involute–involute contact, and two types of involute-tip rounding arc contact. The 6 degrees of freedom system solved for a single couple of gears has been expanded to 6 + 3n degrees of freedom for a planetary transmission with n planets. Furthermore, the coupling of deformations through the gear bodies’ flexibility has been also implemented and assessed. A step-by-step integration of the planetary is presented, using two typical configurations, demonstrating the model capability for transmission simulation of a planetary with distinct pressure angles on each mesh. The model is also put to the test with the simulation of the transmission error of a real transmission system, including the effect of different levels of external torque. The model is assessed by means of quasi-static analyses, and the meshing stiffness values are compared with those provided by the literature.

Abstract This paper presents a planar spur gear planetary transmission model, describing in great detail aspects such as the geometric definition of geometric overlaps and the contact forces calculation, thus facilitating the reproducibility of results by fellow researchers. The planetary model is based on a mesh model already used by the authors in the study of external gear ordinary transmissions. The model has been improved and extended to allow for the internal meshing simulation, taking into consideration three possible contact scenarios: involute-involute contact, and two types of involute-tip rounding arc contact. The 6 degrees of freedom system solved for a single couple of gears has been expanded to 6 ? 3n degrees of freedom for a planetary transmission with n planets. Furthermore, the coupling of deformations through the gear bodies' flexibility has been also implemented and assessed. A step-by-step integration of the planetary is presented, using two typical configurations, demonstrating the model capability for transmission simulation of a planetary with distinct pressure angles on each mesh. The model is also put to the test with the simulation of the transmission error of a real transmission system, including the effect of different levels of external torque. The model is assessed by means of quasi-static analyses, and the meshing stiffness values are compared with those provided by the literature.

Keywords Planetary transmission Á Meshing stiffness Á Transmission error Á Modeling
Nomenclature (x i , y i ) Center position of the i wheel C r Tip rounding arc centre C 10 Nominal position of wheel 1 C 20 Nominal position of wheel 2 FEM Finite Element Method F i Force acting on contact point i K m Meshing stiffness PRi 1,2 Starting and finishing points for each planetring teeth couple R W Nominal primitive radius R i Displacement radius R j Load radius Rmin Minimum radius for the inner gear (tip) SPi 1,2 Starting and finishing points for each planetring teeth couple T ext External torque b k Node displacements of k flank b Pi (R-S) Coupling planet body b R(i-j) Coupling ring b S(i-j) Coupling sun d Geometric overlap k i Auxiliary angle Nowadays, there is a growing demand in the development of more reliable transmission systems, with higher requirements of torque, speed and compactness. In the race to satisfy these demands, engineers and researchers need to improve their understanding of the underlying phenomena involved in gear power transmission. One of the best approaches to accomplish this objective is the development of models intending to realistically reproduce the system behavior. Thereby, while research in this area contributes to improve engineer insight in gear transmissions, it simultaneously provides the industry with simulation tools that lead to more efficient design processes, to better maintenance practices [1] and to the development of other troubleshooting techniques based on vibration analysis [2]. Due to its spatial configuration, planetary transmissions are complicated to model. Nevertheless, being one of the most critical components in several aerospace and energy generation applications, planetary modeling research is booming, partly due to these reliability concerns [3]. One of the main advantages of planetary transmissions is its compactness. For high torques, instead of enlarging wheels size or using advanced materials, the simplest solution is to split the load into a number of paths, so that loadings per unit facewidth remain below nominal values while the torque is multiplied, always maintaining the same wheels size. Planetary transmissions use this approach, being the most compact and lightest possible drives [4].
Initial works on planetary modeling were focused on the extraction of natural frequencies and vibration modes through the use of planar lumped-parameter models [5,6], with constant meshing stiffness. Despite the simplicity that implies the use of a planar model, for spur gears where efforts and displacement can be considered in-plane, a two-dimensional model may be more than sufficient to efficiently and accurately study the transmission behavior [7]. Nevertheless, three dimensional analytical models have also been developed for the study of out-of-plane vibrations, especially when non-spur gears are used [8,9].
The main feature that characterizes the dynamic behavior of gear transmissions is the change in the number of teeth couples simultaneously in mesh. The meshing stiffness is therefore variable, and induces a periodic excitation in the system [10]. Thus, the characterization of this periodic excitation is crucial in order to achieve better simulated results [11]. In a first step to increase modeling realism the static transmission error has been used as excitation to predict dynamic behavior of planetary transmissions [12,13]. Nevertheless, recent studies point that this approach, while remaining relatively valid for ordinary transmissions, may not be applicable to multimesh transmissions such as planetary ones [14]. With a higher degree of accuracy, at a second step evolution there are planetary gear models with time-varying stiffness. They give better off-resonance responses, but they also are used to identify regions of large amplitude vibration near resonances, where damping and other nonlinear phenomena strongly affect the behavior [15,16]. The latest and more advanced planetary transmission models are those based on computational approaches, frequently including FEM techniques in combination with different contact models [17]. In some cases, completely flexible bodies are considered in real time simulation [18]. Depending on the particular application of the model, different emphasize is given to each modeling aspect [11,19].
In this work, the objective is to construct a planetary transmission model capable of obtaining a sufficient degree of accuracy in the simulated meshing stiffness, without impairing the dynamic modeling capabilities. With this aim, analytical solutions are hybridized with finite element models in order to compute the contact forces, making unnecessary the use of mesh stiffness waveform approximations or static transmission error excitation assumptions. The mesh model is based on previous work by the authors [20,21], extended and improved towards the planetary modeling. Coupling through gear body deformations is also given a special attention, due to the multiple meshes per wheel. With respect to the contact point location and geometric overlap modeling approach used in this work, it has been conceived to allow for the almost direct inclusion of additional modeling features, such as tooth profile modifications (with an approach used in [22]) or the use of shifted gears.
In this paper, especially key modeling aspects (such as geometric and contact forces calculation) are described in great detail, with the aim to allow for easily reproducibility of results and implementation by fellow researchers. For this reason and for the sake of shortness, only quasi-static results are presented, although the model is conceived to be part of a complete dynamic planetary transmission model.
A validation of the independent meshing stiffness (sun-planet and planet-ring) is provided by comparison with the ISO 6336-1 norm, as experimental results are scarce and difficult to compare to in the literature. Also comparison between modeling hypothesis adopted is presented, in order to quantify the maximum generated error.
Finally, a step-by-step integration of the planetary transmission is done, presenting results mainly based on the transmission error, which is one of the most representative variables of the planetary transmission behavior, using two typical configurations of a planetary transmission.

Model geometry
The first problem to be addressed in the transmission modeling is the appropriate description of the gear geometry. This definition must be consistent enough to represent with the maximum level of flexibility and robustness the widest possible range of parametric variations of the planetary transmission. In this section two main aspects related to the geometry problem are presented: gear profile definition and potential contact points location. For the first, it must be taken into account that the target is to simulate spur gears with involute flanks, but the model has to be open for the introduction of profile modifications and the use of shifted gears. Regarding the contact point location algorithm, the model must be capable of simulating the behavior of planetary transmissions with any number of planets, allowing the displacement of gear centers.

Gear profile definition
The profile definition of all gears (both external and internal) is made in this work by analytically mimicking the gear cutting process, following a vectorial approach similar to that described by Litvin and Fuentes [23]. For being this the most realistic possible procedure of defining the tooth profile, the model nearly matches the real gear. This analytical definition also implies a great adaptability of the model, allowing the use of shifted gears, also taking into account undercutting conditions. For the planet and sun gears the generating tool used is the rack shown in Fig. 1, defined by four sections. A coordinate system linked to the tool is placed at the center of the rack width space. The parameters that define the reference tool are the gear normal module (m), the pressure angle (u), the addendum (ad), the deddendum (dd) and the radius of the rounding tip (r). In the case of the ring shaping, the tool used is a pinion cutter, and its three section definition is shown on the right in Fig. 1, with the coordinate system linked to the tool placed at the center of the gear cutter. Once the coordinates of every section are defined in the reference system associated with each tool, the tooth profiles are obtained, expressing each of the considered sections in the gear coordinate system, simultaneously imposing the meshing condition. Following this procedure, the straight lines defined by u1 and u4 in Fig. 1 of rack generate the head and root circumferences of the gear, while points in the straight flank of the tool (section u2) create the involute profile. The same is valid for the internal generation, where sections u1 and u3 of the cutter pinion shape the head and root circumferences of the inner gear, while the pinion flank (section u2) define the involute profile of the internal teeth.
The tooth fillet of both external and internal gears are defined by the positions of section u3 and point M of the rack and pinion cutter respectively. In addition to all of this, a tip rounding arc has been added at the top of the teeth, following Vedmar procedure [24], in order to handle occasional contact at these points. Although these contacts occur out of the line of action and should be avoided through appropriate profile modification, the possibility must be taken into account, and the rounding arc is a realistic way of modeling the intersection between external arc and involute flank, where there will always exist some degree of rounding, desired or not. In this work, the profile modifications are not included deliberately, in order to show the consequences of its absence and to show the model performance dealing with contacts out of line of action (OLOA).

Contact point location
Once the wheels geometry is analytically defined, the second step in the modeling process is the location of potential contact points. One common approach to this end is the use of numerical methods that locate the points with smaller separation distances between profiles, later determining which of those pairs present real contact. In this model, being the profiles defined in a purely analytical form, it is possible to define an analytical procedure that provides the exact potential location of the contact points, with very little computational effort.
The meshing between two perfectly shaped and non-deformable gears always takes place along the straight line tangent to both base circles of the involved gears. The contact can be thus defined over this Line of action (LOA), and the separation distances (or geometric overlaps) can be calculated using the analytic formulation of the involute profile. Nevertheless, it must be taken into account that in the actual case of deformable wheels, a non-desirable contact may appear between involute profiles and rounding arcs during initial and final stages of the meshing period. Therefore, the contact point location algorithm must include this contact possibility.
Moreover, in order for the model to reproduce as closely as possible the actual behavior of the planetary transmission, contacts in the counterflank should be considered simultaneously, in order to model situations such as losses of contact, rattle or tooth wedging. In this section the formulation procedure for the location of contact points and the determination of the geometric overlaps is limited to the internal mesh between planet and ring gears, for being more complex than the external mesh case, which can be found thoroughly described in [20]. For internal gearing, with the introduction of tip rounding, the total number of contact possibilities is three: Involuteinvolute contact, Involute-external rounding arc and Involute-internal rounding arc.

Involute-involute contact
The Involute-involute contact takes place along the LOA defined by the tangent to both base circles, which coincides with the direction normal to the profiles, as shown in Fig. 2. As in the real situation wheels supports are not infinitely rigid, gears can move In the case of a planetary transmission, the movement of gear centers is not only defined by the support stiffness: in the planet case the gear center is fixed to the carrier, describing a circular trajectory. Besides, central members often float in planetary systems. For all these reasons the contact point location and geometric overlap calculation cannot be simplified with a fixed center assumption.
In Fig. 2a is presented the graphic construction corresponding with the internal pair at the reference position (contact at the pitch point). The planet gear is defined with the subscript 1 and the ring with the subscript 2, and a fixed reference system with origin O is placed at the system center. The initial position of wheel 2 C 20 (ring) coincides with the origin of the reference system, and the initial planet position C 10 corresponds to the point O 1 , located at the nominal mounting distance d 0 . The problem consists in determining the location of the contact point and the geometric overlap d when both wheels experiment center displacements ðx i ; y i Þ and rotations h i , all measured from the reference position, as shown in Fig. 2b.
Then, the geometric overlap between potential contact points on the involute profiles can be expressed as: According to this formulation, positive values of the distance d InÀIn correspond to situations in which the profiles are overlapped (from a purely geometrical point of view), thereby indicating contact situations. On the other hand, negative values of the overlap correspond to real separation distances between profiles. The obtaining of each term in (1) is made with base on the involute properties as Thus, the geometric overlap for involute-involute contact:

Involute-rounding arc contact
When considering the Involute-rounding arc contact, the normal to the profiles is now a line tangent to the base circle of the involute profile and passing through the center of the tip rounding arc C r , as shown in Fig.  3a. On this account, and to maintain the formulation as close as possible to the Involute-involute contact case, a new base circle (referred to as equivalent base circle) with radius q ri i can be defined. Because of the internal meshing that takes place in the planetary transmission, the symmetry of the formulation is broken, and it is necessary to define a distinct procedure depending on wether the rounding arc is located in the planet or in the ring. Thus, the superscript ri will be used to refer to those parameters affected by the rounding arc. For differentiation, i ¼ 1 when the rounding arc is in the planet, and i ¼ 2 when it is in the ring. For the first case, Involute-external rounding arc (in which the contact takes place between the ring involute profile and the planet rounding arc), the equivalent base radius q r1 1 can be defined, attending to Fig. 3a as Also from Fig. 3a it can be concluded that where arctan2 is a variation of the arctangent function, in which information about the sign is provided, and Similarly to Eq. (1) the new geometric overlap can be formulated as In which Q r1 i P i segments are now: Finally, substituting these expressions in Eq. (4), the geometric overlap for Involute-external rounding arc contact is For the second case, Involute-internal rounding arc, the contact takes place between the planet involute profile and the ring rounding arc, for which the equivalent base radius q r2 2 can be defined now as The formulation is now slightly changed, according to the distinct geometry of the problem. Attending to Fig.  3b it can be concluded that the auxiliary angle is now Correspondingly to Eq. (4), the geometric overlap is In which the new Q r2 i P i segments are: Finally, substituting these expressions in Eq. (9) the geometric overlap for Involute-internal rounding arc contact is

Contact type discrimination
In accordance to the above, there are three possible types of contact between ring and planet, depending on the presence of a non-involute part of the gear tooth in the contact. Thus, three different values of separation distance will be provided for each teeth pair considered during the implementation of the geometric overlap formulation. Nevertheless, depending on the planetary transmission position, only one of the possible contacts can be considered as a potential contact point, and the other two must be eliminated from the calculations. A possible solution to this problem would be to perform the calculation of the whole set of geometrical overlaps for the three possible contact cases. Once the distances were determined, the model could discard the two greatest values, the remaining contact being the valid potential contact point. Despite the simplicity of this approach, it would lead to the complete calculation of three geometrical overlaps for each position and teeth pair.
A much more efficient procedure consists in attending to the pressure angles for each contact, as shown in Fig. 4. In this way, a criterion can be defined for the discrimination of the prevailing contact. The Involute-involute contact will prevail when its pressure angle u T remains between u r1 T and u r2 T . When the contact takes place above the involute limit B, T . This inequality marks the prevailing of the Involute-external rounding arc contact. Correspondingly, a contact below C is always characterized by u T \u r2 T , which indicates an Involute-internal rounding arc type contact.

Multiple and reverse contacts
The procedure for locating the potential contact points and calculating the separation distance between teeth pairs has been established. It is now necessary to note that during the meshing process multiple contacts will occur simultaneously, and this must be taken into account. Furthermore, for the planetary model to be complete, it must also be capable of solving reverse contacts.
With respect to multiple contact, the total number N of potential contact points considered in this work depends on the contact ratio e, and it is obtained from where ceil is the ceiling function, which returns the smallest integer not less than its argument. From (9) it can be appreciated that an extra couple of teeth is added, unnecessary from a purely geometrical point of view. This is justified by the need to consider extreme situations in which great deformations may lead to additional contacts (due to the use of flexible materials or contact ratios near the upper integer). In Fig. 5 are shown the six potential contact points considered for a ring-planet mesh with a contact ratio below 2, with each contact point sequentially named, differentiating between direct (d) and reverse (r) contact. For the sake of simplicity, in the figure it is assumed that the wheels are located at the reference position (pitch point direct contact at 2d). Choosing the contact 2d as the one that defines the angular position of the wheels, to obtain the rest of geometric overlaps (contacts 1d and 3d) it is sufficient the addition (clockwise rotation) or subtraction (counterclockwise rotation) of the angular pitch h ip to the angular positions h i used in Sect. 2.2. All formulation described in said Sect. 2.2 has been derived in order to obtain the geometrical overlap between the tooth profiles for direct contact (positive slope of the LOA). Due to the symmetry of the problem, the inverse contacts can be converted to direct contacts, in order to use the same formulation. This conversion can be made attending to Fig. 6. As an example, reverse contact 1r can be found as a direct contact by simply rotating the wheels the angles h t 2 and h s 1 , corresponding with the internal tooth angular thickness and the external space angular width. To summarize the procedure, in Table 1 are shown the variable and parameter changes that must be introduced in the formulation to obtain the geometric overlaps for all possible contacts, both direct and reverse.

Contact forces calculation
To determine the contact forces in the planetary transmission for a given position, a relationship that links those forces and the deformations produced must be found. Subsequently, a non-linear system of equations subjected to certain conditions can be proposed, from which the contact forces are solved. It is therefore necessary to first obtain the relationship between forces and deformations. For this, a method derived from Andersson work [25] is used. This approach applies a similar procedure as that described by Vijayakar [26] for the formulation of the contact forces. Following previous work by Vedmar [24] for the static case, deformations at the contact point can be obtained as the composition of two different terms; the global deflection (hereinafter refer to as structural) and the local deformation as it can be seen in Fig. 7. The first term (structural) relates to the linear deflections in the gears far from the contact area, resulting from the bending and shear in the teeth and body torsion of the wheel. The second term (local) describes the nonlinear deformation of the contact area.
The calculation of the structural deformation is solved by applying FEM techniques, introducing an individual force at the contact point (Fig. 7a). As the results obtained by the complete finite element model under these loading conditions will only be valid for areas far away from the force application point, it is necessary to introduce a correction for nearby points. This is achieved by superimposing these results to those obtained from a partial finite element model loaded in the opposite direction (Fig. 7b). Finally, the local deformation is found by applying a non-linear closed-form analytical expression derived from the Hertz theory (Fig. 7c). As a summary, the method consists in the application of the superposition principle to the three independent problems shown in Fig.  7, taking into account Saint-Venant's principle for the determination of the boundary depth h, which must be large enough to provide negligible differences at the boundary between subproblems (a) and (b).

Local deformation
Local deformations are formulated using a Weber-Banashek proposal for two-dimensional problems, derived from the Hertz theory. Opposite to Vijayakar model, this proposal avoids the obtaining of the load distribution on the contact area. In order to achieve this, the contact area discretization and the integration of the Boussinesq solution is substituted by a closed   form analytical expression, thus assuming that both the contact area and pressure distribution adopt a predetermined shape. According to Weber-Banashek, the deformation between a point on the surface of a tooth and a point located at a depth h, for a plane strain hypothesis, is given by the expression: where F is the contact force, b is the width of the wheel and 2L is the length of the pressure distribution. In the case of plane stress, the above formulation is modified, resulting in: Regarding the length of the pressure distribution, it is obtained as a function of the load, depending also on geometric characteristics and material properties of the bodies in contact as: where E i , t i and v i are respectively the elastic module, Poisson's ratio and the curvature radius of the body i. Regarding the curvature radii, they are derived from the expressions obtained in Sect. 2.2, taking into account the contact type. The formula presented for the pressure distribution length estimation is correct for any pair of external gears, being the contacting surfaces always convex-to-convex. However, in the planet-ring mesh (Fig. 8), this is only true for the Involute-internal rounding arc contacts, whereas the Involute-involute and Involute-external rounding arc contacts present a greater level of surface conforming, with concave-to-convex surfaces, resulting in the modification of the pressure distribution length as follows: Therefore, it is necessary during the planetary modeling to consider this aspect, applying the adequate formulation for each case in the sun-planet and planet-ring meshes. Attending to the formulation, a sudden drop or rise in the L value is expected when the contact changes from Involute-involute to Involute-tip rounding arc and vice-versa. As this step variation is due to the discrete modeling approach and not due to real changes in the size of the contact zone, a linear transition of the curvature radii has been implemented.

Planetary mesh structural deformation coupling
As the structural deformation calculation approach is thoroughly described in [20] for the external gears case, being completely analogous to the one applied for the planetary transmission ring, only a brief description of the same will be given here. Attention will be focused on the main differences arisen from the introduction of internal gears and planetary arrangement. A two-dimensional FE model is constructed from the geometry obtained following the procedure  Fig. 9), on which nodes a unitary force F i is successively applied in the direction normal to the surface. Thus, for each load case, all b displacements for each k flank are obtained, from k ¼ 1; . . .; K. These values are then stored in K flexibility matrices (one for each flank, including the active one), indexing columns and rows by its corresponding displacement radius R i and load radius R j . Therefore, the flexibility values can be designated as b k R i R j , i.e., deflection of node j located at the flank k when the force is applied at point i of the active load flank. The procedure for obtaining the flexibilities in this model is generalized, so that not only is acknowledged the stiffness of the gear due to the flexibility of an individual tooth, but it is also taken into account the deformation of the wheel body. This approach has the effect of coupling the contact forces problem, introducing a cross effect between adjacent pair contacts. Thus, the dimension of the problem presents a dimension n, with the maximum in N, which is the total number of potential points of contact considered for each mesh. Nevertheless, the method used, which is completely valid for ordinary transmissions, present some new aspects when implementing a planetary transmission. The coupling of deformations between teeth pairs must now be extended to the rest of the gear meshes, given that all planets in the transmission mesh with sun and ring, and the latter two present as much meshes as the system has planets. Therefore, in order to calculate the contact forces, the deformation of the wheel bodies will produce not only coupling between the N potential contact points of each gear mesh, but also there will be a coupling between each gear mesh. This coupling is graphically represented in Fig. 10 for the three planet system. There will be a b PiðRÀSÞ for the coupling through each planet body deformation, as well as a b RðiÀjÞ and b SðiÀjÞ for the coupling via deformation of the ring and sun bodies respectively. The FE model shown in Fig. 9 for the structural deformations obtaining is therefore insufficient for the modeling of planetary transmissions. It is necessary in this case to extend the model to be able to obtain the cross-flexibilities between meshes, as shown in Fig. 11. The new flexibilities are stored in the same fashion, and the only change is the dimension of the problem, now enlarged. Thus, for a transmission with a number of planets P and N potential contact points for each mesh, instead of solving 2P systems of N equations, the new coupling transforms the problem into a single system of dimension 2PN. Of course the existence of negative forces is not possible (complementarity condition), so from the above it is concluded that although the size of the problem will have a maximum number of 2PN equations (36 for the case of 3 planets with 6 points of potential contact), the actual number of active contacts and therefore the final dimension n of the problem will be significantly smaller.

Deformation composition and numerical resolution
The total deformation u Tj in a contact point j can then be formulated as the sum of the local and structural contributions of each wheel: Both local and structural deformations also depend on position and orientation variables, that can be grouped in vector q i : For the sake of simplicity, this dependency has been removed from the equations. Besides, the local terms only depend on the force applied at the point of contact j, whereas the structural terms are also dependent on all contact forces acting on adjacent teeth pairs as follows (for wheel W1): where N is the number of potential contact points considered, F i is the force acting on the contact point i and k j;i is the deformation of the contact point j when a unitary force is applied at the contact point i. This linear relationship is implemented in the model simply by properly assigning the values of flexibility k j;i saved in the flexibility matrices b k , knowing that these values are stored in a sorted manner according to their R i and R j . The problem consists in determining the values of the contact forces F that verify the system: which can be reformulated for a certain position as: The non-linearity associated with the local deformation raises the need to apply an iterative method to obtain the forces F. The use of a fixed point iteration method is proposed, in which the resolution of f ðFÞ ¼ 0 is replaced by the obtaining of a fixed point of the function gðFÞ, problem defined as: The first stage of the contact problem resolution is shown in the diagram of Fig. 12. All the flexibilities are computed and stored, as well as the values of the geometrical overlaps d for all the potential contact points at the current position of the system. The necessary data k are extracted from the flexibility matrices, also depending on the wheels position. A priori, the dimension n of the problem is unknown, so it is assumed to match the number of active contacts from a purely geometric point of view. Thus, among all potential contact points N, only those n contacts in which d i [ 0 are initially considered. As shown in the resolution algorithm presented in Fig. 13, the dimension of the matrix k is consistently reduced from N Â N to n Â n. Subsequently, the system of nonlinear equations (15) is solved through the fixed point iteration method previously introduced. When the difference between the forces value for two consecutive iterations is below a tolerance value Tol, the algorithm continues. During this iterative process, the deformations coupling between different points of potential contact may result in the cancelation of contacts initially considered, or the appearance of new ones. This possibility is taken into account in the algorithm by means of two control loops, as shown in the block diagram of Fig. 13.
In the first loop it is checked the complementary condition, assuring that there are no negative forces which would correspond to a loss of contact among the initial contacts considered. On the other hand, the second loop evaluates the possibility of new forces emerging due to system deformation. To perform this check is therefore necessary to assess the effect of forces on the whole set of potential contacts N. Thus, the total deformations u T ðF k Nx1 Þ are obtained, and these values are subtracted to the geometric overlaps d Nx1 . The resulting differences are stored in the vector Dd k , whose positive components will then mark the potential points with actual contact. If these component indexes does not match for two consecutive iterations, it is necessary to consider the new set of contacts and begin the procedure again. On the contrary, if the indexes coincide for two iterations, the problem is considered solved and the contact forces are provided.
After obtaining the magnitudes of the contact forces it is necessary to define their directions, which are illustrated for the planet-ring pair in Fig. 14. The normal to the profile at the contact point is calculated as follows for direct contact: Whereas for the reverse contact: Finally, the set of contact forces acting on each wheel is transferred to their geometric center, yielding the vector On which the torque is obtained through multiplication by the base circle radii.

Planetary transmission integration
It is true that the resolution of the presented non-linear system may imply excessively long calculation times,  The integration of the complete planetary transmission is next presented step by step. First the meshing stiffness is validated for the planet-ring and sun-planet pairs. A transmission with only one planet is then studied, and finally the full system is implemented. As example, a real gear planetary reducer from agricultural machinery is modeled, whose main parameters are shown in Table 2. The real application consists of two stages with a common ring, reason for which the ring width is larger than the rest of the wheels. It is also to be noted that the gears in question are non-standard.

Transmission error (TE)
The most basic problem that can be formulated is the obtaining of the angular position of a pair of wheels for which the static equilibrium presented in Eq. (21) is fulfilled. The TE problem is graphically shown in Fig.   15 for the planet-ring pair. The problem definition remains the same for the sun-planet pair. Similarly, a generalization of the TE is made for the planetary transmission.
Initially, the pair of wheels in contact is positioned at O 1 and O 2 . Wheel 2 is completely fixed at angular position h 2 . Thus, the one single degree of freedom of the problem is the planet rotation h 1 , to which is applied the external torque T ext . The solution to the problem consists in finding the rotation h 1 for which the geometric overlap generates forces such as the external torque is balanced. The equilibrium position h 1 will then be composed of two terms indicated in the following equation where the first term corresponds to the value of a strictly kinematic rotation, due to the ring angular position, and the second term represents the extra rotation (which implies deformation) generating contact forces. The difference Dh 1 between the theoretical (or kinematic) angular wheel position and the actual position is called transmission error:  The resolution of the problem posed in (21) is made by applying a variation of Newton's method with trust regions, directly implemented in the modeling platform. For the external torque applied the model tries to find the position of the wheels for which the meshing forces satisfy the equation, also considering the DOF fixed in the particular case. As a first application of the model, the TE for both gear pairs has been obtained, and it is presented in Figs. 16 and 17. It is easy to observe the parametric excitation due to the variation in the number of teeth pairs in contact, and how the error amplitude increases as the load does likewise. Another effect of the load increase is the reduction of the single contact area, thus altering the contact ratio. This is because greater torque implies greater deformations, which bring the mating profiles closer, advancing the contact start and retarding its end. Considering the formulation of local deformations described in previous sections, the appearance of a singularity in the transmission error figure can be expected in the form of a sharp step when the linear transition between involute-involute and involute-rounding arc contacts is suppressed. This phenomenon is noticeable in the data shown in dashed lines for the greatest torque in Fig. 17.

Meshing stiffness (K m )
An additional magnitude which can be defined on the basis of the transmission error is the meshing stiffness. It is usually assumed as one lineal spring, acting on the line of action. Thus, the possibility of considering one stiffness for each contact as well as one equivalent single stiffness for the whole mesh can be adopted. The equivalent stiffness is defined in this work as: It should be noted that while this magnitude is a useful parameter facing the linearization of transmission models, it must be used carefully. Including the transmission error in its definition, it is possible that the stiffness values obtained are underestimated, some of the TE being attributable to non-elastic deflections such as kinematic or mounting errors. If so, the dynamic behavior of a linearized model using that Transmission error [rad] Transmission error [rad] Fig. 17 Transmission error for the planet-ring pair false meshing stiffness would not correspond to the reality. For this reason, the subindex elast at Dh 1 elast points that only elastic components of the TE can be included in the stiffness calculation. As in this work the TE only includes deviations due to deformation (the gear centres are fixed and there are not any kind of error implementation), this formulation is direct. In a different case scenario, the meshing stiffness must be defined not based on the rotational TE, but on the real overlap between mating profiles. In order to assess a partial verification of the planetary model in terms of meshing stiffness, the values obtained from the model have been compared to those provided by the ISO 6336-1 Calculation of load capacity of spur and helical gears in its clause 9, using method B. This method is based on studies of the elastic behavior of solid disc spur gears, with its accuracy being thoroughly verified by measurement results. As shown in Figs. 18 and 19 the average value of stiffness obtained from ISO 6336-1 fits the stiffness values for both gear pairs provided by the model.
Besides the mean value, the shape and peak to peak variation of the meshing stiffness are factors which greatly affect the dynamic performance of the transmission model. Thus, a comparison between the proposed model meshing stiffness shape and several other formulations has been performed and it is presented next.
For the usual spur gear transmissions with a contact ratio between 1 and 2, the theoretical meshing stiffness expected shape should resemble a square wave, with the double contact period more or less longer than the single contact period, according to said contact ratio. Regarding the maximum and minimum level of stiffness, the minimum level corresponds to the single contact period, while for the double contact the stiffness reaches its maximum. During the double contact period two pairs of teeth share the load, and thus the stiffness can be represented as the parallel combination of the two individual pair stiffnesses. From this it can be concluded that the maximum value of the meshing stiffness for the double contact period cannot exceed twice the value of the individual pair stiffness (single contact period). Nevertheless, several factors cause the maximum value to be significantly less than this theoretical limit. More importantly, the load is shared among two pairs, thus the individual contact force is reduced, and the individual stiffnesses decreases accordingly, due to the non-linearity of the contact problem. As additional factor, the contacts during the double contact period take place far from the primitive point (center of the flank teeth), where the combination of teeth stiffneses is greatest. As result of this, the theoretical limit for the peak-to-peak amplitude of the meshing stiffness should be 2=3 with respect to the mean value of the maximum and minimum stiffneses. In Fig. 20 the meshing stiffnesses of different external pairs extracted from several works in the literature are shown. All stiffnesses have been normalized with respect to their mean value, and all plots are centered at the 0 y-axis for easy comparison. Although the transmission parameters Meshing  Table 3), all of them consist in spur gear transmissions, for which the shape and limits of the normalized stiffness values can be compared. From Fig. 20 the different approaches can be assessed in terms of the shape and amplitude variation of the meshing stiffness. Pedersen et al. [27] employs a hybrid approach with planar FE models and analytical formulations similar to the presented model, but with a number of simplifications. The coupling of deformations is not considered, and the body deformation is underestimated. The non-linear effect coming from the load level is also neglected. Hui-Ma approach [28] is based on elastic mechanics and beam theory, modeling the teeth as a nonuniform cantilever beam on the root circle. Although using a completely different methodology than Pedersen, again the deformation coupling is not considered. This fact account for the greater peak-to-peak amplitude shown in both meshing stiffness plots, compared to the proposed model results.
Two of the compared models provide similar results in terms of meshing stiffness variation to the ones in the present work. Song He et al. [29] uses a finite element/contact mechanics (FE/CM) analysis code which in essence is based on the same principles that the model presented here. Song-He peak-to-peak variation is almost coinciding with Iglesias results, and the main differences in shape are due to dynamic fluctuations in the simulation (the stiffness is not extracted from a pure quasi-static analysis as in this work) and differences in the transmission data. Ambarisha's approach [30], is very similar to the previous one. In this case, the peak-to-peak variation is again similar to the one described in this work, but the differences in shape are due to changes in the load applied to the gear pair. The simulation is quasi-static, but the effect of load sharing variation between paths affect the torque applied, introducing certain waviness.
The two models with smaller peak-to-peak amplitudes are those of Kiekbusch et al. [31] and Howard et al. [32]. Both are based on pure FE analysis, with line-to-line general contact elements. Although both works are remarkably well done, the fact that the meshing stiffness provided is obtained through the relationship between the transmission error and the applied torque could lead to some degree of stiffness corruption (for example with the existence of counterflank contacts due to small backlash width), even if the rest of the data are correct.
Regarding the meshing stiffness for internal pairs, the scarcity of results in the literature makes very difficult to compare. In Fig. 21 the normalized meshing stiffness values for Pedersen et al. [27], Ambarisha et al. [30] and Iglesias are shown. The peak-to-peak variation differences are consistent with the causes described for external pairs. In this case, the coupling deformation for the ring gear in Iglesias' model is lesser than in Ambarisha's work due to the  Fig. 20 Meshing stiffness comparison for external gears boundary conditions used, which may contribute to the greater differences. Regarding the shape, it is also clear that Pedersen approach (as well as Hui-Ma for the external case) provide an unrealistic meshing stiffness, with a purely kinematic transition between single and double contact periods. In view of the results, the proposed model seems to fit well the contrasted meshing stiffness results found in the literature in terms of mean value, shape and peak-to-peak variation. Nevertheless, an experimental validation is yet necessary to assure the model performance.

Deformation coupling
Among the advantages of the presented model for the contact forces calculation it is included the consideration of coupled deformation between adjacent teeth pairs. When a tooth is loaded, the entire body of the wheel is deformed accordingly, including the profiles on adjacent teeth, thereby affecting the contact ratio and meshing stiffness. The effect of coupling between tooth deformations can be observed in terms of TE in Figs. 22 and 23 for the external and internal pairs respectively for an external torque of 100 Nm. Thus, during double contact the meshing stiffness drops due to the deformation coupling, important aspect overlooked with the use of other approaches. The reason for which the coupling of deformations can not be considered in the majority of the simplest models is the consideration of a single equivalent meshing stiffness. In this paper the meshing stiffness is one of many results of the simulation, and at every moment the contact forces are known for each teeth pair, which represent a valuable piece of information in terms of mechanical design. In Fig. 24 the complete evolution of the load carried by a teeth couple is shown. The load is normalized by the theoretical maximum force transmitted by a single teeth pair. The data thus obtained is called tooth load factor, and provides information about the variation of the load transmitted by one teeth pair. The model is capable of provide valuable insight on how the transmitted torque affects the load factor. Although the differences may seem little with respect to other models, the tooth load factor or tooth load sharing is determinant in efficiency calculations for gear transmissions, and the accuracy with which this magnitude is provided is crucial. It is also important to take into account that the proposed model has been constructed having in mind the introduction of shift factors and profile modifications, cases for which the tooth load factor would differ greatly.

Plane elasticity
In order to justify the plane elasticity assumption made in the contact forces model, the limits of the error were quantified and assessed. In many of the works on spur gears the general rule is to simplify the problem, leading to the use of a plane hypothesis. Wang and Howard [33] conducted an extensive study based on finite element models, in order to quantify the error made with the assumption of a plane elasticity hypothesis. In Fig. 25 the results of the relative error found by Wang is particularized for the external wheels (S for sun and P for planet) used in this work. The x-axis in Fig. 25 corresponds to the gear axial width to pitch radius ratio R, and on the y-axis the meshing stiffness relative error is presented for each hypothesis. The reference for the error percentage is the value found in a fully complex three dimensional FE analysis. Attending to Wang results for the example planetary transmission used in this work, the relative difference between hypotheses will be approximately 6 %, and the difference with the actual stiffness of only 3 %. For real hypothesis validation in this work the maximum error has been estimated through model comparison. The percentage differences between the two scenarios for an external torque of 100 Nm are consistent with Wang study, with a 7.6 % difference for the planet-ring pair and 7.3 % for the sun-planet. These findings point to an error in the estimation of the meshing stiffness for the studied gear couples of less than AE4 %.

Sun-planet-ring assembly
As a second step towards the complete integration of the proposed planetary transmission model, the assembly formed by the sun, ring and a single planet is studied, as shown in Fig. 26. Only rotational degrees of freedom are considered, and an external torque is applied. The equilibrium is found for the fixed ring and different positions of the planet carrier. When the assembly is formed, and attending to the example transmission characteristics, there is one aspect that demands special consideration. The operating pressure angles of the sun-planet and planetring pairs are significantly different. In a planetary transmission, for all wheels to work at nominal pressure angle conditions, it must be met the following relationship between nominal primitive radii: In Fig. 27 this relationship is graphically represented. Sun, planet and ring base circles are presented in blue. Equation (25) could be rewritten in terms of teeth, obtaining the theoretical number of ring teeth for a given sun-planet pair. In this case 16 þ 2 Ã 24 ¼ 64.
Since the actual number of teeth in the ring is 65, it is thus not possible to locate the gears observing the  Fig. 28 is shown the configuration for the actual case, where the planet can be placed in a certain range between the nominal position with respect to the sun (blue dot) or with respect to the ring (red dot). The actual position is located in this case at the 86:4 mm from the system center, and marked with a green dot. Disregarding the actual position and attending to the figure, it is apparent that the operating pressure angles will always be greater than nominal for the sun-planet pair and less than nominal for the planet-ring pair. Two different operating pitch circles will appear (orange for the planet-ring contact and pale blue for the planet-sun contact).
From the modeling point of view, this is an aspect of great importance if the planetary model is supposed to flexibly deal with this kind of non-standard transmissions. All the formulae for geometrical overlap calculation are based on a reference angular position, defined by the contact at the pitch point. Having now two different pitch circles, an extra difficulty is consequently added. From the operating standpoint, a direct consequence of this particular configuration is the appearance of a radial force component on the planet. The new operating pressure angles for each gear pair are calculated for the planetary transmission as u PÀR ¼ 24:45 ; u SÀP ¼ 27:37 , and in Fig. 29 the radial component of the meshing forces on the planet for an external torque of 300Nm is presented. In this figure it can be observed the effect of the contacts out of the LOA, when the tips of the teeth come into contact. For the sun-planet-ring assembly the period fraction with rounding contacts reaches up to 45% of the total meshing period. This tip rounding contact duration length may seem very high, but it is due to the absence of profile reliefs and due to the fact that two different meshes (sun-planet and planet-ring) are considered together. For a single pair of gears, this period fraction would be significantly smaller. Of course, the load level also plays an important role in this aspect. The contacts OLOA take place during the change in the number of pairs in contact, when the tip of the teeth engage or disengage. They occur because of the extra rotation of the tooth profiles due to the gear body deformation, allowing points beyond the  imaginary involute curve (e.g. rounding tip section) to come into contact. To avoid these undesirable contacts, reliefs at the tip or root of the mating profiles must be included.
As to the value of the TE of the sun-planet-ring assembly, it is easy to appreciate the combination of the two gear pairs stiffness. The formulation of the TE is now referred to the sun position, which is: where the first two terms corresponds to the two DOF which determine the rotation of a planetary transmission. These two terms are strictly kinematic, and the last term represents the extra rotation suffered by the sun due to all deformation occurred in the system. The difference Dh S between the theoretical (or kinematic) angular sun position and the actual position is now the TE for the sun-planet-ring assembly: In Fig. 30 it is shown the TE for three meshing cycles and a range of external torques. At the initial carrier angular position, the sun-planet mesh presents one teeth couple in contact, while the planet-ring presents two. At h carrier ¼ 0:025 rad the planet-ring mesh drops to one teeth pair in contact. Later a pronounced step due to a nearly simultaneous entry of a new teeth couple for each mesh occurs. The difference between both entries is noticeable for low torques, while the coupling deformation for higher torques blurs the two consecutive steps. The meshing cycle concludes with the finish of an external teeth couple contact, resuming the initial situation.

Planetary transmission
Once simulated and validated (in terms of stiffness) the performance of the model for external and internal mesh simulation, as well as the combination of sun, planet and ring, the study of the transmission modeling is extended to the full planetary. The three planet example transmission is modeled, and two case configurations are approached: the fixed configuration, with only rotational degrees of freedom, and a floating configuration, in which the sun can move freely.
The static equilibrium for a given position (always considering the ring and carrier angular position known) can be posed for the fixed configuration through a 4 equation system as When considering the floating configuration, it is necessary to determine the position of the center of the sun, so the system of equations must be extended with the balance of forces in the sun:

Fixed configuration
Analyzing the TE for the fixed configuration (which can be determined following Eq. [27] it can be appreciated in Fig. 31 how the resulting TE corresponds to the triple combination of the TE shown in Fig. 30. The different starting and finishing points for each teeth couple are marked as PRi 1;2 (planet-ring) and SPi 1;2 (sun-planet), where the i refers to the planet index and the subscript inform about the number of teeth couples in contact.
The effect of torque increase on the planetary TE is threefold: it changes the mean value, the peak-to-peak amplitude and the shape of the lobes, as presented in Fig. 32.

Floating configuration
In planetary transmissions, one of the main concerns is the equilibrium among planet load levels. Due to the variable stiffness and different manufacture and mounting errors, the planets are often unequally loaded. To avoid this, several methods are available. One of the most simple and easily implemented solutions is to allow one central member to move freely, without supports. To prove the model capabilities, a floating sun configuration has been tested, and the results in terms of TE and sun orbit are presented.
The sun describes and orbit that compensate for the stiffness variation among load paths. For the example transmission, with a load level of 300 Nm the orbit radius is 3 lm, as shown in Fig. 33. This small sun displacement also causes an increment in the TE mean level, as the system tends to the least stiff possible position. This fact is shown in Fig. 34, where the TE for fixed and floating configurations is compared.

Planetary mesh deformation coupling
Among the model virtues, one of the most significant ones is the capacity of deformation coupling consideration. The deformation coupling between adjacent teeth pairs in the same mesh has been studied and presented in Sect. 4.1, as well as the model generalization necessary to take into account coupling between meshes (in Sect. 3.2). In Fig. 35 is shown a comparison of the planetary TE, representing four different cases. In blue is depicted the TE without consideration of mesh coupling. In red is depicted the coupling through planet body deformation, whereas in green color it is shown the coupling through sun deformation. Finally, the TE for both planet and sun couplings is represented in black.
The ring body coupling is not shown, due to its insignificant effect on the TE for the example transmission used. This lack of effect is attributable to the construction and mounting of the real transmission, which demands very specific boundary conditions imposed to the ring FE model. The outer nodes are completely fixed, which greatly impedes gear deformation. For other applications, especially flexible rings are used to improve the load sharing between planets. Those implementations would give very different results in terms of ring deformation coupling effect compared to the presented ones, which are focused solely on the gear torsional deformation.
Regarding the coupling effect and the case comparison, the planet contact forces tend to bring closer the tooth profiles of the opposite planet mesh, which has a stiffening effect on the overall stiffness. This expected behavior is found in red in Fig. 35, with a reduction of the TE mean value. On the contrary, the coupling deformation through the sun body has the opposite effect. This coupling is equivalent to the deformation coupling between teeth pairs in the same gear mesh. Contact forces tend to deform the body, separating the rest of the contacting teeth pairs, which causes an increase of the contacts flexibility. Again, this behavior can be confirmed by means of observing the greater TE mean value for the green case, which corresponds to the sun coupling. Regarding the overall consequences, the resultant TE when all couplings are considered (shown in black) shows a drop in the planetary meshing stiffness, due to the prevailing effect of the deformation coupling through the sun over the planet.

Conclusions
A mathematical tool has been developed to simulate the inner gear cutting process, where the generation of the trochoid is implicitly included, as well as the possibility of generating shifted gears. A tip rounding arc has been added at the top of the teeth, in order to avoid singularities during occasional contact at these points. This model feature provides good results, but the approach applied to calculate the local deformations introduces a step in the transition between different contact types, due to the change in the curvature radii.
The analytical model used to calculate the local deformation has been extended to consider the different conditions of curvature which can be found in the planetary transmission. The problem of the step between contact types has been solved with the inclusion of a linear transition between zones, resulting in a significant improvement of the model performance.
The algorithm used for the calculation of the meshing forces has been broadened and extended with respect to previous versions, expanding the number of DOF to 6?3n for a planetary transmission with n planets. Furthermore, the coupling of deformations through the gear body flexibility has been also considered.
While the coupling through the planet deformation causes a torsional stiffening of the system, the deformation coupling of the n meshes of the sun causes a lessening of the system torsional stiffness, equivalent to the effect of coupling between adjacent teeth. This effect has a great impact, partially due to the proximity of meshes, and therefore it is increased with the number of transmissions paths. The coupling through the ring deformation has a negligible effect on  the transmission behavior, due to the boundary conditions used in this study. The contact forces model has been assessed both in inner and external gears. The procedure followed included the comparison of the average meshing stiffness obtained with the model and those defined by the international norm ISO 6336-1-(2006), and the comparison in terms of shape and peak-to-peak amplitude with other models' results, the proposed model fitting well. The choice of applying a hypothesis of plane stress or strain during the calculation of contact forces has an impact lower than 4% on the average meshing stiffness.
A step-by-step integration of the planetary transmission have been done, and the model capability for transmission simulation of a planetary with distinct pressure angles is demonstrated and shown. The model has also been put to the test with the simulation of the transmission error of a real transmission system, and results for varying external torque demonstrate the model possibilities.