Stress calculation on bevel gears with FEM influence vectors

In order to be able to carry out an optimal gear design with the aim of cost reduction and the careful handling of resources, load capacity is an important criterion for the evaluation of a gear. For the calculation of the flank and root load capacity, a precise loaded tooth contact analysis (LTCA) is necessary. With LTCA software like BECAL, influence numbers are used to calculate the deformation of the gear. These influence numbers are calculated with a BEM-module and considered for calculating the local root stress. This method simplifies the coupling stiffness in tooth width direction with a decay function and neglects the influence of local differences in tooth stiffness. In this publication, this simplification shall be questioned and evaluated. Therefore, a new method for calculating stress with FEM influence vectors is presented. This method enables the calculation of full stress tensors at any desired location in the gear with the efficiency of the influence number method. Additionally, the influence of local stiffness variations in the gear is taken into account. Various gear examples show the influence of material connections at the pinion root and the influence of the rim thickness of a wheel on the root stress. To validate the accuracy and the time efficiency of the new calculation method and to compare the results to current state-of-the-art simulations, a well-documented series of tests from the literature is recalculated and evaluated.


Introduction
The calculation of load distribution with bevel gears takes many hours of modeling and calculation time with conventional FEM solvers. Therefore, common design tools like BECAL [1] use the method of influence numbers and achieve a calculation time of only a few minutes. Here, a load distribution is calculated on the basis of a stiffness matrix and the local contact separation between the tooth flanks. The accuracy of the load distribution calculation depends directly on the quality of the stiffness matrix. The determination of tooth stiffness is already well advanced and executed with a 2D BEM approach. That provides fast and sufficiently accurate stiffness values with the aid of approximated, generalized influence functions and a twodimensional, numerically determined reference value at the force application point. The linear deformation components due to bending, shear and compression are taken into account. Added to this is the non-linear contact deformation, which is iterated as a function of force. The compliance of the wheel body is approximated by the connection of the tooth geometry to an elastic half-space. In reference [2] this simplification is questioned. The authors show that the axial deflection of the wheel body and furthermore the twisting of a thin pinion has to be represented with separate FE influence numbers. The described method provides the possibility to calculate the entire gear deformation with the FE method. This provides the possibility of calculating the stress at any desired point on the tooth surface and in the tooth volume. In Fig. 1 the local root stress can be evaluated on the example of a bevel gear. The necessary steps to enable this calculation will now be described in the following. Fig. 3 Load distribution equation system [2] Fig. 4 Comparison of BECAL and ANSYS contact analysis deformation angle φV (left) and maximum contact pressure (right)

Meshing of bevel gears
In order to calculate influence numbers with the FEM, the bevel gear geometry from the manufacturing simulation must be converted into a three-dimensional FE mesh. For this purpose, the tooth flank is meshed over its surface with a combination of rectangular and triangular elements while maintaining a good aspect ratio. The meshing is then carried out in the tooth thickness direction. Similar to the cylindrical gear meshing, here either the elements of the left and right tooth flanks are merged in the tooth center or transferred in the direction of the gear body connection. To do this, the transition point with the smallest distance to the gear body connection point in the center of the tooth has to be identified. Due to the different face geometry at the heel and toe, this calculation is performed separately at the face of the toe and heel for the left and right flank, respectively, and is subsequently averaged. In order to increase the mesh density in the tooth root and near the flank, a variable gradient is inserted in the tooth subsurface direction. The meshing strategy enables consistent aspect ratio for the entire tooth flank (Fig. 2).

