Gear Teeth Deflection Model for Spur Gears: Proposal of a 3D Nonlinear and Non-Hertzian Approach

Abstract

In this paper, a three-dimensional model for the estimation of the deflections, load sharing attributes, and contact conditions will be presented for pairs of meshing teeth in a spur gear transmission. A nonlinear iterative approach based on a semi-analytical formulation for the deformation of the teeth under load will be employed to accurately determine the point of application of the load, its intensity, and the number of contacting pairs without a priori assumptions. At the end of this iterative cycle the obtained deflected shapes are then employed to compute the pressure distributions through a contact mechanics model with non-Hertzian features and a technique capable of obtaining correct results even at the free edges of the finite length contacting bodies. This approach is then applied to a test case with excellent agreement with its finite element counterpart. Finally, several results are shown to highlight the influence on the quasi-static behavior of spur gears of different kinds and amounts of flank and face-width profile modifications.

1. Introduction

Due to their inherent nature, load and stiffness fluctuations are the main source of excitation and cause of failure in geared transmissions [1,2,3]. Early experimental works from the Japanese school highlighted the interesting dynamic features of those systems, analyzing their torsional behavior including other sources such as profile errors or modifications [4,5,6,7,8]. At first, the dynamic factor (DF) was used to characterize and explain certain types of failures [9] by comparing the dynamic loads in operation with those under static conditions. Vibrations and impacts can certainly be traced back to variations in the input torque depending on the machine characteristics [10,11,12], but the time-varying mesh stiffness (TVMS) was quickly found to be a key player [13]. Indeed, the stiffness of an engaging gear pair changes continuously due its nature and can generate self-excited vibrations which led many researches in the study of the transmission error (TE). Over the years several methods have been proposed, starting from integral approaches [14,15] or discrete ones [16,17]. Others considered the tooth as a trapezoidal beam [18] while also semi empirical formulas have been proposed [19,20]. Experimental methods were also proposed to study the static transmission error (STE) [21,22,23], while others included mounting and manufacturing deviations to estimate the manufacturing transmission error (MTE) [24,25,26,27]. Finite element (FE) models were obviously proposed as well, but its time-consuming nature and difficulty to set up made it applicable to limited aspects, such as tooth root stresses and its structural behavior [28,29,30], crack propagation [31,32,33,34,35], or generally as a validation tool for other proposed models. Hertzian theory [36] of cylinder-to-cylinder contact is generally employed to model the contact between engaging flanks, simplifying several key aspects of the gears, such as the continuously varying curvature and the presence of sharp edges. When contact is not neglected it is generally introduced as an addition to the elastic behavior of the mating teeth, but still under Hertz’s hypotheses [37,38,39,40,41], while other works included non-Hertzian properties, but neglected the flexibility under load of the teeth [42]. Hybrid approaches coupling a FE model with a semi analytical (SA) contact model have also been proposed [43] with great success.

In this paper, a combination of some of the presented approaches will be used to accurately estimate the STE and the contact conditions, including the influence of the micro-geometrical modifications of the teeth flanks both along its height and along the face-width. An SA model [17] will be used to determine the stiffness of the different mating pairs. No a priori assumptions will be made regarding the location of the contact point, as well as the number of mating pairs bearing the load and the load sharing among them. The rigid body kinematics will only be employed as the first tentative guess for a nonlinear iterative scheme, in which, also, the working pressure angle will be dependent on the deformation. A natural equilibrium condition is sought for the location, number, and load intensity acting on the contact points found by a surface-to-surface intersection algorithm. A three-dimensional non-Hertzian contact model is the employed to correctly model the interaction between the meshing flanks. Said model is first compared to Hertzian theory with great agreement, and side and tip mirroring corrections are then introduced to relieve the stresses on the free surfaces of the finite-length bodies in contact allowing accurate representation of the varying curvature and discontinuities of the flanks. After defining the theory for the generation of the tooth geometry tooth profile modifications (TPM) are introduced, as well as the ray tracing approach used to define the contact points [44]. Next the SA model for tooth deformations is briefly described, along with the nonlinear iterative scheme used to overcome the rigid kinematics approach. Next, the non-Hertzian contact model is shown along with some results before introducing the concept and the effects of plane mirroring to correctly model bodies of finite length in contact. The deformations of the contact surfaces are then included as superposition to the previous results to finally obtain the STE. This entire approach is then tested against a 3D FE model showing very similar results both in terms of STE and contact pressures. Next, several results are shown on a test case to highlight the influence of micro-geometrical modifications on the quantities of interest and finally conclusions are drawn and future work is introduced.

2. Model Description

2.1. Tooth Geometry and Contact Points Detection

