Ani Ural, Paul A. Wawrzynek, and Anthony R. Ingraffea

Cornell University, Ithaca, New York

Simulating Fatigue Crack Growth

in Spiral Bevel Pinion

NASA/CR—2003-212529

August 2003

ARL–CR–0531

U.S. ARMY

RESEARCH LABORATORY
The NASA STI Program Office . . . in Profile

Since its founding, NASA has been dedicated to

the advancement of aeronautics and space

science. The NASA Scientific and Technical

Information (STI) Program Office plays a key part

in helping NASA maintain this important role.

The NASA STI Program Office is operated by

Langley Research Center, the Lead Center for

NASA’s scientific and technical information. The

NASA STI Program Office provides access to the

NASA STI Database, the largest collection of

aeronautical and space science STI in the world.

The Program Office is also NASA’s institutional

mechanism for disseminating the results of its

research and development activities. These results

are published by NASA in the NASA STI Report

Series, which includes the following report types:

• TECHNICAL PUBLICATION. Reports of

completed research or a major significant

phase of research that present the results of

NASA programs and include extensive data

or theoretical analysis. Includes compilations

of significant scientific and technical data and

information deemed to be of continuing

reference value. NASA’s counterpart of peerreviewed formal professional papers but

has less stringent limitations on manuscript

length and extent of graphic presentations.

• TECHNICAL MEMORANDUM. Scientific

and technical findings that are preliminary or

of specialized interest, e.g., quick release

reports, working papers, and bibliographies

that contain minimal annotation. Does not

contain extensive analysis.

• CONTRACTOR REPORT. Scientific and

technical findings by NASA-sponsored

contractors and grantees.

• CONFERENCE PUBLICATION. Collected

papers from scientific and technical

conferences, symposia, seminars, or other

meetings sponsored or cosponsored by

NASA.

• SPECIAL PUBLICATION. Scientific,

technical, or historical information from

NASA programs, projects, and missions,

often concerned with subjects having

substantial public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientific

and technical material pertinent to NASA’s

mission.

Specialized services that complement the STI

Program Office’s diverse offerings include

creating custom thesauri, building customized

databases, organizing and publishing research

results . . . even providing videos.

For more information about the NASA STI

Program Office, see the following:

• Access the NASA STI Program Home Page

at http://www.sti.nasa.gov

• E-mail your question via the Internet to

help@sti.nasa.gov

• Fax your question to the NASA Access

Help Desk at 301–621–0134

• Telephone the NASA Access Help Desk at

301–621–0390

• Write to:

 NASA Access Help Desk

 NASA Center for AeroSpace Information

 7121 Standard Drive

 Hanover, MD 21076
Ani Ural, Paul A. Wawrzynek, and Anthony R. Ingraffea

Cornell University, Ithaca, New York

Simulating Fatigue Crack Growth

in Spiral Bevel Pinion

NASA/CR—2003-212529

August 2003

National Aeronautics and

Space Administration

Glenn Research Center

ARL–CR–0531

U.S. ARMY

RESEARCH LABORATORY

Prepared under Contract NAS3–1993
Available from

NASA Center for Aerospace Information

7121 Standard Drive

Hanover, MD 21076

National Technical Information Service

5285 Port Royal Road

Springfield, VA 22100

The Propulsion and Power Program at

NASA Glenn Research Center sponsored this work.

Available electronically at http://gltrs.grc.nasa.gov
NASA/CRó2003-212529 1 

1. Introduction 

 

This project investigates computational modeling of fatigue crack growth in spiral 

bevel gears. Current work is a continuation of the previous efforts made to use the 

Boundary Element Method (BEM) to simulate tooth-bending fatigue failure in spiral 

bevel gears. This report summarizes new results predicting crack trajectory and fatigue 

life for a spiral bevel pinion using the Finite Element Method (FEM). 

Predicting crack trajectories is important in determining the failure mode of a 

gear. Cracks propagating through the rim may result in catastrophic failure, whereas the 

gear may remain intact if one tooth fails and this may allow for early detection of failure. 

Being able to predict crack trajectories is insightful for the designer. However, predicting 

growth of three-dimensional arbitrary cracks is complicated due to the difficulty of 

creating three-dimensional models, the computing power required, and absence of closedform solutions of the problem. 

Previously in this project, tooth-bending fatigue failure in spiral bevel gears was 

investigated using the BEM [1]. These analyses were significant in developing a method 

for predicting three-dimensional, non-proportional, fatigue crack growth incorporating 

moving loads. Prior to the BEM study, there were not many examples of threedimensional crack growth work on gears. Most of the gear studies in the literature are 

two-dimensional analyses. Furthermore, moving loads had not been considered in these 

analyses. 

The BEM predictions were helpful to determine areas in need of further research. 

These were: improving accuracy of simulation by switching from the boundary element 

to the finite element method, decreasing computation time by employing parallel FEM 

analysis, and improving modeling of contact tooth loads. 

Substantial improvement in computation time has been achieved by using parallel 

FEM technologies on PC-clusters. This makes it possible to simulate crack growth in 

large models, like the current gear problem, in a few minutes. It is also possible to obtain 

accurate values of stress intensity factors (SIF) using the FEM and elastic, equivalent 

domain J-integral approaches. A new technique was developed for modeling moving 

tooth loads, which decreases the computation time while introducing a more efficient 

procedure for applying contact loads. 

Another focus of this project was performing three-dimensional contact analysis 

of a spiral bevel gear set incorporating cracks. These analyses were significant in 

determining the influence of change of tooth flexibility due to crack growth on the 

magnitude and location of contact loads. This is an important concern since change in 

contact loads might lead to differences in SIFs and therefore result in alteration of the 

crack trajectory. 

Contact analyses performed in this report showed the expected trend of 

decreasing tooth loads carried by the cracked tooth with increasing crack length. 

Decrease in tooth loads lead to differences between SIFs extracted from finite element 

contact analysis and finite element analysis with Hertz contact loads. This effect became 

more pronounced as the crack grew. 
NASA/CRó2003-212529 2 

2. Three-dimensional FEM Simulations of Fatigue Crack Growth in Spiral Bevel 

Gears 

 The main objective of the current work is to demonstrate that we can improve 

accuracy and efficiency of crack growth simulations in spiral bevel pinion through FEM 

computations along with improved modeling of tooth loads. In the current study, as in 

previous BEM simulations, LEFM theories are used for fatigue crack growth predictions 

in gears. Fatigue crack growth rates were determined using a modified Paris model 

accounting for crack closure. Crack trajectory predictions are made using the maximum 

circumferential stress theory. Moving loads were taken into account by discretizing the 

loading into 15 steps that created a non-proportional load history at the tooth root. 

Accuracy of computations can be improved using the FEM that allows a better 

SIF calculation method, the equivalent domain J-integral method. Solution time for each 

load step can be decreased substantially by employing parallel FEM analysis. Another 

improvement in the accuracy of the results can be obtained by developing a better 

representation of moving contact loads. This is accomplished by approximating the 

moving contact ellipse by interpolation using the shape functions on the surfaces of the 

loaded elements. In the earlier boundary element computations, contact ellipses were a 

part of the geometry models. Since the ellipses overlapped, four different geometry 

models were needed for each crack configuration. In the current approach, contact loads 

are no longer a part of the geometry model. This approach brings a substantial decrease in 

computation time since only one geometry model for each crack configuration is now 

needed. 

 The method used in FEM simulations of spiral bevel gear is composed of the 

following steps: 

1. Create initial geometry model of the spiral bevel pinion using OSM. 

2. Specify boundary conditions and material properties in FRANC3D. 

3. Initiate crack in FRANC3D. 

4. Create a surface mesh composed of triangular elements in FRANC3D. 

5. Create a 3D FEM mesh of the model composed of tetrahedra by J-Mesh using a 

geometry description file written out from FRANC3D. 

6. Calculate the magnitude and location of the nodal contact loads for each load step 

on the loaded elements using the mesh file produced by J-Mesh. 

7. Run FEM analysis on a parallel PC-cluster for 15 discrete load steps. Using the 

displacement results, calculate SIF values at each discrete crack front point for 

each load step. 

8. Calculate the amount of crack extension at each discrete crack front point as a 

result of 15 steps of non-proportional moving load. 

9. Determine new crack front by piecewise least squares fit of the propagated points 

corresponding to each discrete crack front point in FRANC3D. 

10. Remesh surface locally and repeat steps 5 to 10. 
NASA/CRó2003-212529 3 

2.1 Modeling of Moving Tooth Loads 

 Contact in a spiral bevel gear occurs in three-dimensions following a path along 

the tooth surface starting from the fillet of the toe to the top of the heel. In analyzing the 

spiral bevel gear, it is assumed that the contact between the gear and pinion follows Hertz 

theory of elastic contact. Hertzian contact holds as long as the significant dimensions of 

the contact area are small compared to the dimensions of each body and to the relative 

radii of curvature of surfaces [2]. Hertz theory assumes that surfaces are continuous and 

non-conforming, strains are small, each solid can be considered as an elastic half-space 

