Several procedures have been developed to detect the motion of bodies and apply constraints to prevent penetration. These procedures employ Perturbed or Augmented Lagrangian methods, penalty methods, and direct constraints. Contact simulation has often necessitated the use of special contact or gap elements handling the interaction between two surfaces coming into contact. However, integrating these elements and ensuring their accurate behavior can be challenging and computationally intensive. Instead, MSC MARC addresses contact interactions at the level of individual nodes and segments within the existing finite element mesh (MSC.Software, 2021). During the iteration process, A bounding box algorithm is employed to swiftly ascertain whether a node is near a segment. Once contact is detected, it is directly satisfied through multi-point constraints (MPC) that ensure compatibility of displacement fields at contact points rather than relying on contact forces or penalties. This direct constraint approach offers a more localized contact treatment, contrasting with method that uses contact elements where contact inequalities are determined at integration points, such as nodal points or Gauss points.

The constraint method, as implemented in MSC MARC, offers an efficient and direct way to handle contact problems without introducing the complexities and potential inaccuracies associated with special contact elements. The node-to-segment procedure simplifies the process further by only checking nodes against segments, reducing the computational burden and potentially increasing the accuracy of the simulation. In the augmented Lagrangian method, penalty terms and Lagrange multipliers are incorporated into the optimization process. The most fundamental difference is the Augmented Lagrangian’s reliance on penalties and iterative refinement of Lagrange multipliers, while the Constraint method enforces constraints directly without using penalties (Kloosterman, 2002). This direct constraint approach offers advantages in terms of simplicity and computational efficiency.

In this paper, the FE-based contact between the gear teeth is solved by using a node-to-segment numerical procedure. At the start of each iteration, the motion of the bodies is checked to see whether a node has crossed a face by verifying its placement within the contact tolerance. The contact tolerance size significantly impacts the solution’s computational costs and accuracy, as the driving gear initially rotates such that the nodes almost touch a surface. In this case, a biased tolerance zone with a smaller distance on the outside and a more considerable distance on the inside is considered to prevent close nodes from contacting each other and separating again.

When contact occurs, a tie is activated, relating the displacement increments of the contacting node to the displacement increments of the boundary nodes. This concept is depicted in Figure 1.

The displacement increment of the tied node during contact T in the normal direction to the plane, which contains the point of contact P, is a weighted average of the displacement increments of the retained nodes R 1 and R 2 , or: Δ u n T = α   Δ u n R 1 + 1 − α   Δ u n R 2 (4)

Where α is defined by α = PR 1 R 1 R 2

The analysis will reveal the following: Δ v n T = α   Δ v n R 1 + 1 − α   Δ v n R 2 (5) and Δ a n T = α   Δ a n R 1 + 1 − α   Δ a n R 2 (6) where Δ u n , Δ v n and Δ a n are respectively the incremental nodal displacement, velocity, and acceleration components normal to segment R 1 R 2 . n being the normal to the segment. At the end of each iteration, if the reaction force between the node and the surface becomes tensile or positive, it should separate, and the contact constraint should be removed.

However, nodes are checked for separation only if the convergence criteria are met. For nonlinear contact analysis, the Full Newton-Raphson Method is used to solve the equilibrium equations considering the following equation: K _ _ u   δ u _ = F _ − R _ u (7) where K _ _ is the tangent stiffness matrix, u is the nodal displacement vector, δ u _ is the displacement increment, F _ is the applied nodal force vector, and R _ is the internal load generated by the elemental stresses. In the proposed approach, convergence is achieved once F − R u <   ε and ε = 10 − 2 . After performing a convergence study, the convergence criterion is chosen to ensure accurate and efficient results. An iterative penetration-checking procedure is used in which the iteration process is done simultaneously to satisfy the contact constraints and global equilibrium using the Newton-Raphson method. After obtaining the solution δ u i from (7), if a node penetrates a contact surface, a scaling factor is utilized to adjust the displacement change so that the node is relocated to the contact surface. Assuming s represents the proportion of δ u i to prevent a new penetration, the displacement increment is modified accordingly: Δ U i = Δ U i − 1 + s   δ u i (8)

Furthermore, the total displacement is U n = Δ U n − 1 + Δ U i (9)

The FE mesh is generated according to the existing gear pairs geometry using solid elements 3D CHEXA Nastran elements (MSC Software Corporation, 2022). A fine mesh is defined in the areas of the tooth where contact is expected to occur as the nonlinear contact can be captured. A much larger mesh was used for the teeth far from the contact zone. The supporting shafts are modeled as infinitely rigid via a central master node rigidly connected to the inner nodes of the gears. The gear backlash is covered by rotating the driving gear to achieve contact at an initial state. Finally, the obtained FE model is imported into MSC MARC Solver for the STE Computation.