The detailed tooth geometry is obtained by applying Litvin’s vector approach [45] by simulating the meshing interaction between a rack cutter and the gear blank upon which TPM are then applied. Among the possible modifications along the tooth height the main ones are tip and root relief, while along the face-width crowning modifications are possible. Tip and root relief are defined by the diameter at which the modification begins ($dCa$ and $dCf$, respectively) and by the amount of material removed $Ca$ and $Cf$, respectively, as visible in Figure 1. The amount of material removed can follow a linear law or a parabolic one. The crowning is a parabolic removal of material from the tooth face-width along its entire height. The amount of material removed at the beginning and end of the face-width is marked as $C_{\beta In}$ and $C_{\beta Fin}$ in Figure 1 and is evaluated at a distance $zLc$ along the face-width, while the axis of the parabola is defined by the coordinate $zL$. If the helix deviation $C_{H\beta }$ is different from zero, the profile of the face-width is obtained by the summation of the previously defined parabola and a further amount due to the line inclined of angle $C_{H\beta }$.

In this work, the location of the contact point is not defined kinematically by finding the location of the intersection of the line of action (LOA) with the tooth profile at different angular positions, but it is instead sought for by a surface to surface intersection method on the rigid profiles for the first iteration, while on the deformed ones for the successive iterations of the iterative algorithm that will be described later. The surfaces of the teeth are triangulated as visible in Figure 2 and then rotated of a small angle (usually ${10}^{-5}$ rad) towards each other around their center of rotation and resulting intersection points are found through the application of the Möller–Trumbore ray-triangle intersection algorithm in which each triangle edge of each surface is treated as an infinitely long ray and an intersection is sought with the triangular faces of the opposite surface as visible in Figure 3 (for more details on that methodology, please refer to [44,46,47]). If an intersection is found between the ray and the surface of one of the triangles a check on the actual intersection is made by verifying if the location of the intersection point $p_c$ lies between the endpoints of the considered edge $(p_1,p_2)$ by

$$\left|\right|p_c-p_1\left|\right|+\left|\right|p_c-p_2\left|\right|=\left|\right|p_2-p_1\left|\right|$$

(1)

Usually this kind of ray tracing is computationally heavy and almost unusable without complex and specific hardware GPU acceleration, but its speed can be extremely improved by using an octree data structure [46,47] to partition the triangles in which the groups of triangles are recursively subdivided in eight bins, thus greatly reducing the computational effort by a huge amount and making it usable in an iterative process, such as the one being presented without slowing it down. Once the points of intersection between the profiles are known, the location of the contact point $\mathbf{O}={\{x_{cp},y_{cp},z_{cp}\}}^{{}^{\prime }}$ is obtained as the average coordinates of all intersection points.

2.2. SA Deflections

Since only spur gears and no misalignment are considered in this paper, the teeth profiles are discretized in $N_i$ points at the middle of the face-width. Bending and shear properties of each tooth are computed by considering a clamped-free beam model with non-uniform geometry as by Cornell [17] instead of the integral approach of Weber [14]. By applying a load $F_j$ [id = r2] (total load exchanged by the considered flanks) the expression of deflection is

$${\delta }_b^i=\frac{F_{j}co{s}^2{\mu }_{\alpha }}{E}\sum _{i=1}^{N_i}{\delta }_i\left[\frac{l_i^2-l_i{\delta }_i+\frac{1}{3}{\delta }_i^2}{{\overline{I}}_i}+\frac{2.4(1+\nu )+ta{n}^2{\mu }_{\alpha }}{{\overline{A}}_i}\right]$$

(2)

where ${\delta }_i$ is the thickness of the $i^{th}$ slice of the tooth cross-section defined by two consecutive points i and $i+1$ ($i=1,2,\dots ,N_i$) of the profile. E and $\nu$ are, respectively, the Young’s modulus and Poisson coefficient of the material, while ${\overline{A}}_i$ and ${\overline{I}}_i$ are the area and moment of inertia of the slice cross section evaluated at the middle point of the slice along tooth height, $l_i$ is the distance of the point of application of the load $F_j$ from root radius of the tooth. ${\mu }_{\alpha }$ is the working pressure angle which will be in general different from the rigid counterparts $\alpha$ due to teeth deflections and it is computed at the location of the contact point, as visible in Figure 4. Since all triangles considered $\left(\widehat{OAB},\widehat{OBC},\widehat{OCD},\widehat{ODA}\right)$ will have in general a different normal vector the average of them is found by constructing an orthogonal reference frame with origin in the contact point. Two reference frames are defined, one referred to the single triangle and, hence, a local frame (X, Y, Z) and a global one (x, y, z). The vector defining the local X axis is defined by

$${\mathbf{v}}_{\mathbf{x}}={\mathbf{u}}_{\mathbf{j}}-\mathbf{O}$$

(3)