and surfaces are frictionless. 

Under these assumptions, only normal pressure acts between two bodies due to 

Hertz contact producing normal displacements of the surfaces. Since the contact area is 

assumed to be elliptical the region is defined by, 

1 2

2

2

2

+ = b

y

a

x (1) 

with semi-ellipsoidal pressure distribution of the form, 

1/ 2

2

2

2

2

0 1 













 = − − b

y

a

x

p p (2) 

Total load acting on the ellipse and average load are then calculated as, 

0 3

2

p

ab

p

π = , 0 3

2

p p avg = (3) 

In the above equations p0 denotes the Hertz stress which is the maximum stress, a and b

are half of major and minor axes of the ellipse, respectively, Figure 2.1. 

 

Figure 2.1: Representation of Hertz stress over an elliptical area. 

 In order to determine the contact points on the gear tooth, the method described in 

[3] is used to map the three-dimensional contact path into a two dimensional space. The 

following equations were used for the transformation: 

2 1/ 2

, ,1,2

2

, ,1,2 , ,1,2 ( ) RN M = X N M + YN M

η = − α + − α

ξ = − α + − α

( ) cos ( )sin

( ) cos ( )sin

N N M N M

N N M N M

R R Z Z

Z Z R R

 (4) 

b

Z Z

a

R R

2

cos

2

sin 1 2 1 − 2 α = − α =

p0 

pavg 

p(x,y) 

2a 
NASA/CRó2003-212529 4 

XM,YM,ZM = coordinates of mean contact point (corresponds to the contact ellipse centers)

XN,YN,ZN = coordinates of an arbitrary contact point 

X1,2,Y1,2,Z1,2 = coordinates of points corresponding to 1 and 2 in Figure 2.2 in three 

dimensions. 

1 2

2

2

2

= η + ξ

a b 

 

Figure 2.2: Sketch of transformed contact ellipse with mean contact point M and an 

arbitrary contact point N. 

 

 An algorithm was developed to find the nodes within the contact ellipses and to 

calculate the applied load on these nodes. This algorithm: 

1. reads the nodal information from the mesh and geometry files created by 3D FEM 

meshing program, 

2. retrieves the geometry information of the n th contact ellipse, 

3. transforms the nodal points defined in 3D coordinates into contact plane, 

4. checks if the mesh nodes are within the contact area bounded by the n th ellipse 

(n = 1Ö15), eg. Figure 2.3, 

5. calculates nodal contact loads if the node is within the contact area, and 

6. writes out the information of nodes found by this search and corresponding loads. 

Figure 2.3: Nodal points where contact loads are applied for load step 11. 

2.2 Parallel FEM 

 Understanding fracture phenomena in materials is crucial for science and 

engineering applications such as crack propagation in spiral bevel gears. However, crack 

growth simulations are complicated and computationally expensive in realistic structures. 

In order to solve realistic 3D problems, it is necessary to develop fast and accurate 

computational simulation tools. 

2a 

ξ

2b M η

1 N 2 
NASA/CRó2003-212529 5 

 The BEM is an alternative to the FEM for solving large problems. The BEM is 

advantageous since it only requires a surface mesh and has fewer equations to be solved. 

However, it entails a fully dense, nonsymmetric system of equations, and it is difficult to 

extract accurate near-crack-front stresses and displacements. In this respect, it is more 

desirable to use FEM. This brings the challenge of implementing fast and robust parallel 

sparse solvers. 

A parallel PC-cluster in Cornell Theory Center was used for the development of a 

parallel FEM solver for crack growth simulations. The hardware consists of 32 nodes of a 

64 node cluster of Dell PowerEdge 1650 servers with a Giganet interconnect running 

Microsoft Windows 2000 Advanced Server. Each node contains 2 Intel Pentium III 

processors clocked at 1 GHz and 2 GB of RAM. Only one processor per node is used for 

the simulations. 

Parallel FEM simulation of crack propagation makes it possible to simulate large 

problems that were impossible or impractical to solve previously. The objective in 

developing such a system is to be able to solve a real-life problem with multiple crack 

propagation steps with 1,000,000 degrees of freedom in an hour. Such a problem with 

10,000 degrees of freedom takes about 100 hours using a state-of-the-art single processor 

workstation [4]. 

Crack propagation in a spiral bevel gear stands as a large, realistic engineering 

problem that incorporates the computational difficulties mentioned above. In that respect, 

the development of parallel FEM simulations is a significant step in achieving the goal of 

developing a practical and accurate numerical design tool for gears. 

2.3 SIF Calculations

 Complex problems, such as the current gear problem, require computation of 

mixed mode SIFs due to the combined loading that they receive. In the literature, several 

approaches for computing SIFs can be found. Among these methods, the displacement 

correlation method [5-8] and equivalent domain J-integral methods [9-14] are widely 

used in numerical computations. 

 In the previous BEM simulations of the spiral bevel gear, the displacement 

correlation method was used for calculating SIFs. This method only requires 

displacement information on the boundary near the crack front and is computationally 

efficient. Since BEM solves for displacement information only on the boundaries, this 

method is suitable for the BEM analysis. The displacement correlation method makes use 

of numerical solutions for displacements in combination with the analytic solutions of 

SIFs. This method is only as accurate as the computed displacements at the correlation 

points. 

The equivalent J-integral method is adaptable to complex problems and presents 

much more flexibility in the type of applications for which it can be used. We used the 

equivalent domain J-integral approach instead of displacement correlation in the FEM 

simulations of the spiral bevel gear. 

 J-integral is the energy release expression associated with crack growth defined as 

an integral along a contour that encloses a crack tip. J-integral is defined as [15], 
NASA/CRó2003-212529 6 

Γ

∂

∂ = − σ ∫

Γ

n d

x

u J Wn j

i

ij ( )

1

1 (5) 

where W is the strain energy density, nj is the unit outward normal to the contour, σij is 

the traction exerted on the body enclosed by Γ and crack surface. This integral is 

evaluated along an arbitrary contour enclosing the crack tip and is path independent. 

ΓT = Γ′ + Γ1 + Γ2 + Γ

Figure 2.4: Paths of J-integral calculations at a crack tip. 

 

Numerical evaluation of the J-integral involves the use of equivalent domain 

integral approach. In this approach, the integration along the contour is replaced by an 

integral over a finite size domain. The expression for J-integral then becomes, 

∫

Γ ∂

∂ = δ − σ

T

n qdA

x

u J W j

k

i

k kj ij [ ] k=1,2. (6) 

A continuous function q(x1,x2) is introduced in the above equation in order to replace the 

integration contour Γ by ΓT. Function q has a unit value at the crack tip and it is zero 

along Γ', Γ1 and Γ2. The variation of this function within the domain is arbitrary. Using 

divergence theorem equation 6 takes the form, 

∫ ∫

Γ Γ ∂

∂

σ

∂

∂ − ∂

∂ − ∂

∂

∂

∂ − σ

∂

∂ =

T T

qdA

x

u

x x

W dA

x

q

x

u

x

q J W

k

i

ij

k j k j

i

ij

k

k [ ] { [ ]} (7) 

 

Note that the second term of this integral vanishes for elastic problems. SIFs are 

calculated within the vicinity of the 3D crack front by the pointwise values of J-integral. 

3D cracks can be represented by 2D plane strain fields at the crack front points. These 

fields are satisfied only within close proximity of the crack tip, however it is well known 

that finite element calculations of the field quantities are not very accurate close to the 

crack tip. Equivalent domain integral approach circumvents this difficulty by replacing 

the line integral by a finite domain integral. The relation between the J-integral and SIFs 

is given for an isotropic, linear elastic material as follows, 

Γ

Γ2

Γ1

x1

x2

Γ'

n
NASA/CRó2003-212529 7 

( ) 8

1 2 2

1 KI KII G

J +

κ + = and KI KII G

J

8

2( 1)

2

κ + = − (8) 

In the above equations J2 is not path independent. For a pure mode I or mode II problem, 

SIFs can be calculated directly from these equations however under mixed mode loading 

KI and KII cannot be calculated separately. Bui [16] proposed a mixed mode separation 

method by decomposing the elastic displacement field into symmetrical and antisymmetrical parts. This results in separation of J into two components that are given as 

follows [16], 

 

1 2 J = J + J , 2

2 1

I KI E

J − ν = , 2

2 1

II KII E

J − ν = , II

i

I

ui = ui + u (9) 

Examples of numerical investigations on mixed mode problems using this approach 

yielded highly accurate results [17,18]. 

Numerical techniques generally provide good results in an average sense over the 

domain. However, governing equations and/or boundary conditions are only 

approximately satisfied at any one point. This indicates that displacement correlation 

SIFs are not very accurate since they depend on the displacement values at a point. On 

the other hand, J-integral SIFs give better results due to considering much more of the 

results domain over which approximation errors tend to cancel. 

Two models with known analytic solution are investigated to show the higher 

accuracy of the FEM calculations. The first model is a single, internal, flat, pennyshaped crack in an infinite body. The radius of the penny-shaped crack is 0.1. This crack 

