UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 1 VARIABLE AMPLITUDE FATIGUE ANALYSIS USING SURROGATE MODELS AND EXACT XFEM REANALYSIS By MATTHEW JON PAIS A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREME NTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2011 PAGE 2 2 201 1 Matthew Jon Pais PAGE 3 3 To my mother and father who have always suppo rted me in all that I have done PAGE 4 4 ACKNOWLEDGMENTS I would first and foremost like to thank my advisor Dr. Nam Ho Kim for his valuable insights and patience with me as we faced our challenges. He has prepared me to solve real world problems and he deserves many thanks. Dr. Timothy Davis and Nuri Yeralan from the Computer Science Department at the Universit y of Florida for their assistance in the creation and implementation of the exact XFEM reanalysis algorithm for quasi static crack growth and optimization. I would also like to thank the members of the Multidisciplinary Design and Optimization group at the University of Florida for their invaluable comments during my rehearsals. Extra thanks go to Dr. Alexandra Copp and Dr. Felipe Viana for our collaboration duri ng my time as part of the group. Dr. Benjamin Smarslok at the Air Force Research Laboratory at Wright Patterson Air Force Base for providing acceleration and velocity data from aircraft flights which were used to create data for the variable amplitude fatigue analysis. Richard Pippy for the original creation of the wing box finite element model whic h was subsequently modified for use in the analysis of an airplane wing subjected to a variable amplitude fatigue analysis. Dr. Nagaraj Arakere, Dr. Youping Chen, Dr. Timothy Davis, Dr. Raphael Haftka, and Dr. Bhavani Sankar for serving on my committee for their thoughtful insights as I worked toward this dissertation. Extra thanks go to Dr. Haftka for comments which lead to the development of the variable step size kriging algorithm. Dr. Jorg Peters from the Computer Science Department at the University of Florida for his comments which lead to a conference paper discussing enrichment functions for weak discontinuities independent of the finite element mesh. PAGE 5 5 TABLE OF CONTENTS page ACKNOWLEDGMENTS ................................ ................................ ................................ 4 LIST OF TABLES ................................ ................................ ................................ ........... 8 LIST OF FIGURES ................................ ................................ ................................ ...... 10 LIST OF ABBREVIATIONS ................................ ................................ .......................... 15 ABSTRACT ................................ ................................ ................................ .................. 16 CHAPTER 1 INTRODUCTION ................................ ................................ ................................ ... 18 Motivation ................................ ................................ ................................ .............. 18 Fatigue Crack Growth ................................ ................................ ............................ 20 Computational Fracture Mechanics ................................ ................................ ....... 26 Scope ................................ ................................ ................................ .................... 28 Outline ................................ ................................ ................................ ................... 30 2 THE LEVEL SET METHOD ................................ ................................ ................... 32 Level Set Method for Closed Sections ................................ ................................ ... 32 Level Set Method for Open Sections ................................ ................................ ..... 34 Summary ................................ ................................ ................................ ............... 36 3 THE EXTENDED FINIT E ELEMENT METHOD ................................ ..................... 37 General Form of the Extended Finite Element Method ................................ .......... 37 Enrichment Functions ................................ ................................ ............................ 40 Crack Enrichment Functions ................................ ................................ ........... 41 Inclusion Enrichment Functions ................................ ................................ ....... 45 Void Enrichment Function ................................ ................................ ............... 49 Integration of Element with Discontinuity ................................ ............................... 50 XFEM Software ................................ ................................ ................................ ..... 52 Abaq us ................................ ................................ ................................ .......... 52 MATLAB XFEM Code (MXFEM) ................................ ................................ .... 54 Other XFEM Implementations ................................ ................................ ......... 57 Summary ................................ ................................ ................................ ............... 58 4 CRACK GROWTH MODEL ................................ ................................ ................... 60 Stress Intensity Factor Evaluation ................................ ................................ ......... 60 Crack Growth Direction ................................ ................................ .......................... 64 PAGE 6 6 Crack Growth Magnitude ................................ ................................ ....................... 66 Finite Crack Growth Increment for Constant Amplitude Loading ..................... 67 Direct Fatigue Crack Growth ................................ ................................ ........... 67 Classical Paris m odel ................................ ................................ ................ 68 Modified Paris m odel ................................ ................................ ................ 69 Equivalent Stress Intensity Factor ................................ ................................ ... 72 Integration of Fatigue Crack Growth Model ................................ ..................... 73 Example p roblem ................................ ................................ ...................... 77 Cycle Counting Methods ................................ ................................ ........................ 79 Summary ................................ ................................ ................................ ............... 83 5 SURROGATE MODELS FOR HIGHER ORDER INTEGRATION OF FATIGUE CRACK GROWTH MODELS ................................ ................................ ................. 85 Integration of Ordinary Differential Equation ................................ .......................... 85 Surrogate Models ................................ ................................ ................................ .. 87 Kriging ................................ ................................ ................................ ................... 88 Kriging for Integration of Fatigue Crack Growth Law ................................ .............. 90 Constant Integration Step Size ................................ ................................ ........ 90 Variable Integration Step Size ................................ ................................ ......... 92 Example P roblems ................................ ................................ ................................ 95 Center Crack in an Infinite Plate Under Tension ................................ .............. 95 Edge Crack in a Finite Plate Under Tension ................................ .................... 99 Inclined Center Crack in a Finite Plate Under Tension ................................ .. 102 Summary ................................ ................................ ................................ ............. 105 6 AN EXACT EXT ENDED FINITE ELEMENT METHOD REANALYSIS ALGORITHM ................................ ................................ ................................ ....... 106 Reanalysis Methods ................................ ................................ ............................ 106 Exact Reanalysis by Modification of an Existing Chol esky Factorization .............. 107 Example Problems ................................ ................................ .............................. 115 Reanalysis of an Edge Crack in a Finite Plate ................................ ............... 115 Optimization for Finding Crack Initiation in a Plate with a Hole ...................... 117 Summary ................................ ................................ ................................ ............. 119 7 VARIABLE AMPLIT UDE MULTI AXIAL FATIGUE ANALYSIS OF WING PANEL 121 Scaling of Normalized Flight Data from the Air Force Research Laboratory ........ 121 Finite Element Wing Model ................................ ................................ .................. 124 Wing Pressure Distribution ................................ ................................ ............ 125 Finite Element Stress Analysis of Airplane Wing ................................ ........... 133 Surrogate Models for Lift Coefficient, Stress Components, and Surface Area of Pressure Distribution ................................ ................................ ..... 133 Analysis of Cracked Plate Subjected to Variable Amplitude Multi axial Load History ................................ ................................ ................................ ....... 135 Summary ................................ ................................ ................................ ............. 142 PAGE 7 7 8 CONCLUSIONS AND FUTURE WORK ................................ .............................. 145 Conclusions ................................ ................................ ................................ ......... 145 Future Work ................................ ................................ ................................ ......... 147 Future Development of MXFEM ................................ ................................ .... 147 Possible Application for the Exact XFEM Reanalysis Algorithm .................... 147 Possible Applications for Surrogate Integration ................................ ............. 149 APPENDIX A MXFEM BENCHMARK PROBLEMS ................................ ................................ ... 151 XFEM Enrichments ................................ ................................ .............................. 151 Center Crack in a Finite Plate ................................ ................................ ........ 15 1 Edge Crack in a Finite Plate ................................ ................................ .......... 153 Inclined Edge Crack in a Finite Plate ................................ ............................. 154 Ha rd Inclusion in a Finite Plate ................................ ................................ ...... 155 Void in an Infinite Plate ................................ ................................ .................. 158 Other Benchmarks ................................ ................................ ............................... 160 Angle of Crack Initiation from Optimization ................................ .................... 160 Crack Growth in Presence of an Inclusion ................................ ..................... 162 Fatigue Crack Growth ................................ ................................ ................... 163 B AUXILIARY DISPLACEMENT AND STRESS STATES ................................ ....... 165 C CONVERSION OF NORMALIZED FLIGHT DATA TO SCALED FLIGHT DATA 166 D STRESS HISTORIES FOR FLIGHT DATA FROM AFRL ................................ .... 170 LIST OF REFERENCES ................................ ................................ ............................ 177 BIOGRAPHICAL SKETCH ................................ ................................ ......................... 194 PAGE 8 8 LIST OF TABLES Table page 5 1. Accuracy of kriging interpolation integration for fixed increment a for a center crack in an infinite plate of Al 2024 with R = 0.5 ................................ ................. 97 5 2. Accuracy of kriging extrapolation integration for fixed increment N for a center crack in an infinite plate of Al 2024 with R = 0.5 ................................ ................. 98 5 3. Effect of material and initial crack size on estimated cycles to failure for the variable step size algorithm for a center crack in an infinite plate ...................... 99 5 4. Accuracy of kriging interpolation integration for fixed increment a for an edge crack in a finite plate of Al 2024 with R = 0.5 ................................ ................... 100 5 5. Accuracy of kriging extrapolation integration for fixed increment N for an edge crack in a finite plate of Al 2024 with R = 0.5 ................................ ................... 101 5 6. Effect of material and initial crack size on estimated cycles to failure for the variable step size algorithm for an edge crack in a finite plate ......................... 102 5 7. Comparison of Euler and variable step size predictions for the coordinates ( x final y final ) of the right crac k tip marked with a dotted circle in Figure 5 6 and cycle number N fail corresponding to a final crack size of 0.6 m from an initial crack size of 0.187 m for four different materials ................................ ............. 103 7 1. Maximum and minimum values of normalized and scaled flight data with conversion relationship for Flight ID 152 ................................ .......................... 124 7 2. The polynomial coefficients for the chord wise pressure distribution .................. 128 7 3. Area under the normalized pressure distribution for some angles of attack ........ 132 7 4. The stress components for various a ngles of attack used in the lift and drag pressure distributions ................................ ................................ ...................... 133 7 5. Summary of the AFRL data, the number of cycles identified, and the number of function evaluations for each Flight ID ................................ ............................. 138 A 1. Convergence of normalized stress intensity factors for a full and half model for a center crack in a finite plate ................................ ................................ .......... 152 A 2. Convergence of normalized stress intensity factors for an edge crack in a finite plate ................................ ................................ ................................ ................ 154 A 3. Convergence of normalized stress intensity factors for an inclined edge crack in a finite plate ................................ ................................ ................................ 155 PAGE 9 9 A 4. Comparison of the theoretical and MXFEM value for maximum energy release angle as a function of average element size ................................ .................... 161 C 1. Normalized and scaled flight data with conversion for Flight ID 2 ...................... 166 C 2. Normalized and scaled flight data with conversion for Flight ID 8 ...................... 166 C 3. Normalized and scaled flight data with conversion for Flight ID 10 .................... 166 C 4. Normalized and scaled flight data with conversion for Flight ID 12 .................... 166 C 5. Normalized and scaled flight data with conversion for Flight ID 16 .................... 166 C 6. Normalized and scaled flight data with conversion for Flight ID 18 .................... 167 C 7. Normalized and scaled flight data with conversion for Flight ID 22 .................... 167 C 8. Normalized and scaled f light data with conversion for Flight ID 26 .................... 167 C 9. Normalized and scaled flight data with conversion for Flight ID 138 .................. 167 C 10. Normalized and scaled flight data with conversion for Flight ID 140 ................ 167 C 11. Normalized and scaled flight data with conversion for Flight ID 146 ................ 168 C 12. Normalized and scaled flight data with conversion for Flight ID 148 ................ 168 C 13. Normalized and scaled flight data with conversion for Flight ID 150 ................ 168 C 14. Normalized and scaled flight data with conversion for Flight ID 152 ................ 168 C 15. Normalized and scaled flight data with conversio n for Flight ID 154 ................ 168 C 16. Normalized and scaled flight data with conversion for Flight ID 156 ................ 169 C 17. Normalized and sca led flight data with conversion for Flight ID 159 ................ 169 C 18. Normalized and scaled flight data wi th conversion for Flight ID 161 ................ 169 C 19. Normalized and scaled flight data wi th conversion for Flight ID 163 ................ 169 PAGE 10 10 LIST OF FIGURES Figure page 1 1. Representative S N curve for a material subjected to cyclic loading .................... 18 1 2. Representation of the Mode I, Mode II, and Mode III opening mechanisms ......... 22 1 3. An example of constant amplitude loading with characteristic parameters such as the maximum and minimum applied stresses and the stress ratio R ............. 23 1 4. An example of an overload followed by an underload and an underload followed by an overload ................................ ................................ ..................... 25 1 5. Examples of possible traction separation models that can be used as part of a cohesive zone model where Y is the yield stress and c is the critical separation ................................ ................................ ................................ ......... 27 2 2. Example of the signed distance functions for an open section ............................. 35 3 1. The nodes enriched w ith the Heaviside and crack tip enrichment functions ......... 42 3 2. One dimensional bi material bar problem ................................ ............................ 48 3 3. Comparison of the various inclusion enrichment functions and locations of the enriched degree(s) of freedom for a one dimensional bar ................................ 48 3 4. Elements containing a discontinuity and the continuous subdomains for integration ................................ ................................ ................................ ......... 51 3 6. Example problem to show plots generated by MXFEM ................................ ........ 56 3 7. Example of the level set functions output fro m MXFEM ................................ ....... 57 3 8. Example of the stress contours output from MXFEM ................................ ........... 57 4 1. The elements within a specified radius for path indep endence ............................ 63 4 2. Orientations of the critical a nd maximum normal stress planes ............................ 66 4 3. The use of the stress ratio correction f actor M R to condense the various stress ratios R into a single curve corresponding to R = 0 ................................ ............ 70 4 4. The plastic zone and crack sizes which are used in the modified Paris model to account for the effect of crack tip plasticity during load interactions ................... 71 4 5. Plate with a hole subjected to tension to examine the effects of mesh density, crack growth increment, and number of elapse d cycles upon the convergence of crack path ................................ ................................ ................ 78 PAGE 11 11 4 6. Convergence of crack path as a function of the mesh density, crack growth increment, and number of elapsed cycles ................................ ......................... 78 4 7. Close up view of the crack paths for a around the hole ................................ ...... 79 4 8. The strategy used to convert the bi axial stress data into cyclic stress data w hich was suitable for use in a fatigue crack growth model ............................... 82 5 1. Available discrete points available which can be fit using a surrogate mode l ....... 86 5 2. Surrogate model fit to either a complete or partial crack growth history ............... 88 5 3. Comparison of an exact function and the corresponding kriging surrogate model for an arbitrary set of five data points ................................ ...................... 90 5 4. Integration of fatigue crack growth model from N i to N i +1 using either Euler or midpoint approximation ................................ ................................ ..................... 92 5 5. Errors used to adjust the adaptive integration step size algorithm ........................ 94 5 6. An inclined center crack in a square finite plate under uniaxial tension .............. 103 5 7. Comparison of cycle number and integration step size or percent error for surrogate for R = 0.1 and R = 0.5 for aluminum 2024 ................................ ....... 104 6 1. Finite element mesh with enriched nodes ................................ .......................... 109 6 2. Edge crack in a finite plate geometry used to test quasi static crack growth with the reanalysis algorithm ................................ ................................ ........... 115 6 3. Normalized computational time for stiffness matrix assembly as well as factorization and solving of system of linear equations as a function of number of iterations ................................ ................................ ......................... 117 6 4. Normalized computational time for stiffness matrix assembly as well as factorization and solving of system of linear equations as a function of number of mesh density ................................ ................................ .................. 117 6 5. Edge crack in a finite plate geometry used to test optimizatio n with the reanalysis algorithm ................................ ................................ ........................ 118 6 6. Normalized computational time for stiffness matrix assembly as well as factoriz ation and solving of system of linear equations as a function of number of iterations ................................ ................................ ......................... 119 7 1. Coordinate system ( x y z ) and directional names for an airplane ........................ 122 7 2. Example of a wing box with span s chord length c which is a function of y and normalized x coordinate u ................................ ................................ ................ 125 PAGE 12 12 7 3. Wing box model in the coordinate s ystem that was used for this analysis .......... 125 7 4. Comparison of elliptical loading model for various wing taper ratios ................... 126 7 5. E xamples of the changing pressure distribution over an airfoil as a function of the angle of attack ................................ ................................ ........................... 128 7 6. Normalized chord wise coordinates u with maximum pressure coordinate u 0 pressure di stribution P ................................ ................................ ..................... 128 7 7. Comparison of the changes in the maximum pressure location u 0 and subsequently the resulting pressure profile for some angles of attack .......... 129 7 8. The total normalized distribution along the airplane wing for = 5 .................... 129 7 9. Comparison between the angle of attac k and the lift coefficient ......................... 131 7 10. The direction of the uniform pressure distribution defined by the drag ............. 132 7 11. The stre ss histories for the two dimensional stress components for Flight ID 152 ................................ ................................ ................................ .................. 135 7 12. Equivalent stress and identified cycles from equivalent stress for Flight ID 152 137 7 13. The boundary conditions and applied bi axial stresses for the analysis of a panel on an airplane wing box subjected to variable amplitude fatigue ............ 13 9 7 14. The initial and final crack geometries from the stresses identified through the normalized flight data from the AFRL ................................ .............................. 140 7 15. Close up view of the final crack geometry from the stresses identified through the normalized flight data from the AFRL ................................ ........................ 141 7 16. Comparison of the crack length as a function of the cycle number for the 19 flights from the normalized AFRL dat a ................................ ............................. 141 7 17. The Mode I and Mode II stress intensity factor ranges for the 37,007 cycles ... 142 8 1. Possible crack propagation pa ths between dissimilar materials ......................... 148 A 1. Representative geometry for a center crack in a finite plate ............................... 152 A 2. Representative geometry for an edge crack in a finite plate ............................... 153 A 3. Representative geometry for an inclin ed edge crack in a finite plate .................. 155 A 4. Representative geometry for a hard inclusion in a finite plate ............................ 156 A 5. The xx contours for a hard inclusion ................................ ................................ 157 PAGE 13 13 A 6. The xy contours for a hard inclusion ................................ ................................ 157 A 7. The yy contours for a hard inclusion ................................ ................................ 157 A 8. Representative geometry for a vo id in an infinite plate ................................ ....... 158 A 9. The xx contours for a void in an infinite plate ................................ .................... 159 A 10. The xy contours for a void in an infinite plate ................................ .................. 159 A 11. The yy contours for a void in an infinite plate ................................ .................. 159 A 12. Representative geometry for a crack initia ting at an angle for a plate with a hole ................................ ................................ ................................ ................. 160 A 13. Comparison of crack paths to those predicted by Bordas for a hard and soft inclusion ................................ ................................ ................................ .......... 163 A 14. Comparison between Paris model and XFEM simulation for fatigue crack growth for a center crack in an infinite plate ................................ .................... 164 D 1. Stress history for Flight ID 2 ................................ ................................ .............. 170 D 2. Stress history for Flight ID 8 ................................ ................................ .............. 170 D 3. Stress history for Flight ID 10 ................................ ................................ ............ 170 D 4. Stress history for Flight ID 12 ................................ ................................ ............ 171 D 5. Stress history for Flight ID 16 ................................ ................................ ............ 171 D 6. Stress history for Flight ID 18 ................................ ................................ ............ 171 D 7. Stress history for Flight ID 22 ................................ ................................ ............ 172 D 8. Stress history for Flight ID 26 ................................ ................................ ............ 172 D 9. Stress history for Flight ID 138 ................................ ................................ .......... 172 D 10. Stress history for Flight ID 140 ................................ ................................ ........ 173 D 11. Stress history for Flight ID 146 ................................ ................................ ........ 173 D 12. Stress history for Flight ID 148 ................................ ................................ ........ 173 D 13. Stress history for Flight ID 150 ................................ ................................ ........ 174 D 14. Stress history for Flight ID 152 ................................ ................................ ........ 174 PAGE 14 14 D 15. Stress history for Flight ID 154 ................................ ................................ ........ 174 D 16. Stress history for Flight ID 156 ................................ ................................ ........ 175 D 17. Stress hist ory for Flight ID 159 ................................ ................................ ........ 175 D 18. Stress history for Flight ID 161 ................................ ................................ ........ 175 D 19. Stress history for Flight ID 163 ................................ ................................ ........ 176 PAGE 15 15 LIST OF ABBREVIATION S AFRL Air Force Research Laboratory AMD Approximate Minimum Degree Algorith m CTOD Crack t ip o pening d isplacement DOF Degrees of f reedom EPFM Elastic p lastic f racture m echanics FE Function evaluations FEM Finite e lement m ethod GUI Graphical u ser i nterface LEFM Linear e lastic f racture m echanics LSM Level s et m ethod MXFEM MATLAB XF EM Code PUFEM Partition of u nity f inite e lement m ethod XFEM Extended f inite e lement m ethod PAGE 16 16 Abstract 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 VARIABLE AMPLITUDE FATIGUE ANALYSIS USING SURROGATE MODELS AND EXACT XFEM REANALYSIS By Matthew Jon Pais December 2011 Chair: Nam Ho Kim Major: Mechanical Engineering Fatigue crack growth occurs as the result of repeated cyclic loa ding well below the stress levels which typically would cause failure. The number of cycles to failure for high cycle fatigue is commonly of the order of 10 4 10 8 cycles to failure. Fatigue is characterized by a differential equation which gives the crack g rowth rate as a function of material properties and the stress intensity factor. Analytical relationships for the stress intensity factor are limited to simple geometries. A numerical method is commonly used to find the stress intensity factor for a given geometry under certain loading. The use of kriging to assist higher order approximations is introduced. Here stress intensity factor data is fit using a surrogate This surrogate is used to extrapolate for the purpose of integration which enables larger s tep sizes to be taken without a loss in accuracy for the solution of the governing differential equation. Furthermore, i t was observed that for the extended finite element method a small portion of the global stiffness matrix is changed as a result of crac k growth. It is possible to use this small portion to save on both the assembly and solution of the resulting system of linear equations. This results in savings in both the assembly and factorization of the stiffness PAGE 17 17 matrix for repeated simulations reduci ng the computational cost associated with numerical fatigue crack growth. The use of the XFEM reanalysis algorithm allow for the analysis of non proportional mixed mode variable amplitude loading upon an airplane wing to be considered. An airplane wing box was analyzed using Abaqus The Abaqus stress solution was used in coordination with airplane flight data provided by the Air Force Research Laboratory. This stress history is converted into a cyclic loading history through the use of the rainflow counti ng method. The resulting analysis is one where approximately 30,000 cycles elapse. Due to the non proportional loading, each cycle must be modeled independently as the direction of crack propa gation will be cycle dependent creating a solution which would b e very expensive with existing techniques. PAGE 18 18 CHAPTER 1 INTRODUCTION Motivation The n ucleation and propagation of cracks in engineered structures is an important consideration for the design of a structure. In particular fatigue fracture, caused by repeated cyclic loading well below the yield stress of a material can cause sudden, catastro phic failure. The relationship between the applied stress and the number of cycles to failure is typically given by a material specific N curve as shown in Figure 1 1 For fatigue failure once a small crack has formed, the cycli c loadings cause the material ahead of the crack to slowly fail and the crack grows. Initially a crack grows very slowly, maybe on the order of nanometers at a given cycle. Over time the crack growth accelerates. Once the crack reaches a critical length a c a large amount of crack growth occurs rapidly and without warning. There are many incidents caused by the growth of fatigue cracks, many of which resulted in the loss of human lives. Figure 1 1 Re presentative S N curve for a material subjected to cyclic loading. ( Note that there is a point where if the applied stress is sufficiently small that the number of cycles to failure is infinite ) [ 1 ] was ente red into service. In January and May of 1954 two of the Comets disintegrated during flights between New York and London. The failure was caused by fatigue cracks which Number of Cycles to Failure Applied S tress PAGE 19 19 initiated near the front of the cabin on the roof. Over time the crack grew until it rea ched a window, causing sudden catastrophic failure. In 1957 the 7 th President of the Philippines died [ 2 ] along with 24 others when fatigue caused a drive shaft to break, subsequently causing power failure aboard the airplane. In 1968 a helicopter [ 3 ] crashed in Compton, California due to fatigue failure of the blade spindle. Twenty one lives were lost. In 1980 the Alexan der L. Keilland [ 4 ] oil platform capsized killing 123 people. The main cause of the failure was determined to be a poor weld, reducing the fatigue strength of the structure. In 1985 Japan Flight 123 [ 5 ] from Tokyo to Osaka crashed in the deadliest plane crash in history. The incorrect repair of a prior incident where the tail of the airplane impacted the runway in 1978 eventually lead to f atigue failure. There were 4 survivors of the 524 people on the airplane. Fatigue cracks slowly grew until causing sudden rupture 7 years later. In 1988 an Aloha Airlines [ 6 ] flight between the Hawaiian islands of Hilo and Honolulu suffered extensive damage after an explosive decompression caused by th e combined effects of a fatigue crack and corrosion. The plane safely landed in Maui with 94 survivors, 65 injuries and 1 death. In this case, the hours instead of 75,0 00 design hours) as well as operating in a corrosive environment caused by exposure to high levels of humidity and salt. In 1989 a United Airlines flight [ 7 ] from Denver, Colorado to Chicago, Illinois crashed due to a maintenance crew failing to find a crack in a fan disk within the engine 112 deaths occurred. In 1992 a Boeing 747 [ 8 ] crashed in to the Bijlmermeer neighborhood in Amsterdam, Netherlands, killing 43. T he plane crashed when fatigue failure did not PAGE 20 20 allow for the engine to cleanly separate from the wing as designed, leading to the accident. In 1998 an InterCityExpress train [ 9 ] crashed in Eschede, Germany caused by fatigue failure of the train wheels. Of the 287 passengers aboard, there were 101 deaths and 88 injuries. It is the deadlie st train disaster in German history. In 2002 China Airlines Flight 611 [ 10 ] broke apart during a flight killing all 225 people aboard the airplane. Similar to Japan Airlines Flight 123, an incorrect repair procedure allowed fatigue cracks to grow eventually causing failure. In 2005 metal light [ 11 ] from Fort Lauderdale, Florida to Bimini, Bahamas to crash in Miami Beach, FL after metal fatigue broke off the right wing. There were 20 casualties. In 2007 a Missouri Air National Guard F 15C Eagle [ 12 ] crashed due to a structural part not meeting specifications, leading to a series of fatigue cracks to develop and propagate. As recently as July 2009 a Southwest Airlines flight [ 13 ] from Nashville, Tennessee to Baltimore, Maryland had to make an emergency are still undergo ing. The series of accidents through history including those in the 1990s and 2000s show that there is still much work to be done to prevent fatigue failures from taking human lives. Fatigue Crack Growth Fatigue crack growth models are empirical models whi ch are generally created by performing one or a series of experiments and fitting the resulting data to a function of the form [ 14 ] 1 1 PAGE 21 21 where d a / d N is the crack growth rate and is the stress intensity factor ratio, which is a driving mechanism to crack growth. The stress intensity factor range can be calculated as 1 2 where and are the stress intensity factors at the maximum and minimum applied stress during a given loading cycle. The stress intensity factor is used to describe the state of stress at the tip of a crack and depends upon, crack location, crack size, distribution and magnitude of loading, and specimen geometry. There are three modes of fracture [ 15 ] as shown in Figure 1 2 On the left is Mode I or opening mode. In the middle is the Mode II or in plane shear mode. Finally, on the right is the Mode III or out of plane shear mode. In two dimensions only Modes I and II is considered, while Mode III is also considered for three dimensional problems. When multiple modes occur simultaneously, the underlying problem is commonly referred to as being mi xed mode. Mixed mode fractures can be separated into the underlying Mode I, II, and III components through the use of methods such as the domain form of the contour integral [ 16 19 ] For some simplified cases analytical expressions for the stress intensity factors are available [ 20 ] but for a more general case a num erical method such as the finite element method may be used to find the stress intensity factor. Specifically, the use of a method such as the crack tip opening displacement (CTOD) [ 15 ] J integral [ 17 ] or the domain form of the contour integral [ 18 19 ] may be used to extract the mixe d mode stress intensity factors from a finite element solution. Note that Eq. 1 1 does not provide the crack size directly. Rather, it provides the rate of cr ack growth at a given range of PAGE 22 22 stress intensity factor, which also depends on the current crack size. Thus, a numerical integration method should be employed to predict the crack growth by integrating the given fatigue model. In addition, the crack in a co mplex geometry may not grow in a single direction. The direction of future growth is often considered to be governed by the maximum circumferential stress criterion in the finite element framework as closed form solutions for the direction of crack growth are available [ 21 ] Figure 1 2 Representation of the Mode I, Mode II, and Mode III opening mechanisms. There are two main opinions on how to attempt to approximate the solution of the governing differential equation in Eq. 1 1 for constant amplitude loading such as that shown in Figure 1 3 : controlling cycles or crack increment. The first model assumes th at a fixed number of cycles is chosen prior to starting the simulations [ 22 ] Thus, is variable, being initially small and increasing with the iteration number. At each crack gr owth iteration, is evaluated from an analytical formula or numerical simulation which is then used to approximate the solution of the ordinary differential equation governing fatigue crack growth. Since the expression of is unknown as a function of crack size, the differential equation cannot be exactly integrated using any numerical integration technique. Instead, there are only values of at the current and all past PAGE 23 23 simula tion iterations available for use during the integration. An explicit numerical integration method such as the forward Euler method [ 23 ] can be used to a pproximate the current growth increment The use of a higher order integration method such as the midpoint [ 23 ] or Runge Kutta [ 23 ] methods require additional approximations for the crack sizes which require the corresponding function evaluations. Fur thermore, the required function evaluations may have no physical meaning, especially considering that crack growth does not occur at all instances within a loading cycle [ 24 ] The forward Euler method requires small time steps to be accurate; otherwise c rack growth will be under predicted. As the crack increment is based on it is essential for the mesh to be sufficiently refined around the crack tip such that has converged and an accurate representa tion of the localized state of stress around the crack tip. Figure 1 3 An example of constant amplitude loading with characteristic parameters such as the maximum and minimum applied stresses and the stress ratio R The other solution procedure assumes a fixed size of crack increment [ 25 ] at each simulation iteration. Therefore, the number of elapsed cycles is initially large and decreases with the iteration number. Here, the challenge is accurately back calculating the number of elapsed cycles from the fixed crack growth data. The forward Euler max N min R = min / max PAGE 24 24 method is commonly applied to the back calculation of the number of elapsed cycles [ 26 ] but as with the approach of fixing there must be care taken in selecting a sufficiently small value of for the forward Euler method to maintain accuracy. A higher order approximation could be applied her e as well, but discrete values of are not available for all data points required for the evaluation of the slope of the a N curve for these higher order methods. These values require additional function evaluations. For a consta nt amplitude loading simulation involving mixed mode crack growth, it is also imperative that the choice of fixed or be made with care. The path that the crack may take is also influenced by the choic e of or as having too large of a value of either can result in deviation from the crack growth path that would be predicted with the use of smaller growth increments. Once this deviation occurs the pa th of future crack growth is definitely affected, but this deviation is not clear to a user unless a convergence study with respect to crack path is performed. In the literature [ 22 25 ] it is clear that the chosen fixed crack increments influence the crack path under mixed mode loading. For the case of variable amplitude fatigue where the applied loading does not have a constant amp litude, It may not be possible to choose an incremental or and have an accurate solution to the governing differential equation given by Eq. 1 1 As the magnitudes of the applied stress are changing with each cycle, the crack path is changing on a cycle by cycle basis. Attempting to model this behavior considering multiple loading cycles at a time would assume a constan t crack growth direction for many loading cycles and fail to capture the changes in the crack path based on the changes in the applied stress profile. PAGE 25 25 During each loading cycle, plasticity forms around the crack tip due to the stress concentration that the crack introduces into a material. There are many models for the size of this plastic zone [ 14 15 27 ] The main challenge associated with modeling variable amplitude fatigue i s with the interaction of overloads and underloads [ 14 24 27 30 ] and the associated changes in the plastic zone size. An overload can be considered to be a peak load which is an uncha racteristically large magnitude, while an underload is a valley load with an uncharacteristically small magnitude as shown in Figure 1 4 Figure 1 4 An example of an overl oad followed by an underload and an underload followed by an overload. The effect of an overload is an increase in the crack tip plasticity compared to that which would come from the constant amplitude loading portion of Figure 1 4 This increase in plasticity acts to retard future crack growth. Similarly, an underload results in a reduced amount of plasticity compared to the constant amplitude loading portion of the cyclic loading curve. The two cases shown in Figure 1 4 an overload followed by an max N min Overload Underload PAGE 26 26 underload and an underload followed by an overload result in different crack growth paths and must be considered for an accurate analysis of variable amplitude fatigue. Attempts to capture the behavior associate d with the interaction between overload and underloads has been explored through the use of empirical models such as crack growth models [ 14 27 31 34 ] the small time scale model [ 24 ] the state space model [ 29 30 ] i n addition to elastic plastic fracture mechanics (EPFM) analysis in the FEM [ 35 36 ] as well as the extended finite element method (XFEM) [ 37 38 ] Computational Fracture Mechanics With the ever increasing speed of computers numerical methods such a s the finite element method (FEM) are able to model problems with increasing fidelity and increasing complexity, but finite element modeling of fatigue crack growth is still a challenging computational fracture mechanics problem. As the finite element mesh must conform to the geometry the mesh around the crack tip must be recreated [ 25 ] whenever growth occurs. Even if a concept such as crack blocks [ 39 ] is used where a small region around a crack tip is remeshed at each iteration, the computational demands of the remeshing can contribute signif icantly to the simulations especially if a large number of iterations of crack growth are to be modeled. In addition to the loc alized mesh reconstruction around the crack tip, care needs to be taken to ensure accuracy in the crack path, displacement, and stress intensity factors with FEM. It is possible to model crack growth based on the cohesive zone mo del in the classical FEM fr amework [ 40 ] The cohesive zone model assumes that crack growth can be represented as a traction separation model as shown in Figure 1 5 Special cohes ive zone elements can then be used which allow for crack growth without the need to remesh. However, some knowledge about the crack path is needed a priori in PAGE 27 27 order to define which elements in the mesh will be cohesive elements instead of traditional methods T he model is based on the assumption that once the yield stress Y is reached, the material properties start to degrade until some critical separation distance c is reached. At this point the given element/interface in the finite element model can no lon ger sustain loading and fracture occurs. Figure 1 5 Examples of possible traction separation models that can be used as part of a cohesive zone model where Y is the yield stress and c is the cri tical separation. The main challenge in cohesive zone model s is to create an accurate model for the given traction separation relationship for a given material system including the shape of the traction separation model and the critical separation distance [ 40 42 ] It is often the case that the parameters of the traction separation model for a given material change according to the length scale in the model and co rresponding experimental tests In practice, it can be the case that small changes to the parameters in the traction separation model result in substantial differences in the resulting finite element solution [ 40 ] Moody [ 41 ] introduced a series of experiments which can be used to Y Y Y Y c c c c PAGE 28 28 determine the traction separation relationship. Scheider [ 42 ] explored the use of finite element simulations to determine the traction separation relationship. The shape of the softening function that will be used to degrade the material propert ies is also of much importance [ 40 ] The XFEM [ 25 ] alleviates the challenges associated with the mesh conforming to the geometry by allowing discontinuities or other localized phenomena to be represented independent of the finite element mesh. Enrichmen t functions based on linear elastic fracture mechanics (LEFM) are available which are independent of material, unlike cohesive zone models, which are material and possible length scale dependent. Additional functions, referred to as enrichment functions, a re introduced into the displacement approximation through the property of the partition of unity [ 43 ] Additional nodal degrees of enrichment functio ns as well as are used for interpolating within an element and in the calculation of stresses or stress intensity factors using a method such as the domain form of the conto ur integrals. While there is no need to worry about mesh construction or the need t o create a traction separation relationship the XFEM still requires the convergence of the crack path, displacement, and stress intensity factors for accurate crack modeling. Scope The goals and scope of this work are focused on accurate modeling of fatig u e crack growth under constant and variable amplitude loading for complex geometr ies with unknown relationships with computational efficiency and high accuracy In particular the following topics are addressed: PAGE 29 29 Increasing the al lowable step size for a given fixed increment of a or N without a loss in accuracy when approximating the solution to the governing ordinary differential equation for constant amplitude loading The influence that this increased step size has on the conv ergence of the crack path is also considered. A surrogate model is exploited to enable the use of higher order approximations to the calculation of magnitude and direction of crack growth with a single expensive function evaluation A variable step size al gorithm is also introduced to better allocate the av ailable computational resources based on the accuracy of the surrogate model when extrapolating. The fundamental formulation of the XFEM is exploited to enable the modeling of quasi static crack growth wi th reduced computational time through a n exact reanalysis algorithm. When crack growth occurs, the changes to the global stiff ness matrix are limited to a localized number of elements about the crack tip. Here, a supernodal Cholesky factorization is used t o exploit these properties. During the first XFEM simulati on, the full stiffness matrix is formed and a Cholesky factorization is performed. In all subsequent iterations this Cholesky factorization is directly modified to account for the changes to in the stiffness matrix caused by crack growth. Computational savings are realized for both the assembly and factorization of the global stiffness matrix. This reanalysis algorithm is also employed as a means to consider optimization problems in the XFEM framewor k where the location of a discontinuity is iteration dependent. Multi axial variable amplitude fatigue analysis is performed for the case where the stress components are non proportional. In this problem, traditional approaches based upon a or N fail as the crack growth direction becomes cycle dependent. As the crack path is unknown, the likelihood of having an analytical equation available for the stress PAGE 30 30 intensity factors as a function of crack size is small. Due to the variable amplitud e loading, the changing plasticity at the crack tip caused by the varying maximum and minimum applied loads affects the crack growth rate. A test problem of an airplane wing panel subjected to biaxial stresses during 19 flights based upon normalized flight data from the Air Force Research Laboratory (AFRL) is considered. A method to scale and convert the normalized data into equivalent biaxial loading cycles is given. Outline Chapter 2 introduces the level set method for tracking closed and open sections. T his method is used to track the location of discontinuities in the XFEM as they do not conform to the mesh. This includes c racks, inclusions and voids. Chapter 3 introduces the extended finite element method First the general form is considered. Then enri chment functions are introduced for cracks, inclusions and voids. A discussion of the commercial and open source implementations of XFEM is presented Chapter 4 details the domain form of the contour integral for the extraction of mixed mode stress intensi ty factors from a XFEM analysis of a cracked body. Criteria which de termine the direction of crack growth are introduced. Finally, fatigue crack growth models predicting the magnitude of crack growth are given and limitations associated with these models a re discussed. In addition, rainflow counting methods are introduced for the purposes of converting stress histories into equivalent stress cycles. Chapter 5 introduces the use of a surrogate model for increased accuracy in the integration of the ordinary d ifferential equation governing fatigue crack growth allowing larger step sizes to be cons idered without loss of accuracy for the case of constant amplitude loading. A surrogate model may also be used to enable higher order approximations to the crack growt h direction. An algorithm is also presented where the step size dynamically changes based on the PAGE 31 31 surrogate accuracy in order to minimize the number of finite element simulations needed. Chapter 6 introduces and details the use of a reanalysis algorithm to make the repeated simulations of crack growth in a quasi static environment affordable, allowing for additional simulations in a fixed amount of time. The sparse Cholesky factorization is detailed as well as the algorithms which an be used to modify an exi sting sparse Cholesky factorization. The proposed reanalysis algorithm leads to reduced computational cost for repeated XFEM simulations Chapter 7 details the variable amplitude fatigue analysis procedure including the calculation of stress histories from aircraft flight data, the conversion of these stress histories into cyclic loads which are used in fatigue model prediction, and the prediction of crack growth under variable amplitude loading using surrogate models and an exact XFEM reanalysis algorithm Chapter 8 summaries the research as well as suggesting possible areas of future work. PAGE 32 32 CHAPTER 2 THE LEVEL SET METHOD Level Set Method for Closed Sections The level set method was introduced by Sethian and Osher [ 44 ] as a numerical method which can be used to track the evolution of interfaces and shapes. The method is based on evolving an interface subjected to a front velocity given by the physics of the underlying problem which is being modeled. Level set methods have been used in a wide range of engineering applicatio ns in topics such as compressible [ 45 ] and incompressible [ 46 ] flow, computer vision [ 47 ] image processing [ 48 ] manufac turing [ 49 51 ] and structural optimization [ 52 53 ] While the figures given in this chapter are two dimensional, the principles can also applied to three dimensional problems [ 54 55 ] The level set method is used to discretize the domain of interest into discrete points Each of these points is assigned a signed distance value from that point to the nearest intersection with the interface denoted A continuous level set function is introduced where is a point in the domain of interest The level set function can be characterized a s a function of the do main and with a time component as 2 1 Thus, points inside the domain of interest are given negative signs, points outside of the domain of interest are given positive signs, and points on the interface have no sign as their signed distance is zero. An example of the signed distance function for a circular domain is given in Figure 2 1 From Eq. 2 1 it can be noted that at any time t the location of the interface can be found as the locations where PAGE 33 33 2 2 and is commonly referred to as the zero level set of The evolution of the level set function is usually assumed to follow the Hamilton Jacobi equation [ 44 ] where the evolution can be predicted as 2 3 where is the front velocity and is the spatia l gradient of the level set function The solution of Eq. 2 3 is usually approximated using finite differencing techniques [ 23 ] which provide sufficient solution accuracy when the time step is small. When the forward finite difference technique is considered, the derivative of with respect to time can be approximated as 2 4 where is the updated level set value, is the current level set value, is the front velocity vector, and is the elapsed time between and Equation 2 4 can be rewritten in a more convenient form in two dimensions as 2 5 where is the front velocity in the x direction and is the front veloc ity in the y direction. In Eqs. 2 4 and 2 5 the time step is limited by the Courant Friedrichs Lewy (CFL) condition [ 56 ] which ensures that the approximation to the solution of the partial differential equation converges. The CFL condition is given as 2 6 PAGE 34 34 where and are the grid spacing in the x and y directions. In practice, the level set function needs to only be defined in a narrow band [ 57 59 ] around the inter faces of interest or can be represented using the fast marching method [ 60 ] a variant of the level set method with improved computational efficiency. Figure 2 1 Example of a signed distance function for a closed domain. Level Set Method for Open Sections The version of the level set method presented in the previous section is only valid for a closed section. For an open section, the definition of the interior, exterior, and interface as defined in Eq. 2 1 no longer have a physical meaning. Stolarska [ 57 ] introduced a modified version of the level set method which allows for open sections to be tracked with the use of multiple level set functions. A n open section as shown in Figure 2 2 can be described by two level sets and The interface of interest is given as the intersection of and where 2 7 PAGE 35 35 Figure 2 2 Example of the signed distance function s for an open section. (Note: t he interface of interest is given as the region where phi is negative and psi is equal to zero ). An updating algorithm for these two coupled level set function is also given by Stolarska [ 57 ] For the case of the level set function, the update is identical to that presented in Eqs. 2 3 2 5 Two regions are defined with respect to the level set function, and which correspond to the regions which will and will not be updat ed. The level set function is updated in two dimensions a ccording to 2 8 where the crack tip displacement vector is given as and the current crack tip is given by the coordinates The sign of the updated value is chosen to correspond to the location of that node with respect to Figure 2 2 PAGE 36 36 Summary The level set method allows for a closed or open section to be tracked by defining signed distance values at discrete points in the domain of in terest. The value of the level set function at these points is then updated based on the front velocity at each point in the domain using a finite difference technique to approximate the solution to the governing partial differential equation. The CFL cond ition is used to ensure that the solution of the differential equation converges. The level set method seems to be ideal for use in a finite element environment where the nodes of the finite element mesh could be used as the fixed points in the level set a lgorithm. The finite element shape functions could be used to interpolate within an element to identify the values of the level set function if this would be of interest. PAGE 37 37 CHAPTER 3 THE EXTENDED FINITE ELEMENT METHOD General Form of the Extended Finite E lement Method The extended finite element method (XFEM) allows for discontinuities to be represented independent of the finite element mesh by exploiting the partition of unity finite element method [ 43 ] (PUFEM). In this method additional functions, commonly referred to as enrichment funct ions, can be added to the displacement approximation as long as the partition of unity is satisfied, i.e. for all where are the finite element shape functions. The XFEM uses these enrichment functions as a tool to represent a non smooth behavior of field variables, such as stress across the interface of dissimilar materials or displacement across cracks. In general, the enrichment functions introduced into the displacement ap proximation are only defined over a small number of elements relative to the total size of the domain. Additional degrees of freedom are introduced in all elements where the discontinuity is present, and depending upon the type of function chosen, possibly some neighboring elements which are known as blending elements. The additional functions used in the displacement approximation are commonly called enrichment functions and the approximation takes the form: 3 1 where are the classical finite element degrees of freedom (DOF) is the J th enrichment function at the I th node and are the enriched DOF corresponding to the J th enrichment function at the I th node. The enriched degrees of freedom introduced by Eq. 3 1 generally do not have a physical meaning and instead can be considered as a PAGE 38 38 calibration of the enrichment functions which result in the correct displacement approximation. Note that Eq. 3 1 does not satisfy the interpolation property, due to the enriched DOF instead additional calculations are required in order to ca lculate the physical displacement using Eq. 3 1 The interpolation property is important in practice in applying boundary or contact conditions Therefore, it is common practice to shift [ 61 ] the enrichment function such that 3 2 where is the value of the J th enrichment function at the I th node. As the shifted enrichment function now takes a value of zero at all nodes, the solution of the resulting system of equations satisfies and the enriched DOF can be used for additional actions such as interpolation and post processing operations Here, the shifted enrichment functions are referred to with upper case characters, and the unshifted enrichment functions are referred to with lower case characters. The shifted displacement approximation is given by 3 3 wh ere is the J th shifted enrichment function at the I th node. Hereafter, and will be written as and The Bubnov Galerk in method [ 62 ] may be used to convert the displacement approximation given by Eq. 3 3 into a system of linear equations of form 3 4 PAGE 39 39 where is the global stiffness matrix, are the nodal DOF and are the applied nodal forces. By appropriately ordering degrees of freedom, t he global stiffne ss matrix can be considered as 3 5 where is the classical finite element stiffness matrix, is the enriched finite element stiffness matrix, and is a coupling matrix between the classical and enriched stiffness comp onents. The elemental stiffness matrix for any member of may be calculated as 3 6 where is the constitutive matrix for an isotropic linear elastic material, is the matrix of classical shape function derivatives, an d is the matrix of enriched shape function derivatives. The general form of and is given by 3 7 where is the derivative of with respect to and is the derivative of with respect to In practice, is c alculated with the product rule: PAGE 40 40 3 8 Similarly and in Eq. 3 4 are given by 3 9 where and are vectors of the classical and enriched degrees of freedom and 3 10 where and are vectors of the applied forces for the classical and enriched components of the displacement approximation. The vectors and are given in terms of applied tractions and body forces as 3 11 and 3 12 Stress and strain must be calculated with the use of the enrichment functions and enriched degrees of freedom such that the effect of the discontinuity with in a particular element is considered. Therefore the strain and stress may be calculated as 3 13 and 3 14 Enrichment Functions The XFEM has been used to solve a wide range of problems involving disco ntinuities. In general, discontinuities can be described as either strong or weak. A PAGE 41 41 strong discontinuity can be considered one where both the displacement and strain are discontinuous, while a weak discontinuity has a continuous displacement but a discont inuous strain. There exist enrichment functions for a variety of problems in areas including cracks, dislocation, grain boundaries, and phase interfaces [ 63 66 ] Aquino [ 67 ] has also studied the use of proper orthogonal decomposition to incorporate experimental data in to the displacement approximation for cases with no logical choice of enrichment function. Fries [ 68 ] introduced the use of hanging nodes in the XFEM framework with respect to inclusions, cr acks, and fluid mechanics to allow for automated mesh refinement around discontinuities. Crack Enrichment Functions The modeling of cracks in the XFEM has been thoroughly explored [ 63 66 69 70 ] Belytschko [ 71 ] was t he first to study cracks in the XFEM framework based on the element free Galerkin crack enrichment of Fleming [ 72 ] Mo s [ 25 ] introduced the use of the Heaviside enrichment function to simplify the representation of the crack away from the tip. Work has been done in two [ 25 57 71 73 74 ] and three dimensions [ 54 58 60 75 ] for linear elastic [ 25 54 57 60 71 73 75 ] elastic plastic [ 37 76 ] and dynamic [ 77 81 ] fractu re. The common practice is to incorporate two enrichment functions into the XFEM displacement approximation to represent a crack. A Heaviside step function [ 25 ] is use to represent the crack away from the tip and a more complex set of functions is used to represent the crack tip asymptotic displacement field. The Heaviside step function is given as 3 15 PAGE 42 42 It can be noticed that the enrichment given by Eq. 3 15 introduces a discontinuity in displacement across the crack. For a linear elastic crack tip, four enrichment functions [ 72 ] are used to incorporate the crack tip displacement field into elements containing the crack tip: 3 16 where and are the polar coordinates in the local crack tip coordinate system the origin it at the crack tip and is parallel to the crack Note that the first enrichment function in Eq. 3 16 is discontinuous across the crack behind the tip in the element containing the crack tip, acting as the Heaviside enrichment does. Should a node be enriched by both Eqs. 3 15 and 3 16 only Eq. 3 16 is used as shown in Figure 3 1 where the Heaviside and crack tip nodes are denoted by filled circles and squares Figure 3 1 The nodes enriched with the Heaviside and crack tip enrichment functions. Elgeudj [ 37 38 ] id entified crack tip enrichment functions [ 37 ] w hich can be used to capture the elastic plastic fracture behavior of the Hutchinson R ice Rosengren (H R R) singularity [ 82 83 ] The H R R singularity is a model for confined plasticity in the fatigue PAGE 43 43 of power law hardening materials. From a Fourier analysis three basis functions were iden tified, and the basis with the lowest rank was chosen as the XFEM enrichment function. These basis functions are given by 3 17 where is the power law hardening exponent for the given material Comparison between the enrichment functions of Eq. 3 16 and Eq. 3 17 showed very similar predictions in stress intensity factor s between the elastic and elast ic plastic case for several values of These crack tip enrichment functions were then impl emented to model elastic plastic fatigue crack growth [ 38 ] for a material subjected to a combination of overload and underload conditions with a very limited number of loading cycles The goal of this work was the capture the underlying plasticity evolution caused by the interaction of overload and underloads. In benchmark problems the stress intensity factor results do not show a significantly different prediction in stress intensity factor over the traditional crack tip enrichment meth od, but require significantly more computational resources for a given analysis. The similarity between the elastic plastic and linear elastic cases is likely due to the confined plasticity about the crack tip during fatigue crack growth. This was also obs erved by Anderson in the classical FEM [ 84 ] Alternative crack tip conditions have also been explored such as bi material cracks [ 85 ] cohesive cracks [ 86 88 ] branching cracks [ 89 ] cracks under frictional contact [ 38 90 ] fretting fatigue cracks [ 91 ] interfacial cracks [ 76 85 92 93 ] cracks in orthotropic materials [ 94 ] and cracks in piezoelectric materials [ 95 ] Mousavi [ 96 ] introduced a unified framework for the enrichment of homogeneous, intersecting, and branching cracks through the use of harmonic enrichment functions. The XFEM has also been PAGE 44 44 used to study a variety of problems involving cracks including: the effect of cracks in plate s [ 97 98 ] crack detection and identification [ 99 101 ] shape optimization [ 26 ] and optimization with changing crack locati on using a reanalysis technique [ 22 ] Because the XFEM mesh does not need to conform to the domain, a method must be used to track of the location of the cracks. To this end the use of the open segment level set method introduced by Stolarska [ 57 ] and detailed in Chapter 3 is used Two level set functions are used to track the crack, the zero level set of represents the crack body, while the zero level sets of which is orthogonal to the zero level set of represents the location of the crack tip s The two enrichment functions given in Eqs. 3 15 and 3 16 can be calculated in terms of and such that 3 18 Furthermore, the polar crack tip coordinates are given as 3 19 The enriched nodes corresponding to the crack tip enrichment can also be determined through the use of th e le vel set functions defining the crack. Consider an element where the maximum and minimum values of and are given as and Then an element is enriched with the Heaviside enrichment when 3 20 and the crack tip enrichment when 3 21 PAGE 45 45 Therefore, the extended finite element and level set methods complement one another well for the tracking of the location of the cracks. The representation of cracks in three dimensions [ 54 58 60 ] follows a similar methodology. In practice the level sets are defined in only a narrow band about the crack as discussed in Chapter 2 or the fast marching method [ 54 60 ] is used. The convergence rate of XFEM with crack enrichment functions has been an area of interest [ 65 102 107 ] particularly with respect to the challenges presented by the partially enriched or blending elements caused by the crack tip enrichment. No blending issues exis t with the Heaviside function as it vanishes along all element boundaries. It was noticed by Stazi [ 104 ] that the convergence rate for the XFEM was lower than the equivalent tr aditional finite element problem. Chessa [ 108 ] identified that the partially enriched crack tip elements lead to p arasitic terms in the displacement approximation and introduced an enrichment dependent assumed strain model to increase the convergence rate Fries [ 102 ] introduced a linearly decreasing enrichmen t weight function in the blending elements to increase convergence. An area [ 103 109 110 ] instead of single element crack tip enrichment has also been shown to increase convergence. Through the use of these methods the convergence rate of cracked domains with the XFEM has become equivalent to the equival ent traditional finite element problem [ 65 102 107 ] I nclusion Enrichment Functions The modeling of material interfaces independent of the finite element mesh through the element free Galerkin [ 111 ] as well as partition of unity finite e lement method [ 55 61 63 65 66 112 115 ] has been studied. The enrichment function should PAGE 46 46 incorporate the behavior of the weak discontinuity, i e. continuous displacement, but discontinuous strain The Hadama rd condition [ 113 ] given by 3 22 where is the deformation gradient, is the outward normal material interface, and is an arbitrary vector in the plane. The Hadamard condition must be satisfied by the chosen enrichment funct ion Sukumar [ 113 ] first introduced the use of the absolute value enrichment in terms of the level set function which gives the shortest signed distance from a given point to the interface between the two ma terials. Therefore, the enrichment function takes the form: 3 23 The enrichment function is assumed to be nonzero only over the domain of support for the enriched nodes, as with the crack enrichment function. For a bi material boundary value benchmark problem the absolute value enrichment given by Eq. 3 23 led to a convergence rate which was lower than the equivalent tr aditional finite element method problem where the mesh conforms with the material interface. It was hypothesize d that the poor convergence was related to the blending elements containing a partial enrichment. In an attempt to improve the convergence rate a smoothing algorithm was introduced to reduce the effects of the blending elements whic h increased the converge nce rate, but did not equal the traditional finite element method Mo s [ 55 ] studied modeling complex microstructure geometries with the use of level set defined material interfaces and introduced a new enrichment function. The modified absolute value enrichment takes the form: PAGE 47 47 3 24 Note that the enrichment function given by Eq. 3 24 is zero at all nodes and thus, does not need to be shifted such that traditional degrees of freedom are recovered automatically If an interface corresponds to the mesh, then no nodes are enriched as the enrichment function will be zero and the problem will be equivalent to the traditional finite element problem. The same benchmark problem considered by Sukumar [ 113 ] was considered as well as a similar problem in three dimensions. In two dimensions the convergence rate was shown to equal the traditional finite element method. In three dimensions the convergence rate was slightly less than the traditional finite element method. This method is co nsidered the current state of the art for modeling inclusions with the XFEM Pais [ 112 115 ] considered an element based enrichment instead of nodal enrichment where the displacement approximation took the form 3 25 w here is a piecewise linear enrich ment function where 3 26 and are elemental degrees of freedom. Thus, the enrichment function vanishes at all nodes and takes a value of one at the interface locations. The proposed method allows for the number of elem ental degrees of freedom to be equal to the number of dimensions of the problem. The resulting system of equations needs fewer degrees of freedom than either the traditiona l or extended finite element method to represent the PAGE 48 48 same domain. It was found that the convergence rate is comparable to the absolute value enrichment given by Eq. 3 23 due to errors in the prediction of the shear stress distribution in two and three dimensions. A comparison of the enrichment functions for Eqs. 3 23 3 25 i s given for the case of a bi material bar shown in Figure 3 2 A comparison of the enrichment functions and locations of t he enriched degrees of freedom for the absolute value, modified absolute value and element based enrichment functions is given in Figure 3 3 Figure 3 2 One dimensional bi material bar problem. Figure 3 3 Comparison of the various inclusion enrichment functions and locations of the enriched degree(s) of freedom for a one dimensional bar. The problem given in Figure 3 2 can be solved using the absolute value enrichment, modified absolute value enrichment, or element based enrichment all of x L x 0 0 0 1 y Absolute value Element based Modified absolute value Enriched DOF Bar 1 P 0.25L 0.75L Bar 2 PAGE 49 49 which yield equivalent final answers. When additional elements are considered, the smoothing of the abso lute value enrichment presents a challenge not found with the modified absolute value or element based enrichment for recovering the theoretical displacement. Due to the improved convergence rate for the modified absolute value enrichment this method is th e most popular approach in the literature for modeling inclusions with XFEM. For the most part, continued work in the area of modeling inclusions in the XFEM tends to application oriented instead of focusing upon the development of enrichment functions whi ch offer improved performance. Other work on modeling inclusions in the XFEM include a unified model for the representation of arbitrary discontinues and discontinuous derivatives [ 61 ] The imposition of constraints along moving or fixed interfaces was considered by Z ilian [ 114 ] Dirichlet and Neumann boundary conditions for arbitrarily shaped interfaces were presented for a general enrichment function. Instead of Lagrange multiplier [ 116 117 ] or penalty method [ 118 ] approaches, a mixed hybrid method was introduced. Constant boundary tractions, prescribed displacement differences and prescribed interfacial displacement states were applied to a bi material pro blem. Hettich [ 119 120 ] studied the modeling and failure of the inte rface between fiber and matrix in composite materials. Void Enrichment Function Daux [ 89 ] was the first to represent voids with the XFEM Sukumar [ 113 ] later extended the void enrichment to take advantage of the use of the level set function to track the void. Unlike the other enrichment functions presented here, the void PAGE 50 50 enrichment function does not require additional DOF ; instead the displacement approximation for a domain with a hole takes the form 3 27 where takes a value of 0 inside the void and 1 anywhere else. In practice, integration is simply skipped where Additionally, nodes whose support is completely within the void are considered fixed DOF Integration of Element with Discontinuity An area where the XFEM differs from the classical finite element method is on the scheme used to perform the numerical integration of Eq. 3 6 as c hallenges arise in elements which contain a discontinuity. Standard Gauss quadrature [ 23 ] requires that the integrands are smooth, which is not t he case for an element containing a strong or weak discontinuity. The approach introduced by Mo s [ 25 ] was to divide a two dimensional element into a set of triangular subdomains, where the discontinuity was placed along the boundary of one of the subdomains. Integration would then be performed over each subdomain, resulting in a series of integrations over continuous domains. The common number of gauss points for integration in each triangular subdomain with the Heaviside enrichment is 3 while 7 are used for the crack tip enrichment functions given in Eq 3 16 [ 73 ] An example of an element completely cut by a crack as well as containing a crack tip and the associated subdomains for integ ration are shown in Figure 3 4 The creation of the subdomains is straightforward with the use of Delaunay tessellation for the nodal coordinates along with the coordinates of the zero level set of in t he parametric space PAGE 51 51 Figure 3 4 E lements containing a discontinuity and the continuous subdomains for integration. A) element cut by crack, B) element cut by crack divided into four continuous subd omains, C) element containing crack tip, D) element containing crack tip divided into 5 continuous subdomains. Elgu e dj [ 37 38 ] introduced a hybrid integration method for the case of modeling the evolution of plasticity of an elast ic plastic crack on a structured mesh. Elements which were within a given radius around the crack tip were subdivided into 16 square subdomains. Within each square subdomain, 16 gauss points were used in order to track the evolution of plasticity within th e model. As the crack evolved, all gauss points associated with previous plasticity were updated and additional elements were divided into subdomains according to crack growth. Note that 256 gauss points are used for an elastic plastic crack enriched eleme nt, while one would expect an average of 35 gauss A C B D PAGE 52 52 points for the case without plasticity. The implication is that in addition to solving a nonlinear system of equations to capture the plasticity based on the HHR singularity model, more gauss points will be used for a comparable simulation. This leads to a computational problem which is much more demanding if plasticity is to be modeled with XFEM. In three dimensions, it is possible to decompose elements cut by planar cracks into a series of tetrahedrons [ 75 ] in a similar fashion to that of two dimensional elements. Mousavi also explored integration over arb itrary polygons [ 121 ] with an application to XFEM [ 121 ] and through the use of the Duffy transformation [ 122 ] showing excellent accuracy in the presence of singularities Sukumar [ 123 ] presented a method for the integration of an arbitrary polygon based on Schwarz Christoffel conformal mapping Yamad a [ 124 ] presented a hybri d numerical quadrature scheme based on a modified Newton Cotes quadrature scheme. Park [ 125 ] introduced a mapping method for integration of discontinuous enrichments. Other methods are detailed by Belytschko [ 63 ] and Fries [ 66 ] XFEM Software Due to the r elatively short history of the XFEM commercial codes which have implemented the method are not prevalent. There are however, many attempts to incorporate the modeling of discontinuities independent of the finite element mesh by either a plug in or native support [ 66 126 129 ] As most of these implementations are works in progress there are various limitations on thei r practical use. Abaqus In 2009, the Abaqus 6.9 release [ 128 ] introduced basic XFEM functionality to the Abaqus CAE environment. The Abaqus implementation of the XFEM is somewhat PAGE 53 53 different from that which was previously presented in this chapter. The implementation is based on the phantom node method which was introduced by Hansbo [ 130 ] and subsequently modified by Song [ 131 ] and Rabczuk [ 132 ] The fundamental difference between this implementation and the original XFEM is that the discontinuity is described by superimposed element s and phantom nodes. In effect, an element is only defined in an area where an element is continuous Several elements are combined together such that the total behavior of a discontinuous element is described. The method implemented by Abaqus considers a cohesive crack model, which is only enriched with the Heaviside function given by Eq. 3 18 Note that it has been shown repeatedly that the use of only the Heaviside enric hment leads to poor accuracy of the resulting J integral calculation [ 129 ] As a result of this enrichment scheme, all enriched elements must be completely cut by the crack, as no crack tip field is considered. Some limitations with modeling crack growth within Abaqus using the XFEM are: Only the STATIC analysis procedure is allowed Only linear continuum elements are allowed with or without reduced integration No fatigue crack g rowth models are available No intersecting or branching cracks are allowed A crack may not turn more than 90 within a particular element The Abaqus 6.9: Extended Functionality [ 127 ] update allows for energy release rate and stress intensity factors to be evaluated for three dimensional cracked domain. There is currently no method available for the extraction of stress intensity factors in two dimensions. In practice, the modeling of crack growth within Abaqus with the current implementation of the XFEM is challenging. As the formulation is based on the cohesive model the same challenges exist of solving the system of equations as with PAGE 54 54 cohesiv e elements and methods such as viscous regularization must be applied to solve the system of equations. Care must be taken to choose the regularization parameters in a way that has a minimal impact on the resulting solution The Abaqus 6.10 [ 133 ] updated enables parallel processing of some XFEM simulations. The Abaqus 6.10: Extended Functionality [ 36 ] update enables low cycle fatigue predictions. MATLAB XFEM Code ( MXFEM ) Basic functionality of the XFEM was implemented by Pais [ 134 ] in MATLAB for two dimensional plane stress and plane strain probl ems This code has been given the name MATLAB XFEM Code (MXFEM) rectangular domain may be defined with a structured grid of linear square quadrilateral elements with arbitrary loading and boundary conditions. Enrichments pr ovided include the homogeneous [ 25 ] crack, inclusion [ 55 ] and void [ 89 113 ] All discontinuities are tracked using the level set method detailed in Chapter 3 and enrichment functions are calculated from level set values as detailed previously in this chapter The and level set functions track the crack, the level set function t racks the voids, and the level set function tracks the inclusions. Integration of enriched elements is done through subdivision of elements into triangular regions [ 25 73 ] A variety of plotting outputs may be requested including: level set functions, finite element mesh, deformed finite element mesh, elemental and contours of stress, and stress intensity factor history for growing cracks. A graphical user interface (GUI) is available which offers simplified functionality compared to the direct modification of the input file. The GUI writes an input file based on the values of t he GUI and then solves the problem. An example of the GUI is given in Figure 3 5 Examples of some of these PAGE 55 55 plots are given for the geometry shown in Figure 3 6 Figure 3 7 and Figure 3 8 which contains a circular inclusion below the crack and a void above the crack. The domain form of the contour integrals [ 18 19 ] is used to calculate the mixed mode stress intensity factors. The mixed mode stress intensity factors are used in the critical plane approach [ 133 ] or the maximum circumferential stress criterion [ 18 ] to give the direction of crack extension. Crack growth may be modeled using either a constant increment of growth or a fatigue crack growth model, such as Paris model [ 34 ] or a modified Paris model [ 14 27 ] which considers the effects of variable amplitude loading. Four models [ 133 135 137 ] are available for the conversion of the mixed mode stress intensity factors into an equivalent stress intensity factor range. For the case of variable amplitude loading a separate input file may be used to define the cycle d ependent loading parameters for the two dimensional state of stress. Either the forward Euler method or surrogate models as presented in Chapter 4 can be used in the integration of the fatigue model. Surrogate models are used through the use of the SURROGA TES Toolbox [ 138 ] All crack growth problems are solved using the reanalysis algorithm presented in Chapter 5. Additional functionality includes an optimization algorithm for f inding some optimum crack location and the ability to define a variable load history. The code has been prepared in an attempt to facilitate other models for crack growth direction, equivalent stress intensity factor, and fatigue crack growth models. Bench mark problems for the various enrichment functions are provided in Appendix A including: crack, inclusion, and void enrichment functions, fatigue crack growth, and optimization to identify the crack location with maximum energy release rate for a plate wit h a hole. PAGE 56 56 Figure 3 5 MXFEM GUI for automated input file creation. A B Figure 3 6 Example problem to show plots generated by MXFEM A) The geometry being considered, B) Example of t he mesh output from MXFEM where circles and squares denote the Heaviside and crack tip enriched nodes and the circles denote the inclusion enriched nodes. H/4 H/4 H/4 H/4 a r r PAGE 57 57 A B C D Figure 3 7 Example of the level set functions output from MXFEM A) B) C) D) A B C D Figure 3 8 Example of the stress contours output from MXFEM A) B) C) D) Other XFEM Implementations Bordas [ 69 ] implemented the XFEM as a object oriented library in C++. Cenaero [ 139 ] has implemented XFEM functionality into the Morfeo finite element environment. PAGE 58 58 Giner [ 129 ] implemented the XFEM in Abaqus with the traditional crack tip enrichment scheme for modeling two dimensional growth through the use of user defined elements Global Engineering and Materials, Inc. [ 126 140 ] is developing a X FEM based failure prediction tool as an Abaqus plug in which includes a small time scale fatigue model [ 24 ] for variable amplitude loading Sukumar [ 73 74 ] implemented the XFEM in Fortran specifically the finite element program Dynaflow Wyart [ 141 ] discussed the implications of the implementation of the XFEM in general commercial codes. Some a mount of functionality is also available in getfem++ [ 142 ] and openxfem++ [ 143 ] Summary The XFEM allows for strong and weak discontinuities to be represented independent of the finite element mesh by incorporating discontinuous functions into the displacement approximation through the partition of unity finite element method. Additional nodal degrees of freedom are introduced at the nodes of elements cut by a discontinuity. An assortment of enr ichment functions are available for a variety of crack tip conditions and the Heaviside step function is used to model the crack away from the tip. Elastic plastic crack tip enrichment functions are available for materials which power law hardening based o n the HHR singularity model The stress intensity factor results do not show a significantly different prediction in stress intensity factor over the traditional crack tip enrichment method, but require significantly more computational resources for a give n analysis. Inclusions are modeled using the modified absolute value enrichment, which shows convergence equivalent to that of the classical finite element method. Voids can be incorporated through the use of a step function and without additional nodal de grees of freedom. Integration of elements containing PAGE 59 59 enriched degrees of freedom, and therefore a discontinuity require special integration treatment through the use of either a subdivision or equivalent algorithm. As the discontinuities do not correspond to the finite element mesh, some other method must be used to track the discontinuities. The level set method for open sections is used to track any cracks in the domain. The level set method for closed sections is used to track inclusions and voids. The l evel set values are used as part of the definition of the enrichment functions, leading to a symbiotic relationship between the extended finite element and level set methods. A version of the XFEM for the modeling of cohesive cracks has been introduced in Abaqus 6.9, but it is not without its limitations. An implementation within MATLAB had been completed here. The XFEM is also available through the use of plug ins or some open source finite element codes. PAGE 60 60 CHAPTER 4 CRACK GROWTH MODEL Stress Intensity Factor Evaluation The most common way to extract the mixed mode stress intensity fact ors is through the use of the domain form of the interaction integrals [ 16 18 19 85 ] in two [ 25 73 89 ] or three dimensions [ 54 59 75 ] The domain form of the interaction integrals is an extension of the J integral originally introduced by Cherepanov [ 144 ] and Rice [ 17 ] In the domain form, the line inte gral specified by the J integral is converted to an area integral which is much more amenable to use with finite element simulations. As the J integral is used to calculate the energy release rate for a given crack, the interaction integrals are used to ex tract the mixed mode stress intensity factors. This method has been shown to have excellent accuracy when applied to suitable meshes for many crack conditions including homogenous [ 25 ] bi material [ 85 ] and branching [ 89 ] cracks. Other methods [ 63 64 66 ] have also been explored in the literature. For a general mixed mode homogeneous crack in two dimensions the energy release ra te can be expressed in terms of the relationship between the J integral, as 4 1 where is defined by the state of stress as 4 2 where The J integral takes the form PAGE 61 61 4 3 where is the strain energy densit y. Equation 4 3 can be rewritten in the equivalent form using the Dirac delta which is easier to implement in finite element code as 4 4 In order to calculate the mixed mode stress intensity factors, two displacement and stress states are superimposed onto one another. Auxiliary stress and displacement states are superimposed onto the stress and displacement solution from XFEM. The auxiliary stress and displacement states at the crack tip introduced by Westergaard [ 145 ] and Williams [ 146 ] for a homog enous crack and by Sukumar [ 85 ] for a bi material crack are used in the calculation of the mixe d mode stress i ntensity factors. The auxiliary stress and dis placement states for a homogene ous crack are given in Appendix B Hereafter, the XFEM states are given as and while the auxi liary states are given as and The superposition of stress states into Eq. 4 4 leads to 4 5 The expansion of the terms of Eq. 4 5 allows for the J integral to be separated to the auxiliary state XFEM state and interaction state given by 4 6 where is the i nteraction strain energy density given as PAGE 62 62 4 7 The two superimposed st ress states can be expressed using Eq. 4 1 after rearrangement as 4 8 Therefore, from Eqs. 4 6 and 4 8 the interaction state is given as 4 9 T he stress intensity factors for the XFEM state and are given by selecting and followed by and such that and are : 4 10 where is the interaction integral for and and 4 11 where is the interaction integral for and T he interaction integral in Eq. 4 6 is converted from a line integral into an area integral and a smoothing function which takes a val ue of 1 on the interior of the line integral defined by Eq. 4 6 and a value of 0 outside of the integral. Elements for the integration are selected by choosin g a radius from the crack tip. When this radius becomes sufficiently large the integral becomes path independent ; i e. any path larger than the path independent radius will yield equivalent solutions to Eq. 4 6 In practice PAGE 63 63 a radius of three elements about the crack tip is typically sufficient for path independence. An example of the elements contained within a given radius from the crack tip is given in Figure 4 1 The divergence theorem is used to create the equivalent area integral 4 12 which is simpler to implement in the finite element environment. Figure 4 1 The elements within a specified radius f or path independence. A) Radius of three elements, B) Radius of 5 elements. Other methods have also been used to extract the mixed mode stress intensity factors using XFEM. Duarte [ 147 ] used a least squares fit of the localized stress state around the crack tip to extraction the stress intensity factors. Karihaloo [ 64 ] included higher order terms of the asymptotic crack tip expansion in two dimension which allowed for the stress intensity factors to be obtained directly without the use of interaction integrals. No equi valent extension has been performed in three dimensions [ 63 ] Lua [ 126 148 ] introduced nodes along the crack fronts which were used with the crack tip opening displacement [ 15 ] to calculate the mixed mode stress intensity factors. A B PAGE 64 64 Recall that the crack tip opening displacement approach for finding stress intensity factors is only valid for linear elastic cracks. Sukumar [ 60 ] attempted to evaluate the stress intensity factors directly from the enriched degrees of freedom corresponding to the crack tip enrichment function for pure mode I problem, but t he accuracy of the stress intensity factors was found to be insufficient. A similar method was introduced by Liu [ 149 ] with much better results for homogeneous and bi material cracks. Crack Growth Direction The direction of crack propagation is accepted to be a function of the mixed mode stress intensity factors present at a crack tip. While there are several criteria available for both two [ 21 133 150 153 ] and three dimensions [ 133 153 ] in general they only differ in the angle of the initial kink, but then converge to similar crack paths [ 73 ] In two dimensions, these methods tend to give the angle of crack extension, which in general is the direction which will minimize In two dimensions the main criteria for the crack growth direction are the critical plane approach [ 133 ] maximum circumferential stress [ 151 ] maximum energy release rate [ 152 ] and the maximum strain energy density criterion [ 21 ] Other criterion available in the literature are the criterion of energy release rates and the generali zed fracture criterion [ 153 ] The criterion which is most amenable to modeling crack growth with finite elements is the maximum circumferential stress, as the growth dire ction is given in a closed form solution in terms of the mixed mode stress intensity factors. The maximum circumferential stress criterion [ 153 ] is given as the angle given by 4 13 PAGE 65 65 Al ternative, but equivalent formulas are given by Mo s [ 25 ] as 4 14 and Sukumar [ 73 ] 4 15 The critical plane approach [ 133 ] predicts t hat crack growth will occur along the angle specified by as shown in Figure 4 2 for mixed mode I and II loading as 4 16 where the angles and are 4 17 and 4 18 The parameter is the ratio of the Mode I and II crack growth rates as 4 19 where and are the stress intensity factors corresponding to a given crack growth rate The critical plane approach has the benefit of matching the observations of Feng [ 154 ] where the crack paths were different for equivalent magnitude loads in the axial and torsion dir ections. PAGE 66 66 Figure 4 2 Orientations of the critical and maximum normal stress planes. There are two options available for determining the direction of crack growth, the maximum circumferential stress and the critical plane approach. One difference between these two methods is that the critical plane approach considers the effect of the chosen material in addition to the current stress intensity factors through the use of the Mode I and II threshold str ess intensity factors, which may not be available for many materials. Due to the simplicity of the maximum circumferential stress criterion and it's acceptance in the XFEM literature [ 22 57 69 71 74 86 90 129 155 ] it will be the preferred method for determining the crack growth direction. Crack Growth Magnitude A t each cycle of constant amplitude fatigue the amount of crack growth may be on the order of nano meters In practice, the governing differential equation in Eq. 1 1 is solved at discrete points, which will be referred to as simulation iteration in this work. There are two main approaches presented in the literature for the amount of cr ack growth at a given simulation iteration for constant amplitude loading The first method assumes that a known and finite amount of growth will occur at a given iteration The second method assumes that some governing law, such as a fatigue crack growth law can be used to find the corresponding increment of growth at a particular iteration For x y Critical plane Maximum normal stress plane PAGE 67 67 the case of variable amplitude loading with an unknown correction factor model for the given crack geometry, the approaches used for constant amplitude loading are no longer valid and each cycle must be modeled in order to predict the magnitude and path of crack growth. Finite Crack Growth Increment for Constant Amplitude Loading In the literature, the selection of a constant crack growth increment is very popular [ 22 57 69 71 74 86 90 129 155 ] At a given cycle the amount of crack growth is generally very small, approximately on the order of 10 8 Therefore, it is more computationally attractive to choose a small increment of growth to represent m any cycles instead of modeling them independently. The choice of is almost always made a prior i. It has been shown in the literature by Mo s [ 25 ] Pais [ 22 ] and Edke [ 26 ] that the path of crack growth is r elated to the assumed incre ment. As the crack growth increment is reduced, crack path convergence can be observed. The convergence of crack path has also been related to the mesh density in the vicinity of the crack [ 22 26 ] Challenges with this method include selecting a such that the crack path converges to the appropriate path. With a single analysis it is unclear how one can know whether or not the predicted crack path has converged. It is also unclear how to interpolate between the data points since the crack growth curve is nonlinear. Direct Fatigue Crack Growth The slow progressive failure of a structure caused by crack growth under cyclic load ing is called fatigue. There are many fatigue crack growth models of varying complexity [ 14 ] Most fatigue crack growth models are of the f ollowing f orm : 4 20 PAGE 68 68 where d a / d N is the crack growth speed, is the stress intensity fa ctor range, and R is the stress ratio. Many different models attempt to consider the effects of phenomenon [ 14 ] such as crack growth instability when the stress intensity factor appr oaches it s critical value th reshold stress intensity factor changin g crack tip geometry and elastic crack tip c onditions For a variable load history, the interactions between overloads and underloads are a crucial component of any fatigue crack growth model An overload is a load which is substantially greater than the mean, while an underload is a load which is substantially less than the mean. Overloads increase crack tip plasticity which leads to less crack growth in resulting loading cycles. An und erload has the opposite effect and leads to increased growth. The interaction of a n overload and underl oad in subsequent cycles is not well understood. Recent work such as the state space model [ 29 30 ] small time scale model [ 24 ] and further modifications to fatigue models [ 14 27 ] has resulted in better p redictions of fatigue crack growth under variable amplitude loading. Classical Paris m odel The first fatigue model was presented by Paris [ 34 ] and will be referred to here as the classical Paris model The relationship between the crack growth rate and the stress int ensity factor range is given as 4 21 where C is the Paris model c onsta nt is the stress intensity factor range, and m is the Paris model exponent The Paris model constant and exponent depend upon material properties and loading conditions which can be derived from fatigue testing data and are fun ctions of the stress ratio [ 156 ] PAGE 69 69 The classical Paris model given by Eq. 4 21 has limitations in its ability to model fatigue crack growth. Research has shown that for a given material there is a minimum stress intensity factor range which is required to cause crack growth which is commonly referred to as the thresh old stress intensity factor range [ 14 27 157 ] Crack closure [ 14 28 32 ] has been shown to be an important quantity to consider when modeling fatigue failure but is not a part of the classical Paris model. The classical Paris model does not consider the applied stress ratio which is the ratio betw een the minimum and maximum load for a given cycle [ 14 27 31 33 ] Therefore, for the case of variable amplitude loading where the stress ratio is a function of the cycle number, each cycle would need its own value of C and m For the case of loading which is variable in nature, the interaction between over loads and under loads is also not consid ered as no variable is used which considers the intera c tion between loading cycles Modified Paris m odel One fatigue model which a ttempts to capture the behavior which is not part of the classical Paris model is a modified version of the Paris model introduced by Xiaoping [ 27 ] which will be subsequently referr ed to as the modified Paris model. Three additional parameters are added to the classical Paris model; which accounts for the threshold stress intensity factor range, which accounts for the stress rat io effect, and which accounts for the interaction between over loads and under loads. The modified Paris model takes the form 4 22 PAGE 70 70 where and are the same material properties for the classical Paris model. The value of the stress ratio correct ion factor is given as 4 23 The values of and are found by calibrating fatigue data curves from several stress ratios such that they follow the same trend as i s shown in Figure 4 3 Figure 4 3 The use of the stress ratio correction factor M R to condense the various stress ratios R into a single curve corresponding to R = 0 [ 27 ] A) Fatigue data for some stress ratios R B) Condensed fatigue data through the use of M R The value of is calculated as a function of the load interactions as 4 24 where is the plastic zone radius at the beginning of the cycle, is the crack length when the overload occurs, is the plastic zone radius when the overload occurs, is A K M R K B R = 1 .0 R = 0.0 R = 0.3 R = 0.5 R = 1 .0 R = 0.0 R = 0.3 R = 0.5 d a / d N d a / d N PAGE 71 71 the characteristic crack length for the i th cycle, is the chang e in plastic zone size caused by load interactions, and n is a shaping exponent determined by fitting to experimental data of overload and underload cases. These plastic zones are shown in Figure 4 4 where diameters are shown inst ead of radii. The plastic zone sizes follow 4 25 where is a plastic zone size factor, is the maximum stress intensity factor for the i th cycle, is the stress intensity factor when the overload occurs, is the inter action between loads given as is the tensile yield stress of the material, and is the compressive yield stress of the material. Figure 4 4 The plastic zone and crack sizes which are used in the modified Paris model to account for the effect of crack tip plasticity during load interactions. The modified Paris model [ 27 ] has been validated for a variety of over load, under load, and combinations of over load and under load cases as well as variable amplitude a OL a i d OL d y d Overload plastic zone d OL Overload underload plastic zone d OL d Current p lastic zone d y PAGE 72 72 loading [ 27 ] showing improved agreement with experimental results compared to AFGROW [ 35 ] and FASTRAN [ 84 ] models and comparable agreement with the state space model [ 29 30 ] While the modified Paris model proposed by Xiaoping [ 27 ] does not directly contain a cr ack closure model, the effects of crack closure are approximated by the tracking of plasticity associated with the interaction between the load cycles. Equivalent Stress Intensity Factor In order to apply the classical Paris model to mixed mode loading, se veral relationships have been proposed for the calculation of a single effective stress intensity factor range which can be used in Eq. 4 21 Tanaka [ 136 ] proposed the relationship 4 26 based on curve fitting observ ed experimental data. Yan [ 137 ] proposed a different correction from the maximum circumferential stress criterion as 4 27 where is given by Eqs. 4 13 4 15 Finally, a relationship based on energy release rate [ 135 ] is given as 4 28 Liu [ 133 ] propose d a variety of equivalent stress intensity factors for any combination of Mode I, II, and III loading based upon the critical plane model. For the case of mixed mode I and II loading the equivalent stress intensity factor is given as 4 29 w here k 1 k 2 and k h are related to the given mixed mode stress intensity factor values as PAGE 73 73 4 30 4 31 and 4 32 The angle in Eqs. 4 30 and 4 31 is calculated from the angle defining the critical plane as given by Eqs. 4 16 4 19 Finally, the constants A and B are found from the value of s from Eq. 4 19 as 4 33 and 4 34 The choice of an equivalent stress intensity factor model is not a trivial task. The models based on the energy release rate and proposed by Tanaka do not require any informat ion about the direction of crack propagation. Of these models, only the model of Liu considers and accounts for the effects of failure due to shear or tension according to the underlying model parameters. Therefore, this model will be the preferred choice. Integration of Fatigue Crack Growth Model For an arbitrary differential equation of the form the goal is to predict the value of while the current data is available up to The simplest numerical method available which is applicable to the current crack growth problem is the explicit forward Euler method [ 23 ] given as PAGE 74 74 4 35 where h is the step size from to Note that is the slope of at ; therefore a linear approximation is being made from to using the s lope at the This method has been popular in engineering applications because it only requires evaluating the slope at the current step, which often required expensive computational simulation. For example, in the case of crack growth simulation, corresponds to calculating the stress intensity factor with given crack size Therefore, evaluating requires a finite element modeling with the current g eometry of the crack. If a different method, such as the backward Euler method, is used, then is required, which is the stress intensity factor corresponding to the unknown crack size Therefore, ther e exists a fundamental difficulty to use a numerical integration method other than the forward Euler method. The simplest approach to integrating the Paris model for fatigue crack growth is the use of the forward Euler method. Here the stress intensity fa ctor range at the current crack geometry is the only information needed to find the increment of growth between the current and future crack geometries. The growth increment is calculated as 4 36 where is the current increment. The corresponding number of elapsed cycles can be approximated as 4 37 PAGE 75 75 for a fixed increment of crack growth. As the forward Euler method on ly uses the slope at the current point and linearly interpolates to the next crack size it can lead to large inaccuracies for relatively small crack growth increments. Although the physical meaning of is the number of elapsed cy cles, in this model it is treated as a continuous variable as crack growth occurs over tens of thousands of cycles. Choosing a larger results in needing fewer finite element simulations while sacrificing accuracy. Although the f orward Euler method can resolve the issue related to evaluating the slope at unknown crack sizes, it has a drawback in slow convergence; the accuracy of the method is proportional to the step size, h In crack growth analysis, the step size represents the number of fatigue loading cycles between two evaluation points. Since cracks grow slowly throughout the lifecycle of a product, a small step size means tens of thousands of simulations. Therefore, it is highly desired to use a numerical integration method that allows a larger step size, while maintaining accuracy. There are numerous numerical methods that allow larger step sizes for integration, but the midpoint integration method is used as a demonstration tool in this paper. The midpoint method [ 23 ] generally provides a better accuracy than the forward Euler method as it takes the slope at the midpoint between the current and future data points and uses that value to approximate the interval from to such that 4 38 As it is rare for to be evaluated directly it can be approximated as 4 39 PAGE 76 76 The accuracy of the method is proportional to h 2 ; therefore, it allows a larger step size for the same accuracy with the forward Euler method. While the midpoint approximation specifies a crack length for the evaluation at the midpoint location, no direction is specified. For a mixed mode crack growth problem, the direction of crack growth can be approximated based upon the first required function evaluation at This may limit the choice of integration step size however, as in effect an Euler approximation is used for the crack growth direction while a higher order approximation is used for th e crack growth magnitude. In terms of integrating the Paris model is the use of the midpoint method where the growth increment is calculated as 4 40 where is the midpoint between the current and next increment. Similarly, is given as 4 41 Here the slope at the midpoint of the current cycle is use to extrapolate ahead to This leads to a better approximation of the crack size at the next increment as a more accurate approximation of the slope over the entire interval of the chosen growth increment is used. This method requires additional function evaluations at for each crack geometry, effectively doubling the number of function evaluations needed for the given simulation. PAGE 77 77 With regards to the XFEM crack growth models have been limited primarily to the Paris model and modeling in terms of the crack growth increment Belytschko [ 61 ] used it to drive the crack growth for a plate with hole s, but n o specifics about how the classical Paris model was used is given. Gravouil [ 58 ] coupled the level set method to the classical Paris model but no clear definition of the number of elapsed cycles was given. Instead the growth was governed by assuming to be equivalent to the traditional time for the level set method. Sukumar [ 60 ] calculated the increment of crack growth as 4 42 where was taken as for initial crack length for a value of corresponding to the maximum value of along the crack front. Note that here a value for acts very similarly to that of a fixed such that a small value is required to ensure convergence Example p roblem An example is provided here where an XFEM analysis was performed on a plate with a hole as shown in Figure 4 5 The chosen plate dimensions were a width of 4 m and a height of 8 m with a hole with radius 1 m centered at (2,2) m. The edge crack has an initial lengt h of 0.5 m. 1.5 10 10 and 3.8 respectively. The applied stress was 14.5 MPa. The purpose of this study was to atte mpt to assess the effect that the mesh density and crack growth increments and on the predicted crack path as shown in Fi gure 4 6 PAGE 78 78 Figure 4 5 Plate with a hole subjected to tension to examine the effects of mesh density, crack growth increment, and number of elapsed cycles upon the convergence of crack path. A) Geometry, B) Example mesh Fi gure 4 6 Convergence of crack path as a function of the mesh density, crack growth increment, and number of elapsed cycles. A) Mesh density, B) a C) N From Fi gure 4 6 it is clear that the choice o f mesh size is very important as expected. Without a sufficiently refined mesh, the accuracy of the stress intensity factors will be poor leading to an inaccurate prediction of the crack growth direction. The crack path is relatively insensitive to the cho ice of in this choice, but there is deviation between the steps if a close up image of the crack paths near the hole is A B C h = 1/5 h = 1/10 h = 1/20 a = 1/10 a = 1/2 0 a = 1/4 0 N = 50 N = 25 N = 1 a A B PAGE 79 79 examined as in Figure 4 7 The most sensitive case is that of where the cases of = 25, and 50 cycles actually miss the hole as predicted by = 1. This reinforces the observations that the direct solution of a fatigue crack growth model is very sensitive t o the choice of mesh, and crack growth increment. Figure 4 7 Close up view of the crack paths for a around the hole. Cycle Counting Methods In the life estimation of a given structu re it is necess ary to convert the history of a given loading parameter for a given structure into cycles which allow for the application of a fatigue crack growth model. For structures which experience constant amplitude loading with a constant mean stress over the load i ng history, this is a trivial task. However, for a structure where the mean and amplitude of loading are non constant it is more challenging to determine what portion(s) of the given load parameter correspond with a single cycle of loading. This loading pa rameter can be any parameter of interest such as but not limited to acceleration, deflection, force, strain, stress, or torque. Cycle counting methods do not consider the effect of the loading sequence when converting a given loading parameter into an equi valent cyclic history for the same a = 1/10 a = 1/20 a = 1/40 PAGE 80 80 parameter. For many high cycle fatigue applications, the loading sequence does not significantly affect the prediction of damage accumulation from a fatigue damage model instead of a fatigue crack growth model [ 158 ] For the case of variable amplitude fatigue, it has been established that the order of loading has an effect on the crack growth prediction due to the changes in the plastic zone size caused by variable amplitude fatigue however, including the plasticity in the ev aluation of the stress intensity factors has little difference from the LEFM solution [ 37 38 159 ] There are many methods available to attempt to convert a loading history into a cyclic history such as level crossing cycle counting, peak valley cycle counting, range counting, range pair method, and ra inflow cycle counting method [ 160 165 ] Of these methods, the focus here will be upon the use of the rainflow counting method where stress is the loading param eter of interest. To this end, the MATLAB code rainflow by Nieslony [ 166 167 ] which follows ASTM E 1049, Standard Practices for Cycle Counting in Fatigue Analysis is used to convert a given stress history into a series of cyclic loads. For the case of non proportional multi axial stress histories the common practice is to first convert the multi axial stress components (e.g. and for two dimensional stress) into a single equivalent stress In order to convert the multi axial stress components into a single equivalent stress there are me thods based upon approaches such as the critical plane method [ 158 159 168 170 ] or th e Wang Brown method [ 171 172 ] This equivalent stress is then used in coordination with some fatigue damage model s uch as 4 43 PAGE 81 81 where is the amount of damage at cycle i is the number of cycles for a given applied load, is the total number of loading cycles, and is the n umber of cycles to failure at the same applied load. The cumulative damage at cycle i is given as 4 44 It is commonly assumed that once = 1, a given structure is fully damaged, can no longer carry load, and has failed [ 15 ] The fatigue damage model approach serves as the basis for the identification of cycles from multi axial stress data. For a two dimensional analysis, the equivalent bi axial stress [ 173 ] for use in a fatigue damage model can be calculated from 4 45 This equivalent stress can be used to identify the equivalent cyclic loading history. In effect, for a fatigue damage model, once the equivalent stress has been calculated the bi axial stress components (e.g. and ) can be discarded as they will no longer be used. This is not the case here as the biax ial stresses are directly applied. For variable amplitude fatigue, the interaction of the biaxial stresses will affect the magnitude and direction of crack growth and should be applied directly to a cracked plate. If a rainflow counting approach is applied to the individual stress components, it is likely that each component will have a different number of cycles. T he biaxial stress components will be applied to a cracked plate in order to calculate and Before the fatigue crack growth model can be d irectly applied, the stress history is converted into an equivalent cyclic stress history. First, the equivalent biaxial stress is calculated PAGE 82 82 by Eq. 4 45 These cycles will then be superimposed upon the individual stress components. This approach is detailed in Figure 4 8 where the conversion of the biaxial stresses to an equivalent stress, the rainfl ow counting of that equivalent stress, and the superposition of that equivalent stress onto the individual stress components is shown. Figure 4 8 The strategy used to convert the bi axial stress da ta into cyclic stress data which was suitable for use in a fatigue crack growth model. An unexpected result of this procedure was that the ratio of and are typically diff erent at the peak and valley of a given cycle. This creates a need for two function evaluations one at the maximum equivalent stress and one at the minimum equivalent stress for a particular loading cycle, to be conducted in order for the stress Calculate Equivalent Stress Rainflow Counting Data Cycle Data Data Cycle Cycle Cycle D ata xx x y yy vm vm xx x y yy Superimpose Cycles 1 2 5 6 3 4 1 2 5 6 3 4 1 2 5 6 3 4 1 2 5 6 3 4 PAGE 83 83 intensity factor range to be calculated. It should be noted that while the individual stress components are being applied instead of the equivalent stress directly, the same cycles are being followed which were identified according to th e equivalent stres s Summary Crack growth is modeled using mixed mode stress intensity factors which are extracted from XFEM analysis with the use of the domain form of the contour integrals. Four models are given for the direction of crack growth based upon the mixed mode stress int ensity factors. The most popular of these models in the computational fracture mechanics community is the maximum circumferential stress criterion. For cases where shear may be a significant source of failure the Liu model may provide a better match to exp erimental data. Two of the many fatigue crack growth models are introduced, the classical Paris model and a modified Paris model which has additional parameters to consider the effects of the stress ratio R threshold stress intensity factor and changes in the crack growth rate caused by the evolution of the plastic zone ahead of the crack tip under variable amplitude loading. Four models are presented for the calculation of an equivalent stress intensity factor range from the mixed mode stress intensity factors. The model proposed by Liu will be used here based upon it's superior performance when compared to a range of experimental data for failure caused by shear and tension. The challenges with solution of a fatigue crack gro wth model were highlighted by an example of fatigue crack growth in a plate with an edge crack an a hole. The effects PAGE 84 84 of mesh density, crack growth increment, and the elapsed number of cycles were considered with respect to the convergence of the crack pro pagation path. Finally, the rainflow counting method was introduced for the conversion of stress histories into equivalent cycles. As these models are usually used for fatigue damage prediction instead of direct fatigue crack growth modeling, the steps wh ich are used in the identification of cyclic stress histories for the biaxial stresses is detailed. First, the equivalent biaxial stress is calculated. Cycles for the equivalent biaxial stress are calculated according to the rainflow counting method. These cycles are then superimposed upon the biaxial stresses, allowing for the biaxial stress components to be applied directly for the analysis of a cracked plate subjected to variable amplitude biaixal loading. PAGE 85 85 CHAPTER 5 SURROGATE MODELS FOR HIGHER ORDER INTEGRATION OF FATIGUE CRACK GROWTH MODELS Integration of Ordinary Differential Equati on T he objective of this chapter is to use a surrogate model to reduce the number of expensive function evaluations needed in the numerical integration of an ordinary differential equation which follows the general form 5 1 For an unknown function f ( x y ) a common method to approximat e the solution is to discretize f ( x y ) in to discrete points, w here f ( x y ) is approximated using some numerical technique. In this work, finite element simulations are used to approximate f ( x y ). There are numerous methods available to approximate the solution to Eq. 5 1 at some from some current state such as the forward Euler approximation which is 5 2 where h is the integration step size given as x i + 1 x i or the midpoint approximation as 5 3 The midpoint approximation allows for larger step sizes h to be taken compared to the forward Euler app roximation for the same integration accuracy. However, this requires two function evalu a tions one at x i and another at x i +1/2 which is one more than the forward Euler approximation requires. For many engineering problems, no analytical equation is avail able for the right hand side given as In some of these cases can be approximated at PAGE 86 86 discrete locations through the use of finite element analysis simulations. However, the se simulations can be expensive, especially for complex geometries. In the case of fatigue crack growth, the number of analysis needed further increases the cost of the finite element simulations needed to approximate In this c hapter, two cases with various amounts of knowledge about are considered. First, the case where a history of discrete function evaluations is available is considered A surrogate model is then fit to the discrete function evalua tions and is used to interpolate the values of required by the midpoint approximation, replacing finite element simulations at these points. Second, the case where only the discrete points up to ( x i y i ) a re available is consider ed. These two cases are shown in Figure 5 1 where ( x f y f ) is some final value of ( x y ) which ends a given analysis. A surrogate model is fit to the available data points and is used to extrapolate to r eplacing an expensive finite element simulation at that point. The methods introduced allow for the benefits of the midpoint approximation to be realized with the number of function evaluations needed for the forward Euler approximation. Figure 5 1 Available discrete points available which can be fit using a surrogate model. A) Complete history is available to ( x f y f ) B) History is available to ( x i y i ). x f x y x A B y i y f y y f x f x i Data point Data point PAGE 87 87 Surrogate Models Surrogate models are commonl y used to reduce the cost associated with the evaluation of expensive functions. Types of surrogate modes include polynomial response surface [ 174 ] kriging [ 175 177 ] radial basis neural networks [ 178 ] and support ve ctor regression [ 179 ] Kriging is used here instead of other surrogate models based on an initial trial of various surrogate models. The results of these tests showed that all surrogates besides kriging would require additional tuning before they would provide ac ceptable accuracy, while kriging simply worked. While the results here were generated with kriging, the methods presented here should be applicable to any surrogate model which provides sufficient accuracy for interpolation and extrapolation. Here three ca ses are considered one for the case in Figure 5 1 where a compete history of discrete points is available and two where only the history up to some ( x i y i ) is available First, the approach of fitting surrogate models to the stre ss intensity factor and crack growth direction history is tested based on a fixed increment of crack growth or a fixed number of elapsed cycles for a given crack geometry. When a fixed is used the relationship between and is fit using a kriging surrogate. This surrogate is then used to approximate the corresponding cycle number N for a given crack size a such that the cyclic crack growth history can be obtained. Correspondingly, when a fixed is used the relationship between and is fit using a kriging surrogate. This surrogate is then us ed to extrapolate forwar d in t he cycle history so that a higher order approximation to the chosen fatigue crack growth law can be used to find the corresponding Second, a variable integration step size algorithm is introduced f or the automatic generation of a step size based on the accuracy in the surrogate PAGE 88 88 extrapolation of and the crack growth direction. This is an attempt to only perform expensive function evaluations whe n needed. These two cases are shown in Figure 5 2 where K c is the critical stress intensity factor range, a c is the corresponding critical crack size, K i is the stress intensity factor range corresponding to the crack size at cy cle N i Third, the case of mixed mode crack growth is considered and the use of a surrogate model is introduced which enables the use of a higher order approximation to the crack growth direction which would otherwise not be available. Figure 5 2 Surrogate model fit to either a complete or partial crack growth history. A) Complete history, surrogate for a K B) Par t ial history, surrogate for N K Kriging Kriging [ 175 177 ] can be used to represent a function of interest As this function is expe nsive to evaluate, it may be approximated by a cheaper model based on assumptions on the nature of and on the observed values of at a set of data points called an experimen tal design More explicitly 5 4 where is a real dimensional vector and represents both the error of approximation and random errors. a c a N A B K i K K I C N f N i Data point Su rrogate Data point Surrogate K K IC PAGE 89 89 Kriging estimates the value of the unknown function as a combination of basis fu nctions such as a polynomial basis and departures by 5 5 where satisfies for all sample points and is assumed to be a realization of a stochastic process with mean zero, 5 6 and process variance and spatial covariance function given by 5 7 where is the correlation bet ween and is the value of the actual responses at the sampled points is the Gramian design matrix constructed using the basis functions at the sampled points is the matrix of correlations among sample points and is an approximation of the vector of coefficients of Eq. 5 5 Figure 5 3 shows the prediction and the error estimates of kriging. It can be noticed that since the kriging model is an interpolator, the erro r vanishes at sampled data points. It may be possible to use the variance provided by kriging as a means to assess its accuracy or to find an appropriate integration step size. Since many surrogate models do not provide an estimate of the variance, this ap proach has not been used. The goal was to create a general algorithm applicable to any surrogate model PAGE 90 90 Figure 5 3 Comparison of an exact function and the corresponding k riging surrogate model for an arbitrary set of five data point s Kriging for Integration of Fatigue Crack Growth Law Constant Integration Step Size The limitation on the types of methods which can be used in the direct integration of a fatigue model come s from the need to evaluate t he stress intensity factor at some future crack tip position. The idea presented here is to use a surrogate model to fit the available stress intensity factor history and to replace some of the expensive finite element simulations with in expensive functio n evaluations using the surrogate model. Two cases are considered : one in which there is a constant magnitude of crack growth at each integration step, another where a constant number of elapsed cycles occurs for each integration step For t he case of a fixed is calculated a posteriori to the function evaluations using either an Euler approximation as Data point Exact function Kriging surrogate x y PAGE 91 91 5 8 or a midpoint approximation as 5 9 With a traditional approach the function evaluation has not yet been completed and would require an additional expens ive function evaluation. Through the use of surrogate models this extra function eval ua tion at can be avoided by fitting a surrogate model to the history The extra needed function evaluation can then be approximated using interpolation. Four initial data points are provided before the use of kriging. To get the three additional data points, the forward Euler method can be used with a small integration step size (e.g. ) which has sufficient accuracy away from the critical crack length. When a constant number of elapsed cycles is considered the corresponding amount of crack growth for a given ca n be calculated using the forward Euler approximation as 5 10 or the m idpoint approximation as 5 11 In Eq. 5 11 not e that there are two expensive function evaluations, one for and the other at the approximated value for The use PAGE 92 92 of surrogate models allow s for only one expensive function evaluation to be required at A surrogate model is fit to the history up to and including The surrogate model is then used to extrapolate to replacing an expensive finite element simulation as shown in Figure 5 4 For problems with mixed mode loading, two kriging surrogates are used. First, a surrogate is used to approximate the respons e of as a function of either the crack length or the number of elapsed cycles The surrogate approximation is then used for a higher order approximation to Paris model suc h as the midpoint method given by Eqs. 5 9 and 5 11 Second, a surrogate is fit to the h istory for the history of crack growth direction in the global coordinate system. The same extrapolation and interpolation strategy using the surrogate model is used to enable higher order approximations for the crack growth dir ection. Figure 5 4 Integration of fatigue crack growth model from N i to N i +1 using either Euler or midpoint approximation. A) Function evaluations needed for Euler and midpoint approximations, B) A ccuracy of crack size calculation at N i +1 Variable Integration Step Size There are numerical methods which will automatically adjust the integration step size to some allowable error [ 23 ] However, these methods typically use function Data point K ( a i ) K ( a i+ 1/2 ) K N i N i + 1 N N a i +1 Euler Midpoint a N i N i + 1 N N A B PAGE 93 93 evaluations to estimate the error associated with a given step size before choosing upon a step size to use. In the finite element framework, this algorithm would result in needless expensive function evaluations in order to determine the allowable step size. To avoid these additional function evaluations, an algorithm is presented which allows for the step size to be dynamically changed based on the accuracy of the surrogate extrapolation as shown in Figure 5 5 This algorithm should be useful for any surrogate model. For the case of adjusting the step size based on the prior accuracy, the i +1 th step size can be found from the solution for the i th cr ack geometry as 5 12 where is the current step size, is the target allowable error in percent difference between the surrogate extrapolation and XFEM values at i +1, is an exponent which determines how quickly the step size changes, and and are evaluated based on i +1 for the current crack geometry using either a numerical (e.g. XFEM) or surrogate (e.g. kriging) model. The step size is scaled based on the accuracy of the predicted value. Fo r the case where either = or = then = The values of and in Eq. 5 12 were found to be 0.001, while and were defined to be 0.1 based on a parameter study where a range of materials with unique C and m values and a variety of different initial crack sizes were used in the calibration of the variable integration step size algorithm When the error shown in Figure 5 5 is greater than the target error, the integration step size is decreased Similarly, when the error is smaller than the target error, the integration step size is increased. PAGE 94 94 Figure 5 5 Errors used to adjust the adaptive integration step size algorithm. A) Error used to adjust integration step size according to stress intensity factor surrogate, B) Error used to adjust integration step size according to crack growth direction surrogate. A summary of this appro ach follows. A series of XFEM simulations were performed to evaluate as a function of the crack size In addition, the angle of crack growth given in terms of the global coordinate system was evaluated as a function of crack size. Once these easily evaluated functions were available, the crack growth was simulated such that: 1. For crack geometry i fit surrogate for history of and evaluate and at i +1/2 and i +1. 2. Calculate the crack growth increment from 3. Calculate the crack growth increments in the x and y dire ctions from and 4. Update the crack tip location based on and evaluat e and 5. Find the next step size from Eq. 5 12 repeat for next crack geometry. It is noted that the current dynamic step size determination scheme is not recovering the error at step i +1. K i+ 1 (NUM) K i+ 1 (SUR) K i+ 1/2 (SUR) K N i N i + 1 N N i A N i+1 error i+ 1 (NUM) i+ 1 (SUR) i+ 1/2 (SUR) N i N i + 1 N N i B N i+1 error PAGE 95 95 Examp le Problems First a center crack in an infinite plate under uniaxial tension is considered. This case is convenient as i t has a closed for m solution for the crack size and stress intensity factor as a function of applied stress, loading cycle, and material properties. This allow for an estimate of the amount of error that can be expected by the use of kriging to provide data points instead of performing a function evaluation. The effect on the choice of a fixed or is considered for the theoretical model of a center crack in an infinite plate to show the validity of the approach. The presented variable step size algorithm is also used for each example problem. Results are presented for each case as the f inal approximate value normalized by the exact final value of crack size or cycle number. A value of one denotes perfect agreement with the exact solution. The tested materials and the corresponding material properties are aluminum 2024 at stress ratios of 0.1 ( C = 1.60 10 1 1 m = 3.59, = 30 ) and 0.5 ( C = 3.15 10 11 m = 3.59, = 30 ) as well as austenitic ( C = 1.36 10 10 m = 2.25, = 50 ) and martensitic ( C = 5.60 10 1 2 m = 3.25, = 50 ) steel [ 156 ] The resul ts are presented for each case as a function of the ratio between the theoretical solution and the calculated solution for either crack size or cycle number such that a value of one denotes no difference between the calculated and theoretical solution s. A value of 0.98 or 1.02 represents a 2 percent error in the calculated value. Center Crack in an Infinite Plate Under Tension First, a center crack in an infinite plate with initial crack length of 10 mm is considered. For the case of a center crack in an in finite plate under uniaxial tension of 78.6 MPa, the Mode I stress intensity factor [ 20 180 ] is PAGE 96 96 5 13 where is the applied stress, and is the half crack length. By substituting Eq. 5 13 into Eq. 4 21 for Paris model and rearranging terms, the number of cycles for a crack to grow can be calculated using the following integration: 5 14 where is the initial crack size. Rearranging Eq. 5 14 the crack size after cycles can be calculated by 5 15 This allows for the midpoint method to be applied as and can be directly evaluated and used in the integration of Paris model at any cycle number. F or a geometry without known relationship between and a finite element solution would need to be performed at at which time an Euler approximation could be used to approxima te a i +1 allowing a second finite element simulation to be performed and the midpoint method to be applied. An aluminum 2024 plate with a stress ratio of 0.5 is chosen to validate the use of kriging extrapolation for the integration of the fatigue crack gr owth model. The critical stress intensity factor is reached for a crack size of 50 mm given an initial crack size of 10 mm in about 22,300 cycles. For a fixed the results are given in Table 5 1 where t he data is presented as normalized values of the cycle at which the crack length is 50 mm. The kriging surrogate is used to calculate cycle number N for each discrete crack PAGE 97 97 length a as a post processing operation. Comparing the midpoint approximation to th e kriging midpoint approximation provides an indication in the amount of error introduced by using kriging to replace the additional function evaluation needed for the midpoint approximation Thus, all kriging evaluations are based on interpolation. Recall that the common used in the literature is a i /10. For that crack growth increment, the Euler approximation has over 5 percent error. Through the use of the midpoint approximation, the error is reduced to less than 1 percent. Lar ger crack growth increments for the Euler method lead to very large errors, which can be drastically reduced through the use of the midpoint approximation. As interpolation is being used here it would also be possible to apply other higher order approximat ions such as the Runge Kutta method [ 23 ] for the back calculation of elapsed cycles for a fixed Table 5 1 A ccuracy of kriging interpolation integration for fixed increment a for a ce nter crack in an infinite plate of Al 2024 with R = 0.5. The number of function evaluations is denoted as FE. Normalized N ( a = 50 mm) Eule r Euler FE Midpoint Midpoint FE Kriging Midpoint Kriging FE a i /160 1.00 640 1.00 1280 1.00 640 a i /80 1.01 320 1.00 640 1.00 320 a i /40 1.01 160 1.00 320 1.00 160 a i /20 1.03 80 1.00 160 1.00 80 a i /10 1.05 40 0 .999 80 0.999 40 a i /5 1.11 20 0.997 40 0.99 7 20 a i /2 1.30 8 0.981 16 0.981 8 a i 1.66 4 0.935 8 0.944 4 For a fixed the results are given i n Table 5 2 where the data is presented as normalized values of the crack length at cycle number 23,00 0. The kriging surrogate is used to extrapolate a stress intensity factor range allowing the midpoint method to be applied. The first observation from Table 5 2 is that the accuracy of Euler approximations is much more sensitive t o changes in fixed when compared to a PAGE 98 98 fixed For this problem, the function evaluations are equal for each approach, in general midpoint would be twice as expensive as Euler or the kriging midpoint. H owever, the midpoint approximations are largely insensitive to the chosen crack growth increment. There is some loss of accuracy associated with the use of the midpoint approximation method when compared to the exact midpoint approximation. It should be no ted that the kriging surrogate produced improved accuracy when compared to the midpoint approximation method. The kriging assisted midpoint method offers improved accuracy to the midpoint approximation method with less function evaluations. Table 5 2 A ccuracy of kriging extrapolation integration for fixed increment N for a ce nter crack in an infinite plate of Al 2024 with R = 0.5. The number of function evaluations is denoted as FE. Normalized a ( N = 23,000) Euler Euler FE Midpoint Midpoint FE Kriging Midpoint Kriging FE 1 1.00 23,000 1.00 46,000 1.00 23,000 25 0.995 920 1.00 1,840 1.00 920 50 0.990 460 1.00 920 1.00 460 100 0.981 230 0.999 460 1.00 230 500 0.918 46 0.997 92 0.999 46 1000 0.856 23 0.989 46 0.995 23 To assess the applicability of the variable step size algorithm to a range of fatigue problems four different materials were chosen along with initial crack sizes of either 1 or 10 mm and grown to failure at 50 mm. As with the fixed step size approach, three data points are found using the forward Euler approach with = 1 before the variable step size algorithm begins. The results for the materials and crack sizes is given in Table 5 3 where a comparison between the cycle where failure occurs based on the Euler and kriging extrapolation methods with the number of f unction evaluati ons Note that the number of function evaluation s for the Euler method is the same as It can be noted from the above table that the number of function evaluations that PAGE 99 99 the algorithm uses to model the crack growth to failure is generally independent of the number of cycles to failure. The general trend of the variable step size algorithm is that initially there are smaller steps, then the step size increases. When the crack approaches the critical crack size, the integration step size decreases ensuring an accurate solution. T his behavior is shown in Figure 5 7 and while the problem being solved is different the trend of the integration step size curve is the same. Table 5 3 Effect of material and initial crack size on est imated cycles to failure for the variable step size algorithm for a center crack in an infinite plate. Material Kriging FunEval Norm N fail Aluminum, R = 0.1 1 365,740 370, 292 52 1.01 Aluminum, R = 0.1 10 44,301 44,358 49 1.00 Aluminum, R = 0.5 1 185,773 184,286 46 0.992 Aluminum, R = 0.5 10 22,502 23,121 50 1.02 Aust Steel 1 808,426 813,819 61 1.01 Aust Steel 10 285,645 292,556 55 1.02 Mart Steel 1 2,103,878 2,085, 640 67 0.991 Mart Steel 10 346,500 347,733 55 1.00 Edge Crack in a Finite Plate Under Tension For the case of an edge crack in a finite plate under uniaxial tension of 78.6 MPa the Mode I stress intensity factor [ 180 181 ] is 5 16 where is the crack length and is the plate width. The plate width for this problem was W = 0.2 m. As th ere is no analytical expression for the crack length after 10,700 cycles was found to be 30 mm using the forward Euler method where the step size was reduced until convergence in crack size at 10,700 cycles was achieved. For a fixed step size, aluminum 2024 with a stress ratio of 0.5 is chosen to validate the use of kriging extrapolation for the integration of the fatigue crack growth model with PAGE 100 100 a constant step size. The critical stress intensity factor is reached for a crack si ze of 5 0 mm given an initial crack size of 10 mm in about 11,000 cycles. For a fixed the results are given in Table 5 4 as normalized values of the cycle at which the crack length is 50 mm. The same nu mber of function evaluations is shared by all integration methods. The kriging surrogate is used to calculate cycle number N for each discrete crack length a as a post processing operation. Thus, all kriging evaluations are based on interpolation. The midp oint method was used here by evaluating Eq. 5 16 for any crack size a This gives an error estimate for kriging interpolation. Here for a fixed crack growth increment of a i /10 the corresponding Euler approximation has an error of 7 percent, while the midpoint approximation yields less than 1 percent error. Table 5 4 A ccuracy of kriging interpolation integration for fixed increment a for a n edge crack in a finite plate of Al 2024 with R = 0.5. The number of function evaluations is denoted as FE. Normalized N ( a = 50 mm) Euler Euler FE Midpoint Midpoint FE Kriging Midpoint Kriging FE a i /160 1.00 640 0.996 1280 0.996 64 0 a i /80 1.00 320 0.996 640 0.996 320 a i /40 1.01 160 0.996 320 0.996 160 a i /20 1.03 80 0.995 160 0.996 80 a i /10 1.07 40 0.995 80 0.994 40 a i /5 1.14 20 0.991 40 0.991 20 a i /2 1.39 8 0.969 16 0.973 8 a i 1.88 4 0.905 8 0.944 4 For a fixed elapsed numb er of cycles in each function evaluation, the exact results for the midpoint method are not available as there is no explicit value for From Table 5 5 where data is presented as normalized values of th e crack length at cycle number 11,000 it is apparent that once again, the kriging assisted midpoint method allows for larger step sized compared to the forward Euler method. The number of function evaluations is equal for Euler and kriging integration. The kriging surrogate is PAGE 101 101 used to extrapolate a stress intensity factor range allowing the midpoint method to be applied. The exact midpoint method may not be applied to this problem as no closed for solution for a ( N ) is known to the authors for the stress int ensity factor given by Eq. 5 16 In this case, for a step size of 100, the errors for the Euler and kriging assisted midpoint methods are 3.5 and 0.2 percent. The increased errors in the approximation compared to the case of a center crack in an infinite plate can most likely be explained by the increased nonlinearity caused by the edge and finite effects in this geometry. Table 5 5 A ccuracy of kriging extrapolation integration for fixed increment N for a n edge crack in a finite plate of Al 2024 with R = 0.5. The number of function evaluations is denoted as FE. Normalized a ( N = 11 ,000) Euler Euler FE Mi dpoint Midpoint FE Kriging Midpoint Kriging FE 1 1.00 11,000 1.00 22,000 1.00 11,000 25 0.994 440 1.00 880 1.00 440 50 0.987 220 1.00 440 1.00 220 100 0.975 110 0.999 220 0.999 110 500 0.947 22 0.996 44 0.995 22 1000 0.877 11 0.984 22 0.978 11 The variable step size algorithm is again tested with a range of fatigue problems. Four different materials were chosen along with initial crack sizes of either 1 or 10 mm and grown to failure at 50 mm. As with the fixed step size approach, three data points a re found using the forward Euler approach with = 1 before the variable step size algorithm begins. The results for the materials and crack sizes are given in Table 5 6 A comparison between the cycle wh ere failure occurs based on the Euler and kriging extrapolation methods is given in addition to the number of required function evaluations. Note that the number of cycles to failure does not change t he number of kriging function evaluations used by the variable integration algorithm. PAGE 102 102 Table 5 6 Effect of material and initial crack size on estimated cycles to failure for the variable step size algorithm for an edge crack in a finite plate. Material Kriging FunEval Norm N fail Aluminum R = 0.1 1 235,374 240,401 42 1.02 Aluminum R = 0.1 10 20,701 21,473 38 1.03 Aluminum R = 0.5 1 119,557 12 2,705 44 1.03 Aluminum R = 0.5 10 10,516 10,799 38 1.03 Aust Steel 1 594,634 601,308 49 1.01 Aust Steel 10 189,225 193,493 41 1.02 Mart Steel 1 1,415,883 1,459,668 46 1.03 Mart Steel 10 196,924 200,215 41 1.01 Inclined Center Crack in a Finite Plate Under Tension The first test problem which also considers the effect of crack growth direction is that of an inclined center crack in a finite plate subjected to uniaxial tension in the y direction as shown in Figure 5 6 Th e magnitude of the applied stress was material dependent to reduce the number of cycles to failure for the steel samples and therefore, the number of function evaluations for the forward Euler solution with step size = 1. The plate was chosen to be a 2 m x 2 m plate with an initial half crack size of 0.187 m. The crack was grown by XFEM simulations on a structured mesh of element size h = 0.05 to a final half crack size of 0.6 m. The values of the cons tants in the variable step size algorithm were retained from the previous sections. For the variable step size algorithm, different applied stresses were considered for different materials in order to reduce the number of XFEM simulations needed for the Eu ler method The comparison of the number of cycles to failure and the final crack position for the Euler with constant to the variable step size is presented in Table 5 7 Note that the cycle number corresponding to failure in each case is in e xcellent agreement with the reference Euler solution. For the case of austenitic steel, it should be noted that the larger difference between the number of cycles to failure and the final PAGE 103 103 coordinate s of the crack tip can be explained by the last step simply occurring in a way which passed the true failure cycle. Figure 5 6 An inclined center crack in a square finite plate under uniaxial tensi on. A) Initial geometry, B) Final geometry. Table 5 7 Comparison of Euler and variable step size predictions for the coordinates ( x final y final ) of the right crack tip marked with a dotted circle in Figure 5 6 and cycle number N fail corresponding to a final crack size of 0.6 m from an initial crack size of 0.187 m for four different materials. Material R = 0.1 R = 0.5 Austenitic Martensitic Applied Stress MPa 50 50 125 75 x fail Euler 1.5431 1.5432 1.5431 1.5431 Variable 1.5432 1.5433 1.54 5 0 1.5431 y fail Euler 1.0902 1.0902 1.0902 1.0902 Variable 1.0902 1.0902 1.0900 1.0902 N fail Euler 29 519 14 996 68 331 80 597 Variable 29 574 15 025 68 622 80 707 FunEval VAR 356 157 52 245 The integration step size and percent errors as a function of cycle for two stress ratios of aluminum 2024 are given in Figure 5 7 It can be observed that for each material the agreement in final crack tip position and number o f cycles to failure agrees very well to the Euler solution with small step size. Note that the integration step size y x W y x a W A B PAGE 104 104 increases in Figure 5 7 until the target error is exceeded, at which point in time, the integration step size dec reases to limit the amount of error introduced into the solution. This phenomenon can be noted as occurring for aluminum for R = 0.1 around cycle number 11,000 where the integration step size begins to decrease. Also for aluminum at stress ratio R = 0.5 th e reduction in integration step size occurs twice, first around cycle number 4,000 and second around cycle number 10,000. A B C D Figure 5 7 Comparison of cycle number and integration step size or percent error for surro gate for R = 0.1 and R = 0.5 for aluminum 2024. A) Integration step size for aluminum 2024, R = 0.1, B) Integration step size for aluminum 2024, R = 0.5, C) Percent error for aluminum 2024, R = 0.1, D) Percent error for aluminum 2024, R = 0.5. (Note when e rror exceeds target, step size reduces). PAGE 105 105 Summary The ordinary differential equation which governs fatigue crack growth can not easily be solved due to the fact that in general no analytical expression of the stress intensity factor as a function of crack si ze and crack geometry is available. For a general geometry, the stress intensity factor can only be calculated at the current crack size using the extended finite element method. It is common to take a forward Euler approach to the direct calculation of bo th the magnitude and direction of crack growth for a given function evaluation or for the back calculation of the elapsed number of cycles. The forward Euler approach limits the step size and increases the number of function evaluations needed to achieve s ufficient accuracy. A kriging surrogate was used to fit the available history of the stress intensity factor and crack growth direction. The surrogate model allows one to extrapolate ahead of the current data point, enabling the use of the midpoint approxi mation method for the evaluation of both stress intensity factor and crack growth direction. A variable step size algorithm is also presented to allow for the step size to be adjusted for a given allowable error. This method produces comparable accuracy wi th a significantly reduced number of function evaluations for a range of initial crack sizes, applied stresses, materials, and crack geometries. Another area of study is to consider cases of mixed mode crack growth. In this case the path of crac k growth wi ll be curved. Thus, the allowable step size based on the approximation of the solution of the fatigue crack growth law may be limited by the need to have a converged crack path. This is an important interaction which will need to be considered for practica l use of the kriging assisted methods for use in a finite element environment. PAGE 106 106 C HAPTER 6 AN EXACT EXTENDED FINITE ELEM ENT METHOD REANALYSIS ALGORITHM Reanalysis Methods Crack growth is often modeled as a quasi static problem. In this methodology, only a small portion of the XFEM stiffness matrix will change as crack growth occurs. Once Heaviside terms are added to the XFEM stiffness matrix in order to model the discontinuity due to a crack they will remain unchanged in future iterations of crack growth. At each iteration of crack growth, the previous crack tip components will need to be modified to account for the changing crack tip position. The objective of this paper is to maximize recycling computational resources during crack growth simulation by using this property. For the first iteration of the quasi static solution procedure th e Cholesky factorization of the full XFEM stiffness matrix is calculated. In subsequent iterations, only the changed portion of the XFEM stiffness matrix is calculated (i.e. the new Heaviside and crack tip enrichment components). The existing Cholesky fact orization is directly modified using row add and row delete operations. This solution procedure significantly reduces the cost of assembling the stiffness matrix and solving the system of matrix equations in the XFEM framework, allowing for increased conve rgence of crack paths within a computational budget. For many optimization applications involving XFEM, the same basic algorithm may be used to reduce the computational cost associated with the repeated iterations within the optimization cycle. In the case of damage identification [ 99 100 182 ] al l enriched elements may evolve over time (i.e. the Heaviside enrichment is non constant). For shape optimization [ 140 ] the changes to the elements along the shape boundary being modified can be updated along with any enriched degrees of freedom being considered. PAGE 107 107 Reanalysis algorithms [ 183 184 ] have been developed primarily for use in the fields of design and optimization to efficiently solve problems where small perturbations to the finite element domain are made. This may take the form of changing the location of elements, adding ad ditional elements, or a combination of both. These methods can be classified as either being exact or approximate [ 185 ] The exact methods are generally based on the Sherman Mo rrison [ 186 ] inversion formula and consider cases where the modified elements affect a small number of degrees of freedom. The approximate methods are typically iteration based and are used when the modified elements affect a large number of degrees of freedom. Li and Wu [ 184 185 187 ] have introduced algorithms for the reanalysis of structures with added degrees of freedom through the use of an iterative solver. An exact reanalysis algorithm is introduced here based on the direct modification of an existing Cholesky f actorization [ 188 194 ] In this paper Cholesky factorization was chosen instead of a Sherman Morrsion approach due to collaboration with individuals containing ex tensive knowledge of Cholesky factorizations. Exact Reanalysis by Modification of an Existing Cholesky Factorization Recall that the displacement approximation and stiffness matrix for the XFEM with a crack is given as 6 1 and 6 2 PAGE 108 108 It can be noticed from E qs. 6 1 and 6 2 that the stiffness component associated with the traditional finite element approximation is not a function of the crack location, which implies that the component of the stiffness matr ix will be constant at each iteration of crack growth. This implies that the changing portion of the stiffness matrix is limi ted to the enriched portion, which will be small compared to the un enriched portion. F urthermore, it can also be noticed that while the Heaviside enrichment term s are function s of the crack location within an element, that once a Heaviside enrichment or coupling to the classical FE approximation is introduced its value will not change in any future iterations T he stiffness components and will be constant for future iterations of crack growth. An example of the nodes whose corres ponding stiffness components need to be updated to account for crack growth is shown in Figure 6 1 Once the boundary conditions have been applied, the stiffness matrix becomes positive definite and it is possible to compute the C holesky factorization for the XFEM stiffness matrix. The Cholesky factorization for the XFEM stiffness matrix in Eq. 6 2 is given by 6 3 where the components of and can be calculated as 6 4 and 6 5 PAGE 109 109 In the case of modeling crack growth in a quasi static manner, there will be many iterations of cr ack growth. Each iteration will require its own Cholesky factorization and corresponding floating point operations for factoring A B Figure 6 1 Finite element mesh with enriched nodes. A) In itial crack geometry, B) First iteration of crack growth. (Note c ircles denote nodes enriched with the Heaviside function and squares denote nodes enriched with th e crack tip enrichment function, f illed circles and squares denote new enriched nodes for the current iteration while open circles and squares denote previously enriched nodes ). As the XFEM stiffness matrix is largely constant and slowly evolving (i.e. additional Heaviside and changing crack tip stiffness components), it would be computationally a ttractive if the Cholesky factorized matrices, and are directly updated. The modification of the initial Cholesky factorization through sparse row add and row delete operations [ 188 191 194 ] requires fewer floating point operations for subsequent iterations of crack growth. Suppose the i th row and column of is equal to the i th row/column of the identity matrix. This i th row has no effect on the Cholesky factorization of the rest of the matrix. We refer to replacing this i th row/column of with a nontri PAGE 110 110 1 downdate, as described trivial row/column of with the corresponding row/col umn of the identity matrix. An alternative approach would be to append newly enriched degrees of freedom at the bottom/right of the system of equations, and thus the size of and correspondingly and would dynamically grow. One could find the fill reducing ordering based on the initial and then simply add the new rows and columns of and to the bottom right corner of the matrices. This will introduce drastic fill in since the original fill reducing ordering did not account for the new entries. Another option would be to create a matrix that assumes that all nodes in contain both the Heaviside and crack tip enrichment. The size of this matrix is approximately ten times the size of for the case of a two dimensional structure. Any enriched nodes that are inactive hav e their corresponding rows and columns set to zero with unity along the diagonal. Therefore, all the inactive enriched nodes do significantly affect the size of the resulting sparse matrix. If the possible crack growth path is known a priori it would be p ossible to limit the number of initially inactive enriched nodes that are considered. In order to find the fill reducing ordering for this matrix, the connectivity of all possible elements is used to create a matrix with unity i n all matrix locations which would contain information if all possible Heaviside and crack tip nodes were active. The nonzero pattern of is a superset of any that would be seen in the subsequent numer ical factorizations computed during crack propagation, and thus the nonzero pattern of its Cholesky factor will also be a s uperset of any Cholesky factor seen later. A fill reducing ordering is obtain ed for from the approximate minimum degree PAGE 111 111 ordering (AMD) algorithm [ 189 195 196 ] Finding an optimal fill reducing ordering for an arbitrary is computationally impossible, and thus AMD is a heuristic, as are all fill reducing algorithms. However, since AMD is applied to the entire matrix, it will limit the fill in in the matrix and thus it will also limit the upper bound of any fill in in any Cholesky factor computed later. If AMD were simpl y applied to the first then new entries in that AMD did not consider could cause catastrophic fill in in After applying the fill reducing ordering, the permuted stiffness matrix is referred to as The fill reducing ordering is applied to reduce the fill in of subsequent iterations of crack growth as the active enriched nodes evolve. For elements that are initially intersected by the crack, the c orresponding enriched degrees of freedom are set active, while all other nodes are defined inactive. Any rows or columns of and subsequently corresponding to an inactive enriched degree of freedom h ave rows and columns which are set to zero, except for their diagonal component which is unity. In order to update the Cholesky factorization, only the modified components of the updated stiffness matrix (e.g. new Heaviside and crack tip components) are needed. Before updating the previous Cholesky factorization, the fill reducing ordering is applied to yielding The initial Cholesky factorization is 6 6 PAGE 112 112 where the identity row denotes an inactive enriched node and is the XFEM stiffness matrix permuted to a fill reducing ordering. After crack growth has occurred and the previously inactive node becomes active the new stiffness matrix is given as 6 7 The modification of the initial to add a new enriched degree of freedom to the C holesky factorization is accomplished as: 1. Solve the lower triangular system for 2. 3. 4. 5. Perform the rank 1 downdate The total time taken for steps 1 to 3 is proportional to the sum of the number of entries in each column j of ( and ) for which the j th entry in is nonzero. This is an optimal result. The time taken for step 4 is proportional to the number of entries in each column j in for which is nonzero. The time taken for this step is e xactly proportional to the number of entries in that change. This time includes both the symbolic work of figuring out which entries change, and then doing the numerical work to change them. While this seems like an obvious goal of a rank 1 downdate, it is far from trivial to obtain [ 188 191 194 ] PAGE 113 113 Crack tip stiffness components will need to be modified after each iteration of crack growth, therefore it is necessary to also consider the case of a row delete operation to remove old entries into the stiffness matrix associated with the crack tip enrichment function. In this case from th e previous iteration the stiffness matrix and corresponding factorization are given as 6 8 where the second row and column correspond to a enriched crack tip entry in the XFEM stiffness matrix. The corresponding stiffness matrix and factorization with the old crack tip component removed is given as 6 9 The modification of the initial to remove a crack tip enriched degree of freedom from the Cholesky factorization is accomplished as: 1. 2. 3. 4. 5. Perform the rank 1 update Just as the rank 1 update, above, the total time taken by steps 1 to 5 is proportional to the number of entries in that change. PAGE 114 114 The new stiffness components associated with the crack tip enrichment function are added using the row add operation pr eviously discussed. For more details on the row add and row delete algorithms, readers are referred to the work of Davis [ 193 194 ] including proofs of optimal convergence. Our im plementation was accomplished based on the creation of a MEX file which created a link between our MATLAB XFEM code [ 134 ] and CHOLMOD [ 190 ] The proposed reanalysis algorithm for crack growth can be summarized as follows: 1. Calculate the superset matrix for the initial quasi static iteration. 2. Find the fill reducing permutation of with AMD 3. C alculate and from the initial using the fill reducing ordering permutation. 4. Solve for the traditional and enriched degrees of freedom. 5. C alculate and Grow the crack. 6. C alculate the new portions of from prior crack growth 7. Use row delete to remove prior crack tip components from and 8. Use row add to add new Heaviside and crack tip components to and 9. Repeat steps 3 7. For the optimization example, no c rack growth occurs during step 5, and all Heaviside and crack tip components of are removed from and using row delete operations during step 7. While not shown here, this m ethod could easily be applied to the evolution of a void or inclusion as was considered by Mo s and Sukumar [ 55 113 ] The use of row add and row delete operations could also be used to consider cases where boundary conditions are non constant between iterations of crack growth. PAGE 115 115 Example Problems Reanalysis of an Edge Crack in a Finite Plate The first step in evaluating the proposed reanalysis algorithm is to consider a simple geometry that will allow us to both assess the savings from the proposed algorithm as we ll as to study the effect that increasing or decreasing the mesh density locally around the crack tip has on the proposed reanalysis algorithm. To perform this study, an edge crack of initial length 0.25 m in a finite plate of size 2 m x 2 m under unit uni axial tension with thirty increments of 0.05 m crack growth was considered as shown in Figure 6 2 The material properties were elastic modulus E = 1 Pa and Poisson's ratio v = 0.3. In order to study the effect of repeated iterati ons of crack growth a structured mesh of square plane strain elements with characteristic length 0.05 m was used, which was sufficient to achieve convergence in the Mode I and Mode II stress int ensity factors Figure 6 2 Edge crack in a finite plate geometry used to test quasi static crack growth with the reanalysis algorithm. H W a i PAGE 116 1 16 The savings is both the stiffness matrix assembly time as well as the time to factoring the global stiffness matrix or updat ing an existing factorization and then solving the resulting system of equations. The normalized times were calculated as the reanalysis algorithm divided by the traditional approach as a function of crack growth iteration as shown in Figure 6 3 The traditional approach is considered to be calculating or factoring and solving at each iteration from scratch. The comparison between the native MATLAB function backslash and CHOLMOD for factoring and solvin g the system of equations is considered reasonable as they both were written by the same author. The times were normalized in order to facilitate easy viewing of the savings that the proposed reanalysis algorithm yields. The proposed reanalysis reduces the assembly time by approximately 80% and factorization time by approximately 70%. It can be noticed from Figure 6 3 that the general trend for assembly is increased savings with additional iterations. This can be explained as havin g a larger crack increases the bandwidth of while only calculating the new enriched elements is approximately constant for the case of a fixed The increased fill in as the number of enriched element increases subtly increases the cost of updating the existing Cholesky factorization as the number of crack growth iterations increases. The second consideration was the effect that increasing or decreasing mesh density would have on the assembly and facto rization/solving savings. Two iterations of the previous problem were repeated for a range of mesh densities. The number of elements per unit length is plotted against the normalized time in Figure 6 4 where the time from the firs t iteration where the Cholesky update is directly modified is normalized by the time from the iteration where the Cholesky factorization is first performed. Note PAGE 117 117 that both trends show increased savings with increasing mesh density (i.e. number of degrees o f freedom). A B Figure 6 3 Normalized computational time for stiffness matrix assembly as well as factorization and solving of system of linear equations as a function of number of iterations. A) Stiffness matrix assembly, B) Factorization and solving of system of linear equations. A B Figure 6 4 Normalized computational time for stiffness matrix assembly as well as factorization and solving of system of linear equations as a function of number o f mesh density. A) Stiffness matrix assembly, B) Factorization and solving of system of linear equations. Optimization for Finding Crack Initiation in a Plate with a Hole To assess the applicability of the proposed reanalysis algorithm in an optimization f ramework, a simple optimization problem is considered as a second example. A finite plate of size 4 m x 4 m with a hole of radius 0.5 m centered at the middle of the plate is PAGE 118 118 considered. The material properties were repeated from the previous example of E = 1 Pa and v = 0.3. The mesh was a structured mesh of square elements with characteristic length of 0.05 m. The plate is subjected to uniaxial tension. The objective is to find the crack initiation point on the edge of a hole. A crack of length 0.15 m eman ates from the hole at an angle as shown in Figure 6 5 Figure 6 5 Edge crack in a finite plate geometry used to test optimization with the rean alysis algorithm. The optimization problem is solved using the fminbnd function native to MATLAB for 6 10 where for plane stress and for plane strain. As the optimization algorithm changes the angle theta, the factorized stiffness matrix is modified by deleting H W a i PAGE 119 119 the cra ck in the previous iteration and adding a crack corresponding to the new angle. The solution converges after 25 iterations to an angle of = 9.92 10 8 radians as the angle yielding the maximum energy release rate, which is in go od agreement with the expected solution of = 0 radians given the tolerances on the function and were 1 0 9 and 1 0 6 Figure 6 6 shows that for the optimization problem both the assembly and factorization/solving trends are comparable to the results from the quasi static problem. The assembly time no longer decreases as the iteration increases as no additional degrees of freedom are being introduced into the syst em of equations. The factorization and solving time is also largely constant as the matrix size is not changing as no crack growth is occurring. A B Figure 6 6 Normalized computational time for stiffness matrix assembly as we ll as factorization and solving of system of linear equations as a function of number of iterations. A) Stiffness matrix assembly, B) Factorization and solving of system of linear equations. Summary An exact reanalysis algorithm is presented based on the d irect modification of a sparse Cholesky factorization through the use of row modification operations. It has been shown that there is a significantly reduced computational cost for each iteration of PAGE 120 120 crack growth based on locally updating a Cholesky factori zation instead of simply factoring and solving the resulting system of linear equations. A study of the effect on mesh density was performed and it was found that while savings can be achieved over the full range of mesh densities ; the largest savings were achieved when the mesh had an increased level of refinement. The reanalysis algorithm allows for one to consider a smaller crack growth increment, which should lead to a better prediction of the crack growth path within a comparable computational budget. The use of the reanalysis algorithm in an optimization framework was also introduced and shown to reduce the cost of solving optimization problems involving cracks. PAGE 121 121 CHAPTER 7 VARIABLE AMPLITUDE MULTI AXIAL FATIGUE ANALYSIS OF WING PANEL Scaling of Normalized Flight Data from the Air Force Research Laboratory The Air Force Research Labo ratory (AFRL) at Wright Patterson Air Force Base in Ohio has provided normalized flight data for 19 different aircraft flights The normalized flight data includes: normal acceleration, lateral acceleration, longitudinal acceleration, roll acceleration, pi tch acceleration, yaw acceleration, roll rate, pitch rate, yaw rate, airspeed, altitude, angle of attack, flap angle, fuel quantity, and Mach number. The orientations related to the various acceleration and rate variables are given in Figure 7 1 For the finite element model, the x axis is located at the root of the wing at the midpoint of the chord length. The combined flight histories from the 19 flights was chosen to be the data set for this example problem. The flight data is u sed in conjunction with an Abaqus finite element model of an airplane wing box in order to calculate the stress history for a wing panel. At this point, the stress is calculated assuming there is no crack in the panel (nominal stress). This nominal stress history is then used to calculate stress intensity factors at the crack tip of the plate using the MXFEM code where the advantages of the exact reanalysis algorithm are used to solve the problem of an airplane wing panel subjected to variable amplitude mu lti axial loading. In order to simplify the analysis of a wing model, only a subset of the available normalized flight data have been used as a part of this work The parameters which were identified for use in the predictive model are: normal acceleration roll acceleration, airspeed, Mach number angle of attack, and fuel quantity. As the values provided by the AFRL are normalized, they must be scaled according to our estimation of the PAGE 122 122 expected range of values that the chosen parameters would take if they were not normalized The parameters were scaled based on the assumption that the data corresponded to a commercial aircraft flight instead of that of a military aircraft. The conversion the data from normalized to scaled values is summarized in Table 7 1 for Flight ID 152. The conversion for all flights is presented in Appendix C Figure 7 1 Coordinate system ( x y z ) and directional names for an airplane. The normal accele ration is related to the gravity forces (g force) acting upon the aircraft. Aerobatic and military aircraft can experience maximum g force of 9 12 g [ 197 ] Level flight would expect to experience forces equal to about 1 g. A typical banked turn with bank angle of 30 may cause forces of about 1.5 g [ 19 8 ] It is assumed that the range of g force experienced during a commercial aircraft flight is from a lower value of 0.75 g and a maximum value of 1.5 g. The normalized roll acceleration is used to approximate the bank angle associated with a banked turn. From a simple physics model it can be noticed that the z, normal x longitudinal y lateral pitch yaw roll PAGE 123 123 roll acceleration is directly related to the bank radius [ 199 ] For a larger roll acceleration the bank radius will decrease, while for a smaller roll acceleratio n the bank radius will increase. The banking radius is also inversely proportional to the bank angle. Therefore, the general trend can be established that the roll acceleration and bank angle are proportional. It is assumed that a commercial aircraft will not exceed the bank angle of 30 with a minimum bank angle of 0 during level flight. The normalized fuel capacity is related to the mass of fuel by considering the fuel capacity a nd mass of fuel of an aircraft. The Boeing 767 has a maximum fuel capacity of 63,000 L [ 200 ] There are two main types of commercial aircraft fuel in use in the world [ 201 ] In the United States, the aircraft fuel is commonly known as Jet A. In the remainder of the world a different specification is used which is known as Jet A 1. Additional additives are commonly used to convert commercial aircraft fuel into military fuel. The density of Jet A and Jet A 1 are very similar and depending upon the particular mixture will have density between 0.775 to 0.840 kg/L. Therefore, the expected fuel mass for an aircraft would range between 48,825 and 52,920 kg for a fully fueled aircraft. The assumed mass of a fully fueled aircraft is chosen here to be 5 0,000 kg. Furthermore, it is assumed that upon landing an aircraft will have at least 1000 L of fuel in reserve, which corresponds to a minimum fuel mass of 800 kg. The fuel mass is combined with the empty mass of the aircraft, 86,070 kg for a Boeing 767, to yield the total mass of the aircraft and fuel at each AFRL data point. The normalized Mach number M norm is related to the true mach number M by scaling the maximum and minimum values to 0 and 0.8 Mach [ 169 ] The normalized airspeed is related to the true airspeed of the aircraft by simply scaling the PAGE 124 124 minimum values of normalized airspeed to be equal to 0 m/s and the maximum value to 236 m/s [ 169 ] Similar ly, the normalized angle of attack is scaled to such that the maximum and minimum normalized values correspond to 0 and 10 [ 202 ] Table 7 1 Maximum and minimum values of normalized and scaled flight da ta with conversion relationship for Flight ID 152 Variable Norm Min Norm Max Scaled Min Scaled Max Conversion g 0.208 0.957 0.75 1.5 g = 0.644 g norm +0.884 B 1.00 0.624 0 30 B = 18.5 a roll +18.48 m fuel 0.169 0.975 800 50 000 m fuel = 51300 f norm M 0.0191 0.9 77 0 0.8 0 M = 0.835M norm 0.0159 v 0.0222 0.839 0 851 v = 289 v norm 6.42 0.00 0.853 0 10 = 11.7 norm Finite Element Wing Model For the stress anal ysis of an airplane wing a finite element model of a wing box is modeled. A wing box is the main structural support for a wing and supports most of the loading that occurs on an airplane wing. An airfoil is then attached to the wing box in order to produce the desired aerodynamic force. In this work, a finite element model of the wing box geometry originally used as part of a probabilistic optimization problem was used. In this previous work, an elliptical distribution was used to determine the span wise pr essure distribution, but the chord wise pressure distribution was neglected. A figure of a wing geometry showing the chord and span lengths is given in Figur e 7 2 and a figure of the wing in the chosen coordinate s ystem is given in Figure 7 3 The chord wise pressure distribution caused by the lift force is introduced a long with the lift induced drag and the magnitude of the resulting pressure s applied to the wing are given in terms of the chosen scaled variables from the normalized AFRL data. For the purposes of the finite element simulation the effects of lift and drag on the wing are considered to act independently and are separated into unique analysis. PAGE 125 125 Figur e 7 2 Example of a wing box with span s chord length c which is a function of y and normalized x coordinate u Figure 7 3 Wing box model in the coordinate system that was used for this analysi s. Wing Pressure Distribution Pippy [ 203 ] presented a simplified analysis of the normal pressure load caused by normal acceleration through the use of an elliptical loading model based upon the work of Chklovshi [ 204 ] and Lobert [ 205 ] I t has been shown that an elliptical loading model approximates the span wi se lift distribution for many taper ratio s as shown in Figure 7 4 The elliptical loa ding distribution gives the lift on a wing span as u c s x y PAGE 126 126 7 1 where is the lateral coordinate and is the maximum pressure applied to the wing. Note that the wing span here is the length of a si ngle wing and not the total wing span of a given airplane The wing span is assumed to take a value equal t o 18.9 m which corresponds to an aircraft similar to a Boeing 767. Equation 7 1 can be rewritten in terms of the lift for any span wise coordinate as 7 2 The maximum pressure applied to the wing can be found from the lift force which is acting upon the wing. The normalized span wise pressure distribution with range from 0 to 1 over the chosen airplane wing model is given by 7 3 Figure 7 4 Comparison of elliptical loading model for various wing taper ratios. Normalized Span Location Normalized Lift 1.6 0.2 0.1 1.0 Elliptical l oading Taper ratio, 1.0 Taper ratio, 0.5 Taper ratio, 0.0 PAGE 127 127 An estimate of the chord wise pressure distribution i s created based on some theoretical models as shown in Figure 7 5 is made as shown in Figure 7 6 A sixth order polynomial is used to represent the normalized chord wise pressure distribution with range f rom 0 to 1 based upon the normalized coordinate which has range from 0 to 1. An illustration of a the resulting polynomials is given in Figure 7 7 for several angles of attack, which modifies the point of maximum pressure The approach taken here will be to do finite element analysis for several angles of attack, and then to fit a surrogate model relating angle of attack to stress components at a point of interest. An expressi on for the normalized sixth order polynomial is given as 7 4 where the coefficients for various angles of attack are given in Table 7 2 The normalized coordinate is found as 7 5 where and are the coordinates of the leading and trailing edge of the airplane wing as a fu nction of the chord wise coordinate For the coordinate system which was chosen for the finite element model the leading edge follow s the relationship 7 6 and the trailing edge is given as 7 7 PAGE 128 128 Figure 7 5 Examples of the changing pressure distribution over an airfoil as a function of the angle of attack. A) 10 B) 5 C) 0 D) 5 Figure 7 6 Normalized chord wise coordinates u with maximum pressure coordinate u 0 pressure distribution P Table 7 2 The polynomial coefficie nts for the chord wise pressure distribution. Angle of Attack a 0 a 1 a 2 a 3 a 4 a 5 a 6 5 0.0258 7.73 20.0 22.0 11.5 1.83 0.00 0 0.0329 9.60 32.0 49.0 37.6 11.1 0.00 5 0.0253 11.7 47.0 58.7 75.4 25.1 0.00 10 0.0307 16.0 90.1 245 351 252 71 .7 1 0 u 0 1 P u A ngle of attack = 0 Angle of attack = 5 Angle of attack = 5 Angle of attack = 10 A C B D PAGE 129 129 The total pressure distribution for some point ( ) along the airplane wing is 7 8 where is the maximum value of the pressure distribution, is the chord wise pressure distribution given by Eq. 7 4 and is the span wise pressure distribution given by Eq. 7 3 The pressure distribution for = 1 N/m 2 is given in Figure 7 8 Figure 7 7 Comparison of the changes in the maximum pressure location u 0 and subsequently the resulting pressure profile for some angle s of attack Figure 7 8 The total normalized distribution along the airplane wing for = 5 1 0 1 P u 0.15 0.30 = 10 5 PAGE 130 130 The maximum value of the pressure distribution is found based upon a simple mode l of the lift force From the previous conversion of the roll acceleration to a bank angle it is possible to find a corresponding bank radius which is given as 7 9 where is the scaled airspeed and is the scaled normal acceleration from the AFRL data. The current mass of the airplane is calculated as 7 10 where is the empty mass of the aircraft and is the current fuel mass of the aircraft. The lift force [ 206 ] is now given as 7 11 The angle of att ack of the wing has a direct effect on the lift coefficient as shown in Figure 7 9 [ 201 ] Instead of directly calculating the lift from the lift coefficient, here it was decided to scale the lift to account for the changing angle of attack is g iven by 7 12 where for a given angle of attack is found from Figure 7 9 The lift calculated as 7 13 where is the total area under the normalized pressure distribution specified in Eq. 7 8 The area of the normalized pressure distribution was found using exact integration for the chord wise pressure distribution and the rectangle method of integration was used PAGE 131 131 in the span wise direction. The span wise step siz e was reduced until convergence had been achieved. Using the rectangle method 7 14 where The corresponding area for the normalized pressure distribution as a function of angle of attack are summarized in Table 7 3 where convergence is shown. Figure 7 9 Comparison between the angle of attack and the lift coefficient. Drag acts perpendicular to the span of the aircraft wing as shown in Figure 7 10 The drag on t he airplane wing is calculated through the use of the drag coefficient which is approximated as a linear function of the Mach number from 0 to 0.9 Mach [ 207 ] as 7 15 T his allows for the drag for a given wing to be calculated from the lift as 7 16 5 25 Lift Coefficient C L 0.2 1.8 5 15 1.0 PAGE 132 132 and subsequently converted into a uniform pressure with magnitude by dividin g by the surface area of the edge where the drag acts as 7 17 Equation 7 16 is used as the lift and drag can be calculated from the same equation by interchanging the lift and drag coefficient [ 208 ] Table 7 3 Area under the normalized pressure distribution for some angles of attack. Angle of Attack Number of Integration Steps, N Pressure Distribution Area, 10 1000 37.08 m 2 2000 37.06 m 2 5000 37.05 m 2 10000 37.05 m 2 5 1000 39.24 m 2 2000 39.22 m 2 5000 39.21 m 2 10000 39.21 m 2 0 1000 41.28 m 2 2000 41.26 m 2 5000 41.25 m 2 10000 41.25 m 2 5 1000 43.14 m 2 2000 43.12 m 2 5000 43.11 m 2 10000 43.11 m 2 Figure 7 10 The direction of the uniform pressure distribution defined by the drag. (Note that the drag is independent of the angle of attack, unlike the lift model) x y D rag PAGE 133 133 Finite Element Stress Analysis of Airplane Wing L inear elastic finite element analysis was performed on the airplane wing box geometry which was originally created by Pippy [ 197 ] within Abaqus The boundary conditions were created by connecting the edge of the wing box at x = 0 to th e center of gravity of the wing (0.43, 0.86, 0) by rigid elements and fixing all degrees of freedom at the reference point at the center of gravity. Five l oading conditions were considered here, the chord wise lift distribution with angles of attack of 5 0 5 and 10 as well as the drag distribution. All cases were considered with unit magnitude pressure. The two dimensional stresses were taken from ea ch analysis at the point ( 4.32, 7.42, 0) where the crack is located in the panel These values are given in Table 7 4 and were used for the stress calculation for each of the 19 unique flight profiles provided by the normalized d ata from the AFRL. Table 7 4 The stress components for various angles of attack used in the lift and drag pressure distributions. Analysis Dra g 22.3 33.4 32.2 Lift, = 5 135 481 300 Lift, = 0 134 483 398 Lift, = 5 133 484 296 Lift, = 10 131 480 291 Surrogate Models for Lift Coefficient, Stress Components, and Surface Area of Pressure Distribution Several surrogate models are used in the stress analysis of the airplane wing box. A surrogate model is used to allow for the lift coefficient to be evaluated as a function of the an gle of attack which follows Figure 7 9 Another surrogate model is used to fit the surface area of the pressure distribution given by Table 7 3 as a function of the angle of at tack For each case, a kriging surrogate as described in Chapter 5 is used. PAGE 134 134 A kriging surrogate model is fit to each stress component as a function of angle of attack for the chord wise lift distribution in order to allow for th e stress to be evaluated for any value of the angle of attack from the values in Table 7 4 The magnitude of the lift and drag distributions are given by in Eq. 7 13 and in Eq. 7 17 Therefore, the stress for a given set of scaled fli ght data can be calculated as 7 18 where and are the stress components from the kriging surrogate model for the lift pressure distribution and and are the stress components from the drag pressure distribution. The resulting stresses at the crack location and the stresses normalized by are given in Fi gure 7 11 for Flight ID 152. Note that the normalized s tresses are not constant in Fi gure 7 11 This presents a challenge in the fatigue analysis of this structure as the crack growth direction will change with each data point. This means no increment of crack growth across many cycle s can be used. The stresses and normalized stresses for all 19 flights are given in Appendix D A flowchart of the method for the conversion of the normalized flight data into biaxial stress history is summarized as: 1. Perform 5 finite element analysis; four for the chord wise pressure distribution as a function of angle of attack as shown in Figure 7 7 and Figure 7 8 ,and one for the drag pressure distribution as shown in Figure 7 10 to find the biaxial stresses at a point of interest of the wing box. 2. For each normalized variable find the conversion of the maximum and minimum normalized data point to the maximum and minimum scaled value as given in Table 7 1 Then find and use the conversion to create scaled data points from the normalized data points. PAGE 135 135 3. Create surrogate model s for the lift coefficient, area under the chord wise pressure distribution, and the biaxial stresses as a function of the angle of att ack. 4. Follow Eq. 7 9 to Eq. 7 13 for the calculation of w o an d Eq. 7 15 to Eq. 7 17 for the calculation of q o This will use the surrogate models for the lift coefficient and area under the chord wise pressure distribution. 5. Calculate the biaxial stresses from the values of w o q o and the surrogate models for the biaxial stresses as a function of the angle of attack. A B Fi gure 7 11 The stress histories for the two dimensional stress components for Flight ID 152. A) Bi axial stress components, B) Normalized stress components Analysis of Cracked Plate Subjected to Variable Amplitude Multi axial Load History D ue to the fact that the stresses in this problem are non proportional as seen in Fi gure 7 11 as a function of the history (e.g. stress ratios are non constant) this problem creates challenges for approximating the fa tigue crack growth of the problem. As the stress history is non proportional, the ratio of to is non constant, which implies that the angle of crack growth will change as a function of the loading his tory. As the path of crack growth is unknown, it is unclear how to select an analytical relationship between the crack size and stress intensity factor which will correctly consider the unknown crack path. Therefore, it is not possible to use a finite crac k growth increment as the crack growth direction cannot be assumed to be constant over a range of PAGE 136 136 elapsed cycles Therefore, the only choice is to choose the case where the cyclic history is followed Before the fatigue analysis can be started, first the biaxial stresses for each of the 19 flights must be converted into equivalent cyclic stresses. This is done according to the method outlined in Chapter 4 where an equivalent stress is calculated from the biaxial stress components. The rainflow counting algorithm is used to find the cycles from the equivalent stress history. The identified stress history is then superimposed back onto the biaxial stress data, resulting in cyclic biaxial stresses which c an be applie d to a plate and be used with a fatigue crack growth model. An example of the equivalent stress for Flight ID 152 and the identified equivalent cyclic loading history is shown in Figure 7 12 The results of the convers ion from normalized data points to stress cycles for each of the 19 flights is given in Table 7 5 There is an uneven number of data points can be explained as the rainflow counting method being interested in the extrema for the g iven data points. Not all data points which were provided by the AFRL correspond to an extrema. The conversion from biaxial stress components to an equivalent single stress can also result a suppression of the number of extrema present in the resulting rai nflow counting algorithm. Two function evaluations are required for each identified cycle as the ratio of the biaxial stresses changes between the maximum and minimum data point for each cycle, and the calculation of the mixed mode stress intensity factors without two function evaluations is not possible with sufficient accuracy. The initial geometry was that of an edge crack of initial crack size of 5 mm in a square finite plate with sides of length 0.1 m as shown in Figure 7 13 Note that a PAGE 137 137 geometry with does not necessarily agree with that of an actual wingbox panel at the specified location is used here for demonstration purposes.. The material was assumed to be aluminum 7075 T6 where the material properties and model constants are: threshold stress intensity factor of 2.2 Mode II threshold stress intensity factor of 1.0 Paris model constant of 6.85 10 11 Paris model exponent of 3.21, shaping parameters and of 0.30, 0.70, and 0.84 respectively [ 27 133 ] To encourage consistent, small scale crack growth the stresses calculated from the analysis of the wing box were scaled linearly by a factor of 10, which resulted in crack growth between 10 9 and 10 5 cycles for most of the cycles. In order to make the solution of this problem more affordable, the proposed exact reanalysis algorithm of Chapter 6 is used for the 74,014 required function evaluations. A B Figure 7 12 Equivalent stress and ident ified cycles from equivalent stress for Flight ID 152 A) Equivalent stress, B) Identified cycles. A study was performed on the initial edge crack geometry to ensure that the mixed mode stress intensity factors converged according the mesh density. This me sh was found to be that of a structured mesh of square elements with sides of length 1/200 m. The average computer time for each function evaluation was approximately 2.3 PAGE 138 138 seconds with the proposed exact reanalysis algorithm, which corresponds to about 2 da ys for a single analysis. As the computational cost without reanalysis is too high for a complete analysis, 100 iterations without the reanalysis algorithm were performed such that an estimate of the total computer time without the reanalysis algorithm cou ld be approximated for the purposes of comparison. It was found that the mean computer time for a cycle for these 100 iterations without reanalysis was about 20 seconds, which would correspond to about 17 days for a proposed fatigue analysis. Note that as crack growth occurs and more degrees of freedom are added to the system of equations, the cost of assembly, factorization, and solving will increase. Due to this, it is very possible that the true computer time required would be in excess of the predicted time of 17 days. Table 7 5 Summary of the AFRL data, the number of cycles identified, and the number of function evaluations for each Flight ID. Flight ID Number of Data Point Number of Identified Cycles Function Evaluations 2 11,158 1,780 3,560 8 12,893 2,207 4,414 10 6,050 1,409 2,818 12 9,628 1,838 3,676 16 7,544 1,311 2,622 18 7,422 1,126 2,252 22 4,057 797 1,594 26 10,124 1,895 3,790 138 6,038 1,270 2,540 140 10,427 2,275 4,550 146 10,431 2,266 4,532 148 9,766 2 ,290 4,580 150 16,038 3,415 6,830 152 15,120 3,470 6,940 154 11,044 2,754 5,508 156 5,631 1,025 2,050 159 7,616 1,669 3,338 161 9,623 2,180 4,360 163 9,978 2,030 4,060 Total 180,588 37,007 74,014 PAGE 139 139 Figure 7 13 The boundary conditions and applied bi axial stresses for the analysis of a panel on an airplane wing box subjected to variable amplitude fatigue. The initial and final crack paths from the analysis are given in Figure 7 14 For a more detailed view of the crack propagation path, refer to Figure 7 15 Note that initially, the crack tip begins to move down and to the right. This agrees with the observation that the largest stress compon ent is followed by With the loading specified in Figure 7 13 the resultant stress would be approximately perpendicular to the initial direction of crack growth. For the in itial crack has no effect on either the stress intensity factors or the crack propagation path. The subsequent crack propagation path represents a more difficult challenge to analyze in detail as the changing orientation of the crack tip represents a changing crack tip coordinate system and resultant stress acting at the crack tip in the Mode I and II. The general trend is to move in a horizontal motion, which is consistent with being the largest appli ed stress. The crack length as a function of the cycle number is given in Figure 7 16 where the grey lines denote the bounds associated with each unique flight One of the xx xx yy x y a W PAGE 140 140 observations from these plot s is that most of the crack gr owth occurs early in the life of the wing panel. In fact, of the 37,007 loading cycles, crack growth occurs in only about 2,500 of these cycles. This can be explained partially with the assistance of Figure 7 17 which gives the M ode I and II stress intensity factor ranges as a function of cycle number. Recall, that the threshold Mode I and II stress intensity factor ranges for aluminum were identified as taking values of 2.2 and 1.0 From Figure 7 17 it is clear that rarely does the Mode I stress intensity factor range exceed this value. Furthermore, it can be noticed that early in the panel s life, the Mode I stress intensity factor range tends to ta ke larger values than later in the panels life. Later in the lift of the panel, when the Mode II stress intensity range takes larger values, the Mode I values are at a minimum, reducing the net effect of increased It should als o be recalled that the predicted stresses at the crack location were increased by a factor of real application it is unlikely that this level of crack growth would occur at such an accelerated rate. A B Figure 7 14 The initial and final crack geometries from the stresses identified through the normalized flight data from the AFRL. A) Initial crack, B) final crack PAGE 141 141 Figure 7 15 Close up view of the final crack geometry from the stresses identified through the normalized flight data from the AFRL Figure 7 16 Comparison of the crack length as a function of t he cycle number for the 19 flights from the normalized AFRL data. Another cause for the reduce d amount of crack growth later in the panels life may be considered as a function of the and terms in the modified Paris model used in the analysis. The stress ratio R will increase the value of and therefore, accelerate crack growth. It can easily be seen in Figure 7 17 that the value of R will be PAGE 142 142 take a maximum value early in the life of the panel. An observation about what may be perceived as jumps in the crack length around cycle numbers 3500, 4000, 6000, and 9000. Each of these instan ces has a corresponding peak in the curve This peak will act as an overload, increasing the plasticity at the crack tip and slowing crack growth until the point in time when a load is applied which results in the creation of a larger plastic zone. As most of the peaks in the curve exist early in the life of the panel, the crack growth later in the panel's life is retarded as there are few instances where the applied load to exceed the previously developed crack tip plasticity. A B Figure 7 17 The Mode I and Mode II stress intensity factor ranges for the 37 007 cycles. A) Mode I stress intensity factor range, B) Mode II stress intensity factor range. Summary The fatigue analysis of an airplane wing subjected to mixed mode variable amplitude loading was conducted through the use of a simple stress model of an airplane wing box in coordination with normalized flight data provided by the AFRL. First, the normalized flight data was scaled to approximate values related to a commercial aircraft fl ight and a Boeing 767 in particular. Then, the two dimensional stress history from the 180,588 normalized data points was calculated at a particular PAGE 143 143 point of interest along an airplane wing. Next, an equivalent stress was calculated and the rainflow counti ng algorithm was used in order to convert the stress data into cyclic loading data consisting of 37,007 cycles which could be used in a fatigue crack growth model. These cycles were then superimposed onto the and stress histories. Due to the non proportional loading represented by the stress histories, fatigue modeling techniques which assume multiple elapsed cycles for each simulation cannot be used to solve this problem. Instead, each loading cycle was modeled using the exact XFEM reanalysis method based upon the direct modification of an existing Cholesky factorization. It was estimated that the combination of savings in assembly and factorization time for this problem r esulted in the problem being solved in about 86% of the time that would be needed for the traditional solution procedure. At each cycle, the crack growth direction was calculated according to the maximum circumferential stress criterion The fatigue crack growth model which was used was a modified version of the classical Paris model which also considers the effects of the threshold stress intensity factor, stress ratio R and load interactions as a function of the changes in plasticity as a result of overl oads and underloads. The mixed mode stress intensity factors are converted into a single equivalent stress intensity factor for use in the modified Paris model through the use of the methods introduced by Liu. Over the life of the panel, crack growth occur red in about 2,500 of the 37,000 modeled cycles. However, this is not information which would have been available before the analysis was run as no information about the amount or direction of crack growth was known before modeling the response of the pane l to the biaxial stresses. PAGE 144 144 Crack growth in a limited number of cycles was caused by a combination of the stress intensity factor ranges rarely exceeding the threshold values and the crack tip plasticity model in the modified Paris model which retarded crac k growth. It was observed that the stress intensity factor history for the crack showed that early in the panel's life, greater maximum values of and were often realized. This resulted in an increased amount of crack growth early in the life of the panel. L ater in the panels life, crack growth occurred more slowly and in a more controlled fashion, most likely due to both the lower peak values of the stress intensity factor ranges as well as the crack t ip plasticity which had been developed earlier in the panels life. PAGE 145 145 CHAPTER 8 CONCL USION S AND FUTURE WORK Conclusion s Fatigue is the process by which materials fail due to repeated loading well below the levels that they would fail under static loading conditions. It is not uncommon for 10 4 to 10 8 cycles to be needed for failure to occur This presents a challenge for computational simulation as this number of cycles is not feasible for modeling in a finite element environment. With classical finite element methods the mesh corresponds to the domain of interest, which means that the mesh must be recreated locally around the crack tip each time crack growth occurs. Fatigue crack growth has been one of the leading causes of aircraft accidents throughout history and still causes incidents today. Better computational tools for the modeling of fatigue crack growth will help to reduce the risk of accidents from the initiation and propagation of a crack. The goal of this work is to create strategies such that fatigue crack growth can be numerically modeled with greatly reduced computational requi rements. The challenges associated with mesh generation in the classical finite element method are addressed through the use of the extended finite element method (XFEM). In the XFEM, discontinuities are represented independent of the finite element mesh. Even when the challenges associated with mesh generation can be avoided; there is still a large number of function evaluations required in modeling fatigue crack growth. The fundamental structure of XFEM can be exploited to further reduce the cost of the r epeated function evaluations necessary for modeling fatigue crack growth. When crack growth occurs in the XFEM, only a small portion of the finite element stiffness matrix changes as crack growth occurs. These small changes in the finite element stiffness PAGE 146 146 matrix are exploited through the direct modification of an existing Cholesky factorization through the use of row add and row delete operations. The approximate fill reducing ordering for this matrix which has changing, but unknown connectivity is found th rough the use of a dummy stiffness matrix. Surrogate models enable the use of higher order numerical methods in the integration of a fatigue crack growth model for both the magnitude and direction of crack growth with the use of a single expensive finite e lement simulation. Test cases show that the use of kriging to replace expensive function evaluations offers comparable accuracy to classical numerical techniques with significantly fewer expensive function evaluations. A variable integration step size algo rithm is also introduced which uses the difference between the surrogate prediction and the finite element simulations to dynamically change the integration step size to attempt to achieve some target accuracy for each simulation. In the case that the part icular iteration has accuracy less than the target error, then the in tegration step size is increased. Similarly, for an iteration where the accuracy is greater than the target error, the integration step size is decreased. Finally, the XFEM exact reanalys is algorithm is used to predict fatigue crack growth due to variable amplitude biaxial loading in an airplane wing box. Normalized data from the Air Force research laboratory was scaled and then used to create a stress history for the airplane wing in term s of the scaled data. The rainflow counting algorithm was used to convert the discrete stress history into a cyclic stress history which was used as the input for the fatigue analysis. This is a first step towards realizing fatigue PAGE 147 147 crack growth from servic e load histories for variable amplitude multi axial stress histories. Future Work Future Development of MXFEM Additional enrichment functions could be implemented for cracks or inclusions, giving the user more options for learning about the XFEM, of partic ular interest may be allowing only the Heaviside enrichment for a crack and neglecting the crack tip terms, as this comparison is often made in the literature Contact between crack faces could be implemented to prevent for the crack surfaces from entering the materials as would currently occur for compressive loading Modifying the definition of where boundary conditions and loading are applied instead of requiring the direct modification of files other than the input file for the given loading conditions may facilitate the ease of use for simulations Removing the limitation for square elements and a rectangular domain, perhaps with the option to import the nodal coordinates and element connectivity information from a text file which would allow for comple x geometries to be considered Possible Application for the Exact XFEM Reanalysis Algorithm An example of a crack growth problem which could benefit from an application of the proposed exact XFEM reanalysis algorithm is that of elastic crack growth between dissimilar materials. Consider the case of a crack which is propagating towards a material interface as shown in Figure 8 1 Upon reaching the interface the crack may either be arrested, penetrate into the interface, propagate al ong the interface or branch and propagate along the interface [ 209 ] The challenge in modeling this problem is that the crack propagation behavior is a function of the two materials which are considered as well as the angle at which the crack approaches the interface. There is no closed PAGE 148 148 form solution for how the crack will behave upon reaching the interface, and one possible solution procedure would involve the so lution of an optimization problem to find the path with the greatest energy release rate from the available possible propagation paths. Figure 8 1 Possible crack propagation paths between dissimila r materials A) Arrest, B) Penetration, C) Propagate along interface D) Branch along interface The proposed exact XFEM reanalysis algorithm could be used to create an optimization problem where the penetration or growth patterns along the material interf ace are considered. The angle of penetration which maximizes the energy release rate could be compared to the energy release rates for growth along the interface or crack branching at the interface with the maximum energy release rate among these options b eing the most likely crack growth path. This energy release rate can then be converted to a stress intensity value and compared with the fracture toughness of the A C B D Material 1 Material 2 Material 1 Material 2 Material 1 Material 2 Material 1 Material 2 PAGE 149 149 material interface or material into which the crack is penetrating. If the stress intensity i s not sufficiently large to cause fracture, then the case of a crack arresting at the interface may be observed. The major modifications which need to be considered for this problem with the MXFEM code include: implementation and verification of bi materia l [ 85 ] and branching [ 89 ] crack enrichment functions, the implementation of the auxiliary stress and displacement states for the bi material crack [ 85 ] and the formulation of the optimization problem to determine the direction of maximum energy release rate. The analysis may be v erified against the work of He [ 209 210 ] and Zhang [ 211 ] If this implementation was completed it would be possible to fu r ther expend this appr oach to in terface layers [ 212 ] and laminates [ 213 ] Possible Applications for Surrogate Integration An example of a problem which may benefit from the use of the proposed surrogate integration method is that of wear occurring in the joint of multibody systems due to cont act. The Archard wear model [ 214 ] is common ly used to predict the amount of wear such that 8 1 where h is the wea r depth, k is a wear constant, p is the contact pressure between the two bodies in contact, and s is the sliding distance between the two bodies. Note that this problem is quite similar in form to the classical Paris model where a correlation could be made where and Finite element simulations with a contact model are commonly used to calculate the contact pr essure and sliding distance [ 215 ] while the wear constant k usually is determined from experimental data. PAGE 150 150 Wear is a slow process which occurs over a large number of iterations, similar to fatigue. In this way, it wo uld be too expensive to perform a finite element simulation during each loading cycle to find p and s calculate the wear depth, update the mesh, and then perform the next analysis. Ins tead, a forward Euler approach [ 215 ] is often used to approximate the wear over a number of elapsed cycles according to a single finite element simulation for p and s as 8 2 This is the same basic problem which was considered in the case of modeling fatigue crack growth and it is possible that the use of kriging extrapolation and/or the variable step size approach could be applied to a wear analysis and result in a better approximation of the amount of total wear over a number of elapse d cycles with some modifications to account for the differences between the two problems. PAGE 151 151 APPENDIX A MXFEM BENCHMARK PROB LEMS XFEM Enrichments C enter Crack in a Finite Plate The first problem considered in the verification of the MATLAB XFEM code is that of a center crack in a finite plate. The geometry and specifications for the values of and is shown in Figure A 1 The material properties used in the analysis were .3. The full domain was a plate with height 10 m and width 6 m with a center crack of length 2 m. The applied stress was 1 Pa Square plane strain quadrilateral elements with a structured mesh were used. Both a full and half model were considered based on the symmetry in the problem. The theoretical stress intensity factor for a center crack in a finite plate [ 20 ] is A 1 where is the applied stress, and is the half crack length. A compariso n of the convergence as a function of the average element size h is given in Table A 1 where the normalized stress intensity factors are given as A 2 where is given by Eq. A 1 and is the va lue calculated by the XFEM analysis using the domain form of the contour integral. The results for the analysis of the full and half model are given in Table A 1 From this table it can be noticed that an accurate solution is achi eved for a mesh with characteristic element size of 1/20. There is limited PAGE 152 152 improvement with increased mesh density in the vicinity of the crack tip. It can be noticed that there is a difference between the full and half model created by symmetry. These cha nges are most likely due to the slightly different boundary conditions used for the two cases. In the case of a full center crack, the midplane where the crack is located is restricted from moving the vertical direction and the midpoint of the plate is fix ed. For the half model, the left edge are given rollers in the y direction and the bottom left hand corner is fixed. The XFEM model shows good accuracy without any of the cumbersome challenges associated with solving a similar problem in the traditional FE M framework. Figure A 1 Representative geometry for a center crack in a finite plate. Table A 1 Convergence of normalized stress intensity factors for a full and half m odel for a center crack in a finite plate. h Full, Right Tip Full, Left Tip Half, Right Tip 1/5 0.954 0.952 0.962 1/10 0.977 0.977 0.967 1/20 0.989 0.989 0.980 1/40 0.9 96 0.996 0.990 1/80 0.999 0.999 0.996 2a 2W PAGE 153 153 Edge Crack in a Finite Plate The next problem considered in the verification of the MATLAB XFEM code is that of an edge crack in a finite plate. The geometry and specifications for the values of and is shown in Figure A 2 The material properties used in the analysis were was a plate with height 6 m and width 3 m with a center crack of length 1 m. The applied stress was 1 Pa. Square plane strain quadrilateral elements with a structured mesh were used. Figure A 2 Representative geometry for an edg e crack in a finite plate The theoretical stress intensity factor for an edge crack in a finite plate [ 20 ] is A 3 where is the applied stress, and is the characteristic crack length. A com parison of the convergence as a function of the average element size h is given in Table A 2 where the normalized stress intensity factors are given by Eq. A 2 From this H W a PAGE 154 154 table it can be noticed that an accurate solution is achieved for a mesh with characteristic element size of 1/10. There is limited improvement with increased mesh density in the vicinity of the crack tip. This repres ents a decreased mesh density to be needed for a solution of comparable accuracy to the case of a center crack in a finite plate. This can be explained by the value of the Mode I stress intensity factor in this case taking a larger value than in the previo us example, which makes the normalization less sensitive to small differences between the theoretical and calculated values. Table A 2 Convergence of normalized stress intensity factors for an edge crack in a finite plate. h 1/5 0.884 1/10 0.991 1/20 1.001 1/40 1.002 1/80 1.002 Inclined Edge Crack in a Finite Plate The case of an inclined edge crack in a finite plate benchmark problem is used to show accuracy in mixed mode problems. A figure of the geometry is given in Figure A 3 2 m and width 1 m with an edge crac k from (0,1) to (0.4,1.4). The applied stress was 1 Pa. Square plane strain quadrilateral elements wi th a structured mesh were used. A comparison of the convergence as a function of the average element size h is given in Figure A 3 where the normalized stress intensity factors are given by Eq. A 2 The reference solution [ 216 ] which was used for comparison is = 1.927 and = 0.819. This reference solution is commonly used i n the XFEM literature [ 96 ] and is a solution with accepted accur acy. convergence of the Mode I and Mode II stress intensity PAGE 155 155 factors for this problem are comparable to that for the center crack in a finite plate problem. Figure A 3 Representative geometry for an inclined edge crack in a finite plate. Table A 3 Convergence of normalized stress intensity factors for an inclined edge crack in a finite plate. h 1/10 0.958 1.028 1/20 0.980 1.022 1/40 0.987 1.018 1/80 0.990 1.016 Hard Inclusion in a Finite Plate T he inclusion enrichment shown in Figure A 4 is calibrated based on results from the commercial finite element code Abaqus The mat erial properties used in the analysis were c G for The full domain was a plate with height 10 m and width 6 m with a center inclusion of W a H PAGE 156 156 radius 0. 5 m. Through symmetry only the top right hand quarter of the plate was modeled. The applied stress was 1 Pa. Square plane strain quadrilateral elements with a structured mesh with average element size h of 1/30 m were used for comparison against Abaqus The stress contours for the Abaqus and MXFEM solutions are presented in Figure A 5 Figure A 6 and Figure A 7 All contours have color bar s with matching ranges. Abaqus plots from black to grey, while MATLAB plots from black to white, which results is some slight differences between figures. Note that there is also a difference in how the interface is represented as Abaqus plots a black l ine, while the MXFEM code denotes all inclusions with a white line. Some of the differences in contour shape may be attributed to a different number of contours being plotted in the two programs. Figure A 4 Representative geometry for a hard inclusion in a finite plate. 2 H 2 W r PAGE 157 157 A B Figure A 5 The xx contours for a hard inclusion. A) Abaqus B) MXFEM. A B Figure A 6 The x y contours for a hard inclusion. A) Abaqus B) MXFEM. A B Figure A 7 The yy contours for a hard inclusion. A) Abaqus B) MXFEM. PAGE 158 158 Void in an Infinite Plate The void enrichment function is verified against the Abaqus solution for the trad itional FEM solution of a plate containing a hole. The material properties used in the full domain was a plate with height 3 m and width 3 m with a center hole of radius 0 3 m. The applied stress was 1 Pa Square plane strain quadrilateral elements with a structured mesh with average element size h of 1/20 m w as used. The stress contours based on the Abaqus solution are compared to the XFEM stress contours in Figure A 9 Figure A 10 and Figure A 11 All stress components show excellent agreement with the Abaqus values. Any differences can be attributed to either having a structured mesh for the XFEM solution compared to a conforming mesh for the Abaqus solution or the color scale in Abaqus going from black to grey, while it goes from black to white within MATLAB Figure A 8 Representat ive geometry for a void in an infinite plate. H W r PAGE 159 159 A B Figure A 9 The xx contours for a void in an infinite plate. A) Theoretical, B) MXFEM. A B Figure A 10 The x y contours for a void in an i nfinite plate. A) Abaqus B) MXFEM. A B Figure A 11 T he yy contours for a void in an infinite plate. A) Abaqus B) MXFEM. PAGE 160 160 Other Benchmarks Angle of Crack Initiation from Optimization An optimization algorithm is introdu ced using fminbnd in MATLAB to determine the angle of crack initiation which results in the maximum energy release rate. The example of a 8 m x 8 m plate with a hole with 1 m radius subjected to uniaxial tension as shown in Figure A 12 is used to show that the optimization implementation has been completed with accuracy Square plane strain quadrilateral elements with a structured mesh with an average element size of h equal to 1/20 m was used. For this problem the location of max imum stress corresponds t o the angle of crack initiation for a 0.15 m crack which will result in the maximum energy release rate. For a plate under uniaxial tension this angle is 0 or radians Figure A 12 Representative geometry for a crack initiating at an angle for a plate with a hole. The optimization problem used to benchmark the implementation of the XFEM optimization subroutine is given as H W a i PAGE 161 161 A 4 where is the energy release rate, and are the Mode I and II stress intensity factors, and are the lower and upper bounds on the angle and is given as A 5 The results are pr esented in Table A 4 for two sets of bounds and Resu lts are presented as the optimum angle identified by fminbnd with the function tolerance set to 1E 9 and the angle tolerance set to 1E 6. The normalized energy release rate is A 6 where is the energy release rate corresponding to an initiation angle of 0 or radians and is the energy release rate from the optimization with the XFEM For this example the theoretical maximum energy release rate is given as = 3. 5166 E 7 J/m 2 Based upon the function and angle tolerances, the solutions represent an excellent agreement between the known solution and the MXFEM optimization result. Table A 4 Comparison of the theoretical and MXFEM value for maximum energy release angle as a function of average element size. h 1/20 1.715E 5 1.001 3.142 1.000 PAGE 162 162 Crack Growth in Presence of an Inclusion Bordas [ 69 ] presented the case of edge crack growth under mixed mode loading in a plate with either a hard or soft inclusion. This problem acts to verify the implementation of mixed mode crack problems, the inclusion enrichment function, and the use of both the cr ack and inclusion enrichment functions in a single framework. The ratio of 0.3 for mat erial 2. For the case of a hard inclusion, material 2 is for the inclusion and material 1 is for the plate. For the case of a soft inclusion, material 1 is for the inclusion and material 2 is for the plate. The full domain was a plate with height 8 m and w idth 4 m with an edge crack of length 0.5 m centered on the left edge. The circular inclusion of radius 1 is placed at the horizontal midpoint and the vertical quarter point. T here were 34 iterations of growth where at each iteration the amount of growth w as equal to = 0.1 m and in the direction given by the maximum circumferential stress predicted by a unit tensile load which was applied to the top edge of the domain in the y direction Square plain strain quadrilateral elements with a structured mesh with average element size h of 1/20 m were used. It is important to note that no boundary conditions were given by Bordas so here it is assumed that edge crack boundary conditions are to be used. This may explain the slight differen ces between the results of Bordas and those presented in Figure A 13 Note that the scales of the x and y axis are different, which result in the circular soft inclusion appearing to take a shape more similar to that of an ellipse Also, what appears to be a large initial crack growth increment corresponds to the two data points corresponding to the initial crack geometry. PAGE 163 163 A B Figure A 13 Comparison of crack paths to those predicted by Bordas for a ha rd and soft inclusion. A) Hard inclusion, B) Soft inclusion. Fatigue Crack Growth The implementation of the Paris model [ 34 ] for fatigue crack growth was verified within the XFEM framework. The Paris model predicts that crack growth occurs according to the differenti al equation A 7 where is the crack growth rate, is the Pa ris model constant, is the Paris model exponent, and is the stress intensity factor range under cyclic loading. Here the geometry is chosen to mimic a center crack in an infinite plate where is A 8 where is the applied stress range and is the half crack length. 10, and Paris model exponent of 3.8. The full domain was a plate with height 0.4 m and width 0.6 m with a center crack of length 0.02 m. The applied stress was 78.6 MPa. Square plane strain quadrilateral elements with a structured mesh with average element size h of 1/500 m PAGE 164 164 was used. Through symmetry, only the right half of the geo metry was modeled. Fatigue crack growth was measure to 2450 cycles with step sizes of equal to 50 at each simulation. A comparison of the theoretical and simulated crack growth is given in Figure A 14 From this analysis it is apparent that the XFEM and Paris model solutions are in good agreement. As the crack size increases, there begins to be a difference between the two solutions. It is likely that the choice of a smaller number of elapsed cycles at e ach iteration would reduce the differences between the two predictions for large crack sizes. Figure A 14 Comparison between Paris model and XFEM simulation for fatigue crack growth for a center crack in an infinite plate PAGE 165 165 APPENDIX B AUXILIARY DISPLACEME NT AND STRESS STATES The auxiliary displacement and stress states given by Westergaard [ 145 ] and Williams [ 146 ] for that are used in the interaction integrals [ 16 18 19 ] to calculate the stress intensity factors are given in terms of polar coordinates from crack tip and shear modulus Poisson's rat io and Kosolov constant as 10 1 10 2 10 3 10 4 10 5 10 6 10 7 10 8 10 9 10 10 PAGE 166 166 APPENDIX C CONVERSION OF NORMAL IZED FLIGHT DATA TO SCALED FLIGHT DATA Table C 1 Normalized and scaled flight data with conversion for Flight ID 2. Variable Norm Min Norm Max Scaled Min Scaled Max Conversion g 0.0186 0.908 0.75 1.5 g = 0.843 g norm +0.734 B 0.711 0.574 0 30 B = 23.2 a roll +16.5 m fuel 0.00 0.971 800 50,000 m fuel = 51500 f norm M 0.0191 0.921 0 0.80 M = 0.887 M norm 0.0170 v 0.0222 0.822 0 236 v = 295 v norm 6.56 0.0178 0.853 0 .001 10 = 12.0 norm 0.212 Table C 2 Normalized and scaled flight data with conversion for Flight ID 8. Variable Norm Min Norm Max Scaled Min Scaled Max Conversion g 0.0609 0.914 0.75 1.5 g = 0.769 g norm +0.797 B 0.485 0.490 0 30 B = 30.8 a roll +14.9 m fuel 0.00 0.990 800 50,000 m fuel = 50500 f norm M 0.0191 0.954 0 0.80 M = 0.856 M norm 0.0163 v 0.0222 0.962 0 236 v = 252 v norm 5.59 0.00 0.833 0.001 10 = 1 2.0 norm + 0.001 Table C 3 Normalized and scaled flight data with conversion for Flight ID 10. Variabl e Norm Min Norm Max Scaled Min Scaled Max Conversion g 0.0793 0.337 0.75 1.5 g = 2.91 g norm +0.519 B 0.196 0.195 0 30 B = 76.7 a roll +15.0 m fuel 0.00 0.986 800 50,000 m fuel = 50700 f norm M 0.191 0.579 0 0.80 M = 1.43 M norm 0.0273 v 0.0222 0.663 0 236 v = 369 v norm 8.19 0.00 0.844 0.001 10 = 11.9 norm + 0.001 Table C 4 Normalized and scaled flight data with conversion for Flight ID 12. Variable Norm Min Norm Max Scaled Min Scaled Max Conversion g 0.0599 0.847 0.75 1.5 g = 0.952 g norm + 0.693 B 0.324 0.312 0 30 B = 47.2 a roll +15.3 m fuel 0.00 0.993 800 50,000 m fuel = 50 4 00 f norm M 0.0191 0.867 0 0.80 M = 0. 944 M norm 0.0 180 v 0.0222 0.828 0 236 v = 293 v norm 6.51 0.00 0.723 0.001 10 = 1 3.8 norm + 0.001 Table C 5 Normalized and scaled flight data with conversion for Flight ID 16. Variable Norm Min Norm Max Scaled Min Scaled Max Conversion g 0.0248 0.965 0.75 1.5 g = 0.797 g norm + 0.730 B 0.639 0.634 0 30 B = 31.2 a roll + 19.9 m fuel 0.00 0.996 80 0 50,000 m fuel = 50 2 00 f norm M 0.0191 0.987 0 0.80 M = 0. 827 M norm 0.01 58 v 0.0222 0.919 0 236 v = 2 64 v norm 5.85 0.00 0.794 0.001 10 = 1 2. 6 norm + 0.001 PAGE 167 167 Table C 6 Normalized and scaled flight data with conversion for Flight ID 18. Variable Norm Min Norm Max Scaled Min Scaled Max Conversion g 0.145 0.880 0.75 1.5 g = 0.732 g norm + 0.856 B 0.817 0.759 0 30 B = 19.0 a roll +15.5 m fuel 0.00 0.941 800 50,000 m fuel = 531 00 f norm M 0.0191 0.941 0 0.80 M = 0. 868 M norm 0.0 166 v 0.02 22 0.825 0 236 v = 295 v norm 6.54 0.00 0.721 0.001 10 = 1 3.9 norm + 0.001 Table C 7 Normalized and scaled flight data with conversion for Flight ID 22. Variable Norm Min Norm Max Scaled Min Scaled Max Conversion g 0.0838 0 .848 0.75 1.5 g = 0. 981 g norm +0. 668 B 0.488 0.505 0 30 B = 30.2 a roll + 14.7 m fuel 0.00 0.998 800 50,000 m fuel = 501 00 f norm M 0.0191 0.947 0 0.80 M = 0.86 2 M norm 0.016 5 v 0.0222 0.821 0 236 v = 2 96 v norm 6.57 0.00 0.819 0.001 10 = 1 2.2 norm + 0.001 T able C 8 Normalized and scaled flight data with conversion for Flight ID 26. Variable Norm Min Norm Max Scaled Min Scaled Max Conversion g 6.88E 4 0.586 0.75 1.5 g = 1.28 g norm + 0.749 B 0.463 0.291 0 30 B = 3 9 .8 a roll +18.4 m fuel 0.00 0.992 800 50,000 m fuel = 50 4 00 f norm M 0.0191 0.617 0 0.80 M = 1.34 M norm 0.0256 v 0.0222 0.709 0 236 v = 344 v norm 7.65 0.0178 0.858 0.001 10 = 11.9 norm + 0.001 Table C 9 Normalized and scaled flight data with c onversion for Flight ID 138. Variable Norm Min Norm Max Scaled Min Scaled Max Conversion g 0.0269 0.882 0.75 1.5 g = 0.825 g norm +0.772 B 0.569 0.651 0 30 B = 24.6 a roll +14.0 m fuel 0.00 0.913 800 50,000 m fuel = 54700 f norm M 0.0191 0.961 0 0.80 M = 0. 8 49 M norm 0.0162 v 0.0222 0.836 0 236 v = 290 v norm 6.45 0.00890 0.869 0.001 10 = 1 1.6 norm + 0.001 Table C 10 Normalized and scaled flight data with conversion for Flight ID 140. Variable Norm Min Norm Max Scaled Min Scaled Max Conversion g 0.00730 1.00 0.75 1.5 g = 0.756 g norm +0.744 B 0.607 0.756 0 30 B = 22.0 a roll +13.4 m fuel 0.00 0.934 800 50,000 m fuel = 53500 f norm M 0.0191 0.951 0 0.80 M = 0.886 M norm 0.0163 v 0.0222 0.838 0 236 v = 290 v norm 6.43 0.0124 0.799 0.0 01 10 = 1 2.7 norm 0.157 PAGE 168 168 Table C 11 Normalized and scaled flight data with conversion for Flight ID 146. Variable Norm Min Norm Max Scaled Min Scaled Max Conversion g 0.0775 0.913 0.75 1.5 g = 0.757 g norm +0.809 B 0.709 0 .727 0 30 B = 20.9 a roll +14.8 m fuel 0.00 0.914 800 50,000 m fuel = 54700 f norm M 0.0190 0.991 0 0.80 M = 0. 823 M norm 0.0156 v 0.0222 1.00 0 236 v = 242 v norm 5.37 0.00 0.879 0.001 10 = 1 1.3 norm + 0.001 Table C 12 Normalize d and scaled flight data with conversion for Flight ID 148. Variable Norm Min Norm Max Scaled Min Scaled Max Conversion g 0.0136 0.954 0.75 1.5 g = 0.775 g norm +0.761 B 0.775 0.611 0 30 B = 21.6 a roll +16.8 m fuel 0.00 0.961 800 50,000 m fuel = 52000 f norm M 0.019 1.00 0 0.80 M = 0. 815 M norm 0.0155 v 0.0222 0.983 0 236 v = 246 v norm 5.46 0.00 0.908 0.001 10 = 1 1.0 norm + 0.001 Table C 13 Normalized and scaled flight data with conversion for Flight ID 150. Variable Norm Min Norm Max Scaled Min Scaled Max Conversion g 0.169 0.919 0.75 1.5 g = 0.689 g norm +0.866 B 0.920 0.745 0 30 B = 18.0 a roll +16.6 m fuel 0.00 0.956 800 50,000 m fuel = 52300 f norm M 0.0190 0.954 0 0.80 M = 0.856 M norm 0.0163 v 0.0222 0.837 0 236 v = 290 v nor m 6.44 0.00 1.00 0.001 10 = 1 0.0 norm + 0.00 1 Table C 14 Normalized and scaled flight data with conversion for Flight ID 152. Variable Norm Min Norm Max Scaled Min Scaled Max Conversion g 0.208 0.957 0.75 1.5 g = 0.644 g n orm +0884 B 1.00 0.624 0 30 B = 18.5 a roll +18.5 m fuel 0.00 0.975 800 50,000 m fuel = 51300 f norm M 0.0191 0.977 0 0.80 M = 0. 835 M norm 0.0159 v 0.0222 0.839 0 236 v = 289 v norm 6.42 0.00 0.853 0.001 10 = 1 1.7 norm + 0.001 Table C 15 Normalized and scaled flight data with conversion for Flight ID 154. Variable Norm Min Norm Max Scaled Min Scaled Max Conversion g 0.011 0.922 0.75 1.5 g = 0.804 g norm +0.759 B 0.696 0.513 0 30 B = 24.8 a roll +17.3 m fuel 0.00 0.993 800 50,00 0 m fuel = 50400 f norm M 0.0191 0.979 0 0.80 M = 0.8 34 M norm 0.01 59 v 0.0222 0.978 0 236 v = 247 v norm 5.49 0.00 0.931 0.001 10 = 1 0.7 norm + 0.001 PAGE 169 169 Table C 16 Normalized and scaled flight data with conversion for Flight ID 156. Variable Norm Min Norm Max Scaled Min Scaled Max Conversion g 4.20E 4 0.900 0.75 1.5 g = 0. 833 g norm +0.750 B 0.652 0.697 0 30 B = 22.2 a roll +14.5 m fuel 0.00 0.927 800 50,000 m fuel = 5 39 00 f norm M 0.0219 0.928 0 0.80 M = 0.8 83 M norm 0.0 19 3 v 0.025 6 0.821 0 236 v = 2 97 v norm 7.61 0.0373 0.831 0.001 10 = 12.6 norm 0.469 Table C 17 Normalized and scaled flight data with conversion for Flight ID 159. Variable Norm Min Norm Max Scaled Min Scaled Max Conversion g 0.07 65 0.917 0.75 1.5 g = 0. 755 g norm +0.808 B 0.530 0.492 0 30 B = 29.3 a roll +15.6 m fuel 0.00 0.946 800 50,000 m fuel = 5 28 00 f norm M 0.0192 0.993 0 0.80 M = 0.8 21 M norm 0.01 58 v 0.0225 0.985 0 236 v = 2 46 v norm 5.53 0.00 0.734 0.001 10 = 1 3. 6 norm + 0.001 Table C 18 Normalized and scaled flight data with conversion for Flight ID 161. Variable Norm Min Norm Max Scaled Min Scaled Max Conversion g 0.0182 0.878 0.75 1.5 g = 0. 872 g norm +0.734 B 0.492 0.471 0 30 B = 3 1.2 a roll +1 5.3 m fuel 0.00 0.936 800 50,000 m fuel = 534 00 f norm M 0.0201 0.9808 0 0.80 M = 0.8 32 M norm 0.016 7 v 0.0236 0.966 0 236 v = 25 1 v norm 5. 92 0.00 0.826 0.001 10 = 1 2.1 norm + 0.001 Table C 19 Normalized and scaled flight data with conversion for Flight ID 163. Variable Norm Min Norm Max Scaled Min Scaled Max Conversion g 0.00360 0.900 0.75 1.5 g = 0. 837 g norm +0.747 B 0.518 0.725 0 30 B = 24.1 a roll +12.5 m fuel 0.00 0.946 800 50,000 m fuel = 529 00 f norm M 0.0209 0.946 0 0.80 M = 0. 865 M norm 0.01 81 v 0.0244 0.843 0 236 v = 2 89 v norm 7.05 0.00890 0.879 0.001 10 = 1 1.5 norm 0 .1 01 PAGE 170 170 APPENDIX D STRESS HISTORIES FOR FLIGHT DATA FROM AFR L Figure D 1 Stress history for Flight ID 2. A) Bi axial str ess, B) Normalized stress. Figure D 2 Stress history for Flight ID 8. A) Bi axial stress, B) Normalized stress. Figure D 3 Stress history for Flight ID 10. A) Bi axial stress, B) Normalized stress. PAGE 171 171 Figure D 4 Stress history for Flight ID 12. A) Bi axial stress, B) Normalized stress. Figure D 5 Stress history for Flight ID 16. A) Bi axial stress, B) Normalized stress. Figure D 6 Stress history for Flight ID 18. A) Bi axial stress, B) Normalized stress. PAGE 172 172 Figure D 7 Stress history for Flight ID 22. A) Bi axial stress, B) Normalized stress. Figure D 8 Stress history for Flight ID 26. A) Bi axial stress, B) Normalized stress. Figure D 9 Stress history for Flight ID 138. A) Bi axial stress, B) Normalized stress. PAGE 173 173 Figure D 10 Stress history for Flight ID 140. A) Bi axial stress, B) Normalized stress. Figure D 11 Stress history for Flight ID 146. A) Bi axial stress, B) Normalized stress. Figure D 12 Stress history for Flight ID 148. A) Bi axial stress, B) Normalized stress. PAGE 174 174 Figure D 13 Stress history for Flight ID 150. A) Bi axial stress, B) Normalized stress. Figure D 14 Stress history for Flight ID 152. A) Bi axial stress, B) Normalized stress. Figure D 15 Stress history for Flight ID 154. A) Bi axial stress, B) Normalized stress. PAGE 175 175 Figure D 16 Stress history for Flight ID 156. A) Bi ax ial stress, B) Normalized stress. Figure D 17 Stress history for Flight ID 159. A) Bi axial stress, B) Normalized stress. Figure D 18 Stress history for Flight ID 161. A) Bi axial stress, B) Normalized stress. PAGE 176 176 Figure D 19 Stress history for Flight ID 163. A) Bi axial stress, B) Normalized stress. PAGE 177 177 LIST OF REFERENCES 1. "de Havilland Comet." Wikipedia: The Free Encyclopedia Wikimedia Found ation, Inc. 25 Oct. 2011. Web. 25 Oct. 2011. 2. "1957 Cebu Douglas C 47 Crash." Wikipedia: The Free Encyclopedia Wikimedia Foundation, Inc. 17 Mar. 2011. Web. 25 Oct. 2011. 3. "Los Angeles Airways Flight 417." Wikipedia: The Free Encyclopedia Wikimedia F oundation, Inc. 8 Oct. 2011. Web. 25 Oct. 2011. 4. "Alexander L. Kielland (platform)." Wikipedia: The Free Encyclopedia Wikimedia Foundation, Inc. 25 Aug. 2011. Web. 25 Oct. 2011. 5. "Japan Airlines Flight 123." Wikipedia: The Free Encyclopedia Wikimedia Foundation, Inc. 22 Oct. 2011. Web. 25 Oct. 2011. 6. "Aloha Airlines Flight 243." Wikipedia: The Free Encyclopedia Wikimedia Foundation, Inc. 3 Oct. 2011. Web. 25 Oct. 2011. 7. "United Airlines Flight 232." Wikipedia: The Free Encyclopedia Wikimedia Fou ndation, Inc. 8 Oct. 2011. Web. 25 Oct. 2011. 8. "El Al Flight 1862." Wikipedia: The Free Encyclopedia Wikimedia Foundation, Inc. 23 Sep. 2011. Web. 25 Oct. 2011. 9. "Eschede Train Disaster." Wikipedia: The Free Encyclopedia Wikimedia Foundation, Inc. 1 Oct. 2011. Web. 25 Oct. 2011. 10. "China Airlines Flight 611." Wikipedia: The Free Encyclopedia Wikimedia Foundation, Inc. 14 Oct. 2011. Web. 25 Oct. 2011. 11. "Chalk's Ocean Airways Flight 101." Wikipedia: The Free Encyclopedia Wikimedia Foundation, Inc Web. 12. ACC Public Affairs. F 15 Eagle Accident Report Released 2008, available at http://www.af.mil/news/story.asp?storyID=123081718 13. "Southwest Airlines Flight 2294." Wikipedia: The Free Encyclopedia Wikimedia Foundation, Inc. 18 Apr. 2011. Web. 25 Oct. 2011. 14. Beden S, Abdullah S, Ariffin A. Review of fatigue crack propagation models in metallic components. European Journal of Scientific Research 2009; 28 (3):364 397. 15. Ande rson TL. Fracture Mechanics: Fundamentals and Applications CRC Press: Boca Raton, FL, 2004. PAGE 178 178 16. Nikishkov GP, Atluri SN. Calculation of fracture mechanics parameters for an arbitrary three dimensional crack by the 'equivalent domain integral method'. Inte rnational Journal for Numerical Methods in Engineering 1987; 24 (9):1801 1821. 17. Rice JR. A path independent integral and the approximate analysis of strain concentration by notches and cracks. Journal of Applied Mechanics 1968; 35 (2):379 386. 18. Shih C, Asaro R. Elastic plastic analysis of crack on bimaterial interface: Part I small scale yielding. Journal of Applied Mechanics 1988; 55 (2):299 316. 19. Yau J, Wang S, Corten H. A mixed mode crack analysis of isotropic solids using conservation laws of el asticity. Journal of Applied Mechanics 1980; 47 (2):335 341. 20. Murakami Y, editor. Stress Intensity Factors Handbook Pergamon Press: New York, 1987. 21. Sih GC. Strain energy density factor applied to mixed mode crack problems. International Journal of F racture 1974; 10 (3):305 321. 22. Pais M, Kim NH, Davis TA. Reanalysis of the extended finite element method for crack initiation and propagation. 51st AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference Orlando, FL, 2010. 23. Chapra S, Canale R. Numerical Methods for Engineers McGraw Hill: New York, 2002. 24. Lu Z, Liu Y. Small time scale fatigue crack growth analysis. International Journal of Fatigue 2010; 32 (8):1306 1321. 25. without remeshing. International Journal for Numerical Methods in Engineering 1999; 46 (1):131 150. 26. Edke MS, Chang KH. Shape design sensitivity analysis (DSA) for structural f racture using extended FEM (XFEM). International Journal of Pure and Applied Mathematics 2008; 49 (3):365 372. 27. Xiaoping H, Moan T, Weicheng C. An engineering model of fatigue crack growth under variable amplitude loading. International Journal of Fatigu e 2008; 30 (1):2 10. 28. Dabayeh AA, Topper TH. Changes in crack opening stress after underload and overloads in aluminum alloy. International Journal of Fatigue 1995; 17 (4):261 269. PAGE 179 179 29. Ray A, Patankar R. Fatigue crack growth under variable amplitude loadi ng: Part I model formulation in state space setting. Applied Mathematical Modelling 2001; 25 (11):979 994. 30. Ray A, Patankar R. Fatigue crack growth under variable amplitude loading: Part II code development and model validation. Applied Mathematical Modelling 2001; 25 (11):995 1013. 31. Huang XP, Moan T. Improved modeling of the effect of R ratio on crack growth rate. International Journal of Fatigue 2007; 29 (4):591 602. 32. Ibrahim FK, Thompson JC, Topper TH. A study of effect of mechanical variables on fatigue crack closure and propagation. International Journal of Fatigue 1986; 8 (3):135 142. 33. Kujawski D. A fatigue crack driving force parameter with load ratio effects. International Journal of Fatigue 2001; 23 (S1):S239 S246. 34. Paris P, Gomez M, A nderson W. A rational analytic theory of fatigue. The Trend in Engineering 1961; 13 :9 14. 35. Shatil G, Ellison EG, Smith DJ. Elastic plastic behavior and uniaxial low cycle fatigue life of notched specimens. Fatigue and Fracture of Engineering Materials a nd Structures 1995; 18 (2):235 245. 36. Zhang JZ, Halliday MD, Bowen P, Poole P. Three dimensional elastic plastic finite element modelling of small fatigue crack growth under a single tensile overload. Engineering Fracture Mechanics 1999; 63 (3):229 251. 37 Elguedj T, Gravouil A, Combescure A. Appropriate extended function for X FEM simulation of plastic fracture mechanics. Computer Methods in Applied Mechanics and Engineering 2006; 195 (7 8):501 515. 38. Elguedj T, Gravouil A, Combescure A. A mixed augmente d Lagrangian extended finite element method for modelling elastic plastic fatigue crack growth with unilateral contact. International Journal for Numerical Methods in Engineering 2007; 71 (13):1569 1597. 39. Maligno AR, Rajaratnam S, Leen SB, Williams EJ. A three dimensional (3D) numerical study of fatigue crack growth using remeshing techniques. Engineering Fracture Mechanics 2010; 77 (1):94 111. 40. Elices M, Guinea GB, Gomez J, Planas J. The cohesive model: Advantages, limitations and challenges. Engineeri ng Fracture Mechanics 2002; 69 (2):137 163. 41. Moody NR, Reedy ED, Kent MS. Physical basis for interfacial traction separation laws Sandia Report SAND2002 8567, Sandia National Laboratories, 2002. PAGE 180 180 42. Scheider I. Derivation of separation laws for cohesive models in the course of ductile fracture. Engineering Fracture Mechanics 2009; 76 (10):1450 1459. 43. Babuska I, Melenk JM. The partition of unity method. International Journal for Numerical Methods in Engineering 1998; 40 (4):727 758. 44. Osher S, Sethian J. Fronts propagating with curvature dependent speed: Algorithms based on Hamilton Jacobian formulations. Journal of Computational Physics 1988; 79 (1):12 49. 45. Mulder W, Osher S, Sethian J. Computing interface motion in compressible gas dynamics. Journal of Computational Physics 1992; 100 (2):209 228. 46. Sussman M, Smereka P, Osher S. A level set approach for computing solutions to incompressible two phase flow. Journal of Computational Physics 1994; 114 (1):146 159. 47. Osher S, Fedkiw R. Level Set Method s and Dynamic Implicit Surfaces Springer Verlag: New York, 2002. 48. Osher S, Paragios N. Geometric Level Set Methods in Imaging, Vision, and Graphics Springer Verlag: New York, 1993. 49. Adalsteinsson D, Sethian J. A level set approach to a unified mode l for etching, deposition, and lithography I: Algorithms and two dimensional simulations. Journal of Computational Physics 1995; 120 (1):128 144. 50. Adalsteinsson D, Sethian J. A level set approach to a unified model for etching, deposition, and lithograph y II: Three dimensional simulations. Journal of Computational Physics 1995; 122 (1):348 366. 51. Adalsteinsson D, Sethian J. A level set approach to a unified model for etching, deposition, and lithography III: Redeposition, reemission, surface diffusion, a nd complex simulations. Journal of Computational Physics 1997; 138 (1):193 223. 52. Wang S, Wang M. Radial basis functions and level set method for structural topology optimization. International Journal for Numerical Methods in Engineering 2006; 65 (12):206 0 2090. 53. Wang S, Lim K, Khoo B, Wang M. An extended level set method for shape and topology optimization. Journal of Computational Physics 2007; 221 (1):395 421. 54. dimensional non planar crack growth by a co upled extended finite element and fast marching method. International Journal for Numerical Methods in Engineering 2008; 76 (5):727 748. PAGE 181 181 55. F. A computational approach to handle complex microstructure geometries. Co mputer Methods in Applied Mechanics and Engineering 2003; 192 (28 30):3163 3177. 56. Courant R, Friedrichs K, Lewy H. Uber die partiellen differenzengleichungen der mathematischen physik. Mathematische Annalen 1928; 100 (1):32 74. 57. Stolarska M, Chopp DL, sets in the extended finite element method. International Journal for Numerical Methods in Engineering 2001; 51 (8):943 960. 58. planar 3D crack growth by the extend ed finite element and level sets Part II: Level set update. International Journal for Numerical Methods in Engineering 2002; 53 (11):2569 2586. 59. planar 3D crack growth by the extended finite element and level sets Part I: Mechanical model. International Journal for Numerical Methods in Engineering 2002; 53 (11):2549 2568. 60. Sukumar N, Chopp DL, Moran B. Extended finite element method and fast marching method for three dimensional fatigue crack propagation. Engine ering Fracture Mechanics 2003; 70 (1):29 48. 61. International Journal for Numerical Methods in Engineering 2001; 50 (4):993 1013. 62. Reddy JN. An Introduction to the Finite Element Method McGraw Hill: New York, 20 05. 63. Belytschko T, Gracie R, Ventura G. A review of extended/generalized finite element methods for material modeling. Modeling and Simulation in Materials Science and Engineering 2009; 17 (4):1 24. 64. Karihaloo BL, Xiao QZ. Modelling of stationary and growing cracks in FE framework without remeshing: A state of the art review. Computers and Structures 2003; 81 (3):119 129. 65. Yazid A, Abdelkader N, Adbelmadjid H. A state of the art review of the X FEM for computational fracture mechanics. Applied Mathem atical Modelling 2009; 33 (12):4269 4282. 66. Fries T P, Belytschko T. The extended/generalized finite element method: An overview of the method and its applications. International Journal for Numerical Methods in Engineering 2010; 84 (3):253 304. 67. Aquino W, Brigham JC, Earles CJ, Sukumar N. Generalized finite element method using proper orthogonal decomposition. International Journal for Numerical Methods in Engineering 2009; 79 (7):887 906. PAGE 182 182 68. Fries T P, Byfut A, Alizada A, Cheng KW, Schrder A. Hanging nodes and XFEM. International Journal for Numerical Methods in Engineering 2010: n/a. 69. Bordas S, Nguyen PV, Dunant C, Nguyen Dang H, Guidoum A. An extended finite element library. International Journal for Numerical Methods in Engineering 2006; 71 (6):70 3 732. 70. Sukumar N, Dolbow JE, Devan A, Yvonnet J, C hinesta F, Ryckelynck D, Lorong unity finite elements. International Journal of Forming Processes 2005; 8 (4):409 427. 71. Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. International Journal for Numerical Methods in Engineering 1999; 45 (5):601 620. 72. Fleming M, Chu YA, Moran B, Belytschko T. Enriched element free Galerkin methods for crack tip fields. International Journal for Numerical Methods in Engineering 1997; 40 (8):1483 1504. 73. static crack growth with the extended finite element method, Part I: Computer implementation. International Journal of Solids and Structures 2003; 40 (26):7513 7537. 74. i static crack growth with the extended finite element method, Part II: Numerical applications. International Journal of Solids and Structures 2003; 40 (26):7539 7552. 75. three di mensional crack modelling. International Journal for Numerical Methods in Engineering 2000; 48 (11):1549 1570. 76. Nagashima T, Miura N. Elastic plastic fracture analysis for surface cracks using X FEM. 8th International Conference on Computational Plastici ty Barcelona, Spain, 2005. 77. Zi G, Chen H, Xu J, Belytschko T. The extended finite element method for dynamic fractures. Shock and Vibration 2005; 12 (1):9 23. 78. Prabel B, Marie S, Combescure A. Using the X FEM method to model the dynamic propagation a nd arrest of cleavage cracks in ferritic steel. Engineering Fracture Mechanics 2008; 75 (10):2985 3009. 79. FEM explicit dynamics for constant strain elements to alleviate mesh constraints on internal or external bou ndaries. Computer Methods in Applied Mechanics and Engineering 2008; 197 (30 32):349 363. PAGE 183 183 80. conserving scheme for dynamic crack growth using the extended finite element method. International Journal for Numer ical Methods in Engineering 2005; 63 (5):631 659. 81. Zamani A, Eslami MR. Implementation of the extended finite element method for dynamic thermoeleastic fracture initiation. International Journal of Solids and Structures 2010; 47 (10):1392 1404. 82. Hutchi nson JW. Singular behavior at the end of a tensile crack in a hardening material. Journal of Mechanics and Physics of Solids 1968; 16 (1):13 31. 83. Rice JR, Rosengren GF. Plane strain deformation near a crack tip in a power law hardening material. Journal of Mechanics and Physics of Solids 1968; 16 (1 12). 84. Anderson MR. Fatigue Crack Initiation and Growth in Ship Structures. Ph.D. Thesis Department of Naval Architecture and Offshore Engineering Technical University of Denmark, 1998. 85. Sukumar N, Huang bimaterial interface cracks. International Journal for Numerical Methods in Engineering 2004; 59 (8):1075 1102. 86. Eng ineering Fracture Mechanics 2002; 69 (7):813 833. 87. Unger JF, Eckardt S, Knke C. Modelling of cohesive crack growth in concrete structures with the extended finite element method. Computer Methods in Applied Mechanics and Engineering 2007; 196 (41 44):408 7 4100. 88. Zi G, Belytschko T. New crack tip elements for XFEM and applications to cohesive cracks. International Journal for Numerical Methods in Engineering 2003; 57 (15):2221 2240. 89. intersecting cracks with the extended finite element method. International Journal for Numerical Methods in Engineering 2000; 48 (12):1741 1760. 90. T. An extended finite element method for modeling crack growth with frictional contact. Computer Methods in Applied Mechanics and Engineering 2001; 190 (51 52):6825 6846. 91. Baietto MC, Pierres E, Gravouil A. A multi model X FEM strategy dedicated to fric tional crack growth under cyclic fretting fatigue loadings. International Journal of Solids and Structures 2010; 47 (10):1405 1423. 92. cracking of thin films with the extended finite element method. Engineering Fracture Mechanics 2003; 70 (18):2513 2526. PAGE 184 184 93. inte rfacial crack in a four point bend test. Engineering Fracture Mechanics 2005; 72 (17):2584 2601. 94. Asadpoure A, Mohammadi S. Developing new enrichment functions for crack simulation in orthotropic media by the extended finite element method. International Journal for Numerical Methods in Engineering 2006; 69 (10):2150 2172. 95. FEM to the fracture of piezoelectric materials. International Journal for Numerical Methods in Engineering 2008; 77 (11):1535 1565. 96. Mousavi SE, Grinspun E, Sukumar N. Harmonic enrichment functions: A unified treatment of multiple, intersecting and branched cracks in the extended finite element method. International Journal for Numerical Methods in Engineering 2011; 85 (10):1306 1322. 97. Re issner plates with the extended finite element method. International Journal of Solids and Structures 2000; 37 (48 50):7161 7183. 98. Bachene M, Tiberkak R, Rechak S. Vibration analysis of cracked plates using the extended finite element method. Archive of Applied Mechanics 2009; 79 (3):249 262. 99. Rabinovich D, Givoli D, Vigdergauz S. XFEM based crack detection scheme using a genetic algorithm. International Journal for Numerical Methods in Engineering 2007; 71 (9):1051 1080. 100. Rabinovich D, Givoli D, Vig dergauz S. Crack identification by 'arrival time' using XFEM and a genetic algorithm. International Journal for Numerical Methods in Engineering 2009; 77 (3):337 359. 101. Waisman H, Chatzi E, Smyth AW. Detection and quantification of flaws in structures by the extended finite element method and genetic algorithms. 10th US National Congress on Computational Mechanics Columbus, OH, 2009. 102. Fries T P. A corrected XFEM approximation without problems in blending elements. International Journal for Numerical Methods in Engineering 2007; 75 (5):503 532. 103. Laborde P, Pommier J, Renard YS, M. High order extended finite element method for cracked domains. International Journal for Numerical Methods in Engineering 2005; 64 (3):354 381. PAGE 185 185 104. Stazi FL, Budyn E, Ches sa J, Belytschko T. An extended finite element method with higher order elements for curved cracks. Computational Mechanics 2003; 31 (2 3):38 48. 105. Ventura G, Gracie R, Belytschko T. Fast integration and weight function blending in the extended finite el ement method. International Journal for Numerical Methods in Engineering 2008; 77 (1):1 29. 106. Tarancon J, Vercher A, Giner E, Fuenmayor J. Enhanced blending elements for XFEM applied to linear elastic fracture mechanics. International Journal for Numeric al Methods in Engineering 2008; 77 (1):126 148. 107. Chanine E, Laborde P, Renard Y. Crack tip enrichment in the XFEM using a cutoff function. International Journal for Numerical Methods in Engineering 2008; 75 (6):629 646. 108. Chessa J, Wang H, Belytschko T. On the construction of blending elements for local partition of unity enriched finite elements. International Journal for Numerical Methods in Engineering 2003; 57 (7):1015 1038. 109. Ventura G, Moran B, Belytschko T. Dislocation by partition of unity. I nternational Journal for Numerical Methods in Engineering 2005; 62 (11):1463 1487. 110. robustness study of the X FEM for stress analysis around cracks. International Journal for Numerical Methods in Engineering 2005; 64 (8):1033 1056. 111. Krongauz Y, Belytschko T. EFG approximation with discontinous derivatives. International Journal for Numerical Methods in Engineering 1998; 41 (7):1215 1233. 112. Pais M, Kim NH, Peters J. Discussion on modeling weak discontinuities independent of the finite element mesh. 10th US Nation al Congress on Computational Mechanics Columbus, OH, 2009. 113. level sets in the extended finite element method. Computer Methods in Applied Mechanics and Engineering 2001; 190 (46 47):6183 6200. 114. Zilian A, Fries T P. A localized mixed hybrid method for imposing interfacial constraints in the extended finite element method (XFEM). International Journal for Numerical Methods in Engineering 2009; 79 (6):733 752. 115. Pais M, Kim NH. Modeling failure in composite materials with the extended finite element and level set method. 50th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference Palm Springs, CA, 2009. PAGE 186 186 116. Dolbow JE, Merle R. Solving thermal and phase change problems with the extended finite element method. Computational Mechanics 2002; 28 (5):339 350. 117. Ji H, Dolbow JE. On strategies for enforcing interfacial constraints and evaluating jump conditions with the extended finite element method. International Journal for Numerical Methods in Engineering 2004; 61 (14):2508 2535. 118. Chessa J, Smolinski P, Belytschko T. The extended finite element method (XFEM) for Stefan problems. International Journal for Numerical Methods in Engineering 2002; 53 (8):1959 1977. 119. Hettich T, Hund A, Ramm E. Modeling of failur e in composites by X FEM and level sets within a multiscale framework. Computer Methods in Applied Mechanics and Engineering 2008; 197 (5):414 424. 120. Hettich T, Ramm E. Interface material failure modeled by the extended finite element method and level se ts. Computer Methods in Applied Mechanics and Engineering 2006; 195 (37 40):4753 4767. 121. Mousavi SE, Sukumar N. Generalized gaussian quadrature rules for disconinuities and crack singularities in the extended finite element method. Computer Methods in Ap plied Mechanics and Engineering 2010; 199 (49 52):3237 3249. 122. Mousavi SE, Sukumar N. Generalized Duffy transformation for integrating vertex singularities. Computational Mechanics 2010; 45 (2 3):127 140. 123. Natarajan S, Mahapatra DR, Bordas S. Integrat ing strong and weak discontinuities without integration subcells and example applications in an XFEM/GFEM framework. International Journal for Numerical Methods in Engineering 2010; 83 (3):269 294. 124. Yamada T, Nagashima T. A simple and efficient hybrid n umerical quadrature scheme for structured extended FEM. 10th US National Congress on Computational Mechanics Columbus, OH, 2009. 125. Park K, Pereira JP, Duarte C, Paulino GH. Integration of singular enrichment functions in the generalized/extended finite element method for three dimensional problems. International Journal for Numerical Methods in Engineering 2008; 78 (10):1220 1257. 126. Global Engineering and Materials, Inc. Computational Software Development 2008, available at http://www.gem consultant.com/Demo.aspx 127. Dassault Systemes. Abaqus 6.9 Extended Functionality (6.9 EF) Overview November 2009 2009, available at http:// www.simulia.com/download/Webinars/6.9_EF/ PAGE 187 187 128. Dassault Systemes. Abaqus 6.9 Overview June 2009 2009, available at http://www.simulia.com/download/Web inars/Abaqus69_overview/Abaqus69Overvi ew new.html 129. Giner E, Sukumar N, Tarancon J, Fuenmayor J. An Abaqus implementation of the extended finite element method. Engineering Fracture Mechanics 2009; 76 (3):347 368. 130. Hansbo A, Hansbo P. A finite elem ent method for the simulation of strong and weak discontinuities in solid mechanics. Computer Methods in Applied Mechanics and Engineering 2004; 193 (33 35):3524 3540. 131. Song J H, Areias PMA, Belytschko T. A method for dynamic crack and shear band propag ation with phantom nodes. International Journal for Numerical Methods in Engineering 2006; 67 (6):868 893. 132. Rabczuk T, Zi G, Gerstenberger A, Wall WA. A new crack tip element for the phantom node method with arbitrary cohesive cracks. International Jour nal for Numerical Methods in Engineering 2008; 75 (5):577 599. 133. Liu Y, Mahadevan S. Threshold stress intensity factor and crack growth rate prediction under mixed mode loading. Engineering Fracture Mechanics 2007; 74 (3):332 345. 134. Pais, M. 2D MATLAB XFEM Codes 2010, available at http://www.matthewpais.com/2Dcodes 135. Rhee H, Salama M. Mixed mode stress intensity factors solutions for a warped surface flaw by three dimensional finite element analys is. Engineering Fracture Mechanics 1987; 28 (2):203 209. 136. Tanaka K. Fatigue crack propagation from a crack inclined to the cyclic tension axis. Engineering Fracture Mechanics 1974; 6 (3):493 507. 137. Yan X, Zhang Z. Mixed mode criteria for the materials with different yield strengths in tension and compression. Engineering Fracture Mechanics 1992; 42 (1):109 116. 138. Viana, FAC. SURROGATES Toolbox User's Guide, Version 2.1 2010, available at http://sites.google.com/site/fchegury/surrogatestoolbox 139. Cenaero. Morfeo 2010, available at http://www.cenaero.be/Page_Generale.asp?DocID=21686&la=1&langue= EN 140. Shi J, Chopp DL, Lua J, Sukumar N, Belytschko T. Abaqus implementation of extended finite element method using a level set representation for three dimensional fatigue crack growth and life predictions. Engineering Fracture Mechanics 2010; 77 (14) :2840 2863. PAGE 188 188 141. Wyart E, Duflot M, Coulon D, Martiny P, Pardoen T, Remacle J F, Lani F. Substructuring FE XFE approaches applied to the three dimensional crack propagation. Journal of Computational and Applied Mechanics 2008; 215 (2):626 638. 142. Renard, Y, Pommier, J. Getfem++ 2010, available at http://download.gna.org/getfem/html/homepage 143. Bordas, S. Downloads 2010, available at http://people.civil.gla.ac.uk/~bordas/codes.html 144. Cherepanov GP. The propagation of cracks in a continuous medium. Journal of Applied Mathematics and Mechanics 1967; 31 (3):503 512. 145. Westergaard I. Bearing pressures and cracks. Journal of A pplied Mechanics 1939; 6 (1):49 53. 146. Williams M. On the stress distribution at the base of a stationary crack. Journal of Applied Mechanics 1957; 24 (1):109 114. 147. Duarte C, Hamzeh OH, Liszka TJ, Tworzydlo WW. A generalized finite element method for t he simulation of three dimensional dynamic crack propagation. Computer Methods in Applied Mechanics and Engineering 2001; 190 (15 17):2227 2262. 148. Lua J, Shi J, Lu Z, Liu Y. Life prediction of aerospace structures by combined XFEM and advanced fatigue mo dels. 51st AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference Orlando, FL, 2010. 149. Liu XY, Xiao QZ, Karihaloo BL. XFEM for direct evaluation of mixed mode SIFs in homogenous and bi materials. International Journal for Nume rical Methods in Engineering 2004; 59 (8):1103 1118. 150. Bilby BA, Cardew GE. The crack with a kinked tip. International Journal of Fracture 1975; 11 (4):708 712. 151. Erdogan F, Sih GC. On the crack extension in plates under plane loading and transverse sh ear. Journal of Basic Engineering 1963; 85 :519 525. 152. Nuismer RJ. An energy release rate criterion for mixed mode fracture. International Journal of Fracture 1975; 11 (2):708 712. 153. Richard HA, Fulland M, Sander M. Theoretical crack path prediction. F atigue and Fracture of Engineering Materials and Structures 2005; 28 (1 2):3 12. 154. Feng M, Ding F, Jiang Y. A study of loading path influence on fatigue crack growth under combined loading. International Journal of Fatigue 2006; 28 (1):19 27. PAGE 189 189 155. Guidaul t P A, Allix O, Champaney L, Cornuault C. A multiscale extended finite element method for crack propagation. Computer Methods in Applied Mechanics and Engineering 2008; 197 (5):381 399. 156. Sun CT. Mechanics of Aircraft Structures John Wiley and Sons: Hob oken, 2006. 157. Tanaka K, Nakai Y, Yamashita M. Fatigue growth threshold of small cracks. International Journal of Fatigue 1981; 17 (5):519 533. 158. Dassault Systemes. Dassault Systemes Accelerates Eploration of Real World Product Performance with Abaqus 6.10 EF 2010, available at http://www.3ds.com/company/news medi a/press releases detail/release/dassault systemes accelerates exploration of real world/single/3233/?cHash=88be0d34424cf37fc71fb63619e5849d 159. Langlais TE, Vogel JH, Chase TR. Multiaxial cycle counting for critical plane methods. International Journal of Fatigue 2003; 25 (7):641 647. 160. Laborde P, Pommier J, Renard Y, Salan M. High order extended finite element method for cracked domains. International Journal for Numerical Methods in Engineering 2005; 64 (3):354 381. 161. Khosrovanah AK, Dowling NE. F atigue loading history reconstruction based on the rainflow technique. International Journal of Fatigue 1990; 12 (2):99 106. 162. Rice RC, editor. Fatigue Design Handbook Society of Automotive Engineers: Warrendale, 1997. 163. Rychlik I. A new definition o f the rainflow cycle counting method. International Journal of Fatigue 1987; 9 (2):119 121. 164. Rychlik I. Simulation of load sequences from rainflow matrices: Markov method. International Journal of Fatigue 1996; 18 (7):429 438. 165. Collins JA. Failure of Materials in Mechanical Design: Analysis, Prediction, Prevention John Wiley and Sons: New York, 1993. 166. Nieslony A. Determination of fragments of multiaxial service loading strongly influencing the fatigue of machine components. Mechanical Systems and Signal Processing 2009; 23 (8):2712 2721. 167. Nieslony, A. Rainflow for MATLAB 2010, available at http://www.mathworks.com/matlabcentral/fx_files/3026/3/content/ind ex.html 168. Karolczuk A, Macha E. A review of critical plane orientation in multiaxial fatigue failure criteria of metallic materials. International Journal of Fracture 2005; 134 (3 4):267 304. PAGE 190 190 169. Papadopoulous IV, Davoli P, Foria C, Fillippini M, Bern asconi A. A comparative study of multiaxial high cycle fatigue criteria for metals. International Journal of Fatigue 1997; 19 (3):219 235. 170. Pitoiset X, Preumont A. Spectral methods for multiaxial random fatigue analysis of metallic structures. Internati onal Journal of Fatigue 2000; 22 (7):541 550. 171. Wang CH, Brown MW. Life prediction techniques for variable amplitude multiaxial fatigue Part I: theories. Journal of Engineering Materials and Technology 1996; 118 (3):367 370. 172. Wang CH, Brown MW. Life prediction techniques for variable amplitude multiaxial fatigue Part II: comparison with experimental results. Journal of Engineering Materials and Technology 1996; 118 (3):371 374. 173. Newman JC. FASTRAN II A fatigue crack growth structural analysis program NASA 104159, NASA, 1992. 174. Myers R. Classical and Modern Regression with Applications Duxbury Thomson Learning: Pacific Grove, 2000. 175. Kleijnen J. Kriging metamodeling in simulation: A review. European Journal of Operational Research 2009; 182 (3):707 716. 176. Martin J, Simpson T. Use of kriging models to approximate deterministic computer models. AIAA Journal 2005; 43 (4):853 863. 177. Stein M. Interpolation of Spatial Data: Some Theory for Kriging Springer: New York, 1999. 178. Cheng B, Ti tterington D. Neural networks: A review from a statistical perspective. Statistical Science 1994; 9 (1):2 54. 179. Smola A, Scholkopf B. A tutorial on support vector regression. Statistics and Computing 2004; 14 (3):199 222. 180. Tada H, Paris P, Irwin G. Th e Stress Analysis of Cracks Handbook ASME Press: New York, 1987. 181. Murakami Y. Stress Intensity Factors Handbook Pergamon Press: New York, 1987. 182. Waisman H, Chatzi E, Smyth AW. Detection and quantification of flaws in structures by the extended fi nite element method and genetic algorithms. International Journal for Numerical Methods in Engineering 2009; 82 (3):303 328. 183. Abu Kassim AM, Topping BHV. Static renalysis: A review. Journal of Structural Engineering 1987; 113 (5):1029 1045. PAGE 191 191 184. Wu B, Li Z. Static reanalysis of structures with added degrees of freedom. Communications in Numerical Methods in Engineering 2006; 22 (4):269 281. 185. Wu B, Lim C, Li Z. A finite element algorithm for reanalysis of structures with added degrees of freedom. Finite Elements in Analysis and Design 2004; 40 (13 14):1791 1801. 186. Sherman J, Morisson W. Adjustment of an inverse matrix corresponding to a change in one element of a given matrix. Annals of Mathematical Statistics 1950; 21 :124 127. 187. Li Z, Lim C, Wu B. A comparison of several reanalysis methods for structural layout modifications with added degrees of freedom. Structural and Multidisciplinary Optimization 2008; 36 (4):403 410. 188. Chen Y, Davis TA, Hager WW, Rajamanickam S. Algorithm 887: CHOLMOD, supern odal sparse Cholesky factorization and update/downdate. ACM Transactions on Mathematical Software 2009; 35 (3):1 14. 189. Davis TA. Direct Methods for Sparse Linear Systems Society of Industrial and Applied Mathematics: Philadelphia, 2006. 190. Davis, TA. CHOLMOD: Supernodal Sparse Cholesky Factorization and Update/Downdate 2010, available at http://www.cise.ufl.edu/research/sparse/cholmod/ 191. Davis TA, Hager WW. Modifying a sparse Choles ky factorization. SIAM Journal on Matrix Analysis and Applications 1999; 20 (3):606 627. 192. Davis TA, Hager WW. Multiple rank modifications of a sparse Cholesky factorization. SIAM Journal on Matrix Analysis and Applications 2001; 22 (4):997 1013. 193. Dav is TA, Hager WW. Row modifications of a sparse Cholesky factorization. SIAM Journal on Matrix Analysis and Applications 2005; 26 (3):621 639. 194. Davis TA, Hager WW. Dynamic supernodes in sparse Cholesky update/downdate and triangular solves. ACM Transacti ons on Mathematical Software 2009; 35 (4):1 23. 195. Amestoy M, Davis TA, Duff IS. An approximate minimum degree ordering algorithm. SIAM Journal on Matrix Analysis and Applications 1996; 17 (4):886 905. 196. Amestoy P, Davis TA, Duff IS. Algorithm 8 37: AMD, an approximate minimum degree ordering algorithm. ACM Transactions on Mathematical Software 2004; 30 (3):381 388. PAGE 192 192 197. "g Force." Wikipedia: The Free Encyclopedia Wikimedia Foundation, Inc. 21 Oct. 2011. Web. 25 Oct. 2011. 198. Aerospace Web. Bank Angle and G's 2010, available at http://www.aerospaceweb.org/question/performance/q0146.shtml 199. "Banked Turn." Wikipedia: The Free Encyclopedia Wikimedia Foundation, Inc. 15 Oct. 2011. Web. 25 Oct. 2011. 200. Boeing. 767 Family Technical Information 2010, available at http://www.boeing.com/commercial/767family/specs.html 201. "Jet Fuel." Wikipedia: T he Free Encyclopedia Wikimedia Foundation, Inc. 23 Oct. 2011. Web. 25 Oct. 2011. 202. "Angle of Attack." Wikipedia: The Free Encyclopedia Wikimedia Foundation, Inc. Web. 203. Pippy R. Approximate Probabilistic Optimization of a Wingbox Model Using Exact Capacity Approximate Response Distribution (ECARD). M.S. Thesis Mechanical Engineering University of Florida, 2008. 204. Choklovshi, T. Pointed tip wings at low Reynolds numbers 2011, available at http://www scf.usc.edu/~tchklovs/Proposal.htm 205. Lobert G. Spanwise lift distribution of forward and aft swept wings in comparison to the optimum distribution form. Jornal of Aircraft 1981; 18 (6):496 498. 206. NASA. Inclination effects on lift 2010, available at http://www.grc.nasa.gov/WWW/K 12/airplane/incline.html 207. Rocket and Space Technology. Drag Coefficient 2011, available at http://www.braeunig.us/space/cd.htm 208. "Lift (force)." Wikipedia: The Free Encyclopedia Wikimedia Foundation, Inc. 24 Oct. 2011. Web. 25 Oct. 2011. 209. He MY, Hutchinson JW, Becher PF. Crack deflection at an interface between dissimilar elastic materials. International Journal for Solids and Structures 1989; 25 (9):1053 1067. 210. He MY, Hsueh CH, Becher PF. Deflection versus penetration of a wedge loaded crack: Effects of branch crack length and penetrated layer width. Composites: Part B 2000; 31 (4):299 308. 211. Zhang Z, Suo Z. Split singularities and the competition between crack penetration and debond at a bimaterial interface. International Journal of Solids and Structures 2007; 44 (13):4559 4573. PAGE 193 193 212. He JW, Li N, Jia M, Wen SP, Jia R P. Effect of interface layer on fatigue crack growth in an aluminum laminate. Fatigue and Fracture of Engineering Materials and Structures 2004; 27 (3):171 176. 213. Kao Walter S, Sthle P, Chen SH. A finite element analysis of a crack penetrating or deflec ting into an interface in a thin laminate. Key Engineering Materials 2006; 312 :173 178. 214. Archard JF. Contact and rubbing of flat surfaces. Journal of Applied Physics 1953; 24 (8):981 988. 215. Mukras S, Kim NH, Sawyer WG, Jackson DB, Bergquist LW. Numer ical integration schemes and parallel computation of wear prediction using finite element method. Wear 2009; 266 (7 8):822 831. 216. Sutradhar A, Paulino GH, Gray LJ. Symmetric Galerkin Boundary Element Method Springer Verlag: Heidelberg, 2008. PAGE 194 194 BIOGRAP HICAL SKETCH Matthew Jon Pais graduated from the University of Missouri in 2007 magna c um laude and h onors s cholar with a Bachelor of Science in Mechanical Engineering. As an undergraduate student he participated in a National Science Foundation sponsored Research Experience for Undergraduates (REU) at the U niversity of Missouri at Rolla and worked as an undergraduate research assistant which resulted in the undergraduate honors thesis, Hydroxyapatite Reinforced Dental Composites He joined the University o f Florida in 2007 after receiving the Alumni Fellowship where he worked with Dr. Nam Ho Kim on his Ph.D. in Mechanical Engineering He spent the summer of 2010 at the Idaho National Laboratory where he modeled failure by liquid metal embrittlement. He rece ived his Ph.D. from the University of Florida in the fall of 2011. 