where ${\mathbf{u}}_{\mathbf{j}}={\{x_j,y_j,z_j\}}^{{}^{\prime }}$ and $j=A,B,C,D$ alternatively. An accessory vector is defined by

$${\mathbf{v}}_{\mathbf{a}}={\mathbf{u}}_{\mathbf{k}}-\mathbf{O}$$

(4)

where ${\mathbf{u}}_{\mathbf{k}}={\{x_k,y_k,z_k\}}^{{}^{\prime }}$ and $k=B,C,D,A$ alternatively. The vectors defining the Y and Z axis are then obtained as

$$\begin{array}{c}{\mathbf{v}}_{\mathbf{z}}={\mathbf{v}}_{\mathbf{x}}\times {\mathbf{v}}_{\mathbf{a}}\hfill \\ {\mathbf{v}}_{\mathbf{y}}={\mathbf{v}}_{\mathbf{z}}\times {\mathbf{v}}_{\mathbf{x}}\hfill \end{array}$$

(5)

The versors are then obtained dividing the vectors with their norm

$${\mathbf{n}}_{\mathbf{x}}=\frac{{\mathbf{v}}_{\mathbf{x}}}{\left|\right|{\mathbf{v}}_{\mathbf{x}}\left|\right|},{\mathbf{n}}_{\mathbf{y}}=\frac{{\mathbf{v}}_{\mathbf{y}}}{\left|\right|{\mathbf{v}}_{\mathbf{y}}\left|\right|},{\mathbf{n}}_{\mathbf{z}}=\frac{{\mathbf{v}}_{\mathbf{z}}}{\left|\right|{\mathbf{v}}_{\mathbf{z}}\left|\right|}$$

(6)

Finally, the averaged versor defined with respect to the global reference frame and normal to the surface in the contact point is simply obtained by

$${\overline{\mathbf{n}}}_{\mathbf{z}}=\{{\overline{n}}_{z,1},{\overline{n}}_{z,2},{\overline{n}}_{z,3}\}=\overline{\{{\mathbf{n}}_{\mathbf{z},\widehat{\mathbf{OAB}}},{\mathbf{n}}_{\mathbf{z},\widehat{\mathbf{OBC}}},{\mathbf{n}}_{\mathbf{z},\widehat{\mathbf{OCD}}},{\mathbf{n}}_{\mathbf{z},\widehat{\mathbf{ODA}}}\}}$$

(7)

From the components of the ${\overline{\mathbf{n}}}_{\mathbf{z}}$ versor the actual pressure angle ${\mu }_{\alpha }$ can be found by

$${\mu }_{\alpha }=ta{n}^{-1}\left(\frac{{\overline{n}}_{z,2}}{{\overline{n}}_{z,1}}\right)$$

(8)

Additionally, the compliance of the fillet region of the tooth and its foundation is included through an analytical expression based on the theory of elastic rings as [48,49]

$${\delta }_f^i=\frac{F_{j}co{s}^2{\mu }_{\alpha }}{bE}\left[L^{*}{\left(\frac{u}{s_f}\right)}^2+M^{*}\left(\frac{u}{s_f}\right)+P\ast \left(1+Q\ast ta{n}^2{\mu }_{\alpha }\right)\right]$$

(9)

where b is the facewidth of the gear, u is the distance between the tooth base and the point of application of the load on the tooth centerline, while $s_f$ is the chordal resisting section measured at tooth root. The coefficients $L^{*}$, $M^{*}$, $P^{*}$, $Q^{*}$ are computed as

$$X^{*}(h_{f,i},{\theta }_f)=A_i/{\theta }_f+B_i{h}_{f,i}^2+C_i{h}_{f,i}{\theta }_f+D_i/{\theta }_f+E_i{h}_{f,i}+F_i$$

(10)

where $h_{f,i}$ is the ratio between the root circle and hub radii, while ${\theta }_f$ is half of the tooth root angular thickness. The coefficients $A_i$, $B_i$, $C_i$, $D_i$, $E_i$, $F_i$ are given in

Table 1. Polynomial coefficients for Equation (10).

	${\mathit{A}}_{\mathit{i}}$	${\mathit{B}}_{\mathit{i}}$	${\mathit{C}}_{\mathit{i}}$	${\mathit{D}}_{\mathit{i}}$	${\mathit{E}}_{\mathit{i}}$	${\mathit{F}}_{\mathit{i}}$
$L^{*}$	$-5.547\times {10}^{-5}$	$-1.9986\times {10}^{-3}$	$-2.3015\times {10}^{-4}$	$4.7702\times {10}^{-3}$	0.0271	6.8045
$M^{*}$	$60.111\times {10}^{-5}$	$28.100\times {10}^{-3}$	$-83.431\times {10}^{-4}$	$-9.9256\times {10}^{-3}$	0.1624	0.9086
$P^{*}$	$-50.952\times {10}^{-5}$	$185.50\times {10}^{-3}$	$0.05380\times {10}^{-4}$	$5.300\times {10}^{-3}$	0.2895	0.9236
$Q^{*}$	$-6.2042\times {10}^{-5}$	$9.0889\times {10}^{-3}$	$-4.0964\times {10}^{-4}$	$7.8267\times {10}^{-3}$	−0.1472	0.6904