is centered in a block with dimensions 10 x 5 x 5. 

Figure 2.5: Flat, internal penny-shaped crack subjected to remote tension. 

For this problem Mode II and Mode III SIFs are zero. The analytical solution for 

the Mode I SIF is given by [19], 

π = σ

a KI 2 (10) 

σ

σ

a y 

x 

z 
NASA/CRó2003-212529 8 

Table 2.1 shows the comparison of solutions from the BEM, and FEM and the analytical 

solution. BEM results shown in Table 2.1 are calculated by the displacement correlation 

method using linear quadrilateral elements with 12 subdivisions along the radius and 

along the quarter of the crack front. FEM calculations are carried out using the equivalent 

domain J-integral method with 36,583 tetrahedra. Note that an elastic modulus of 10,000 

and a Poissonís ratio of 0.0 are assigned to the model, and a uniform, far-field tension of 

unity is applied. 

Table 2.1: Average Mode I SIF values along the crack front calculated by the BEM and 

the FEM compared with the analytical solution. 

Solution method Penny-shaped % Error

Analytical 0.357 - 

BEM 0.348 2.6 

FEM 0.352 1.4 

 

The second example is a single, internal, flat, elliptical crack in an infinite body 

with an elastic modulus of 10,000 and a Poissonís ratio of 0.0. The orientation of the 

crack is the same as in Figure 2.5. An elliptical crack with 0.1 minor and 0.2 major axis 

dimensions is embedded in a 10 x 5 x 5 block. A uniform far-field tension of unity is 

applied. The analytic solution of KI involves elliptic integrals whereas KII and KIII are 

zero. Detailed discussion of the analytical solution can be found in [19]. Comparison of 

solutions from the BEM, and FEM and the analytical solution are given in Table 2.2. 

BEM calculations are done using the displacement correlation method with linear 

quadrilateral elements with 5 subdivisions along the axes and 12 subdivisions along the 

quarter of the crack front. The FEM results are computed using J-integral method with 

13,464 tetrahedra. 

Table 2.2: Maximum and minimum Mode I SIF values along the crack front calculated 

by the BEM and the FEM compared with the analytical solution. 

Solution method Max KI % Error Min KI % Error 

Analytical 0.463 - 0.327 - 

BEM 0.435 5.9 0.304 7.0 

FEM 0.442 4.5 0.333 1.8 

 

2.4 Determining New Crack Front

 Crack extension due to one non-proportional load cycle is computed using the 

approach developed in [1]. As previously mentioned, a modified Paris model was used 

that accounts for crack closure effects. A crack is assumed to advance when its SIF is 

large enough to overcome closure and is larger than the SIF of the previous load step. 

The new crack front is determined as follows: 
NASA/CRó2003-212529 9 

1. For every point along the crack front: 

a) Calculate amount of crack growth, dai corresponding to each load step i= 

1,Ö15. 

da

K

K K da

eff

i

I

i

i I

∆

− =

−1

 

n C Keff dN

da = (∆ )

∆Keff = KI − Kop

max 

3

3

2

0 1 2

max

A A R A R A R

K

Kop = + + + for R≥ 0 

α





 





σ

π = − α + α

1/

0

2 max

0 ) 2 (0.825 0.34 0.05 ) cos( S A

0

2 max

1 (0.415 0.071 )

σ = α − α

S A

2 0 1 3 A = 1− A − A − A

2 1 A3 = A0 + A1 −

b) Calculate the angle of crack growth corresponding to each load step using 

maximum circumferential stress theory 

 















θ = ± + − ( ) 8

4

1

4

1 2 tan 1 2

i

II

i

I

i

II

i

i I

K

K

K

K

c) Calculate the final coordinates of the crack front and trajectory angles by 

approximating the contributions from each load step by a straight line. 

d) Determine number of cycles, N, necessary to have a reasonable amount of 

crack advance compared to pinionís geometry. 

2. Update geometry using the new crack front points. 

3. Repeat 1 and 2 for each crack growth step. 

2.5 Results

 Three-dimensional finite element analyses were performed for an OH-58 spiral 

bevel pinion for 54 crack growth steps, Figure 2.6. In the previous BEM-based analysis, 

only 13 growth steps were performed. One crack growth step consists of 15 static finite 

element analyses in order to simulate moving loads on the pinion tooth. After the 54th

crack growth step the analysis was stopped because it was concluded that propagating the 

crack further would not provide additional insights into the simulations and results. All 

predictions were based on the methodology outlined in section 2. 
NASA/CRó2003-212529 10 

 Figure 2.6: A spiral bevel pinion model with surface FE mesh (a) whole pinion (b) closeup view of pinion teeth. Note the location of the initial crack in the middle tooth. 

 Figures 2.7 through Figure 2.9 show Mode I, II and III stress intensity factors for 

the initial crack configuration for the first eleven load steps. Load steps 1-4 correspond to 

double tooth contact, load steps 5-11 correspond to single tooth contact and load steps 

12-15 again correspond to double tooth contact. Since a crack is not assumed to advance 

when its SIF is smaller than the SIF of the previous load step, only the first 11 steps 

practically contribute to crack growth calculations. In these figures crack front position 

corresponds to the points near the crack front at which SIF values are evaluated. In the 

initial step of crack growth simulations the crack front was discretized such that there 

were 222 points at which SIF values were calculated. This is compared to only 61 points 

previously used with the BEM model [1]. 

Figure 2.10 shows the crack trajectory predictions by the FEM on the tooth 

surface and cross section of the tooth for several crack growth steps including the initial 

and final configurations of the crack. Figure 2.11 is a close-up view of the last step of 

crack growth from the toe end of the tooth. 

Figure 2.7: Initial Mode I SIF distribution for load steps one through eleven. Note that 

orientation is from heel to toe. 

Initial 

Crack 

0

2000

4000

6000

8000

10000

12000

14000

16000

18000

0 20 40 60 80 100 120 140 160 180 200 220

crack front position

KI (psi in1/2)

Load step 1 Load step 3 Load step 5 Load step 7 Load step 9 Load step 11
NASA/CRó2003-212529 11 

Figure 2.8: Initial Mode II SIF distribution for load steps one through eleven. Note that 

orientation is from heel to toe. 

Figure 2.9: Initial Mode III SIF distribution for load steps one through eleven. Note that 

orientation is from heel to toe. 

-5000

-4000

-3000

-2000

-1000

0

1000

2000

3000

4000

5000

6000

0 20 40 60 80 100 120 140 160 180 200 220

Crack front position

KII (psi in1/2)

Load step 1 Load step 3 Load step 5 Load step 7 Load step 9 Load step 11

-7000

-6000

-5000

-4000

-3000

-2000

-1000

0

1000

2000

3000

4000

5000

6000

7000

0 20 40 60 80 100 120 140 160 180 200 220

crack front position (in)

KIII (psi in1/2)

Load step 1 Load step 3 Load step 5 Load step 7 Load step 9 Load step 11
NASA/CRó2003-212529 12 

Initial Step: Step 0: 

Step 15: 

Step30: 

Step 45: 

Final Step: Step 54: 

Figure 2.10: Crack trajectory predictions by FEM on tooth surface and tooth cross section. 

Midpoint of crack where cross 

section is taken. 
NASA/CRó2003-212529 13 

Figure 2.11: Close-up view of the crack after the last step of crack growth. 

2.6 Comparison of FEM and BEM Results

This section compares FEM and BEM results in terms of accuracy and efficiency. 

BEM results presented in this section are from reference [1]. One of the most significant 

improvements in the FEM approach over the BEM simulations is the difference in the 

solution time of the problem. This comparison is summarized in Table 2.3. 

Table 2.3: Comparison of performance between FEM and BEM simulations. 

Results FEM BEM 

Computer power 32 nodes#

PC-cluster, 1 GHz 

Single DEC-Alpha 

Workstation, 175 MHz 

Number of unknowns 

(coefficient matrix) 

~1-3 million 

(sparse, symmetric) 

6700 

(dense, unsymmetric) 

Solution time per load cycle 

(1 load cycle=15 load steps) 

5-15 hours 

(20min-1 h/load step) 

~60 hours 

(4 models) 

Typical time for one crack step ~8-18 hours*

 ~65 hoursÜ

#

Only one processor per node is used. 

*

Solution time plus 3 hours post-processing and remeshing 

Ü

Solution time plus 5 hours post-processing and remeshing 

The performance comparison in Table 2.3 clearly shows that there is a very 

significant improvement in computation time with the development of a parallel FEM 

solver. The solution time for BEM simulations was about 65 hours for one crack growth 

step whereas FEM simulations took much less time for each crack growth step even 

though the models have many more degrees of freedom. The number of unknowns for the 
NASA/CRó2003-212529 14 

FEM was about a million at the initial step and increased to more than 3 million as the 

crack front advanced. 3D volume meshing of one model ranged from 20 minutes for the 

smallest model to about an hour for the largest models. 

Improvement in the computation time for each crack growth step allowed us to 

increase solution resolution by taking smaller crack increment steps and to simulate more 

