Script-Driven ABAQUS Workflow for Loaded Tooth Contact Analysis of Spiral Bevel Gears

Spiral bevel and hypoid gearsets carry torque between intersecting or offset axes through flanks that are only locally conjugate. The mismatch is deliberate: cutter geometry, machine settings and tooth depth are chosen so that the unloaded contact stays inside a small ellipse near the mid-face and expands into a controlled pattern as torque rises. Because the contact pattern, the load sharing between teeth and the transmission error are all governed by micrometre-scale flank geometry, loaded tooth contact analysis (LTCA) by finite elements is only as trustworthy as the geometric description fed into the solver. Mesh density and solver settings matter, but a flank that is displaced by a few micrometres will move the contact ellipse and distort the predicted transmission error long before numerical error becomes visible.

Where the conventional modelling chain loses accuracy

The usual route to an LTCA model has four stages: flank coordinates are computed from the cutting theory; the point data are imported into a CAD kernel and lofted into a solid; the gear pair is assembled in its mounting position; and the assembly is transferred to a finite element package for meshing and pre-processing. Each hand-off introduces its own tolerance. Spline fitting inside the CAD kernel approximates the computed points, and neutral-format exchange between CAD and FE environments resamples the surfaces again. The result is that the solid reaching the mesher is not exactly the surface defined by the cutting kinematics.

The penalty is not limited to geometric error. Studies of parameter sensitivity — flank modification, material data, boundary conditions, mounting misalignment — require the model to be rebuilt every time an input changes. Repeating the four-stage chain for each variant turns a parametric study into a manual routine, and it also prevents the analysis from being embedded in a digital manufacturing loop where flank corrections must be evaluated repeatedly against measured tooth surfaces.

An alternative is to move the geometry generation into the finite element kernel itself. The ABAQUS scripting interface exposes the full set of CAE operations to Python and adds several hundred library functions beyond the standard language, so that part creation, assembly, meshing, interaction definition and job submission can all be driven by a single script. Building the flank from the cutting equations inside that kernel removes the neutral-format exchange entirely and makes the model regenerable from a parameter list.

Spiral bevel gear pair with curved teeth and offset pinion

Cutting kinematics as the geometric source

The flank is defined by the generating process rather than by a drafting operation. A cutter head mounted on a cradle represents one tooth of an imaginary generating gear; the workpiece and the generating gear rotate at a fixed ratio while the blade sweeps through the tooth space, and the envelope of the blade positions forms the flank. The wheel of the pair discussed here is cut by the formate method, in which no generating roll is used and the flank shape is set directly by the cutter profile and the machine settings. The pinion is cut by the duplex helical method, which superimposes an additional helical motion on the cradle rotation.

For the pinion, the rotation of the workpiece follows the cutter-head rotation through the blade-group ratio, and the cradle rotation adds the generating motion:

$$\psi_1 = R_b\,\theta,\qquad R_b=\frac{N_t}{N_w}$$

$$\psi = R_a\,\phi + \psi_1,\qquad R_a=\frac{N_p}{N_w}$$

where θ is the cutter rotation angle, φ the cradle adjustment angle, Nt the number of blade groups, Nw the workpiece tooth number and Np the number of teeth of the generating gear. The helical motion of the duplex process is introduced through a vertical slide whose position varies linearly with the generating rotation:

$$X_b = X_{b0} – H\,R_a\,\psi$$

with Xb0 the vertical position without helical motion and H the helical motion coefficient, expressed in mm·rad−1. The blade profile is described in its own coordinate system by Gaussian coordinates s1 and θ1, with α1 the blade profile angle, r1 the tip radius, s01 the distance from the reference point of the root blade and R1 the radius of the circular blade profile. Upper and lower signs in the profile expressions refer to the outer and inner blades respectively.

The generating surface and its unit normal are carried from the cutter frame into the machine frame by a chain of homogeneous transformations, one for each machine setting: radial cutter position and cutter swivel, cutter tilt, vertical and horizontal offsets, machine root angle, and the cradle rotation itself:

$$\mathbf{r}_1(\phi,s_1,\theta_1)=M_{st}(\phi)\,\mathbf{r}_t(s_1,\theta_1),\qquad \mathbf{n}_1(\phi,\theta_1)=T_{st}(\phi)\,\mathbf{n}_t(\theta_1)$$

$$M_{st}(\phi)=M_{sk}M_{km}M_{mj}M_{ji}M_{it}$$

Here Tst is the leading third-order submatrix of Mst, which transforms normals without translation. The contact between the blade family and the workpiece is then fixed by the meshing equation, which states that the velocity of the generating surface relative to the workpiece has no component along the common normal:

$$f_1=\mathbf{n}_1\cdot\frac{\partial \mathbf{r}_1}{\partial \phi}=0$$

Solving f1 = 0 together with the surface parameterisation gives, for each cradle position, the line of contact and hence the flank point coordinates. The fillet between flank and root is generated by the tip arc of the blade. With tip arc radius ρ1 and parameter λ1, the fillet surface and its normal follow the same transformation:

$$\mathbf{r}_{c1}(\lambda_1,\theta_1,\phi)=M_{st}(\phi)\,\mathbf{r}_{ct}(\lambda_1,\theta_1),\qquad \mathbf{n}_{c1}(\lambda_1,\theta_1,\phi)=T_{st}(\phi)\,\mathbf{n}_{ct}(\lambda_1,\theta_1)$$

For the formate wheel the flank is obtained without a generating roll, as a direct offset of the blade tip path in the machine frame:

$$\mathbf{r}_2=\mathbf{r}_m-s_2\,\mathbf{t}_2-X_G\,\mathbf{p}_2-E_{m2}\begin{bmatrix}0\\1\\0\end{bmatrix}-X_{b2}\begin{bmatrix}0\\0\\1\end{bmatrix},\qquad \mathbf{n}_2=R(Z,q_2)\,R(Z,-\theta_2)\,\mathbf{n}_0$$

where s2 is the blade profile parameter, t2 the unit tangent at the considered point, p2 the gear axis direction, XG the horizontal setting, Em2 the vertical offset and Xb2 the sliding base.

From flank equations to a meshable solid

The equations define an extended surface, only part of which belongs to the actual tooth. The first practical step is therefore to bound the parameter domain: the flank is divided into m intervals along the face width and n intervals along the tooth height, producing an m × n grid whose physical extent is limited by the blank geometry. Each grid point is first expressed in a local frame tied to the pitch cone and then rotated into the design coordinate system, using the pitch cone angle δm, the cone distance Am and the distance Zp from the cone apex to the crossing point:

$$X_1 = X\cos\delta_m – Y\sin\delta_m + A_m\cos\delta_m – Z_p,\qquad Y_1 = X\sin\delta_m + Y\cos\delta_m + A_m\sin\delta_m$$

These coordinates are linked to the computed position vector r1 and the pinion axis direction p by

$$X_1=-\mathbf{r}_1\cdot\mathbf{p},\qquad Y_1=\left|\mathbf{r}_1\times\mathbf{p}\right|$$

The numerical solution proceeds iteratively. Starting from initial values of cradle angle and cutter phase, the meshing equation is solved by Newton iteration in a rotational projection scheme, and the converged surface parameters are converted into three-dimensional points. The loop continues until m × n points have been obtained for each flank. The point sets are then lofted into a family of curves, the curves into the flank surface, and the closed tooth slot is cut from the gear blank by a Boolean subtraction. Repeating the slot around the axis produces the full gear body, which is the input to the finite element stage.

Assembling the finite element model through the script interface

The pair used to demonstrate the workflow is a 7 × 43 spiral bevel gearset with a 25.4 mm pinion offset and a 90° shaft angle, cut with the settings listed below. The geometric data describe the blank and the mounting position; the machining data describe the cutter head and the machine motions that produce the flanks.

Table 1. Blank and mounting geometry of the 7 × 43 gear pair
Parameter Pinion Wheel
Number of teeth 7 43
Module (mm) 6.860
Face width (mm) 43.512 40.000
Pinion offset (mm) 25.400 —
Shaft angle (°) 90.000
Mid-point spiral angle (°) 49.000 38.770
Hand Left Right
Outer cone distance (mm) 151.744 150.465
Outer pitch diameter (mm) 58.945 294.980
Pitch cone angle (°) 11.200 78.588
Face cone angle (°) 14.299 78.932
Root cone angle (°) 10.862 75.433
Table 2. Cutter-head and machine settings for the two members
Parameter Pinion (duplex helical) Wheel (formate)
Radial cutter position (mm) 102.805 117.492
Cutter tilt (°) 3.582 —
Cutter swivel (°) 196.300 —
Vertical offset (mm) 23.114 —
Machine root angle (°) 11.427 75.433
Horizontal setting (mm) −1.988 9.652
Sliding base (mm) 4.813 —
Roll ratio 6.070 —
Angular cutter position (°) 65.581 70.977
Helical motion coefficient (mm·rad−1) −1.475 —
Cutter radius (mm) 111.882 114.300
Inner / outer blade profile angle (°) 25.000 / 20.000 25.000 / 20.000
Point width (mm) 2.536 3.300

The script instantiates the two generated bodies, positions them at the standard mounting distance and defines the contact pair, with the pinion flank as master surface. To keep the model within available computing resources, the wheel is represented by a segment of eight teeth; the pinion is modelled over the corresponding sector. Both members are meshed with C3D8R brick elements. Coupling constraints tie the flank surfaces to reference points on the respective axes, so that rotation and torque can be applied through a single node. The wheel reference point carries the resisting torque, the pinion reference point the rotational boundary condition, and the analysis steps are defined to output contact pressure, contact area and the rotation of both members.

Accuracy of scripted geometry compared with imported CAD