. Finally, by using the superposition principle, the total deflection of a tooth pair j contacting at point i can be defined as

$${\delta }_j^i={\left({\delta }_b^i\right)}_p+{\left({\delta }_f^i\right)}_p+{\left({\delta }_b^i\right)}_g+{\left({\delta }_f^i\right)}_g$$

(11)

where the subscript p indicates the deformation of the driving pinion, while g of the driven gear. Those deformations are computed for each point of the profiles and will later be applied to the 3D flanks in the procedure to obtain the equilibrium contact point. Hence, the total stiffness of the engaging teeth pair j contacting in point i can be expressed as

$$k_j^i=\frac{F_j}{{\delta }_j^i}$$

(12)

2.3. Nonlinear Algorithm

As detailed in [50] during engagement the elastic deformation of the meshing teeth pairs causes a relative sliding between the contacting flanks causing a subsequent shift of the contact point where the load should be applied. Since the contact point changes, the stiffness of the engaged pair changes thus also altering the load sharing characteristics. Furthermore, due to the deflection of the gear body teeth pairs not originally in rigid contact can touch, while for the same reason the opposite might happen. For those considerations an iterative approach has been applied starting from the rigid contact conditions and then updating the contact point, the load of each teeth pair and the number of pairs in contact. At the $k^{th}$ iteration using the updated contact point for the $j^{th}$ pair the load sharing coefficient is computed by [51]

$$C_{k,j}=\frac{k_j}{{\sum }_{i=1}^N{k}_i}\left(\frac{1+{\sum }_{i=1}^N{k}_l{\tilde{E}}_{jl}}{F}\right)$$

(13)

where ${\tilde{E}}_{jl}={\delta }_j-{\delta }_l$ and $k_j$ is defined in Equation (12) while $F=T/r_b={\sum }_{j=1}^N{F}_j$, where T is the total torque to be transmitted and $r_b$ is the base radius of the pinion. A natural equilibrium condition is reached when the contact points of the different engaging pairs are in a stable position, as well as the load sharing coefficients, meaning

$$\frac{x_{k,j}-x_{k-1,j}}{x_{k,j}}<{ϵ}_x\wedge \frac{y_{k,j}-y_{k-1,j}}{y_{k,j}}<{ϵ}_y\wedge \frac{z_{k,j}-z_{k-1,j}}{z_{k,j}}<{ϵ}_z\wedge \frac{C_{k,j}-C_{k-1,j}}{C_{k,j}}<{ϵ}_C$$

(14)

where $x_{k,j}$, $y_{k,j}$, $z_{k,j}$ are the coordinates of the contact point of the $j^{th}$ engaging pair at the $k^{th}$ iteration and ${ϵ}_x$, ${ϵ}_y$, ${ϵ}_z$, and ${ϵ}_C$ are tolerance values generally equal to ${10}^{-2}\%$. Once equilibrium is reached a detailed contact model, described in the next paragraph, is used to study the contact between the so obtained deformed profiles.

2.4. Non-Hertzian Contact Model

To account for the continuously changing curvature of the flanks, the effects of the profile modifications and discontinuities, such as the gear edges and tip, a detailed numerical rough frictionless non-Hertzian contact model is implemented. The contact problem is usually stated as the Hertz–Signorini–Moreau problem [52,53,54]

$$\mathbf{h}\ge 0,{\mathbf{p}}_{\mathbf{n}}\ge 0,\mathbf{h}·{\mathbf{p}}_{\mathbf{n}}=0,$$

(15)

The first condition enforces that no interpenetration can occur between the contact bodies and, therefore, the gap function $\mathbf{h}$, which measures the distances between the surfaces, can only be positive or equal to 0 in the contact area. The second condition imposes that the contact is non-adhesive whereas the third condition enforces that the normal pressures $p_n$ can only be different from 0 inside the contact area where $\mathbf{h}=0$ and null everywhere else. The gap function $\mathbf{h}$ is expressed as

$$\mathbf{h}=h_0+\mathbf{g}+\delta$$

(16)

where $h_0$ is the indentation between the profiles imposed as a rigid body motion along the contact normal, $\mathbf{g}$ is the initial separation of the contacting surfaces and represents its topography, while $\delta$ represents the elastic deformation of the surfaces due to the applied normal pressure ${\mathbf{p}}_{\mathbf{n}}$ and can be expressed as [55]