steps of propagation. Figures 2.12 to 2.14 show the comparison of crack area, depth and 

crack front length by FEM and BEM. These plots indicate that the FEM-predicted fatigue 

growth rate is marginally lower when we compare the crack area and the crack front 

length and more cycles are required to achieve the same amount of crack growth when 

compared with previous BEM-based results. Figures 2.15 to 2.17 show the SIF 

comparisons between FEM and BEM for load steps 3, 7, and 11. 

Table 2.4: Comparison of crack growth amount between FEM and BEM simulations. 

Results FEM BEM 

Number of steps 54 13 

Total number of cycles 383,000 311,000 

Final crack front length 1.55 in 1.42 in 

Final crack area 0.319 in2

 0.186 in2 

Final crack depth 0.268 in 0.188 in 

Figure 2.12: Comparison between crack area versus number of cycles predicted by BEM 

and FEM analyses. 

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0 50000 100000 150000 200000 250000 300000 350000 400000

Number of cycles

Crack area (in2

)

BEM-311,000 cycles 

FEM-383,000 cycles 
NASA/CRó2003-212529 15 

Figure 2.13: Comparison between crack depth versus number of cycles predicted by 

BEM and FEM analyses. 

Figure 2.14: Comparison of crack front length between BEM and FEM analyses. 

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

0 50000 100000 150000 200000 250000 300000 350000 400000

Number of cycles

Crack front length (in)

BEM-311,000 cycles 

FEM-383,000 cycles 

0

0.05

0.1

0.15

0.2

0.25

0.3

0 50000 100000 150000 200000 250000 300000 350000 400000

Number of cycles

crack depth (in)

BEM-311,000 cycles 

FEM-383,000 cycles 
NASA/CRó2003-212529 16 

Figure 2.15: Comparison of initial Mode I, II and III SIFs computed by the FEM and the 

BEM for load step 3. 

Figure 2.16: Comparison of initial Mode I, II and III SIFs computed by the FEM and the 

BEM for load step 7. 

-7000

-2000

3000

8000

13000

18000

0 0.04 0.08 0.12 0.16 0.2

Crack front position (in)

K (psi in1/2)

Mode I,FEM Mode I,BEM ModeII,FEM

ModeII,BEM Mode III,FEM Mode III,BEM

Load Step 3

-7000

-2000

3000

8000

13000

18000

0 0.04 0.08 0.12 0.16 0.2

Crack front position (in)

K (psi in1/2)

Mode I,FEM Mode I,BEM Mode II,FEM

Mode II,BEM Mode III,FEM Mode III,BEM

Load Step 7
NASA/CRó2003-212529 17 

Figure 2.17: Comparison of initial Mode I, II and III SIFs computed by the FEM and the 

BEM for load step 11. 

Figure 2.18: Comparison of typical SIF history computed by the FEM and the BEM for a 

load cycle at the midpoint of the initial crack front. 

-4000

-2000

0

2000

4000

6000

8000

10000

1 3 5 7 9 11 13 15

Load step

K (psi in1/2)

Mode I,FEM Mode I,BEM Mode II,FEM

Mode II,BEM Mode III,FEM Mode III,BEM

-7000

-2000

3000

8000

13000

18000

0 0.04 0.08 0.12 0.16 0.2

Crack front position (in)

K (psi in1/2)

Mode I,FEM Mode I,BEM Mode II,FEM

Mode II,BEM Mode III,FEM Mode III,BEM

Load step 11
NASA/CRó2003-212529 18 

Figure 2.19: Comparison of typical non-proportional load cycle at the midpoint of the 

initial crack front computed by the FEM and the BEM. 

Figures 2.18 and 2.19 present the comparison of SIF variation over one load cycle 

and non-proportional load history. In Figure 2.18, it is seen that KI predictions are lower 

in the FEM analysis than the BEM analysis over most of the steps. As seen in Figure 

2.19, the ratio of KII to KI over 15 load steps shows a similar trend as observed in the 

BEM analysis with a slight difference in values. 

Figure 2.20: Comparison of crack trajectory predicted by the BEM, the FEM and an 

experiment on the tooth surface and tooth cross section at the midpoint of crack. 

-1.4

-1.2

-1

-0.8

-0.6

-0.4

-0.2

0

0.2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

Load step

KII/KI

BEM

FEM

BEM Experiment BEM

FEM

FEM 

Heel 

Toe 

Experiment

Experiment 

BEM 

FEM
NASA/CRó2003-212529 19 

 Crack trajectory comparisons in Figure 2.20 show that, on the tooth surface, crack 

trajectory predicted by FEM is similar to the BEM results. The deviation of crack path at 

the toe end in the BEM predictions was also observed in FEM predictions. On the heel 

end, FEM results exhibit a less steep propagation angle than the experimental trajectory. 

On the tooth cross-section it is observed that the predicted crack trajectory by FEM is less 

steep than BEM and experimental results. FEM results also exhibit a shallower ridge 

formation than the experiment. Although the FEM-predicted crack path still does not 

exactly represent the crack trajectory observed in experiments, it captures the global 

behavior of the crack propagation trajectory. 

3. Three-dimensional Finite Element Contact Analysis 

The meshing of a gear and pinion produces loads on the tooth surface that change 

in magnitude and position in time. In section 2.1 analyses assuming Hertzian contact 

theory were described. This assumption, however, is a simplification for a number of 

reasons: 1) the teeth and rims that transmit the load are flexible, 2) assumed elliptical 

contact area on the tooth surface may not be accurate since the curvatures of the teeth 

surfaces in contact are not constant over the contact zone, 3) the center of contact, where 

the maximum contact stress occurs, may differ from the theoretical contact point, and 4) 

most importantly for the current investigation, the tooth flexibility changes as the crack 

propagates in a tooth. The change in flexibility due to cracks may shift the contact area, 

and change the magnitude of the contact loading. 

In this section a series of three-dimensional finite element analyses that model the 

contact conditions between the gear and pinion explicitly are described. The purpose of 

this study is to assess what influence the details of the load transfer through the contact 

conditions have on the computed stress intensity factors, and through them the predicted 

crack trajectory. 

Three-dimensional finite element contact analyses in spiral bevel gears have been 

undertaken by several researchers to study the stresses induced at the root and on the 

surface of the tooth. Earlier work by Bibel et al. [20] utilized gap elements to simulate the 

contact boundary conditions and to evaluate the stresses on the tooth surface. Further 

study on this subject eliminated the need of using gap elements through the use of contact 

capabilities of commercial finite element programs. Preliminary results showing the 

moving contact patches over the tooth surface are reported in Bibel and Handschuh [21]. 

These analyses also incorporate automatic meshing of the gear and pinion via the 

addition of user subroutines to a commercial finite element program. The same authors 

performed additional analyses in order to compare experimentally measured tooth 

bending stresses of spiral bevel gears to FEM results [22]. Other studies on the contact 

analysis of spiral bevel gears, including the development of a methodology combining a 

surface integral and a finite element solution, are described by Vijakaar [23]. Numerical 

examples carried out using this approach show the change in the contact area as a result 

of the misalignment of the gears. 
NASA/CRó2003-212529 20 

3.1 Approach 

The goal of the present study is to perform contact analyses incorporating cracks 

to determine the effect of changing tooth flexibility on the magnitude and location of the 

contact load and therefore on the stress intensity factors. A commercial finite element 

program, ABAQUS [24], was used to perform the contact analysis in conjunction with the 

software developed by the CFG to calculate fracture-related parameters. A general 

outline of the approach is as follows: 

1. Create initial geometry model of spiral bevel pinion and gear with OSM using 

the geometry information provided by NASA Glenn Research Center. 

2. Specify boundary conditions, material properties, and contact surfaces in 

FRANC3D. 

3. Create a surface mesh composed of triangular elements in FRANC3D. 

4. Create a 3D FE mesh of the model composed of tetrahedra using the J-mesh 

program. 

5. Transform the FRANC3D model into an ABAQUS input file. 

6. Perform contact analysis in ABAQUS and determine the location and magnitude 

of the contact loads. 

7. Calculate SIFs using the results of FEM contact analysis and also using the results 

from analysis with Hertzian contact loads. Compare these values and investigate 

their effects on the crack trajectory predictions. 

 

 Geometries of a one-tooth sector of the spiral bevel pinion and gear were created 

separately in OSM using the information provided by NASA. This information regarding 

the geometric positions of points defining the gear set was determined analytically using 

the method developed by Litvin et al. [3]. The one-tooth gear and pinion models were 

created in a common coordinate system with the same axis of rotation (z-axis). 

Additional teeth of the gear and pinion were created by copying and rotating the original 

tooth by an angle equal to 360o divided by number of gear teeth. The resulting models 

were further rotated to orient them in the proper meshing position. At the meshing 

orientation there was no interference between the gear and pinion teeth. 

The contact analysis was performed using the contact pair approach in ABAQUS.

Contact pairs, named as master and slave surfaces, define surfaces which can potentially 

come into contact. If one of the contact surfaces is stiffer or has a coarser mesh, that 