A single-stage transmission is considered, where the elasticity of the gears is taken into account, and the quasi-static behavior of the gear pair is evaluated. The meshed bodies are automatically detected, and the driven and driving wheel are identified. The material is assigned, and a new contact interaction is defined. The next step is to set the relative boundary conditions at the wheel’s central nodes 6 DOFs. The applied torque is assigned to the driving gear where the rotational DOF along the axial direction is left unconstrained. While the driven one will get the requested rotation angle α of the wheel to span, at minimum, an entire meshing cycle. Then the static transmission error is computed along the line of action. The discretization of the angular position of the driving gear must be finely defined to describe the sudden change in the STE correctly. Thus, a fixed-time stepping is used with α = 360 / z / n , where z is the number of teeth and n is the number of discrete rotational positions. For spur and helical gears, n is set to 30 and 45 for every mesh period, respectively.

The accuracy of the approach can depend on several factors, e.g., the choice of contact parameters and the mesh density. As with any numerical simulation technique, it is essential to carefully validate and verify the results simulations to ensure that they accurately capture the behavior of the modeled system. For this reason, some numerical test-cases for standard spur and helical gears are considered in the next section.

A convergence study on the FE mesh of the tooth is conducted to establish an efficient multibody simulation. Three mesh cases along the involute profile were examined, corresponding to mesh element sizes of 100, 70, and 50 µm. The goal was to achieve an accurate STE response estimation with minimal computational expense. This procedure is validated using a spur gear pair with the characteristics outlined in Table 1. The mesh generation, consisting of approximately 170,000 nodes and around half a million degrees of freedom, is illustrated in Figure 2A. A detailed view of the teeth region can be seen in Figure 2B.

For the nonlinear quasi-static analysis, an output torque of T = 50 Nm was applied. Figure 3A presents the associated STE for different tooth mesh detail levels. The shapes of the curves are nearly identical, indicating a smoother trend with a more refined mesh. All three curves exhibit comparable PPTE and shapes, reflecting similar amplitude and harmonic distributions of the periodic excitations. Notably, as the mesh size reduces from 100 µm to 70 µm to 50 µm, the simulation time over two mesh periods increases sequentially, requiring 1177 s, 2030 s, and 3031 s respectively. Consequently, the 100 µm mesh size was selected as the optimal choice for spur gear cases in terms of involute profile refinement.

For the selected case, inflations of 10 and 50 layers for the bulk and along the tooth width were compared, corresponding to 1 mm and 0.2 mm lengths. Figure 3B displays consistent results, even though the wall time notably increased to 13449s. This increase in time is attributed to the DOFs nearing 1.5 million in this instance. Therefore, the final mesh chosen for the spur gear case was set with a 100 µm mesh element size along the involute profile and a 1 mm length along the face width.

A similar procedure was applied to the helical gear pair described in Table 2 to derive an accurate yet cost-effective solution, considering the helix angle. Then a 0.5 mm element length along the face width is adopted.

In ROMAX^© (Romax Technology, 2021), the microgeometry analysis is based on a thin strip model where the contact teeth are divided into many uncoupled springs across the face width. Figure 4 gives a simplified overview of how this works.

The first tooth is modeled as a series of uncoupled springs. These springs capture the profile of the micro geometry modifications. In the example below (for simplicity), only one gear is modeled with springs. However, in the ROMAX^© software, both teeth are modeled with springs to capture the difference in pinion and wheel gear tooth stiffness.

Considering the backlash, misalignment and torque to calculate the deformed profile, the two gears are bought into contact. This is repeated over the roll angle steps. The microgeometry analysis assumes that the misalignment is constant with the roll step.

The image also shows the absolute value of the TE for one roll step (the Sum of the effects of backlash and deformation). This is calculated at each role step, and the peak-to-peak and harmonics are calculated.

On the other hand, VIBRAGEAR^© presents two different ways to calculate the flexibility matrix:

• An analytical model based on the Ritz method, used to calculate the flexibility matrices associated with the contact lines from an analytical model of a thick plate of variable thickness using a Reissner-Mindlin theory.

A spur and helical gear pair configurations were tested to develop different test conditions in order to better understand the general trend between the results. The macro-geometry parameters and the operating torque are reported respectively in Tables 1, 2 for spur and helical pairs.

Figures 5A, B demonstrate that all software yield comparable peak-to-peak transmission error (PPTE) values for the helical and spur gear cases, respectively, with results reported in Table 3.

Another comparative criterion is the shape of the curve. ROMAX^© and VIBRAGEAR^© exhibit similar behavior, with a sudden transition observed from two teeth in contact to one. In contrast, the Multibody FE approach presents a more gradual transition. This difference is attributed to commercial software not directly calculating displacement at a contact point but interpolating it from surrounding nodes, as shown in Figure 6B. Consequently, this approach fails to capture the tip contact at the transition phase, a capability inherent to the multibody method, as shown in Figure 6A. This phenomenon affects the energy distribution between the harmonics and should be considered for a more accurate dynamic response.

