# Contact simulation of tooth flanks using Isogeometric Analysis

## Subtitle: Leveraging the power of high order, increased smoothness and curved geometries

Dipl. Ing. **C. Karampatzakis**, Aristotle University of Thessaloniki, Thessaloniki;

Dr. **A. Mantzaflaris**, Centre Inria d'Université Côte d'Azur, Sophia Antipolis;

Prof. **C. Provatidis**, National Technical University of Athens, Athens;

Prof. **A. Mihailidis**, Aristotle University of Thessaloniki, Thessaloniki.

## Zusammenfassung

Die Finite-Elemente-Modellierung (FE) einer Zahnflanke beeinflusst sowohl die Genauigkeit der Ergebnisse als auch die benötigte Rechenzeit und Ressourcen. Bei der FE-Analyse (FEA) wird die Evolvente durch lineare Segmente modelliert, wodurch Abweichungen von ihrer tatsächlichen Geometrie entstehen. Um genauere Ergebnisse zu erzielen, muss das FE-Netz möglichst dicht erzeugt werden. Dadurch erhöht sich die Anzahl der Freiheitsgrade und damit auch der Rechenaufwand. Dieser inhärente Nachteil der FEA wirkt sich ungünstig auf die Berechnung der Druckverteilung aus.

Die isogeometrische Analyse (IGA) ist eine aktuelle Alternative zur FEA. Danach werden anstatt von linearen Segmenten B-Splines verwendet. Diese Technologie wird auch in Computer Aided Design (CAD)-Systemen verwendet, um sowohl die Geometrie als auch das Lösungsfeld zu modellieren. Somit wird beim Übergang vom CAD zum geometrischen Analysemodell kein geometrischer Fehler eingeführt. Ein inhärentes Merkmal von B-Splines ist die Kontinuität zwischen benachbarten Elementen. Darüber hinaus ist das glatte Vektorfeld der Oberfläche-Normalen an jedem Punkt im Inneren von den Elementen bekannt. Dies ist besonders für Kontaktalgorithmen von Vorteil, da dadurch die Notwendigkeit spezieller Glättungstechniken, wie sie häufig in der FEA verwendet werden, entfällt.

In der aktuellen Studie wird ein Stirnradpaar simuliert und der Pressungsverlauf an der Kontaktfläche berechnet. Die Ergebnisse werden hinsichtlich Genauigkeit und Rechenaufwand mit denen der FEA verglichen.
**Abstract**

The Finite Element (FE) modelling of a tooth flank affects both the accuracy of the results as well as the computational time and resources required. In FE Analysis (FEA), the involute curve is modelled by linear segments, thus deviations from its actual geometry are introduced. In order to achieve accurate results, the mesh must be generated as dense as possible. This increases the number of degrees of freedom and the computational cost along with it. This inherent drawback of FEA unfavourably affects the calculation of the pressure distribution of the mating tooth flanks.

Isogeometric Analysis (IGA) is a recent alternative to the FEA. It uses B-Splines, a technology used by Computer Aided Design (CAD) systems, to model the geometry as well as the solution field. Thus, no geometric error is introduced in the transition from CAD to analysis. An inherent characteristic of B-Splines is the continuity between adjacent elements. Furthermore, the smooth normal vector field of the surface is known at every point in the interior of elements. This is particularly advantageous for contact algorithms because this alleviates the need for special smoothing techniques such as those often used in FEA.

In the current study, a spur gear pair is simulated and the pressure at the contact area is calculated. The results are compared to those obtained by FEA in terms of both accuracy and computational cost.

**1. Introduction**

Gears are one of the most prominent parts of mechanical transmission systems. Power is transmitted from one gear to its mating one through the contacting teeth. Thus, contact is a central aspect of a gear's operation. It also plays a crucial role in a gear's wear behaviour. Some of the major failure modes of gears, such as scuffing, pitting and spalling are surface contact induced [1]. As such, it is of great importance for design engineers to accurately predict the contact pressure of a mating gear pair.