should be selected as the master surface. In this formulation contact direction is always 

normal to the master surface. In ABAQUS, contact can occur in the form of small-sliding 

or finite-sliding. Analyses performed in this work used the finite-sliding formulation 

where surfaces are allowed to undergo finite separation and sliding. This formulation 

defines the master surface as composed of element faces. Continuity of the surface 

normals was attained by smoothing the transitions between the elements. Contact surface 

tractions were calculated assuming that the surfaces are perfectly hard, which means 

there will be no penetration of surfaces into one another. Contact tractions were defined 

by the local basis system formed by the normal to the master surface. 

The analyses used both FRANC3D and ABAQUS. FRANC3D was used to generate 

the complicated geometry of the gear, to specify the analysis attributes, and to define the 
NASA/CRó2003-212529 21 

crack. Also, it was used as a post-processor to calculate the SIFs. ABAQUS was utilized 

for finite element contact analysis of the model. In order to perform contact analysis on 

models with cracks, features that complement FRANC3D and ABAQUS were needed. In 

this respect FRANC3D was extended to include a capability of specifying the master and 

slave contact surfaces and to include a translator to convert the information regarding the 

geometry, mesh, boundary conditions, material properties and contact surfaces to an 

ABAQUS input file. Also, equivalent domain J-integral implementation was added to 

FRANC3D to obtain more accurate SIF values. 

Contact analyses involve nonlinear geometric analysis. The gear/pinion model is 

meshed with linear tetrahedra. Volume meshing is done using J-mesh as before. The 

mesh is finer on the contacting tooth surfaces of both gear and pinion compared to the 

rest of the model. In performing the analyses, it was assumed that there was no friction 

between the contacting surfaces. 

3.2 Analysis 

 

Finite element contact analysis of a spiral bevel gear and pinion was performed on 

a four-tooth model composed of two pinion teeth and two gear teeth. This model is 

shown in Figure 3.1 and is adapted from [22]. Only four teeth of the gear and pinion were 

explicitly modeled since this is the minimum number of teeth that is needed to simulate 

double tooth contact. We chose to model only a small portion of the gear and pinion to 

decrease the computation time and effort. This model was adequate to simulate the 

contact conditions and also demonstrate how the stress intensity factors change with 

increasing crack length. 

The model was designed such that load applied at a known distance from the axis 

of rotation transmits the required amount of torque. The gear model was extended to 

include a solid, "wedge shaped" part. This portion was a modeling artifact added so that 

the gear could be allowed to rotate freely about its axis of rotation and come into contact 

with the pinion as a result of the applied torque. As seen in Figure 3.1a, the pinion teeth 

were fixed with zero displacements in x,y,z directions on the rim surfaces. Displacement 

components of the gear were fixed in x,y,z directions along the line of axis of rotation. 

Gear rim surfaces were only constrained in the z direction; x and y degrees of freedom of 

the gear were not fixed in order to allow for free rotation in that plane. Dimensions of the 

gear sector were chosen such that the bending effects caused by the application of load is 

minimized. Early studies showed that when smaller portions of the gear were modeled 

contact loads changed with the change in the location of the applied load due to the 

excessive deformation occurring in the gear hub. 

Full design torque for the pinion was 3099 in-lb. In order to obtain the torque that 

is transferred by the gear this value was multiplied by the gear ratio, giving 11580 in-lb. 

A point load was applied at a distance 3.58 in. from the gear axis of rotation to produce 

this amount of torque. The location of the applied load is shown in Figure 3.1a. 
NASA/CRó2003-212529 22 

 (a) 

(b) 

 

 (c) 

 

Figure 3.1: Contact model of spiral bevel gear set. (a) whole model (b) close up view of 

gear and pinion teeth in mesh (c) FE mesh of the model. Note that two teeth of pinion 

and gear are explicitly modeled. Note also the refined mesh on the pinion and gear teeth 

surfaces where contact is expected to occur. 

Pinion tooth 1 

Pinion tooth 2 

Gear tooth 1 Gear tooth 2 

Axis of rotation 

ux, uy, uz =0 

ux, uy, uz =0 

x 

y 

z 

uz =0 

uz =0

ux, uy, uz =0 

Pinion 

Gear 

P 
NASA/CRó2003-212529 23 

The original configuration of the gear and pinion corresponds to a double tooth 

contact load transfer. In order to obtain a single tooth contact configuration the gear set 

was rotated equivalent to a 12o pinion rotation. During the analyses, first a very low 

displacement was applied to the gear in order to ensure contact between gear and pinion 

and to fix the model in space. Then the load was increased to its full design value. 

In order to obtain reasonable contact loads, an adequately refined mesh was very 

important. Since the gear has a doubly curved geometry, a fine discretization was needed 

to model the curved edges. In our studies, we observed that coarse meshes lead to 

inaccurate contact stresses and stress concentrations. To this end, we discretized the 

surfaces where contact is expected to occur fine enough to obtain reasonable results 

without excessive computational cost. 

3.3 Results 

Contact analyses were carried out on four models which have: no crack, an initial 

crack, and two crack sizes, crack #1 and crack #2. Crack #1 is within the fillet region of 

the pinion tooth. Crack #2 is a larger crack extending to the surface of the tooth. These 

four geometries were analyzed for double tooth contact and single tooth contact 

configurations that correspond to 0° and 12° pinion rotation, respectively. 

The number of nodes and elements for the analyzed models ranged between 

12,000-16,000 and 52,000-70,000, respectively. Solution time of a finite element contact 

analysis is about 30 minutes on a 1 GHz Intel Pentium III machine. 

 

(a) (b) 

 (c) 

Figure 3.2: Pinion teeth with a) initial crack b) crack #1 c) crack #2. Note that in all the 

models, crack is in tooth 1. 

Pinion tooth 1 

Pinion tooth 2 

Initial crack Crack #1 

Crack #2 
NASA/CRó2003-212529 24 

3.3.1 Comparison of contact load and area calculated by Hertz theory and threedimensional FEM contact analysis 

Results of contact analyses and how they compare to the Hertz theory are 

summarized in Table 3.1. Total contact load calculated by finite element analysis 

compared well with the total load predicted by Hertz theory. However, distribution of this 

total contact load between two adjacent teeth shows some discrepancies. There is nearly a 

10% difference in the values between contact analysis and Hertz theory predictions on 

each tooth. 

We were unable to attain single tooth contact conditions using this model. At all 

rotations, there is always at least a very little portion carried by the adjacent tooth. This 

might be due to the inaccuracies introduced by the geometric model of the gears. The 

closest case to single tooth loading was obtained when a 12o pinion rotation was 

introduced to the model. This case was assumed to correspond to load step 11 during 

which the highest contact load occurs. 

Contact loads calculated by finite element analyses decrease on tooth 1, where the 

crack is located, and increase on tooth 2 as the crack grows. These results indicate that as 

the crack becomes larger, the flexibility of the pinion tooth is enhanced, therefore the 

contact loads taken up by the two adjacent teeth are redistributed. A larger portion of the 

total contact load has to be carried by the uncracked neighboring tooth. This trend is in 

agreement with the expected behavior. 

Table 3.1: Comparison of contact loads and areas between FEM contact simulations and 

Hertzian stress calculations for two contact configurations with different crack lengths. 

Note that in the first configuration the Hertz contact load given below corresponds to load 

step 3 for tooth 1 and load step 13 for tooth 2, and in the second configuration it 

corresponds to load step 11. 

Model Contact Load (lb) Contact Area (in2) 

 Tooth #1 Tooth #2 Tooth #1 Tooth #2 

Configuration # 1 (Rotation 0o)

No crack 1800 1399 0.0129 0.0126 

Initial crack 1794 1404 0.0129 0.0126 

Crack #1 1792 1406 0.0123 0.0126 

Crack #2 1743 1457 0.0177 0.0126 

Hertz Solution 1584 1602 0.0114 0.0146 

Configuration # 2 (Rotation 12o) 

No crack 2899 318 0.0193 0.0018 

Initial crack 2885 322 0.0193 0.0018 

Crack #1 2881 326 0.0193 0.0018 

Crack #2 2860 348 0.0252 0.0018 

Hertz Solution 3196 - 0.0216 - 
NASA/CRó2003-212529 25 

 

 (a) 

 

 (b) 

Tooth 1

Tooth 2

Tooth 1

Tooth 1

Tooth 2

Tooth 2
NASA/CRó2003-212529 26 

(c) 

Figure 3.3: Contact pressures in on pinion tooth surface for a double tooth contact 

configuration with (a) initial crack (b) crack #1 (c) crack #2. Note that the values on the 

contour bars are in psi. 

For earlier steps of the crack growth, when the crack is still in the fillet of the 

pinion tooth, as in crack #1, there is no significant change in the contact load and area on 

the tooth as seen in Table 3.1. However, when the crack extends into the tooth surface, as 

in crack #2, its effect on the flexibility of the tooth is more pronounced. This results in a 

more significant decrease in the contact load on tooth 1. 

There is no significant difference in the centroid of contact loads calculated by 

finite element contact analysis for different crack lengths. However, locations of contact 