Figure 7 presents an example of the results obtained from the nonlinear FE method near the meshing area. Each image displays the equivalent Von Mises stress distribution for the spur gear pair. When the number of engaging teeth changes from one to two, a reduction in the maximum stress value is noticeable from 1,005 MPa in Figure 7A to 695 MPa in Figure 7B. This variation is due to the time-varying nature of the mesh stiffness. Additionally, this approach allows for a closer look at root stresses, where fluctuations in stress levels can be identified during the contact treatment phase.

Another case study is performed between pairs of spur gears with a different number of teeth, i.e., a reduction ratio R ≠ 1 . The macro-geometry parameters are outlined in Table 4. Several analyses were performed for different input torques. The results presented in Figure 8 are consistent with previous findings, as all three approaches produced similar results. Additionally, an increase in PPTE values was observed with increasing torque. However, differences between the methods in STE form were influenced by the torque level.

Similar to the earlier observations, ROMAX^© and VIBRAGEAR^© demonstrated a sudden shift of the STE in this study, with the flexible multibody method indicating a more gradual transition. As shown in Figures 8A–C, this effect becomes increasingly pronounced as torque levels rise. Such a pattern is linked to the theoretical contact lines proposed by ROMAX^© and VIBRAGEAR^©, which restrict potential contact points to no-load contact points for the corresponding conjugate gear. However, due to load-induced deflection, the tips of the teeth can initiate contact, a phenomenon referred to as “corner contact” in the literature (Langlois et al., 2016).

Figures 8D–F validate this behavior, demonstrating the contact status from the proposed method at the same rotational position across three different torque levels. Figure 8F displays two teeth in contact where the highest torque (200 N m) is applied, while Figure 8D shows a single tooth in contact when subjected to a torque of 80 N m. This increase in the operational contact ratio contributes to the changing shape of the STE curve observed with the present approach.

On the other hand, Figure 9 is showing a comparison between the methods from a peak-to-peak prospective. It is observed that PPTE results obtained with the proposed approach are located in between the ones estimated by ROMAX^© and VIBRAGEAR^©.

The latter methods present a linear increase against the applied load, whereas the multibody simulation shows a nonlinear trend. This aspect is related to the following points:

• The proposed approach presents a contact resolution based on a nonlinear contact analysis. However, an analytical resolution of the contact equations is addressed with ROMAX^© and VIBRAGEAR^©.

In this validation phase, the STE results obtained from the multibody finite element approach were compared with those generated by ROMAX^© and VIBRAGEAR^© for standard gear configurations. Despite the higher computational cost of the multibody approach, the necessity of using the proposed method for studying lightweight gears is summarized as follows:

• As lightweight gears are most likely compliant compared to standard ones, a tooth is submitted to higher deflection under load, leading to corner contact at early torque stages.

Gathering all these aspects justify the use of the presented method for studying lightweight gears as to be discussed in Section 4.

This section examines the impact of holes in the gear blank on the STE response, using a driven gear with eight holes as an example. Table 5 provides the parameters of the holes. The gear pair achieves a mass reduction of 21%. The torque applied is 50 N m, the same as in the case of standard gears discussed in section 3.

Figure 10A illustrates the STE fluctuations of the standard gear pair throughout a complete rotation. The frequency spectrum, displayed in Figure 10B, highlights the mesh harmonic and its multiples.

The presence of the holes significantly alters the curve shape, as shown in Figure 10C, compared to the standard configuration. This material discontinuity introduces additional low-frequency harmonic components. An angular frequency equal to the number of holes in the lightweight gear corresponds to the eighth order of the gear pair rotation.

A similar effect was discussed by Shweiki et al. (Shweiki et al., 2019), where their experimental STE investigation of a lightweight gear with three slots demonstrated a low-frequency contribution corresponding to the third order of the gear pair rotation, equivalent to the number of slots.

Figure 10D corroborates these observations, revealing an eighth-order harmonic with an amplitude surpassing the meshing frequency. Additionally, the harmonic content manifests the presence of sidebands that can be expressed as H s b = H z + H h with H z and H h denote the meshing and hole harmonics, respectively.

Figure 11 illustrates the impact of material removal from the blank on the local and global deformation of the lightweight gear with holes. The absence of material in the direction of the meshing forces results in more significant deflections than when the meshing teeth are close to a web region with no gaps. Furthermore, Figure 11B demonstrates that stresses are more widely distributed, with certain regions exhibiting higher values, such as the hole’s circumference near the meshing area, where concentrated stresses reach 400 MPa, as opposed to 100 MPa in Figure 11A. This discrepancy indicates a stiffer state of the gear and explains the STE curve fluctuations.