$$\delta =\mathbf{C}·{\mathbf{p}}_{\mathbf{n}}$$

(17)

where $\mathbf{C}$ is a matrix of the influence coefficients which introduces the elasticity of the contacting surfaces. Its components $C_{ij}(i,j=0,1,\dots N)$ relate the displacement ${\delta }_i$ at a point i due to the application of a unit pressure at point j were proposed by Kalker and Van Randen [56] and later fully derived and corrected by Boedo [57]. A linearly varying normal pressure $p(\xi ,\eta )$ is imposed on the half-space region as

$$p\left(\xi ,\eta \right)=\frac{P_0\eta }{\gamma }$$

(18)

In which $P_0$ is the pressure in the apex $(0,\gamma )$ of the triangle in the local coordinate frame $(\xi ,\eta )$ with $\gamma >0$ as visible in Figure 5.

The displacement Equation (17) is then solved, and a closed-form solution is obtained

$$\begin{array}{cc}I(a,b,c,t)=\hfill & \frac{(a+bc)(t-c){[{(a+bt)}^2+{(t-c)}^2]}^{1/2}}{4\pi {(1+b)}^2\left|(t-c)\right|}+\hfill \\ \hfill & \frac{t^2-c^2}{4\pi }arcsinh\left(\frac{a+bt}{t-c}\right)+\frac{(a+bc)(2c-ab+b^{2}c)}{4\pi {(1+b^2)}^{3/2}}ln(A+B)\hfill \end{array}$$

(19)

where

$$A=\frac{2b(a+bt)+2(t-c)}{{(1+b^2)}^{1/2}}$$

(20)

$$B=\frac{2(t-c){[{(a+bt)}^2+{(t-c)}^2]}^{1/2}}{|t-c|}$$

(21)

While parameters a,b,c and t are dependent on the shape and size of the triangle, as described in [57]. Equation (19) is used in several limits in order obtain the displacements $\delta$ depending on the non-dimensional distance $\overline{y}=y/\gamma$, since it is undefined when $t=c$. The dimensional pressure-displacement coefficients are obtained as

$$w_{ij,k}=\delta \frac{1-\nu }{G}$$

(22)

where $G=\frac{E}{2(1+\nu )}$ is the shear modulus of the material, E is the Young’s modulus and $\nu$ is the Poisson coefficient and $k=1,2,3$ since the values of $w_{ij,k}$ for the entire triangulation need to be obtained three times, considering as the apex $\gamma$ a different vertex of each triangle. Once this is done, the influence coefficient relating node i to node j is found by

$$C_{ij}=\sum _{j=1}^{j=n}\sum _{k=1,2,3}{w}_{ij,k}$$

(23)

where n is the number of nodes in the triangulation and the summation is carried out by choosing the appropriate k in which node i is the vertex $\gamma$ in the local coordinate frame. The process is repeated for the pinion p and the driven gear g and the final influence coefficient matrix is defined as

$$\mathbf{C}={\mathbf{C}}_{\mathbf{p}}+{\mathbf{C}}_{\mathbf{g}}$$

(24)

The surfaces displacements can then be obtained from

$${\mathbf{h}}_{\mathbf{s}}={\mathbf{C}}_{\mathbf{s}}\mathbf{p}$$

(25)

where $s=p,g$. To solve the problem stated in Equation (15) and satisfy all the conditions a two nested sub-iterative processes are needed, as detailed in [50,58] to remove the residual traction stresses and to obtain the correct total transmitted force based on the computed force at the $k^{th}$ iteration $f_k$ and the value at the previous one $f_{k-1}$. The first one is used to identify the correct contact footprint devoid of tension forces given a certain initial indentation $h_0$, while the second is used to estimate the correct indentation to impose to the profiles so that the resulting transmitted force is within tolerance of the applied load. To validate the proposed approach and ellipsoid to ellipsoid contact is compared to the results from Hertz theory [52]. The two contacting ellipsoids ($E_1=E_2=210$ GPa, ${\nu }_1={\nu }_2=0.3$) have major radii $R_1=100$ mm and $R_2=20$ mm and minor radii $r_1=40$ mm and $r_2=5$ mm, respectively, pressed together by a normal load F = 10,500 N on a $3\times 3$ mm contact plane with equilateral triangles with a side length $l=0.06$ mm. The obtained pressure distribution is visible in Figure 6 in which the white line represents the predicted contact area from Hertz’s theory. The visual comparison shows a very good agreement and a numerical comparison is shown in Figure 7 where the obtained values for a varying load between 1000 ÷ 15,000 N show an error for both the maximum pressure and the contact area always lower than $1.6\%$. The influence of the mesh size is also visible in the same figure for the maximum load and shows a trend quickly approaching lower errors as the length of the triangles sides decreases.