loads determined by these analyses were shifted towards the toe compared to the Hertz 

theory predictions. 

Tooth 2

Tooth 2

Tooth 1
NASA/CRó2003-212529 27 

 

(a)

(b) 

(c) 

Figure 3.4: Contact pressures on pinion tooth surface for an almost single tooth contact 

configuration with (a) initial crack (b) crack #1 (c) crack #2. Note that the values on the 

contour bars are in psi. 

Tooth 1

Tooth 1

Tooth 1
NASA/CRó2003-212529 28 

 Contact areas shown in Figures 3.3 and 3.4 are almost elliptical, showing the 

same trend as the Hertz contact predictions. Also, the distribution of the contact loads is 

close to an ellipsoid shape predicted by Hertzian contact. Maximum contact loads 

calculated by contact analysis were higher than the predictions by the Hertz theory. 

 Results of contact analyses are very sensitive to the mesh size of the contacting 

surfaces. If a coarse mesh is used, computed contact pressure is not continuous on the 

surface since only certain nodes are in contact. In order to overcome this issue, the mesh 

on tooth surfaces was refined several times until a continuous contact pressure 

distribution was obtained. The slightly wavy pattern observed in the contact areas in 

Figures 3.3 and 3.4 are an artifact of the mesh size, shape and orientation. Meshes on the 

gear and pinion tooth surfaces are built such that there is compatibility of the mesh nodes 

between these two contact surfaces. This is important since the contact is calculated 

based on the interaction of mesh nodes between contacting surfaces. If the nodes of one 

of the surfaces fall between the two nodes defining an edge of the other surface, then the 

contact load is taken up by either of the nodes that might lead to different contact 

patterns. The discontinuous contact load pattern observed on crack #2 (Figure 3.3c and 

3.4c) for both configurations might be due to the surface mesh of the tooth. At this crack 

step, the crack was extended to the surface of the pinion tooth. This distorted the 

regularity of the mesh on the pinion surface. Since the compatibility of the meshes on 

both surfaces is very important for the contact analysis, the loss of this compatibility 

might lead to the formation of the contact areas as seen in the above figures. However, 

this issue does not bear a high significance for the calculation of the SIFs. For the 

accuracy of the SIF values, the dominating variable is the total contact load transferred 

rather than the exact distribution. 

 Figures 3.5 and 3.6 show the comparison of location of contact loads by Hertzian 

analysis and FE contact analysis for both single-tooth and double-tooth contact at initial 

crack. Gray squares in the figures define the extreme points of the Hertzian contact 

ellipses. Black points are the mesh nodes where contact occurs in the FE contact analysis. 

The location of contact loads is almost the same at crack #1 and crack #2. 

Figure 3.5: Comparison of location of contact loads by Hertzian analysis and FE contact 

analysis for double-tooth contact at initial crack. Hertzian loads on tooth 1 and tooth 2 

correspond to load step 3 and load step 13, respectively. 
NASA/CRó2003-212529 29 

Figure 3.6: Comparison of location of contact loads by Hertzian analysis and FE contact 

analysis for double-tooth contact at initial crack. Hertzian load on tooth 1 corresponds to 

load step11. 

3.3.2 Comparison of SIFs extracted from three-dimensional FEM contact analysis and 

finite element analysis with Hertz contact loads 

 Change in contact loads due to cracking is expected to have an influence on the 

stress intensity factors and therefore on the crack trajectory. In order to investigate this 

issue, SIFs were calculated using the results of three-dimensional finite element contact 

analysis and finite element analysis using Hertzian contact loads. The latter set of 

analyses were carried out using a two-tooth pinion model where the Hertzian contact 

loads, corresponding load step 3 and load step 11 of the fifteen discretized moving load 

steps described in section 2.1, were applied to the model. 

 SIFs were first calculated by the displacement correlation technique for both 

cases; however, contact analysis results did not compare well with the Hertzian contact 

load results. A 30 % difference in SIFs was observed between the two set of analysis. We 

used linear elements for contact calculations, which is required for this type of analysis, 

and extracted SIFs from these results by the displacement correlation technique. This 

approach is not accurate enough with the linear elements. 

The equivalent domain J-integral method is known to give more accurate results, 

and this method was implemented in FRANC3D. However, SIFs obtained by linear 

elements were still somewhat less accurate than the quadratic elements used for the finite 

element analysis with Hertz contact loads. Improved results with the linear elements were 

obtained as the size of the radius determining the domain of the J-integral calculations 

increased. 

Comparisons between Hertz contact loads and contact loads calculated by threedimensional finite element analysis were given in the previous section. In this section, we 

present the effects of the difference in contact loads calculated by two different methods 

on the SIFs. Figures 3.7 through 3.12 compare mode I, II and III SIFs extracted from the 

results of the analyses with Hertz contact conditions and three-dimensional finite element 

contact analyses. 

SIF values were calculated for the models mentioned in the previous section, 

namely, a pinion with an initial crack, and with cracks corresponding to crack #1 and 
NASA/CRó2003-212529 30 

crack #2. Two loading conditions were investigated, one corresponding to double tooth 

contact and the other corresponding to single tooth contact, assumed to represent load 

step 3 and load step 11, respectively. 

The difference between SIFs calculated by three-dimensional finite element 

contact analyses and analyses using Hertz contact loads increased as the crack length 

increased due to the enhanced flexibility of the tooth. The change in tooth flexibility as 

the crack grew, resulted in lower contact loads and therefore lower SIFs. 

For the gear/pinion model with an initial crack under double-tooth contact, the 

difference was not expected to be that significant since the initial crack was very small 

and did not introduce significant additional flexibility to the tooth. This behavior is 

shown in Figure 3.7. At crack #1, as shown in Figure 3.8, the change in the SIFs became 

more pronounced. An even more amplified effect of the change in tooth flexibility was 

observed for crack #2, as plotted in Figure 3.9. 

Figure 3.7: Comparison of mode I, II and III SIFs calculated by three-dimensional finite 

element contact analysis and finite element analysis with Hertzian contact analysis for a 

pinion with initial crack and at a double tooth contact configuration. 

-3000

-2000

-1000

0

1000

2000

3000

4000

5000

0 0.2 0.4 0.6 0.8 1

normalized crack front

K (psi in1/2)

KI_contact

KI_Hertz

KII_contact

KII_Hertz

KIII_contact

KIII_Hertz
NASA/CRó2003-212529 31 

Figure 3.8: Comparison of mode I, II and III SIFs calculated by three-dimensional finite 

element contact analysis and finite element analysis with Hertzian contact analysis for a 

pinion with crack #1 and at a double tooth contact configuration. 

Figure 3.9: Comparison of mode I, II and III SIFs calculated by three-dimensional finite 

element contact analysis and finite element analysis with Hertz contact loads for a pinion 

with crack #2 and at a double tooth contact configuration. 

-1000

0

1000

2000

3000

4000

5000

0 0.2 0.4 0.6 0.8 1

normalized crack front

K (psi in1/2)

KI_contact

KI_Hertz

KII_contact

KII_Hertz

KIII_contact

KIII_Hertz

-1000

0

1000

2000

3000

4000

0 0.2 0.4 0.6 0.8 1

normalized crack length

K (psi in1/2)

KI_contact

KI_Hertz

KII_contact

KII_Hertz

KIII_contact

KIII_Hertz
NASA/CRó2003-212529 32 

The analyses were repeated for an almost single tooth contact case on all three of 

the crack configurations. The same trends were observed as in the double tooth contact as 

seen in Figures 3.10 through 3.12. However, the difference between the contact analyses 

and analyses with Hertz contact loads were higher for this case. This might be due to the 

fact that the total load applied to the tooth was less than the total load applied in analysis 

with Hertzian contact. As described in section 3.3.1, and shown in Table 3.1, there was at 

least a very little portion carried by the adjacent tooth decreasing the amount of total load 

carried by the cracked tooth. 

Figure 3.10: Comparison of mode I, II and III SIFs calculated by three-dimensional finite 

element contact analysis and finite element analysis with Hertzian contact loads for a 

pinion with initial crack and at a single tooth contact configuration. 

Figures 3.13 through 3.15 show the variation of the ratio of KII to KI along the 

crack front for both double-tooth and single-tooth contact. These graphs are plotted using 

the results obtained by FE contact analysis for all three models with initial crack, crack 

#1 and crack #2. 

-6000

-4000

-2000

0

2000

4000

6000

8000

10000

12000

14000

0 0.2 0.4 0.6 0.8 1

normalized crack front

K (psi in1/2)

KI_contact

KI_Hertz

KII_contact

KII_Hertz

KIII_contact

KIII_Hertz
NASA/CRó2003-212529 33 

Figure 3.11: Comparison of mode I, II and III SIFs calculated by three-dimensional finite 

element contact analysis and finite element analysis with Hertzian contact loads for a 

pinion with crack #1 and at a single tooth contact configuration. 

Figure 3.12: Comparison of mode I, II and III SIFs calculated by three-dimensional finite 

element contact analysis and finite element analysis with Hertzian contact loads for a 