Various methods are used by engineers to estimate the pressure developed at the contacting regions. The oldest among them is the analytical solution developed by Heinrich Hertz [2] for the contact between two elastic cylindrical bodies. This theory regards the contacting bodies as half-spaces of constant curvature. This is almost never the case in real-world applications, but it is a good enough first estimate. Analytical solutions have also been provided for various other common indenter shapes [3].
The elastic foundation method [4; 5] offers a compromise between solution time and accuracy of results. This method approximates one of the contacting surfaces as a series of linear springs, independent of one another. As the influence of adjacent spring elements is neglected, it provides only approximate solutions of the contact pressure, but does that in less time than more complex methods.

The Finite Element Method (FEM) has been widely used to simulate contact phenomena. It is more thorough than all aforementioned methods, but this comes at the expense of increasing complexity and solution time. Some shortcomings of the finite element method for contact modelling are that, since the most commonly used basis functions are linear, geometric error is introduced in the discretization process. Furthermore, there is only $C^0$ continuity between adjacent elements, regardless of polynomial degree, and thus the surface's normal and tangent vectors are ambiguous in that areas. This has resulted in the creation of various smoothing techniques that try to alleviate this flaw [6], at the expense of increased computational overhead.

Isogeometric Analysis (IGA) [7] is a recent descendant of the FEM. Their core difference is that IGA utilizes splines as basis functions, both for geometry as well as for the solution field, instead of Lagrange polynomials. Splines have been the basis of choice for most Computer Aided Design (CAD) software for many years. This unified approach has many favorable characteristics in comparison to the classical FEM [8]. Since the same basis is used for design and analysis, no geometric error is introduced in the transition between the two. In addition, splines are invariant under refinement, so remeshing or increasing the degree of a patch does not alter it geometrically or parametrically. Furthermore, since Spline patches typically have inter-element continuity greater than $C^0$, the normal and tangent vector fields are well defined and are smooth on all points of a patch except the corners.

Isogeometric gear modelling has been the focus of some previous studies [9-11]. In these studies, the elasticity equations have been solved using IGA, but the contact has not been modelled in an Isogeometric setting. An equivalent concentrated load has been applied instead. In [12], the contact between a 2D pair of gears was modelled using a hybrid IGA - meshless formulation. In this study, a pure isogeometric formulation is used to simulate contact between a pair of mating gear teeth. All the presented IGA routines have been imple-
mented within the framework of Geometry plus Simulation Modules (G+Smo) [13], an open source isogeometric analysis library written in C++.

Section two provides a short summary of isogeometric analysis and Splines. In section three the method used to model contact is presented and in section four the benchmarks performed to verify the soundness of the implementation are presented. Section five describes the problem setup while the results of the simulations are presented and discussed in the last section.

## 2. Isogeometric Analysis

In this section the Isogeometric paradigm is briefly reviewed. The seminal paper [7] presents the isogeometric method in detail.

### 2.1 Splines

Splines are the primary mathematical entities used by modern CAD systems to represent arbitrary geometries. The most prevalent type of polynomial splines are basis splines (B-Splines) and their rational counterpart, namely non uniform rational B-Splines (NURBS). All aspects regarding the theory and implementation of NURBS have been discussed in [14] and thus the topic is only briefly presented here.

A univariate polynomial B-Spline is computed as the sum of the basis functions which are defined by its knot vector, evaluated at parameter value $\xi$, each weighted by a coefficient. The knot vector is a non-decreasing sequence of values which essentially discretizes the parameter domain into elements. The coefficients of the B-Spline can be d-dimensional, and thus can be interpreted as points in $\mathbb{R}^d$. A B-Spline curve $C(u)$ of polynomial degree $p$ defined by control points $P_i \in \mathbb{R}^d$, is given as:

$$C(\xi) = \sum_{i=1}^{n} N_{i}^{p}(\xi) P_{i}$$

The basis functions $N_i^p(\xi)$ are evaluated using the Cox-de Boor formula [15]

$$N_{i}^{p}(\xi) = \frac{\xi - \xi_{i}}{\xi_{i+p} - \xi_{i}} N_{i}^{p-1}(\xi) + \frac{\xi_{i+p+1} - \xi}{\xi_{i+p+1} - \xi_{i+1}} N_{i+1}^{p-1}(\xi)$$

where $\xi_i$ is the i-th entry of the knot vector. The same principles are extended to represent surfaces (and volumes) by means of the tensor product between univariate bases. By defin-
ing the bivariate basis functions as $N_{i,j}^{p,q}(\xi, \eta) = N_i^p(\xi)N_j^q(\eta)$ we can now define a B-Spline surface in the same manner:

$$S(\xi, \eta) = \sum_{i=1}^{n} \sum_{j=1}^{m} N_{i,j}^{p,q}(\xi, \eta) \mathbf{P}_{i,j}$$

This framework is flexible and easy to use, but it lacks the ability to exactly represent conic sections, which are of great interest in engineering applications. This is achieved with the introduction of rational basis functions. Every basis function, or equivalently every control point, is equipped with a “weight” coefficient $w_{i,j}$. The rational basis functions are then defined as follows:

$$R_{i,j}^{p,q}(\xi, \eta) = \frac{N_i^p(\xi)N_j^q(\eta)w_{i,j}}{\sum_{i=1}^{n} \sum_{j=1}^{m} N_i^p(\xi)N_j^q(\eta)w_{i,j}}$$

A NURBS surface is calculated in the same manner as before but using the rational basis functions instead.

$$S(\xi, \eta) = \sum_{i=1}^{n} \sum_{j=1}^{m} R_{i,j}^{p,q}(\xi, \eta) \mathbf{P}_{i,j}$$

Since the Isogeometric Analysis follows the isoparametric concept [16], the same basis is used to approximate the solution field as well. So, the displacement field $\mathbf{u}$ will take the following form:

$$\mathbf{u}(\xi, \eta) = \sum_{i=1}^{n} \sum_{j=1}^{m} R_{i,j}^{p,q}(\xi, \eta) \mathbf{u}_{i,j}$$

## 2.2 Linear Elasticity

Using this discretization, one can now proceed with the formulation of the problem in the same manner as with classical FEM. Assuming plane strain conditions the constitutive equation for a linear elastic isotropic body is $\sigma(\mathbf{u}) = 2\mu\epsilon + \lambda \text{trace}(\epsilon)I$, with Lame's parameters:

$$\lambda = \frac{E\nu}{(1 + \nu)(1 - 2\nu)}, \quad \mu = \frac{E}{2(1 + \nu)}$$

Strain is given by the infinitesimal strain equation:

$$\epsilon(\mathbf{u}) = \frac{1}{2} [\nabla \mathbf{u} + (\nabla \mathbf{u})^T]$$

The weak form of the equilibrium equation $-\text{div} \sigma(\mathbf{u}) = \mathbf{f}$ subject to: $\mathbf{u} = \mathbf{u}_d$ on $\Gamma_D$, $\sigma(\mathbf{u}) \cdot \mathbf{n} = \mathbf{t}$ on $\Gamma_N$ is then given as:
$$ \int_{\Omega} \sigma(\mathbf{u}) : \varepsilon(\mathbf{v}) d\Omega = \int_{\Omega} \mathbf{f} \cdot \mathbf{v} d\Omega + \int_{\Gamma_N} \mathbf{t} \cdot \mathbf{v} d\Gamma $$

By substituting the spline discretization of $\mathbf{u}$ in the weak form, one can retrieve the discrete system for linear elasticity.

### 3. Contact Modelling

The review of De Lorenzis et al. [17] provides a very concise summary of the available algorithms for contact modelling in an isogeometric setting. The algorithms are grouped into three categories: a) collocation approaches, b) Gauss Point to Surface (GPTS) algorithms and c) mortar methods. The methods have been presented in order of increasing complexity as well as fidelity of the results. In the current study, the GPTS algorithm was preferred as the review noted that it is relatively straightforward to implement, and it provides reasonably accurate results.

#### 3.1 Gauss Point to Surface Algorithm

The Gauss Point to Surface algorithm has been named as such, because the contact inequalities are evaluated and enforced at the integration (Gauss) points of the master-worker interface. Not all integration points are taken into account when integrating the interface equations, but only those in the region where contact between the two bodies has occurred.

In order to distinguish between the regions of the body where contact occurs, it is necessary to define a gap function. Such a function takes a point on the perimeter of body 1, hereafter “worker” body, and returns its distance to the point of the body 2, hereafter “master” body, which is closest to it. Thus, a closest point projection must be performed. In our implementation, a Newton iterative scheme has been employed for the task of closest point projection. The quantities associated with the master point that are evaluated on the result of this projection, are hereafter denoted with a subscript $(\cdot)_p$. Once the two points have been identified, the gap is simply the difference between their current position vectors $\mathbf{d} = \mathbf{x}^{(1)} - \mathbf{x}_p^{(2)}$. The normal gap is the inner product of the gap with the master surface’s normal vector $g_N = \mathbf{d} \cdot \mathbf{n}_p^{(2)}$.
<figure><img src="image_1.png" /><figcaption>Figure 1: Master and worker bodies and the kinematic quantities.</figcaption></figure>

