UFDC Home  myUFDC Home  Help 



Full Text  
ANALYTICAL AND NUMERICAL EVALUATION OF SUBSURFACE STRESSES IN ANISOTROPIC (SINGLECRYSTAL) CONTACTS By ERIK C. KNUDSEN A THESIS PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE UNIVERSITY OF FLORIDA 2003 Copyright 2003 by Erik C. Knudsen This thesis is dedicated to my grandmother, Madeline Rosenthal. I will forever miss her and always cherish the memories. ACKNOWLEDGMENTS I thank my advisor, Nagaraj Arakere, for his support and guidance while I did this research. Greg Swanson, Greg Duke, Gilda Batista, and JeffRayburn were a big help when I was at the Marshall Space Flight Center. My labmates (Jeff, Shadab, and Srikant) were an invaluable resource. My parents have been a source of wisdom and comfort all my life and none of this would have been possible without them. TABLE OF CONTENTS Page A C K N O W L E D G M E N T S ................................................................................................. iv LIST OF TABLES ............. ................... ........... .. .. ............... .. vii LIST O F FIG U RE S ......................................................... ......... .. ............. viii ABSTRACT .............. ..................... .......... .............. xii CHAPTER 1 IN TR OD U CTION ............................................... .. ......................... .. R e search O bje ctiv e s............... ................................ .................. .............. ............. ..... Introduction to N ickelBase Superalloys.................................. ........................ 3 C oordinate Transform nations ............................................................................. 4 The [110] O orientation ................................................. .... .. ............... The [213] O orientation ................................................... ........ ............... 2 A N A L Y TIC A L SO LU TIO N ........................................................... .....................14 3 FINITE ELEMENT FORMULATIONAPPLIED STRESS MODEL......................22 4 RESULTS OF THE APPLIEDSTRESS MODEL ..............................................33 Comparison of Contour Plots: AppliedStress Model and Analytical Solution.........34 Numerical Comparison of the Analytical and AppliedStress Models ...................39 5 C O N T A C T M O D E L ......................................................................... ...................42 Examination of the Normal Stresses at the Surface................................ ........ 44 Comparison of Contact Model Stresses to the Analytical Solution .........................48 6 DEFORM ATION M ECHANISM S.................................. .................................... 56 7 L ITE R A TU R E R E V IE W ........................................ ............................................67 8 CON CLU SION S .................................. .. .......... .. .............73 9 RECOMMENDATIONS FOR FUTURE WORK ............................................... 74 v L IST O F R E F E R E N C E S ...................................... .................................... ....................75 B IO G R A PH IC A L SK E TCH ...................................................................... ..................77 LIST OF TABLES Table p 11 D direction cosines .............. ...................... ........................... .... .... 9 12 Direction cosines for the [110] orientation ............................................................9 51 Comparison of Hertzian and numerical contact pressures and patch sizes.............47 6 1 S lip S y stem s ....................................................... ................ 6 0 71 N um ber of elem ents used ........................................................................... .... ... 68 LIST OF FIGURES Figure page 11 Visual description of the contact problem ............................................. ...............2 12 D am per and blade in contact................................ ....... ................ ............... 2 13 Profile of the contact problem .......................... .................. ............... 14 D am aged turbine blade............................................................... ....................... 3 15 Specim en and m material coordinate axes .......................................... ............... 5 16 A ligned coordinate system s ............................................... ............................ 6 17 Rotation of specimen coordinate axes about the z axis.........................................7 18 Twodimensional view of rotation (z axis is into the page).......................... ........8 19 Schem atic of the [213] orientation .............................................. ........ ............... 10 110 First rotation about y, y" axis ...................... ............................... ..................... 111 S econ d rotation ................................................. ................. 12 112 Second rotation about z' axis ............................................................................12 21 Profile view of contact problem with normal and tangential tractions .................. 15 22 E lliptic distribution ........ ....... ............................................................ .. .... .. ... ... 16 23 Applied forces used to derive the equilibrium equations................. ............. ...16 31 C ylinder in contact w ith slab ........................................................... ......... ......23 32 Faceon view of contact problem with applied load shown.............................. 23 33 Elastic modulus of PWA 1480 as a function of crystal orientation at room tem perature ..................................... ................................ ........... 25 34 C contact patch dim ensions................................................ ............................. 25 35 FE model; applied stress (no contact elements).................................................26 36 Constraints ......... ...... ........ ............................28 37 Elliptic contact load ......... ..... .......................... ................ .... 28 38 A applied pressure............ .............................................................. .......... ....... 29 39 Visualization of crystallographic orientation ................................ ..................... 32 41 Results were taken from the elements at the midplane ........................................33 42 Schematic of casting coordinate system.................................... ................34 43 Contour plots of the ox stress (case 0). ................................................ .... ........... 35 44 Contour plots of the Gy stress (case 0). ...................................... ...............35 45 Contour plots of the z stress (case 0). ................................................ .... ........... 36 46 Contour plots of the Txy stress (case 0). ............................................... ............... 36 47 Contour plots of the ox stress (case 33). ........................................ .... ........... 37 48 Contour plots of the Gy stress (case 33). ........................... ....................... 37 49 Contour plots of the cz stress (case 33). ............................................... .. ........... 38 410 Contour plots of the Txy stress (case 33). ........................ ....................................... 38 411 N odes used to extract results ........................................................... ......... ......39 412 ox stress vs depth .............. ............ ........ ............... .. .. ........... 40 413 Gy stress vs depth ......... .... ........................... ............. 40 414 ,z stress vs depth .......................................... ............. .... ....... 40 415 Txy stress vs depth ........................................ ......... .............. .. 41 416 Tyz stress vs depth ........................................ ........ .............. .. 41 417 Txz stress vs depth ............................... ..... .. ... ......... .............. .. 41 51 C contact m odel ................................................................42 52 Densely meshed regions of the contact model .................................... ..................43 53 Contact and target elem ents .............................................................................. 43 54 Nodes used to obtain the contact patch and peak pressure at the surface ................45 55 Y stress for case 17 .......................................... ................... ........ 45 56 Y stress for case 21 ............................................ ................. ........ 46 57 Y stress for case 25 .......................................... ................... ........ 46 58 Y stress for case 29 .......................................... ................... ........ 47 59 Dimensions of the brick element used in the contact model.............. .......... 49 510 Sigm a X com prison case 17 ............................................................................49 511 Sigm a X com prison case 25 ............................................................................50 512 Sigm a Y com prison case 17 ............................................................................50 513 Sigm a Y com prison case 25 ............................................................................51 514 Sigm a Z com prison case 17......................................................... ............... 51 515 Sigm a Z com prison case 25 ............................... ........................ ............... 52 516 Tau XY com prison case 17 ............................................................................. 52 517 Tau XY com prison case 25 .............................................................................. 53 518 Sigm a Y Z com prison case 17........................................... .......................... 53 519 Tau Y Z com prison case 25 ......... ............... .................................. ............... 54 520 Tau X Z com prison case 17.......................................................... ............... 54 52 1 T au X Z com prison case 25 ......................................................................... ... ... 55 61 Exam ple slip system ....................................................... ................. 57 62 Four octahedral planes and slip directions. ................................... ............... 57 63 Cube slip planes and slip directions ....................................................................... 58 64 Tau 1 contour plot. ....................................... ... .... ........ ......... 61 65 Tau 2 contour plot. ......................................... ... .... ........ ......... 61 66 Tau 3 contour plot. ......................................... ... .... ........ ......... 62 67 Tau 4 contour plot. ......................................... ... .... ........ ......... 62 68 Tau 5 contour plot. ......................................... ... .... ........ ......... 63 69 Tau 6 contour plot. ......................................... ... .... ........ ......... 63 610 Tau 7 contour plot. ......................................... ... .... ........ ......... 64 611 Tau 8 contour plot. ......................................... ... .... ........ ......... 64 612 Tau 9 contour plot. ......................................... ... .... ........ ......... 65 613 Tau 10 contour plot. ....................................... ... .... ........ ..... .... 65 614 Tau 11 contour plot. ....................................... ... .... ........ ..... .... 66 615 Tau 12 contour plot. ....................................... ... .... ........ ..... .... 66 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 ANALYTICAL AND NUMERICAL EVALUATION OF SUBSURFACE STRESSES IN ANISOTROPIC (SINGLECRYSTAL) CONTACTS By Erik C. Knudsen December 2003 Chair: Nagaraj K. Arakere Major Department: Mechanical and Aerospace Engineering The Space Shuttle main engine alternate high pressurefuel turbo pump (SSME AHPFTP) turbine blades commonly fail from subsurface contact stresses imparted by a frictiondamping device designed to reduce vibratory stresses. Frettingfatigue failure is a major cause of failure in turbine components. Fretting occurs when components such as blade dampers and dovetail joints are subjected to vibration, which results in contact damage. When combined with a mean stress in the component, fretting can lead to a reduction in high cycle fatigue (HCF) life. Fretting (along with wear and corrosion) initiates cracks where the blade and damper are in contact. Although the turbo pump is only operating for ten minutes, it is cycling at a very high rate, which leads to rapid accumulation of fatigue cycles. The Marshall Space Flight Center (MSFC) has developed finite element (FE) models to investigate this problem. These models have very fine meshes and also incorporate nonlinearities such as friction. The sophistication of these models is expensive from a computational standpoint. Since the blades are made of anisotropic material, the properties of the blade change with the crystallographic orientation of the material. Manufacturing tolerances limit the precise orientation of the material, thus there are numerous orientations to examine. When one factors in the time it takes to run the models and also examine the results, many hours would have to be put in to give a thorough analysis. To speed the analysis up and also make things a bit simpler, it would be helpful to avoid such large FE models (that take hours to solve and many more hours to analyze). Ideally, one should build a contact model, to get a true sense of the stresses and contact patch size. If possible, much less complicated models should be built and compared to larger contact models. If there is reasonable agreement between the two, then the less sophisticated FE models can be used outright and the results arrived at much faster. This was done here. Results for a simple, appliedstress, FE model are compared to those for the contact model. Concurrently, analytical (or closedform solutions) are also being developed for this problem. Getting results using analytical methods is almost always faster and easier then using numerical approaches. Numerical (or FE) results are also compared to those for analytical solutions. Results for numerical vs analytical approaches have fairly good agreement with the normal stress components, but the shear stresses do differ. The most probable source of that disagreement lies in the mesh density of the FE models. CHAPTER 1 INTRODUCTION Extensive work has been done in the nonconformal classical area of contact mechanics of elastic solids dating back to Hertz's seminal work in the area. Johnson (1985) gives a complete discussion of his results. Relatively little work has been focused on contact mechanics of nonconformal anisotropic solids. Interest in this field is increasing because of applications of single crystal materials used in high temperature environments. Research Objectives The primary goal of the research presented is to examine the subsurface stresses in facecentered cubic (FCC) single crystal (SC), orthotropic, turbine blade damper and dovetail attachment contacts using the finite element (FE) method. * Develop a simple FE model that has no contact elements (appliedstress problem) and compare the results to the more complex contact FE model. * Develop a contact FE model that captures the contact patch dimensions and stresses in the blade, for fully anisotropic contacts; and compare the results to Hertzian solutions. * Develop an analytical solution for anisotropic subsurface stress fields. * Compare results of both the contact and appliedstress models to analytical solutions; and determine if the analytical approach is better than the FE approach for looking at the subsurface. Figures 11 and 12 show what the contact problem looks like. The damper shown in Figure 11 is touching the axial face of the turbine blade. Notice that the damper is upsidedown in Figure 12. Figure 13 shows a profile view of the contact problem. Essentially, the damper and blade are modeled as a halfcylinder (damper) touching a rectangular slab (blade). The damper touches the blade platform at two separate locations. Only the circled location in Figure 13 is analyzed here, although this analysis could be applied to other parts of the damper. The resulting contact load is one that has a normal and tangential component. The work presented in this thesis is for normal loading only. Ultimately, cracks form and damage the blade as shown in Figure 14. Blade Coordinate System Figure 11. Visual description of the contact problem Figure 12. Damper and blade in contact Figure 13. Profile of the contact problem ,,4a I i~iGISZ~iil~c1ft ,~'i~m~ir.s Figure 14. Damaged turbine blade Source: N. Arakere, personal correspondence, 2003 Introduction to NickelBase Superalloys The material chosen for the blade is a precipitation strengthened singlecrystal Nickelbase superalloy, Pratt and Whitney (PWA) 1480/1493. Generally speaking, when working at high temperatures and high tensile stresses where creep dictates the mode of failure, it is desirable to have a material with a close packed structure, such as FCC (which happens to be the structure of this alloy); a high melting temperature; and coarse grains (for isotropic materials). The closepacked structure makes diffusion more difficult. Since grain boundaries can be the weak link at high temperatures, it is a good idea to either eliminate them or make sure one does not have too fine a grain size (Hummel 1998). The blades in the SSME AHPFTP are made from a monograin material. The absence of grain boundaries gives improved creep resistance (which is necessary since the blades operate at high temperatures). Also, with no grain boundaries, there is no need for elements, such as Boron, that are added to increase grain boundary strength, but can also lower the melting point of the alloy (Donachie and Donachie 2002). Nickelbase superalloys can be used at a higher fraction of their melting point than just about any other alloy. Also, many elements can be dissolved into Nickel. This means that lots of additives can be used for hot corrosion and oxidation resistance, solid solution strengthening and y' formation (precipitation hardening). For example, the addition of Al provides increased oxidation resistance and also forms the compound Ni3Al (y'). Chromium and Tungsten are usually added as solid solution strengtheners and Cr also helps combat corrosion (Donachie and Donachie 2002). The precipitation hardening (y') comes from heattreating an alloy that has limited solid solubility and decreasing solid solubility with decreasing temperature. The heat treatment forms a fine dispersion of precipitates (small particles) that inhibit the motion of dislocations and make the material stronger. The highest strength is achieved when many small, round, closely spaced particles are coherently spread throughout the alloy (Hummel 1998). Since these alloys have only one grain, they are orientation dependent. The orientation of the crystal lattice plays a pivotal role in the properties of the material. Typically the blades are cast in the <001>, or low modulus, direction. This orientation gives the blades great thermalmechanical fatigue resistance. It should be pointed out that one cannot cast the blades perfectly. Casting tolerances allow a 15 deviation from the stacking line (the line between the <001> direction and the airfoil stacking line) (Arakere and Swanson 2000). Coordinate Transformations To account for the orientation dependency of these materials, one must understand how to transform from one coordinate system to another. From there quantities such as the resolved shear stress (RSS) on different slip planes can be determined. Let us first begin with a few definitions. Figure 17 shows the material coordinate and specimen coordinate systems. The primed axes (x' y' z') are defined as the specimen coordinate system. The unprimed axes (x y z) are defined as the material (crystal) coordinate system. Coordinate transformations are necessary because material properties are functions of the orientation. This is not the case for an isotropic material. Since this material has cubic symmetry, there are three independent elastic constants as well as a tensor matrix that relates stress to strain (Dame and Stouffer 1996). y <010> y x<100> z <001> z Figure 15. Specimen and material coordinate axes One way to relate one coordinate system to another is through the use of direction cosines. Once one has the cosines of the angles between the specimen and material axes, the proper coordinate systems can be constructed in whatever FE code one chooses. The model is initially built in the global (x y z) coordinate system. The material is cast with a certain crystallographic orientation which must be accounted for to obtain the proper results (stresses, strains, etc.). The orientation of the material is always given, and from there direction cosines of the angles between the material and specimen coordinate system can be determined using methods described by Lekhnitskii (1963) and also Magnan (2002). The r1101 Orientation This example will deal with the [110] orientation. Let us align the x' axis with the [110] direction. To do that imagine the specimen axes (they are primed) are initially aligned with the material, unprimed, coordinate system. y <010> y' / / z' x <100> z <001> Figure 16. Aligned coordinate systems The orientation of the x' axis is known (defined as [110]). One way to find the y' and z' axes is to arbitrarily define the y' axis (it just has to be normal to the x' axis, i.e. the dot product between the two will be zero). Now that the x' axis is defined and the y' axis is arbitrarily chosen, the z' axis is obtained by taking the cross product of the x' and y' axis. The downside of employing such a method is that there is an infinite number of y' vectors normal to the x' axis. Thus, there are no unique vectors, or directions, that define where the y' and z' axes could wind up. But, the method used here does not involve an arbitrary assignment of the y' and z' axes. Continuing on, the y' axis will also be rotated about the z axis and has the direction [110]. <110> 0 y y <010> x'<110> 0 x <100> z <001> Figure 17. Rotation of specimen coordinate axes about the z axis Notice in Figure 17 the orientation (or coordinates relative to the material coordinate system) of the y' vector was not arbitrarily picked. Both the specimen and material coordinate systems are orthogonal (i.e. the angle between the x, y, z axes are 900). The entire specimen (primed) coordinate system was rotated about the z axis. The next step is to write the coordinates of the new axes with respect to the original axes. They are as follows x'= xcos O + y sin y'= xsin0 + ycosO (11) Z'= z Putting Equation (11) into matrix form x' cosO sinO 0 x y' = sin0 cosO 0 y z' 0 0 1 z y x'  x 1 Figure 18. Twodimensional view of rotation (z axis is into the page) Trigonometry gives that 0 (x y z) axes. 45 degrees (measured counterclockwise) from the tan(O)= S= 450 (13) No further transformations are needed here to obtain the direction cosines for this case. It is convenient to enter the direction cosines into a table. (12) ~I Table 11. Direction cosines x y z x' a 13Pi y y' a2 12 72 z' 13 03 73 Where at is the cosine of the angle between the x' and x axes, P1 is the angle between the x' and y axes, etc. Thus, Table 11 becomes Table 12. Direction cosines for the [110] orientation x y z x' 0.707 0.707 0 y' 0.707 0.707 0 z' 0 0 1 The [213] Orientation The last example involved only one rotation. Frequently, one will have to perform two or even three rotations to get the desired crystal orientation. This example will look at the [213] orientation. Let us look at Figure 16 again. Once more, the x' axis will be aligned with the [213] direction as shown in Figure 19. To do this transformation correctly, two rotations will be executed. First a rotation about the yy" axis by 0 degrees is performed. Then a second rotation about the z' axis by (p degrees takes place and the transformation is complete. y <010> P W x <100> z <001> Figure 16. Reprinted y <010> x <100> z <001> Figure 19. Schematic of the [213] orientation In Figure 111 one can see the first rotation from the projection (denoted by the dashed line) of the [213] vector into the xz plane. Since the dimensions of the projection are known, the angle 0 can be readily determined. zoo' '"'~r.: ~... ~ 3 tan =  2 S= tan' (14) 2 0 = 56.3090 y,y" <010> x<100> z x": projection ofx' into xz plane z <001> Figure 110. First rotation about y, y" axis For clarity, the rotation is redrawn in Figure 110. Writing the coordinates of the transformed coordinates in terms of the original set of axes we have x" cos 0 sinO x y' = 0 1 0 y (15) z"1 sinO 0 cosOI z There is just one more step which entails a rotation about the z" axis by an amount of p degrees. Figure 111 illustrates this. y <010> x <100> 2 Figure 111. Second rotation To make this second rotation clear, it is redrawn in Figure 112. y,y x": x' projected into xz plane Figure 112. Second rotation about z' axis Once more this transformation is put in matrix form y' = sin cos 0 y" z' 0 0 1 z" From looking at Figure 111, the value of p is z <001 y z',z" (16) tan l = tan (17) =15.5050 The only thing left to do is solve for the table of direction cosines. La A, / 0 cos sin 0 cos 0 0 sin0 a2 A 72 = sine cos 0 0 1 0 (18) a3 P3 73 0 0 1 sin 0 cos6O Finally, we note that counterclockwise rotations are positive and clockwise rotations are negative. Following this convention, theta is negative and phi is positive. Substituting and solving for the table of direction cosines a2 2 2 = 0.148 0.964 0.222 (19) a,3 p3 73 0.832 0 0.555 Of course one could perform a third rotation about the x' axis, perhaps by an angle y. If that were to happen, we would have al A 1 0 0 cos sino 0 cosO 0 sino a2 = c sinin/ sin 0 cos 0 0 1 0 (110) a3 /73 73 0 sin/ cosi/l 0 0 1 sinO 0 cosO CHAPTER 2 ANALYTICAL SOLUTION A closedform solution to this contact problem was formulated by Dr. Arakere (University of Florida) which utilizes complex potential stress functions developed by Lekhnitskii (1963). This is a twodimensional method with a generalized plain strain condition. The load and resulting stress functions have no zdependency. The shape of the contact patch was assumed and the applied load equations were taken from Hertz's work for linecontact problems (such as a cylinder touching a cylinder) with isotropic materials. The subsurface stresses can be solved for any general traction distribution over the surface. In chapter two it was stated that the normal pressure distribution for Herzian (nonconformal) anisotropic bodies is the same as for isotropic ones. From Hertz the normal load distribution is N(x)= po 1 x a2 (21) For a tangential load we have T(x) = p, po 1 x/a (22) Where po is the peak pressure, a is the contact halfwidth, and [f is the coefficient of friction. In Figure 13, a profile view of the contact problem was shown. That figure is reprinted in Figure 21 with a normal and tangential traction applied. If both materials are isotropic, then a solution exists for this contact problem.. Essentially all one needs are the material properties and applied loads. But since the blade is made out of an orthotropic material, there is no closedform solution available. p Q Figure 21. Profile view of contact problem with normal and tangential tractions Once again, there is only a normal load applied to the top of the cylinder (i.e. Q = 0). Hertz has proved, for nonconforming bodies in contact, that the normal pressure distribution on the block above is elliptic, Figure 22, for isotropic materials (Johnson 1985). In his paper of anisotropic contact problems, Willis (1966), states that this same distribution can be applied when anisotropic materials are used. What follows is an overview of Lehknitiskii's work applied to a nonconformal contact problem such as a cylinder touching a flat substrate. Equations that determine the subsurface stresses are developed with the main assumption is knowing the contact patch dimensions beforehand. Figure 22. Elliptic distribution Although Hertz's equations are used to describe the normal contact pressure, the analytical solution is derived in such a fashion that any loading can be applied to the surface of the anisotropic substrate (Figure 23). Where Hertz's equations are used can be seen in Equations 227 and 228. Lekhnitskii presents a stress function approach, F(x,y) and y(x,y), to obtain the subsurface stresses. The stress functions are used to satisfy the equations of equilibrium which come from Figure 23. Surface Traction Anisotropic HalfSpace Figure 23. Applied forces used to derive the equilibrium equations Equation 23 shows the relations for equilibrium. U is the applied body force. Notice that if no body forces are present Equation 23 becomes homogeneous. a0 x 8) a x 17 + o 0 (23) y The relations listed in Equation 24 are the stress functions used to satisfy the equilibrium equations. 2 F 02F 2F x = +U 2 = +U (24) By 2 a2 X axB all a'i yz z y This approach to solving these kinds of problems was initially developed by Airy in 1882 when he was solving twodimensional stress problems. A good primer to the concept of stress function can be found in Timoshenko (1982). The stress functions must satisfy the above equations (which become homogeneous when there are no body forces). Lekhnitksii integrates three stressstrain equations that come from the generalized Hooke's law. The aforementioned integration gives rise to three displacement equations. These equations must also satisfy the leftover three equations from the generalized Hooke's law. They are substituted into the remaining three stressstrain relations that give rise to another set of three equations where the coefficients of the left and righthand sides are equated. Those equations are differentiated to eliminate certain variables from the generalized displacement equations that will be determined in the future. The result of that differentiation is listed in Equation 25 and Equation 26. "2U 02U 82U L4FL = + L3 12 22) a2 + ( 6 )c (A 11 + 12) (25) ~sC2 ax@ OU OU L3F+L, = 20+ Aa34 Ba35 +(f,1 + 9, ) (,15 + f2,) (26) L2, L3, L4 are differential operators of the second, third, and fourth orders. A and B describe the bending in the xz and yz planes. 0 is the rotation about the zaxis. In the special case of plane deformation; A = B = 0 = 0. Equation (26) can be written as follows F = F + Fo (27) y = y + Vo (28) where F' and y' are the general solutions to a system of homogeneous differential equations. They are homogenous because of the absence of body forces. L4F' + L3Y' = 0 (29) L3F' + L2y' = 0 (210) To obtain the solution to the set of differential equations above, one of the functions, such as y', is eliminated. The next step is to apply the operator L2 on Equation 29 and L3 on Equation 210 and subtract what remains. The result is an equation of the sixth order (L4L2L32)F = 0 (211) To get the solution for y', the exact same procedure is followed. Since Equation (211) is of the sixth order, it can be broken down into six linear operators of the first order D6 D5 D4 D3 D2 D1F' = 0 212) Where Dk = (a/ay)Ltka/ax and [tk are the roots to the characteristic equation 2( 3 (/()= 0 (213) where 12(11) =5512 2.451 + P55 (214) 13(1) = 1531 (114 + 356).2 +(1325+ 146) B 24 (215) 14() = 11.4 21316.3 +(2.P12+ P66)2 22P261 + P22 (216) and ai3aj3 Pij = aij  a33 (217) The term aj are the elastic constants for this orthotropic material. See Equation 313. Since Equation 212 is of the sixth order, it can be resolved into the integration of six equations of the first order D1F' = (2, D2(P2 = (P3, D3(P3 = (p4, D4(4 = (P5, D5(P5 = (6, D6(6 = 0 (218) In order to integrate these first order equations, it is assumed the general integral is a function of the argument x + [ty. For example, the integral for (p6 is denoted by fV6(x + t6y). Following this convention, ps5 satisfies the nonhomogeneous equation D5(p5 = fV6(x + t6y) (219) Integrating (219) the following expression is obtained (p5 = fsIV(x + k5sy) + fV(x + X 6y)/ 16 [15 (220) Using the same method (p4, (p3, (P2, and F' will be obtained. Using a change of notation, expressions for F' and y' are determined 6 F'= FFk(x +ky) (221) k=1 (222) 6 'k= Yk +k Y) k=l Lekhnitskii goes on to prove four theorems that are necessary for the theory governing a body bounded by a cylindrical surface. He shows that a certain set of four equations cannot have real roots and that ultimately tlk is always complex or imaginary. Now the stress functions take the following form F = 2 Re[F,(zi) + F,(z,) + F (z3)] + F (223) (224) I= 2Re F(z,)+ A2F'2 (z,) + F'3() + Yo A3 Where I3(tl) I3(112) I3(3) 1 = (2 = (3 = 1(3) 12G12) 12(1 2) 14( 3) (225) Introducing the complex variable zk (k = 1, 2), new functions are introduced 1 k (zk)= F'k(zk) 3(z3)= F'3(z3) k= 1,2 ( A3 A (z)+ A2(z)+ 3 (z)= 0 ( 2m J x z /0/1 (Z) + //2 (Z)+ 303 (Z) =  ( 2;d xz A10(z)+., 202(z)+ +3(z) =( 0 (= Equations 23035 are the resulting stress functions. cx= 2.Re (112. 1(z) + 12 .2(z) + 132 3 ,3'. )) ( Ty = 2.Re(1(z) + 2(z 2) + )3 '3', )) 226) 227) Z28) 229) Z30) (231) 21 txy = 2Re(Gi'1(z) + 12'2(z) + 3 1'3' )) (232) Txz = 2Re(Gl l 'lI(z) + 22 '2(z) + ,'3 '3' ) (233) tyz = 2Re( l4'(z) + 2'2() + (z)) (234) 1 7z  (al3x+ a23'Sy + a34yz + a35'cxz + a36.'xy) a33 (235) zi= x+ rliy (236) CHAPTER 3 FINITE ELEMENT FORMULATIONAPPLIED STRESS MODEL Our first goal with this model is to evaluate the subsurface stresses in the substrate by decoupling the contact mechanics from the stress analysis. To do this, an anisotropic substrate is modeled in 3D, subject to normal loading over a contact width chosen arbitrarily. One of the difficulties in trying to get reliable results is that FE models need very fine meshes and must account for nonlinearities such as friction and plasticity. The numerical modeling procedure would be simplified considerably if the contact problem is reduced to a stress analysis problem. This can be done if the size of the contact patch and the normal contact pressure distribution is determined from Hertzian assumptions. These models are less complex than contact models and will run faster. They are easier to troubleshoot and modify should there be a problem or design change. Hertz developed closedform solutions for isotropic contact problems. His theory is grounded on four assumptions: the bodies in contact are continuous and non conforming, the strains are small, each solid is considered to be an elastic half space, and there is no friction on the surfaces (Johnson 1985). These assumptions were built into the models to analyze the anisotropic material. The applied stress model is a rectangular slab (composed of the superalloy) with an applied load. The load was determined from the Hertzian assumptions stated above. Since the contact is modeled as a cylinder touching a rectangular slab, this is a line contact problem. Hertz determined the contact patch dimensions, maximum contact pressure, and profile of the contact load (which happens to be elliptic). It has been shown that the pressure distribution for contact problems of a similar type (nonconforming) involving anisotropic materials of is also parabolic (Willis 1966). The true contact patch size and load distribution were not known when these models were first being built. The idea is to build a simple model and compare it a more sophisticated one that uses contact elements to see how well Hertz's equations will predict the load distribution and contact patch size for the anisotropic material. Figure 31. Cylinder in contact with slab I P Figure 32. Faceon view of contact problem with applied load shown Figure 32 shows the nonconformal contact problem of a cylinder on a flat substrate and in Figure 34 shows a closeup view of the contact region. Hertz (1881) determined this contact width and also the maximum contact pressure. The equations are as follows: fi1 R K + 1 (31) Ri R2 where R is the equivalent radius of curvature. E* = + v2 (32) 2P a E (33) S= 2P (34) Where P is the load per unit length, a is the contact halfwidth, and po is the maximum contact pressure. The subscripts 1 and 2 refer to the two bodies that are in contact with each other. The elastic modulii, El and E2, can be easily defined if the elastic bodies are isotropic. The problem inherent here is the contact of isotropic indenter (damper) with anisotropic substrate (blade). For the anisotropic substrate the modulus will vary as a function of crystallographic orientation. There is no procedure to determine the effective elastic modulus, E*, for an anisotropic contact, since El and E2 vary with the orientation of the material. See Figure 23; the modulus is 126 GPa for the [001] orientation and increases to 310 GPa for the [111] orientation (Sims et al. 1987). %.\ "\ \ \ ) 32.1 0 s Io Iis a6 30 :s 40 50 Figure 33. Elastic modulus of PWA 1480 as a function of crystal orientation at room temperature Source: Modified from Donachie and Donachie, 2002 2a Figure 34. Contact patch dimensions The contact dimension, 2a, is known here and is given by Equation 33. The contact takes place over an area that is rectangular (width x length). The width is 2a and the length is 0.5 inches as shown in Figure 33. Figure 35 shows the FE model used for the applied/load stress conditions. The model was built using ANSYS FE software and is composed of solid45 elements. These are 8node elements used for 3D models. This element has three degrees of freedom at each node and also has the capability to analyze creep, swelling, large deflection and large strains (ANSYS 1999). Figure 35. FE model; applied stress (no contact elements) Since a plane stress or plane strain condition cannot be enforced with this material, a 3D model must be built. The analytical solution (discussed in chapter two) uses a 'generalized plane strain' condition that generates equations for the subsurface stresses that are functions ofx and y only. When Lekhnitskii (1963) developed the closedform solution used in this analysis, his boundary conditions were prescribed on an infinitely long cylinder. This is not a plain strain condition is the conventional sense since an orthotropic material is involved. The orthotropy induces shear coupling in the z direction. One cannot introduce infinitely large elastic bodies into a numerical analysis. The dimensions used in the FE model drawn in Figure 35 were sufficiently large so the rectangular slab could be treated as an elastic half space. That is, the stresses are uniform throughout the thickness (z dimension) of the material. In most nonconformal contact problems, the contact patch is small compared to the dimensions of the bodies in contact. Consequently, stress gradients in the region affected by the contact are very large, necessitating an extremely refined FE mesh to accurately capture the subsurface stresses. The densely meshed region in Figure 35 is where the contact would take place. It is in this region that the results were extracted and where the stresses are the highest. A coarser mesh is used for the rest of the substrate since the stresses are not as high and to reduce the overall size of the model. The less elements there are the better, since the model will run faster on the computer. The base of the model, or slab, is constrained so the block cannot move in the normal, or ydirection. It is important to allow the elements to have some flexibility, thus one does not want to overconstrain the model. That is why only two of the bottom corners of the slab were constrained in the xdirection. This ensures that the model cannot rotate about the y axis, or have any rigid body motion. Now that the model is properly constrained, only entering the material properties and applying the loads remain. The load was applied using a macro (algorithm written in a computer language such as FORTRAN). The Hertzian contact pressure distribution for normal, frictionless, loading is shown above in Figure 37. The loading curve takes the form of an ellipse when the normal load P is applied to the cylinder. This distribution was applied over the densely meshed region in Figure 38. Figure 36. Constraints Figure 37 Elliptic contact load 111111 *^J~j kl ''il 17P I Figure 38. Applied pressure (x)= p 1a (36) With the geometry and loading in hand, the model was drawn and then meshed with the solid45 elements. Finally, the material properties had to be entered along with the crystal orientation of the material. Determining the elastic constants in the specimen coordinate system can be achieved by using the following equations taken from Lekhnitskii (1963). The relationship between the stresses in the material and specimen coordinate systems is: {u} = [Q]{cr'} (37) where anx 073x 7'Y '7 ' S. and (38) TZ 7zx 7T T X TX)? TX)? Recall that aL is the cosine of the angle between the x' and x axis, 01 is the cosine of the angle between the x' and y axes, etc. a2 32 2aaa2 2aaa3 2a2a1 #12 2 32 2#3#2 2#1A3 2f2A1 7= 2 Y22 32 27372 2yl73 2721 ArY A22 A73 (3l273 + 3 2) (A213 +A3y1) (A1/2 + )2y1)) ryl, 7,2a2 ,3a3 (Y203 + y3a2) (Y1a3 + y3 1) (y1a2 +2C) a 1 a2l2 a3 3A (a2A3 + 1a3 2) (a1i3 +a3)1) (a1f2 +a2A1) The strain vector transformations are similar. The [Q] matrix is different because there is a factor of two in the relationship between engineering and tensor shear strain components(i.e. yxy = 2 yxy, yzx = 2 yzx, and 7yz = 2 yyz). {'} =[Q ]{e} (310) a12 a2 a a3a2 a1a3 a2a1 2 2 2 Q= 7 2 73 7Y32 7173 7271 (311) 2Ay, 2027, 2A373 (y073 + 73y2) (A,y3 + 37y,) (Ay2 + P7y,) 2y a, 22a,2 273a3 (y2a3 + y3a2) (y1a3 + y3a) (1,a2 + ya ) 2aAP, 2a2,2 2a3A3 (a2,3 + a 2) (afl + a3/1) (af2 + a2/,) Hooke's law for a homogeneous anisotropic body using Cartesian coordinates is: {} = [a ]o} (312) Since this superalloy is FCC and has cubic symmetry, there are only three elastic independent elastic constants needed to define the shear modulus, elastic modulus and Poisson ratio. a11 a12 a12 0 0 0 a12 a11 a12 0 0 0 [, aL12z a12 a111 0 0 0 a = a(313) 0 0 0 a44 0 0 0 0 0 0 a44 0 0 0 0 0 0 a44 1 1 vyx vV all = a44 a12 (314) Ex Gyz Ex Eyy The equation below relates the stresses to the strains in the material specimen coordinate system: {}c = L[a'J{ 1' (316) The elastic constants can be transformed from one coordinate system to another by: 6 6 [a ]= T [airk 1=iria^ a, 31 m=l n=l (317) (i, j=l,2,......,6) With these relationships in hand, the stresses, strains, and elastic constants can be determined in either the material or specimen coordinate systems. Consider a block of the Nibased superalloy shown in the Figure 39. The desired orientation of the crystal lattice is <110> with the x' axis aligned with the <110> direction. One way to look at the problem is to imagine a cylindrical piece of material, shown in the Figure 38, cut from the block. That cylinder now has the desired material properties in the <110> direction. Using equation 317 the elastic constants in any orientation can be determined if one has the direction cosines between the specimen and material coordinate system. ANSYS allows one to input values of the stiffness (units of psi, or Pa) or compliance matrix (units of 1/psi or 1/Pa). Once those values are entered into ANSYSthe finite element model now has the material properties in the desired orientation. <110> y v A x <010> <100> <001>, z, z' Figure 39. Visualization of crystallographic orientation CHAPTER 4 RESULTS OF THE APPLIEDSTRESS MODEL Contour plots of the normal and shear stresses were obtained from the densely meshed region at the midplane of the model (see Figure 41). To make the analysis simpler, the analytical solution does not account for the boundary conditions on the free ends, or surfaces. Since the analytical solution uses a generalized plane strain condition, it is important to extract the results from nodes that are not near the edges, or faces, of the model. midplane of Densely meshed region the model Figure 41. Results were taken from the elements at the midplane How the material orientations are selected is best seen in Figure 42. The cones represent variations in the primary orientation. The inner and outer cones represent a 7.5 degree and fifteen degree variation, respectively. The red dots indicate some possible positions of the zaxis (33 total) that come from the intersection of the xz and yz planes. The third plane is the theta (0) plane, or a rotation about the z axis. Since this material exhibits cubic symmetry, theta only has to vary between 0 and 80 degrees. Following the convention of Arakere and Swanson (2000), theta is broken up into nine degree increments (0, 10, 20, 30, 40, 50, 60, 70, 80). Multiplying the possible values for theta by the 33 primary axis cases, there is a total of 297 material orientations. z 29^^ j 2 27l2n Casing Coo e System correspondence, 2002 i1"' .^^ " 9 Plane Casting Comdinate System Figure 42. Schematic of casting coordinate system Source: G. Duke, personal correspondence, 2002 Case 0 corresponds to the angles A, y, and 0 all being equal to zero as shown in Figure 42. Another example would be case 17. Here, A = 15, y and 0 are zero degrees. Two sets of results will be presented here: case 0 and case 33. Case 33 is not on the chart but A = 30, y and 0 are zero degrees. Comparison of Contour Plots: AppliedStress Model and Analytical Solution The first set of results is a contour plot for case 0. This is where the material and specimen coordinate systems are coincident with each other. Right away, one notices the symmetry between all of the plots. The ANSYS plots have the black background. The analytical plots for the corresponding stress are placed next to the ANSYS contour plots. xx A B Figure 43. Contour plots of the ox stress (case 0). A) ANSYS contour plot. B) Analytical contour plot A B Figure 44. Contour plots of the Gy stress (case 0). A) ANSYS contour plot. B) Analytical contour plot Figure 45. Contour plots of the cz stress (case 0). A.) ANSYS contour plot. B) Analytical contour plot Figure 46. Contour plots of the Txy stress (case 0). A). ANSYS contour plot. B). Analytical contour plot The next set of results is for case 33. Again, the two plots look very similar from a graphical standpoint. Figure 47. Contour plots of the cx stress (case 33). A). ANSYS contour plot. B). Analytical contour plot Figure 48. Contour plots of the Gy stress (case 33). A.) ANSYS contour plot. B). Analytical contour plot A B Figure 49. Contour plots of the cz stress (case 33). A.) ANSYS contour plot. B). Analytical contour plot A B Figure 410. Contour plots of the Txy stress (case 33). A.) ANSYS contour plot. B). Analytical contour plot I00 IOOU rmo Numerical Comparison of the Analytical and AppliedStress Models Contour plots establish a pattern of how the stresses vary, but it is a good idea to also look at the numbers and compare the magnitudes of the stresses. The following plots are for case 33. These results were taken from the elements at the midplane of the model in the densely meshed region (see Figure 411). In Figures 416 and 417, the shear stresses between the two models do not agree very well. The likely cause will be discussed in chapter four. results were obtained from these nodes densely meshed region y mx z midplane Figure 411. Nodes used to extract results Sigma X vs Depth 00 1 0( Depth into Material (in.) Depth into Material (in.) p11 SCZgnm3 r lS S ' m g .rnj$,$ Figure 412. cx stress vs depth Sigma Y vs Depth ' C"J5 L J 1 ~_U Depth into Material (in.) Figure 413. oy stress vs depth Sigma Z vs Depth 0005 001 0i Depth into Material (in.) Sigma ZANSYS mSigma Z Analytical Figure 414. cz stress vs depth 0.00E+00 5.00E+04 1.00E+05 ,. 1.50E+05 2.00E+05 2.50E+05 3.00E+05 0.00E+00 5.00E+04 1.00E+05 ' 1.50E+05 2.00E+05 2.50E+05 3.00E+05 0.OOE+00 5.00E+04 1.00E+05 S1.50E+05 2.00E+05 i 2.50E+05 Tau XY vs Depth Tau XY ANSYS Tau XY Analytical Depth into Material (in.) Figure 415. Txy stress vs depth 1400 1200 1000 . 800  S600 400 200 0 Tyz stress vs depth Tau XZ vs Depth 0 0.005 0.01 Depth into Material (in.) Sigma XZANSYS mTau XZ Analytical 0.015 Figure 417. Txz stress vs depth 20000 15000 10000 CL 5000 01 5000 ) 005 0 01 0 C Tau YZ vs Depth 0 0005 001 0O':I 50 S100 Sigma YZANSYS a 150 Tau YZ Analytical 200 250 Dpeth into Material (in.) Figure 416. ~4~___ CHAPTER 5 CONTACT MODEL Figure 51. Contact model The contact model is a halfcylinder resting on a rectangular block. The same assumptions are in place for the applied stress model: no friction and normal loading only. The volumes are comprised of SOLID45 elements. CONTA174 and TARGE170 elements were used to establish the contact between the halfcylinder and block. The primary aim of running the contact model is to determine the dimensions of the contact patch (the contact halfwidth, a). The applied stress and analytical solutions relied on an assumed value for that number. The contact model will take much longer to run on the computer versus the applied stress or analytical solution, but the computations are necessary to get a handle on the contact patch size. Figure 52. Densely meshed regions of the contact model In Figure 52 above one can see the densely meshed region where the contact elements were placed on the block and where the target elements were applied on the halfcylinder. Figure 53 shows the resulting contact pair created in ANSYS. Figure 53. Contact and target elements Notice the curved shape of the target elements. That is because they are attached to the curved profile of the halfcylinder. In order to create contact models in ANSYS, a contact pair of elements must be createda contact element and a target element. ANSYS has general guidelines as to what line, surface, or volume these elements should be applied. Perhaps the most critical feature is the mesh size. For example, a large target element size and very fine contact element will not work. The sizes of the contact and target elements should be fairly close to one another. It is possible to get a solution to converge, but the results will most likely be incorrect. That is why there is a densely meshed region in both the bottom part of the halfcylinder and on the surface of the block shown in Figures 52. Another important characteristic of this contact model is that displacement boundary conditions had to be used instead of applied forces. Applied displacements are more stable than applied forces in ANSYS. Examination of the Normal Stresses at the Surface The analysis was run for four cases: case 17, 21, 25, and 29. To obtain the size of the contact patch and peak pressure at the surface, the same procedure in chapter three was followed. A group of nodes selected at the midplane of the densely meshed region of the block (see Figure 54). Below are plots of the stress in the ydirection for the selected nodes. Notice that they have an elliptic profile. The normal contact pressure should be elliptical (even though orthotropic materials are involved). If the correct distribution is not apparent, then the size of the elements may be too big and additional refinement is needed so the proper load is applied to the surfaces of the bodies in contact. results were obtained from these nodes densely meshed region midplane Figure 54. Nodes used to obtain the contact patch and peak pressure at the surface Case 17 Sigma Y  Sigma Y mils Figure 55. Ystress for case 17 50000 0 50000 CL 100000 150000 200000 250000 15 20 5 46 Case 21 Sigma Y 50000 0  50000 CL 100000 150000 200000 250000 Sigma Y mils Figure 56. Ystress for case 21 Case 25 Sigma Y 50000 0 50000 C 100000 150000 200000 250000 *Sigma Y mils Figure 57. Ystress for case 25 15 20 ) 5 10 15 20 ^ I Case 29 Sigma Y 50000 0 50000 5 10 15 20 25 50000 5 100000 Sigma Y 150000 200000 250000 mils Figure 58. Ystress for case 29 Table 51 lists the contact patch (taken from the above plots) and also the peak pressures from the four cases that were run in ANSYS. These numbers are compared to peak pressure and contact patch when Equations 32, 33, and 34 are used. The cylinder was composed of an isotropic steel (E = 30 x 106 psi and v = 0.3), the effective radius, R, was 0.1 in. Table 51. Comparison of Hertzian and numerical contact pressures and patch sizes Case 17 21 25 29 ANSYS Contact HalfWidth (in.) 0.003 0.003 0.002 0.002 ANSYS Peak Pressure (psi) 195000 206000 194770 200240 Using Equations 32, 33, and 34 the peak (Hertzian) contact pressure and contact halfwidth are 180100 psi and 0.002774 inches (from the 79 lbf normal load). Looking at the table above it would appear the true contact halfwidth is between 0.002 and 0.003 inches. A finer mesh would help to determine if that is the case. The consequence of a better mesh (more elements) is computational cost. ANSYS does have a limit on the number of elements that can be used. Approaching that limit would tax ANSYS and the abilities of the computer to the point that one run could take many hours to finish. Comparison of Contact Model Stresses to the Analytical Solution The following plots compare stresses between the contact model and analytical solution taken at the midplane straight down (y direction) into the rectangular block. The stresses between the applied stress model and analytical solution were already compared and the results below are similar in that the cx, y oz, and Txy stresses match well but the remaining two shear components do not. It should be pointed out that the Txz and Tyz are not always zero for the analytical solution. As the orientation changes, those two stress components will not be zero. For case 21, Txz is nonzero and in case 29 Tyz is nonzero. Even though Lehnitkskii (1963) uses what he refers to as a generalized plane strain condition, there is still shear coupling in the z direction. The elements in the densely meshed regions of both models had finer meshes in the x and y directions than the z direction. Figure 59 illustrates the size of the brick element used in the contact model. Brick elements of a slightly smaller size were used in the applied stress model (the zdimension was 0.005 in.). One can see that there is a large difference between the length in the x and y directions and the z direction of the brick element. The difference between the z dimensions was only 0.002 inches, but that is enough to throw off at least one of the shear stress components. For example, the Txy component from the contact model does not agree with the analytical solution. The issue here is probably one of convergence. One way that might help is to use a finer mesh. Four of the six stress components from the applied stress model match up very with the analytical solution and that is most likely due to the smaller element size that model uses. 0001 in 0. 001 in. 4 0.07n. 0.007 m. Ly x z z .001 in. Figure 59. Dimensions of the brick element used in the contact model The reason why these elements had such dimensions is computational cost. A smaller dimension in the z direction would have resulted in more elements. Ideally one would like to have many cube elements or small tetrahedrons. But that would have resulted in contact models having a large number of elements. A model that big would take a very long time to solve and the results files (where the stresses are stored) also would have been large. Sigma X vs Depth (case 17) 0.00E+00 5.00E+04 'C 1.00E+05 1.50E+05 2.00E+05 Sigma XANSYS  Analytical mils Figure 510. Sigma X comparison case 17 50 Sigma X vs Depth (case 25) 0.00E+00 5.00E+04 C. 1.00E+05 1.50E+05 2.00E+05  Sigma XANSYS Analytical mils Figure 511. Sigma X comparison case 25 Sigma Yvs Depth (case 17) Sigma Y ANSYS  Analytical mils Figure 512. Sigma Y comparison case 17 0.00E+00 5.00E+04 S1.00E+05 1.50E+05 2.00E+05 5 10 51 Sigma Y vs Depth (case 25) O.OOE+00 5.00E+04 'i. 1.00E+05 1.50E+05 2.00E+05  Sigma Y ANSYS Analytical mils Figure 513. Sigma Y comparison case 25 Sigma Z vs Depth (case 17) 0.00E+00 2.00E+04 4.00E+04 6.00E+04 i 8.00E+04 1.00E+05 1.20E+05 1.40E+05 1.60E+05  Sigma ZANSYS  Analytical mils Figure 514. Sigma Z comparison case 17 ) 5 10 1 s  S5 1 r 52 Sigma Z vs Depth (case 25) O.OOE+00 2.00E+04 4.00E+04 6.00E+04 '. 8.00E+04 1.00E+05 1.20E+05 1.40E+05 1.60E+05  Sigma ZANSYS  Analytical mils Figure 515. Sigma Z comparison case 25 Tau XY vs Depth (case 17) Tau XY ANSYS  Analytical mils Figure 516. Tau XY comparison case 17 S5 1 /____ 20000 10000 0 10000 1 ' 20000 30000 40000 50000 60000 5 0~ * 53 Tau XY vs Depth (case 25) 0 5000 10000 15000 '. 20000 25000 30000 35000 40000 mils Figure 517. Tau XY comparison case 25 Tau YZ vs Depth (case 17) 3000 2000 1000 0 1000 2000 Tau XY ANSYS  Analytical  Tau YZ ANSYS  Analytical mils Figure 518. Sigma YZ comparison case 17 1) 5 54 Tau YZ vs Depth (case 25) Tau YZANSYS  Analytical mils Figure 519. Tau YZ comparison case 25 Tau XZ vs Depth (case 17)  Tau XZ ANSYS  Analytical mils Figure 520. Tau XZ comparison case 17 3000 2500 2000 1500 '. 1000 500 0 500 1000 5 1 20000 15000 10000 5000 0 5000 S1 10 Tau XZ vs Depth (case 25) Tau XZANSYS Analytical mils Figure 521. Tau XZ comparison case 25 20000 15000 10000 5000 0 5000 10 1 CHAPTER 6 DEFORMATION MECHANISMS Plastic deformation in metals results from a phenomenon called slipthat is the movement of dislocations throughout the material. Typically slip occurs over planes with high atomic density in the closepacked directions. When slip occurs atoms are being forced to relocate to new lattice positions within the crystal structure. Atoms located on densely packed planes will have a shorter distance to travel when compared to atoms situated on other, less densely packed, planes. This is not always the case, and at higher temperatures dislocations can move along planes with lower atomic densities because of two things: cross slip and climb. Cross slip (phenomenon that entails dislocations moving onto an intersecting slip plane) is a thermally activated process. Climb (also a thermally activated process) is associated with diffusion (movement of atoms into vacancies, etc) (Hummel 1998). Since the material in question is a Nickel based superalloy, the planes with the highest atomic density are the octahedral planes. A slip system requires the definition of the slip plane and the slip direction. For example, (111)[01 1] is a slip system where (111) is the slip plane and [0 1 1] is the slip direction. The slip direction always lies in the slip plane (Figure 62). The aforementioned notation comes from the Miller indices. Let us look at the (111)[0 1 1] system graphically (Figure 61). The (111) plane is encased by the cube and the vector <011> is the slip direction drawn from the origin. With the four octahedral planes there are 12 closepacked directions the planes can move and also another set of 12 directions in the nonclosepacked slip direction. Thus, 57 there are 24 slip systems on the octahedral planes. At higher temperatures, slip can occur on the cube planes. These planes are not as densely packed as the octahedral planes and they, along with the octahedral planes, are shown below in Figures 62 and 63. Counting the cube slip systems, there are a total of 30 slip systems (Dame and Stouffer 1996). Figure 61. Example slip system DIi. I 010 eoPYd]i.z Pla 2 Pnay 0, t, T' Secondaiv T11, T17 t11 Pl, 3 010 iD . Sec..day Oe 010 Figure 62. Four octahedral planes and slip directions. Source: Modified from Dame and Stouffer 1996 11W Flamel PHM 001010 100 001 010 Plane 3 S100 1 1001 1 Figure 63. Cube slip planes and slip directions. Source: Modified from Dame and Stouffer 1996 With the 30 slip systems defined one can now establish equations to acquire local stresses and strains within the material. From Dame and Stouffer (1996) we have S1 0 1 1 0 1 S0 1 1 1 1 01 1 S1 1 0 0 1 1 4 1 0 1 1 0 1 co 5 1 1 0 0 1 1 O r6 1 0 1 1 1 1 0 ^ (61) r7 6 1 1 0 0 1 1 ,y T8 0 1 1 1 1 0 o" r9 1 0 1 1 0 1 yz r10 0 1 1 1 1 0 11 1 0 1 1 0 1 12 1 1 0 0 1 1 59 13 1 2 1 1 2 1 14 2 1 1 1 1 2 15 1 1 2 2 1 1 16 1 2 1 1 2 1 o, 17 1 1 2 2 1 1 o 718 1 2 1 1 1 1 2 2) I (62) 19 3= 1 1 2 2 1 1 C ,X 20 2 1 1 1 1 2 o 21 1 2 1 1 2 1 o z22 2 1 1 1 1 2 23 1 2 1 1 2 1 z24 1 1 2 2 1 1 The six slip systems below are for the cube slip planes. '25 0 0 0 1 1 0 oxx S26 0 0 1 1 0 S27 0 0 0 1 0 1 C ( '28 V2 0 0 0 1 0 1 x 29 0 0 0 0 1 1 O 30 0 0 0 0 1 1 OC cI'c30 are the local shear stresses associated with that slip system. For example, TI is associated with slip system number one. If any of these values ('lC30) are above a threshold stress called the critical resolved shear stress, then that slip system will become active. It is possible for more than one slip system to become active depending on the orientation and loading of the material. It might be of some use to investigate how the local shear stresses ('C30) vary in the material. In order to calculate these local shear stresses, all six stress components (Ox, O, Oz, C xy,'zx, Tyz) are needed for the equations listed above. Then at each node all of the local shear stresses can be tabulated. Table 61. Slip Systems Slip Number Slip Plane Octahedral Slip a/2<110>{111} 1 [1 1 1] 2 [ 1 1] 3 [1 1 1] 4 [1 11] 5 [1 11] 6 [1 1 1] 7 [1 1 1] 8 [1 1 1] 9 [1 1 1] 10 [1 1 1] 11 r1 1 1] 12 [1 1 1] Octahedral Slip a/2<112>{111} 13 [1 1 1] 14 [1 1 1] 15 [1 1 1] 16 [1 1 1] 17 [1 1 1] 18 [1 1 1] 19 [1 1 1] 20 [1 1 1] 21 [1 1 1] 22 [1 1 1] 23 [1 1 1] 24 [1 1 1] Cube Slip a/2<110>{100} 25 [1 0 0] 26 [1 0 0] 27 [0 1 0] 28 [0 1 0] 29 [0 0 1] 30 [0 0 1] Source: Modified from Dame and Stouffer 1996 L Slip Direction 12 Primary Slip Directions [1 0 1] [0 1 1] [1 1 0] [1 0 1] [1 1 0] [0 1 1] [1 1 0] [0 1 1] [1 0 1] [0 1 1] [1 0 1] [1 1 0] 12 Secondary Slip Directions [1 2 1] [2 1 1] [1 1 2] [1 2 1] [1 1 2] [2 1 1] [1 1 2] [2 1 1] [1 2 1] [2 1 1] [1 2 1] [1 1 2] 6 Cube Slip Directions [0 1 1] [0 1 1] [1 0 1] [1 0 1] [1 1 0] r1 1 0] 61 Once again, nodes from the densely meshed region at the midplane of the rectangular slab were used to obtain the required stresses. Local shear stresses from the primary planes ('1'c12) are given below for cases 17 and 25. In all cases the peak stresses occurred a few mils below the surface, and then decrease as a function of depth (y direction) into the material. Also, in comparing the local shear stresses one notices how similar the contours are between the two cases. 5 10 15 A Figure 64. Tau 1 contour plot. A.) 0 5 10 15 20 A Figure 65. Tau 2 contour plot. A.) Case Case 17. B.) Case 25 5 10 15 20 17. B.) Case 25 2 0 0 5 10 15 20 dd A Figure 66. Tau 3 contour plot. A.) Case 17. 8 62 1 44 4 1. 0 0 5 10 15 20 ee A Figure 67. Tau 4 contour plot. A.) Case 17. 4 a 0 0 5 dd B.) Case 25 8 6 4 2 0 0 5 ee B.) Case 25 10 15 20 B 4 4  0 5 10 15 20 ff A Figure 68. Tau 5 contour plot. A.) Case 17. :1 1 0 5 10 15 20 gg A Figure 69. Tau 6 contour plot. A.) Case 17. 0 5 ff B.) Case 25 8 4 0 5 gg B.) Case 25 10 15 20 B 10 15 20 B 64 8  0 0 0 5 10 15 20 0 5 hh hh A Figure 610. Tau 7 contour plot. A.) Case 17. B.) Case 25 4.0 a 0 5 10 15 20 A Figure 611. Tau 8 contour plot. A.) Case 17. 10 15 20 B 4 4. 0 5 B.) Case 25 15 20 8 0 8 6 6 4 4 4 2 1 2 0 0 0 0 5 10 15 20 0 5 ii ii A Figure 612. Tau 9 contour plot. A.) Case 17. B.) Case 25 8 6 o 0 0 10 15 20 0 5 kk kk A Figure 613. Tau 10 contour plot. A.) Case 17. B.) Case 25 10 15 20 B 15 20 8 6 41 21 0 0 5 10 15 20 11 A Figure 614. Tau 11 contour plot. A.) Case 17. 0 5 10 15 20 mm A Figure 615. Tau 12 contour plot. A.) Case 17. 8 6 4 4 2 0 0 0 5 2 11 0 0 5 mm B.) Case 25 10 15 20 B CHAPTER 7 LITERATURE REVIEW This section gives a discussion of the limited work done in the area of anisotropic contact mechanics. Relevant work dealing with analytical and numerical solutions is highlighted. Sinclair and Beisheim (2002) developed a procedure to analyze stresses of dovetail attachments in turbine blades using a submodeling routine. These stresses are difficult to obtain because of the nature of the contact problem (conforming) and also the large magnitude and gradient of the contact stresses. These parameters require a high number of elements to achieve accurate results. A submodeling routine was implemented in order to get converged peak stresses without considerable computational cost. The first step is to run a global analysis that includes the highly stressed contact region. Three grids labeled coarse, medium, and fine are used to mesh the model. The FE model is run with each of the three grids and the results from each grid are compared to check for convergence. The stress component in question, such as a peak stress, is said to have converged if f max mmax < (71) < Es (71) 0" max where Es is the error level (0.01 is excellent, 0.05 is good, and 0.1 is satisfactory). The superscripts f and m denote fine and course grids. The next step is the submodeling routine. This is where the highly stressed region, or submodel, is broken out and analyzed separately using the information obtained from the global analysis. Distances, or boundaries, away from peak stress component are chosen such that the stress component in question has converged. Then the submodel is divided up with its own coarse, medium, and fine grids. Displacements were used to check for convergence since they converge faster than stresses. Boundary displacements are said to be converging if Uc _Um > zm Ulf r (72) uf tU, U1i (72) where i denotes the node in question. Cubic spline equations were fit to the nodal displacements. This was done so even more displacements could be ascertained between nodes for all of the grids. Using the fitted splines, the three grids are run and checked to see if the submodel boundary conditions have converged using two expressions that calculate the boundary condition error and the discretization error. The latter is built in to all FE analysis. Even though Sinclair and Beisheim (2002) claim that their numerical approach requires 'modest' computing resources, a very high number of elements was used for their analysis. Table 71. Number of elements used Grid Spatial Contact Total Global Course 2127 478 2614 Medium 14990 1943 16993 Fine 108870 7765 116635 Outer Ring Course 540 0 540 Medium 4140 0 4140 Fine 34147 0 34147 Submodel Coarse 9216 1152 10368 Medium 73728 4608 78336 Fine 589824 18432 608256 Source: Modified from Sinclair and Beisheim (2002) On the average, contact models in this work, composed of roughly 30,000 elements, took approximately two hours for the computer to complete its run for each crystal orientation. Over 600,000 elements were used in the submodeling routine to get accurate results in the contact region. Only had normal tractions applied to the model, just like the models used in this analysis. In chapter three the applied stress results were compared to analytical solution. The normal stress components agreed fairly well, but the shear stresses did not. A potential source of this problem was the size of the elements used, or the mesh size. Smaller elements, and consequently bigger FE models, are needed to get the correct stressesas shown in Table 71. Using that many elements resulted in a total estimated error of 2.7%. When looking at contact problems (such as those involving indenters) the contact patch is a function of the material properties, geometry, and load. In most cases, one has a handle on what the geometry of the indenter and halfspace is and also what external force is being applied. When isotropic materials are used the contact patch is only dependant on two material properties: such as the elastic modulus and the poisson ratio. Anisotropic materials are different in that the contact patch is related to an equivalent, or composite, modulus since the material properties are direction dependent. Vlassak et al. (2003) looked at calculating the indentation modulus of anisotropic materials. Analytical solutions were developed for indenters of arbitrary shape being pressed into anisotropic halfspace. A theorem stating that the solution to this contact problem is one that maximizes the load on the indenter for a certain depth. One of the principle assumptions was to assume the contact patch had a certain shape, then the best approximate solution exists in the RayleighRitz sense. For axisymmetric indenters, a limited family of Green's functions, using the first term in the Fourier expansion, can be used to obtain a much simpler solution. This solution is denoted as the "equivalent isotropic solution." First, the Green's function used to obtain displacements at the surface of the half space is defined. The next step is to determine the load on the punch, or indenter. With the load function available, it can be differentiated and the modulus of the halfspace is found. The Green's function gives the displacement in the direction of a load applied to the halfspace. w(y) = I, akk :1 a (73) Equation 74 is the force that will establish the elliptical contact area. F(A) u(x, y)dxdy (74) a2 b2 1 where po = and u(x,y) is the displacement from the indenter. 7;iba(e, p) Now the solution for a conical indenter is produced so the modulus for this case can be found. For a conical indenter the displacement is u(x,y) =U Cx2 + y (75) where U is the rigid body displacement and C is the cotangent of the cone angle. Using Equation 74 the force on the indenter is now U CCx + y2 dxdy F(A) JJf (76) r a2 b2 After a change of variables and also switching to polar coordinates the force is 2Ua Ca2E(e) F(a, b, () = (77) The resulting contact stiffness for a conical indenter is sdF 21 S = 1(78) dU 4 a(e,p)(1e2)1/4 where U is the approach of the two bodies and A = nab. The modulus of the halfspace is eqv = 1/4 (79) a(e, p)(1 e2)1 A less complicated expression for conical indenters can be found by assuming the contact patch is circular. Now the orientation angle, p, is arbitrary and the eccentricity drops out a(O,() = h(O)dO = 7h (710) 0 ho is the first term of the Fourier expansion of the function h(O), and modulus now becomes qM (711) zho which is the same solution for a flat circular punch for an isotropic halfspace; the equivalent isotropic solution. The modulus in Equation 711 is compared to the modulii one gets from when a sphere and cone are the indenters and the Green's function is not truncated using a Fourier expansion. When the modulii are plotted against crystal orientation, the modulii of the cone and sphere do not vary much (the difference between the two is less than 0.03%). Using the equivalent isotropic solution, the difference is approximately 0.1%. Green and Zerna (1954) looked at the twodimensional contact problem for anisotropic materials. Their work is somewhat confusing because they retain the notation used for isotropic materials. Their solution is not applicable for coupled inplane and antiplane deformation. Fan and Keer (1994) reexamine the twodimensional contact problem using the analytic function continuation approach. They arrive at a compact solution that may be applied to other problems such as examining subsurface stresses in contact problems involving sliding and sticking. Willis (1966) examined the Hertzian contact for anisotropic halfspaces. He assumed the contact patch is elliptical and the contact pressure will be of the form Hertz used when he was looking at this problem for isotropic materials. Willis does not use potential theory to solve this problem. He used Fourier transforms and the solution lies in solving contour integrals. Evaluating those integrals can be done numerically. CHAPTER 8 CONCLUSIONS The 2D analytical solution developed for subsurface stresses in anisotropic half spaces is an efficient and practical method for performing design iterations. Contact dimensions and pressure distributions obtained by Hertzian isotropic assumptions are a good approximation for anisotropic contacts also. The 3D FE model, or applied stress model, for stress analysis in the anisotropic substrate (single crystal), using contact dimensions and pressure profile from Hertzian solutions is an excellent representation of the contact problem. This is an accurate and efficient technique for determining subsurface stresses for anisotropic contacts. The 3D contact model for anisotropic substrates, while numerically intensive, is the best way to solve problems involving complex geometries. Just like the applied stress model, effects of crystal orientation were included. Additional mesh refinement might produce more accurate stress solutions in the edge of contact. Overall, most of the results had excellent agreement with the applied stress model and analytical solution. Contour plots for resolved shear stresses on slip planes are presented in the contact region. These results are required for determining fatigue life (crack initiation) in the substrate. CHAPTER 9 RECOMMENDATIONS FOR FUTURE WORK Effects of contact friction and tangential forces must be incorporated in the FE models. The analytical solution can already account for tangential traction forces of the contact. Relationships between shear stress amplitude on the slip planes and cycles to fatigue crack nucleation need to be determined. To effectively model the fatigue damage process, the ability to evaluate plastic stress and strain amplitudes on the contact is required. This will entail developing constitutive relations for crystal plasticity, which will prove to be a complex endeavor. LIST OF REFERENCES ANSYS Elements Reference, ANSYS Release 5.6; ANSYS, Inc. November 1999. Arakere N, Swanson G. Effect of crystal orientation on analysis of singlecrystal, Nickel base turbine blade superalloys, National Aeronautics and Space Administration Technical Publication 2000; Feb: 610. Beisheim JR, Sinclair GB. Threedimensional finite element analysis of dovetail attachments. American Society of Mechanical Engineers Turbo Expo. 2002 Jun 36; Amsterdam, the Netherlands. Dame L, Stouffer D. Inelastic deformation of metals. New York: John Wiley and Sons; 1996. p. 346, 387435. Donachie M, Donachie S. Superalloys. Materials Park, Ohio: American Society for Metals International, 2002. p. 2539. Fan H, Keer LM. Twodimensional contact on an anisotropic halfspace. Transactions of the American Society of Mechanical Engineers 1994; 61: 250255. Green AE, Zerna W. Theoretical Elasticity. Oxford, UK: Clarendon Press; 1954. Hummel R. Understanding materials science. New York: Springer; 1998. p.74101, 102 122. Johnson KL. Contact mechanics. Cambridge, UK: Cambridge University Press; 1985. p. 84106 Lekhnitskii SG. Theory of elasticity of an anisotropic body, San Francisco: HoldenDay; 1963. Magnan S. Threedimensional stress fields and slip systems for single crystal superalloy notched specimens, Master's Thesis, University of Florida, Gainesville 2002, p. 4152. Sims CT, StoloffNS, Hegel WC, editors. Superalloys II. New York: John Wiley and Sons; 1987. Timoshenko SP, Goodier JN. Theory of elasticity. Tokyo: McGrawHill; 1983. p. 1534. 76 Vlassak JJ, Ciavarella M, Barber JR, Wang X. The indentation modulus of elastically anisotropic materials for indenters of arbitrary shape. Journal of the Mechanics and Physics Solids 2003;51: p. 17011721. Willis JR. Hertzian contact of anisotropic bodies. Journal of the Mechanics and Physics Solids 1966; 14: 163176. BIOGRAPHICAL SKETCH Born on February 18, 1979 in West Chester, Pennsylvania I moved to Naples, Florida when I was four. I graduated with a B.S. in mechanical engineering from Villanova University and am now seeking a doctorate in mechanical engineering. I plan on becoming a professor someday. 