pinion with crack #2 and at a single tooth contact configuration. 

-3000

-1000

1000

3000

5000

7000

9000

0 0.2 0.4 0.6 0.8 1

normalized crack length

K (psi in1/2)

KI_contact

KI_Hertz

KII_contact

KII_Hertz

KIII_contact

KIII_Hertz

-6000

-4000

-2000

0

2000

4000

6000

8000

10000

12000

0 0.2 0.4 0.6 0.8 1

normalized crack front

K (psi in1/2)

KI_contact

KI_Hertz

KII_contact

KII_Hertz

KIII_contact

KIII_Hertz
NASA/CRó2003-212529 34 

Figure 3.13: Ratio of KII to KI along the crack front for double-tooth and single-tooth 

contact at initial crack. 

Figure 3.14: Ratio of KII to KI along the crack front for double-tooth and single-tooth 

contact at crack #1. 

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

normalized crack front

KII/KI

Single-tooth contact

Double-tooth contact

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

normalized crack front

KII/KI

Single-tooth contact

Double-tooth contact
NASA/CRó2003-212529 35 

-1

0

1

2

3

4

5

6

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

normalized crack front

KII/KI

Single-tooth contact

Double-tooth contact

Figure 3.15: Ratio of KII to KI along the crack front for double-tooth and single-tooth 

contact at crack #2. 

By performing these contact analyses, we were able to show that the change in the 

stiffness of the tooth due to a growing crack has significant effects on the SIFs and 

therefore on the crack trajectory and fatigue life. If LEFM analyses of section 2 are 

repeated under these new contact loading conditions, the change in the crack trajectory 

and the fatigue life of the pinion, as a result of taking into account the change in the tooth 

flexibility, can be determined. 

4. Summary 

 This report summarized new results for predicting crack trajectory and fatigue life 

for a spiral bevel pinion using the Finite Element Method (FEM). The predictions 

presented in this report are based on LEFM theories combined with the FEM, 

incorporating plasticity induced fatigue crack closure and moving loads. 

Previously in this project, three-dimensional computational simulation of crack 

growth in spiral bevel gears was done using BEM [1]. As a result of these analyses, a 

method for predicting three-dimensional, non-proportional, fatigue crack growth 

incorporating moving loads was developed. The BEM studies showed the need for 

improving the efficiency and the accuracy of the computations by employing the state-ofthe-art computational capabilities. 

 In the current work, three-dimensional finite element modeling of the spiral bevel 

pinion was done using OSM/FRANC3D. A three-dimensional FEM mesh of the pinion 
NASA/CRó2003-212529 36 

model was created by a 3D mesher. Moving contact loads on the tooth surface were 

calculated by interpolation using the shape functions on the surfaces of the loaded 

elements. The analyses were carried out using a parallel FEM solver, which calculates 

SIFs using equivalent domain J-integral method. Fatigue life predictions were made 

based on a modified Paris model incorporating crack closure. 

In this report, we show that we can simulate fatigue crack growth in a spiral bevel 

gear more accurately and efficiently using the FEM along with a better representation of 

moving loads. Another very significant improvement was the decrease in solving time of 

the problem by employing parallel PC-clusters. This reduces the computation time to a 

few minutes while substantially increasing the resolution of the solution. 

To obtain a more detailed understanding of the contact between a cracked pinion 

tooth in mesh with an uncracked gear tooth, three-dimensional contact analyses were 

performed on a spiral bevel gear set incorporating a crack. The goal in carrying out these 

analyses was to capture the redistribution of contact loads due to crack growth. Results of 

these analyses showed the expected trend of decreasing tooth loads carried by the cracked 

tooth with increasing crack length. We also showed that this decrease in contact loads 

had an impact on the SIF values and therefore would also affect the crack trajectory and 

fatigue life predictions. 

 Predicting crack trajectories more efficiently and accurately is very important for 

gear designers. Development of practical and accurate numerical tools for evaluating the 

gear performance and life enables for new and better gear designs. The work presented 

here provides such a tool capable of predicting crack trajectory and fatigue life in a spiral 

bevel gear. 
NASA/CRó2003-212529 37 

Appendix: Abbreviated sample ABAQUS input file for contact analysis 

*HEADING 

FRANC3D Model 

** FRANC3D Version 2.0 Generated Input File 

** NUMBER OF NODES: 13029 

** NUMBER OF ELEMENTS: 58712 

** 

** NODE DEFINITION 

*NODE 

1, 3.815870, 0.262049, 1.500000 

2, 3.585600, 0.038081, 1.500000 

3, 3.824640, 0.040620, 1.500000 

. 

. 

. 

13028, 5.109482, 0.491277, 0.734872 

13029, 5.110899, 0.478949, 0.730531 

** END NODE DEFINITION 

** 

** ELEMENT DEFINITION 

*ELEMENT, TYPE=C3D4, ELSET=PID1 

1,1904,1855,2100,1903 

*ELEMENT, TYPE=C3D4, ELSET=PID1 

2,1904,2100,1953,1952 

*ELEMENT, TYPE=C3D4, ELSET=PID1 

3,125,107,106,846 

. 

. 

. 

*ELEMENT, TYPE=C3D4, ELSET=PID2 

58711,7717,12943,7703,13029 

*ELEMENT, TYPE=C3D4, ELSET=PID2 

58712,6078,7717,7703,13029 

** END ELEMNT DEFINITION 

** 

*SOLID SECTION, ELSET=PID1, MATERIAL=MAT1 

*MATERIAL, NAME=MAT1 

*ELASTIC, TYPE=ISO 

 3e+07, 0.3 

*DENSITY 

 7800 

*SOLID SECTION, ELSET=PID2, MATERIAL=MAT2 

*MATERIAL, NAME=MAT2 

*ELASTIC, TYPE=ISO 

 3e+07, 0.3 

*DENSITY 

 7800 

** 

**Definition of contact for tooth 1 

** 

*ELSET, ELSET=MASTER_1 

26454, 

26457, 

26458, 

. 
NASA/CRó2003-212529 38 

. 

. 

57531, 

58366, 

58367, 

*ELSET, ELSET=SLAVE_1 

1, 

2, 

5, 

. 

. 

. 

22208, 

24437, 

24438, 

*SURFACE DEFINITION, NAME=MASTER_1 

26454,S4 

26457,S4 

26458,S4 

. 

. 

. 

57531,S1 

58366,S1 

58367,S1 

*SURFACE DEFINITION, NAME=SLAVE_1 

1,S2 

2,S4 

5,S4 

. 

. 

. 

22208,S1 

24437,S1 

24438,S1 

*SURFACE INTERACTION, NAME=NOFRIC 

*FRICTION 

 0 

*CONTACT PAIR, INTERACTION=NOFRIC 

MASTER_1, SLAVE_1 

** 

** Definition of contact for tooth 2 

** 

*ELSET, ELSET=MASTER_2 

26455, 

26456, 

26459, 

. 

. 

. 

57320, 

57753, 

57754, 

*ELSET, ELSET=SLAVE_2 

3, 

4, 

9, 
NASA/CRó2003-212529 39 

. 

. 

. 

24532, 

24540, 

24548, 

*SURFACE DEFINITION, NAME=MASTER_2 

26455,S4 

26456,S4 

26459,S4 

. 

. 

. 

57320,S1 

57753,S1 

57754,S1 

*SURFACE DEFINITION, NAME=SLAVE_2 

3,S2 

4,S4 

9,S4 

. 

. 

. 

24532,S1 

24540,S1 

24548,S1 

*SURFACE INTERACTION, NAME=NOFRIC1 

*FRICTION 

 0 

*CONTACT PAIR, INTERACTION=NOFRIC1 

MASTER_2, SLAVE_2 

** 

** Application of displacement to ensure contact 

** between gear and pinion 

** 

*STEP, NLGEOM, INC=100 

**FRANC3D load case 1 

*STATIC 

1.0, 1.0 

*BOUNDARY 

2209,3,3, 0 

2211,3,3, 0 

2213,3,3, 0 

2215,3,3, 0 

2217,3,3, 0 

2219,3,3, 0 

2221,3,3, 0 

2223,3,3, 0 

2232,2,2, -0.01 

2247,3,3, 0 

. 

. 

. 

33,1,1,0 

33,2,2,0 

33,3,3,0 

2322,1,1,0 
NASA/CRó2003-212529 40 

2322,2,2,0 

2322,3,3,0 

2415,1,1,0 

2415,2,2,0 

2415,3,3,0 

*END STEP 

** 

** Application of full design torque 

** 

*STEP, NLGEOM, INC=100 

**FRANC3D load case 1 

*STATIC 

1.0, 1.0 

*BOUNDARY, OP=NEW 

2209,3,3, 0 

2211,3,3, 0 

2213,3,3, 0 

2215,3,3, 0 

2217,3,3, 0 

2219,3,3, 0 

2221,3,3, 0 

2223,3,3, 0 

2247,3,3, 0 

. 

. 

. 

33,1,1,0 

33,2,2,0 

33,3,3,0 