The proposed contact model has been derived from the elastic half-space theory and, hence, implies that in any transverse section a plane state of deformations is respected. However, when one or both the contacting bodies have finite length, it is evident that the end faces are to be treated as free boundaries, but in a plane state of deformations two shear stresses and a normal stress would be present at the free faces. To remove those unrealistic stresses a correction based on mirror pressure distribution and an iterative computation was proposed by Heteńyi [59,60]. However, the iterative part is time consuming as the number of iterations could theoretically go to infinity as the accuracy required increases [61] and for this reason a correction factor was proposed [62] of the form

$${\psi }_g=1.29-\frac{1}{1-\nu }(0.08-0.5\nu )$$

(26)

For a body a of finite length with both free ends, the influence coefficient matrix with the corrected mirroring planes is obtained as

$${\mathbf{C}}_{\mathbf{corr},\mathbf{a}}={\mathbf{C}}_{\mathbf{a}}+{\psi }_g{\mathbf{C}}_{\mathbf{mirror},\mathbf{left}}+{\psi }_g{\mathbf{C}}_{\mathbf{mirror},\mathbf{right}}$$

(27)

If the two bodies in contact have different finite length the mirroring operation is carried out accordingly at the free edges of each body. In order to verify the non-Hertzian results and the efficiency of the mirroring correction a model of a crowned roller bearing contacting with a cylindrical inner race has been compared with the results from [61] and the geometrical and material data are reported in

Table 2. Crowned roller and inner race geometrical and material properties, from [61].

Crowned Roller
Diameter	$d_1$	15 mm
Crowning radius	$r_1$	1114 mm
Fillet radius	$f_1$	$1.006$ mm
Width	$b_1$	16 mm
Shear modulus	$G_1$	$79.3$ GPa
Poisson ratio	${\nu }_1$	$0.3$
Inner race
Diameter	$d_2$	$58.5$ mm
Width	$b_2$	25 mm
Shear modulus	$G_2$	$79.3$ GPa
Poisson ratio	${\nu }_2$	$0.3$

. A mesh of right triangles with equal sides of $0.02$ mm was defined and the resulting non-Hertzian pressure distribution obtained by the application of the current method can be seen in Figure 8 for a load F = 33,800 N where a distinct pressure peak can be distinguished at the edge of the crowned roller, while the central part has a trend similar to the ellipsoid to ellipsoid contact seen in Figure 6. The obtained results agree extremely well with the numerical and FE values detailed in [61], although minor differences are present due to the different discretization of the contact plane.

2.5. Application to Gear Contact

Since the equilibrium contact point is known through the algorithm detailed in Section 2.3, the mean tangent plane to both gear profiles is taken as the plane where contact will lie. From this line the initial separation $h_0$ is obtained as the normal distance between the mean common tangent at the equilibrium contact point and the deformed profiles. In the ellipsoid–ellipsoid contact, as well as the crowned roller and inner race contact shown earlier, the imposed rigid body indentation $h_0$ was intended as a vertical displacement of either body towards the other. In the pinion–gear contact instead, in order to respect the meshing kinematics, a rotation is imposed as a rigid body rotation of the pinion towards the gear. Therefore, at each iteration it is needed to estimate again the initial separation ${\mathbf{g}}_{\mathbf{k}}$ obtained through a tentative rigid body rotation ${\theta }_{0,k}$ for the $k^{th}$ iteration computed as

$${\theta }_{0,k}={\theta }_{0,k-1}+\frac{{\theta }_{0,k-1}-{\theta }_{0,k-2}}{f_k-f_{k-1}}(F_j-f_{k-1})$$

(28)

Evidently also gears have finite face-width and, therefore, the mirroring operation is carried out as described earlier its effect is visible in Figure 9. The data of the tested gears are listed in

Table 3. Gear pair parameters.

Parameter	Pinion p	Gear g
Number of teeth z [-]	28	28
Module m [mm]	3.175	3.175
Pressure angle ${\alpha }_p$ [${}^{\circ }$]	20	20
Facewidth b [mm]	6.35	6.35
Hub radius [mm]	20	20
Young modulus E [MPa]	210,000	210,000
Poisson coefficient $\nu$ [-]	0.3	0.3

and it is clear that if no mirroring is applied two important pressure peaks are present at the free edges of the contacting flanks, which disappear as expected when the proper measures are taken. However, if one of the two bodies has a face-width greater than the other those peaks correctly reappear, although with lower values. Furthermore, gears have another free side at the tip of the tooth which must be taken into account. For this reason, a tip mirror plane is accounted for when part of the contact plane of one of the two gears lies outside of the flank and in this case the tip edge is taken as the mirror line, as visible in Figure 10.

3. Application to Case Studies

The final expression for the STE, considering the cumulative effect of the SA deflections and the contact analysis, can be given as

