Spiral bevel gears transmit motion between intersecting shafts in aerospace and automotive drivetrains. Their flanks are generated on computer-controlled machines, then carburized and quenched to produce a hard, wear-resistant case. The same thermochemical cycle that creates the case also distorts the tooth flank. Because the loaded contact pattern is sensitive to small geometric changes, distortion can shift load away from the intended region, raise local contact pressure, and alter transmission error. Predicting that distortion before hard finishing is therefore a process-planning problem as well as a materials problem.
The Coupled Nature of Carburizing and Meshing
Carburizing introduces a carbon gradient, a hardness gradient, and a non-uniform volumetric response during cooling. The resulting distortion is not a simple uniform expansion; it depends on tooth geometry, case depth, cooling path, and the constraint imposed by the surrounding web, hub, and shaft. In parallel, the loaded meshing state is a nonlinear contact problem: the contact footprint, pressure distribution, and approach depend on the actual flank geometry after distortion.
These two effects should not be treated as independent post-processing steps. A distorted flank changes the contact, and the contact state determines where mechanical and thermal deformations must remain consistent. In the coupled formulation, the carburizing-induced deformation and the loaded contact deformation are solved on the same finite-element mesh, so the instantaneous contact ellipse, contact force, and distortion field are mutually consistent. The term “meshing” here refers to the loaded contact solution used to evaluate the final flank state; it does not imply that carburizing occurs during gear operation.
The practical reason for coupling is compensation. Machine tool settings define the unloaded flank before heat treatment. If the mapping from those settings to the final coupled distortion is known, settings can be screened or corrected before expensive cutting and heat-treatment trials. A purely geometric correction of the cold flank may be insufficient because the loaded contact region is where the coupling has the greatest effect.
Generating the Tooth Flank from Machine Settings
Spiral bevel gears do not have a standard involute profile. Their flanks are generated as the envelope of a cutter surface moving according to the machine kinematics. A flank point can be written as
\mathbf{p}(\phi,\mu,\theta)
=
\mathbf{M}_{b1}(\phi_P,X_D,\gamma_m,X_B,E_M,S_r,\sigma,\zeta,R_a)\,
\mathbf{r}_1(\mu,\theta),
\qquad
\mathbf{n}_p\cdot\mathbf{v}_{b-1}=0.
$$
Here, \(\mu\) and \(\theta\) are cutter surface parameters, \(\phi\) is the cradle angle, and \(\mathbf{r}_1\) is the cutter position vector. The transformation \(\mathbf{M}_{b1}\) contains the machine settings: axial cradle position \(X_D\), root cone mounting angle \(\gamma_m\), bed position \(X_B\), vertical offset \(E_M\), radial cutter position \(S_r\), cutter tilt \(\sigma\), cutter swivel \(\zeta\), and roll ratio \(R_a\). The envelope condition \(\mathbf{n}_p\cdot\mathbf{v}_{b-1}=0\) enforces tangency between the cutter and the workpiece; it is solved iteratively to obtain discrete flank points.
The validation case is a non-orthogonal aerospace spiral bevel gear pair with a 52° shaft angle. Its basic geometry is listed below.
| Parameter | Pinion | Gear |
|---|---|---|
| Number of teeth | 31 | 38 |
| Normal module | 5.3 mm | 5.3 mm |
| Face width | 32 mm | 32 mm |
| Pressure angle | 22.5° | 22.5° |
| Outer cone distance | 208.81 mm | 208.81 mm |
| Pitch cone angle | 23.167° | 28.833° |
| Root cone angle | 22.4° | 27.733° |
| Face cone angle | 24.267° | 29.6° |
| Spiral angle | 35° | 35° |
| Hand | Right | Left |
| Addendum | 5.09 mm | 3.44 mm |
| Dedendum | 4.65 mm | 6.30 mm |
| Shaft angle | 52° | 52° |
The corresponding duplex helical grinding settings are summarized below. These are machine inputs, not nominal CAD dimensions, and they control the generated flank curvature and contact position.
| Setting | Symbol | Pinion | Gear |
|---|---|---|---|
| Roll ratio | \(R_a\) | 2.352220 | 2.073214 |
| Radial cutter position | \(S_r\) | 152.07209 mm | 163.35610 mm |
| Vertical offset | \(E_M\) | 1.95331 mm | 0 |
| Bed position | \(X_B\) | −37.13876 mm | −3.95200 mm |
| Axial cradle position | \(X_D\) | −15.90914 mm | 0 |
| Root cone mounting angle | \(\gamma_m\) | 34.34° | 27.73° |
| Cutter tilt | \(\sigma\) | 17.46° | 0.63° |
| Cutter swivel | \(\zeta\) | 106.62° | 14.83° |
| Angular cutter position | \(\phi\) | 61.19° | 49.83° |
| Helical correction coefficient | \(\chi_H\) | 38.95265 | 0 |
A three-dimensional finite-element model of the gear pair, including the teeth, web, hub, and shaft, was built with 18,000 C3D8 solid elements. This mesh is used for both the contact solution and the thermal deformation field, which is necessary for the coupling.
Loaded Contact Analysis with DC-FFT
Numerical loaded tooth contact analysis (NLTCA) projects the two flanks onto a common tangent plane, interpolates them onto a rectangular grid, and uses a discrete-convolution fast Fourier transform (DC-FFT) to compute elastic influence coefficients. The contact constraints are
[C]\{F\}+\{D\}-\Theta\{e\}-[I]\{Y\}=\{0\},
\qquad
\sum_{K=1}^{N}F_K=P.
$$
In this expression, \([C]\) is the compliance matrix, \(\{D\}\) is the initial gap, \(\Theta\) is the rigid-body approach, \(\{Y\}\) contains slack variables for non-contact points, and \(P\) is the external load. The unilateral contact inequalities are solved by an improved simplex method. The active set gives the instantaneous elliptical contact region; nine contact stations along the meshing path are used in the reported example.
When contact and thermal effects are evaluated separately, the uncoupled NLTCA gives loaded contact deformations on the order of \(10^{-3}\) mm: the pinion reaches 0.004645 mm at station 9 and a minimum of 0.0002905 mm at station 5, while the gear reaches 0.0055463 mm at station 1 and a minimum of 0.0021249 mm at station 8. These values describe contact deformation only; they do not include the carburizing–meshing coupling and therefore are not the final heat-treatment distortion.
Carburizing Model for 9310 Steel
The aerospace gear material in the example is 9310 steel. Among the elements listed for the grade, Ni is the largest at 3.13 wt% and P the smallest at 0.008 wt%. In the carburizing atmosphere, the main surface reactions are
2\mathrm{CO}+2\mathrm{H}_2\rightleftharpoons \mathrm{CH}_4+\mathrm{CO}_2,
\qquad
\mathrm{CO}+3\mathrm{H}_2\rightleftharpoons \mathrm{CH}_4+\mathrm{H}_2\mathrm{O}.
$$
Active carbon atoms are absorbed at the surface and diffuse inward. The non-steady diffusion problem is described by Fick’s second law with a concentration- and temperature-dependent diffusion coefficient:
\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[-B(0.4-C)].
$$
Here, \(C\) is the carbon mass fraction, \(t\) is carburizing time, and \(x_i\) is the diffusion direction. The surface carbon flux is driven by the difference between the furnace carbon potential \(C_e\) and the surface carbon content \(C_s\), with a transfer coefficient \(\beta=\beta_0\exp(-E/RT)\). The parameters used in the example are given below.
| Parameter | Value |
|---|---|
| \(D_{0.4}\), diffusion coefficient at 0.4 carbon mass fraction | 25.5 mm2/s |
| \(Q\), carbon diffusion activation energy | 14 kJ/mol |
| \(B\), diffusion–concentration constant | 0.8 |
| \(\beta_0\), carbon transfer constant | 0.00347 mm/s |
| \(E\), activation energy for surface transfer | 34 kJ/mol |
The carbon field creates a graded case rather than a uniform surface layer. In the coupled model, the thermal deformation associated with this case is carried by the same finite elements used in the contact solution, so the local stiffness and contact geometry are not treated as independent inputs.
Coupled Semi-Analytical Solution
The coupled problem is expressed as a functional map from machine settings and contact position to flank displacement. At each instantaneous contact ellipse center, the force is split into a thermal contribution \(\mathbf{F}_e\) and a contact contribution \(\mathbf{F}\):
\mathbf{F}_{\Omega}=\mathbf{F}_{e}+\mathbf{F},
\qquad
\mathbf{u}:=\mathbf{u}(\phi_i,\theta,t)
=\mathbf{C}_{\Omega}\mathbf{F}_{\Omega}
=\Theta_{\Omega}\mathbf{e}
+\mathbf{I}_{\Omega}\mathbf{Y}
-\mathbf{D}.
$$
Here, \(\mathbf{u}\) is the flank deformation, \(\mathbf{C}_{\Omega}\) is the compliance of the local contact region, \(\Theta_{\Omega}\mathbf{e}\) represents the rigid-body approach, \(\mathbf{I}_{\Omega}\mathbf{Y}\) accounts for slack, and \(\mathbf{D}\) is the initial gap. The corresponding weak form balances the coupled contact force against the external load:
\int_{\Phi}\left(\mathbf{e}^{T}\mathbf{F}_{\Omega}-P\right)\delta\mathbf{u}\,d\Phi=0.
$$
The semi-analytical character comes from the influence-coefficient representation of the contact and thermal response. The local contact inequalities are handled by linear programming with the improved simplex method, while the global nonlinear equilibrium is updated by Newton–Raphson iteration. Once the instantaneous contact centers are solved, the same functional relation is evaluated over the full tooth flank. This is the step that makes the method a full-flank distortion predictor rather than a single-point contact calculation.
Predicted Distortion of the Pinion and Gear
All deformations below refer to the tooth flank after cooling to room temperature. The R–Z–U frame is defined with R radial, Z axial, and U along the tooth normal. Positive and negative signs indicate the direction of displacement; total deformation is the vector magnitude of the three components, not an arithmetic sum. Contact stations 1 to 9 run from the start to the end of meshing.
| Part | Component | Station 1 | Station 9 | Additional extremum | Observed trend |
|---|---|---|---|---|---|
| Pinion | R | −0.29639 | −0.41332 | — | Magnitude increases along the path |
| Pinion | Z | −0.21098 | −0.18488 | — | Magnitude decreases |
| Pinion | U | −0.02134 | +0.023334 | −0.000882 at station 5 | Sign reversal along the path |
| Pinion | Total | 0.364532 | 0.453224 | Minimum at station 1; maximum at station 9 | Increases |
| Gear | R | −0.44409 | −0.52285 | — | Magnitude increases along the contact path |
| Gear | U | +0.040989 | +0.020768 | — | Decreases |
| Gear | Total | 0.462923 | 0.532417 | Minimum at station 1; maximum at station 9 | Increases |
Pinion. In the R direction, the contact-center deformation is negative and grows in magnitude from −0.29639 mm to −0.41332 mm. Across the whole gear structure, the R distortion increases from the tooth small end to the large end. In the Z direction, the behavior is opposite in trend: the center deformation becomes less negative, from −0.21098 mm to −0.18488 mm. The shaft, hub, and web show much smaller Z distortion than the tooth flank.
The pinion U component is the most distinctive because it changes sign. At station 1 it is −0.02134 mm; at station 5 it is −0.000882 mm; at station 9 it is +0.023334 mm. A single scalar distortion value would miss this reversal. Over the full pinion, the larger U distortions are near the concave flank at the small and large ends, while the shaft, web, and hub remain smaller. The total pinion distortion rises from 0.364532 mm at station 1 to 0.453224 mm at station 9. The largest total values occur near the tooth tip at the large end. On the shaft, the middle region has the smallest distortion and the lower region the largest; the shaft–hub connection varies little, while the tooth flank changes strongly from the large end to the small end.
Gear. The gear R distortion is larger in magnitude than the pinion R distortion. It changes from −0.44409 mm at station 1 to −0.52285 mm at station 9, and the pattern follows the contact path. Among the structural regions, the shaft is the least distorted in R. In Z, the gear differs from the pinion: the center deformation decreases from station 1 to station 9, and along the large-to-small-end direction the whole structure’s Z distortion increases and changes from negative to positive. From tooth tip to root, the flank Z distortion increases, and the web, hub, and shaft show relatively large positive values.
The gear U distortion decreases from +0.040989 mm at station 1 to +0.020768 mm at station 9. From the upper shaft to the lower shaft, the U distortion increases; from the large end to the small end, the flank distortion also increases, consistent with the loaded contact footprint. The total gear distortion rises from 0.462923 mm to 0.532417 mm. It increases downward along the shaft and from the small end to the large end along the flank, with large values near the tooth tip at the large end.
Accuracy and Scope
The U-direction values at the instantaneous contact ellipse centers were compared with an earlier adaptive data-driven prediction. The maximum relative error for the pinion was 0.87% at contact station 2; the maximum for the gear was 0.80% at station 8. All contact stations were within 1%. This supports the consistency of the coupling implementation and the contact solution.
The comparison is numerical, not a full metrological validation. It covers the U direction at the contact centers of one aerospace gear pair and one material condition. Absolute accuracy over the entire flank would require experimental distortion measurement after carburizing and quenching, including the actual cooling path and fixturing. Extrapolation to other modules, case depths, quenching schedules, or structural stiffnesses should therefore be checked rather than assumed.
Manufacturing Implications
The main value of the semi-analytical route is not a single distortion number but a traceable map from machine settings to coupled distortion. That map can be used to screen machine-setting changes before cutting trials, to estimate how a flank correction propagates into the loaded contact pattern, and to decide where hard-finishing allowance should be placed. Because the method keeps the contact and thermal fields on the same mesh, it also exposes why a correction based only on the cold, unloaded flank can miss the region that matters most under load.
The same contact solution can be extended to contact stress, contact pressure, and transmission error, since these quantities are derived from the same loaded equilibrium. For production planning of high-precision spiral bevel gears, the useful role of this kind of model is to reduce the number of full commercial finite-element trials while retaining the coupling between carburizing distortion and meshing contact. Its limitation is equally clear: the prediction is only as good as the machine-kinematic model, the carburizing parameters, and the assumed thermal boundary conditions, so it should be used as a physically grounded screening and compensation tool rather than a substitute for final experimental verification.