The advantage of generating the solid inside the solver becomes measurable when the same point cloud is processed by both routes: once through a solid modeller and a neutral-format import, once through the direct script. Three comparisons are informative. First, the surface mesh of the imported model is irregular and contains distorted elements, and merging the meshes of the circularly arrayed teeth at a given tolerance leaves gaps; the scripted model meshes regularly and merges cleanly at the same tolerance. Second, when a single flank from each route is overlaid, the four corner vertices coincide and the discrepancy only appears in the interior of the surface, at edge midpoints, where corresponding points differ by up to about 91 nm. Third, the assembled single-tooth solids show deviations of roughly 0.07–1.93 µm at key points. The error grows with geometric complexity, which is precisely where the imported route is weakest.

Table 3. Comparison of the two modelling routes from an identical point set
Check Script-generated model CAD model imported into the solver
Surface mesh quality Regular; elements merge within tolerance Irregular, distorted elements; merging unsatisfactory at equal tolerance
Flank interior deviation Reference Up to ≈91 nm at edge midpoints
Single-tooth key-point deviation Reference ≈0.07–1.93 µm
Corner vertices Reference Coincident with reference

These magnitudes are small in absolute terms but not in the context of local conjugate gearing, where the designed mismatch between the two flanks is itself of the order of tens of micrometres across the face. A deviation that is systematic rather than random, and that is largest in the interior of the flank, biases the predicted contact ellipse instead of merely adding scatter to it.

Contact pattern and transmission error under load

With the model assembled, static analyses were run at wheel torques of 10, 300, 600 and 900 N·m on both the drive and coast sides. The predicted contact pattern starts as a narrow, elongated ellipse at 10 N·m, consistent with a locally conjugate pair whose elastic deformation is still small. As torque increases, the ellipse lengthens and spreads across the face, and at the highest load the contact approaches the tooth edges — the expected consequence of elastic deflection enlarging the instantaneous contact ellipse beyond the region that the mismatch design had reserved for it. The peak contact stress remains in the central part of the flank at all four load levels, and no abnormal contact concentration appears in the computed patterns.

Transmission error was extracted from the rotation histories of the two members:

$$\varepsilon=(\varphi_2-\varphi_{20})-\frac{Z_2}{Z_1}(\varphi_1-\varphi_{10})$$

where φ1 and φ2 are the instantaneous rotation angles of pinion and wheel, the subscript 0 denotes the initial position, and Z1, Z2 are the tooth numbers. At 10 N·m the finite element curve lies slightly below the curve obtained from the numerical tooth contact analysis, with comparable amplitude and period. The vertical offset is the elastic compliance introduced by the applied torque. The horizontal phase difference has a different origin: the finite element assembly does not begin meshing exactly at the design reference point, so its initial instantaneous ratio differs from the nominal one, and the numerical curve was shifted along the abscissa to align the two for comparison.

Increasing the torque changes the curve in three ways. The whole curve moves downward, because the elastic rotation of the driven member grows with load. The amplitude decreases, because at high torque the load is shared by more than one tooth pair, so the mesh entry and exit impacts are softened and the motion becomes smoother. The position of the maximum migrates toward the reference point at mid-face, since a locally conjugate pair has a low contact ratio at that point and therefore deforms most there under load. At equal torque, the coast side shows a larger transmission error amplitude than the drive side, which indicates smoother running on the drive side under load.

The agreement at light load is the main validation of the modelling route: a model whose flank geometry had been corrupted by format exchange would not reproduce the phase and period of the numerical transmission error curve, nor would it show the correct direction of the load-induced offset. The load-dependent trends — growth of the contact ellipse, downward shift and amplitude reduction of the transmission error, migration of its maximum — are consistent with the mechanics of a locally conjugate pair rather than with numerical artefacts.

What the workflow enables, and where it stops

The practical value of generating the model inside the solver is regenerative rather than one-off. Once the flank equations, the blank parameters and the machine settings are exposed as script variables, a change in cutter profile, roll ratio or mounting offset requires editing those variables and rerunning the script; the mesh, contact definitions and boundary conditions are rebuilt automatically. This makes sensitivity studies of flank modification, material parameters and boundary conditions feasible at a fraction of the manual effort, and it opens the possibility of closing the loop between measured flank deviations and corrected machine settings, which is the point at which LTCA becomes a manufacturing tool rather than a verification exercise.

The limits of the present implementation should be stated with equal clarity. The wheel is represented by an eight-tooth sector, which is adequate for contact and transmission error over the analysed arc but suppresses any whole-body dynamic response of the complete wheel. The analyses are static, with torque applied through coupled reference points, so the transmission error reported here is a quasi-static quantity driven by tooth compliance and the imposed load, not a prediction of gear whine under running speed. Accuracy is controlled by the density of the m × n surface grid and by the tolerance used when lofting points into surfaces; both must be refined until the contact pattern and transmission error stop changing. Within these bounds, the direct scripted route removes the dominant geometric error source of the conventional CAD-to-FEA chain and produces an LTCA model whose fidelity is limited by the cutting theory and the discretisation, not by data exchange.

Scroll to Top