$$STE=max\left[{\left({\theta }_g{r}_{b,g}-{\theta }_p{r}_{b,p}\right)}_j\right]$$

(29)

where the rotations of the $j^{th}$ engaged pair of the pinion (p) and driven gear (g) are measured at the point of maximum displacement inside the contact area. Since the analyses in the present paper will deal with quasi-static conditions no dynamic effect is included. Additionally, no manufacturing, mounting, or other possible errors are considered at this point, hence the STE is equal to the global transmission error (TE). Results from the proposed SA approach are firstly compared against a 3D FE model from Ansys which adopts quadratic 20-node solid elements (SOLID186) and a very refined mesh ($0.01$ mm for a total of 441,562 nodes, 394,002 elements and 1,311,456 degrees of freedom) for the contacting flanks which are modeled as pairs of frictionless contact-target elements (CONTA174-TARGE170) in surface-to-surface contact using the augmented Lagrangian algorithm. The hub of the pinion gear is grounded through rigid body elements to a central node which is only free to rotate around its axis and where the torque is applied and the displacement recorded. Similarly, the hub of driven gear is constrained against all possible motion. The FE mesh and an example of the obtained contact pressures are visible in Figure 11, while the relevant data are listed in Table 3. A sequence of 20 nonlinear static analyses with automatic substepping is performed at different rotation angles to cover the whole mesh cycle which corresponds to $\frac{2\pi }{z}$, as done in [63,64]. Since the transmission ratio of the studied gears is equal to one, between each analysis the FE mesh of the gears is rigidly rotated of the same angle according to their rotation directions, hence the motion is not caused by the applied torque but imposed by the angular displacement of the contacting gears. A comparison of the different approaches is visible in Figure 12 for three different torque levels. The overall agreement of the results is evident, while a small difference in the double contact zone is visible which amounts to a maximum of $0.{002}^{\circ }$ at the highest load. The agreement is extremely good also in terms of maximum contact pressure, as visible in Figure 13, for both the single and double contact regions, as well as the pressure peaks occurring when tip-corner contact is present. However, it must be noted the the calculation in Ansys on a laptop equipped with an Intel Core i7-7700HQ CPU and 32Gb ram takes around 4 h with a so refined mesh at the surface of the engaging teeth, whereas the proposed SA approach employs only 3 min on the same machine, yielding a significant time saving without even considering the time spent to properly generate the FE model.

The STE and load sharing coefficients of the same gear pair under a load of 100 Nm are visible in Figure 14. In the same figure, it is possible to appreciate the dependency of the load sharing factor and contact ratio with respect to the torque level. In particular, it is evident that the increase in the contact ratio as the torque increases which is not considered by the standards. The actual contact ratio, which can be extracted by the load sharing characteristics, ranges from the value computed according to literature [65] 1.638 to the highest 1.85 due to the maximum torque, that corresponds to an increase of about 13%. A detailed pressure map of the flank at the highest torque displaying the pressure peaks due to tip contact is shown in Figure 15. Under these conditions the maximum pressure reaches 2648 MPa which could cause surface damage since it is higher than the design value in the single contact zone (1336 MPa).

To mitigate this problem, the effect of the linear tip relief is now analyzed. Firstly, an amount of material removed $Ca$ is analyzed for a fixed modification diameter $dCa=94.245$ mm as visible in Figure 16 varying from $Ca=0.008÷0.040$ mm. Evidently, this modification slightly affects the mesh-in and mesh-out condition of the engaging flanks, effectively reducing the contact ratio since the single contact zone of the mesh cycle grows. Furthermore, the load sharing characteristics are not greatly affected since after a certain threshold increasing the amount of material removed effectively means reducing the working part of flanks due to a lower $dNa$ [66], thus the curves of the STE and load sharing factor (LSF) converge to a new stable condition. Instead, modifying the length of the modification by reducing the start diameter of the modification $dCa$ has a huge impact, as visible in Figure 17. The start of the modification is changed in the range $dCa=94.245÷90.245$ mm and evident improvements are visible on both the STE and LSF. Indeed, the minimum value of the STE increases as the length of the modification increases, effectively reducing the peak to peak transmission error (PPTE), while simultaneously smoothing the transition between the single and double contact reducing the possibility of impacts on the flanks due to abrupt discontinuities. This kind of profile modification also has beneficial effects on the pressure distribution on the flanks as visible in Figure 18, where a linear tip relief was applied with $Ca=0.032$ mm and $dCa=94.245$ mm. The contact area is reduced since the flanks are not able to come into contact up to the tip circle, but are limited by the presence of the modification. However, the pressure peak is greatly reduced, displaying a value of 1780 MPa which, while higher than the corresponding value in the single contact zone, is significantly reduced compared to the case where no TPM was applied.