With the use of the normal gap, it is easy to see that penetration has occurred where $g_N < 0$. For two contacting bodies, the normal gap between them must be zero where contact occurs and greater elsewhere, i.e. $g_N \ge 0$. This constraint is enforced with the penalty method. A penalty parameter $\epsilon$ must be specified and with it the contact pressure is identified as $p = \epsilon\langle -g_N \rangle$. The system's total energy $\Pi = \Pi_{elastic} + \Pi_{external}$ is augmented with an extra contact penalty term, integrated on the worker side of the contact interface $c^{(1)}$:

$$ \Pi_{contact} = \frac{1}{2} \int_{c^{(1)}} \epsilon \langle -g_N \rangle^2 da $$

Following the procedure described in detail in [18] but in a two-dimensional setting, we arrive at the variation of the contact penalty integral, which is calculated as:

$$ \begin{aligned} \delta\Pi_{contact} &= \int_{C^{(1)}} \epsilon g_N \delta g_N da = - \sum_{I \in N^{(1)}} \mathbf{R}_I^{(1)} \delta \mathbf{q}_I^{(1)} - \sum_{J \in N^{(2)}} \mathbf{R}_J^{(2)} \delta \mathbf{q}_J^{(2)} \\ \mathbf{R}_I^{(1)} &= - \sum_{i \in G^{(1)}} \epsilon w_i J^{(1)}(u_i) N_I^{(1)}(u_i) [\mathbf{x}^{(1)}(u_i) - \mathbf{x}_p^{(2)}] \\ \mathbf{R}_J^{(2)} &= \sum_{i \in G^{(1)}} \epsilon w_i J^{(1)}(u_i) N_J^{(2)}(u_p) [\mathbf{x}^{(1)}(u_i) - \mathbf{x}_p^{(2)}] \end{aligned} $$

where the index $i$ denotes the integration point belonging to the set $G^{(1)}$ of integration points of the worker side (1) within the contact area $c^{(1)}$. $w_i, J(u_i)$ are the quadrature weight and Jacobian at point $i$.

Since the residuals depend on the current configuration of the bodies, the problem is non-linear, and thus needs to be solved iteratively, in this case with the Newton method. This requires the computation of the contributions of the contact residuals to the tangent stiffness matrix. The process is also presented in [18], but here it has been adapted to the two-dimensional case.
## 3.2 GPTS-2p & 2hp

In the standard GPTS method described above, the contact integral is integrated only on the worker side of the interface. Notice that even the residual of the master body, $R_j^{(2)}$, is integrated on the worker side. This induces an inherent bias between the master and worker sides, in the sense that the quantities of the master side are misrepresented in the contact integral, as they are not evaluated at the master surface's integration points. To alleviate this imbalance, two modifications to the standard method have been proposed.

The first alternative, GPTS two-pass (GPTS-2p), is named as such because the integral is calculated twice, each time switching the role of worker and slave, and the two results are averaged.

$$ \delta\Pi_{contact} \approx \frac{1}{2} \int_{c^{(1)}} \epsilon[\delta\mathbf{x}^{(1)} - \delta\mathbf{x}^{(2)}] \cdot [\mathbf{x}^{(1)} - \mathbf{x}^{(2)}]da + \frac{1}{2} \int_{c^{(2)}} \epsilon[\delta\mathbf{x}^{(2)} - \delta\mathbf{x}^{(1)}] \cdot [\mathbf{x}^{(2)} - \mathbf{x}^{(1)}]da $$

As will be shown later, this approach does not fully address the issue, but tries to lessen its extent by dispersing the error on both sides of the interface.

The second variation of the method, GPTS two-field (also known as two-half-pass), similarly splits the integral in two, evaluated with the roles of master and worker switched. This time only the terms associated with the current integration side, are considered on the respective integral. This is equivalent to introducing two independent pressure fields, one for each side of the interface.

