Spiral bevel gears transmit torque between intersecting shafts through curved teeth that roll and slide under high contact pressure. In aerospace transmissions, their flanks are carburized to produce a hard, wear-resistant case over a tough core. The subsequent quench and temper inevitably distort both the thin, asymmetric gear body and the generated tooth surface. The resulting error is not a uniform shrinkage that can be corrected by a single offset. It changes the contact ellipse, the loaded transmission error, the backlash, and the edge-load distribution. A useful prediction method therefore has to connect three items that are often treated separately: the machine-generated tooth surface, the carbon diffusion field produced by carburizing, and the loaded contact state during meshing.
Why heat-treatment distortion is a contact problem
For spiral bevel gears, the tooth surface is not defined by a simple analytic profile. It is the envelope of a cutter family swept by a multi-axis generator. The surface geometry is therefore controlled by machine settings such as roll ratio, radial cutter position, vertical offset, sliding base, axial setting, root-cone mounting angle, cutter tilt, cutter swivel, and modified-roll coefficients. Small changes in these settings alter the curvature and the contact pattern. Heat treatment adds another layer: carburizing changes the near-surface carbon profile, and the subsequent thermal cycle introduces local strain and stress. Because the gear body, web, hub, and shaft are not symmetric, the distortion field varies strongly from the tooth flank to the shaft.
A full commercial finite-element simulation of carburizing, quenching, and loaded contact can resolve these effects with high fidelity, but it is expensive to repeat inside a machine-setting optimization loop. The semi-analytical route replaces repeated full-field heat-treatment runs with a parameterized mapping from machine settings and process time to tooth-flank deformation. The contact and diffusion problems are still solved numerically; the “semi-analytical” character comes from representing the generated surface and its sensitivity to machine settings analytically, while the local contact and thermal contributions are assembled through a finite-element update.
The coupled physics: carbon diffusion plus loaded contact
During carburizing, carbon enters the steel surface from the furnace atmosphere and diffuses inward. The basic diffusion equation is
$$
\frac{\partial C}{\partial t}
= \frac{\partial}{\partial x_i}\left( D(T,C)\,\frac{\partial C}{\partial x_i} \right),
\qquad
D(T,C)=D_{0.4}\exp\left(-\frac{Q}{RT}\right)\exp\left[-B(0.4-C)\right],
$$
where \(C\) is the carbon mass fraction, \(t\) is time, \(x_i\) is the diffusion direction, \(T\) is temperature, \(D_{0.4}\) is the diffusion coefficient at a carbon mass fraction of 0.4, \(Q\) is the diffusion activation energy, \(R\) is the gas constant, and \(B\) is a constant describing the concentration dependence. The surface mass-transfer coefficient is written as \(\beta=\beta_0\exp[-E/(RT)]\), with \(\beta_0\) a constant and \(E\) an activation energy. The boundary condition is governed by the difference between the furnace carbon potential \(C_e\) and the surface carbon content \(C_s\), so the carbon profile is both temperature- and time-dependent.
At the same time, the gear pair is loaded. The contact is unilateral: a point either remains open with zero force or closes with a positive contact force. For \(N\) candidate contact points, the discrete complementarity relation can be written as
$$
[C_c]\{F\}+\{D\}-\Theta\{\mathbf 1\}-[I]\{Y\}=\{0\},
\qquad
\sum_{K=1}^{N}F_K=P,
$$
where \([C_c]\) is the contact compliance matrix, \(\{F\}\) is the contact-force vector, \(\{D\}\) is the initial gap, \(\Theta\) is the rigid-body approach, \(\{Y\}\) is a slack vector, and \(P\) is the external normal load. The inequality form is converted into a linear-complementarity problem and solved by linear programming at each instantaneous elliptical contact center.
The coupling appears because the thermal contribution changes the local gap and contact pressure, while the contact pressure changes the local deformation. The total contact force at a surface point is therefore not obtained by adding two independent analyses. It is represented as
$$
F_\Omega(\phi_i,\theta,t)=F_e(\phi_i,\theta,t)+F(\phi_i,\theta,t),
$$
where \(F_e\) is the force contribution associated with the thermal/carburizing deformation and \(F\) is the mechanical contact contribution. The variables \(\phi_i\), \(\theta\), and \(t\) denote the machine-setting parameter, the tooth-surface parameter, and process time. This coupling is the central difference between a conventional “distortion plus contact” workflow and a contact-consistent heat-treatment prediction.
Solution sequence for the whole tooth surface
The numerical loaded tooth contact analysis (NLTCA) projects the two tooth surfaces onto a common normal grid. The grid is zero-padded, and the elastic influence coefficients are computed with the discrete-convolution and fast Fourier transform (DC-FFT) method. The contact solver then determines the instantaneous contact area, the pressure distribution, and the rigid-body approach under the applied load. The grid density \(\Delta L\) and the rotation step \(\Delta\phi_{CP}\) are controlled variables; they determine the resolution of the contact ellipse and the smoothness of the loaded contact path.
After the contact forces are obtained, the thermal and mechanical contributions are assembled into a finite-element update. The Newton–Raphson step has the usual form
$$
(K_\Omega)^i_{t+\Delta t}\,\Delta \mathbf U_e
= (P_\Omega)_{t+\Delta t}-(R_\Omega)^i_{t+\Delta t},
$$
where \(K_\Omega\) is the global stiffness matrix, \(P_\Omega\) is the equivalent nodal load vector, \(R_\Omega\) is the internal reaction vector, and \(\Delta \mathbf U_e\) is the displacement increment. The deformation is resolved in three local directions: radial \(R\), axial \(Z\), and tooth-normal \(U\). The vector resultant is
$$
\mathbf u = u_R\,\mathbf e_R+u_Z\,\mathbf e_Z+u_U\,\mathbf e_U,
\qquad
u_{\mathrm{total}}=|\mathbf u|.
$$
The sign of each component is retained because a positive radial displacement and a negative axial displacement can partially cancel in the resultant. Reporting only a scalar shrinkage would lose this directional information.
Demonstration case and baseline contact results
The demonstration uses a non-orthogonal spiral bevel gear pair representative of a helicopter transmission. The shaft angle is 52°, the pinion has 31 teeth, the gear has 38 teeth, the normal module is 5.3 mm, the face width is 32 mm, the pressure angle is 22.5°, the outer cone distance is 208.81 mm, and the spiral angle is 35°. The pinion is right-hand and the gear is left-hand. The material is 9310 steel, and the finite-element model uses 18,000 eight-node brick elements (C3D8). Nine instantaneous contact-ellipse centers along the meshing cycle are used as sampling positions.
| Parameter | Pinion | Gear |
|---|---|---|
| Number of teeth | 31 | 38 |
| Normal module / mm | 5.3 | 5.3 |
| Face width / mm | 32 | 32 |
| Pressure angle / (°) | 22.5 | 22.5 |
| Outer cone distance / mm | 208.81 | 208.81 |
| Spiral angle / (°) | 35 | 35 |
| Hand | Right | Left |
| Shaft angle / (°) | 52 | |
| Material | 9310 steel | |
| Finite-element mesh | 18,000 C3D8 brick elements | |
Without coupling, the NLTCA baseline gives separated contact and thermal deformation values. For the pinion, the contact deformation reaches 0.004645 mm at position 9, is 0.004326 mm at position 1, and falls to 0.0002905 mm at position 5. For the gear, the contact deformation is 0.0055463 mm at position 1 and 0.0021249 mm at position 8, with an overall decrease from position 1 to position 9. These values are useful as a reference, but they do not include the feedback between carburizing-induced distortion and the loaded contact state. The coupled prediction is therefore the relevant result for manufacturing compensation.
Coupled deformation results
Selected values at the first and ninth contact-ellipse centers are summarized below. The first position is near the beginning of meshing and the ninth is near the end of the sampled cycle.
| Component | Position 1 | Position 9 | Observed trend |
|---|---|---|---|
| Pinion \(R\) / mm | −0.29639 | −0.41332 | More negative |
| Pinion \(Z\) / mm | −0.21098 | −0.18488 | Less negative |
| Pinion \(U\) / mm | −0.02134 | +0.023334 | Sign reversal; minimum −0.000882 at position 5 |
| Pinion total / mm | 0.364532 | 0.453224 | Increases |
| Gear \(R\) / mm | −0.44409 | −0.52285 | More negative |
| Gear \(U\) / mm | +0.040989 | +0.020768 | Decreases |
| Gear total / mm | 0.462923 | 0.532417 | Increases |
The pinion does not behave like a uniformly shrinking body. In the radial direction, the contact-center deformation becomes more negative from −0.29639 mm at position 1 to −0.41332 mm at position 9, and its magnitude increases from the tooth toe toward the heel. The axial component moves in the opposite sense: it changes from −0.21098 mm to −0.18488 mm, so its magnitude decreases as meshing proceeds. The tooth-normal component is the most position-sensitive: it starts at −0.02134 mm, passes through −0.000882 mm at position 5, and ends at +0.023334 mm. This sign reversal means that a fixed offset or a single scalar shrinkage cannot represent the distortion field.
The vector resultant still grows from 0.364532 mm to 0.453224 mm because the radial and axial contributions dominate. Spatially, the smaller deformations are located on the shaft, hub, and web, while the tooth flank and the upper part of the shaft show larger values. The total deformation is largest near the tooth-tip/heel region. Along the shaft, the middle region is least affected and the lower region most affected. The shaft–hub connection shows only modest change, whereas the flank changes strongly from the toe to the heel.
For the gear, the radial component also becomes more negative, from −0.44409 mm to −0.52285 mm, and the deformation follows the contact path. The axial component decreases overall from the first to the ninth position and can change sign across the tooth, so it is not a scaled copy of the radial field. The tooth-normal component falls from 0.040989 mm to 0.020768 mm. The total deformation nevertheless rises from 0.462923 mm to 0.532417 mm. The largest total values are concentrated near the tooth-tip/heel region; the shaft shows a gradient from top to bottom, and the flank deformation varies strongly across the tooth.
These results have a practical consequence for compensation. If a manufacturer corrects only the radial shrinkage, the axial and normal components remain uncompensated. A correction based on the average distortion over the whole flank can also be misleading because the sign of \(U\) changes along the meshing path. The machine-setting adjustment should therefore be evaluated against the coupled vector field, not against a single dimensional error.
Validation and scope
The tooth-normal deformation was compared with an earlier adaptive data-driven prediction for the same gear pair. The maximum relative difference is 0.87% for the pinion at position 2 and 0.80% for the gear at position 8. All sampled positions are within 1%. This supports the numerical consistency of the semi-analytical formulation and its ability to reproduce the coupled trend at the contact centers.
The comparison is a model-to-model check, not an experimental validation. It does not replace distortion measurements from a quench trial on the same lot, fixture, and furnace cycle. The reported case is also limited to one non-orthogonal pair in 9310 steel. Transfer to other gear sizes, case depths, quenching paths, and fixture constraints requires recalibration of the thermal and contact inputs. The formulation treats the carburizing contribution through carbon diffusion and thermo-mechanical loading; it does not present a full phase-transformation model for the quenching sequence. For final qualification, commercial finite-element analysis and physical distortion measurements remain the reference.
Practical implications
The main value of the semi-analytical framework is that it links machine settings to a contact-consistent distortion field at design-iteration speed. It can be used to screen machine-setting changes before a full heat-treatment simulation, to identify the dominant direction of distortion for a given flank region, and to support pre-compensation of spiral bevel gears before grinding or lapping. Because the contact solution also provides pressure and approach, the same framework can be extended to contact pressure, contact stress, and loaded transmission error without changing the underlying contact model. For spiral bevel gears that must meet both geometric accuracy and physical performance targets, this coupling-aware view is more useful than treating heat-treatment distortion and loaded contact as separate problems.