2322,1,1,0 

2322,2,2,0 

2322,3,3,0 

2415,1,1,0 

2415,2,2,0 

2415,3,3,0 

*CLOAD 

2799,2,-1000 

*END STEP 

*STEP, NLGEOM, INC=100 

**FRANC3D load case 1 

*STATIC 

1.0, 1.0 

*CLOAD 

2799,2,-3230 

*FILE FORMAT, ASCII 

*NODE FILE, FREQUENCY=100 

U 

*EL FILE, FREQUENCY=100, POSITION=NODES 

S,E,PE 

*RESTART,WRITE, FREQUENCY=100 

*NODE PRINT, FREQUENCY=100 

U 

*EL PRINT, FREQUENCY=100, POSITION=NODES 

*CONTACT PRINT,FREQUENCY=100 

CFN,CAREA,XN 

*END STEP 

**END OF INPUT 
NASA/CRó2003-212529 41 

References 

[1] Spievak, L. E., Wawrzynek, P. A., Ingraffea A. R., and Lewicki, D. G. 

ìSimulating Fatigue Crack Growth in Spiral Bevel Gears.î Engineering Fracture 

Mechanics 68, no. (2001):53-76. 

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

1985. 

[3] Litvin, F. L., and Zhang, Y. ìLocal Synthesis and Tooth Contact Analysis of 

Face-Milled Spiral Bevel Gearsî NASA Contractor Report 4342/AVSCOM 

Technical Report 90-C-028, 1991. 

[4] Carter, B., Chen, C-S , Chew, P., Chrisochoides, N., Gao, G.R., Heber, G., 

Ingraffea, A.R., Krause R., Myers, C., Nave, D., Pingali, K., Stodghill P., 

Vavasis, S., Wawrzynek P.A. ìParallel FEM Simulation of Crack Propagationó

Challenges, Status, and Perspectives.î Lecture Notes in Computer Science 1800, 

(2000):443-449. 

[5] Chan S.K., Tuba I.S. and Wilson W.K. ìOn the finite element method in linear 

fracture mechanics.î Engineering Fracture Mechanics 2, no.1, (1970): 1-17. 

[6] Barsoum R.S. ìOn the use of isoparametric elements in linear fracture 

mechanics.î International Journal for Numerical Methods in Engineering 10, 

(1976): 25-37. 

[7] Shih C.F., deLorenzi H.G. and German M.D. ìCrack extension modeling with 

singular quadratic isoparametric elements.î International Journal of Fracture 12, 

(1976): 647-651. 

[8] Ingraffea A.R. and Manu C. ìStress intensity factor computation in three 

dimensions with quarter-point elements.î International Journal for Numerical 

Methods in Engineering 15, (1980): 1427-1445. 

[9] deLorenzi H.G. ì3D elastic-plastic fracture mechanics with ADINA.î 

International Journal of Computers and Structures 13, (1981): 613-621. 

[10] Li F.Z., Shih C.F. and Needleman A. ìA comparison of methods for calculating 

energy release rates.î Engineering Fracture Mechanics 21, 1985: 405-421. 

[11] Shih C.F., Moran B. and Nakamura T. ìEnergy release rate along a threedimensional crack front in a thermally stressed body.î International Journal of 

Fracture 30, (1986): 79-102. 

[12] Moran B. and Shih C.F. ìA general treatment crack tip contour integrals.î 

International Journal of Fracture 35, (1987): 295-310. 

[13] Moran B. and Shih C.F. ìCrack tip and associated domain integrals from 

momentum and energy balance.î Engineering Fracture Mechanics 27, no.6, 

(1987): 615-642. 

[14] Dodds Jr., R.H. and Vargas, P.M. ìNumerical evaluation of domain and contour 

integrals for nonlinear fracture mechanics: formulation and implementation 

aspects.î Report, University of Illinois at Urbana-Champaign, Department of 

Civil Engineering, 1988. 

[15] Rice J.R. ìA path independent integral and the approximate analysis of strain 

concentration by notches and cracks.î Journal of Applied Mechanics, 35 (1968): 

379-386. 
NASA/CRó2003-212529 42 

[16] Bui H.D. ìAssociated path independent J-integrals for separating mixed modes.î 

Journal of the Mechanics and Physics of Solids 31, no.6, (1983): 439-448. 

[17] Nikishov G. P. and Atluri S.N. ìCalculation of fracture mechanics parameters for 

an arbitrary three-dimensional crack, by the ëEquivalent Domain Integralí 

method.î International Journal for Numerical Methods in Engineering 24, no.9 

(1987): 1801-1821. 

[18] Chen K.L., and Atluri S.N. ìComparison of different methods of evaluation of 

weight functions for 2D mixed-mode fracture analysis.î Engineering Fracture 

Mechanics 34, no.4 (1989): 935-956. 

[19] Murakami Y. (Editor-in-Chief) Stress Intensity Factors Handbook. Pergamon 

Press, (1987). 

[20] Bibel, G.D., Kumar, A., Reddy, S., and Handschuh, R. ìContact stress analysis of 

spiral bevel gears using finite element analysis.î Journal of Mechanical Design

117, no.2A (1995): 235-240. 

[21] Bibel, G.D., and Handschuh, R. ìMeshing of a Spiral Bevel Gearset with 3D 

Finite Element Analysis.î NASA Technical Report ARL-TR-1224 (1996). 

[22] Handschuh, R.F., and Bibel, G.D. ìComparison of experimental and analytical 

tooth bending stress of aerospace spiral bevel gears.î NASA Technical Report 

ARL-TR-1891 (1999). 

[23] Vijakaar, S. ìA combined surface integral and finite element solution for a threedimensional contact problem.î International Journal for Numerical Methods in 

Engineering 31, no. (1991): 525-545. 

[24] ABAQUS/Standard, Version 6.1.1, Hibbitt, Karlsson & Sorensen Inc. 
This publication is available from the NASA Center for AeroSpace Information, 301–621–0390.

REPORT DOCUMENTATION PAGE

2. REPORT DATE

19. SECURITY CLASSIFICATION

 OF ABSTRACT

18. SECURITY CLASSIFICATION

 OF THIS PAGE

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources,

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this

collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson

Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188), Washington, DC 20503.

NSN 7540-01-280-5500 Standard Form 298 (Rev. 2-89)

Prescribed by ANSI Std. Z39-18

298-102

Form Approved

OMB No. 0704-0188

12b. DISTRIBUTION CODE

8. PERFORMING ORGANIZATION

 REPORT NUMBER

5. FUNDING NUMBERS

3. REPORT TYPE AND DATES COVERED

4. TITLE AND SUBTITLE

6. AUTHOR(S)

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

11. SUPPLEMENTARY NOTES

12a. DISTRIBUTION/AVAILABILITY STATEMENT

13. ABSTRACT (Maximum 200 words)

14. SUBJECT TERMS

17. SECURITY CLASSIFICATION

 OF REPORT

16. PRICE CODE

15. NUMBER OF PAGES

20. LIMITATION OF ABSTRACT

Unclassified Unclassified

Final Contractor Report

Unclassified

1. AGENCY USE ONLY (Leave blank)

10. SPONSORING/MONITORING

 AGENCY REPORT NUMBER

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546–0001

and

U.S. Army Research Laboratory

Adelphi, Maryland 20783–1145

Available electronically at http://gltrs.grc.nasa.gov

August 2003

NASA CR—2003-212529

ARL–CR–0531

E–14091

48

Simulating Fatigue Crack Growth in Spiral Bevel Pinion

Ani Ural, Paul A. Wawrzynek, and Anthony R. Ingraffea

Spiral-bevel gears; Cracks; Finite element analysis; Boundary element analysis; Fracture

mechanics; Stress intensity factors; Contact analysis

Unclassified -Unlimited

Subject Category: 37 Distribution: Nonstandard

WBS–22–708–90–01

NAS3–1993

L162211A47A

Cornell University

Cornell Fracture Group

Rhodes Hall

Ithaca, New York 14850

Project Manager, David Lewicki, Structures and Acoustics Division, NASA Glenn Research Center, organization code

5950, 216–433–3970.

New results for predicting crack trajectory and fatigue life for a spiral bevel pinion using the Finite Element Method

(FEM) are reported. The predictions presented are based on linear elastic fracture mechanics combined with the FEM,

incorporating plasticity induced fatigue crack closure and moving gear tooth loads. The analyses were carried out using a

parallel FEM solver, which calculates stress intensity factors using equivalent domain J-integral method. Fatigue life

predictions were made based on a modified Paris model incorporating crack closure. To obtain a more detailed understanding of the contact between a cracked pinion tooth in mesh with an uncracked gear tooth, three-dimensional contact

analyses were performed on a spiral bevel gear set incorporating a crack. The goal in carrying out these analyses was to

capture the redistribution of contact loads due to crack growth. Results of these analyses showed the expected trend of

decreasing tooth loads carried by the cracked tooth with increasing crack length. It was also showed that this decrease in

contact loads had an impact on the stress intensity factor values and therefore would also affect the crack trajectory and

fatigue life predictions.