$$ \delta\Pi_{contact} \approx \int_{c^{(1)}} \epsilon \delta\mathbf{x}^{(1)} \cdot [\mathbf{x}^{(1)} - \mathbf{x}^{(2)}]da + \int_{c^{(2)}} \epsilon \delta\mathbf{x}^{(2)} \cdot [\mathbf{x}^{(2)} - \mathbf{x}^{(1)}]da $$

## 4. Benchmarks

To verify the integrity of our implementation of the aforementioned algorithms in G+Smo, two benchmark problems were solved. First the contact patch test which was introduced in [19] and second the classical two-dimensional Hertzian contact [2].

### 4.1 Patch Test

The contact patch test consists of two dissimilar rectangular patches, with the smaller of the two resting on top of the other. An ambient vertical traction of -100 N/mm is applied on all the exposed upper horizontal surfaces of both patches. The lower patch is fixed along the vertical direction, on its bottom side (Figure 2a). The expected solution is to recover a uniform stress field all throughout the two patches.
<figure><img src="image_2.png" /><figcaption>Figure 2: From left to right: a) Setup of patch test, stress field with b) standard method, c) two-field method.</figcaption></figure>

In agreement with the observations of previous studies, the standard and 2p methods fail to pass the patch test exactly. While the bulk mostly behaves as expected, oscillations of the stress are observed within the contact region, as is shown on Figure 2b. These are accredited to the master-worker bias that was discussed previously, and on the dissimilarity of the elements meshes on the interface. The two-field method, on the other hand, does recover a truly uniform stress field on both patches and thus passes the patch test, as it can be seen in Figure 2c.

## 4.2 Hertzian Contact

Hertzian theory provides analytical solutions to the cases of contacting cylinders and contacting spheres, under the assumption of small deformations. For the purpose of this study, the contact between an infinitely long cylinder and a half-space is considered. Due to the nature of the problem, it is examined at a two-dimensional section perpendicular to the cylinder's axis. Exploiting the symmetry of the geometry, only a quarter of the disk is modelled, and the relevant symmetry boundary conditions are applied. The quarter-disk is loaded with a uniformly distributed load on its top edge whereas the plane is fixed along the vertical direction on its lower side, as it can be seen in Figure 3a.
<figure><img src="image_3.png" /><figcaption>Figure 3: From left to right: a) Setup of Hertzian contact benchmark b) resulting von Mises stress field c) comparison of analytic and calculated contact pressure.</figcaption></figure>

Figure 3c shows the resulting contact pressure calculated with the two-field method compared to the one predicted by the analytical formulas. Both axes have been normalized by the analytical values. The oscillation observed close to the edge of the contact area agrees with what has been observed in previous studies. Due to the smooth nature of the basis used to represent the geometry as well as the solution field, the kink between the contacting area and the free surface cannot be captured exactly, thus oscillations occur. This problem can be mitigated by further refining the region close to the edge of the contact zone, thus allowing for a more accurate representation of the transition. Another possibility is to perform adaptive knot insertion to the basis, similar to the rp-FEM method presented in [20]. In this way the continuity of the basis can be lowered to $C^0$ at the crossing point, so that the kink can be represented exactly. Furthermore, a modified version of the GPTS algorithm that aims to improve on the aspect of oscillations has been proposed in [21]. For the present study however, this is not critical for the analysis since the calculated maximum contact pressure, as well as the contact width are in good agreement with the analytical solution.

# 5. Problem setup

## 5.1 Geometry and positioning

The generation of the NURBS geometry for standard involute gears is performed automatically with the established engineering parameters as input. First, two point sets are created by following the procedure described in [22]. The two sets correspond to the involute and fillet regions of the tooth. Then a Bezier spline of user-specified degree, is fitted to each of these sets of points, by means of minimizing the least square error. As is evident the accuracy of the resulting geometry is affected by the choice of degree. These two splines are then
spliced into one, and extended radially from the root to the hub diameter. The resulting spline forms one side of the bivariate patch.

This side is subsequently mirrored, and the bulk of the patch is created by assuming a circular geometry along the other direction. This justifies the need for NURBS basis, so that the circular parts of the gear are represented precisely. The resulting patch constitutes a single tooth of the gear, centered along the Y axis. Once the patches for both gears have been initialized, they need to be appropriately positioned in their mating configuration. To specify the desired position, the distance along the line of action, between the desired contact point and the pitch point is specified and then the rotation angles are computed accordingly.

