UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 1 CRITICAL FLAW SIZE IN SILICON NITRIDE BALL BEARINGS By GEORGE ARTHUR LEVESQUE A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTO R OF PHILOSOPHY UNIVERSITY OF FLORIDA 2009 PAGE 2 2 2009 George Arthur Levesque PAGE 3 3 To My Father PAGE 4 4 ACKNOWLEDGMENTS This project was funded by the United States Air Force under Versatile Affordable Advanced Turbine Engines (VAATE) Contract# F3361503D 2353 003 and managed by Robert Wolfe, Manager, Materials Technology, Timken Company, Canton, O hio Robert Wolfe is thanked for extensive discussions on all aspects of critical flaw size determination. William Ogden and Herb Chin of Pratt & Whitney Co. Ne lson Forster and Vaughn Svendsen of the Air Force Research Labs and Michael J. OBrien of the Aerospace Corporation are thanked for their insightful contributions. This work has been done while supported by the Timken Company, Canton, OH, and Air Force Res earch Labs, Dayton, OH. At the University of Florida there were two collaborations worth noting. Dr. John J. Mecholsky and Karthik Gopalakrishnan have contributed an experimental endeavor on the Brazilian disc test. Also, Dr. Nam Ho Kim and Sriram Pattabhiraman have contributed in the endeavor of an uncertainty analysis of rolling element survival. In addition, my labmates and my advisor Dr. Nagaraj Arakere are thanked for their support throughout my graduate career PAGE 5 5 TABLE OF CONTENTS P age ACKNOWLEDGMENTS ...................................................................................................... 4 LIST OF TABLES ................................................................................................................ 8 LIST OF FIGURES .............................................................................................................. 9 L IST OF ABBREVIATIONS .............................................................................................. 14 NOMENCLATURE ............................................................................................................ 15 ABSTRACT ........................................................................................................................ 19 CHAPTER 1 INTRODUCTIO N ........................................................................................................ 21 1.1 General Background ............................................................................................. 21 1.1.1 Bearing Design Issues ................................................................................ 21 1.1.2 Issue s With Silicon Nitride .......................................................................... 22 1.2 Objective And Scope Of Research Work ............................................................ 23 1.3 Literature Review .................................................................................................. 24 1.3.1 Rolling Contact Fatigue Loading In Ball Bearings ..................................... 24 1.3.2 Limitations Of LundbergPalmgren Theory ................................................ 26 1.3.3 Hertzian Fracture Of Brittle Materials ......................................................... 27 1.3.4 Computational Simulation Of Surface Flaws ............................................. 28 1.3.5 MixedMode Fracture Testing .................................................................... 30 1.4 State Of The Art From Literature Survey ............................................................. 32 1.5 Outline Of Thesis .................................................................................................. 33 1.6 Figures .................................................................................................................. 37 2 DETERMINATION OF SURFACE FLAW GEOMETRY............................................ 45 2.1 Introduction ........................................................................................................... 45 2.2 Analysis ................................................................................................................. 46 2.2.1 Stress State Analysis .................................................................................. 47 2.2.2 Incremental Crack Growth Analysis ........................................................... 53 2.3 Results On Crack Shape Determination .............................................................. 55 2.4 Conclusions .......................................................................................................... 57 2.5 Figures .................................................................................................................. 59 PAGE 6 6 3 GENERAL MODELING ISSUES OF SURFACE FLAWS SUBJECT TO ROLLING CONTACT FATIGUE ................................................................................ 70 3.1 Introduction ........................................................................................................... 70 3.2 Com puting Stress Intensity Factors ..................................................................... 71 3.3 Parametric Studies ............................................................................................... 73 3.4 Physical Effects .................................................................................................... 74 3.5 Conclusions .......................................................................................................... 75 3.6 Figures .................................................................................................................. 76 4 EXPERIMENTAL ANALYSIS ..................................................................................... 84 4.1 A Novel Fracture Toughness Test ....................................................................... 84 4.2 V Ring Test ........................................................................................................... 85 4.2.1 Experimental Procedure to Determine Growth Regime ............................ 85 4.2.2 Specimen Analysis ..................................................................................... 86 4.2.3 Stress Intensity Factor Determination Of Individual Test Specimens ....... 86 4.2.4 Calculating A Keff ......................................................................................... 87 4.3 Brazilian Disc Test ................................................................................................ 87 4.3.1 Experimental Procedure ............................................................................. 88 4.3.2 Analytical Expressions for MixedMode Stress Intensity Factors ............. 89 4.3.3 Stress Intensity Factor Calculation Via Finite Element Modeling ............. 89 4.3.4 MixedMode Fracture Criteria ..................................................................... 92 4.3.5 Analysis Procedure ..................................................................................... 96 4.3.6 Results And Discussion .............................................................................. 97 4.4 Conclusions .......................................................................................................... 99 4.5 Figures ................................................................................................................ 100 5 CFS EVALUATION: COMPUTATIONAL INVESTIGATION ................................... 117 5.1 Introduction ......................................................................................................... 117 5.2 Analytical Procedure........................................................................................... 120 5.2.1 Cra ck Geometries Induced By Oblique Interactions Of Silicon Nitride Balls ................................................................................................................. 120 5.2.2 Finite Element Modeling Of A ThreeDimensional Surface Flaw ............ 122 5.3 Mixed Mode Stress Intensity Factors Due to Rolling Contact Fatigue ............. 124 5.3.1 Effect Of The Variation Of Load Ellipticity On Stress Intensity Factors .. 126 5.3.2 Comparison Of C Cracks And Semi Elliptical Cracks ............................. 127 5.3.3 The Effect Of Traction ............................................................................... 128 5.3.4 C ritical Flaw Size Evaluation Results ....................................................... 129 5.4 Conclusions ........................................................................................................ 130 5.5 Figures ................................................................................................................ 13 2 PAGE 7 7 6 EMPIRICAL STRESS INTENSITY FACTOR EQUATIONS FOR SURFACE CRACKS UNDER RCF ............................................................................................ 142 6.1 Introduction ......................................................................................................... 142 6.2 Analysis ............................................................................................................... 143 6.3 Crack Geometry .................................................................................................. 144 6.4 Orientation Of Load ............................................................................................ 145 6.5 Load .................................................................................................................... 146 6.6 Results ................................................................................................................ 148 6.7 Conclusions ........................................................................................................ 149 6.8 Figures ................................................................................................................ 157 7 UNCERTAINTY ANALYSIS FOR FATIGUE FAILURE PROBABILITY OF SILICON NITRIDE ROLLING ELEMENT ................................................................ 161 7.1 Introduction ......................................................................................................... 161 7.2 Modeling Rolling Contact Fatigue Orientation Effects ...................................... 165 7.3 Uncertainty Analysis ........................................................................................... 167 7.3.1 Uncertainty Quantification ........................................................................ 167 7.3.2 Surrogate Modeling Uncertainty Propagation ....................................... 169 7.3.3 Uncertainty Analysis Monte Carlo Simulation Using Surrogate ........... 173 7.3.4 Parameter Study ....................................................................................... 176 7.4 Conclusions ........................................................................................................ 178 7.5 Figures ................................................................................................................ 179 8 DISCUSSION AND CONCLUSIONS ....................................................................... 186 8.1 Discussion ........................................................................................................... 186 8.2 Conclusions ........................................................................................................ 188 8.3 Future Work ........................................................................................................ 188 APPENDIX A THE NORMAL INDENTATION OF BRITTLE SPHERES ....................................... 190 B BUILDING A CCRACK MODEL .............................................................................. 195 LIST OF REFERENCES ................................................................................................. 216 BIOGRAPHICAL SKETCH .............................................................................................. 223 PAGE 8 8 LIST OF TABLES Table P age 1 1 Comparison of material properties between fully densified silicon nitride and a common ball material predecessor M50 steel. .................................................. 37 1 2 Mea surements for a sense of scale of the problem. ............................................. 41 1 3 RCF tests on ceramics and their failure modes in literature. ................................ 41 2 1 Estimated c crack geometries produced by different values of from stress state analysis. ......................................................................................................... 62 4 1 Vring test conditions for 28.575 mm (1 in.) balls. ............................................... 105 4 2 The sample number, turning angle, crack dimensions and fracture load of each specimen examined. ................................................................................... 106 4 3 The mean Kc values in MPa std). ..................................................................................................................... 106 6 1 The xd distances (nondimensinoalized w.r.t. the crack semi width) where KI was maximal. ........................................................................................................ 150 6 2 Coefficients for all parametric cases for KI, KII, and KIII ...................... 151 6 3 Coefficients for all parametric cases for KI, KII, and KIII .................... 153 6 4 Coefficients for all parametric cases for KI, KII, and KIII .................... 155 7 1 Input parameters and their distributions for uncertainty analysis ....................... 181 7 2 Cross validation errors (PRESSRMS) for different surrogate models ................ 182 7 3 Probability of failure values for different sample sizes ........................................ 182 7 4 Effects of improvements pertaining to crack size on probability of failure ......... 183 7 5 Effects of improvements pertaining to fracture toughness on probability of failure .................................................................................................................... 183 7 6 Effects of combined improvements pertaining to crack size and fracture toughness on probability of failure ....................................................................... 183 PAGE 9 9 LIST OF FIGURES Figur e P age 1 1 grains are apparent. ............................................................................................... 37 1 2 A) Example of a cone crack in silicon nitride B) Experimental image of oblique spher e interaction with contact patch and resultant c crack shown.. ...... 38 1 3 Before and after of Timken test #7 of the V ring single ball test. .......................... 39 1 4 Image of steel raceway spall from RCF. ............................................................... 40 1 5 The elliptical contact patches between the ball and raceways. ............................ 42 1 6 F EA simulation of a steel raceway section being subjected to contact with a sphere. .................................................................................................................... 43 1 7 Fracture stress as a function of dominant defect dimension. ............................... 44 1 8 Image of a spalled c crack clearly showing its face as it extends into the material. .................................................................................................................. 44 2 1 A) Example of a cone crack in silicon nitride B) Experimental image of obl ique sphere interaction with contact patch and resultant c crack shown. ....... 59 2 2 Diagram of coordinate system for the equations of stress for oblique colliding spheres. .................................................................................................................. 60 2 3 Surfaces representing constant values of the max principal stress field. ............ 61 2 4 A) C one crack as generated from the Hertzian stresses without traction. B) A c crack generated for v =0.3 and =0.1. ................................................................ 62 2 5 Initial configuration of crack and load at the beginning of the crack growth analysis in the FE model. ....................................................................................... 63 2 6 Incremental images of crack growth for =0.1. ..................................................... 64 2 7 Incremental images of crack growth for =0.5. ..................................................... 65 2 8 Cracks produced for = 0.1, 0.2, 0.3, 0.4, and 0.5 and v =0.28. ........................... 66 2 9 Cross sections of multiple natural c cracks on 12.7 mm diameter balls that (after manufacturing) had a typical ...................................... 67 2 10 Cross sections of artificially produced c cracks. ................................................... 68 PAGE 10 10 2 11 For =0.1, s uperposed images of half cracks p roduced by the stress stateand the incremental growth. .......................................................................... 69 3 1 A partial cone crack geometry generated from some experimental cross sections of the crack created with a spline technique. .......................................... 76 3 2 General J integral equation reference coordinate system. ................................... 76 3 3 Crack tip coordinate systems for CTOD. ............................................................... 77 3 4 A c crack specimen with subsurface crack face highlighted above. In tension below, indicating opening of the crack faces. ........................................................ 78 3 5 Elements from the cra ck tip in fig 3 3 reoriented to show that different parts of the crack tip undergo all three of the different types of deformation. ................... 79 3 6 Visual illustration of the Matlab code which is used t o create, submit, post process, and organize the desired FEA results. ................................................... 80 3 7 An illustration of a ball in a ball bearing raceway system to indicate the location and orientation of the contact p atches. .................................................... 80 3 8 Relative orientation of an elliptical load near a surface crack. ............................. 81 3 9 KI for =60 near r =1b load for varying aspect ratio. ............................................ 82 3 10 A) worst case orientation of crack and load but if the traction direction is reversed as in B) the crack tends to close. ........................................................... 82 3 11 KI as a function of the coefficient of friction on the surface for po=540 ksi (3.7 GPa), Crack angle = 30 b =250 mm, a =75 mm, a/b =0.3. ................................... 83 4 1 On the left, the designed fixture and a tested ball are displayed. On the right, the platens and ball are placed into the fixture in the orientation for testing. ..... 100 4 2 Closeup of indentation flaw, with measurement taken in the direction of growth due to applied load. Arrows indicate direction of tensile field. ................ 101 4 3 The fracture toughness of tested specimens versus the final crack length they became, indicating a lack of interdependence. ........................................... 101 4 4 The fracture toughness of tested specimens versus the Vickers indentation load indicating a lack of interdependence. .......................................................... 102 4 5 FEA of the ball in the platens undergoing compression, clearly indicating the tensile field along the equator of the ball. ............................................................ 102 4 6 Pendulum device for impacting balls and inducing cracks ................................. 103 PAGE 11 11 4 7 Images of V ring single ball test rig. Lower image shows contact path and cracks oriented within them. ................................................................................ 104 4 8 Diagram illustrating measurement parameters to classify a c crack as observed on the surface of a cracked specimen. ............................................... 105 4 9 Schematic showing the flaw orientation on a Brazilian disk specimen as well as the coordinates and components used in the analytical stress equations. ... 107 4 10 Nondimensionalized KI and KII from empirical equat ions for different angles of loading orientation for b =.02m a =.0025m. ...................................................... 108 4 11 BD test model showing a) global uncracked model with square partition indicating submodel boundaries b) cracked me shed submodel. ........................ 109 4 12 The displacement field of Sample 1 for =20, to illustrate the presence of KII at the face and KIII at the depth. ........................................................................... 109 4 13 KI for all analyzed specimens as a function of position along the crack front. ... 110 4 14 KII for all analyzed specimens as a function of position along the crack front. .. 111 4 15 KIII for all analyzed specimens as a function of position along the crack front. 112 4 16 Keff using Tanakas relationship of three SIFs as a function of crack front position. ................................................................................................................ 113 4 17 Keff using the strain energy density release rate using all three SIFs as a function of crack front position where a) includes KIII and b) neglects it. ........... 114 4 18 Plot of 2 2 Mi IIIIIIKKBKK where B=1. ................................................... 115 4 19 Fractured specimens at various values of ....................................................... 115 4 20 Fracture surface of BD specimen showing the crack front growing as it links up with newly formed cracks from indents along the surface. ............................ 116 5 1 Images of artificially induced c cracks a) ball surface view, and b) a subsurface crosssection extending into the material ........................................ 132 5 2 The elliptic contact patches on th e ball surface and the band which they remain in. .............................................................................................................. 133 5 3 A) The coordinate system used in the equations to plot the C crack and B) the rendered 3D crack with notation for the SIF plots to be discussed. ............. 1 34 5 4 Illustration of A ) displacement correlation variables around an opening two dimensional crack and B) the coordinate systems near the crack tip. ............... 135 PAGE 12 12 5 5 Elliptical Hertzian contact load coordinate system and variables for elliptical load simulation. .................................................................................................... 135 5 6 The cross section of the general worst case orientation with tra ction for crack tip deformation ..................................................................................................... 136 5 7 KI stress intensity factor for a circular load ( ac crack with a variation in the max pressure magnitude. ........................................................... 136 5 8 A plot of the three SIFs as a function of moving load of Po=2.8 GPa. ............... 137 5 9 The effect of the ellipticity of the load, on the SIFs on the semi minor axis of the ellipse, for Po=2.8 GPa (412 Kpsi). ................................................................ 138 5 10 Superimposed SIFs for a c crack and penny crack of similar dimensions under elliptical load. ............................................................................................. 139 5 11 The max KI found on the crack front of a semi elliptical flaw as a function of varying friction coefficients. .................................................................................. 140 5 12 The nondimensionalized crack s emi width plotted against a) the max KIeq and b) the max KI ....................................................................................................... 141 6 1 Model configuration displaying orientation of load and defining variables and dx .................................................................................................................. 157 6 2 KI for a single load geometry with and without contact defined. ......................... 158 6 3 KI for r=1 b along the crack front for all variations of crack angle and depth. ..... 159 6 4 KI for r=2 b along the crack front for all variations of crack angle and depth. ..... 159 6 5 KI for r=3 b along the crack front for all variations of crack angle and depth. ..... 160 7 1 Before and after of a partial cone crack placed under rolling contact fatigue. ... 179 7 2 Elliptical Hertzian contact load coordinate system and variables for elliptical load orientation simulation. .................................................................................. 180 7 3 The elliptic contact patches on the ball surface and the band which they remain in. .............................................................................................................. 180 7 4 Surrogate models for the contribution of lateral position to the stress intensity factor ..................................................................................................................... 181 7 5 Surrogate models for the contribution of crack orientation to the stress intensity factor ...................................................................................................... 182 PAGE 13 13 7 6 Cumulative distribution function of the limit state function g = Keq Kth ............. 183 7 7 Variation of probability of failure with percentage reduction of crack size. ........ 184 7 8 Variation of probability of survival with crack size. .............................................. 184 7 9 Variation of probability of survival with crack size for a Kth of 4.50.54 MPa m. ................................................................................................................. 185 PAGE 14 14 LIST OF ABBREVIATION S Abbreviation Definition 3D Three dimensional BD Brazilian disc BEA Boundary element analysis CFS Critical flaw size EHD Elastohydrodynamic FE Finite el ement FEA Finite element analysis FPI Fluorescent penetrant inspection MCS Monte Carlo simulation NDE Non destructive evaluation RCF Rolling contact fatigue SIF Stress intensity Factor \ VAATE Versatile Affordable Advanced Turbine Engine s w.r.t. With respect to PAGE 15 15 NOMENCLATURE a Semi elliptical crack depth a` Axis dimension of ellipse in x direction a* Revolved ellipse axis dimension in the x direction for a c crack ac Radius of contact patch size ia Coe fficients of a polynomial for KI b Semi elliptical crack semiwidth b` Axis dimension of ellipse in z direction b* Axis dimension of revolved ellipse in the z direction for a c crack ib Coefficients of a polynomial for KII dc Distance between crack and load edge ic Coefficients of a polynomial for KIII ds Length element along the contour e 21 e E Youngs Modulus 12, EE Youngs Moduli for bodies 1 and 2, respectively E Effective Youngs Modulus, 1 12 1211vv EE f Coefficient of Friction g Limit state function G Strain energy release rate, 22 21 2(1)III IIIv GKKK EE h x coordinate of crack at surface on x z plane. PAGE 16 16 Ij Index function JR Rices contour integral J 3 432 3 xyzz ruu k z coordinate of crack at surface on x z plane. Keq 0.5 22 21 2(1)eq III IIIvKEKKK EE eqK */`eqeqOKKPa KI Opening stress intensity factor *iK i oK a p Q where ,, iIIIIII KII Sliding stress intensity factor KIII Tearing stress intensity factor Kth Mixed mode t hreshold SIF for crack growth m Ball mass (when single ball is discussed) m Effectiv e ball mass, (for contacting balls m = 1 1211 mm N A harmonic potential function (Love, 1927) N Number of Monte Carlo simulation samples Nf Number of simulation samples that fail NC Number of simulation samples t hat lie in the contact patch O Coordinate origin po Max value of pressure in elliptical pressure dome PAGE 17 17 op Max Pressure to induce cracking P Total load PF Probability of failure Q Shape factor for elliptical crack 12, RR Radius of sphere 1 and 2, respectively r Radial cylindrical coordinate, 22xy for circular contact rd Distance from crack tip where u, v, and w are measured R Effective ball radius for a ball on ball interaction, 1 1211 R RR T A harmonic potential function (Love, 1927) with partial derivatives Ti Components of traction vector Friction coefficient u 2 222221 114 2 rzrzz u, v, w Crack face opening, sliding, and shearing displacements ui Displacement vector components w Strain energy density Poissons ratio V Velocity of moving ball zV Velocity component of moving ball normal to the contact surface xd x displacement of pressure distribution from coordinate origin xD Lateral position of crack to contact patch center yd y displacement of pressure distribution from coordinate origin PAGE 18 18 Aspect ratio of load ( b/a) Arbitrary counterclockwise path around the crack tip Material density Position along crack front pl Angle of principal plane Angle of crack inclination towards the vertical max Half the total amount which the crack subtends a circle. c Critical max periphery stress to induce cracking Elastic potential P oissons ratio 12, Poissons ratios for bodies 1 and 2 respectively 123,, First, second, and third principal stresses where 123 c Contact periphery stress to indu ce cracking max Contact periphery stress max f Contact periphery stress with friction N ij Stress components from normal load T ij Stress components from traction l oad PF Standard deviation of the probability of failure N ij Shear components of stress from normal load T ij Shear components of stress from traction load z, Parameters to sample points on a sphere PAGE 19 19 Abst ract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy CRITICAL FLAW SIZE IN SILICON NITRIDE BALL BEARINGS By George Arthur Leves que December 2009 Chair: Nagaraj K. Arakere Major: Mechanical Engineering Aircraft engine and bearing manufacturers have been aggressively pursuing advanced materials technology systems solutions to meet main shaft bearing needs of advanced military air craft engines. Ceramic silicon nitride hybrid bearings are being developed for such high performance applications. Though silicon nitride exhibits many favorable properties such as high compressive strength, high hardness, a third of the density of steel, low coefficient of thermal expansion, and high corrosion and temperature resistance, they also have low fracture toughness and are susceptible to failure from fatigue spalls emanating from preexisting surface flaws that can grow under rolling contact fat igue (RCF). Rolling elements and raceways are among the most demanding components in aircraft engines due to a combination of high cyclic contact stresses, long expected component lifetimes, corrosive environment, and the high consequence of fatigue failure. The cost of these rolling elements increases exponentially with the decrease in allowable flaw size for service applications. Hence the range of 3D nonplanar surface flaw geometries subject to RCF is simulated to determine the critical flaw size (CFS) or the largest allowable flaw that does not grow under service conditions. PAGE 20 20 This dissertation is a numerical and experimental investigation of surface flaws in ceramic balls subjected to RCF and has resulted in the following analyses: Crack Shape Determ ination: the nucleation of surface flaws from ball impact that occurs during the manufacturing process is simulated. By examining the subsurface Hertzian stresses between contacting spheres, their applicability to predicting and characterizing crack size and shape is established. It is demonstrated that a wide range of cone and partial cone cracks, observed in practice, can be generated using the proposed approaches. RCF Simulation: the procedure and concerns in modeling nonplanar 3D cracks subject to RCF u sing FEA for stress intensity factor (SIF) trends observed from parametrically varying different physical effects are plotted and discussed Included are developments in contact algorithms for 3D nonplanar cracks, meshing of nonplanar cracks for SIFs, para metric studies via MATLAB and other subroutines in python and FORTRAN. Establishing Fracture Parameters: the fracture toughness, Kc, is determined by using numerical techniques on experimental tests namely the Brazilian di sc test and a novel compression test on an indented ball. The fatigue threshold for mixed mode loading, Keff, is determined by using a combination of numerical modeling and results from the V ring single ball RCF test CFS Determination: the range of 3D nonplanar surface flaw geometries s ubject to RCF are simulated to calculate mixed mode SIFs to determine the critical flaw size, or the larg est allowable flaw that does not grow under service conditions. The CFS results are presented as a function of Hertzian contact stress, traction magnit ude, and crack size. Empirical Equations: accurate empirical equations (response functions) for the KI, KII, and KIII SIFs for semi elliptical surface cracks subjected to RCF as a function of the contact patch diameter, angle of crack to the surface, max pressure, position along the crack front, and aspect ratio of the crack are developed via parametric 3D FEA. Statistical Probability of Failure: since the variability in mechanical properties for brittle materials is high a probabilistic investigation of variations in flaw size, flaw orientation, fracture toughness, and Hertzian load on failure probability is conducted to statistically determine the probability of ball failure for an existing flaw subjected to the service conditions. PAGE 21 21 CHAPTER 1 INTRODUCTION 1.1 General Background 1.1.1 Bearing Design Issues Advanced rotor systems for military aircraft engines are consistently being designed to operate at higher speeds and efficiencies. The bearings of these aircraft engines are often designed to their limits to meet the increased demands of thrust, speed, and weight. Developments in materials consistently push these limits and the case with hybrid bearings is no different. By combining steel raceways and silicon nitride (Si3N4) balls, researchers have establis hed that: Si3N4 balls can improve bearing fatigue life by up to 4 times than their M50 steel replacements and when they reach the end of their fatigue lives they fail in a similar fashion (Miner et al. 1996). Hybrid bearings display superior thermal chara cteristics to their all steel components (Miner et al. 1996). Hybrid bearings survived longer (by about 5 times) oil starvation periods than their all steel predecessors (Miner et al. 1996). Si3N4 exhibits favorable corrosion resistance (Klemm, 2002). T hough this hybrid system is an improvement over its all steel predecessors, the research and problems from the implementation are justifiably pursued considering the benefits that the balls provide. Silicon nitride (Si3N4) exhibits favorable material prope rties for application in hybrid high speed ball bearings such as high compressive strength, high hardness, third of the density of steel, low coefficient of thermal expansion, and high corrosion and temperature resistance. Its corrosion resistance and high hardness are good matches for new raceway materials that are very hard and corrosion resistant. Also, its corrosion resistance is beneficial for harsh environments, PAGE 22 22 like a NASA shuttle turbo pump (Miner et al. 1996). Comparison is drawn between it and it s steel replacement as rolling elements (see Table 11) and notes the high modulus, hardness and compressive strength that allow the ball to survive the loads and number of cycles that the engines are designed for. Silicon nitride is a ceramic that has gone through much development over the past few decades. Silicon nitride balls used in rollingelement bearings is an engineered ceramic product made from phase Si3N4 needle grain size held by a glassy phase. Balls for aerospace applications are made by sintering and hot isostatic pressing of Si3N4 grains with the glassy phase forming roughly 1% of the volume fraction with ball diameter and surface finish held to submicron levels (see Fig 11). The final ball dimension is achieved in a lapping process that produces (unavoidable) low velocity ball collisions resulting in surface cracks initiated by the radial tensile stress field at the contact periphery (Wang and Hadfield, 2000, Levesque and Arakere, 2008) (see Fig 12). 1.1.2 Issues W ith Silicon Nitride Even though there are many beneficial properties of silicon nitride ( Si3N4), it also has low fracture toughness (Piotrowski and OBrien, 2006). Hertzian fractures resulting from light collisions in the manufacturing proces s are cone cracks for normal interactions and, the more common, c cracks (or partial cone cracks) for oblique interactions (see Fig 12). These surface flaws would be subjected to rolling contact fatigue (RCF) if placed into service as a rolling contact el ement. These surface flaws often propagate to a fatigue spall and are the leading cause of ball failure (Hadfield, 1993b) (see Fig 13 ). Ball failure quickly results in raceway failure as the higher ball PAGE 23 23 hardness and rapid rotational speeds machine away th e raceway from a simple raceway spall (see Fig 14 ) to a damaged surface. In addition, brittle materials generally exhibit much variation for similar interactions. For example, there is a significant variation in the depth of a crack from a controlled imp act (Lawn 1967, Wang and Hadfield 2000). Also, there is a noted variation in the fracture toughness specimen to specimen (for example 4.85 0.36 MPa and OBrien, 2006). The material variations have a direct dependence on how each specimen is manufactured. While the flaws (illustrated in Fig. 12) can be readily induced and are generally present, methods for their detection prior to service is still in development. Development of a nondestructive evaluation (NDE) method is complicated because hairs) long on the surface and are difficult to find under a microscope. In addition, the material is nonconductive and only slightly translucent which complicates optical inspection procedures. From a fracture mechanics and structural integrity perspective, the larg est allowable surface flaw that does not propagate under RCF loading is of design significance, and is termed the critical flaw size. 1.2 Objective A nd Scope O f Research Work The main objective of the study presented in this thesis is utilize and develop numerical and analytical techniques to investigate general three dimensional (3D) nonplanar surface flaws subjected to the physical effects of RCF to determine their influence on the possibility of rolling element survival after the service conditions.. The scope of this study includes: (1) determining a surface flaws shape from known nucleation conditions, (2) simulating the range of possible flaws under RCF in finite PAGE 24 24 element analysis to determine how physical effects decrease fatigue life, (3) developing empirical equations to make these SIF results transferrable to future engineers (4) numerically analyzing experimental data to determine how three stress intensity factors (SIFs) are best combined into Kc or Keff parameters, and (5) exploring statistical methods to determine how variations in the RCF system affect ball survival probabilities. 1.3 Literature Review 1.3.1 Rolling Contact Fatigue Loading I n Ball Be arings The study of stresses resulting from contact has a long history. After the solutions for a point contact (Boussinesq, 1878) the stresses of two contacting elastic spheres were developed (Hertz, 1882a). For a thorough derivation of these theories an d other elastic contacts see Love, 1944. Since then there has been much development in the field of elastic contact (Johnson, 1987). In the case of ball bearings, the contact patch between the ball and raceway are ellipses (see Fig 15 ) and, in general, t he contact between any two objects that have two radii of curvature have elliptical contact patches. The stresses beneath an elliptical contact can be calculated by integrating derived equations (Hills et al., 1993) and has been analyzed previously in FEA (see Fig 16 ) (see Arakere et al., 2009). As both the elliptical loads (as each the inner and outer raceways have elliptical contacts on opposing ends of the ball) circumnavigate the ball, each point near the contact band experiences two cycles of fatigue and, as a result of possibly high operating speeds (see Table 1 2) high cycle fatigue threatens survival of these surfaces. This is rolling contact fatigue in action. The differences between rolling contact fatigue and standard fatigue have been previousl y listed in a literature review on rolling contact fatigue (Sadeghi et al., 2009) as: PAGE 25 25 1 The state of stress in nonconformal contacts where RCF occurs is complex and multi axial and governed by the Hertzian contact theory. 2 Contrary to most classical fatigue phenomena, rolling contact fatigue is typically a mu ltiaxial fatigue mechanism. 3 Contrary to classical fatigue, the loading history at a point below the surface is nonproportional, i.e., the stress components do not raise and fall with time in the same p roportion to each other. [For example], the stress history for a point located at the depth where the orthogonal shear stress xz is maximal [T]here is a complete reversal of the shear stress xz, while th e normal stresses x z always remain compressive. Also, the [maximal] peaks of the two normal stresses do not coincide [in time] with the peaks for the shear stress. 4 There is a high hydrostatic stress component present in the case of nonconformal contacts, which is absent in classical tensioncompression or bending fatigue. 5 The principal axes in non conformal contacts constantly change in direction during a stress cycle due to which the planes of maximum shear stress also keep changing. Thus, it is difficult to identify the planes where maximum fatigue damage occurs. 6 The phenomenon of RCF occurs in a very small volume of stressed material, because the contact stress field is highly localized. Typical bearing contact widths are of the order of 200 7 [In ductile materials] the evolution of RCF damage leading to a fatigue spall involves a threestage process: (i) shakedown, (ii) steady state elastic response, and (iii) instability. Localized plastic deformation and development of residual stresses are precursors to fatigue damage, and therefore the ability to compute the 3D elastic plastic stress fields that accounts for cyclic loading and traction effects, and acknowl edge microstructural changes, are necessary tool for quantifying raceway fatigue damage. Rolling contact fatigue in metal bearings is manifested as a flaking off of metallic particles from the surface of raceways and/or rolling elements. This process commences as a crack below the surface and is propagated to the surface, eventually forming a pit or spall in the loadcarrying surface. Lundberg and Palm gren (year) postulated that the maximum orthogonal shear stress is a subsurface crack initiation parameter (though not all re searchers accept that this i s the significant stress initiating failure) and they postulated that fatigue cracking commences at weak points subsurface PAGE 26 26 (Lundberg and Palmgren, 1947). However, modern (beyond 1990s) metallic bearing failures are frequently associated with surface initiated flaws. The LundbergPalmgren (LP) theory is widely used for fatigue life estimation in metal ball and roller bearings. The damage indicator in the LP theory, and life, can be written as: 3 0 0, 1 ln P C L Life Z V N s Parameter Damageh c e (1 1 ) Where N is the number of cycles endured with a probability of survival s 0 is the maximum value of the subsurface orthogonal shear stress, V is the volume of the subsurface material subjected to stress, z0 is the depth at which 0 occurs. L is the fatigue life in millions of revolutions, which 90% of the bearings will achieve; L10 or rating life. C is the basic dynamic capacity of the bearing as provided in the bearing manufacturers catalogs, calculated according to LP methods. P is the contact load. 1.3.2 Limitations O f L undbergPalmgren Theory The fatigue damage process in cerami c rolling elements is very different from metal bearings discussed above. Metals are weaker in shear than tension. In contrast to metallic materials Si3N4 material is weaker in tension than compression or shear The LP methods do not nor cannot account fo r this significant difference in material behavior. On the other hand, the fracture mechanics approach proposed is fully capable of modeling such material characteristics along with the complex systems design, lubrication and tribological interaction in be aring contact under steady state and adverse operating conditions. PAGE 27 27 The LP approach does not deal with fracture mechanics issues in predicting the subsurface initiated flaws. Also, it does not apply to this application as pre existing surface flaws in the ball grow under RCF and they are not initiated subsurface as the result of material cleanliness. It is akin to the stress life approach or the S N approach to fatigue life evaluation. Surface cracks and initial flaws also cannot be effectively addressed by the LP lifing, because of its probabilistic approach. A fracture mechanics methodology that evaluates the critical flaw size in balls and raceways is required. Failure typically occurs, for hybrid bearings, when the stress field at the crack tip equals th e fracture properties of the material. 1.3.3 Hertzian Fracture O f Brittle Materials Normal contact of spheres made of brittle materials has been analyzed extensively and is known to generate cone cracks (Chen et al. 1995, Lawn, 1994) (Fig. 12a). Normal i ndentation of a hard sphere into a brittle material typically results in a ring crack on the surface, which then propagates into a frustum of a cone. Auerbach (1891) and Hertz (1895) first investigated this problem in the 1890s and many have revisited the problem since. Lawn (1994) provides a comprehensive review of extensive investigations made subsequently, to elucidate a qualitative and quantitative description of the initiation and growth of cone cracks in brittle materials. Mackerle (2001) has provided a biography of indentation simulations and with an emphasis on ceramic materials (Mackerle, 2002). In contrast to these bodies of work, cracks generated by oblique contact have received much less attention. Partial cone or c cracks are the observed resul t of the oblique impact of brittle materials (Frank and Lawn, 1967) (Fig. 12b). These c cracks are considered the most damaging surface defect that limits ball life in hybrid bearings under RCF (Evans, 1983, PAGE 28 28 and Hadfield et al. 1993a). (See Fig 17 ). The se cracks, when subjected to RCF often undergo spallation which occasionally gives a clear image of the crack front and its interior shape (see Fig 18 ). C cracks are more commonly observed in sphereto sphere interactions and they not only have nonplanar crack faces but also possess non planar crack tips, as opposed to the axisymmetric cone cracks. This makes their shape more difficult to describe and, therefore, more difficult to analyze in any linear elastic fracture mechanics (LEFM) based analysis. 1. 3.4 Computational Simulation O f Surface Flaws Surface cracks are among the most common flaws (Hadfield 1993a) and most critical in structural components (Evans, 1993) and are directly exposed to RCF in multiple systems. The stress intensity factors (SIFs) are important for analysis of growth and fracture nucleating from these flaws. Due to the complexity of the three dimensional (3D) subsurface Hertzian stress field at the ball raceway interface and the steep stress gradients from the edge of contact to th e surface crack, exact solutions are intractable. Accounting for these complexities entails the use of a comprehensive 3D finite element contact/fracture mechanics analysis. Previous researchers have only provided a limited amount of approximations for SIF s for this system. In earlier works, where the investigated geometries had similar levels of complexity to the current work, FEA was not always used. Zalounia (1993) used the Fourier transforms method to obtain SIFs for an edge crack penetrating the interface in a coated solid subjected to contact loading. Keer and Bryant (1983) used a two dimensional linear elastic theory to simulate a crack in a rail wheel to examine the effects of contact friction, lubricant pressure and friction between the crack face s on the KI and KII stress intensity factors. Karapetian and Hanson (1994) derived weight PAGE 29 29 equations for the simpler geometry of a submerged circle crack subjected to point loads on the crack face. Later, these equations would be applied by Kida (2000) (wit h some untested approximations) to the case of a subsurface penny crack subjected to RCF at the surface. Hasebe and Qian (1995) visited the problem of an indenter of different geometries contacting a plane strain infinite plate with an angled surface flaw in multiple notable works by using complex stress functions derived by the rational mapping function (Hasebe, Qian, 1999, Qian, 1997, Qian, 1996). The results presented were in graphical form and not readily applicable to a 3D problem. Noda and Miyoshi ( 1996) used the body force method to come up with an integral equation that could be integrated with a polynomial stress distribution for a half penny crack in a semi infinite body. Pommier (1998) in a related work on semi elliptical cracks which are normal to the surface provided a set of equations for KI for a polynomial stress distribution that could be explained by the equation: mn xxxyab (1 2 ) w here the stress state is normal to the crack and n and m are constants whose sum is less than 4. Kaneta et al. (1991) investigated the penny crack by using the body force method. Their mode I SIFs occasionally went negative, due to a neglect of crack closure, when the contact patch started to go over the load. Kaneta et al. (1989) also visited the problem of an inclined subsurface penny crack subjected to line contact loading to investigate subsurface crack growth under RCF. FEA would also be used for analysis for a similar family of crack problems. Komvopoulos and Cho (1997) used two d imensional (2D) FEA to simulate a subsurface PAGE 30 30 plane strain crack that was parallel to the surface under sliding asperity contact. Zhang et al. (1999) also used FEA but for the analysis of multiple (two) vertical surface cracks under contact loading in 2D. K ojima (1999) used 2D FEA to analyze the angled surface crack under contact loading with a viscous lubricant to penetrate the crack. For the case of cone cracks with an internal circularly symmetric load, FEA was used to find the worst crack angle and the critical flaw size as determined by what causes the SIFs to reach the Keff to initiate fatigue (Warrier, 2000). Fletcher and Beynon (1999a, 1999b) applied their work to calculating SIFs for the specific case of contact loading for inclined cracks. The ir method for calculating SIFs involved edge Greens functions presented graphically in the works of Rooke et al. (1991) which must then be integrated based on the amount of crack opening. The method uses an approximation of an infinitely wide crack to cal culate SIFs for the deepest point of a semi circular surface crack. Problematically, this point may not be where the highest SIFs occur along the crack front depending on the depth of the crack (Newman and Raju, 1981). 1.3.5 Mixed M ode Fracture Testing Mixed mode fracture has been investigated by many researchers for finding a single parameter which is effective at determining fracture toughness, or Kc, or the mixed mode fracture toughness of a given material. In reality, there has been much divergence in t he equations provided for brittle materials which are often empirical in nature and both material and experiment dependent (Qian and Fatemi, 1996). While there are a variety of tests for fracture parameters this thesis contains analyses with three: the Br azilian disk test, the V ring test, and a novel fracture test from Michael J. OBrien (Piotrowski and OBrien, 2006). PAGE 31 31 The Brazilian disk (BD) test, first proposed by Carniero and Barcellos in 1953 has since been used to measure the tensile strength and fr acture toughness of brittle materials like rocks, concrete and ceramics. The BD test consists of a circular disk that, like most mixed mode fracture tests, have inserted planar flaws, by the process known as pre cracking, in the middle of these samples. Th ere are two widely used approaches for producing precracks in a specimen to determine its fracture toughness, Kc, either from surface flaws produced by indentations (where small flaws are needed) (like Ayatollahi and Aliha, 2005) or by using the chevron n otch, or a sharp edged notch, for larger cracks since there is stable crack growth during the initial precracking (like Fowell, 1995). For silicon nitride, the high loads that are required to grow a centrally placed flaw under testing conditions require a larger crack to be capable by common load frames and since the material is very tough to cut or grind (a favorable material property for the harsh conditions it will face) a Chevron notch is not favorable. The BD test consists of a precracked circular disk loaded in compression that results in a tensile stress at the center (perpendicular to the loading direction) and can spl it the disk along the load axis (Canerio, 1953) With the crack size visible at the end of the experiment, and the final load to fr acture known from the load frame, fracture toughness can be calculated (Piotrowski and OBrien, 2006). The advantage of this test is that it can be performed under a range of modemixities (from pure mode I to a mix of mode I and II for through cracks) to evaluate a mixed mode fracture toughness along with simple specimen geometry and minimal requirements for implementation. Brazilian disk specimen can also be used to measure interfacial fracture toughness in biomaterials (Banks Sills, 1999). PAGE 32 32 There have b een quite a number of Brazilian disc studies conducted on fracture toughness via multiple approaches on different materials. Zhou et al. (2006) tested Polymethylmethacrylate (PMMA) with a chevron notch to determine its mode I fracture toughness in the BD t est. Tong, et al. (2007) used the Brazilian disk test to determine the interfacial fracture toughness in bone cement interfaces. Awaji and Sato (1978) employed this test for examining the fracture toughness of graphite, plaster and marble in mixedmode loading using machined central cracks. Petrovic and Mendiratta (1976) have examined mixedmode fracture toughness in hot pressed silicon nitride using controlled surface flaws produced by Knoop indentation in a four point bending test. Khandelwal, et al. (199 5) have found the high temperature mixed mode fracture toughness of hot isostically pressed, PY6 silicon nitride using bend bars which contained cracked indentations produced using the Vickers hardness indenter and tested using the four point bend test. Be yond the BD test, there is a group at Bournemouth University that uses a four ball test rig to test rolling contact fatigue and it has become a centerpiece of their publications (Wang and Hadfield, 2000). In this test a single ball is rotated axially in a collet and contacts three balls below which rotate in a cup. Surface flaws are induced in the upper ball and are carefully placed in the contact path such that on each upper ball rotation the flaw goes through three RCF cycles as a result of touching the t hree balls below. This test and others are summarized in table 13. 1.4 State O f T he Art F rom Literature Survey C cracks or partial cone cracks are commonly found on the surface of silicon nitride balls and, if subjected to rolling contact fatigue can result in spallation failure. The Hertzian fracture process has been largely ignored for these partial cone cracks even PAGE 33 33 though they are more common than the symmetrical cone crack. Numerical simulation of 3D nonplanar flaw for SIF calculation is also a relativ ely unanalyzed field (especially near an elliptical bearing load). A Keff parameter, while researched by many, has left no consensus for a general case or material type and, even in specific analyses, empirical equations are often provided. Probability of bearing survival has only been looked at by empirical equations for steels rather than a fracture mechanics approach for pre existing flaws on the ball surface. 1.5 Outline O f Thesis In chapter 2, subsurface Hertzian stresses between contacting spheres ar e examined, using an analytical stress solution, to investigate their applicability to predicting and characterizing crack size and shape. Also, these cracks are incrementally developed through an iterative crack growth procedure using a 3D finite element analysis. Comparisons are then made to experimental images of the flaws in silicon nitride. By varying the initial conditions during the contact interaction of the balls it can be demonstrated that a wide range of cone and partial cone cracks, observed in practice, can be generated using both the analytical and numerical fracture mechanics approaches. Furthermore, an expression is presented for the impact velocity that induces a cone crack from a maximum radial stress criterion at the contact periphery. In chapter 3, the researchers discuss the modeling techniques that were developed to simulate the resultant types of flaws. Included are developments in contact algorithms for 3D curvilinear cracks, meshing of curvilinear cracks in 3D, parametric studies via MATLAB and more. In addition, SIF trends observed from parametrically varying different physical effects are plotted and discussed. PAGE 34 34 In chapter 5 the range of threedimensional non planar surface flaw geometries subject to RCF is simulated to calculate m ixed mode stress intensity factors to determine the critical flaw size, or the larg est allowable flaw that does not grow under service conditions. The cost of nondestructive evaluation methods for silicon nitride balls scales exponentially with decreasing CFS and increasing ball diameter and can become a significant fraction of the overall manufacturing cost. Stress intensity factor variability is analyzed for variations of the location and orientation of the load relative to the crack, the geometry of load, and full slip traction. The modeling techniques utilized in the creation of a three dimensional FEA model is discussed and the maximum tensile contact periphery stress is examined for effect on crack driving force under RCF. The CFS results are presented as a function of Hertzian contact stress, traction magnitude, and crack size. In chapter 4 Kc and Keff are determined via the experimental/numerical analysis of three tests: Brazilian disc test, a novel test from Michael J. OBrien, and the V ring test. Mixed mode fracture toughness parameters are tested for applicability via the Brazilian disc test on silicon nitride. Specimens are precracked, measured, loaded to fracture and then simulated individually in finite element analysis to calculate the speci fic stress intensity factors under failure. Avoiding empirical equations, the max strain energy release rate criterion is found to be the easiest and acceptably accurate for mixed mode conditions. In a novel test from Michael J. OBrien, full silicon nitr ide balls are tested by indenting them and compressing them in a load frame until the tension developed at the equator propagates the cracks extending from the Vickers indentation. Empirical equations are then used to calculate the fracture toughness of ea ch specimen. The V  PAGE 35 35 ring test is the only test utilized to calculate a Keff. Herein, to opposing inner rings squeeze a rotating ball at a certain cyclic rate and load. An assortment of flaws are measured and placed in the contact path to interact with bot h the upper and lower rings. The test is stopped when a fatigue spall trips an accelerometer or a certain number of cycles is reached. E ach tests dimensions and orientation are modeled to calculate SIFs associated for tests that fail and do not fail to de termine what SIFs are sufficient for a ball failure. Chapter 6 contains empirical equations for the KI, KII, and KIII stress intensity factors for semi elliptical surface cracks subjected to rolling contact fatigue as a function of the contact patch diamet er, angle of crack to the surface, max pressure, position along the crack front, and aspect ratio of the crack. The equations were developed from SIFs calculated by parametric 3D finite element evaluation for a range of contact patch radii (1 b 2 b and 3b ) and angles of the crack to the surface (0, 45 and 60). The comprehensive empirical curve fits presented are accurate to within 0.5% of FE simulations. The results are of particular relevance to hybrid silicon nitride ball bearings which are susceptible to failure from fatigue spalls emanating from preexisting surface cracks, due to crack growth driven by RCF Chapter 7 contains an uncertainty analysis of rolling element survival under rolling contact fatigue. Previous chapters focused on a determinist ic categorization of allowable flaw size. This chapter determines ball survivability if there is a variation in tolerable crack size, fatigue threshold, and crack orientation on the ball. First we create a surrogate modeling based on the physics of the S IFs under rolling contact fatigue to account for variations in crack size, max pressure, and orientation and then conduct a PAGE 36 36 Monte Carlo Simulation to determine the effect of crack size and fatigue threshold on ball survivability. PAGE 37 37 1.6 Figures Table 1 1. C omparison of material properties between fully densified silicon nitride and a common ball material predecessor M50 steel. Parameter M50 Steel Silicon Nitride Specific Gravity 7. 8 3.25 Youngs modulus (GPa) 207 310 (ksi) 30 45 Vickers hardness (GPa) 8 14 18 Coeff. of Thermal Expansion (x10^ 6/K) 12 3 (x10^ 6/F) 6.7 7.1 Specific Heat (J/kg K) 450 800 (Btu/hr ft F) 0.11 0.19 Thermal Conductivity (W/m K) 30 20 (Btu/hr ft F) 17.3 11.5 Upper use temp erature (K) 600 1300 (F) 620 1800 Fracture Toughness (MPa m^0.5) 18 6 *Note: Toughness from (Rescalvo, 1998) and (Piotrowski, 2006). All other properties (Galbato, 1992). Figure 1grains are apparent [Image courtesy of The Timken Company ]. PAGE 38 38 Figure 12. A) Example of a cone crack in silicon nitride (Mecholsky, 2008), B) Experimental image of oblique sphere interaction with contact patch and resultant c crack shown. Image in proportion but scales are not shown ( Image courtesy of The Timken Company ). (a) (b) PAGE 39 39 Figure 13 Before and after of Timken test #7 of the V ring s ingle ball test. PAGE 40 40 Figure 14 Image of steel raceway spall from RCF (Branch et al., 2009). PAGE 41 41 Table 1 2. Measurements for a sense of scale of the problem (Miner et al., 1996). Feature Measurement Size of Balls 14.29 25.4 mm (9/16 1 in) Crack Width on Surface 300 Grain Length Speed 11000 rpm Radial Load 9000 lb Axial Load 1300 lb p o 3800 MPa Table 1 3. RCF tests on ceramics and their failure modes in literature. Publication Test Failure Parker and Zaretsky, 1975 Five ball test on 12.7mm dia. HIPed Si3N4 at 4.3 6.2 GPa, 9600rpm, and 30 contact angle Spallation without edge cracking Lucek and Crowley, 1978 Disc on rod on hot pressed ceramic 4.1 5.5GPa Non catastrophic spallation failure. Hadfield and Stolarski, 1995 Dis c on rod test of varying ceramics and lubricants Non catastrophic spallation failure. Morrison et al. 1984 Si 3 N 4 hybrid bearings with M50 steel between 1.952.44 GPa Spallation with no cases of ball fracture. Fujiwara et al. 1989 Ball on plate machin e to test Si3N4 raceways up to 6.4 GPa at 1400rpm. Spalling is most common long term failure but cavein and peeling also occurred in short tests. Lucek, 1990 Hot pressed Si 3 N 4 rods at 6.4GPa up to 8600 rpm Spallation occurred Hadfield, 1993a,b,c Four ball test at 6.4 GPa at 5000 rpm with different lubricants Non catastrophic spallation Burrier, 1996 Ball on rod test of 11 Si 3 N 4 materials at 5.93 GPa at 3600rpm Large life dependency on finer microstructure with primary failure of spallation. PAGE 42 42 Figu re 15 The elliptical contact patches between the ball and raceways. The contact patch size and the subsurface stress distribution are a function of raceway geometry, conformity between the ball and raceways, radial and thrust loads, etc. PAGE 43 43 Figure 16 FEA simulation of a steel raceway section being subjected to contact with a sphere from a distance and with a cutaway section revealing the subsurface stresses. PAGE 44 4 4 Figure 17 Fracture stress as a function of dominant defect dimension (Evans, 1983). Figure 18 Image of a spalled c crack clearly showing its face as it extends into the material. Other markers are indicating possible dimensions of measurement for further classification and analysis. PAGE 45 45 CHAPTER 2 DETERMINATION OF SURFACE FLAW GEOMETRY 2 .1 Introduction In this chapter, the goal is to obtain the possible range of crack shapes in three dimensions that result from oblique interaction for future finite element (FE) analysis. The analysis of crack shapes produced by oblique interactions has di rect and immediate application. The performance of hybrid silicon nitride ball/steel raceway bearings has been shown to be much superior to steel bearings (Miner et al. 1996 and Tanimoto et al. 2000). However, silicon nitride balls are sensitive to surfa ce defects and can fail from fatigue spalls emanating from preexisting c cracks, due to crack growth driven by RCF (Hadfield et al. 1993b). Silicon nitride (Si3N4) exhibits favorable material properties for application in hybrid high speed ball bearings such as high compressive strength, high hardness, a third of the density of steel, low coefficient of thermal expansion, and high corrosion and temperature resistance. However, it also has low fracture toughness (Piotrowski and OBrien, 2006). The final b all dimension is achieved in a lapping process that produces (unavoidable) low velocity ball collisions resulting in surface cracks initiated by the radial tensile stress field at the contact periphery (Wang and Hadfield, 2000). A qualitative and quantitat ive description is presented on the range of possible partial cone crack shapes in 3D as functions of initial conditions (coefficient of friction and contact patch size) during oblique spherical contact. The emphasis is not on simulating the dynamics of oblique contacts, but simulating the interaction forces via the coefficient of friction and contact patch size, and then focusing on the effect of resulting subsurface Hertzian stress fields in inducing brittle fracture. The resulting 3D geometry PAGE 46 46 of crack faces is predicted using the following two approaches: 1) An analytical description of 3D subsurface stress state is used to predict the crack front based on the max principal stress trajectory and 2) A 3D finite element contact stress analysis is used in conjunction with a numerical 3D fracture analysis program, FRANC3D, to grow the partial cone crack incrementally. Results from the two methods are compared and show good agreement. Qualitative features are examined from experimentally produced C cracks in silicon nitride balls. By varying the initial conditions of the contact interaction it is shown that a wide range of partial cone cracks, observed in practice, can be generated using both the analytical and numerical fracture mechanics approaches. The resu lts generated are of immediate engineering relevance to the hybrid ball bearing industry towards evaluating critical flaw size for defining limits for non destructive evaluation methods for silicon nitride ball quality control and for developing a fracture mechanics based life prediction methodology for hybrid bearings. 2.2 Analysis Here the stress state and fracture mechanics approaches mentioned earlier are discussed to yield the range of partial cone crack shapes produced during oblique interaction. The stressstate analyses have been done for normal indentation for some time. The cone crack path due to a spherical indenter was shown by Frank and Lawn (1967), to a first approximation, to be orthogonal to the maximum tensile principal stress distribution o f the statically loaded elastic half space. Later, Lawn (1967) showed that traction forces induced by a sliding spherical indenter strongly influenced the quasi static stress field, and generated an early partial cone crack profile for brittle materials in two dimensions. After the works of Frank and Lawn (1967), there was a continuing use of the max principal stress trajectory in cone crack analyses. As a result of this history, PAGE 47 47 this body of work will extend the stress state analysis to describe the crack shapes in threedimensions and discuss a couple variations. Crack shape analysis has also been done using an incremental growth approach as endorsed by Kocer and Collins (1998). In later work on cone cracks, Kocer emphasizes that the study of stress in an uncracked body to determine a crack shape had not been analytically supported to be reminiscent of the resultant crack shape (as once a crack exists, the stress field changes significantly and the crack grows in compliance with this new stress field) and s o believed that the crack shape should be acquired based upon an iterative crack growth analysis (or an experimental approach). An incremental crack growth analysis is conducted, and in the case of the partial cone cracks, an iterative crack growth study c an only be done correctly in threedimensions. Doing such a fracture mechanics analysis using 3D FEA requires the ability to incrementally grow cracks and remesh, resulting in a computationally intensive endeavor. T hese two analyses are compared with experimental images that were obtained through multiple sources. The surface dimension of c cracks typically seen in Si3N4 balls ranges from 50 400 m, depending on the ball diameter (Wang and Hadfield, 2000). Experimental examination of crack shape of these v ery small cracks that are embedded in a nontranslucent material like Si3N4 which is difficult to cut or grind, and hence poses considerable difficulty. 2.2.1 Stress State Analysis To attempt to generate the crack shape from its uncracked stress state, the stress state of obliquely interacting spheres must first be obtained. Oblique impact with contact friction gives rise to tangential tractions potentially leading to micro slip at the interface. PAGE 48 48 The complex interaction between tangential traction and micro slip and their time history dependence has been visited by Mindlin and Deresiewicz (1953) and analyzed by Maw et al. (1976, 1981). However, these surface tractions are not closedform and more difficult to incorporate into a 3D subsurface stress field sol ution and are not worth the possible loss in accuracy (Andersson, 1996). Also, balls are typically well lubricated during the lapping process where these cracks originate. With this support a sliding contact analysis is conducted. To investigate cone and c crack size and shape as a function of normal and traction force distributions, a comprehensive 3D analytical subsurface stress solution is used for efficiency. The Hertz theory of quasi static impact holds for evaluating the variation of contact size and contact pressure during impact under familiar criteria of ball velocities relative to wave speeds (Johnson, 1987). Applying superposition for normal and tangential stresses, the Hertz normal pressure is given by 2()1oprpr (2 1) where po is the peak pressure and r is the cylindrical coordinate. For a contact in full slip the tangential traction distribution is ) ( *r p where is the coefficient of friction. Hamilton and Goodman (1966) presented the 3D subsurfa ce stress distribution for a spherical sliding contact. Solutions are given elsewhere (Hamilton, 1983, Sackfield and Hills, 1983, and Hills et al. 1992). The Cartesian coordinate system is used in the configuration displayed in Fig. 22. To acquire a sol ution, the harmonic potentials, N and T, are chosen to satisfy the necessary boundary conditions, namely: zz() zxprra f (2 2a) PAGE 49 49 0 zzzxra (2 2b) So they may be substituted into the relations for displacements u, v, and w, (Love, 1921) namely: ,,,,,2(12) 22xxzxxzzxxzuNzNTTzT (2 3a) ,,,,2(12) 2yyzxyxyzvNzNTzT (2 3b) ,,,,22(1) (12)zzz zxxzzwNzNTzT (2 3c) where commas denote derivatives and is the Poissons ratio. Solutions are given elsew here (Hamilton, 1983, Sackfield and Hills, 1983, and Hills et al. 1992) and are reproduced below for completeness. 3 22 43 3 2 24 2 2 22 2421 12 1 3 1 1arctan21 12 1 1 1N xx oz vxy pru z u x xu vu vv v uuu ru uuz (2 4a) 3 22 43 2 2 24 2 2 22 2421 12 1 3 1 1arctan21 12 1 1 1N yy oz vxy pru z u y yu vu vv v uuu ru uuz (2 4b) 3 2 24212 1N xy oxyzu vJ p uuz (2 4c) 3 42 N zz oz p uuz (2 4d) 2 2421N zx oxzu p uuz (2 4e) PAGE 50 50 2 2421N yz oyzu p uuz (2 4f) 21 1arctan 1T T T yy xx zz ooou xv fpfpfp u u (2 4g) 25 23 2 2 242313 2arctan 881 41 1 12T yy ouuyu vx fp u u u uuz J vz y (2 4h) 2 2421T zz oxzufp uuz (2 4i) 3 2 2421T yz oxyzu fp uuz (2 4j) 23 2 2 242311 arctan 2 21 1T zx ouxu z fp uu u uuz (2 4k) 25 23 2 2242 211 2arctan 8 81 411 1 arctan 1221T xy ouu xu vy fp u u uuuz yuJ vz uux (2 4l) where 3 432 3 xyzz J ruu (2 5) and 2 2222221 114 2 urzrzz (2 6) PAGE 51 51 The formulation used here, for spheres in sliding contact, is nondimensionalized such that th e only inputs required are a Poissons ratio ( ), coefficient of friction ( ), and magnitude of po. The stress has been normalized with respect to (w.r.t.) po and the contact patch is normalized w.r.t. itself and always is at r=1. Fig. 2 3 is an example contour plot. At points close to the surface, coordinate planes and axes, and at r=1, the limits of these functions were evaluated for efficiency and to avoid numerical difficulties. With a 3D stress field well described, it is po ssible to generate the shape of the c cracks according to the max principal stress trajectory for brittle materials by Frank and Lawn (1967). When max principal stress trajectory was used by Frank and Lawn, they left the phrase largely undefined. Generally, trajectory refers to some curve that intersects a family of surfaces at the same angle. However, for any given angle (say 90 degrees) there are an infinite number of surfaces that can cut through the contours of max principal stress values. It can b e seen that they implied to generate a line that cut through the contours at 90 degrees and originated at the surface, near the edge of the contact patch, and thus followed something close to the direction of maximum change in the field. Doing so computati onally in a 3D field can quickly become difficult as gradients of stress become demanding depending on the level of refinement necessary (especially near the surface) (no matter if the stress state is being evaluated using equations or the finite element m ethod). Also, as the point of crack origin is not always in agreement to what is seen in experiment (Kocer and Collins, 1998) some perturbations about the contact periphery are necessary to see the changes in crack shape for a given stress field. PAGE 52 52 A trajec tory approach, like Frank and Lawn (1967) was extended to three dimensions, where the 3D contours of the third (minimum) principal stress are plotted according to 33(1,0,0)(,,) xyz where is a chosen small value that force s the contour value of interest to be near the periphery of the contact. This method functions well because it is orthonormal at the intersections of contours of the max principal stress and has a similar result to searching out the direction of max change in the max principal stress. However, it should be noted that the choice of epsilon is significant since the stresses on the periphery of the contact is quite rapid and two close adjacent contours will diverge from each other as distance from the contact periphery increases. The cracks displayed in Fig. 24 are examples of this procedure. The method is effective in generating cone cracks for normal contact and partial cone cracks for oblique contact. The generated cone angle of 63 to the vertical agrees w ith experimental observations (Kocer and Collins, 1998). Briefly, similar crack shapes we re also generated using two other techniques. For example, a maximum t ensile stress approach, where the 3D stress field is searched for points which are maxima in the z direction according to 1(,,) 0 xyz z and these points are then connected to form smooth, continuous plane originating at the contact periphery. Also, the plane of maximum change where the surfaces are produced by local maxima is investigated. [When investigating local maxima, the investigator is forced to choose a coordinate variable in which the max principal stress undergo maximum change. Here it is noted that the surfaces produced by 10 r (for cone PAGE 53 53 cracks) and 10 x (for c cracks) agree best with experimental observation.] These two methods were seen to be similar when the derivatives were taken with respect to the mentioned variables. Since the contact patch radius has been normalized w.r.t. itse lf, the crack shapes generated are attributable to all cracks that are generated by contacting spheres with the prescribed Poissons ratio (which produces little effect when analyzing the parametric variation) and coefficient of friction but its size must be scaled to the crack generating contact patch. The generated crack surface is smoothed by a moving average on a discrete number of evaluation points [by the smooth() function in Matlab (Mathworks, 2005)] and then patched for a 3D plane [by the patch() function (Mathworks, 2005)]. Since the surfaces produced by the above procedure are based on stress distributions alone, they can differ from the physical cracks whose extension is limited by the physics of dynamic brittle fracture. A stressanalysis based criterion to determine the crack depth can be explored in future work For this work, surfaces were curtailed to show common cracking proportions observed through experiments. 2.2.2 Incremental Crack Growth Analysis Also generated was the 3D partial cone crack using the finite element method and an incremental crack growth procedure. Starting with an initial small penny shaped crack, the crack was then incrementally grown using the crack growth software FRANC3D/NG developed by the Fracture Analysis Consul tants (2005) after being analyzed by the commercial finite element software ABAQUS (Dessault Systmes, 2007). While FRANC3D/NG is normally used for fatigue crack growth, its abilities to incrementally grow and remesh cracks to simulate c crack creation ar e used. PAGE 54 54 In 3D FEA, a small vertical half penny flaw ( a b silicon nitride block ( E =310 GPa, v =0.28). The body was meshed with a combination of quadratic hexahedral, pyramidal, and tetrahedral elements (so that the crack had collapsed hexahedral quarter point elements to preserve the singularity and avoid warning elements), which totaled around 25,000 elements. This part was then loaded according to Fig. 25 with an applied Hertz pressure distribution, as analyzed with the s tress analysis technique and is fixed far field. Far field boundary conditions were simulated on this small cracked block using a submodeling technique (Dessault Systmes, 2007) where displacements are interpolated from a larger block undergoing the same load orientation and is fixed at its furthest side. To grow a crack after the displaced state is calculated, stress intensity factors are calculated by the displacement correlation equations: 12 22 4(1)IE K uu vr (2 7a) 12 22 4(1)IIE K vv vr (2 7b) 12 22 4(1)IIIE K ww vr (2 7c) The computed stress intensity factors are inserted into the maximum stress crack turning angle criteria and then the crack is extended an average length as specified by the user according to a relative extension power law written as: n nodei nodeimean meanK aa K (2 8) PAGE 55 55 where the exponent, n, is iteratively chosen to avoid high gradients in the stress intensity factors along the crack front and the average length, meana was often chosen to be about 10 um. Incremental remeshing w as stopped once the relative proportions of cracks observed in experiment were reached. Images of incremental growth are shown for =0.1 and =0.5 in Fig. 26 and for =0.5 in Fig. 27. 2.3 Result s O n Crack Shape Determination Three dimensional crack shapes resulting from the stress state analysis exhibit the following characteristics: (a) the shape of the cone crack is expectedly axisymmetric and exhibits the experimentally observed angle of appr oximately 63 to the vertical (Frank and Lawn, 1967) and (b) when spheres collide obliquely the resulting max tensile stress is in the wake of the contact and is the site of c crack origination. From the stress state examination, it is observed that a c cr ack is unlikely to circumvent more than 180 about the contact patch since the other side of the contact region is in compression. Research on silicon nitride confirms that their circumference is between a third and a fourth of a complete circle (Hadfield et al. 1993b). Also, c cracks tend to depart from the contact patch radii as friction increases (see Fig. 21b). Furthermore, a plot of the c cracks produced as a result of different friction coefficients (Fig. 28) shows that the range of c crack angles observed (Zhao et al. 2006) are accounted for and that higher coefficients of friction increase the radius of the c crack on the surface and steeper cracks. The results are summarized in table 21. C ross sectional images produced by the methods are compared to those examined in experiment. Firstly, an array of cracks from Zhao et al. (2006) is shown which is a collection of naturally occurring (those which result from the manufacturing PAGE 56 56 process) c cracks in Si3N4 with a range in steepness (Fig. 29). A ls o display ed are images of artificially produced cracks (created in laboratory settings on finished balls) in Fig. 2 10 ( courtesy of The Timken Company ). Reflecting on the incremental growth results it is observed that the range of crack shapes produced by the two analyses are quite similar and agree with the range of shapes that are observed in experiment. The incremental growth analysis may produce small perturbations from the actual crack shape that can expectedly result from an incremental remeshing te chnique. The incremental growth analysis also tends to have more curvature as it extends into the material as opposed to the stress state analysis. Also, the incremental results seem to conform to the contact patch more than the stress results that may b e the result of the incremental approach. A superposed image of the incrementally grown crack was developed by importing the global coordinates of the created crack into Matlab and performing Delaunay triangulation (Mathworks, 2005) to obtain a series of triangles that were plotted as a triangular surface plot. The result is reminiscent of the actual mesh (see Fig.2 11). Comparing the two shapes, they are in excellent agreement on the midsection where increments were comparably smaller but were somewhat divergent on the surface. As a result, it can be seen the benefits of both approaches. The incremental crack growth analysis is well supported by a history of LEFM based analysis but regions of significant crack turning may take significant amounts of com putational effort as the size of the meana must be decreased. The incremental growth analysis also produces a meshed part that the user can apply under any other type of loading for other analyses. The stress state analysis is relati vely fast to implement and runs without user PAGE 57 57 interaction. In addition, their good agreement may leave both as viable options for future analysts depending on the analysts time and needs. 2.4 Conclusion s This analysis of oblique contacting spheres lead to the following conclusions: 1 Three approaches are presented that operate on the 3D subsurface stress field and effectively predict potential crack shapes for both the axisymmetric cone crack and the oblique resultant c crack for Hertzian fracture. 2 An iter ative, 3D crack growth simulation can also be used to generate the potential crack shapes. 3 B oth methods to agree well with each other and experimental images. 4 The methods can predict a wide range of crack shapes observed in oblique interactions by varying the friction coefficient at the contact. 5 With knowledge derived from this investigation it is possible to approximate a maximum periphery stress to induce cracking, an equation to determine the velocity required to induce cracking, and an approximation of the original flaw size from which cracks may nucleate. 6 An equation was derived for calculation of the maximum tensile periphery stress that drives crack nucleation, including the effect of friction. While it was not originally intended for these analyses to be deterministic of crack shapes from the initial conditions of colliding spheres, results have indicated the expected range of crack shapes in three separate ways: experimental reflection, a stressstate analysis and an iterative growth technique. Wi th this family of crack shapes, one conducting an analysis on these types of cracks should consider the range of steepness and proportion produced by this analysis. Since the range of crack shapes produced is observed in experiment, an examination of the u ncracked stress field can yield accurate predictions of crack shapes for brittle materials and propose the method to be used in other instances where a stress field is known (even by computational methods) and the stress field is torsion free PAGE 58 58 (Frank and Lawn, 1967) as long as the researchers note the possible variation of traction in the (small) contact region. Also, these results are in excellent agreement with the shapes produced by an iterative remeshing technique that requires the user to have an adequate mesh and a small increment of crack advance. As an added benefit, the iterative analysis creates a mesh suitable for analysis of c cracks under RCF. PAGE 59 59 2.5 Figures Figure 21. A) Example of a cone crack in silicon nitride (Mecholsky, 2008), B) Experimental image of oblique sphere interaction with contact patch and resultant c crack shown. Image in proportion but scales are not shown ( Image courtesy of The Timken Company ). (a) (b) PAGE 60 60 Figure 22. Diagram of coordinate system for the equ ations of stress for oblique colliding spheres. PAGE 61 61 Figure 23. Surfaces representing constant values of the max principal stress field. PAGE 62 62 Figure 24. A) The cone crack as generated from the Hertzian stresses without traction with gradients of nondimensional max principal stress. B) A c crack generated for v =0.3 and =0.1. Table 2 1. Estimated c crack geometries produced by different values of from stress state analysis Crack Radius Angle to ball surface 0.1 1.07 48 0.2 1.12 54 0.3 1.17 63 0.4 1.22 72 0.5 1.28 77 Note: Crack radius is on the surface (nond imensionalized w.r.t. the contact patch) and the estimate of angle in degrees is measured down from the surface PAGE 63 63 Figure 25. Initial configuration of crack and load at the beginning of the crack growth analysis in the FE model. PAGE 64 64 Figure 26. Incremental images of crack growth for =0.1. a) half width=0.04a depth=0.07 a b) half width=0.07 a depth=0.14a c) half width=0 .09 a depth=0.16 a d) half width=0 .11 a depth=0.32 a e) half width= 0.12 a depth=0.41 a f) half width =0.15 a depth=0.58a PAGE 65 65 Figure 27. Incremental images of crack growth fo r =0.5. a) half width=0.11a depth=0.17 a b) half width=0.17 a depth=0.35a c) half width=0.29 a depth=0.69 a d) half width=0.39 a depth=1.07a e) half width=0.523a depth= 1.46a f) half width=0.59 a depth=1.59 a PAGE 66 66 Figure 28. Cracks produced for = 0.1, 0.2, 0.3, 0.4, and 0.5 and v =0.28. Note the increasing surface radius and steepness with the increase of the friction coefficient. PAGE 67 67 Figure 29. Cross sections of multiple natural c cracks on 12.7 mm diameter, grade 5 balls that (after manufacturing) h et al. (2006) for details. PAGE 68 68 Figure 210. Cross sections of artificially produced c cracks. Images in proportion but scales are not shown [Image courtesy of The Timken Company] PAGE 69 69 Figure 2 11. For =0.1, fro m the front a) incremental growth crack, b) stress state crack. Superposed images of half cracks produced by the stress state (in blue) and the incremental growth (in red) at c) the cracks midsection and d) from an az i muthal angle of 45 and a elevation angle of 18 as defined by matlab (Mathworks, 2005). PAGE 70 70 CHAPTER 3 GENERAL MODELING ISSUES OF SURFACE FLAWS SUBJECT TO R OLLING CONTACT FATIGUE 3.1 Introduction The partial cone crack has a very complicated shape as illustrated in chapter 2. It is not planar, like a semi elliptical flaw, nor axisymmetric, like a cone crack. Modeling the shape accurately in order to calculate accurate SIFs for the geometry has been a focus of this research. A few of techniques to acquire the desired shape have been used where t he first technique was dependent upon the complex use of spl ines to generate a 3D surface (see Fig 31). In general, the simulation of a semi elliptical crack is an advanced application of FEA code. In addition to the fact that the part has highly curved surfaces, it is also a closed geometry (two faces of the same body are immediately adjacent to each other) and this complicates meshing. The proper mesh density for a surface crack under RCF is really dependent on the purpose of the model. For the purpose of calculating accurate stresses, the user must not only have accurate enough meshes to account for the high stress gradients near the contact periphery but also to account for the 1/ r singularity in the stress field because of the presence of the crack in the material. For the purpose of calculating SIFs, the mesh can be coarser but if the crack faces close than the contact calculated between them needs to be exceptionally accurate and this demands more refinement. The models boundary conditions are also quite complicated. There is contact between a ball and raceway (of different curvatures) whose interfacial traction is unknown. The contact patch is very small relative to the size o f the entire cracked body and an entire cracked ball is not simulated, but a large portion of the ball whose edge PAGE 71 71 effects do not contribute to the overall analysis. In addition, the contact itself has high stress gradients near the contact periphery which should be properly accounted for with mesh density if stress results are required. At the same time, many look into the problem of cracks in rolling contact fatigue and attribute the growth mechanism to fluid penetrating the crack and the crack face shutti ng down, resulting in an interior pressure on the crack faces. This is neglected here due to the complex modeling issues of the fluid structure interaction problem. 3.2 Computing S tress I ntensity F actor s For this analysis, computing stress intensity facto rs ( SIFs ) from the crack tip is an additional difficulty. By s imply applying the cracked modeled specimen to tensile loading ( or mode I loading) (see fig 34 ) will result in all three modes being present at the crack tip (see fig 35 ). This is due to the c omplex shape of the crack which can turn even the simplest loading into complex stress intensity factors. Whichever method of SIF calculation is chosen, it must be able to calculate for all three modes of crack tip displacement. T here are a few methods f or SIF computation, but the main focus herein is on J integral computation and displacement correlation. ABAQUS has a built in J integral decomposition tool in order to extra ct SIFs where they are needed, h owever, this tool is not complete. Its formulation i s based on energy methods and, even though ABAQUS is completely capable of solving a multiphysics problem, it omits multiple physical contributions (like heat, dynamic behavior, etc.) in SIF calculation. It even neglects a component of the traction term seen in PAGE 72 72 i Riu JwdyTds x (3 1) which arises due to crack face contact (see fig 3 1). So even if the displacement of every node is computationally accurate, the final SIFs will be as though no contact was defined. Since the model is often in a state of compression somewhere, crack face contact needs to be properly accounted for in an external SIF calculation. As a result, displacement correlation became easier to apply in this regard. No matter which physical interactions are influencing the crack opening, if the displacements are correct, so are the calculated SIFs. The theory for all three modes of crack opening can be written as: 12 22 4(1)I dE K uu r (3 2) 12 22 4(1)II dE K vv r (3 3) 12 22 4(1)III dE K ww r (3 4) N ote that rd is the radius from where the relative displacements are measured from the crack tip, and u v and w correspond to the displacements in the x y and z directions for the local coordinate system at the tip of the crack (see Fig 33). As of this writing, displacement correlation is not available in ABAQUS (v.6.71) and so an external code must be written. ABAQUS Python libraries have been used to construct a post processing routine. Since the crack face and the crack tip are nonplanar, the coordinate system for displacement correlation must be calculated at each point SIFs PAGE 73 73 are extracted along the crack front. To illustrate the complexity of the post processing routines, the following is done for the computation of SIFs via this method: N odes near the crack tip (I have used the quarter point nodes of the first element adjacent to the crack tip) are chosen and cataloged with their neighbors on the opposing faces. The radius for each of these quarter point elements is measured to the crack t ip. A coordinate system is created that has components tangent to the tip line, locally normal to the face, and in the direction of the crack plane. These are created by using Macros for establishing the relative positions of nearby nodes and cataloging th eir normalized directions for later scripting. The distances of the opposing nodes measured in the new coordinate system are substituted into the above equations. The resultant SIFs are output to a text file for later manipulation and are saved with the same file name as the job which the data has originated from. For ease of use and flexibility, the displacement correlation routine was encoded via a Python code subroutine (using the ABAQUS Python library) that would run on job completion and output all thr ee SIFs to an external text file for array manipulation. 3.3 Parametric Studies Even with the SIF calculation code written, there is still much work to be done in terms of the number of models to be run and how the resulting data will be organized for any parametric study. For example, a given model with a certain geometry may need to have 10 load orientations analyzed. This means that 10 jobs should be created for one .inp files which refer to different load files. After each analysis is completed, the post processing code must save that resulting data so it can be easily found. Then this data must be read into software that can convert the units, manipulate the arrays, and create useful plots. This entire process should be automated and tell the user when it has PAGE 74 74 completed to remove user downtime. This has been done by using MATLAB in the structure below in Fig 36 3.4 Physical Effects The contact patch between a ball and raceway can vary in ellipticity based on the radii of curvature of the ball/raceway p air. There are two contact patches on opposite sides of each ball but that of the inner raceway is smaller and therefore has a higher max pressure in its distribution. See fig 37 for a visual description. If there exists a flaw on the ball surface and it is on the path where a contact patch will pass, th e n a crack will alternately see contact patches approach it and move away. See fig 38 Once the models and the tools are set up, analysis may commence to see what physical effects drive up the SIFs. For example, there is the magnitude of the normal load, the shape of the load (i.e. the amount of ellipticity), the orientation of the load (on which side of the crack its on and how far away from the crack it is), the amount of traction that the load exhibits on the ball surface, the effect of friction between the two crack faces, if there is a pressure acting upon the crack faces, the crack depth, the crack width, and the amount of curvature on the crack face. For example, if the depth of crack is considered f or analysis for geometry in 3 7 the below variation for KI is obtained. While most of these effects have been analyzed, the one that seems to exhibit the largest effect upon the SIFs is the amount of traction that exists between the ball a nd the raceway. T he location of the load is even more important when traction exists at the contact interface. Even if the location of the load is limited to immediately next to the crack on the surface, there are still two possible directions of traction, resulting in four PAGE 75 75 possible load orientations. In reality, the most effective orientation for crack opening is illustrated below. See fig 310. For a lubricated ball bearing contact the EHD viscous friction coefficient is in the range of 0.05 to 0.09. The effect of this surface traction on the SIFs is very significant as demonstrated in a plot below (see fig 3 1 1 ) which displays an increase of more than 250% as the coefficient of friction increases from 0 to 0.07. 3.5 Conclusions This chapter has discussed some of the main concerns with modeling a small nonplanar surface flaw with a nonplanar crack tip. The benchmarks achieved include: Meshed models of nonplanar surface crack geometries were created (see Appendix B for a step by step illustration). Stress intensity fact ors are calculated in a geometry that can undergo compression with nonplanar crack face contact. Parametric studies were conducted in an organized fashion to analyze multiple effects. Trends were observed and are discussed in depth in chapter 5. PAGE 76 76 3.6 Fi gures Figure 31. A partial cone crack geometry generated from some experimental cross sections of the crack created with a spline technique. Figure 32. General J integral equation reference coordinate system about the local crack tip. Ti is the tr action term which is equivalent to the forces in the normal direction along the integral path which would be applied to the interior region if it was independent of its surrounding substrate in order to maintain the identical state of equilibrium. PAGE 77 77 Figu re 33. Crack tip coordinate systems for CTOD (crack tip opening displacement). PAGE 78 78 Figure 34. A c crack specimen with subsurface crack face highlighted above. In tension below, indicating opening of the crack faces. Displacements exaggerated. PAGE 79 79 Fi gure 35 Elements from the crack tip in fig 33 reoriented to show that different parts of the crack tip undergo all three of the different types of deformation. (Displacements exaggerated.) PAGE 80 80 Figure 36 Visual illustration of the M atlab code which is used to create, submit, post process, and organize the desired FEA results. Figure 37 An illustration of a ball in a ball bearing raceway system to indicat e the location and orientation of the contact patches. PAGE 81 81 Figure 38 Relative orientation of an elliptical load near a surface crack indicating, cd, a variable defining the change in distance between the crack and load in a single pass. cd PAGE 82 82 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 2 3 4 5 6 7 8 9 10 x 104 KI for =60 near r=1a load for varying aspect ratiosKI* a/b=1 a/b=0.8 a/b=0.6 a/b=0.4 a/b=0.2 Figure 39 KI for =60 near r=1b load for varying aspect ratio showing a semi elliptical crack steeply angled to the surface will produce the SIFs in MPa po=1 MPa and for the geometry in fig 37 along the crack tip. a) b) Figure 310 A) the location and orientati on of load, C and traction, T to have the most amount of crack opening but if the traction direction is reversed as in B) the crack tends to close. PAGE 83 83 Figure 31 1 KI as a function of the coefficient of friction on the surface for po=540 ksi (3.7 GPa), Crack angle = 30 b =250 mm, a =75 mm, a/b =0.3. PAGE 84 84 CHAPTER 4 EXPERIMENTAL ANALYSIS This chapter focuses on the bodies of experimental work that has been analyzed by us or ha s contributed to this work, in the pursuit of calculating a fatigue threshold, Kth. A mixed mode fatigue threshold has not been agreed upon in literature (see Chapter 1) and is necessary for computing what size flaws are allowed into service. The three tests discussed here are: 1 A novel test for fracture toughness (Piotrowski and OBrien, 2 006) for a fracture toughness parameter, Kc. 2 Brazilian disc (BD) testing for finding a mixedmode fracture toughness, Kc. 3 A Vring test for finding a mixed mode fatigue threshold parameter, Kth. 4.1 A Novel Fracture Toughness Test An experimental analysi s conducted by the Aerospace Corporation has contributed to figuring out how to treat these mixed mode SIFs for a single fracture parameter. In this new test, a Vickers indented ball is placed in the compression of opposing conformal platens (see Fig 41). This compression results in bulging at the equator and a tensile field along the surface where the cracks emanate from the Vickers indentation (Fig. 42). With the magnitude of the tensile field determined from FEA and the flaw size known from measurements, SIFs can be calculated by combining the empirical equations of Newman and Raju (1981) with equations to account for the residual stress field induced by the indent, as below: 2 11(/)10.1(1sin)appKaQMf (4 1) 11.130.09(/) M ac (4 2 ) PAGE 85 85 1 222 4(/)cossin fac (4 3) 1.6511.464(/) Qac (4 4) And to account for the residual stress field: 123/2 10.016(/)/resid cKEHPc (4 5) (which ignores the correction from Smith and Scattergood (1992) for the crack shape dime nsions). The SIFs are additively combined and the residual stress correction makes for a reasonable percentage of the final calculation, like one test would indicate: KIC = 3.19 + 1.19 = 4.38 MPa (4 6) Angled indentations are currently the subject of analysis to determine how the mixed modes of deformation contribute to the resulting fracture parameter. The plots resulting from this work are in fig 43 and 44 which indicate a lack of dependence on indentation load and crack size. The test is featur ed in Piotrowski and OBrien (2002) and concludes that the mode I fracture toughness (for NBD 200 material) is 4.850.36 MPa a fracture toughness of 5.481.34 MPa 4.2 V Ring Test 4.2.1 Experimental Procedure to Determine Growth Regime The singleball, v ring test consists of tw o opposing inner rings that squeeze a sandwiched ball as it rotates (see Fig. 45 ). The lower ring is rotated at a certain speed to drive rotation while the upper ring induces vertical load. The two contact paths on the sides of the ball are perfect locati ons for c crack initiation and testing. Balls are cracked with a rigid pendulum device (see Fig. 4 6) PAGE 86 86 After placing the ball into the device the crack is subjected to RCF. After a sufficient number of cycles, the ball is presumed to survive and testing st ops. If there is spallation before this number of cycles is achieved, an accelerometer is tripped and the test stops. The v ring test has been mentioned in prior work for similar testing (Miner et al., 1996). In order to determine what flaw size is so sma ll that it will not grow under the expected life of the part, the following steps must be followed: Run flaws of known sizes and geometries under well understood RCF conditions. Analyze specimens to determine which flaws grew to failure. Simulate each of t hese experiments in FEA to determine each samples maximum SIFs. Combine the SIFs into one parameter, Keff, that will bracket the possible values of Kth. 4.2.2 Specimen Analysis Once testing is stopped images are taken by whitelight microscope images and destructive grinding to determine depth and geometry. Once all tests are completed, the data is compiled in a plot of applied stress versus crack size on the surface. [This data was not available for publication.] (See fig. 71 for an example image.) 4.2 .3 Stress Intensity Factor Determination Of Individual Test Specimens The contact patch size was calculated for the prescribed load configuration using standard Hertzian contact theory. This contact patch was simulated as a pressure distribution in the ways mentioned in the simulation chapter of this document (chapter 3) The SIFs for each case of crack size and load magnitude were stored for later use. [This data was also unavailable for publication.] PAGE 87 87 4.2.4 Calculating A Keff In simple fracture tests, the mode I stress intensity factor may be the only nonzero SI F. In this case it is straightforward to use the fracture toughness or the Kth, for either static or fatigue load experiments, to determine what loading conditions will lead to failure. However, in these rolling contact fatigue tests all three modes are present. How these three modes of fracture should be combined to determine a fatigue threshold has been the subject of study for many and has resulted in a wide variety of answe rs which have no common ground (Richard et al., 2009). For Kth we have continued to use the critical strain energy release rate (CSERR) discussed below. 4.3 Brazilian Disc Test Herein, the BD test is implemented to determine mixed mode fracture toughness of silicon nitride. S ilicon nitride has material properties that are favorable to hybrid bearing applications including a low coefficient of thermal expansion, corrosion resistance, and a third the density of common bearing metals (Jahanmir, 1994). T he material also has a low fracture toughness of about 6.0 MPa OBrien, 2006). In the course of the manufacturing process, silicon nitride balls are subjected to a lapping process, to achieve the final submicron accuracy in ball diameter and circularity, where the balls are susceptible to low velocity ball collisions. These low velocity interactions induce small surface cracks on the surface of the ball that, under operating conditions in a ball bearing, are subjected to severe rolling contact fatigue (RCF) conditions (Levesque and Arakere, 2008). Under mixed mode loading crack growth can lead to spallation and, as a result, these surface cracks are the biggest contributors to ball (and therefore bearing) failure (Levesque and Arakere, 2008, Wang, 2000). If finite element analysis (FEA) can simulate these flaws under service PAGE 88 88 conditions, then the nature of the effective crack driving force resulting from stress intensity factors (SIFs) KI, KII and KIII must be properly understood in order to determine whether or not a crack is in the regime of growth. Herein, the focus is on defining an effective mixed mode fracture parameter for silicon nitride. 4.3.1 Experimental Procedure Brazilian disc NBD 300 silicon nitride specimens are 25.4mm (1 in.) in diameter and 2 mm thic k. They are cut from bar stock discs are precracked using about 20 collinear Vickers indentations performed on a Zwick Vickers Hardness Testing Machine using 20 kg load such that the penny cracks formed by individual indents link up to for m one long penny crack. Multiple indents are used because the load required to propagate a penny crack from a single indent is very high, causing the discs to fracture at the support platens rather than propagate the center crack. The dimensions of the cra cks induced are recorded of table 1. With the specimens prepared in this fashion, they are placed into a conformal fixture in the plane of loading with the cracks oriented vertically or slightly off the load line as illustrated in fig. 49. The specimens are then loaded at a rate of 100 lbs/min (faster loading can result in fracture occurring at the platens instead of growing the created center crack) until fracture and then the test is stopped. The load at which the specimen cracked and dimensions of the crack are recorded for each specimen at a number of different angles to the load line (see Table 42). This information will be used to calculate SIFs for the specimens at fracture in FEA. PAGE 89 89 4.3.2 Analytical Expressions for MixedMode Stress Intensity Fact ors Prior research has used the equations for stress intensity factors (SIFs) for small (with respect to the specimen thickness) semicircular flaws (formed by Vickers indentation) loaded under mixedmode conditions approximated by (Marshall, 1994): 1 2 2 1sinIKYc (4 7) 1 2 2sincosIIKYc (4 8) where is the applied stress, cab is the equivalent semicircular crack length, a the length of the semi minor axis of the critical flaw and 2b the length of the semi major axis of the elliptical crack. Y1 and Y2 are geometric constants which are equal to 1.65 and 1.55 respectively for indented specimens (Petrovic and Mendiratta, 1976). These equations were implemented to approximate the SIFs for semi elliptical flaws from equations for angled throughcracks in a biaxial stress field with a geometric constant. E ven though experimental observations indicate crack closure at angles around 35 ( as there is observed specimen shattering at angles higher than this ), these equations do not discern crack closure and will always give non zero KI valu es ( excluding =0) Also, the equations (though attributed to a semi elliptical crack) only gives SIF values for a single point along the crack front instead of giving SIFs along the entire crack front as could be expected and also neglect KIII. For b = 0. 02m and a = 0.0025m, KI and KII are plotted as a function of the angle that the specimen in loaded at in figure 410. 4.3.3 Stress Intensity Factor Calculation Via Finite Element Modeling N umerical methods are resorted to f or calculating stress intensity f actors for this centrally cracked, three dimensional (3D) BD. The finite element method (FEM) is PAGE 90 90 utilized to model each of the cracked specimens under the load which induced cracking to calculate the critical SIFs. Firstly, the disk is modeled uncracked and subjected to a uniform pressure distribution over the small contact region on the top and bottom edges of the disk similar to as observed in experiment. The stresses at the center of the disk are compared to those calculated analytically before proceedin g. The stresses at the center of the disc can be calculated according to: 2 (0,0)xP BD (4 9) 6 (0,0)yP BD (4 10) (0,0)0xy (4 11) where y is in the direction of load, P is the total l oad, B is the specimen thickness, and D is the diameter of the disc (Mitchell, 1961) and the full stress state for an uncracked disc is given by: 2244 122()()1 2xPRyxRyx BrrR (4 12) 33 44 122()()1 2yPRyRy BrrR (4 13) 22 44 122()()xyPRyxRyx Brr (4 14) 222 1() rRyx (4 15) 222 2() rRyx (4 16) PAGE 91 91 where variables are defined in figure 49. The 3D FE model stress results on the mid plane of the disc showed excellent agreement with the analytical expressions (Eqs. 412 4 13). The FE mesh for a 3D disc model with a center crack with varying angles to the load necessitates the use of submodeling as the disc is ten times wider than the crack but does not require the same mesh density Also, f or simpli fying t he remeshing process for cracks of different angles and aspect ratios a smaller square submodel was used (See fig. 4 9). The boundaries of this model are displaced according to a submodeling technique, where the displacements are applied to the boundary o f a small cracked block to simulate its being a part of a much larger half space (Dessault Systmes, 2008). B enchmark studies were conducted t o ensure that the boundary of the submodel is sufficiently far away such that the crack presence does not affect i ts displacements. The SIFs calculated for an analysis of mixed mode Kc parameters was a two part process. Firstly, the meshed cracks are analyzed with FRANC3D/NG by M integral decomposition where material is still treated as isotropic (Banks Sills et al., 2005). The SIFs are computed as a function of position along the crack front, and not just at one location as in Eqs. ( 4 7 4 8 ). The SIFs calculated from the BD FEM simulations neglect the presence of indents and the resultant residual stress induced in the specimen. To correct for this the effect of the residual stress is superposed on the mode I SIF using the well known equation (Lawn, 1993, Anstis, 1981): 1/2 3/20.016resid Ic dEP K Hc (4 17) PAGE 92 92 where E is the Youngs modulus, H is the hardness (and equal to 15.5 GPa according to indentation testing), P is the indentation load (196 kN), and cv is the distance from the center of the last indent to the end of the large semi elliptical crack (here 0.47 mm is used as measured in sample S3 for =10). So KI is the result of superposition of mode I SIF contributions due to the str ess field from the applied load computed via FEA and from the residual stress field: FEAresid IIIKKK (4 18) 4.3.4 Mixed M ode Fracture Criteria Once these SIFs are calculated, there are multiple equations that could be used to calculate an effective single mixed mode fracture parameter. For this analysis, equations that have a basis in physical behavioral characteristics are relied upon rather than the multitude of empirical equations that have been proposed in other works (Qian and Fatemi, 1996). The noncoplanar strain energy release rate (NCSER) criterion assumes that the crack growth starts at a critical value of the strain energy release rate and it grows in a direction along which it is maximal. It is defined by (Hussain et al., 1974): 2 1/2 2 1112 221 1II II NCESRI IIKK KKbbb KK (4 19) 2 11 2 2413cos 3cos b (4 20a) PAGE 93 93 12 2 232sincos 3cos b (4 20b) 2 22 2 24(95cos) 3cos b (4 20c) Griffith pr oposed a maximum normal stress (MNS) criterion where fracture nucleates when the maximum stressed point on the surface reaches a specific tensile stress and the direction of the crack growth is normal to it (Griffith, 1920). The fracture condition is given by (Jayatilaka, 1979, Chen et al., 1986): At the crack depth, 2 221 12124 2MNS I I IIIKKvKvK (4 21a) On the surface, 23 coscossin 222MNS I IIKKK (4 21b) where is the angle between the initial crack plane and the plane of propagation (shown in fig. 49). Other methods of crack propagation, in mixedmode conditions, have been proposed with energy considerations. Sih developed a theory based on the strain energy distribution near the crack tip and defined a parameter called the strain energy density factor, S, given by (Qian and Fatemi, 1996): 2 22 11 12 22 332I III II IIISaKaKKaKaK (4 22) Sih (1974) proposed that crack propagates in a radial direction along which S is min imum and the crack growth starts at a critical value of S. At the depth, PAGE 94 94 1/2 221 12MSEI IIIKKK v (4 23a) At the surface, where (3)/1 vv and KIII=0 22 11 1214 162(16)sin 1cos1cos3cos1 412 1MNS I III IIK aKKKa K vv (4 23b) where, for Poissons ratio of 0.25, 111 1coscos 16 a (4 24a) 121 sin2cos1 16 a (4 24b) 221 11cos1cos3cos1 16 a (4 24c) 331 4 a (4 24d) Sihs minimum strain energy density (MSED) hypothesis is the most accepted theory for failure in mixedmo de loading for brittle materials (Chen et al., 2002). The fourth fracture criterion used for the calculation of KIC is based on the assumption that the crack starts propagating at a critical value of the strain energy release rate in the same plane as that of the initial crack (Erdogan and Sih, 1963) 1/2 21II GCI IK KK K (4 25) If the strain energy release rate is considered to be the most critical component to drive crack growth than is evaluated according to: PAGE 95 95 22211 2III IIIGKKK E (4 26) (where is the shear modulus) to calculate an critical strain energy release rate (CSERR) SIF parameter as if there was only one mode of crack tip displacement according to: GC cKEG (4 27a) 2 2 2) 1 (III II I cK v K K K (4 27b) The "coplanar method" becomes identical to using a CSERR, G or Gc. The equations used are the same rewritten as a single SIF parameter, [where only two parameters ( KI and KII) are present]. By ignoring the tearing mode of crack tip displacement 0IIIK and then: 221IIIGKK E (4 28) a nd to substitute and rewrite: 1/2 2 2222 21 1II GC III IIII IK KEKKKKK EK (4 29) and this is the same as the CSERR criterion. These equations have also been written by Qian as : At the crack depth, 221 1GI IIIKKK v (4 30a) And at the specimen surface PAGE 96 96 2 22 211 cos1cos3sin3cos1 122G I II IIK KKK v (4 30b) For good measure, a couple of empirical equations are included for comparison. From Tanaka (1996), 1/4 4441 88 1TIII IIIKKKK v (4 31) Also, from Lucht (2009) 2 2 Mi IIIIIIKKBKK (4 32) where B =1 reduces to more familiar equations. 4.3.5 Analysis Procedure The best mixed mode criteria will be determined by analyzing each of the theories by using SIF s that are calculated as a function of position on the crack front. Previous work has dealt with through cracks and only had two points where SIFs were necessary to be calculated and others have used empirical equations that only provide a single SIF for the entire crack front (Banks Sills, 1998) Here FEA provides us the advantage of calculating SIFs along the entire crack front. The calculated SIFs are found in Figs 413 4 15. Having these SIFs as a function of crack front position complicates matters. SIF combinations can be high at either the surface or the depth of the crack. On the surface, SIFs can be affected by the residual stress present from the multiple indents but there is no KIII component (which is the component that has the least amount of analysis in prior work). Only a battery of analyses can be conducted to see th e variation of the mixed mode parameters between the experimental samples We compare each PAGE 97 97 method where the crack meets the surface, the depth of the crack and where the max value of the parameters is calculated on the crack front FEA requires the specimen dimensions a s an input, load at fracture, and crack dimensions. The specimen dimension is well defined and has negligible error. The load at fracture and the crack shape are measured at the fracture of the specimen in experiment. While the load at fracture is measured by the load frame, its error is also quite negligible. The crack dimensions, however, are only determined by analyzing the fracture surface as in fig 4 20 By prescribing an average err or of 50 um to the crack depth, than the resulting error in the SIFs would be: ,,,, ,,150 50 2IIIIIIxyxy IIIIIIm K amErrorinK a (4 33) w hich is significant in the light of the specimen dimensions in table 4 2 An additional amount of variation comes from the material itsel f. Just because specimens with the same size cracks of the same type of material are subjected to the same type of loads does not guarantee failure as a result of the variation that occurs within the material (and between materials of different grades, as well). In previous studies on fracture toughness of a different grade of silicon nitride (TSN 03 NH where only KI was present in testing) had a variation of 4.850.36 MPa OBrien, 2006). 4.3.6 Results A nd Discussion Stress intensity f actors were calculated for each experimental study by simulating the measured crack dimensions under the load that induced failure. The SIFs were tabulated for reference. (See Table 41) The SIFs tabulated are those that occur where PAGE 98 98 the crack meets the surfa ce of the specimen. T hese are plotted since this is where the angle of crack growth, was measured on the experimental specimens. Also, at this location KIII=0 which makes these Kc calculations easier and transferrable to future research ers In addition, experimental knowledge indicates that the point where crack growth is likely to init iate is at this location Due to the geometry of the experiment, as the angle of the crack increases, KI decreases while KII increases. KIC is clearly measured for =0 and is very close to the 4.6 MPa KI loadin g seems to induce the smallest of the calculated Kc parameters but with KII contribution the mixed mode Kc ends up at an average 5.891.13 MPa deviation among tests). Figure 419 shows the fractured specimens as a function of the angle The fracture of the specimen at =90 sted entirely of pure mode I loa d ing, occurred along the initial plane of the fl aw. However, at all other angles the crack deviated from its initial plane. At the same time, it can be noticed that the inclined crack deviates from its initial plane to propagate in a direction normal to the maximum principal stress. It was also seen that in some of the samples, the cracks did not encompass all the Vickers indents on the sample. The di sk shattered to several pieces at =55 =50 4 19(D)). Note that these angles correspond to 35 and 40 respectively, from the loading axis. FEA demonstrated that there this observation was due to crack closure around 31 without accounting for the residual opening stre ss produced by the indent Examples of cracks on the fracture surfaces are shown in Figure 43. PAGE 99 99 4.4 Conclusions 1 A novel test to calculate Kc from Michael J. OBrien (Piotrowski and OBrien, 2006) was used to calculate fracture toughness from the cracks propagating from an indent as 4.850.36 MPa TSN 03 NH. 2 Vring testing was used to experimentally calculate a CFS. 3 Vring tests were simulated to calculate Kth boundary. 4 Measured KC and KIC through BD test simulations. 5 Investigated applicability of mixed mode parameters through FEA simulation of BD testing. 6 Characterized fracture toughness for silicon nitride as a function of modemixity based on multiple criteria. 7 The max strain energy release rate criterion seems to be most effective, as its standard deviation is lowest for all the visited methods and it predicts 5.891.13 MPa PAGE 100 100 4.5 Figures Figure 41. On the left, the designed fixture an d a tested ball are displayed. On the right, the platens and ball are placed into the fixture in the orientation for testing. PAGE 101 101 Figure 42. Close up of indentation flaw, with measurement taken in the direction of growth due to applied load. Arrows indi cate direction of tensile field. Figure 43. The fracture toughness of tested specimens versus the final crack length they became, indicating a lack of interdependence. PAGE 102 102 Figure 44. The fracture toughness of tested specimens versus the Vickers inden tation load indicating a lack of interdependence. Figure 45. FEA of the ball in the platens undergoing compression, clearly indicating the tensile field along the equator of the ball. PAGE 103 103 Figure 46. Pendulum device for impacting balls and inducing cr acks PAGE 104 104 Figure 47. Images of V ring single ball test rig. Lower image shows contact path and cracks oriented within them. PAGE 105 105 Table 4 1. V ring test conditions for 28.575 mm (1 in.) balls. (Miner et al., 1996) Test Property Measurement Hertzian Stress 3800 MPa (550 ksi) Shaft Speed 7800 rpm Stress Rate 6*10 6 cycles/hour Lubricant Flow Rate 3.4 kg/min (7.5 lb.min) Lubricant Temperature 366K (200F) Raceway Surface Finish 0.05 Figure 48. Diagram illustrating measurement parameters to classify a c crack as observed on the surface of a cracked specimen. PAGE 106 106 Table 4 2. The sample number, turning angle, crack dimensions and fracture load of each specimen examined. 90 Sample (deg) 2b (mm) a (mm) P (N) 0 S1 0 3.01 0 .26 8940 0 S2 0 3.10 0.39 7659 10 S1 12 5.00 0.23 7200 10 S2 18 2.68 0.49 10640 10 S3 15 1.94 0.49 9341 20 S1 26 0.65 0.28 10100 20 S2 28 0.63 0.23 12342 20 S3 20 0.50 0.19 13950 30 S1 19 1.38 0.32 9870 30 S2 21 2.70 0.44 9010 30 S3 38 0.49 0.23 14540 Table 4 3. The mean Kc values in MPa the standard variation std) Location KNCESR KMNS KMSE KGC KTAN KL Max 6.191.16 7.041.59 6.721.47 5.891.13 7.552.29 7.571.28 Surface ( K III =0) 5.181.49 2.253.08 6.611.57 4.961.79 6.672.09 7.041.3 2 Depth ( K II =0) 3.251.54 6.201.35 6.671.63 5.721.13 7.142.15 6.891.49 PAGE 107 107 Figure 49. Schematic showing the flaw orientation on a Brazilian disk specimen as well as the coordinates and components used in the analytical stress equations. PAGE 108 108 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 (in radians)K/ KI KII Figure 410. Nondimensionalized KI and KII (with resptect to normal and shearing stresses respectively) from empirical equations for different angles of loading orientation for b =.02m a =.0025m. PAGE 109 109 Figure 411. BD test model showing a) global uncracked model wit h square partition indicating submodel boundaries b) cracked meshed submodel for sample 2 for =20 Figure 412. The displacement field of Sample 1 for =20, with a scale factor of 100 and translucency applied to view the crack superposed with elements of the circular contour regions to illustrate the presence of KII at the face and KIII at the depth. PAGE 110 110 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 1.5 2 2.5 3 3.5 4 0 to 1 along the crack frontKI (MPam) =0 S1 =0 S2 =10 S1 =10 S2 =10 S3 =20 S1 =20 S2 =20 S3 =30 S1 =30 S2 =30 S3 Figure 413. KI for all analyzed specimens as a function of position along the crack front. PAGE 111 111 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 8 6 4 2 0 2 4 6 8 0 to 1 along the crack frontKII (MPam) =0 S1 =0 S2 =10 S1 =10 S2 =10 S3 =20 S1 =20 S2 =20 S3 =30 S1 =30 S2 =30 S3 Figure 414. KII for all analyzed specimens as a function of position along the crack front. PAGE 112 112 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 1 2 3 4 5 6 7 0 to 1 along the crack frontKIII (MPam) =0 S1 =0 S2 =10 S1 =10 S2 =10 S3 =20 S1 =20 S2 =20 S3 =30 S1 =30 S2 =30 S3 Figure 415. KIII for all analyzed specimens as a function of position along the crack front. PAGE 113 113 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 2 4 6 8 10 12 0 to 1 along the crack frontKTanaka (MPam) =0 S1 =0 S2 =10 S1 =10 S2 =10 S3 =20 S1 =20 S2 =20 S3 =30 S1 =30 S2 =30 S3 Figure 416. Keff using Tanakas relationship of three SIFs as a function of crack front position. PAGE 114 114 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1 2 3 4 5 6 7 0 to 1 along the crack frontKeff (MPam) =0 S1 =0 S2 =10 S1 =10 S2 =10 S3 =20 S1 =20 S2 =20 S3 =30 S1 =30 S2 =30 S3 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 1 2 3 4 5 6 7 0 to 1 along the crack frontKeff (MPam) =0 S1 =0 S2 =10 S1 =10 S2 =10 S3 =20 S1 =20 S2 =20 S3 =30 S1 =30 S2 =30 S3 Fi gure 417. Keff using the strain energy density release rate using all three SIFs as a function of crack front position where a) includes KIII and b) neglects it. PAGE 115 115 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1 2 3 4 5 6 7 8 9 0 to 1 along the crack frontKeff from Lucht (MPam) =0 S1 =0 S2 =10 S1 =10 S2 =10 S3 =20 S1 =20 S2 =20 S3 =30 S1 =30 S2 =30 S3 Figure 418. Plot of 2 2 Mi IIIIIIKKBKK where B =1. Figure 419. Fractured spe cimens at various values of PAGE 116 116 Figure 419. Continued Figure 420. Fracture surface of BD specimen showing the crack front growing as it links up with newly formed cracks from indents along the surface. PAGE 117 117 CHAPTER 5 CFS EVALUATION: COMPUTATIONAL INVESTIGAT ION 5.1 I ntroduction Silicon nitride balls, used in hybrid ball bearings, are susceptible to failure from fatigue spalls emanating from pre existing partial cone cracks that can grow under rolling contact fatigue (RCF). The range of 3D nonplanar surface flaw geometries subject t o RCF are simulated to calculate mixed mode stress intensity factors to determine the critical flaw size (CFS), or the largest allowable flaw that does not grow under service conditions. The cost of nondestructive evaluation (NDE) methods for silicon nitr ide balls scales exponentially with decreasing CFS and increasing ball diameter and can become a significant fraction of the overall manufacturing cost Stress intensity factor variability is analyzed for variations of the location and orientation of the l oad relative to the crack, the geometry of load, and full slip traction. The modeling techniques utilized in the creation of a 3D FEA model is discussed and the maximum tensile contact periphery stress is examined for effect on crack driving force under RC F. The CFS results are presented as a function of Hertzian contact stress, traction magnitude, and crack size. The evaluation of a critical flaw size (CFS) is of immediate engineering relevance to the hybrid ball bearing industry towards defining refinemen t necessary for non destructive evaluation methods for silicon nitride ball quality control. Also, a CFS can be applied for developing a fracture mechanics based life prediction methodology for hybrid bearings. The performance of hybrid silicon nitride (Si3N4) ball/steel raceway bearings has been shown to be superior to their all steel predecessors (Miner et al. Tanimoto et al. ). However, silicon nitride balls are sensitive to surface defects and can fail from PAGE 118 118 fatigue spalls emanating from pre existing c c racks or partial cone cracks (see Fig. 51), due to crack growth driven by RCF (Hadfield et al. 1993a, Levesque and Arakere, 2008). Silicon nitride exhibits favorable material properties for application in hybrid high speed ball bearings such as high compr essive strength, high hardness, a third of the density of steel, low coefficient of thermal expansion, and high corrosion and temperature resistance (Jahanmir, 1994). Unfortunately, they also have low fracture toughness (4 6 MPa m ) (Piotrowski and OBrien, 2006). Silicon nitride ball manufacturing involves a lapping process to achieve the submicron ball finishing dimension and surface characteristics that produces (unavoidable) low velocity ball collisions resulting in surface cra cks nucleated by the radial tensile stress field which is maximal at the contact periphery (Levesque and Arakere, 2008, Wang and Hadfield, 2000). Partial cone or c cracks are the result of the oblique impact of brittle spheres (Frank and Lawn, 1967) and th ese cracks are considered the most damaging surface defect that limits ball life in hybrid bearings under service conditions (Evans, 1993, Hadfield et al., 1993b). C cracks are more commonly observed in sphereto sphere collisions and they not only have non planar crack faces but also possess nonplanar crack tips, as opposed to the axisymmetric cone cracks that result from the normal collisions of spheres. Figure 51 illustrates the 3D nonplanar crack features of a c crack. This increased complexity makes their shape more difficult to describe and, therefore, more complicated to analyze in any linear elastic fracture mechanics (LEFM) based analysis. While the flaws (illustrated in Fig. 51) can be readily induced and are generally present, methods for th eir detection prior to service is still in development. Development PAGE 119 119 of a non destructive evaluation (NDE) method is complicated because these cracks are surface and are diff icult to find under a microscope. In addition, the material is nonconductive and only slightly translucent which complicates inspection procedures. From a fracture mechanics and structural integrity perspective, the larg est allowable surface flaw that does not propagate under RCF loading is of design significance, and is termed the critical flaw size. The goal of this research is to present a systematic procedure to compute the CFS based on fracture mechanics principals, RCF loading, and ball material properties. This analysis is driven by the cost of NDE method for silicon nitride balls, which scales up very steeply with decreasing CFS and increasing ball diameter. Thus the cost associated with NDE can become a significant fraction of the overall manufacturing cost of the silicon nitride ball. Prior work has presented a qualitative and quantitative description of the range of possible partial cone crack shapes in three dimensions depending on the initial conditions of the interaction (coefficient of frict ion and contact patch size) during oblique spherical contact (Levesque and Arakere, 2008). With this range of possible nonplanar shapes characterized, they may be analyzed under RCF conditions that are seen in service, to determine which physical effects h ave the largest impact on crack tip displacement. Rolling elements are manufactured with silicon nitride because of the severe RCF conditions the material can withstand. These conditions are severe for any pre existing surface flaws that exist. For example in the case of a main shaft jet engine bearing inner ring rotating at 15,000 rpm, and as a worst case scenario considering that the c crack lies along the ball track (see Fig. 52), the crack can experience in excess of PAGE 120 120 90,000 fatigue cycles per minute ( 1,500 Hz), or 5.4 million cycles/hour. Furthermore, as the c crack passes through the contact patch, the crack tip experiences a mixedmode loading at every point along the crack front. T he effects of the location and orientation of the load relative to the crack, the geometry of load, the geometry of the crack and the respective effects on crack tip stress intensity factors (SIFs) and, therefore, are presented There is no published work concerning an elliptical Hertzian contact load inducing RCF for nonplanar c cracks, which is what is most often seen in service applications. Herein, the discussed effects on the SIFs are analyzed in the ball bearing system and the process of extracting them. 5.2 Analytical Procedure 5.2.1 Crack G eometries I nduced B y O blique I nteraction s Of Silicon Nitride Balls The c cracks produced during oblique interaction have multiple geometric features that need to be understood before modeling is undertaken. It can be shown that the crack size scales with the ball radius, R, and that the velocity needed to induce cone cracks is very low and can be shown to be proportional to R3 (see Appendix I). Also, it is noteworthy that the angular extend of the c cracks have been observed to be roughly 90 to 120 (Levesque and Arakere, 2008, Wang, 2000). However, cracks that are the result of normal indentation of spheres have a weak correlation about the depth to which they extend into the material and exhibit considerable variation even fo r similar indentation conditions (Johnson, 1987). Avoiding generalities, the range of the geometry of cracks produced by ball onball collisions was established in a paper using a stress state and numerical iterative growth and com pared with experimental images (Levesque and Arakere, 2008). The angles to which the crack extends subsurface have already been bracketed within an established range in the same reference. The resultant shape PAGE 121 121 from these oblique interactions is complex and other prior works charac terized the approximately modeled crack geometry with sets of parametric equations (Wang, 2000). However, the complexity of the shape leads to a difficulty in generalizing results and in finding which geometry of crack will be the most severely affected by RCF. With this information in mind, the c crack geometry can be modeled as a partial cone, or a partial frustum, whose dimensions are limited by inequalities which are equations of an offset cylinder. In other words, the shape is an axially revolved qua rter ellipse with limitations to its extension in the rand direction (considering cylindrical coordinates, see Fig. 53). 2 *222 **1 zaxyhk ab (5 1) Where, 111 tan 22range rangey x (5 2) *22 1 *111 sectan 2sin() 22rangey haxy ha xha (5 3) where a* and b* are the axes of the ellipse on the x and z axis, respectively, h and k are the x and z coordinates of the center of the quarter ellipse, respectively and is used to denote the fraction of a circle which the crack is seen to be on the surface of the specimen ( see Fig. 53). PAGE 122 122 5.2.2 F inite Element Modeling Of A ThreeDimensional Surface Flaw The shape of the crack is critical to understanding the difficulties of generating an adequate mesh for finite element analysis (FEA). In order to accurately model a crack f or SIFs, quadratic elements are utilized to account for the large gradients of stress that occur throughout the model as a result of the Hertzian contact loading. To accurately model the crack tip stress, collapsed hexahedral elements are used at the tips with the midside nodes moved to the quarter point (Andersson, 1996). Since hexahedral elements are utilized at the tip, C1 continuity of the shape functions is maintained, which is desirable especially if SIFs are intended to be extracted in the region jus t surrounding the crack tip. It is noteworthy to mention that others have developed methods to extract SIFs from tetrahedral based meshes, which simplifies meshing and remeshing in the instance of crack growth but there is still much variability in the met hod (Rajaram et al., 2000). However, it would be difficult to create a mesh that would be made of all hexahedral elements and avoid distorted warning elements considering that the geometry is curvilinear and closed. Through benchmark analyses types of mesh es and their SIFs were tested against previously established solutions (like that of Newman and Raju, 1981). As a result, a tube of elements can be created which were hexahedral, about the crack tip region, but were tied to a bordering tetrahedral region t hat meshed the majority of the cracked body and maintain SIF accuracy. SIF extraction from this mesh is feasible through a few methods (including stress matching, and virtual crack extension) but the most appropriate approach for the analyses is via crack tip opening displacement (CTOD) correlation since it accounts for the effects of all physical phenomena without a need to change its formulation (as in PAGE 123 123 contour integral, energy based methods). The equations for displacement correlation can be written as: 12 22 4(1)IE K uu r (5 4a) 12 22 4(1)IIE K vv r (5 4b) 12 22 4(1)IIIE K ww r (5 4c) where E and are the Youngs modulus and Poissons ratio, u v and w are the three displacements in the directions normal to the crack face, in the dire ction of the continuing crack plane normal to the tip, and tangent to the crack tip respectively, and r is the distance from the crack tip, where these displacements are measured. Figure 4a illustrates the CTOD coordinate system. Figure 4b shows the 3D dis placement vectors at the tip of the nonplanar c crack mesh. The FEA model represents a computationally intense effort due to the dense mesh (for capturing the Hertz stress gradients) and implementation of quadratic elements (for moving the midside node of the tip elements to the quarter point, preserving the 1/ singularity at the crack tip), and accounting for contact accuracy. Mesh refinement for Hertz contact problems has been visited in other works (Gu, et al., 2005) usually for stress calculations. H ere the displacements are of primary importance, especially those in the crack tip vicinity. A mesh refinement study has been performed in the tensile region on the contact periphery and established guidelines for accurate determination of SIFs. The mesh density utilized is about 50 elements along the crack to the tip where the hexahedral region has 47 rings of elements in the contour region, PAGE 124 124 meaning that for this length of the longest grains in TSN 03 NH silicon nitride). Contact of the crack faces is defined by using a surface based contact definition with no friction where a direct, iterative sol ver is used to achieve accurate displacement and stress field solutions. The direct iterative solver must be used as the accuracy in displacement results required for adequate SIF calculations is high. The direct iterative solver is the most accurate sol ver available in ABAQUS (Dessault Systmes, 2007). The mesh refinement previously established is based upon (a) refinement that the Hertz contact requires for adequate resolution and (b) refinement in the vicinity of the crack tip necessary for accurate di splacements needed for SIF computation. Loading is done via FORTRAN user subroutines DLOAD and UTRACLOAD for the normal and traction loads. The equations to describe the pressure distributions can be written as: 22 22(,)dd oxxyy pxyp ab (5 5) and (,)(,) fxypxy (5 6) where xd and yd are the distances from the global coordinate system to the load center, a and b are the dimensions of the ellipse along these dimensions, and is the friction coefficient for the moving load in a full slip interaction (see Fig. 5 5). 5.3 Mixed M ode S tress I ntensity F actor s D ue to R olling C ontact F atigue For generalization, the trends observed in the SIFs for each of the possible physical param eters should be investigated and quantified. For example, if the max PAGE 125 125 Hertzian pressure is increased, the SIFs are mostly observed to directly scale with the value of Po. The maximum continuous operating Hertzian contact pressure, Po, typically encountered in turbine engine mainshaft hybrid ball bearings is in the range of 1.7 2.8 GPa ( 250412 Kpsi), for achieving reasonable L10 fatigue life for the bearing raceway material. SIF results are presented at the upper end of that range, at 2.8 GPa, to obtain cons ervative estimates on CFS for ball material. For the orientation in Fig. 56 without traction, the linear effect of max pressure on SIFs is demonstrated in Fig. 57. This information may be applied by others in the field who require SIF calculations for si milar phenomena but different magnitudes of pressure in their contact ellipse. In the two dimensional case, the width of the contact, the angle of the crack, and the depth of the crack are enough parameters to define an RCF model (Bogdanski and Trajer, 2 005). However, for 3D analysis, the model must be characterized by the load aspect ratio and specific crack dimension, making generalization of results difficult. An analysis of SIFs as a function of approaching load is illustrated in Fig. 58. There are m ultiple static steps to analyze these c cracks under a RCF cycle (or single passage of the load rolling over the crack). As the load approaches the KI SIF increases until the load is very close to the crack tip, resulting in crack closure. At this point, KI rapidly decreases to zero as the contact patch moves directly over the crack. As the contact patch advances further the crack tip opens due to the influence of the maximum tensile stress on the trailing edge of the load. While KI reaches a second maximum in the trailing periphery edge, it is not as high as that of the approaching edge because the crack angle (illustrated in Fig. 56) is such that the depression from the load has a tendency to cause the crack to close as shown in Fig. 58. As the load mov es away PAGE 126 126 from the crack all the SIFs eventually decay to zero. In figure 58, the relative orientation of crack tip and load allows for influence of the load at 3a but KI and KIII decay quickly thereafter. For design application, an effective stress intensi ty parameter is used that embodies all modes of crack deformation as: 0.5 22 21 2(1)eq III IIIv KEGEKK K EE (5 7) This representation of the three modes of crack tip displacement was chosen due to its physical basis in fracture mechanics (since it utilizes the st rain energy release rate, G ) (Anderson, 2005). The most severe orientation of a purely normal load (i.e. no traction load) is such that the periphery stress is near the crack on the side where the crack is oriented (see Fig. 5 6). However, the distance to the center of the approaching load ( dx ) where max SIFs are calculated is affected by the crack angle. 5.3.1 Effect O f T he V ariation O f L oad E llipticity O n S tress I ntensity F actor s The ellipticity of the contact patch also influenc es the induced SIFs. For a generalized Hertzian elliptical contact seen in a ball bearing the ellipticity ratio ( ba ) can vary depending on the raceway curvatures, speed, radial and thrust loads. A nearly circular contact patch is o ften seen in RCF test setups like the ball on rod test, the four and fiveball tests, v ring test and more. Also, if ba is very small, the contact can be reminiscent of a line contact. Since the crack is moving with the rotating bal l, the relative orientation of crack and load can vary greatly. The effect of ellipticity of load on the SIFs from the semi minor axes is studied in Figures 5 9. Again, SIF results are presented for a peak Hertz stress Po =2.8 GPa (412 kpsi) The circular case yields comparatively PAGE 127 127 higher SIFs. This is due to the nature of the stresses that circular loads can produce at their periphery (Johnson, 1987). 5.3.2 Comparison O f C C racks A nd S emi Elliptical C racks The range of possible crack shapes and sizes has b een established by Levesque and Arakere (2008). The complex nonplanar 3D geometry of the c crack requires parameters such as crack depth and angle, which are difficult to measure, and a function that represents the position of the surface and inequalities to indicate the end of the crack (see Eqs. 51 5 3). Also, some NDE inspection procedures can only detect the width of the crack on the surface. Much of this issue could be resolved by comparing c cracks to a semi elliptical flaw that is fully classified by three parameters (depth, width, and angle of inclination). The semi elliptical crack also yields higher SIFs, for the same peak Hertz stress, Po, leading to a conservative analysis and then the feasibility of simulation becomes much more accessible to ot her workers in the field. With this in mind, comparison of SIFs for the same orientation of elliptical contact load of both a c crack and a semi elliptical surface crack of similar width, depth and angle toward the surface, was completed. See the results i n Fig. 510. When the SIFs are compared it can be notice d that considerable divergence in their values even though they occur in the same relative orientation to the load. However, the maxima and minima of each, though at different locations, are quite s imilar. As a result of the complex geometry and loading, all three modes of crack tip deformation are present. Modeling and experimental efforts need to collaborate to determine the best mixed mode parameter for crack growth of brittle materials under thes e types of conditons for future design purposes. PAGE 128 128 5.3.3 The Effect O f Traction Under elastohydrodynamic (EHD) lubrication conditions in a ball bearing, the effective friction coefficient at the contact is typically in the 0.05 0.09 range. Although the coefficient of friction appears to be low, it will be show n that it has a significant effect on the maximum SIFs reached. The direction of traction also becomes important. When traction is incorporated, it is seen that the worst case load moves immediately n ext to the crack and that the SIFs increase substantially (see Fig. 56). If the traction direction is reversed, the worst case load jumps to the opposite side of the crack where the tensile stress field has the most effect. (It should be noted that this position will not have as high SIFs as that on the opposing side since the normal load only contributes to crack closure on that side of the crack.) During RCF, the elliptical contact patch induces both normal and traction forces onto the ball surface. S ince the bearings of interest are lubricated, the interfacial contact is in fullslip (a hydrodynamically lubricated contact) and the amout of traction on the surface is limited by the shearing properites of the lubricant that is in application. For the mo st part, it is fou nd that =0.05 0.09 is an acceptable range for high speed bearing applications. The traction force is parametrically varied and the max stress intensity factor values along the front are plotted in Fig. 511 below. It is quite apparent that the slope of the incr easing SIFs is quite steep as the traction force is increased. For silicon nitride, the fracture toughness ( Kc= 6.0 MPa reached when =0.2 (even though steel on silicon nitride has a friction coefficient of about =0.15) and indicates that the fricti on coefficient a required variable to evaluate a critical flaw size for ball bearing inspection. The reason for the strong influence of the friciton coefficient is that the surface flaw is easily effected by surface behavior. In PAGE 129 129 addition, there are high no rmal loads that are directly contributing to a traction load according to (where, for example, if there is a total load in the contact of P=2500 lbf than, for =0.07, the total tration load, f is 175 lbf is contributing to crack opening). 5.3.4 Critical F law Size Evaluation R esults If the operating conditions of a bearing are known, and the material properties (especially the Kth for initiating crack growth) are established then a CFS (or flaw size so small that it experiences crack tip loading below the Kth threshold value under the operating conditions) can be found These models can be generalized for exactly this application. For example, it has been established that the SIFs scale linearly with the max pressure, Po, of the elliptical contact patch. It has also been established that a circular load will produce higher SIFs than the commonly found elliptical dimensions. In addition, it has been noticed that similar magnitudes of maximum SIFs across the crack front of c crack and semi elliptical flaws occur under i dentical loading conditions. It has also been establis hed that traction has a significant effect on the SIFs for the RCF of surface flaws. With all these facts in mind, models can be ran that simulate a decreasing crack size near a circular load and plot the change in the SIFs as a function of the changing cr ack size. By nondimensionalizing the SIFs with respect to the max contact pressure ( */IIOcKKPa ), the results can be generalized. The crack dimensions chosen were based on what is commonly observed in experiments. In this case, cracks have b een modeled that have an aspect ratio of 0.3 and are angled 30 from the surface. Figure 12 displays the results. To utilize this graph, all that is needed is the effective threshold value Kth of the material, under mixed mode loading, and the PAGE 130 130 CFS can then be looked up for design purposes after dimensionalizing the data for specific applications. For example, for Po=2.8 GPa, =0.07, and ac=0.5mm, if Kth=2.3 MPa Kth can also be determined via RCF testing of c cracked balls for specific contact loading conditions. For example, if a ball with a crack of a =0.35 mm, and has loading conditions of Po=2.8 GPa, =0.07, ac=.5mm, than Kth=2.4 MPa 5.4 Conclusions The evaluation of a CFS has important engineering relev ance towards defining limits for NDE methods for silicon nitride ball quality control, for reducing the cost of ball inspection, and for developing a fracture mechanics based life prediction methodology for hybrid bearings. To this end, it has been analyzed that complex surface cracks present in silicon nitride balls used in hybrid bearings subject to RCF, using 3D FEA and fracture mechanics principals. While the relative position of load affects the SIFs at any given instant, the max SIFs, used in a calculation, are observed in trends and were utilized in the determination of the smallest allowable CFS. These observations include: Surface cracking can be induced in silicon nitride ball collisions with low velocities of ~16 cm/sec Surface crack dime nsions that result from collisions can be shown to scale linearly with ball radius, R, and the velocity to induce cracking with the cube of the ball radius, R3. C cracks are a complex shape that can be modeled through iterative growth, or solid modeling, w here the shape has been characterized by equations. Maximum Hertzian stress Po, scales the SIFs for ccracks in balls subject to RCF in a nearly linear fashion. Surface traction between the ball and the raceway and its direction has a significant detrimen tal effect on SIFs and crack driving force. PAGE 131 131 The max SIFs are seen while the crack is on the edge of a circular load when compared to approaching elliptical loads. The semi elliptical crack is simpler to model and does exhibit similar magnitudes of SIFs which are necessary for Kth calculations to determine a CFS. CFS can be determined from plot of SIFs versus crack dimension (Fig. 12) by looking up the appropriate mixed mode effective threshold SIF value, Kth, which can be determined via RCF testing of c cracked silicon nitride balls for specific contact loading conditions. PAGE 132 132 5.5 Figures Figure 51. Images of artificially induced c cracks a) ball surface view, and b) a subsurface crosssection extending into the material [C ourtesy of The Timken Company ] PAGE 133 133 Figure 52. The elliptic contact patches on the ball surface and the band which they remain in. Cracks inside and very close to this region (assuming a fixed contact angle and attitude of ball rotation) will experience RCF. PAGE 134 134 Figure 53. A) The coordinate system used in the equations to plot the C crack and B) the rendered 3D crack with notation for the SIF plots to be discussed. PAGE 135 135 Figure 54. Illustration of A) displacement correlation variables around an opening two dimensional crack and B) the coordinate systems are different for each point near the crack tip at a few element borders. Figure 55 Elliptical Hertzian contact load coordinate system and variables for elliptical load simulation. PAGE 136 136 Figure 56. The cross section of the gen eral worst case orientation with traction for crack tip deformation. Here the traction direction acting on the silicon nitride surface is indicated and considerations of relative movement must be based on which surface is driving each contact which is dif ferent for inner and outer raceway interactions. Figure 5 7 KI stress intensity factor for a circular load ( ac crack a variation in the max pressure magnitude. The effect of a scalar increase in pressure is a near scalar increase in KI. PAGE 137 137 Figure 58. A plot of the thr ee SIFs as a function of moving load of Po=2.8 GPa. When abs( dx / a )=1, the load is just on the edge of the crack and the coordinate system is oriented where the crack breaks the surface. The maximum SIFs along the crack front are p lotted and they occur around dx / a = 1.7 for this geometry crack and load. PAGE 138 138 Figure 59. The effect of the ellipticity of the load, on the SIFs on the semi minor axis of the ellipse, for Po=2.8 GPa (412 Kpsi) PAGE 139 139 Figure 510. Sup erimposed SIFs for a c crack and penny crack of similar dimensions under elliptical load ( a b Po= 3.7 GPa (540kPsi)) where the surface crack is tangent to the semi minor edge of the contact patch. PAGE 140 140 0 0.05 0.1 0.15 0.2 0.25 1 2 3 4 5 6 7 KI, KIeq (MPam) KI KIeq Figure 511. The max KI found on the crack front of a semi elliptical flaw that is 30 to the surface and has a` =250 b` =75 as a function of varying friction coefficients on an elliptical contact patch in full slip for Po=3.7 GPa. PAGE 141 141 Figure 512. The nondimensionalized crack semi width plotted against a) the max KIeq and b) the max KI where nondimensional units are on the left and in MPa for po=2.8 GPa and ac=0.5 mm on the right The SIFs are nondimensionalized by dividing the calculated values by poc. PAGE 142 142 CHAPTER 6 EMPIRICAL STRESS INTENSITY FACTOR EQUATIONS FOR SURFACE CRACKS UNDER RCF 6.1 Introduction Reflecting on prior works, there are a few discrepancies. Real RCF produces a complex stress field that is not easily characterized by many of the above methods. In many cases KI is sometimes supplied as a negative number, which, while physically impossible, has ofte n gone unmentioned, i.e. effects of crack closure have been ignored. If KI is supplied, KII or KIII may not be supplied even though all modes of crack tip deformation are present and contributing to mixed mode crack growth. Also, the problem requires a 3D analysis which rules out much of the prior work. Beyond this, previous results were often case specific. The pre existing surface flaws present in silicon nitride balls are typically partial cone cracks or c cracks and their 3D crack geometry was described by Levesque and Arakere (2006). The c cracks have nonplanar crack faces but also possess nonplanar crack tips, making their shape more difficult to describe and, therefore, more difficult to analyze in any linear elastic fracture mechanics (LEFM) based analysis. In this work, comprehensive results for the calculation of SIFs for semi elliptical cracks under RCF are presented. In another chapter a comparison of SIFs for partial cone or c cracks and semi elliptical cracks subject to circular and elliptical RCF will be presented that shows that use of semi elliptical crack with circular contact leads to conservative estimates for critical flaw size, required for defining limits for non destructive evaluation methods for silicon nitride ball quality control. The results generated are of immediate engineering relevance to the hybrid ball bearing industry towards evaluating critical flaw size and for developing a fracture mechanics based life prediction methodology for hybrid bearings. The empirical PAGE 143 143 curve fits will also be of relevance to other areas of component design where contact initiated fatigue damage is important such as gears, roller bearings, and railway wheels. 6.2 Analysis A detailed investigation of SIFs for a semi elliptical flaw under rolling contact fatigue can only be properly done in 3D. This implies a 3D FEA analysis, as any analytical method would quickly become intractable under this complex stress state. We have modeled the penny crack for different aspect ratios, orientations of load, elli pticity of load, and angles to the surface. The resulting trends are plotted and are fitted to equations for ease of use. The FEA model has been created using either FRANC3D/NG, as developed by the Fracture Analysis Consultants or ABAQUS. FRANC3D/NG makes it easier for the user to create meshes of cracked bodies with limited control of the mesh density. The mesh density is critical for both the accurate calculation of the crack tip opening displacements and for capturing the stress gradients of the Hertzian contact. As alluded to, SIFs are calculated by using a displacement correlation technique as described by the equations: 12 22 4(1)IE K uu vr (9 2) 12 22 4(1)IIE K vv vr (9 3) 12 22 4(1)IIIE K ww vr (9 4) However, S IFs could also be determined via J integral decomposition through ABAQUS or the M integral in FRANC3D/NG (Banks Sills) but were not as data obtained PAGE 144 144 provided similar results and is intended to be used in other work reflecting on crack closure under RCF and neither method is currently formulated for accounting for the traction encountered on the crack faces during closure. Quadratic elements are required by the Displacement correlation method for calculating SIFs since the insertion of a quarter point element greatly increases the accuracy of the model without a dramatic increase in mesh density. To reduce the size of the problem a technique referred to as submodeling is utilized where displacements are applied to the boundary of a small cracked block to sim ulate its being a part of a much larger half space. The contact load is applied to the surface with FORTRAN user subroutines DLOAD (and UTRACLOAD) for ABAQUS (see Fig 61). This load is applied up till the edge of the crack but not over the crack, as this would cause the applied pressure to be inaccurate (Fujimoto). Also, crack closure would result and the model would then require a contact algorithm on an already large quadratic, fracture model that would be accurate enough to yield crack tip displacements 6.3 Crack Geometry Surface cracks can have a variety of shapes. Cone cracks have received much attention in the literature (Mackerle). Partial cone cracks are reviewed in separate paper (Levesque and Arakere) and have also been the subject of some analy ses [Hadfield, 1993a, Hadfield, 1993b, Wang and Hadfield, 2000, Zhao, 2006]. Penny cracks seem to be the simplest of which to mesh and simulate in a 3D FEA analysis. However, if attention were limited to the family of possible half penny surface flaws, th ere remain only two features that can be altered, aspect ratio and angle towards the surface. Five aspect ratios have been analyzed ( a/b=0.2, 0.4, 0.6, 0.8, and 1) under a few different angles (0, 45, and 60). PAGE 145 145 6.4 Orientation O f Load As a surface flaw is subjected to rolling contact fatigue, a load is set in a direction across the free surface of the cracked body. In application, this pass of a load will happen repeatedly, subjecting the part to cyclic loading. While at each instant, the loads orientat ion generates a different set of SIFs as it moves along the surface, there is one orientation (for a given load size and magnitude) that produces the highest SIFs along the crack front. It is this orientation that is of most interest for a design analysis that would determine if the SIFs reach a critical value to induce fracture, Kc, initiate fatigue, Keff, or will direct a growth analysis, effK In some orientations of load, the crack faces interpenetrated when contact was not defined, indicating crack closure occurred. The interpolation resulted in a negative KI stress intensity factor but KII and KIII remained identical. Comparing KI along the front of a crack for an orientation with and without contact defined results in fig 62. Re flection indicates that simply setting the KI=0 when KI<0 is conservative, under conditions of crack closure. It must be noted that the coefficient of friction for contact elements between the crack faces was set to zero. For non zero friction between the crack faces it is possible that closure might couple modes of deformation. Contact is computationally very difficult to implement in these problems because traction encountered on the crack faces during closure is not satisfactorily incorporated into the J integral decomposition in ABAQUS and necessitates the use of displacement correlation for evaluating K solutions. However, this procedure has been successfully implemented and will be discussed in a future work on modeling concerns for these types of problems. PAGE 146 146 While the tensile region is maximal on the periphery of the contact patch in an uncracked body (according to the stress solution) this does not necessarily mean that this same position is where the maximum SIFs will occur. So, to find the worst cas e load orientation the distance between the contact patch edge and the crack center at the surface was perturbed and SIF s calculated for these orientations. The SIFs of multiple contact patch p ositions relative to the crack were compared and the position t hat yielded the highest SIF was chosen for curve fitting. These distances for each model are given in Table 1. It was noticed that this distance tends to decrease as load size is increased and as crack aspect ratio decreases (and the crack tip is generally closer to the surface) for the more vertical cracks (as this allows their tip to be more generally located in the small tensile region about the contact). 6.5 Load Rolling contact fatigue can have multiple different shapes of pressure distribution, but if the two contacting bodies can be characterized with two radii of curvature each, as in a Hertzian contact, then the contact patch will be elliptical. For bodies whose radii of curvature are identical the contact patch is circular and its pressure distri bution is 2(,)1oPxypr (9 5) It is important to note that the radial tensile stress region at the edge of contact is both small and quickly decays with increasing depth into the material. Also, it clearly shows that the stress state is quite far from anything that can be considered as far field. The absence of a clearly identifiable far field stress seen by the crack makes the curve fits more complicated than those employed by the widely referenced paper by Newman and Raju (1981) PAGE 147 147 Emp irical equations for circular loading of penny cracks at a few angles are provided. The range of size of circular loads has been nondimensionalized with respect to the crack half width and has been analyzed for the r =1 b 2 b and 3 b cases, which should cover the majority of cases which would be of engineering interest. The data obtained from FEA was fit to the equations: 234 12345 Ioa Kpaaxaxaxaxgf Q (9 6) 234 12345 IIoa Kpbbxbxbxbxgf Q (9 7) 234 12345 IIIoa Kpccxcxcxcxgf Q (9 8) In this way, the results can be easily fit with very little error (~0.5%) and can be easily nondimensionalized. It is noted that when /1 ab Q is approximated by the equation: 1.6511.464 a Qb (9 9) Also, 210.1(1sin) g (9 10) 1 2 2 22cossin a f c (9 11) While the fit uses components of the Newman and Raju fit (1981), it is more involved because of the complex 3D subsurface stress field and the steep stress gradients at the edge of contact to crack tip. The c oefficients for these cases are provided in Table 2 4. PAGE 148 1 48 6.6 Results As a consequence of running these analyses for different aspect ratios and angles a few notable trends can be seen as a result. Firstly, the location along the crack front where the SIFs are the highest can change based on crack depth and angle. Generally, the shallower cracks can have their SIFs toward the crack center while the deeper cracks tend to have their SIFs highest near the surface. For a single load size the worst case crack i s nonobvious. See figs. 63 6 5. By reflecting on the results, it is observed that the =45 a/b=1 case has the highest SIFs for the r=2 b and 3 b cases. In these situations, the highest SIFs occur near the free surface. For the r=1 b case, it is observed th at the =60, a/b =0.2 case has the highest SIFs and this occurs at the deepest portion of the crack. Neither of these observations is intuitive. In fact, they display that it is not always the steepest or the longest of cracks which will have the highest S IFs under a given load, but rather an interesting combination thereof which varies with load size. The variation of SIFs across the crack front also brings many issues in terms of crack growth. For example, if the highest SIFs are produced by a certain geometry crack and these SIF values are largest near the free surface, if this crack grows it will grow wider on the surface rather than into the depth. Then the geometry may become one which yields relatively lower SIFs than its former shape. This means the crack which produces the highest SIFs, for a given geometry and load, may not be the worst geometry crack under RCF because its shape evolution effects its SIFs and so the worst crack geometry is something that can be revisit ed in future work. The above graphs contain every model created for each of the three types of load. The trends between models of the same angle but different depths is quite smooth PAGE 149 149 and is reflective of the transitions of the stress state as a crack is longer. The trends between crac ks of different angles is less obvious as the change in angle places the cracks in a different part of the stressed state. Quite surprisingly, the trends in r =3 b seems smooth between all 15 models graphed even though the load orientation may not be identic al for all models and the tip location may be quite different. This trend is observed because the larger contact patch produces lower stress gradients. 6.7 Conclusions Three dimensional FEA was performed on a parametric variation of semi elliptical flaws under circular RCF patches to provide a set of comprehensive empirical equations for mixed mode SIFs applicable in the hybrid ball bearing industry and other contact fatigue applications. Modeling concerns have been addressed and a modeling framework has been provided for analyzing these types of geometry and load. The parametric variation of load size, crack angle, and crack depth displayed a several interesting trends including the type of cracks that provide the highest SIFs for the given load geometry which is neither the deepest crack nor the steepest crack all the time. The comprehensive and accurate (~0.5%) empirical equations for the KI, KII, and KIII SIFs presented are of immediate engineering relevance to the hybrid silicon nitride ball bearing in dustry towards evaluating critical flaw size, and for developing a fracture mechanics based life prediction methodology. The empirical curve fits will also be of relevance to other areas of component design where contact initiated fatigue damage is importa nt such as gears, roller bearings, and railway wheels. PAGE 150 150 Table 6 1. The xd distances (nondimensinoalized w.r.t. the crack semi width) where KI was maxim al r = 1 b 2 b 3 b a/b =1.0 1.3 1.0 0.875 a/b =0.8 1.3 0.875 0.875 a/b =0.6 1.1 0.875 0.875 a/b =0.4 0.95 0.625 0.75 a/b =0.2 0.95 0.625 0.625 a/b =1.0 1.0 0.625 0.5 a/b =0.8 1.25 0.5 0.5 a/b =0.6 1.25 0.5 0.5 a/b =0.4 1.25 0.5 0.375 a/b =0.2 1.25 0.5 0.375 a/b = 1.0 0.0 0.0 0.0 a/b =0.8 0.0 0.0 0.0 a/b =0.6 0.0 0.0 0.0 a/b =0.4 0.0 0.0 0.0 a/b =0.2 0.0 0.0 0.0 PAGE 151 151 Table 6 2. Coefficients for all parametric cases for KI, KII, and KIII coefficients angle aspect ratio load size ( a ) 5 4 3 2 1 0 K I 1 1 0.00057 0.00152 0.00233 0.00027 0.00209 2 0.00160 0.00449 0.00743 0.00082 0.00529 3 0.00326 0.00928 0.01603 0.00172 0.01017 0.8 1 0.00631 0.00882 0.00516 0.00059 0.00240 2 0.01813 0.02534 0.0 1571 0.00173 0.00654 3 0.03734 0.05212 0.03330 0.00359 0.01313 0.6 1 0.01528 0.02003 0.00842 0.00108 0.00250 2 0.04472 0.05859 0.02521 0.00318 0.00712 3 0.09423 0.12337 0.05377 0.00671 0.01478 0.4 1 0.01711 0.02316 0.00882 0.00120 0 .00232 2 0.05121 0.06928 0.02663 0.00361 0.00688 3 0.10771 0.14566 0.05629 0.00760 0.01439 0.2 1 0.00905 0.00734 0.00285 0.00021 0.00181 2 0.02741 0.02232 0.00860 0.00063 0.00546 3 0.05815 0.04745 0.01823 0.00135 0.01153 K II 1 1 0 .00005 0.00059 0.00135 0.00011 0.00043 2 0.00011 0.00167 0.00400 0.00030 0.00130 3 0.00027 0.00369 0.00867 0.00063 0.00286 0.8 1 0.00167 0.00220 0.00033 0.00009 0.00027 2 0.00517 0.00691 0.00083 0.00031 0.00078 3 0.01093 0.01453 0.00185 0.00065 0.00167 0.6 1 0.00260 0.00398 0.00072 0.00020 0.00015 2 0.00788 0.01206 0.00224 0.00062 0.00042 3 0.01667 0.02550 0.00476 0.00132 0.00086 0.4 1 0.00058 0.00043 0.00017 0.00006 0.00007 2 0.00190 0.00114 0.00055 0. 00018 0.00018 3 0.00423 0.00215 0.00124 0.00036 0.00037 0.2 1 0.00502 0.00535 0.00209 0.00023 0.00001 2 0.01555 0.01661 0.00647 0.00071 0.00003 3 0.03324 0.03555 0.01384 0.00152 0.00006 K III 1 1 0.00072 0.00206 0.00129 0.00047 0.0 0003 2 0.00215 0.00613 0.00383 0.00136 0.00009 3 0.00457 0.01299 0.00804 0.00290 0.00019 0.8 1 0.00349 0.00719 0.00367 0.00026 0.00002 2 0.01045 0.02144 0.01095 0.00071 0.00007 3 0.02215 0.04527 0.02312 0.00145 0.00015 PAGE 152 152 Table 6 2. Continued coefficients angle aspect ratio load size ( a ) 5 4 3 2 1 0.6 1 0.00218 0.00543 0.00275 0.00041 0.00001 2 0.00651 0.01619 0.00823 0.00121 0.00004 3 0.01380 0.03416 0.01738 0.00250 0.00009 0.4 1 0.00440 0.00 286 0.00054 0.00079 0.00000 2 0.01361 0.00898 0.00172 0.00242 0.00001 3 0.02883 0.01913 0.00369 0.00510 0.00001 0.2 1 0.01197 0.01271 0.00448 0.00119 0.00001 2 0.03696 0.03926 0.01385 0.00369 0.00004 3 0.07848 0.08338 0.02941 0.0 0783 0.00009 PAGE 153 153 Table 6 3. Coefficients for all parametric cases for KI, KII, and KIII coefficients angle aspect ratio load size (a) 5 4 3 2 1 45 K I 1 1 0.03069 0.06817 0.05008 0.01249 0.00783 2 0.05798 0.12787 0.10044 0.02374 0.01392 3 0.06165 0.13609 0.11087 0.02560 0.02027 0.8 1 0.25896 0.37265 0.16790 0.02668 0.00796 2 0.47901 0.68372 0.32295 0.04879 0.01569 3 0.50815 0.72296 0.34653 0.05168 0.02271 0.6 1 0.01294 0.01710 0.00755 0.00106 0.00249 2 0.32412 0.36314 0.13671 0.01950 0.01478 3 0.12658 0.08750 0.02082 0 .00246 0.02104 0.4 1 0.08439 0.06357 0.00203 0.00421 0.00718 2 0.14478 0.15116 0.04963 0.00780 0.01204 3 0.09795 0.10360 0.03560 0.00500 0.02300 0.2 1 0.02247 0.00565 0.02538 0.00281 0.00998 2 0.01594 0.02200 0.01337 0.00007 0.01610 3 0.03883 0.04039 0.01662 0.00126 0.02068 K II 1 1 0.02901 0.05987 0.07188 0.01255 0.04438 2 0.02804 0.06717 0.09026 0.01562 0.05196 3 0.02127 0.06061 0.08914 0.01507 0.04856 0.8 1 0.15705 0.17951 0.14463 0.01446 0.04219 2 0.08488 0.07200 0.10265 0.00823 0.04500 3 0.00015 0.02817 0.05834 0.00234 0.03989 0.6 1 0.00480 0.00766 0.00272 0.00063 0.00019 2 0.97142 1.49346 0.60414 0.10156 0.03866 3 1.02226 1.53021 0.62327 0.10318 0.03278 0.4 1 0.08691 0.12303 0.132 84 0.00258 0.02877 2 0.01543 0.05833 0.02924 0.00684 0.02684 3 0.05186 0.12077 0.00423 0.00754 0.01685 0.2 1 0.17706 0.30827 0.18665 0.00855 0.01485 2 0.06102 0.06311 0.02082 0.00388 0.00863 3 0.12743 0.16371 0.03198 0.00820 0.00 413 K III 1 1 0.02890 0.10038 0.09167 0.00947 0.00268 2 0.02893 0.11510 0.10462 0.01713 0.00316 3 0.02652 0.10900 0.09448 0.01920 0.00290 0.8 1 0.09734 0.25130 0.16807 0.00512 0.00214 2 0.03831 0.21799 0.16516 0.01420 0.00229 3 0.01420 0.18514 0.14208 0.01526 0.00202 PAGE 154 154 Table 6 3. Continued coefficients angle aspect ratio load size ( a ) 5 4 3 2 1 0.6 1 0.00117 0.00032 0.00057 0.00048 0.00001 2 0.10361 0.00956 0.05369 0.01667 0.00096 3 0.04763 0.1030 9 0.09433 0.00746 0.00105 0.4 1 0.25381 0.26703 0.08972 0.03561 0.00028 2 0.39896 0.42784 0.15003 0.04753 0.00017 3 0.56860 0.61293 0.22239 0.05780 0.00015 0.2 1 0.17767 0.13319 0.05089 0.05330 0.00111 2 0.16738 0.15268 0.04809 0.03087 0.00064 3 0.17233 0.16379 0.05305 0.02531 0.00044 PAGE 155 155 Table 6 4 Coefficients for all parametric cases for KI, KII, and KIII coefficients angle aspect ratio load size (a) 5 4 3 2 1 60 K I 1 1 0.00022 0.0021 9 0.00166 0.00151 0.01305 2 0.00154 0.00912 0.00541 0.00388 0.02207 3 0.00407 0.01589 0.00991 0.00569 0.02850 0.8 1 0.00580 0.00104 0.00646 0.00132 0.01405 2 0.00649 0.02848 0.01860 0.00504 0.02342 3 0.02152 0.05393 0.03099 0.00769 0.02959 0.6 1 0.00992 0.02657 0.02304 0.00145 0.01485 2 0.04943 0.07986 0.04148 0.00617 0.02394 3 0.07570 0.11359 0.05482 0.00858 0.02964 0.4 1 0.19388 0.20215 0.05007 0.01112 0.01417 2 0.18599 0.22079 0.07396 0.01166 0.02164 3 0. 12169 0.15525 0.05367 0.00837 0.02604 0.2 1 0.13338 0.10915 0.00607 0.00692 0.01262 2 0.08585 0.09914 0.02550 0.00562 0.01888 3 0.01541 0.03383 0.00914 0.00297 0.02266 K II 1 1 0.00062 0.00185 0.00699 0.00085 0.01279 2 0.00082 0.0 0312 0.00757 0.00122 0.00995 3 0.00176 0.00287 0.00670 0.00122 0.00437 0.8 1 0.00520 0.02522 0.00543 0.00325 0.01226 2 0.01770 0.04188 0.00250 0.00464 0.00800 3 0.02757 0.04846 0.00790 0.00505 0.00245 0.6 1 0.04532 0.07520 0.00876 0.00680 0.01044 2 0.04845 0.08291 0.01570 0.00701 0.00509 3 0.04929 0.07602 0.01641 0.00641 0.00023 0.4 1 0.15883 0.19031 0.11165 0.00839 0.00678 2 0.11770 0.17533 0.04090 0.00955 0.00119 3 0.31135 0.40073 0.13295 0.02029 0.0031 7 0.2 1 0.36183 0.47099 0.22470 0.02001 0.00731 2 0.05244 0.05544 0.05487 0.00289 0.00191 3 0.15749 0.19656 0.04810 0.00870 0.00180 K III 1 1 0.00218 0.00133 0.00088 0.01215 0.00111 2 0.00378 0.00180 0.00232 0.01097 0.00098 3 0 .00263 0.00168 0.00087 0.00567 0.00056 0.8 1 0.03250 0.03575 0.01271 0.01721 0.00131 2 0.04629 0.04433 0.01722 0.01513 0.00107 3 0.03760 0.02715 0.01109 0.00829 0.00060 PAGE 156 156 Table 6 4. Continued coefficients angle aspect ratio load size ( a ) 5 4 3 2 1 0.6 1 0.06033 0.05899 0.01928 0.02105 0.00187 2 0.06898 0.05935 0.02020 0.01772 0.00148 3 0.05150 0.03334 0.01189 0.01070 0.00087 0.4 1 0.37372 0.38473 0.14244 0.05317 0.00015 2 0.60388 0.64698 0.23694 0.05732 0.00054 3 0.59932 0.63865 0.23462 0.04821 0.00062 0.2 1 0.49753 0.49000 0.18141 0.06698 0.00012 2 0.82548 0.87869 0.31396 0.07222 0.00066 3 0.82666 0.88589 0.31707 0.06345 0.00078 PAGE 157 157 6 8 Figures Figure 61. Model configuration displaying orientation of load and defining variables and dx PAGE 158 158 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 2 0 2 4 6 8 10 x 103 KI with and without contact defined phiKI* without contact with contact Figure 62. KI for a single load geometry with and without contact defined. PAGE 159 159 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 I yg p KI* a/b=1.0 =0 a/b=0.8 =0 a/b=0.6 =0 a/b=0.4 =0 a/b=0.2 =0 a/b=1.0 =45 a/b=0.8 =45 a/b=0.6 =45 a/b=0.4 =45 a/b=0.2 =45 a/b=1.0 =60 a/b=0.8 =60 a/b=0.6 =60 a/b=0.4 =60 a/b=0.2 =60 Figure 63. KI for r=1 b along the crack front for all variations of crack angle and depth. 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 KI for r2b load for varying aspect ratiosKI* a/b=1.0 =0 a/b=0.8 =0 a/b=0.6 =0 a/b=0.4 =0 a/b=0.2 =0 a/b=1.0 =45 a/b=0.8 =45 a/b=0.6 =45 a/b=0.4 =45 a/b=0.2 =45 a/b=1.0 =60 a/b=0.8 =60 a/b=0.6 =60 a/b=0.4 =60 a/b=0.2 =60 Figure 64. KI for r=2 b along the crack front for all variations of crack angle and depth. PAGE 160 160 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 I o 3b oad o ayg aspect atosKI* a/b=1.0 =0 a/b=0.8 =0 a/b=0.6 =0 a/b=0.4 =0 a/b=0.2 =0 a/b=1.0 =45 a/b=0.8 =45 a/b=0.6 =45 a/b=0.4 =45 a/b=0.2 =45 a/b=1.0 =60 a/b=0.8 =60 a/b=0.6 =60 a/b=0.4 =60 a/b=0.2 =60 Figure 65. KI for r=3 b along the crack front for all variations of crack angle and depth. PAGE 161 161 CHAPTER 7 UNCERTAINTY ANALYSIS FOR FATIGUE FAILURE PROBABILITY OF SILICON NITRIDE ROLLING ELEMENT 7.1 Introduction Ball and roller bearings are widely used in a variety of industrial machinery to allow relative motion and support load in rotating shafts. In conventio nal ball bearings, with metal raceways and balls, subsurface originated spalling and surface originated pitting have been recognized as the dominant modes of failure due to rolling contact fatigue (RCF). Bearing subsurface material is subjected to RCF cyc les that induces a complex triaxial stress state, nonproportional loading, high hydrostatic stress component, and changing planes of maximum shear stress during a loading cycle, and eventually leads to a subsurface crack. Spalling occurs when subsurface cr acks propagate towards the s urface to form a surface spall ( Sadeghi et al ., 2009) This mechanism is the dominant mode of failure in rolling element bearings that have smooth surfaces and operate under elastohydrodynamic lubrication (EHL) conditions. Sur face originated pitting occurs in cases where surface irregularities in the form of dents or scratches are present. Here, cracks initiate at the surface stress concentrators and thereafter propagate at a shallow angle to the surface (Bower, 1988). This me chanism of failure is more common in gears where substantial sliding occurs between the contacting surfaces. Aircraft engine manufacturers have been aggressively pursuing advanced materials to meet mainshaft bearing requirements of advanced engines for m ilitary, commercial and space propulsion. These requirements include bearings with extended life, superior corrosion resistance, surface durability and tribological performance. Hybrid bearings utilize silicon nitride balls and steel raceways to achieve c onsiderably PAGE 162 162 longer fatigue lives, to have similar thermal behavior, to last up to five times longer in oil starvation conditions (Miner et al. 1996, Tanimoto et al., 2000) and perform well under corrosive conditions (Klemm, 2002). Silicon nitride balls h ave many physical properties which allow for bearing technology advancement including low density and high compressive strength but also have low fracture toughness of 46 MPa m (Piotrowski and OBrien, 2006). This low fracture toughness, in combination with unavoidable ball to ball collisions in the manufacturing (lapping) process often results in ring or c cracks (partial cone cracks) (Cundill, 1996). These flaws can grow under RCF when placed in service and result in a spall on t he ball surface, (Fig 71) (Levesque and Arakere, 2008). The failures caused by surface cracks of silicon nitride ball bearings under rolling contact have been addressed by ring or c cracks and wedge effect models. In the c crack model the Hertzian contact stresses occurring around the circumference of the bearing contact area is thought to dominate the crack growth. The lubricant is assumed to have no effect on the crack growth. In the wedge effect model the fluid pressure is thought to penetrate into the crack by contact pressure. When the fluid is pressurized by the maximum contact pressure at the contact center it is thought to cause crack growth. The mechanisms of these two models have been investigated separately. A surface crack is likely to be affected by the stresses used in both models as a ball passes over the crack The stress intensity factor (SIF) at the crack tip is also adversely affected by the presence of frictioninduced traction forces at the contact. However, effects of the fluid pressure are not well understood or characterized yet and the dominant mode of failure is thought to be crack growth driven by SIFs arising from the contact patch passing over the c crack. Partial cone or c cracks are considered the most damaging PAGE 163 163 surface de fect that limits ball life in hybrid bea rings under service conditions ( Eva ns, 1983, Hadfield et al., 1993b) The fatigue damage process in ceramic rolling elements is very different from metal balls. RCF in metal bearings is manifested as a flaking off o f metallic particles from the surface of raceways and/or rolling elements. As described earlier, this process commences as a crack below the surface and is propagated to the surface, eventually forming a pit or spall in the loadcarrying surface. Bearing f atigue life estimation is still largely based on the seminal probabilistic life model by Lundberg and Palmgren (LP) proposed in 1945. Despite many improvements to the LP model current probabilistic bearing life prediction methods are based on the ISO stan dard set up in 1989 and continue to rely on extensions to the LP model, are empirical in nature, and include variables that are obtained from extensive experimental testing. The LP theory states that for bearing rings subjected to N cycles of repeated load ing the probability of survival S is given by, 0 01 lnec hNV A Sz (7 1) where, 0 is the maximum orthogonal shear stress in the contact, z0 is the corresponding depth at which this stress occurs and V is the stressed volume of material. T he parameters A, c and h are material characteristics that are determined experimentally. The parameter e is the Weibull slope for the experimental life data plotted on a Weibull probability paper. Metals are weaker in shear than tension. In contrast to metallic materials, Si3N4 material is weaker in tension than compression or shear. The LP methods cannot account for this difference in material behavior. Furthermore, c cracks already exist on PAGE 164 164 the ball surface and hence subsurface crack initiation and subsequent growth to a spall is not pertinent to ball failure. A fracture mechanics methodology that evaluates the critical flaw size (or the largest allowable flaw size such that growth does not occur under the service conditions) in balls is required, since failure typically occurs, for hybrid bearings, when the effective SIF at the crack tip equals the mixedmode fracture toughness of the silicon nitride ball material. The larg est allowable surface flaw that does not propagate under RCF loading is of design significance, and is termed the critical flaw size (CFS). A systematic procedure to compute the CFS based on fracture mechanics principles, RCF loading, and ball material properties has been recently presented by Levesque and Arakere (2009). Nondestructi ve evaluation (NDE) techniques are being developed and used on each ball in order to determine if it is acceptable to enter service (Wang, 2000). The resolution to which each ball is inspected, as determined by the CFS, has a strong effect on the cost of each ball and has been analyzed deterministically in prior work (Levesque and Arakere, 2009). The cost of NDE method for silicon nitride balls scales up very steeply with decreasing CFS and increasing ball diameter. Thus the cost associated with NDE can become a significant fraction of the overall manufacturing cost of the silicon nitride ball. Brittle materials inherently exhibit considerable variation in material mechanical and fractureproperties. For example, there is a significant variation in the d epth of a crack from a controlled impact (Lawn, 1967, Wang and Hadfield, 2000). Also there is a noted variation in the fracture toughness specimen to specimen (for example 4.85 0.36 MPa PAGE 165 165 7.2 Modeling Rolling Contact Fat igue Orientation Effects Finite element analysis techniques have been developed to analyze 3D nonplanar surface flaws for stress intensity factor calculation (SIF) (Levesque and Arakere, 2008). Mesh density must be accurate enough to account for the strong gradients in stress from the contact patch and element types must allow SIF extraction SIF calculation is done by the displacement correlation method, which can be written as: 12 22 4(1)I dE K uu r (7 2a) 12 22 4(1)II dE K vv r (7 2b) 12 22 4(1)III dE K ww r (7 2c) It is noted that rd is the radius from where the relative displacements are measured from the crack tip, and u v and w correspond to the displacements in the x y and z directions for the local coordinate system at the tip of the crack. depth angled 30 to the surface. These specific dimensions were chosen since they have been experienced to be a set of dimensions that have been observed in experiment (Wang 2000) and have exhibit high stress intensity factors when compared to other geometries (Levesque and Arakere, 2009). Loading is done via FORTRAN user subroutines DLOAD and UTRACLOAD for ABAQUS (Dessault Syst ms, 2007) for the normal and traction loads. The equations to describe the pressure distributions can be written as: PAGE 166 166 22 22(,) ``dd oxxyy pxyp ab (7 3) and (,)(,) fxypxy (7 4) where xd and yd are the distances from the global coordinate system to the load center, a and b ar e the dimensions of the ellipse along these dimensions, and is the friction coefficient for the moving load in a full slip interaction (see Fig. 72). For the specific cases that were analyzed, were used, which is believed to be representative of lubricated bal l bearing contacts. The stress intensity factors are calculated in all three modes. For implementation in determining what combination of these three modes of crack displacement have met a fatigue threshold, or Kth, the critical strain energy release rat e (CSERR) (Levesque et al. 2009) is utilized as adapted from Anderson (2005) and can be written as. 222(1)eqIII IIIKKKvK (7 5) For simplification, the SIFs (calculated from displacement correlation) as calculated as a function of crack tip p osition a re combined them into a single parameter (via CSERR) and have then chosen the highest value along the crack front as the representative value to determine if crack growth can occur for each orientation. For an uncertainty analysis, the highest Keq for a randomly oriented crack relative to the moving contact patch must be calculated. Here a few assumptions are made. The contact patch passes on the ball surface in a band and does not change direction fo r the entire sequence (see Fig 73 ). The load magnitude and the traction magnitude and direction do not change relative to the elliptical contact. As the load passes around PAGE 167 167 the ball, there exists only one position where the max SIF is reached and this orientation is the nearest to the contact patch and on the opposite side of the traction direction. This maximum position can be classified with two variables: the lateral position of the crack relative to the center of the contact patch, xD, and the angle that the crack makes on the surface relati ve to the contact patch tangent, (Herein, the distance of the crack to the center of the contact load in the yD direction is always 0 for maximum SIFs.) (see Fig. 7 3) Parametric studies for these two parameters have been analyzed and their effect s se parated as independent parameters in Figs 74 and 75. 7.3 Uncertainty Analysis In this section, uncertainty analysis of estimating the failure of silicon nitride ball bearings is presented using surrogate modeling and Monte Carlo simulation. The fatigue f ailure model in the previous section is based on the effective stress intensity factor, which depends on the size, direction and location of a crack; fracture toughness of the material; and applied loads. However, it is difficult to determine if a bearing may fail or not because there are many uncertainties in the practical operating environment. For example, an existing crack may not fail if its location is out of the contact path with bearing track or its size is too small. Thus, the failure of a bearing can only be evaluated in terms of probability. The goal of this section is to evaluate the probability of fatigue failure of silicon nitride ball bearings under representative operating conditions. 7.3.1 Uncertainty Quantification The first step in uncertainty analys is is to quantify the uncertainty of the input parameters. Table 1 summarizes input parameters that are used in uncertainty analysis. Some parameters are considered to be deterministic, while the others are random. Even if all bearings do not ha ve the same diameter due to manufacturing tolerances, it PAGE 168 168 is considered a deterministic value because its uncertainty is relatively small compared to others and its contribution to the uncertainty in fatigue failure is also small. In general, the uncertaint y in the applied pressure is large. In fact, it is known that the applied load is the largest source of uncertainty. However, a deterministic value of the maximum applied pressure is used in uncertainty analysis due to the following two reasons: first, it is difficult to characterize probabilistic distribution of the applied load, and second, it will provide conservative estimate of the probability of failure. Major uncertainties are related to the crack configuration. Based on inspection data (Cundill, 1 996), a new silicon nitride ball bearing can have many micro cracks on the surface. Levesque and Arakere (2008) showed that these cracks can be modeled using half penny shaped cracks with a specific size, aspect ratio and orientation. Since the manufacturi ng process does not have any preference in the orientation of nucleated cracks and ball orientation in the assembly can vary often, it is assumed that they are uniformly distributed between 0 and 180 (see fig 72) For the same reason, the location of cr ack is also uniformly distributed on the ball surface. Based on sample measurements by (Wang, 2000), it is reported that the initial crack size can be between 200 and 350 microns. Since no probabilistic distribution information is available, it is assumed that the initial crack size is uniformly distributed between the minimum and maximum sizes. Piotrowski and OBrien (2006) have computed the fracture toughness value for the silicon nitride ball bearings when an applied force generates a tensile hoop stre ss. The fracture toughness value, computed from 16 experiments shows a uniform variation within 12% of the nominal value. Thus, the uncertainty of the fatigue threshold is PAGE 169 169 modeled using a uniform distribut ion with the nominal value of 2.8 MPa he suggested range of 24 MPa 7.3.2 Surrogate Modeling Uncertainty Propagation In the uncertainty propagation stage, the input uncertainties are propagated through the governing physics of the system to yield uncertainty in output, which is the effective stress intensity factor. Traditionally, Monte Carlo simulation (MCS) is often employed for this purpose where many samples of input random variables are generated accord ing to their distribution types and samples of output variable are produced by solving the governing equation with each set of input variables. This stage is computationally intensive because it involves 3D finite element modeling and volume integrals. For example, for given values of input parameters, the computation of eff ective stress intensity factor takes about 144,000 sec in desktop computer with 2 processors. In order to identify the effect of input uncertai nty on the output uncertainty numerous repetitions of this calculation are required, which can quickly become imp ractical. Thus, the key issue in uncertainty analysis is how to effectively propagate the input uncertainty to the output uncertainty. There are many methods available in uncertainty propagation that use approximation. First order reliability method (FORM ) (Allen and Camberos, 2009), second order reliability method (SORM) (Allen and Camberos, 2009), importance sampling method (Cano et al., 1996), and MCS are a short list of available methods. Except for MCS, all other methods use approximation in either ou tput variable itself or its distribution type, which inevitably involves error especially when the governing system equation is complex. Herein, surrogate modeling techniques (FAC Viana, 2009) are utilized to approximate the relation between input parameters and output. Instead of PAGE 170 170 black box type surrogate model, physics based surrogate modeling is employed to represent the relationship more effectively. Once the surrogate model is obtained, uncertainty in output variable can easily be calculated using MCS b ecause the computational cost of function evaluation using the surrogate model is negligible. It is observed that the effective stress intensity factor for cracks on ball bearings can be characterized with four parameters: applied contact pressure, lateral position of crack with respect to contact patch, orientation of crack with respect to contact path, and initial crack size. Thus, the uncertainty propagation can be written in a general form of: 0(,,,)eq DKfpxa (7 6) where p0 is the applied pressure at the center of contact patch, xD is the normalized lateral position between the center of contact patch and the center of crack, is the angle between a line perpendicular to the contact path and the orientation of crack, and a is the length of surface crack. Although the explicit expression of function f is unknown (or, sometimes it is an implicit function), it can be evaluated fo r given input parameters using finite element analysis and volume integrals. The idea of surrogate modeling is to approximate the function f using a simple analytical function (Queipo et al., 2005). The general procedure of surrogate modeling is to gener ate several samples, called design of experiment, and to fit a function using these samples. First, a number of samples are chosen based on different criteria, such as full factorial design or Latin hypercube sampling. The larger the sample size, the bette r the quality of approximation. However, each sample means expensive finite element analysis. In order to have a reasonable accuracy, the number of samples for PAGE 171 171 the case of four variables will be around 100. The number of required samples exponentially incr eases as the number of input variables increases. Once all samples are available, they are used to fit a function. In general, the functional form is fixed with unknown coefficients, which are to be found by minimizing the error between the function and samples. In general, surrogate modeling does not require detailed knowledge on the physical problem; it can be considered as a black box. Herein however, the surrogate model is simplified by observing the physical behavior of the system, in which the effec tive stress intensity factor can be expressed as a function of the four parameters in the following form: 0 2 1() 2700MPa2.5250 eq Dp f a K fx (7 7) First, the stress intensity factor, Keq, is proportional to the contact pressure because the material response i s linear elastic and stress increases linearly with the applied load. In addition, Keq is proportional to the square root of crack size, which is basically identical to the traditional definition of stress intensity factor. Since these two relations are si mple and explicit, there is no need to introduce approximation using surrogate modeling. On the other hand the effects of lateral position and orientation are not straightforward, and thus, require approximation as in Eq. (7 7) An important point in the above equation is that the effects of these two parameters are decoupled. This can be easily explained that if two cracks are located in different lateral positions, the failure mechanism will still be the same except that the magnitude of contact periphery stress will be di fferent In the above equation. I t is assumed that the aspect ratio of the PAGE 172 172 semi elliptical crack, the cracks angle to the surface and that the elliptical contact patch aspect ratio remain constant. After simplifying the relationship betw een input variables and output, the initial surrogate model with four variables can now be simplified to two surrogate models with a single variable, which is much more computationally efficient to build and more accurate. To simulate the variation in SIFs due to orientation, seven equally spaced simulations are generated for the orientation angle and five for the lateral position within their ranges and are interpolated between for computational accuracy. The next step is to choose a surrogate model. The difficulty is that there is no single surrogate that outperforms all others. Depending on functional behavior, one surrogate performs better than the others. In general, however, the functional behavior is unknown a priori. One of the best practices is to build multiple surrogates using the sampled data and choose the best one. It is generally accepted that cross validation (Myers and Montgomery, 2002) is a good tool to choose the best surrogate. This procedure is relatively inexpensive under the assumption that obtaining a sample requires expensive finite element analysis, but surrogate modeling can be finished without requiring intensive computation. In this chapter four different surrogates are considered: Kriging, Radial Based Neural Network (RBNN), Support Vector Regression (SVR), and fourth order polynomial response surface (PRS). Different surrogates have different characteristics. For example, the Kriging always pass the sampled data points, while PRS does not pass the data point exactly. However, th at does not mean that the former is more accurate than the latter. The accuracy of a surro gate should be measured at data points that are not used in fitting the surrogate. PAGE 173 173 Figures 7 4 and 7 5 show the approximation of f1(xD) and f2( ), respectively, usi ng four different surrogate models. Since all four surrogates are close to each other, it is difficult to tell which surrogate is the best. Without having additional test points, crossvalidation can be used to find the best surrogate. In cross validation, one data point is dropped in fitting the surrogate and the error is measured at the dropped data point, which is called a prediction error. If this procedure is repeated for all data points, the root mean square of prediction errors (PRESSRMS) can be used as an indicator of accuracy. PRESSRMS is a well established parameter to compare the effectives of surrogates. The smaller the PRESSRMS value is, the more effective the surrogate is. Table 7 1 compares the PRESSRMS values for all four surrogate models. It shows that Kriging is a good model to fit the lateral position, while PRS is good for the orientation. 7.3.3 Uncertainty Analysis Monte Carlo Simulation Using Surrogate Once surrogate models are selected, it is used to evaluate the uncertainty of stress intensity factor according to the uncertainty in input parameters. Since function evaluation in the surrogate model is very fast, MCS can be used for that purp ose. The procedure is to generate many samples of input parameters according to their distributi on types and to apply them to the surrogate model to generate samples of stress intensity factors. From the data in Table 7 1, it is relatively easy to generate samples of uniformly distributed crack sizes and orientations using a random number generator. However, generating samples of lateral positions is not straightforward because they are uniformly distributed on a sphere. The method to uniformly distribute points on a sphere is based on Archimedes theorem, which is stated as the axial projection of any measurable region on a sphere on the right circular cylinder circumscribed about the sphere preserves area. The PAGE 174 174 physics behind this theorem is explained in the Shao and Badler (1996). According to the theorem, two independent uniformly distributed rand om variables, z ~ U [ 1, 1] and ~ U [0, 2 ] are sampled based on their distribution types. Each combination of ( z ) corresponds to a sample of a point on the surface of the sphere, whose coordinates are given by 22(,,)1cos,1sin, xyzzzz (7 8) Th is method uniformly distributes points on the surface of a unit radius sphere. The normalized lateral position xD in Eq. (7 8) is equivalent to x in Eq. (7 6) if the coordinate system is set such that the center line of contact patch is on the yzplane. Fo r given dimensions of ball bearing in Table 1, the crack will be located within the contact patch if 'Db x d (7 9) During MCS, a crack is located randomly on the surface of the ball bearing, whose lateral position is calculate d from Eq. 7 8. If the crack lies within the contact patch according to Eq. 79, the effective stress intensity factor is calculated using Eq. 77. The objective is to calculate the probability of fatigue failure of the silicon nitride ball bearing under input uncertainties described in Table 7 1. Here, the failure mode is defined when the effective stress intensity factor is larger than the fatigue threshold. For that purpose, a limit state function is first defined as eqthgKK (7 10) and the probability of failure is defined by Prob(0)FPg (7 11) PAGE 175 175 In MCS, the probability of failure is calculated by counting the number of samples that are failed. When the total number of random samples is N the probabilit y of failure can be calculated by 11 ()N F Fj jN PIg NN (7 12) where Ij is an index function whose value is one when the argument is positive and zero otherwise, and NF is the number of failed samples. Due to the random nature of MCS, differen t sets of samples yield different values of probability of failure; exact PF can only be calculated when the number of samples is infinity. The standard deviation of the probability of failure can be estimated using (1)FF PFPP N (7 1 3) Note that the standard deviation is inversely proportional to the number of samples. Table 7 2 shows the results of uncertainty analysis using MCS and surrogate modeling. In order to see the effect of samples, two different sets of samples are used: o ne with 105 and the other with 106 samples. Samples of the uncertain input parameters are randomly generated according to their dist ributions and output samples of effective stress intensity factors are calculated using Eq. (77). Since the fracture toughn ess is also random, the samples of effective stress intensity factors are compared with that of fracture toughness to determine the number of failure cases. In the table, NC is the number of samples in which cracks are located within the contact patch, and NF is the number of samples that fail. It is clear that the uncertainty in the probability, PF, decreases as the number of samples increases. It can be concluded that for given PAGE 176 176 uncertainties in input parameters, the chance that a ball bearing may fail is about 0.5%. In practice, the failure probability will be lower because the maximum applied load is used. Figure 76 shows the plot of cumulative distribution function (CDF) of the limit state function g = Keq Kth for the case of 105 samples. In order t o emphasize the failure probability, 1 CDF is plotted in logscale. The ordinate value at g = 0 is the value of failure probability. The discontinuity in the slope near g = 2.5 is due to the fact that those cracks lying outsid e the contact patch do not fail and hence, do not contribute the stress intensity factor calculation. As can be observed in the plot, the slope of the CDF curve near g = 0 is high, which means that the failure probability can be improved significantly by slightly increasing shiftin g the limit state function. This can be achieved in various ways, such as using a material with higher Kth or reducing initial crack size that can reduce Keq. 7.3.4 Parameter Study Even if the accuracy of the failure probability depends on uncertainty quantification of input parameters, it can still provide the possibility of improving the failure probability. An important question is how much the failure probability can be improved by modifying the input uncertainty. For example, the initial distribution o f crack size was uniform between 200 and 350 microns. If the manufacturing technology is improved such that the largest initial crack size is reduced by 10%, then it would be beneficial to estimate how much the probability of failure can be improved. In pr actice, it is not feasible to change the lateral position and the orientation of cracks as they are random by nature. However, the initial crack size and the fracture toughness can be modified by using different manufacturing technology or different materi als. PAGE 177 177 Table 7 3 compares the probability of failure when the initial distribution of crack size is reduced from U [200, 350] (original) to U [200, 335] (improved). Both cases use the same number of samples, 106. Note that reducing uncertainty in initial cra ck size by 10% improves the probability of failure by 20%. A similar trend can be observed when the fracture toughness is increased from U [2.464, 3.136] to U [2.531, 3.360]. The effect of this improvement on the probability of failure is tabulated in Table 7 4, in which the probability of failure is improved by 20%. It is noted that improvements on either the distribution of crack size or the distribution of fracture toughness has similar effect on the probability of failure. Table 7 5 shows the combined improvements when both the initial crack size and the fracture toughness are improved. It turns out that the probability of failure can be improved by 44% due to this change. Figure 76 plots the variation of the probability of failure with the percentage re duction of crack size. It is noted that a higher reduction in the crack size leads to lower probability of failure. For a 40% reduction in crack size, original materials probability of failure could improve by 71%, while for a tougher material, the improv ement in probability of failure is about 77%. In order to isolate the effects of location and orientation of crack, the probability of failure is computed for a deterministic crack size. Figure 77 plots the variation of the probability of survival as a function of the crack size. Low crack sizes have a very high probability of survival. As the crack size increases, the probability of survival decreases and the decrease is less for a tougher material. Piotrowski and OBrien (2006) has computed the fracture toughness for silicon nitride ball bearings when the applied force produces hoop stress. Assuming a very PAGE 178 178 tough material that its fracture toughness values conform to the range Piotrowski and OBrien (2006) mentioned at 12% unifor m variation about 4.5 MPa m an attempt is made to find the variation of probability of survival as a function of crack size. Figure 78 plots that variation. Tougher material allows the presence of higher crack lengths. Figures 76 7 7 and 78 could be used in good conjunction to optimize the combination of material and NDE techniques cost with the safety of the ball bearings. 7.4 Conclusions This represents the first uncertainty analysis of probability of survival of ceramic roll ing elements involving a statistical investigation on ball survival rates dependent on fatigue threshold, crack size, and crack orientation in a rolling contact fatigue system. The main conclusions of the work are as follows: Finite element models were us ed to calculate stress intensity factors for surface cracks under rolling contact fatigue. An empirical equation was developed for a mixedmode fracture parameter to account for the variation in crack position, load magnitude, and crack size. Statistical t echniques have been developed to determine the probability of ball failure under service conditions. As the crack size increases, the number of possible orientations for failure increases gradually. Decreases in the upper bound of allowable crack sizes has a similar effect as increasing the fatigue threshold of the material. For high fatigue thresholds, the possibility of ball failure is low but increases steadily for larger crack sizes. PAGE 179 179 7.5 Figures Figure 71. Before and after of a partial cone crack placed under rolling contact fatigue [Image courtesy of The Timken Company]. PAGE 180 180 Figure 72. Elliptical Hertzian contact load coordinate system and variables for elliptical load orientation simulation. Figure 73. The elliptic contact patches on the ball surface and the band which they remain in. Cracks inside and very close to this region (assuming a fixed contact angle and attitude of ball rotation) will experience RCF. PAGE 181 181 Table 7 1 Input parameters and their distributions for uncertainty analysis Para meter Type Value (or distribution) Diameter of ball, d Deterministic 25 mm Width of contact patch, b Deterministic 8 mm Pressure, p 0 Deterministic 2.7 GPa Friction coefficient, Deterministic 0.07 Fracture toughness, K th Ran dom U [2.464, 3.136] MPa Crack length, a Random U [200, 350] m Crack orientation, Random U [0, 180] Crack position, x D Random Uniformly distributed on the sphere Figure 74. Surrogate models for the contribution of lateral position to the stress intensity factor PAGE 182 182 Figure 75. Surrogate models for the contribution of crack orientation to the stress intensity factor Table 7 2. Cross validation errors (PRESSRMS) for different surrogate models Model PRESSRMS Lateral position Angular position Kriging 0.1187 0.3921 RBNN 0.2597 0.1895 SVR 0.6224 0.4388 PRS 0.2308 0.1658 Table 7 3. Probability of failure values for different sample sizes N N C N F P F PF 100,000 67,816 557 5.57E 03 2.35E 04 1,000,000 679,338 5,990 6.00E 03 7.71E 05 PAGE 183 183 Figure 76. Cumulative distribution function of the limit state function g = Keq Kth Table 7 4. Effects of improvements pertaining to crack size on probability of failure Case N F P F PF Original 5,990 6.00E 03 7.71E 05 Improved 4,548 4.50E 03 6.73E 05 Table 7 5. Effects of improvements pertaining to fracture toughness on probability of failure Case N F P F PF Original 5,990 6.00E 03 7.71E 05 Improved 3,391 3. 40E 03 5.81E 05 Table 7 6. Effects of combined improvements pertaining to crack size and fracture toughness on probability of failure Case N F P F PF Original 5,990 6.00E 03 7.71E 05 Improved 2,492 2.50E 03 4.98E 05 PAGE 184 184 Figure 77 Variation of probabil ity of failure with percentage reduction of crack size. Figure 78 Variation of probability of survival with crack size. PAGE 185 185 Figure 79 Variation of probability of survival with crack size for a Kth of 4.5 0.54 MP a m. PAGE 186 186 CHAPTER 8 DISCUSSION AND CONCLUSIO NS 8.1 Discussion The fatigue damage process in ceramic rolling elements is very different from metal balls. Silicon nitride balls, used in hybrid ball bearings, are susceptible to failure from fatigue spalls emanating from preexisting c cracks that can grow under RCF. In contrast, RCF in metal bearings is manifested as a flaking off of metallic particles from the surface of raceways and/or rolling elements. This process commences as a subsurface initiated crack below the surface and propagates to the sur face, eventually forming a pit or spall (Sadeghi et al., 2009). Bearing raceway fatigue life estimation models are still mostly based on the seminal probabilistic life model by Lundberg and Palmgren (LP) proposed in 1945, and includes Hertzian subsurface s tress computations and empirical material variables obtained from extensive experimental testing. Metals are weaker in shear than tension. In contrast to metallic materials, Si3N4 material is weaker in tension than compression or shear. The LP methods ca nnot account for this difference in material behavior. Furthermore, failure in ceramic hybrid bearings is a consequence of preexisting c cracks on the ball surface growing to a spall. Fatigue life estimation for ceramic balls therefore requires a fracture mechanics approach that evaluates the critical flaw size in balls, since failure typically occurs when the effective SIF at the crack tip equals the mixed mode fracture toughness of the silicon nitride ball material. This approach is largely deterministic and radically different from the empirical lifing approaches used for metal bearing raceways and elements. This thesis presents a comprehensive numerical and experimental life prediction methodology for ceramic balls subject to RCF using fracture mechani cs principles. This PAGE 187 187 approach breaks new ground in the lifing of ceramic elements subject to RCF. The results are of immediate interest to the hybrid ball bearing and aerospace industry because of applications to military, commercial and space propulsion sy stems. The procedures and results generated have both scientific and technological relevance. Contributions that are of scientific and intrinsic interest to fracture and contact mechanics are described below: 1 3D nonplanar geometry of surface cracks nucleating from the tensile stress field at the contact patch periphery, during collision of ceramic balls, was determined from a stressstate analysis and iterative crack growth via FEA. The stress state analysis utilized a closedform subsurface stress field of contacting spheres with friction to compute the 3D trajectory of the max principle stress plot to generate the crack front. 2 First derivation of tensile stress field at the circular contact periphery including friction. 3 Determination of effective fra cture toughness, under mixedmode loading conditions with KI, KII and KIII active, for silicon nitride. 4 Computation of mixedmode SIFs for 3D nonplanar surface cracks subject to RCF. 5 Establishing the importance of surface traction towards dramatically inc reasing crack driving force during RCF. 6 Derivation of the tensile stress field at the elliptical contact periphery and recognizing that it controls the maximum SIF during RCF. Furthermore, the maximum SIF experienced by a crack due to an elliptic rolling c ontact is that generated by a circular contact enclosed by the ellipse, with the same max Hertz stress. 7 Computation of the critical flaw size, as a function of contact stress, traction condition and material fracture properties. Contributions that are of engineering and industrial interest for the design of ceramic rolling elements are described below: 1 Derivation of minimum collision velocity required to initiate a surface crack between identical brittle spheres. Inclusion of friction at the contact decreas es the velocity required to initiate a surface crack. 2 Development of a technique that can be used to predict the shape of any surface flaw from a known contact interaction. The cone and partial cone crack shapes PAGE 188 188 simulated are reminiscent of the range of surface crack geometries observed in manufactured silicon nitride balls. 3 Rendering and meshing of nonaxisymmetric nonplanar 3D surface cracks in FEA and analysis of the part subject to Hertzian rolling contact. Parametric studies were made possible by user c reated code to automate the creation of input files, loading and boundary conditions, and the post processing of the displaced states. Appendix B serves as an educational tool for future modelers. 4 Simulating the growth of 3D cracks via iterative remeshing methods in FEA. 5 Development of accurate curve fits for computing SIFs for surface flaws subject to rolling contact. 6 Development of probabilistic analysis to determine the survivability of the ball subject to a range of RCF conditions. 7 Integration of resul ts from numerical and experimental investigations into a design methodology for silicon nitride (ceramic) balls subject to RCF. 8.2 C onclusions The range of presented analyses combined in this document results in a complete treatment of surface flaws in brittle materials subject to RCF From the determination of flaw shape from nucleation, to the effects of rolling contact fatigue it sees in servic e, to examining if failure is imm inent, the possible analyses on flaw shapes under rolling contact fatigue is treated. The main contributions of this work are developed analysis tools for crack shape determination, the analysis of fatigue test experiments for measuring a Kth, and transferring the effect of each parameter on SIFs to the engineers who will design hybrid bearings to tolerate their service applications 8.3 Future Work There remain a few issues that would improve this field. Firstly, the establishment of a Kth and Kc for all materials under a range of mixed mode behavior is requisite. As previously discussed, many have tried and the result is a myriad of empirical equations that are often material and experiment specific. Secondly, an PAGE 189 189 examination of the failure process of the crack as it grows to spallation is necessary as it has been visited by W ang (2000) and Kida (2002) but have lead to some disagreement as to the roles of wear between the crack faces and the affect of the lubricant, as well as, how the crack turns branches and eventually spalls to failure. Furthermore, reliable, nondestructive detection of surface flaws in a nonconducting opaque material for flaw morphology (including crack depth) has yet to found and could find immediate application if automated to find and categorize flaws. PAGE 190 190 APPENDIX A THE NORMAL INDENTATION OF BRITTLE SPHE RES In chapter 1, it became clear that the crack nucleates and grows from the point in the trailing contact periphery of the oblique interaction, as that is where stress is at a maximum. With this knowledge, the contact of spheres can be reflected upon to develop predictive equations concerning crack nucleation. There have been a few studies that have normally loaded silicon nitride spheres until cracking was detected by acoustic emission (Ichikawa et al. 1995a, Ichikawa et al. 1995b, and Ohgushi et al. 1996). Interestingly, they have found a range of the maximum Hertz contact pressure to induce cracking to be in the range of 1418 GPa. With this number it can be approximated that the velocity to impart collisions for normal interactions. It is noted that the maximum principal stress from a normal spherical contact on the surface contact periphery is max12 3rr ov p (A1) Which, for a Poissons ratio of 0.28, the normal contact equation yields, 0 max6.8 p (A2) By substitut ing 16 GPa as the max pressure (the median of those pressures that have been observed to cause cracking) a maximum periphery stress of 2.3 GPa is calculated. To find the normal velocity to induce this stress (without the influence of tracti on) take: 4 1 5 2 5 3 45 34 24 3z omV E p R (A3) PAGE 191 191 (Johnson, 1987) and substitute eqn. A 2 for op to get: 2531 222 *515 50c zRm V E (A4) Substituting, for a normal collision of silicon nitride bodies with a density of 2,600 kg/m3 and radius of 6.35 mm (half inch diameter), a velocity of only16 cm/s is obtained, which would need to be verified by experiment. It is not surprising then that surface cracks are generated during low velocity collisions generated during the ball lapping operation. The impact velocity required to cause subsurface yielding during collision between a hard steel ball and a medium hard steel ball with an yield strength of 1 GPa is 14 cm/sec (Johnson, 1987), and for balls made of brittle material exhibiting minimal yielding such as Si3N4 the impact energy is absorbed in the fracture surface generation. It is also noted that the maximum periphery stress may be derived for sliding friction as well for use in maximum contact periphery stress criteria. First t he limits of the stress equations are taken as z and y approach 0 and x approaches 1 and then use the resulting components in the stress matrix to find the equation that represents the maximum root of the eigenvalue problem. The result is 1 max12 1,0,0 1 324fovv p (A5) which is expectedly a function of the contact pressure, the Poissons ratio, and the friction coefficient. Comparing Eq. ( A6) with Eq. ( A2) it is noted that the po required to induce cracking is reduced with a nonzero friction coeff icient. Another approximation can be made from these results concerning crack nucleation. Taking the periphery stress to initiate cracking, 2.3 GPa, and the fracture PAGE 192 192 toughness of 6 MPa Ka an initial flaw size of 0 these larger cracks could nucleate. How the radius of the ball affects the radius of crack can also be derived in a similar manner. For example, contact patch radius is 1 3 *3 4cPR a E (A6) Rewriting *34 3cEa P R (A7) And the max pressure is 22 3 2c o substitutingo cEa P pp aR (A8) For cone cracking to start at the periphery max max12 6.8 3oov pp (A9) Substituting this and rearranging max *13.6cR a E (A10) So if the contact patch size is roughly the same as the crack radius than it scales linearly with ball radius ( R ) assuming that the two interacting balls have the same radius. NOTE: P=total load an d 1 22 12 1211 vv E EE and 1 1211 R RR PAGE 193 193 Before determining the critical flaw size, the size of the cracks can be discussed practically. As discussed above, the fracture toughness of Si3N4 is low (about 6 MPa (Piotrowski and OBrien), and the resultant velocity to induce cone cracking of normally colliding balls is low. The velocity to induce cracking, zV has been derived (Levesque and Arakere) for normal ball collisions as: 2531 222 *515 50c zRm V E (A4) and determines that 16 cm/sec is a sufficient velocity to crack normally interacting 6.38mm (0.5 in) diameter balls. Substituting for the ball mass, 34 3 mVR in Eq. (A4) yields: 2115 3 222 *515 253c zRV E (A11) Equation ( A11) shows that the velocity to induce cracking is proportional to R3. The effect that the radius of the ball has on the size of the crack can also be calculated since contact patch radius is given by: 1 *3 3 *4 3 43c cEa PR aP ER (A12) The max Her tzian contact pressure for a ball on ball interaction is 22 3 2c o substitutingo cEa P pp aR (A13) From prior work (Levesque and Arakere, 2008), cracking will occur when a maximum stress value is reached and for Hertzian contacts the max stress value in the PAGE 194 194 half space is located just outside the contact region. For cone cracking to start at the periphery, the value of the maximum principal stress that must be exceeded is determined by eqn. A 3: max max12 6.8 3oov pp (A3) Substituting (A 13) and rearranging max *3.4cR a E (A15) Therefore, the contact patch size, and therefore the crack radius, scales linearly with ball radius ( R ), assuming that the two interacting balls have the same radius. PAGE 195 195 APPENDIX B BUILDING A CCRACK MOD EL First a note about scale The following model was created in um. This is so the smallest elements (once meshed) are above ABAQUS criteria to avoid flagging an error. The material properties are scaled to conform to this unit system. A crack discussion Multiple geometries can be modeled currently with the exception of self intersecting flaws and flaws that approach a surface at an angle. First Create Part C crack Figure B 1. Sketch module for creating a revolved part. Revolved Shape Create datum plane for next step PAGE 196 196 Figure B 2. Created revolved part. Cut from top to limit extension Figure B 3. Cutting the part from above PAGE 197 197 Figure B 4. The final crack shape. Building Gen III crack Use experimentally obtained cross section image and import vi a View>>>Image/Movie Options. Than trace with spline tool (in proper location for radius on face) and proceed as for Gen I. Create a block to be cracked A block of dimensions about 4a *8a*5a will be cut by crack. Retain edges. Note: cut the block with the crack. PAGE 198 198 Figure. B 5. Merge/cut instances window. Figure B 6. Merging a crack with a block in the Assembly module. PAGE 199 199 Half symmetry constructed first It should be noted that the images presented are of a half symmetry model. However, the creation of a f ull symmetry model is done ONLY by creating two half symmetry models and merging them (as the partitioning is not robust enough to perform on a full symmetry model right away). Partition time N ow partition the model so that it may be feasibly meshed in the future. S wept Hex elements are needed in the contour region so these cells need to have six sides and no more than 90 degrees of a circle. These partitions are done ON THE PART since this part will be used to create another part and these partitions are to be maintained. Inner and outer circles The correct dimensions of these circles were found through trial and error. In general, the tube can be about .3a of the crack and the inner tube a fourth or less of that but do not try to be so small that meshing within this tube will not run. PAGE 200 200 Figure B 7. Partitioning the surface. Partition the cell Sweep these circles along the crack tip to create tubular cells. Figure B 8. Sweeping to partition the cell. One more cell partition to go With the proper fores ight it is known that the modeler needs to have access to the crack faces. The following partitions are for this purpose. PAGE 201 201 Figure B 9. Partitioning to get access to the crack face. Set Declaration There is still quite a bit of work to do as current progress only has a geometry with no seam, crack, contact interaction, or way to get SIFs. For this manipulation, the following cells are to be sets. THESE SETS ARE CREATED ON THE PART SO THEY TRANSFER WHEN AN ORPHAN MESH IS CREATED PAGE 202 202 Figure B 10. Blown apart view of the different cells. In the interaction module, add a seam In order for the part to have crack faces and allow crack opening a seem must be inserted. Figure B 11. Creating the crack seam. Define the crack Accurate crack tip displacements and SIFs are desired so a crack needs to be declared. This will give greater control in the mesh module. Use the following settings: Q vectors, .777, 0, .63 Midside node parameter =0.25 PAGE 203 203 Degenerate Element control at Crack Tip/Line: Collapsed element s ide, single node This will give you the Q vectors that wi ll need to be changed (since the model is nonplanar), collapsed, quarter point elements. Figure B 12. Defining the crack. Loading and BCs Create a User defined load on the top surface. A symmetry condition if you have a half symmetry model. Submodel BCs on the remaining four/five faces. (More on this later). PAGE 204 204 Figure B 13. Loads and boundary conditions defined on the assembly. *DLOAD FORTRAN based code has been created which defines normal and traction pressures based on the global coordinate system. PAGE 205 205 Figure B 14. Defining the user defined load. Meshing A Short Course on How to Crash ABAQUS Small Problem. the brown areas, which can now be observed in the mesh module, are unmeshable. To solv e this we will turn to virtual topology. PAGE 206 206 Figure B 15. The cracked part in the mesh module. Virtual topology Virtual topology is a tool that is used to merge small faces and edges so this partitioned geometry can be feasibly meshed. Basically, ABAQUS will spline a couple of surfaces together instead of having two separate surfaces. This is important since a hexahedral mesh requires (among other things) six sides in the cell. PAGE 207 207 Figure B 16. Virtual topology on crack integration cells. BETTER!! Reassign Element Types Pink can be Tet meshed. This has been assigned via the element type selection for these (two) cells. Green can be hex structured and were easily allowed after the virtual topology was completed. Yellow will be swept wedges (though this as signment is converted by the crack definition to be collapsed quarter point hex element. PAGE 208 208 Figure B 17. Assembly after virtual topology is applied. Seeding the part Notice Biasing toward the crack, a fine mesh on the crack faces (for contact) and a given number of element tubes throughout the structure mesh. See B 18. PAGE 209 209 Figure B 18. Seeding the part for meshing. More on Seeding Note that all lines parallel to the crack tip receive equal seeding to ensure rings of elements for contour integration. Seedi ng is an iterative process where the seeds must be tuned such that where two lines meet, ABAQUS may allow a reasonable number of elements. It may be better for you to fix the number of seeds in the contour integral region so that you have a specific numbe r of contour integrals that can be taken and that the mesh on the inner most and outer most sections of the tube do not have improper aspect ratios for error in this region (remember: elements should be desirably square and anything with a .1 aspect ratio of element sides is not going to behave well). PAGE 210 210 Figure B 19. All edges that exist near the crack front have the same number of seeds and nodes because this is the contour region and planes of nodes are needed. Mesh it Take a break. this is almost halfwa y. As stated above, the seeding and meshing are currently (v. 6.71) an iterative process for this complex geometry and this particular geometry (because it is both complex and closed) will require much time. In terms of procedures to be completed this i s around the halfway mark but in terms of effort exerted this may be where the user will spend the majority of their time. PAGE 211 211 Figure B 19. Image of meshed crack. Create an orphan mesh You need to move away from an eight sided, modifiable object to a mul ti sided unmodifiable orphan mesh in order to: properly assign q vectors, create sets for the top and bottom crack faces (individually), and properly assign contact. Proceedure: Create a job>>Write inp then File>>Import Model>>select inp Or in the mesh m odule>>create mesh part. T he former method is used here PAGE 212 212 Figure B 20. Orphan mesh of crack. Reassign the Q vectors Each crack tip node must have an assigned qvector for the contour integral. This can be done individually but is best repeated with the aid of a macro. PAGE 213 213 Figure B 21. Assigning the Q vectors for each node on the crack tip. Assign Contact Assigning contact in a closed geometry is difficult to begin with. In addition to this, it must be defined between crack faces, which share a line of nodes. On top of all this it must be accurate enough to be verified through a displacement matching algorithm (~10^ Contact Settings Direct solver. PAGE 214 214 NLGEOM on. Initial step size dropped to .1 Small Sliding. Direct normal contact, with frictionless s liding, and geometric properties assigned. Node to surface. Contact Schematic Figure B 22. Meshed crack face illustration of how crack contact is defined. Please refer to the above figure which illustrates the only way you will successfully have the cra ck faces defined in contact and be accurate. (Note: self contact will not work because of the issues related to calculating a vector normal to the surface at the crack tip and surface to surface contact will not work because the crack tip nodes will be sh ared among the two surfaces which results in errors upon submission.) Post Processing Scripts were used to: Immediately plot results in Excel. PAGE 215 215 Write all three ABAQUS SIFs to a text file with corresponding nodal distance along the crack front. Compute SIFs via Displacement Correlation. Clean up the folders. Please view individual scripts for details. Running the model 1 Write the DLOAD code for proper load orientation. 2 Submit the global model for analysis. 3 Submit the submodel for analysis. (The submodel will know to reference the global modal only by the fact that it has a boundary condition stating that the surfaces in the material will have a global model and that the models attributes, MB3 on the model>>Edit Attributes, reference the Job name where the glob al model is subjected to the load.) The global model should share the same coordinates in space (relative to the global coordinate system of each model) AND should have the same user subroutine referenced so that incompatibilities between the user defined load and the submodel displaced boundary condition dont try to tear apart certain elements. 4 Run the script to either generate Ks (File>>Run Script) or extract the necessary data. 5 Save and close. 6 To post process these SIFs, look for the text file in the working directory that has the same job namethese are formatted to be imported into MATLAB for other processing. (Note that the SIFs that are calculated are likely in units of TPa need to be multiplied by 1000 to convert to MPa ion. Note: the entire procedure above could be automated for parametric study. For the desired purposes of this author each parametric study was submitted by MATLAB through the command prompt, where predecided job names were decided for each run to lat er distinguish results. (Job names are important in automation as, not only will the global and submodels need to be able to identify each other, but the post processing script must also access the proper submodel and generate a relevant name of the SIF o utput.) PAGE 216 216 LIST OF REFERENCES Aliabadi, M. H. and Rooke, D. P. 1991. Numerical Fracture Mechanics Southampton UK, Computational Mechanics Publications. Allen, M.S., Camberos, J.A., Comparison of Uncertainty Propagation / Response Surface Techniques for T wo Aeroelastic Systems, 50th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference Palm Springs, CA, May 2009 Andersson, M., 1996. Stress distribution and crack initiation for an elastic contact including friction, International Journal of Solids Structures, 33 3673 3696. Anstis, G.R., Chantikul, P., Lawn, B.R., Marshall, D.B., 1981. A critical evaluation of indentation techniques for measuring fracture toughness: I, direct crack mechanisms, J. Am. Ceram. Soc., 64, 533538. Aue rbach, F., 1891. Absolute Hardness, Ann. Phys. Chem ., 43 61 100. Awaji, H., Sato, S., 1978. Combined mode fracture toughness measurement by the disk test, J of Engg Mat. And Tech. 100 175 182. Ayatollahi, M. R. and Aliha, M. R. M., 2005. Cracked Brazili an disc specimen subjected to mode II deformation. Engineering Fracture Mechanics 72 493 503. Banks Sills, L., Hershkovitz, I., Wawrzynek, P.A., Eliasi, R., Ingraffea, A.R., 2005. Methods for calculating stress intensity factors in anisotropic materials: Part I z = 0 is a symmetric plane, Eng. Fract. Mech. 72, 15, 2328 2358. Bank Sills, L., Travitzky, N., Ashkenazi, D., Eliasi, R., 1999. A methodology for measuring interface fracture properties of composite materials, Int. J. of Fract., 99, 3, 143161. Bogdanski, S., Trajer, M., 2005, A dimensionless multi size finite element model of a rolling contact fatigue crack, Wear 258, 1265 1272. Burrier, H. I. 1996. Optimising the structure and properties of silicon nitride for rolling contact bearing perfor mance. Tribology Transactions 39 (2), 276285. Cano, J.E., Hernndez, L.D., Moral, S., (1996). Importance sampling algorithms for the propagation of probabilities in belief networks, Inte rnational Journal of Approximate Reasoning, 15 1 77 92. Carniero, F., Barcellos, A., 1953. Concrete Tensile Strength. Union of testing and research laboratories for Materials and Structures, 13. Chen, S.Y., Farris, T.N., Chandrasekar, S. 1995. Contac t mechanics of Hertzian cone cracking, International Journal of Solids Structures, 32 329 340. PAGE 217 217 Chen, X.M., Jiao, G.Q., Cui, Z.Y., 1986. Application of CombinedMode Fracture Criteria to Surface Crack Problems, Engineering Fracture Mechanics 24 1,127 14 4. Desault Systmes, 2007. ABAQUS Theory Manual, Dessault Systmes, Providence, RI. Erdogan, F., Sih, G.C., 1963. Crack Extension in Plates Under Plane Loading and Transverse Shear, J. Basic Eng. Trans. A.S.M.E., Ser. D ., 85, 519527. Evans, A. G., 1983. Progress in Nitrogen Ceramics edited by P. L. Riley, Martinus Nijhoff Publishers, The Netherlands. 595. FAC Viana, "SURROGATES Toolbox Users Guide," http://fchegury.googlepages.com, Accessed October 20, 2009. Frank, F.C., Lawn, B.R., 1967. On the theor y of Hertzian Fracture, Proc. Royal Soc. London Ser. A 299, 291306. Fracture Analysis Consultants, 2005. www.fracanalysis.com Ithica, NY. Frank, F.C., Lawn, B.R., 1967. On the theory of Hertzian Fracture, Proc. Royal Soc. London Ser. A, 299, 291306. Freiman, S.W., Gonzalez, A.C., Mecholsky, J.J., 1979. Mixed Mode Fracture in Soda Lime Glass, J. Am. Ceram. Soc., 62 3, 206208. Fowell, R. J., 1995. Suggested method for determining mode I fracture toughness using Cracked Chevron Notched Brazilian Disc (CCNBD) specimens. International Journal of Rock Mechanics and Mining Science & Geomechanics Abstracts 32 1, 5764. Fujiwara, T., Yoshioka, T. Kitahara, T., Koizumi, S., Takebayashi, H. and Taka, T. 1989. Study on load rating property of silicon nitride for rolling bearing material. Journal of JSLE International Edition 10, 81 86. Galbato, A.T. Cundill, R.T. Harris, T.A. 1992. Fatigue life of silicon nitride balls, Lubrication Engineering, 48 11, 886. Gerstle, W. H. 1986. Finite and boundary element modeling of crack propagation in two and threedimensions using interactive computer graphics. Dissertation Abstracts International Part B: Science and Engineering. 47 2, 241. G riffith, A.A., 1920. The Phenomena of Rupture and Flow in Solids, Philosophical Transactions, Series A 221 163 198. Gu, P., Nordlund, P., Asaro, R.J., Uang, C.M., 2005, Crack face contact, frictional sliding and mesh design flexibility, Comm. Numer. M eth. Engng., 21 209217. PAGE 218 218 Hamilton, G. M., Goodman, L. E., 1966. The stress field created by a circular sliding contact, Journal of Applied Mechanics, 33 371. Hamilton, G.M., 1983. Explicit equations for the stresses beneath a spherical sliding contact, Proc. Inst. Mech. Eng ., 217, 2, 281. Hadfiled, M. and Stolarski, T.A. 1995a. The effect of the test machine on the failure mode in lubricated rolling contact of silicon nitride. Tribology International 28 (6), 377382. Hadfiled, M. and Stolarski, T.A. 1995b Observations of lubricated rolling contact fatigue on silicon nitride rods. Ceramics International 21 (2), 125 130. Hadfield, M., Stolarski, T.A., Cundhill, R.T., Horton, 1993a. Failure modes of ceramics in rolling contact, S. Proc. Royal Soc. London Ser A, 443, 607621. Hadfield, M., Stolarski, T., Cundill, R.T., Horton, S., 1993b. Failure modes of ceramic elements with ring crack defects, Tribology International 26, 157164. Hadfield, M., Stolarski, T.A., Cundill, R.T. and Horton, S. 1993c. Failure modes of pre cracked ceramic elements under rollingcontact. Wear 169(1), 6975. Hertz, H., 1896. Hertzs Miscellaneous Papers McMillan, London. Hills, D.A., Nowell, D., Sackfield, A.,1992. Mechanics of Elastic Contacts, (Butterworth Heinemann, Oxford), 2 03208. Hussain, M. A., Pu, S L., Underwood, J., 1974. Strain energy release rate for a crack under mixed mode I and II loading. Fracture Analysis A.S.T.M. STP 560 2 28. Ichikawa, M., Takamatsu, T., Shindou, N., Okabe, N. and Abe, Y., 1995. Ring crack in itiation load of silicon nitride bearing balls, JSME Int. J., Ser. A, Mechanical and Materials Engineering 38 2, 226230. Ichikawa, M., Takamatsu, T., Matsuo, T., Okabe, N. and Abe, Y., 1995. Intra ball and inter ball variability of ring crack initiatio n load of silicon nitride bearing balls, JSME Int. J., Ser. A, Mechanical and Materials Engineering, 38 2, 231235. Jahanmir, S.,1994. Friction and wear of ceramics Marcel Dekker, Inc., New York, Basel, Hong Kong. Jayatilaka, A.S., Fracture of Engineerin g Brittle Materials 90 107, Applied Science Publishers Ltd., 1979. Johnson, K.L., 1987. Contact Mechanics Cambridge Press, Cambridge, 355361. Khandelwal, P., Majumdar, B.S., Rosenfield, A.R., 1995. Mixed mode high temperature toughness of silicon nitrid e, J. Material Science, 30 395 398. PAGE 219 219 Klemm, H. 2002. Corrosion of silicon nitride materials in a gas turbine environment. Journal of the European Ceramic Society 22, 14 15, 2735 2740. Kocer, C., Collins, R. E., 1998. Angle of Hertzian cone cracks, Journa l of the American Ceramic Society, 81, 7, 1736 1742. Lawn, B.R., 1994. Indentation of ceramics with spheres: a century after Hertz, Journal of the American Ceramic Society, 81 1977 1794. Lawn, B., Fracture of Brittle Solids. Cambridge University Press, C ambridge, UK (1993). Lawn, B.R., 1967, Partial cone crack formation in a brittle material loaded with a sliding spherical indenter, Proceedings of the Royal Society 299, 307. Levesque, G., Arakere N.K., 2008. An investigation of partial cone cracks in silicon nitride balls International Journal of Solids and Structures 45 6301 6315. Love, A.E.H., 1927. Treatise on the Mathematical Theory of Elasticity Dover Publications, New York. Lucek, J. W. 1990. Rolling wear of silicon nitride bearing materials. ASME paper No. 90GT 165. Lucek, J. W. and Cowley, P. E. 1978. Investigation of the use of ceramic material in aircraft engine bearings. Dept. of Navy, Code AIF 52032A Washington. L undberg, G.,Palmgren, A., 1952, Dynamic Capacity of Roller Bearings, Acta Polytechnica, Mechanical Engineering Series Royal Swedish Academy of Engineering Sciences, Vol. 2, No. 4. Lundberg, G., and Palmgren, A., 1947, Dynamic Capacity of Rolling Bearings, Acta Polytechnica, Mechanical Engineering Series Royal Swedish Academy of Engineering Sciences, Vol. 1, No. 3, pp 152. Mathworks, Inc., 2005. MATLAB Function Reference Manual, Version 7.0.4.365 (R14) Mathworks, Inc., Cambridge, MA. Mackerle, J., 2001. F inite element and boundary element analysis on indentation problems: a biography (19972000), Finite Elements in Analysis and Design, 37 811819. Mackerle, J., 2002. Ceramics and ceramic matrix composites: finite element and boundary element analyses: a b iography (1998 2000), Finite Elements in Analysis and Design, 38, 567577. Marshall, D. B., 1984. Mechanisms of Failure from Surface Flaws in Mixed Mode Loading, J.Am. Ceram. Soc. ,67, 2, 110 116. PAGE 220 220 Marshall, D.B., Lawn, B.R., Mecholsky, J.J., 1980. Effect of Residual Contact stresses on Mirror/Flaw Size Relation, J. Am. Ceram. Soc., 63, 5 6, 358360. Maw, N., Barber, J.R., Fawcett, J.N., 1976. The oblique impact of elastic spheres, Wear 38, 101114. Maw, N., Barber, J.R., Fawcett, J.N., 1981. The role of el astic tangential compliance in oblique impact, Journal of Lubrication Technology ASME Transaction, 103, 7480. Mecholsky, J. 2008. Private Communication, Department of Materials Science and Engineering, Univesity of Florida. Mi, Y., Aliabadi, M.H., 1995. An automatic procedure for mixed mode crack growth analysis, Commun. Numer. Meth. Engng., 11, (2), 167 177 Mindlin, R.D., Deresiewicz, H., 1953. Elastic spheres in contact under varying oblique forces Journal of Applied Mechanics 75, 327344. Miner, J. R., Dell, J., Galbato, A., and Ragen, M. A., 1996. F117PW 100 hybrid bearing ceramic technology insertion, ASME Journal of Engineering for Gas Turbines and Power 118, 434442. Mitchell, N.B., The Indirect Tension For Concrete, Materials Research and Standards, ASTM October 1961, 780788. Morrison, F. R. McCool, J. I. and Yonushonis, T. M. 1984. The loadlife relationship for M50 steel bearings with silicon nitride ceramic balls. J. of ASLE, Lubrication Engineering, 40 153159. Myers, R. H. and Montgomery, D. C. (2002), Response Surface Methodology, 2nd Ed., Wiley, New York, NY Newman, J.C. and Raju, I.S., 1981, An empirical stress intensity factor equation for the surface crack, Engineering Fracture Mec hanics, 15, pp. 185. Ogden, William, 2008, Private Communication. Pratt & Whitney, Co. Ohgushi, K., Ichikawa, M., 1996. Fracture mechanics study of ring crack initiation in ceramics by sphere indentation, JSME International Journal, Series A: Mechanics and Materials in Engineering 39 4, 489495. Parker R. and Zaretsky, E. V. 1975. Fatigue life of highspeed ball bearings with silicon nitride balls. Trans. of the ASME, Journal of Lubrication Technology (July), 350 357. Paris, P. and Erdogan, F. 1963. A critical analysis of crack propagation laws. ASME J our. Basic Eng. 85 528534. PAGE 221 221 Petrovic, J. J., Mendiratta, M. G., 1977. Correction of MixedMode Fracture from Controlled Surface Flaws in Hot Pressed Si3N4, J. Am. Ceram. Soc., 60, 9 10, 463. Petrovic, J.J., Mendiratta, M.G., 1976. Mixed mode fracture from controlled surface flaws in hot pressed silicon nitride, J Amer Ceram Soc 59, 3 4, 163167. Piotrowski, A. E., OBrien, M. J., 2006. A novel test method to measure th e fracture toughness of ceramic balls used in bearings Fatigue and Fracture in Engineering of Materials and Structures 29 558572 Qian, J. and Fatemi, A., 1996. Mixed mode fatigue crack growth: a literature survey. Engineering Fracture Mechanics 55, 6, 1 10. Queipo, N.V., Haftka, R.T., Shyy, W., Goel, T., Vaidyanathan, R. Tucker, P.K. (2005). Surrogate based Analysis and Optimization, Progress in Aerospace Sciences 41 1 28. Rajaram, H., S. Socrate and D. M. Parks, 2000, Application of domain integr al methods using tetrahedral elements to the determination of stress intensity factors, Engineering Fracture Mechanics 66 5, 455 482. Rice, J.R. 1968, A Path Independent Integral and the Approximate Analysis of Strain Concentration by Notches and Cracks J. of Appl. Mech., 35 379386. Richard, H.A., Fulland, M., and Sander, M. 2005. Theoretical crack path prediction, Fatigue and Fracture in Engineering of Materials and Structures, 28 3 12. Sackfield, A., Hills, D.A., 1983. A note on the Hertz contact pr oblem: a correlation of standard formulae, Journal of Strain Analysis for Engineering Design 18 2, 195 197. Sadeghi, F., Jalalahmadi, B., Slack, T.S., Raje, N., Arakere, N.K., (2009). A review of rolling contact fatigue, Journal of Tribology 131, (Accepted for publication). Sih, G.C., 1974. Strain Energy Density Factor Applied to Mixed Mode Crack Problems, Int. J. Fract., 10 3, 305 331. Shao, M., Badler, N. (1996), Spherical Sampling by Archimedes Theorem, MS CIS 9602, Department of Computer and Inform ation Science, University of Pennsylvania, PA Smith, S.M., Scattergood, R.O., 1992. Crack shape effects for indentation fracture toughness measurements, J. Am. Ceram. Soc. 75 305315. Tanaka, K., 1996. Fatigue Propagation from a Crack Inclined to the Cyc lic Tensile Axis, Eng. Fract. Mech. 55 6, 969. PAGE 222 222 Tanimoto, K., Kajihara, K., and Yanai, K., 2000. Hybrid ceramic ball bearings for turbochargers, SAE Paper 2000 01 1339, 114. Tong, J., Wong, K.Y., Lupton, C., 2007. Determination of interfacial fracture toughness of bone cement interface using sandwich Brazilian disks, Engineering Fracture Mechanics 74 12,1904 1916. Wang L. Snidle RW. Gu L. 2000, Rolling contact silicon nitride bearing technology: A review of recent research. Wear 246 159173. Wang, Y., 2000, Failure modes of silicon nitride rolling elements with ring cracks, PhD Thesis, Bournemouth University, UK. Wang Y., Hadfield, M., 2000. The influence of ring crack location on the rolling contact fatigue failure of lubricated silicon nitride: experimental studies, Wear 243, 157166. Warrier, S.G., Jarmon, D.C., and Chin, H.A. 2000. Finite element analysis of the critical flaw size in hybrid silicon nitride bearing ball. ASME Paper: 2000GT 65. Zhao, P., Hadfield, M., Wang, Y., Vieillard, C., 2006. Subsurface propagation of partial ring cracks under rolling contact Part 1. Experimental Studies, Wear 261, 382389. Zhou, J., Wang, Y. and Xia, Y., 2006. ModeI fracture toughness measurement of PMMA with the Brazilian disk test, J. of Materials Science 41, 17, 57785781. PAGE 223 223 BIOGRAPHICAL SKETCH George Arthur Levesque IV was born and raised in Dartmouth, M assachusetts After completing Dartmouth High School in 2001, he attended Bridgewater State College in Bridgewater, M assachusetts and earned a Bachelor of Science degree in p hysics (with h onors), m athematics, and c riminal j ustice in 2005. Fall 2005, he attended the University of Florida direct PhD program in the Department of Mechanical and Aerospace Engineering as an Alumni Fellow under Dr. Nagaraj Arakere. 