Spur gears are widely used mechanical transmission components in modern equipment, relying on tooth surface contact to achieve power transmission. The mechanical performance of the tooth surfaces directly impacts the stability and reliability of gear transmission systems, as well as their service life. Currently, finite element analysis (FEA) serves as an effective numerical method for analyzing tooth surface contact in gear transmission processes. However, the accuracy of contact calculations based on FEA highly depends on mesh quality. During mesh generation, the smooth tooth surfaces obtained from geometric design are approximately discretized into low-order continuous analysis elements for numerical simulation. Overly coarse meshes introduce geometric approximation errors and reduce computational accuracy, while overly fine meshes lead to an exponential increase in computational cost. Therefore, exploring higher-precision and more efficient contact calculation methods for spur gears holds significant practical value.
Isogeometric analysis (IGA), proposed by Hughes et al., offers a new approach to addressing gear contact problems. IGA utilizes the high-order continuous NURBS basis functions from CAD models as shape functions, treating NURBS curve segments, surface patches, or solid blocks as mesh elements and NURBS control points as nodes. This avoids the mesh generation process in FEA. Consequently, IGA not only achieves consistency between geometric and analysis models but also naturally meets the high-order continuity requirements for contact boundaries, facilitating accurate contact force transmission and rapid convergence in contact calculations. Compared to FEA, IGA exhibits unique advantages in computational accuracy and convergence speed. This paper explores parametric solid modeling and contact calculation methods for spur gears within the isogeometric framework from a CAD/CAE integration perspective. First, geometric design parameters are given to drive the automatic parametric solid modeling of spur gears under specified error conditions. Then, linear elastic tooth surface contact isogeometric calculations based on the node-to-surface method and contact search strategies are investigated. Finally, the influence of gear geometric design parameters on contact calculations is explored, aiming to provide a new alternative solution for contact design and strength verification of spur gears and lay the foundation for the integration and intelligence of gear design and analysis.
The basic parameters of spur gears include the number of teeth $z$, module $m$, pressure angle $\alpha$, addendum coefficient $h_a^* = 1$, dedendum coefficient $c^* = 0.25$, and tooth width $B$. The gear profile curve consists of four segments: the addendum circle arc $AB$, the involute curve $BC$, the root transition curve $CD$, and the dedendum circle arc $DE$. The parametric equations for these segments are derived as follows.
The addendum circle arc $AB$ is expressed as:
$$
\begin{cases}
x = r_a \sin \psi \\
y = r_a \cos \psi
\end{cases}
$$
where $\psi$ is the variable parameter with range $0 \leq \psi \leq \varphi_a$, $r_a$ is the addendum radius given by $r_a = \frac{1}{2} m (z + 2h_a^*)$, and $\varphi_a$ is the central angle corresponding to the addendum tooth thickness, calculated as $\varphi_a = \frac{\pi}{2z} + 2(\text{inv} \alpha_a – \text{inv} \alpha)$, with $\alpha_a = \arccos(\frac{r_b}{r_a})$ being the pressure angle at the addendum circle, and $r_b = r \cos \alpha$ being the base circle radius, where $r = \frac{1}{2} m z$ is the pitch circle radius.
The involute curve $BC$ is parameterized as:
$$
\begin{cases}
x = r_K \sin \varphi_K \\
y = r_K \cos \varphi_K
\end{cases}
$$
where $r_K = \frac{r_b}{\cos \alpha_K}$, $\varphi_K = \theta + \text{inv} \alpha_K – \text{inv} \alpha$, $\theta = \frac{\pi}{2z}$, and $\alpha_K$ is the variable parameter with range $\alpha_C \leq \alpha_K \leq \alpha_B$. The values of $\alpha_C$ and $\alpha_B$ are determined from the gear geometry.
The root transition curve $CD$ is given by:
$$
\begin{cases}
x = \left( r – \frac{a}{\sin \psi} \right) \sin(\varphi + \psi – \alpha) + \rho \cos(\varphi + \psi – \alpha) \\
y = \left( r – \frac{a}{\sin \psi} \right) \cos(\varphi + \psi – \alpha) – \rho \sin(\varphi + \psi – \alpha)
\end{cases}
$$
where $\varphi = \arctan(\frac{a}{b})$, $a = h_a^* m + c^* m – \rho$, $b = \frac{\pi m}{4} + h_a^* m \tan \alpha + \rho \cos \alpha$, $\rho = \frac{c^* m}{1 – \sin \alpha}$, and $\psi$ is the variable parameter with range $\alpha \leq \psi \leq \frac{\pi}{2}$.
The dedendum circle arc $DE$ is expressed as:
$$
\begin{cases}
x = r_f \sin \psi \\
y = r_f \cos \psi
\end{cases}
$$
where $r_f = \frac{1}{2} m (z – 2h_a^* – 2c^*)$ is the dedendum radius, and $\psi$ ranges from $\arctan(\frac{x_D}{y_D})$ to $\frac{\pi}{z}$.
NURBS curves, surfaces, and solids provide a mathematical framework for parametric modeling. A NURBS curve of degree $p$ is defined as:
$$
\mathbf{C}(\xi) = \frac{\sum_{i=0}^{n} N_{i,p}(\xi) \omega_i \mathbf{Q}_i}{\sum_{i=0}^{n} N_{i,p}(\xi) \omega_i} = \sum_{i=0}^{n} R_{i,p}(\xi) \mathbf{Q}_i
$$
where $N_{i,p}(\xi)$ are the B-spline basis functions, $\mathbf{Q}_i$ are control points, $\omega_i$ are weights, and $R_{i,p}(\xi)$ are the NURBS basis functions. Similarly, a NURBS surface is given by:
$$
\mathbf{S}(\xi, \eta) = \sum_{i=0}^{n} \sum_{j=0}^{m} R_{i,p;j,q}(\xi, \eta) \mathbf{Q}_{i,j}
$$
and a NURBS solid by:
$$
\mathbf{V}(\xi, \eta, \zeta) = \sum_{i=0}^{n} \sum_{j=0}^{m} \sum_{k=0}^{l} R_{i,p;j,q;k,r}(\xi, \eta, \zeta) \mathbf{Q}_{i,j,k}
$$
For spur gears, the addendum and dedendum circle arcs are exactly represented using quadratic NURBS curves, while the involute and root transition curves are approximated using least-squares fitting with cubic NURBS curves. The fitting process minimizes the objective function:
$$
f = \sum_{i=1}^{M-1} \left\| \mathbf{q}_i – \mathbf{C}(\xi_i) \right\|^2
$$
where $\mathbf{q}_i$ are data points, and $\mathbf{C}(\xi_i)$ is the NURBS curve evaluated at parameter $\xi_i$. The maximum relative fitting error $E_f$ is defined as:
$$
E_f = \max_{0 \leq i \leq M} \frac{\min \left\| \mathbf{q}_i – \mathbf{C}(\xi) \right\|}{\max(x_{\text{max}} – x_{\text{min}}, y_{\text{max}} – y_{\text{min}})}
$$
For a spur gear with $m=2$, $z=19$, and $\alpha=20^\circ$, the involute and root transition curves are fitted with 5 control points each, achieving $E_f < 0.001$.
The spur gear plane model is constructed by decomposing the single-tooth region into five patches. For instance, patch① has its right boundary as the fitted involute curve, left boundary as the symmetric curve, and upper and lower boundaries as circular arcs. The control points for intermediate columns are derived based on symmetry and geometric constraints. The 3D solid model is generated by offsetting the 2D surface model along the $z$-axis by the tooth width $B$, and the full gear model is obtained by rotating the single-tooth solid around the gear center.
For frictionless contact analysis using IGA, consider two deformable bodies in contact, denoted as slave (s) and master (m). The current configuration at time $t$ is $\Omega^t$, with boundary $\Gamma^t = \partial \Omega^t$. The reference configuration at $t=0$ is $\Omega^0$. The contact interface $\Gamma_C^t$ satisfies $\Gamma_C^t = \Gamma_s^t \cap \Gamma_m^t$. The spatial coordinates are related by $\mathbf{x} = \mathbf{X} + \mathbf{u}(\mathbf{X}, t)$, where $\mathbf{u}$ is the displacement. The normal contact gap $g_N$ is defined as:
$$
g_N = (\mathbf{x}^s – \mathbf{x}^m) \cdot \mathbf{n}^m
$$
where $\mathbf{x}^s$ is a point on the slave surface, $\mathbf{x}^m$ is its closest projection on the master surface, and $\mathbf{n}^m$ is the unit outward normal at $\mathbf{x}^m$. The contact force $\mathbf{t}_c$ is decomposed into normal and tangential components: $\mathbf{t}_c = t_N \mathbf{n}^m + t_T \mathbf{a}^m$, with $t_T = 0$ for frictionless contact. The normal contact constraints are:
$$
g_N \geq 0, \quad t_N \leq 0, \quad g_N t_N = 0
$$
The penalty method is used to enforce these constraints, where the normal contact force is $t_N = \lambda g_N$, with $\lambda$ being the penalty parameter. The total potential energy of the contact system is:
$$
\Pi(\mathbf{u}) = \Pi_{\text{int}}(\mathbf{u}) + \Pi_{\text{ext}}(\mathbf{u}) + \Pi_c(\mathbf{u})
$$
where $\Pi_{\text{int}}$ is the internal elastic energy, $\Pi_{\text{ext}}$ is the external work, and $\Pi_c$ is the penalty term:
$$
\Pi_c(\mathbf{u}) = \frac{1}{2} \int_{\Gamma_C^t} \lambda g_N^2 d\Gamma
$$
The variation of the potential energy yields the weak form:
$$
\delta \Pi = \delta \Pi_{\text{int}} + \delta \Pi_{\text{ext}} + \delta \Pi_c = 0
$$
where
$$
\delta \Pi_c = \int_{\Gamma_C^t} \lambda g_N \delta g_N d\Gamma
$$
For 2D contact, the linearization of $\delta \Pi_c$ involves terms like $\delta g_N$, $\Delta g_N$, and $\Delta \delta g_N$, which are derived using geometric relations. For 3D contact, similar derivations apply with additional parameter directions.
Contact detection involves global and local searches. The local search uses the closest point projection method. For a slave point $\mathbf{x}^s$, the projection on the master surface $\mathbf{S}(\xi, \eta)$ is found by solving:
$$
\begin{cases}
\mathbf{S}_{,\xi} \cdot (\mathbf{x}^s – \mathbf{S}) = 0 \\
\mathbf{S}_{,\eta} \cdot (\mathbf{x}^s – \mathbf{S}) = 0
\end{cases}
$$
using Newton’s method. The iterative equations for parameters $\xi$ and $\eta$ are:
$$
\xi^{i+1} = \xi^i – \frac{g h_{,\eta} – h g_{,\eta}}{g_{,\xi} h_{,\eta} – g_{,\eta} h_{,\xi}}, \quad \eta^{i+1} = \eta^i – \frac{h g_{,\xi} – g h_{,\xi}}{g_{,\xi} h_{,\eta} – g_{,\eta} h_{,\xi}}
$$
where $g = \mathbf{S}_{,\xi} \cdot (\mathbf{x}^s – \mathbf{S})$ and $h = \mathbf{S}_{,\eta} \cdot (\mathbf{x}^s – \mathbf{S})$.
Numerical examples validate the proposed method. For static analysis, a 2D spur gear tooth with $z=19$, $m=2$, $\alpha=20^\circ$ is fixed at the base and subjected to a distributed load $P_0 = -10 \text{N/mm}^2$ at the top. Material properties are $E = 207 \text{GPa}$, $\nu = 0.25$. The IGA and FEA results are compared in Table 1.
| Quantity | Unit | IGA | FEA |
|---|---|---|---|
| Max $U_1$ | mm | $1.0138 \times 10^{-5}$ | $1.1320 \times 10^{-5}$ |
| Min $U_2$ | mm | $-1.5500 \times 10^{-4}$ | $-1.5527 \times 10^{-4}$ |
| Max $U$ | mm | $1.5500 \times 10^{-4}$ | $1.5527 \times 10^{-4}$ |
| Max $S_1$ | MPa | 0.5537 | 0.5636 |
| Min $S_2$ | MPa | -11.1710 | -10.5517 |
| Max $S_{12}$ | MPa | 3.1295 | 2.7768 |
The convergence of maximum displacement magnitude with degrees of freedom is shown in Figure 1, indicating faster and smoother convergence for IGA compared to FEA.
For 3D static analysis, the 2D model is extruded to a tooth width $B=1 \text{mm}$. The results are summarized in Table 2.
| Quantity | Unit | IGA | FEA |
|---|---|---|---|
| Max $U$ | mm | $1.5629 \times 10^{-4}$ | $1.5527 \times 10^{-4}$ |
| Max Mises Stress | MPa | 9.8602 | 9.8602 |
Contact analysis is performed for a spur gear pair with $z_1=19$, $z_2=38$, $m=2$, $\alpha=20^\circ$. The driving gear is rotated by $\theta = \pi/1600 \text{rad}$ while the driven gear is fixed. The penalty parameter is $\lambda = 1 \times 10^7$. The 2D contact results are shown in Table 3.
| Quantity | Unit | IGA | FEA | Relative Error |
|---|---|---|---|---|
| Max $U$ | mm | 0.0324 | 0.0324 | 0% |
| Max Mises Stress | MPa | 2215.7302 | 2212.8679 | 0.13% |
For 3D contact analysis, the tooth width is $B=1 \text{mm}$, and $\lambda = 1 \times 10^8$. The results are in Table 4.
| Quantity | Unit | IGA | FEA | Relative Error |
|---|---|---|---|---|
| Max $U$ | mm | 0.0324 | 0.0324 | 0% |
| Max Mises Stress | MPa | 2200.2902 | 1776.9370 | 23.82% |
The convergence of maximum contact stress with degrees of freedom for 2D contact is plotted in Figure 2, showing superior convergence for IGA.
In conclusion, this paper presents a NURBS-based parametric solid modeling and isogeometric contact analysis method for spur gears. The approach enables high-precision geometric modeling and efficient contact calculations, addressing limitations of traditional FEA. Future work will focus on incorporating friction into spur gear contact analysis, considering manufacturing errors for high-precision simulations, and applying IGA to large deformation elasto-plastic problems like stamping formation. The integration of design and analysis for spur gears paves the way for intelligent and efficient gear transmission systems.