As contact is a highly localized phenomenon and as has been shown in previous studies [17], a further refinement in the vicinity of the contact point is required. Since the contact point is known a priori from the positioning routine, a zone centered on this point is further refined on both gears by means of knot insertion. An example of the created element mesh created by this automatic procedure can be seen in Figure 4.

The specific gear pair under consideration, is the well-known FZG-A, as given in [23]. The key parameters for the gear pair are summarized in the following table.

<table><caption>Table 1: Parameters of gear pair.</caption><thead><tr><th></th><th>Module</th><th>Teeth</th><th>Profile shift coef.</th><th>Centre distance</th><th>Tip diameter</th></tr></thead><tbody><tr><td>Gear</td><td rowspan="2">4.5</td><td>24</td><td>-0.5</td><td rowspan="2">91.5</td><td>112.5</td></tr><tr><td>Pinion</td><td>16</td><td>0.8532</td><td>88.77</td></tr></tbody></table>

## 5.2 Boundary Conditions

The boundary conditions applied to the pair of gears are as follows. The gear tooth is fixed in all degrees of freedom along its innermost edge as well as on the fictitious edges where it would normally be joined with the rest of the gear. On the pinion, a prescribed rotation is applied on the same edges as on the gear. An applied torque would have been a more appropriate loading for this case, but this would require a constraint on the pinion so that it only rotates around its center. Unfortunately, this sort of constraint mechanism is not yet supported in G+Smo. A schematic of the problem setup can be seen in the following figure.
<figure><img src="image_4.png" /><figcaption>Figure 4: Left: Setup of gear teeth and applied boundary conditions. Right: example of one of the employed NURBS meshes.</figcaption></figure>

**6. Results and Discussion**

The two-field algorithm was applied to meshes of consecutive element subdivision steps, to investigate their convergence behavior. For every such mesh, elements of different types of bases were used, to compare their relative performance. NURBS bases of degree 2 and 3 were examined, as well as linear Lagrange elements, the elements most commonly employed by the Finite Element Method.

<figure><img src="image_5.png" /><figcaption>Figure 5: Error in maximum contact pressure vs. total solution time plot.</figcaption></figure>

As can be seen in Figure 5 both the IGA and FEM models, converge to the same maximum pressure value, so there is good agreement between the two methods. The advantage of
IGA that has been observed in our experiments is that for a given amount of computation time, superior accuracy results with respect to FEM are attained. It is clear from the data, that the second order NURBS basis has converged to less than 1% within the final value already from the simulation that took 2 seconds with 8468 degrees of freedom (DOF), whereas the simulation with linear basis functions requires approximately 9 seconds and 67544 DOF to reach the same result. The contact pressure profile recovered from the two aforementioned simulations is compared in Figure 6b. There is very good agreement between the two both in terms of the maximum pressure, as well as the contact area width.

The cost per degree of freedom for the linear Lagrange basis is much lower than for the IGA bases. Indeed, the higher degree and larger support of NURBS basis functions increases the computational cost of the system's assembly and solution. However, this increased cost translates to improved solution accuracy. In addition, due to the higher degree and continuity of NURBS, quantities that are derived from the solution field, like strains and stresses, can be evaluated and are continuous throughout the domain, which is not the case with the FEM. The high smoothness of the acquired solution with quadratic NURBS basis functions can be observed in Figure 6a.

<figure><img src="image_6.png" /><figcaption>Figure 6: From left to right a) resulting von Mises stress field, b) contact pressure comparison between linear and NURBS bases.</figcaption></figure>
Concluding, NURBS based isogeometric analysis is a promising alternative to the estab-
lished practices. It offers a significantly more precise description of the geometry and attains
results of higher accuracy in less time, compared to the linear basis. Naturally there is still
room for improvement. As the cost per degree of freedom is higher for IGA, using a basis
that allows for local refinement might be beneficial. Candidates are the Truncated Hierar-
chical B-Splines (THB-Splines) and T-Splines. On the aspect of contact formulation, the mor-
tar method has been shown to deliver some advantages over the GTPS algorithms, mainly in
the sense of robustness and absence of the observed oscillatory behavior in the contact
pressure, at the expense of increased computational effort.

References