Similar considerations can be applied when the amount of material removed $Ca$ is applied on a parabolic tip relief, as visible in Figure 19. Due to the fact that the removal rate is lower than the linear one, the effects of this parameter for a fixed $dCa=94.245$ mm are even lower than before and can hardly benefit the transmission both in terms of STE and LSF. Indeed, the different curves are almost overlapping, showing only a minimal increase in the single contact region as expected. Additionally, increasing the length of modification has similar effects to what has been shown before as visible in Figure 20. Indeed, the PPTE decreases as the length increases and the transition of the load between the teeth pair is smoothed, but the overall effect is lower than the linear modification due to the lower removal rate of the material earlier mentioned. However, a gear pair modified with parabolic tip relief $dCa=94.245$ mm and $Ca=0.032$ mm has significant benefits in terms of pressure distribution along the flank as visible in Figure 21. This kind of modification is even more effective at reducing the pressure peak around the tip and root regions since the transition between the pure involute and the modified portion of the tooth is smoother and the lower removal rate offers more material for the flanks to come into contact with. The maximum pressure achieved is indeed only 1426 MPa which is just 90 MPa higher than the single contact zone thus reducing any risk of potential damage.

Lastly, a combination of parabolic tip relief and face-width crowning will be discussed. The STE and LSF have been computed for different levels of symmetric crowning with $C_{\beta }=C_{\beta In}=C_{\beta Fin}$ in combination with a parabolic tip relief with $dCa=94.245$ mm and $Ca=32\mathsf{\mu }$m and the results are visible in Figure 22. The axis of the parabola is fixed at $zL=b/2$ (Figure 1), while the amount of material removed is evaluated at $zLc=0.8b$ as commonly done in the industry. Since the symmetric crowning modification does not influence the mesh-in or mesh-out phases the LSF is unaffected while also the STE has minimal changes. Indeed, the only difference is the contribution of the elastic deflections ${\mathit{\delta }}_{\mathit{p}}$ and ${\mathit{\delta }}_{\mathit{g}}$ of the contact surface that increase since the contact are is more and more localized towards the axis of the parabola. This localization is also evident in Figure 23 where the pressures experienced by the flank are visible. Due to this effect the side edges are completely unloaded, however the pressure the gears are subjected towards b/2 are increased up to 1539 MPa in the single contact region and up to 1762 MPa in the tip and root zones. This kind of modification is indeed useful mostly when gears with different face-width need to be designed or when shafts misalignment due to deflections must be compensated, since otherwise side pressure peaks would appear as discussed earlier in Figure 9.

4. Conclusions

Aim of this paper has been to establish an accurate model to determine gear deflections under load and the contact condition that are present during meshing. For this purpose, on the basis of well established SA foundations a nonlinear iterative scheme was implemented seeking for a natural equilibrium condition between the location of the contact point, which slides due to the deflections, the load intensity, since the change in contact position alters the stiffness as well, and, ultimately, the number of the engaged teeth pairs. The contact between the deformed flanks has been studied in detail using a non-Hertzian contact model to account for the variable curvatures of the flanks, as well as the profile modifications. Appropriate measures have been included to avoid overestimations of the pressures at the free edges of the finite length gears by including a mirroring of the contact plane. A mirroring operation has also been carried out when one of the tooth tips comes in contact with the mating flank of the opposing gear. The contact model has been firstly validated against Hertz hypotheses and later for a non-Hertzian contact, giving extremely accurate results. Then, the whole proposed model has been applied to a test case and the comparison of the results against a FE counterpart has shown very good agreement. The same test case has been studied under a variety of combinations of profile modifications, namely linear and parabolic tip relief, as well as crowning. For all those tests the STE and LSF have been analyzed along with the resulting pressure distributions on the flanks of the teeth, showing different behaviors. The linear tip relief has been found to be able to effectively modify the STE and reduce the PPTE, but the parabolic one was more effective at reducing the pressure increase when contact is close to the tip of one of the gears. Indeed, the proposed SA model produces results with the same level of accuracy of a refined FE model, with some significant advantages. The setup of the model for the proposed approach is almost automatic by just specifying the parameters of the generating tool and of the wanted gear pair, while complex operations need to be carried out to set up the FE counterpart correctly. Additionally, the saving on the computational time is huge since the proposed method generates the results in minutes, while several hours are required if the FE approach is chosen instead. The present results also highlight the strong influence between the contact ratio and the applied load, the micro-geometrical modifications and the STE. A correct balance of these aspects can improve not only the endurance of the studied gears, but also their dynamic behavior. Further research will be carried out to further increase the prediction capability of this approach firstly by extending the method to model also helical gears by introducing the appropriate modifications where needed. Next, the influence of other components of the transmission such as shafts, bearings and the casing will be included to correctly estimate the contact under misaligned conditions.