Calculation of load distribution
The BECAL load distribution is calculated with local deformation influence numbers eik.
In [3] the components of the load distribution equation system are described as the following. The influence numbers are representing the deformation at the location i as result of the load at location k. The resulting deviation of the contact line position f i and the resulting pinion torque TRi are considered as inhomogeneous links. The deformation leads to a small torsion φV around the pinion axis, with the sections i of the flanks (in normal direction) approaching around φV rwi. In [4] the calculation of the deformation influence factor is extended to all teeth in contact, all entries of the influence factor matrix can be filled. The equation system for calculating the load distribution has the following appearance (Fig. 3): The findings in [4] allow the calculation of the load distribution for an elastic gear using the FEM influence numbers. In order to evaluate performance and accuracy Fig. 4 presents a comparison of the influence number method with a detailed FEM contact simulation (ANSYS). For this purpose, the simulated deformation angle φV and the pressure distribution for 49 contact positions are evaluated. Compared to the influence number method, the contact simulation is a higher class simulation model that includes the moving contact line under load and the curvature in face width direction. To achieve good pressure results with a commercial FEM contact solver a very fine mesh density is needed. Therefore, a convergence study was performed with the criterion of 1% pressure deviation between the last refinements. The ANSYS contact problem results in 3.8 million nodes and 3.5 million hexahedron elements. For calculating 49 mesh positions the ANSYS contact calculation takes 48 h on a state of the art workstation. The BECAL calculation with FEM influence numbers only needs up to 15 min (FEM EFZ). A very good correlation of the simulation models becomes clear. Considering the two FE calculation methods (FEM EFZ and FEM contact), the overall deformation behavior (' V < 1.6%) and the max. pressure The extremely time-reducing model simplifications carried out in BECAL do not pose a problem in the calculation of the load distribution on this bevel gear. In the next step, the influence numbers can now be extended to other tooth areas in order to determine the local tooth deformation.

Root stress calculation
The FEM influence numbers offer the possibility of calculating the deformation at any location in the tooth, from which the local stresses can then be derived. In contrast to the deformation calculation for the load distribution, the volume deformation must be calculated as a three-dimensional deformation vector. Therefore, it is necessary to switch from the calculation of influence numbers to the calculation of influence vectors (Fig. 5). These vectors represent the local deformation during load application at the tooth flank. While influence numbers always point in normal direction, influence vectors can be oriented arbitrarily.
This method can now be used to calculate the tooth root deformation by superposing the deformation components with the face load of the respective meshing position. For this purpose, the local tooth forces E f k are multiplied by the influence matrix e using the line load and the width of contact deformation.
The tooth deformation E u k is now determined over all mesh positions and written down in addition to the position vector E x k in the element. This forms the basis for local strain and stress calculation at any location in the calculation model. Fig. 6 shows an example of the deformation of a gear due to an applied line load.
With the help of the description in the commercial FE solver [5] and the explanations in the literature [6], the calculation path for the local stress calculation on the hexahedron is explained in the following. First, the shape function is used to calculate the deformation in the element in the image domain. For a hexahedron the first form function is for example: The shape function can then be used to determine the fraction of rigid body displacement and the deformation separately at the integration points at r; s; t =˙1= p 3: x .r; s; t/ y .r; s; t/ z .r; s; t/ By setting up the Jacobian matrix J i .r; s; t/ and D i .r; s; t/, the strain " i at the integration point i can be determined: @x @r @y @r @z @r @x @s @y @s @z @s @x @t @y @t @z @t 3 7 5 I D i .r; s; t/ = 2 4 @u @r @v @r @w @r @u @s @v @s @w @s @u @t @v @t @w @t " xx " xy " xz " yx " yy " yz " zx " zy " zz Subsequently, the extrapolation from the integration point i to the nodes k is performed: The element stress k at the node can be calculated from the element strains using the material-dependent stiffness tensor C , where is the double contraction [7]: By averaging the stress tensors of the adjacent elements, the nodal stress tensors can be calculated from the element stresses. The result of the calculation is the stress tensor for each node in the tooth root. From this the principal stress tensor H can be determined: By saving the maximum values of the first principal stress σI and minimum values third principal stress σIII over the entire tooth contact, the cumulative maximum and minimum root stress can now be determined analogously to the cumulative pressure distribution. Both are presented in Fig. 1. The result can now be used to calculate the local root load capacity according to ISO 10300 [13]. Since the stress result, based on the FEM, depends on the mesh discretization, a convergence analysis must be performed for each calculation run. Fig. 7 (right) shows the result of a convergence analysis on an example gear pair. To simplify the convergence analysis, a global parameter for the mesh density is used. The global mesh parameter modifies the resolution in tooth width, height and subsurface direction. In order to compare the root stress calculation with a commercial FEM solver, the identical FE problem is calculated in ANSYS [6] and compared with the BECAL results. The comparison gives a good correlation with an overall maximal deviation less than 2.5% (Fig. 7 left). The calculation in BECAL show higher stresses in terms of magnitude and is for the studied example mathematically on the safe side.

Validation
Now the two methods for calculating the root load capacity are available and can be compared. As a comparison of the methods, a validation with the well-documented experimental results for the root load capacity from Wirth [8] will be carried out in the following section. The Wöhler curves are determined for the Tooth Root test bevel gear set F0 and hypoid gear sets F15 (15 mm Hypoid offset) and F31-75 (31.75 mm Hypoid offset) in a comprehensive test  program with a calculated 1% failure probability. The gear pairs with a mean normal module between nmn = 2.2 mm and 2.46 mm and an average spiral angle of between βm = 33°a nd 44°will not be presented further here, but reference will be made to the explanations by Wirth [8]. With the aid of the respective gear pair drawings and the assembly drawing, it was possible to model the wheel body with the available wheel body variants from BECAL. For the pinion, a gear body variant with adjacent shaft shoulder with heel constraint was selected, and for the wheel a wheel body with hub and heel constraint. The constraints are visualised in the following graphics with red areas. Fig. 8 shows the maximal tooth root stress in tooth width direction and the overall tooth root safety factor for different calculation variants. Fig. 8 (upper section) shows the comparison between the BECAL calculation method BEM and FEM using the example of the hypoid pinion F31-75. A significantly lower maximum stress with an deviating location is shown for the root stress calculation with FEM. Looking at the tooth root stress curves in the tooth width direction, Fig. 8 (lower section), it is apparent that the pinion body and its clamping at the heel also have an influence on the root stress curve and the root stress maximum. The realistic connection to the pinion shaft with a shaft shoulder (bottom section left) or a transition into the shaft (bottom section center) should therefore be aimed for. A constraint near the tooth root is not recommended due to the additional stiffening effect (bottom section right). Fig. 9 (upper section) shows the comparison between the BECAL calculation method BEM and FEM using the K Fig. 10 Recalculation of the test results [8] with the target SF1 = 1 example of the hypoid wheel F31-75. An comparable stress is shown for the root stress calculation with FEM. Looking at the tooth root stress curves in the tooth width direction in Fig. 9 (lower section), it is apparent that the wheel body and its clamping at the toe have an small influence on the root stress curve and the root stress maximum. The higher rim thickness (B = 16 mm / bottom section left) increases the wheel stiffness and provides slightly lower root stress compared to the smaller wheel thickness (B = 23.5 mm / bottom section right).
Using the experimentally determined pinion torques an exact load capacity calculation must have the result SF1 = 1. Values above SF1 = 1 are on the uncertain side and values below SF1 = 1 are on the conservative but safe side. Therefore, a comparative calculation can now be performed. The results in Fig. 10 consistently show a conservative simulation result for BEM and FEM. The FEM is consistently closer to the experimental results. In addition to the failure-relevant pinion load capacity, the simulated load capacity of the wheel is plotted on the right in Fig. 10. At the wheels, unexpected tooth root fractures occurred in the project despite conservative design. Therefore, to achieve pinion failure, the wheel root had to be shot peened. The comparison of the simulation results of pinion and wheel load capacity using the FEM shows a similar load capacity of both gears here. Therefore, the results show better correlation with the wheel failure observed in the test.

Volume capacity
In the field of volume capacity, there are many simulations and experimental studies, e.g. Wirth [8], Bauer [10] and others [14][15][16][17][18][19][20][21][22]. Currently, there is no subsurface flank fracture ISO Standard for bevel gears available. The knowledge gained from cylindrical gears ISO 6336-4 TS [9] has not yet been transferred to bevel gears. The design of test gears proved to be difficult [10]. In Bauer's series of tests, only two of eight gear pairs failed with flank failure. In these cases, the failure occurred in combination with pitting and the experiment was not fully successful. The results from Bauer [10] shall now be simulated to investigate the rea-sons that failures occurred only sporadically and as mixed failures. The presented stress calculation method is now applied to the tooth volume. Here, the stress tensor profile under the tooth contact for all mesh positions is the center of interest. This time-dependent local stress calculation provides the basis for various verifications, in the following the Von Mises equivalent stress is used as an example. More complex time dependent stress calculation methods such as shear stress intensity hypothesis SIH [12] are also possible. For his flank fracture load capacity Bauer [10] uses an approach using the equivalent stress amplitude according to Von Mises σVA(MISES) and calculates the maximal amplitude in subsurface direction y for all contact positions x: This can now be applied on the local calculation model and describes the maximal equivalent stress amplitude for all contacting teeth k, for each contact position j and every location in the tooth volume i: This equivalent stress amplitude is calculated from the local stress tensor. In Fig. 11 the maximum value of the equivalent stress over all contacting teeth and mesh positions is visualized with the maxima at one tooth. The tooth is cut in the middle, near the heel and the toe to visualize the subsurface stress.
The quality of the stress profile shall now be investigated. For this purpose, an FEM contact analysis was carried out and the stress profile was evaluated as well. The results are presented in Fig. 12. The depth profile shows good correlation. However, a considerable deviation at the surface is noticeable, which can be attributed to the load introduction at the influence vector model. The load applied is not as uniform as determined by the ANSYS contact solver. Therefore, when using the influence vector method the stress evaluation at the surface is currently not recommended. However, the relevant failure mode on the flank surface is pitting load capacity. Stress evaluation is therefore not necessary in this area. The deviation at the maximum stress amplitude is lower than 3%. The contact calculation was performed on the converged model for the tooth root stress at 49 mesh positions. The BECAL calculation time was 15 min, the calculation time of the FEM contact model was 4.5 h. The second side of the volume capacity approach is the calculation of the local fatigue strength. In the following, the fatigue strength is calculated with literature approaches from the local hardness [11]. The local hardness can be determined with the aid of a complex and time-consuming FEM heat treatment simulation or alternatively with a simplified calculation approach. Real gear sets have significant differences in hardness between tooth tip, flank and tooth root zone. Furthermore, there are differences of the core hardness within the tooth between the different zones. For this first simple approach no residual stress is considered, the hardness HHV(y) is calculated using the distance to the surface y according to Lang [11] and depends on the case hardening depth (CHD) and the heat treatment parameters (a,b).
H HV .y/ = H HV;core +.H HV;surface −H HV;core / f .y / (12) f y = 10 .a+b y / y (13) This approach is applied at all nodes of the three-dimensional tooth, therefore all following parameters are local and referred to echo node (i). The result is the hardness HHV(i) calculated with the distance to all free surfaces (y) (Fig. 13). This is a uniform approach. For the gears used in Bauer's experiment, the case hardening depth was determined with CHD = 0.3 mm .H HV;core = 450H V; H HV;surface = 700H V /.
The alternating strength σW(i) is calculated from the hardness HHV(i). This approach is related to Bauer [10] but neglects the proposed material correction factor, as it is not statistically validated. The local flank fracture load capacity In the experiment, for the case hardening depth of CHD = 0.3 mm the global minimum is calculated with SFB1 = 0.8. The minimum of local safety factor against flank failure is presented in Fig. 14.
Due to reusability of the influence vectors this method allows rapid variation of load and case hardening depth. To illustrate the possibilities of the simulation in Fig. 15 a simulation matrix is presented with the load levels 400/600/800 Nm and case-hardening depth of 0.3/0.6 mm and therefore the flank fracture load capacity SFB1, the pitting load capacity SH1 and the root load capacity SF1 can be analyzed. This variation offers the possibility to investigate the damage tendency of the gear. In the presented example, the pitting load capacity is achieved at a pinion torque of approx. 600 Nm. At a case hardening depth of CHD = 0.3 mm the failure mode changes from pitting to flank fracture when increasing the load from 400 to 600 Nm. At a case hardening depth of 0.6 mm, no load can be shown to have a flank fracture load capacity lower than the pitting load capacity. In Bauer's [10] experimental investigations, mixed damages of pitting and flank fracture occurred at the load level of 800 Nm and CHD of 0.3 = mm.

Conclusion
The presented method for stress calculation with FEM influence vectors allows the consideration of the tooth geometry and gear body constraint in the load distribution calculation of bevel gears. Due to reusability of the influence vectors this method allows rapid variation of input parameters and is therefore well suited in the design process of gears. It could be shown that this method allows a local evaluation of the root load capacity and shows compared to the BEM method a better correlation with experimental results for bevel and hypoid gears. In the next step, the method was extended to the tooth volume. It could be shown that the local stress tensors were calculated and, in connection with the load capacity, made accessible to a verification under the surface. The method was only used to illustrate the possibilities of the simulation at an example. In the next step, existing approaches must be incorporated, new approaches must be designed and the accuracy must be assessed on the basis of verified experimental results. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4. 0/.