[1] Vorgerd, J., Tenberge, P. and Joop, M.: Scuffing of cylindrical gears with pitch line velocities up to 100 m/s. Forschung im Ingenieurwesen, 86, (2022), pp. 513–520

[2] Hertz, H.: Über die Berührung fester elastischer Körper. Journal für die reine und an-gewandte Mathematik, 92, (1882), pp. 156–171

[3] Popov, V. L., Heß, M. and Willert, E.: Handbook of Contact Mechanics: Exact Solutions of Axisymmetric Contact Problems. Springer Berlin Heidelberg 2019

[4] Johnson, K. L.: Contact Mechanics. Cambridge University Press 1985

[5] Mukras, S. M. S.: Computer Simulation/Prediction of Wear in Mechanical Components. Advances in Tribology, 2020, (2020)

[6] Neto, D. M., Oliveira, M. C. and Menezes, L. F.: Surface Smoothing Procedures in Computational Contact Mechanics. Archives of Computational Methods in Engineering, 24, (2017), pp. 37–87

[7] Hughes, T. J. R., Cottrell, J. A. and Bazilevs, Y.: Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 194, (2005), pp. 4135–4195

[8] Cottrell, J. A., Hughes, T. J. R. and Bazilevs, Y.: Isogeometric analysis: toward integration of CAD and FEA. Wiley 2009.

[9] Cheng, Q., Yang, G. and Lu, J.: An analysis of gear based on Isogeometric Analysis. Vibroengineering Procedia, 2, (2013), pp. 17–22

[10] Yusuf, O. T., Zhao, G., Wang, W. and Onuh, S. O.: Simulation Based on Trivariate Nurbs and Isogeometric Analysis of a Spur Gear. Strength of Materials, 47, (2015), pp. 19–28
[11] Beinstingel, A., Keller, M., Heider, M., Pinnekamp, B. and Marburg, S.: A hybrid analytical-numerical method based on Isogeometric Analysis for determination of time varying gear mesh stiffness. *Mechanism and Machine Theory*, 160, (2021)

[12] Greco, F., Rosolen, A., Coox, L. and Desmet, W.: Contact mechanics with maximum-entropy meshfree approximants blended with isogeometric analysis on the boundary. *Computers & Structures*, 182, (2017), pp. 165–175

[13] Mantzaflaris, A.: An Overview of Geometry Plus Simulation Modules. *Lecture Notes in Computer Science*, 11989., (2020), pp. Z453-456

[14] Piegl, L. and Tiller, W.: The NURBS Book. Springer Berlin Heidelberg 1995

[15] de Boor, C.: Package for Calculating with B-Splines. *SIAM Journal on Numerical Analysis*, 14, (1977), pp. 441–472

[16] Irons, B. M.: Engineering applications of numerical integration in stiffness methods. *American Institute of Aeronautics and Astronautics Journal*, 4, (1966), pp. 2035–2037

[17] De Lorenzis, L., Wriggers, P. and Hughes, T. J. R.: Isogeometric contact: a review. *GAMM-Mitteilungen*, 37, (2014), pp. 85–123

[18] Lu, J.: Isogeometric contact analysis: Geometric basis and formulation for frictionless contact. *Computer Methods in Applied Mechanics and Engineering*, 200, (2011), pp. 726–741

[19] Taylor, L. R. and Papadopoulos, P.: On a patch test for contact problems in two dimensions. *Nonlinear Computational Mechanics*, (1991), pp. 690–702

[20] Franke, D., Düster, A., Nübel, V. and Rank, E.: A comparison of the h-, p-, hp-, and rp-version of the FEM for the solution of the 2D Hertzian contact problem. *Computational Mechanics*, 45, (2010), pp. 513–522

[21] Temizer, İ., Wriggers, P. and Hughes, T. J. R.: Contact treatment in isogeometric analysis with NURBS. *Computer Methods in Applied Mechanics and Engineering*, 200, (2011), pp. 1100–1112

[22] Lopez, M. A. and Wheway, R. T.: A Method for Determining the AGMA Tooth Form Factor from Equations for the Generated Tooth Root Fillet. *Journal of Mechanisms, Transmissions, and Automation in Design*, 108, (1986), pp. 270–279

[23] DIN 51354-1: Prüfung von Schmierstoffen; FZG-Zahnrad-Verspannungs-Prüfmaschine; Allgemeine Arbeitsgrundlagen, (1990)