- Email: [email protected]
Three dimensional computational analysis of fatigue crack propagation in functionally graded materials
Computational Materials Science 52 (2012) 246–252
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Three dimensional computational analysis of fatigue crack propagation in functionally graded materials Baris Sabuncuoglu a,⇑, Serkan Dag a, Bora Yildirim b a b
Department of Mechanical Engineering, Middle East Technical University, Ankara 06800, Turkey Department of Mechanical Engineering, Hacettepe University, Ankara 06800, Turkey
a r t i c l e
i n f o
Article history: Available online 30 June 2011 Keywords: Functionally graded materials Fracture mechanics Fatigue crack propagation Finite element method Stress intensity factors
a b s t r a c t This article proposes a new finite elements based three dimensional method developed to study the phenomenon of fatigue crack propagation in functionally graded materials (FGMs). The particular problem examined in detail is that of an initially-elliptical crack located in a functionally graded medium, subjected to mode I cyclic loading. The crack is modelled by employing three dimensional finite elements; and the stress intensity factors (SIFs) around the crack front are computed by the application of the displacement correlation technique (DCT). Fatigue crack propagation calculations are based on the Paris– Erdogan crack growth law. The developed procedure makes it possible to generate the crack front profiles corresponding to given numbers of loading cycles. Proposed methods are validated by examining the behaviours of both stationary and propagating cracks. Numerical analyses conducted for an initiallyelliptical crack lying in a graded medium demonstrate that the proposed model can effectively capture the evolution of the crack front morphology. Hence, the methodology presented in this article can be used along with fracture toughness data to assess the fatigue lives of functionally graded components containing cracks that possess arbitrary front profiles. Ó 2011 Elsevier B.V. All rights reserved.
1. Introduction Functionally graded materials (FGMs) are multiphase composites that possess spatial variations in the volume fractions of the constituents. These variations are intentionally introduced so as to take advantage of desirable properties of the constituents, such as the high thermal resistance of ceramics and superior fracture toughness of metals. FGMs have been used successfully in a number of technological applications including thermal barrier coatings, solid oxide fuel cells, high performance cutting tools, and biomedical materials. The spatial variations in the volume fractions of the constituents of FGMs require that certain physical properties be represented by continuous functions of spatial coordinates for analytical and computational purposes. Fracture mechanics of FGMs has been studied extensively in order to develop an understanding of the influence of the continuous gradation on the cracking behaviour. Although great bulk of this work is based on the two-dimensional fracture mechanics theory, in recent years there have been attempts to build computational models for three dimensional fracture mechanics analysis of functionally graded materials. Among the computational methods ⇑ Corresponding author. Present address: Interuniversity Microelectronics Center (IMEC), Kapeldreef 75, B-3001 Leuven, Belgium. Tel.: +44 (0) 1509 227520. E-mail addresses: [email protected] (B. Sabuncuoglu), [email protected] (S. Dag). 0927-0256/$ - see front matter Ó 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2011.06.010
proposed for three dimensional fracture analysis of FGMs, we can mention the domain integral method [1], the displacement correlation technique (DCT) [2,3], the method of enriched finite elements [4], the interaction integral approach [5], and the boundary domain integral method [6]. Note that in all these studies on three dimensional computational fracture mechanics of FGMs, the focus is on stationary cracks, and there seems to be no work on fatigue behaviour of FGMs. Present study is directed towards developing a three dimensional computational analysis technique to study the phenomenon of fatigue crack propagation in functionally graded structures. Fatigue crack growth computations are carried out by considering an embedded crack in the examined functionally graded medium. The initial front profile of the embedded crack is assumed to be in the shape of an ellipse. Elliptical cracks or voids are known to be the source of various structural failures especially under repeated loading conditions. Fatigue crack growth calculations involving elliptical cracks require the computation of the stress intensity factors (SIFs) around the crack border. We use the displacement correlation technique [2,3] in conjunction with the finite element method to evaluate the stress intensity factors. The subsequent front profiles of the initially-elliptical cracks subjected to cyclic loading are generated by adopting the fatigue crack growth analysis methods proposed by Paris and Erdogan [7], and Joseph and Erdogan [8]. The main outcome of this study is the finite elements based three dimensional computational analysis tool, which can be
B. Sabuncuoglu et al. / Computational Materials Science 52 (2012) 246–252
247
used to predict the fatigue lives of functionally graded components. Presented numerical results illustrate the influence of material property gradation on the fatigue propagation characteristics of cracks located in cyclically loaded ceramic–metal FGMs.
2. The computational procedure In this section, we provide the details of the three dimensional computational analysis methodology developed to examine fatigue crack growth in functionally graded materials. The geometry of the particular problem considered is depicted in Fig. 1. An elliptical crack resides in an elastic functionally graded medium, whose material properties are functions of the x-coordinate. The major and minor radii of the elliptical crack are denoted by c and a, respectively. The crack is loaded through the normal stress r applied over the remote bounding planes. Since, the material is assumed to possess property gradation in x-direction, one-quarter of the medium is modelled by considering symmetries about xyand xz-planes. The finite element model for the problem depicted in Fig. 1 is constructed by using the general purpose finite element analysis software ANSYS [9]. Fig. 2 shows the quarter model generated and close-up view of the region around the crack front. The dimensions were taken such that the bounding planes would have no effect on the behaviour of the elliptical crack. In order to simulate the square-root singularity of the strain field near the crack border, the volume around this border is modelled by utilizing collapsed 20-noded isoparametric brick elements, whereas the domain excluding this volume is discretized by means of brick and tetrahedral elements. Smooth spatial variations of the material properties of the graded medium are taken into account in the finite element model by computing the properties of each finite element at its centroid. The properties are uniform over each element; however they vary from element to element according to the position of the centroid. This is the so called homogeneous finite element approach and known to lead to highly accurate numerical results so long as there is appropriate mesh refinement in the finite element model [2,3]. If the displacement field for the cracked functionally graded medium is known, the mode I stress intensity factor distribution around the crack front can be evaluated by using one of the available computational methods. In the present study, we use the displacement correlation technique in the evaluation of the SIFs. In this technique, the expression of the stress intensity factor is derived by considering a point P on the crack front. Fig. 3 shows such
Fig. 2. (a) The quarter model and (b) the crack front region.
a point P lying on the border of the crack. The origin of a Cartesian coordinate system is also located at point P. This coordinate system comprises tangential (t), normal (n), and binormal (b) directions. Fig. 4 depicts the deformed shape of the crack in the normal plane passing through point P. The normal plane is formed by n and b
Fig. 1. (a) An elliptical crack in an infinite volume and (b) a point P on the elliptical crack front.
248
B. Sabuncuoglu et al. / Computational Materials Science 52 (2012) 246–252
crack front is specified as Dq0. The number of cycles corresponding to this increment is computed using
DN 0 ffi
Dq0 : C 0 DK nI00
ð5Þ
The increments at other points around the crack front are then calculated by substituting Eq. (5) into Eq. (4): n
Dqi ffi Fig. 3. Crack front and the local coordinate system at point P.
C i DK Iii Dq0 ; C 0 DK nI00
ð6Þ
where the subscript i stands for a generic point on the crack front. A sufficiently small value needs to be used for Dq0 to be able to obtain accurate results through Eqs. (5) and (6). After calculating the crack growth increments for a specified number of points on the crack front, the new crack front profile is generated by connecting the new locations by splines. The procedure is then repeated to form the crack front profile corresponding to the next specified increment.
3. Verification of the developed procedure Fig. 4. Deformed shape of the crack surface in the normal plane.
axes. The nodes 1–3 in Fig. 4 are considered to be at the edge of a collapsed 20-noded brick element. Under mode I loading conditions, the asymptotic displacements of these nodes in b-direction can be found from [3]
ub ðnÞ ¼
4ð1 ðmP Þ2 Þ EP
rffiffiffiffiffiffiffi jnj KI; 2p
ð1Þ
where KI is the mode I stress intensity factor; EP and mP are the modulus of elasticity and the Poisson’s ratio computed at point P, respectively; and the superscript P implies that the corresponding material property should be calculated at point P. The mode I SIF at P is defined as [2,3]
) ub ðjnjÞ KI ¼ lim pffiffiffiffiffiffi : jnj 4ð1 ðmP Þ2 Þ jnj!0 pffiffiffiffiffiffiffi P 2pE
(
ð2Þ
Then, using Eqs. (1) and (2) and the displacements of the nodes 2 and 3, it is possible to express the mode I stress intensity factor in the following form:
# 3=2 3=2 R3 ub2 R2 ub3 pffiffiffiffiffipffiffiffiffiffi : KI ¼ R2 R3 ðR3 R2 Þ 4ð1 ðmP Þ2 Þ pffiffiffiffiffiffiffi P 2pE
"
ð3Þ
R2 and R3 in Eq. (3) are the distances of nodes 2 and 3 with respect to the origin. The fatigue crack growth calculations for the configuration shown in Fig. 1 are based on the mathematical models proposed by Paris and Erdogan [7], and Joseph and Erdogan [8]. Suppose that the stress applied to the graded medium r varies cyclically in time with a constant stress range Dr and a constant mean stress. Then, the growth rate at a point on the crack front can be calculated using Paris–Erdogan law [7,8], which is expressed as:
dq ¼ C DK nI ; dN
ð4Þ
where q is the radius of curvature at the considered point on the crack front, N is the number of cycles, DKI is the stress intensity factor range, and C and n are material parameters. In the finite element implementation, first, the increment for a particular point on the
We validated the methods described in the previous section on problems involving stationary and propagating cracks. For the stationary cracks, the application of the displacement correlation technique is verified by considering two different problems, the first of which is the problem of an elliptical crack in a homogeneous medium, and the second is that of a surface crack in a functionally graded plate. As for the propagating cracks, the validation study is conducted by examining fatigue growth in a homogeneous medium. In what follows below, we present the results generated for these validation problems. By considering the elliptical crack problem shown in Fig. 1, Irwin [10] derived a simple closed form expression for the mode I stress intensity factor. In this derivation, the material was assumed to be homogeneous. We first focus on this problem solved by Irwin [10] and provide comparisons of the mode I SIFs computed by the DCT developed in the present study to those obtained through Irwin’s formula. These comparisons are given in Table 1. In the calculation of these results, the medium shown in Fig. 1 is assumed to be homogeneous and made of Ti–6Al–4V. The material properties are taken from Refs. [3,11] and provided in Table 2. The minor and major radii of the crack are respectively set as 0.5 cm and 1.5 cm. The normalized mode I stress intensity factor provided in Table 1 at various locations around the crack front is defined as
KI K In ¼ pffiffiffiffiffiffi : r pa
ð7Þ
Examining the results given in Table 1, it is seen that the normalized mode I stress intensity factors computed by the DCT are in excellent agreement with those obtained through Irwin’s for-
Table 1 Comparisons of the normalized stress intensity factors computed by the displacement correlation technique to those obtained by the formula given by Irwin [10]. 2Ø/p
Irwin [10]
DCT
% Difference
1 (minor axis) 0.875 0.75 0.5 0.375 0.25 0.125 0 (major axis)
0.8979 0.8902 0.8671 0.7752 0.7075 0.6293 0.5540 0.5184
0.8977 0.8898 0.8675 0.7755 0.7068 0.6270 0.5523 0.5180
0.02 0.05 0.04 0.03 0.10 0.37 0.31 0.09
249
B. Sabuncuoglu et al. / Computational Materials Science 52 (2012) 246–252 Table 2 Properties of Ti–6Al–4V [3,11]. Modulus of elasticity (GPa)
Poisson’s ratio
Paris–Erdogan law coefficient (C)
Paris–Erdogan law exponent (n)
105.8
0.298
5.2 1012
3.17
Table 3 Comparisons of the normalized mode I stress intensity factors computed by the displacement correlation technique to those given by Walters et al. [1]. a/h = 0.5, a/ c = 2, E(h)/E1 = 5.
mula. Hence, it is possible to compute the mode I stress intensity factors to within a high degree of accuracy by employing the displacement correlation technique developed in the present study. To be able to demonstrate the applicability of the developed method for functionally graded materials, we provide additional comparisons by considering the static surface crack problem studied by Walters et al. [1]. The problem geometry is shown in Fig. 5. The figure depicts a semi-elliptical surface crack located in an isotropic FGM plate, that is under the effect of remote normal stress r at the ends y = ±L. The modulus of elasticity of the plate is represented by an exponential function, and given by
EðxÞ ¼ E1 expðbxÞ;
ð8Þ
where b is a nonhomogeneity constant. Poisson’s ratio is assumed to be constant and equal to 0.25. L and B are taken as sufficiently large by Walters et al. [1]; and as a result the stress field around the crack front is not affected by the boundaries at y = ±L and those at z = ±B. Table 3 tabulates the normalized mode I stress intensity factors computed by using the displacement correlation technique developed in the present study and those given by Walters et al. [1], which are evaluated by means of a three dimensional domain integral method. The normalized mode I SIF is defined by
Fig. 5. (a) Surface cracked FGM plate problem considered by Walters et al. [1] and (b) plan view of the semi-elliptical surface crack.
2Ø/p
Walters et al. [1]
DCT
% Difference
1 (major axis) 0.875 0.75 0.5 0.375 0.25 0.125 0 (minor axis)
0.580 0.595 0.629 0.679 0.685 0.677 0.656 0.596
0.587 0.600 0.631 0.686 0.693 0.684 0.664 0.613
1.21 0.84 0.32 1.03 1.17 1.03 1.22 2.85
KI K In ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi :
r
pa 1þ1:464ðc=aÞ1:65
ð9Þ
The results obtained by the developed displacement correlation technique are in very good agreement with those provided by Walters et al. [1], which is indicative of the high level of accuracy achieved through the application of this technique. Percent difference at the point where the crack front intersects the free surface, i.e. at / = 0, is relatively larger compared to the percent differences computed at the other points. The degree of singularity at this point is different from 0.5 and depends on the Poisson’s ratio [12]. Both in the article by Walters et al. [1] and in the present study, the special singular behaviour at the free surface is not considered, and the results are obtained by assuming that there is a square-root singularity at all points around the crack front including the point of intersection at the free surface. The change in the singularity at this single point could be the main reason behind the difference between the results computed for / = 0. However, although the percent difference is relatively larger at this point, it is less than 3% and the agreement is still acceptable. The finite element mesh used in the calculation of the results tabulated in Table 3 is shown in Fig. 6. Next, we consider fatigue crack propagation in a homogeneous medium subjected to cyclic loading for further verification. The initial geometry of the crack is as depicted in Fig. 1; and the medium is assumed to be subjected to a cyclically varying stress with a constant stress range given as Dr = 100 MPa and a constant mean stress of 50 MPa. The material is homogeneous Ti–6Al–4V. The initial minor and major radii of the elliptical crack are taken as 0.5 cm and 1.5 cm, respectively. Fatigue growth behaviour of this crack is examined by using two different methods, the first of which is the method elucidated in the previous section. The second method is based on the assumption that the initially-elliptical crack retains
Fig. 6. Finite element model used in the solution of the semi-elliptical crack problem considered by Walters et al. [1].
250
B. Sabuncuoglu et al. / Computational Materials Science 52 (2012) 246–252
its elliptical form during propagation. Hence, in the second method it is possible to generate the crack front profile at each step by considering only two points on the crack border. In our calculations pertaining to this second method, we used the points on the minor and major axes to generate the crack front data. In fact, it is previously shown that, for homogeneous materials it is sufficient to use just two points along the crack front for accurate prediction of the crack shape and the fatigue life [13]. The crack front profiles generated by these two different approaches are shown in Fig. 7. The circular symbols and the full lines in this figure represent the crack front profiles formed by using the fatigue crack growth analysis procedure described in Section 2, whereas the dashed lines are for the results obtained through the numerical procedure based on the invariant elliptical crack front assumption. The number of cycles corresponding to each crack front profile is provided in Table 4. As can be seen from Fig. 7, the results obtained using the two different approaches are in very good agreement, which is a validation of the proposed finite elements based crack growth analysis algorithm.
Table 4 Cycle numbers corresponding to the crack front profiles. Profile
Number of cycles (N)
1 2 3 4
178,943 307,326 409,383 494,869
4. Fatigue crack propagation in functionally graded materials The numerical results regarding fatigue crack propagation in FGMs are generated for the initially-elliptical crack depicted in Fig. 1. Two different types of gradation profiles are considered in the analyses. In the first case, material properties are assumed to vary in the direction of the minor axis as shown in Fig. 8a, whereas in the second case the variation is along the major axis as illustrated in Fig. 8b. In both examples, the total length H is taken as 15 cm, and the initial major and minor radii are set as 1.5 cm and 0.5 cm, respectively. The applied stress range is 100 MPa and the mean stress is 50 MPa. The medium is 100% metal at x = H/3 and 100% ceramic at x = 2H/3. The particular metallic and ceramic components used are Ti–6Al–4V and zirconia (ZrO2). Material properties of Ti–6Al–4V are given Table 1 and the properties of zirconia are tabulated in Table 5. The spatial variation of each of the material properties is represented by the following relation:
AðxÞ ¼ Acr exp
2 x Am ln 3 H Acr
;
Fig. 8. (a) Material property gradation in the direction of the minor axis and (b) material property gradation in the direction of the major axis.
ð10Þ
where A symbolizes any of the material parameters E, m, C, and n; and the subscripts cr and m refer to the 100% ceramic and 100% metallic material properties, respectively. Notice that A(H/3) = Am and A(2H/3) = Acr. Three different types of finite elements are used to simulate fatigue crack propagation in functionally graded materials. The volume around the crack border is modelled by utilizing collapsed 20-noded isoparametric brick elements, whereas the domain excluding this volume is discretized by means of brick and tetrahe-
Fig. 7. Crack front profiles for an elliptical crack in a homogeneous medium.
dral elements. Note that remodelling is required to compute stress intensity factors, for the crack front morphology changes under repeated loading. The number of finite elements used at a given step in the finite element computations is different from those used at the other steps. Extensively refined finite element meshes are used throughout the computations so as to assure the generation of convergent results. For example, in the initial configuration of the problem in which the direction of elastic gradation coincides with the minor axis, a total of 71,514 finite elements are used; while for FGMs with elastic gradation in the direction of the major axis total initial element number is 75,635. For each case, 960 collapsed singular elements are utilized in the modelling of the crack front region. The element number is determined through a mesh refinement study. After successive mesh refinements, the mesh resulting in convergent results is chosen to compute the stress intensity factors. Since the examined medium is functionally graded, the material properties vary from element to element. In our analyses, we followed the homogeneous finite element method, and computed the material properties of a given finite element at its centroid. The approach of utilizing homogeneous finite elements leads to the generation of highly accurate results so long as there is sufficient mesh refinement in the model. The validation study we describe in Section 3 indicates that, highly accurate fracture analyses can be conducted through the use of the homogeneous finite element method employed in the present study.
251
B. Sabuncuoglu et al. / Computational Materials Science 52 (2012) 246–252 Table 5 Properties of ZrO2 [3,14]. Modulus of elasticity (GPa)
Poisson’s ratio
Paris–Erdogan law coefficient (C)
Paris–Erdogan law exponent (n)
116.4
0.333
2.7 1016
19
Results illustrating the crack front profiles in FGMs subjected to cycling loading are presented in Figs. 9 and 10. The numbers of cycles corresponding to the crack front profiles are provided in Tables 6 and 7. From Figs. 9 and 10, it can be observed that the elliptical shape of the crack is distorted in functionally graded materials due to the existence of the material property gradation. In both cases, the crack grows more rapidly towards the 100% ceramic surface. This is the expected result since the exponent of the Paris–Erdogan law for ceramic materials is much larger compared to that for metallic materials. Due to the brittle nature of ceramics, crack growth in these materials occurs at much larger rates. Moreover, crack growth in the direction orthogonal to the direction of gradation is seen to be significant when the orthogonal direction coincides with the minor axis and suppressed when the orthogonal direction coincides with the major axis. The results given in Figs. 9 and 10 and Tables 6 and 7 indicate that the method presented in this article could prove useful in keeping the track of crack front morphology and computing the fatigue lives of functionally graded structures.
Fig. 10. Crack front profiles generated by considering property variations along the major axis.
Table 6 Numbers of cycles corresponding to the crack front profiles shown in Fig. 9. Number of cycles
1 2 3 4
5359 8529 10,472 11,863
Table 7 Numbers of cycles corresponding to the crack front profiles shown in Fig. 10.
5. Concluding remarks The main objective in the present study is to develop a computational method capable of generating useful data regarding fatigue crack propagation behaviour of functionally graded materials. For this purpose, we considered the general configuration of an initially-elliptical crack in an elastic functionally graded medium. The crack and the graded medium are modelled by means of a three dimensional finite element approach. The square-root singular behaviour of the strain field around the crack front is integrated into the finite element model by employing singular finite elements. The smooth spatial variations in the material properties of FGMs are taken into account by specifying material properties at the centroid of each finite element. Mode I stress intensity factors around the crack front are computed through the application of the displacement correlation technique. Fatigue crack propagation algorithm is based on the Paris–Erdogan law [7] and can be considered as an extension of the technique proposed by Joseph and Erdogan [8]. The comparisons of SIFs given in Section 3 display the high level of accuracy achieved in the applications of the three dimensional displacement correlation technique for both homoge-
Crack profile
Crack profile
Number of cycles
1 2 3 4
5349 8512 10,442 11,768
neous and functionally graded materials. Moreover, results on fatigue crack growth provided in Section 3 illustrate that, the elliptical crack front form, that occurs during propagation in homogeneous materials, can be captured quite accurately by means of the proposed computational algorithm. In Section 4, we present results pertaining to fatigue crack growth in functionally graded materials, which demonstrate the distortion of the elliptical shape of the crack due to the influence of the material property gradation. Under cyclic loading, the crack is seen to grow much more rapidly towards the 100% ceramic plane; and crack growth increments are generally relatively small towards the 100% metallic plane. Furthermore, crack growth rates in the direction orthogonal to the direction of gradation are larger when the orthogonal direction coincides with the minor axis than when the orthogonal direction coincides with the major axis. It is concluded that the computational method proposed in this study is an effective way of generating crack front profiles and computing fatigue lives of functionally graded structures. This method is deemed to be among the first in the literature dealing with the phenomenon of fatigue in functionally graded materials. It can further be extended to study the propagation of mixed-mode cracks in FGMs. This requires the computation of the modes I–III stress intensity factors as well as the utilization of a relation that correlates crack growth rates with mixed-mode stress intensity factor ranges. References
Fig. 9. Crack front profiles generated by considering property variations along the minor axis.
[1] M.C. Walters, G.H. Paulino, R.H. Dodds Jr., International Journal of Solids and Structures 41 (2004) 1081–1118. [2] O. Inan, S. Dag, F. Erdogan, Materials Science Forum 492–493 (2005) 373–378.
252
B. Sabuncuoglu et al. / Computational Materials Science 52 (2012) 246–252
[3] B. Yildirim, S. Dag, F. Erdogan, International Journal of Fracture 132 (2005) 369–395. [4] A.O. Ayhan, International Journal of Solids and Structures 46 (2009) 796– 810. [5] H. Yu, L. Wu, L. Guo, H. Wu, S. Du, International Journal of Solids and Structures 47 (2010) 2178–2189. [6] Ch. Zhang, M. Cui, J. Wang, X.W. Gao, J. Sladek, V. Sladek, Engineering Fracture Mechanics 78 (2011) 585–604. [7] P. Paris, F. Erdogan, ASME Journal of Basic Engineering 85 (1963) 528–534. [8] P.F. Joseph, F. Erdogan, International Journal of Fracture 41 (1989) 105–131.
[9] [10] [11] [12] [13] [14]
ANSYS, ANSYS Basic Analysis Procedures Guide, Release 5.4, ANSYS Inc., Canonsburg PA, USA, 1997. G.R. Irwin, ASME Journal of Applied Mechanics 29 (1962) 651–654. R.O. Ritchie, D.L. Davidson, B.L. Boyce, J.P. Campbell, O. Roder, Fatigue and Fracture of Engineering Materials and Structures 22 (1999) 621–631. A.O. Ayhan, H.F. Nied, International Journal for Numerical Methods in Engineering 54 (2002) 899–921. J.M. Martinez-Esnaola, A. Martin-Meizoso, International Journal of Fracture 109 (2001) L17–L22. J. Acala, M. Anglada, Materials Science and Engineering A 232 (1997) 103–109.
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Three dimensional computational analysis of fatigue crack propagation in functionally graded materials Baris Sabuncuoglu a,⇑, Serkan Dag a, Bora Yildirim b a b
Department of Mechanical Engineering, Middle East Technical University, Ankara 06800, Turkey Department of Mechanical Engineering, Hacettepe University, Ankara 06800, Turkey
a r t i c l e
i n f o
Article history: Available online 30 June 2011 Keywords: Functionally graded materials Fracture mechanics Fatigue crack propagation Finite element method Stress intensity factors
a b s t r a c t This article proposes a new finite elements based three dimensional method developed to study the phenomenon of fatigue crack propagation in functionally graded materials (FGMs). The particular problem examined in detail is that of an initially-elliptical crack located in a functionally graded medium, subjected to mode I cyclic loading. The crack is modelled by employing three dimensional finite elements; and the stress intensity factors (SIFs) around the crack front are computed by the application of the displacement correlation technique (DCT). Fatigue crack propagation calculations are based on the Paris– Erdogan crack growth law. The developed procedure makes it possible to generate the crack front profiles corresponding to given numbers of loading cycles. Proposed methods are validated by examining the behaviours of both stationary and propagating cracks. Numerical analyses conducted for an initiallyelliptical crack lying in a graded medium demonstrate that the proposed model can effectively capture the evolution of the crack front morphology. Hence, the methodology presented in this article can be used along with fracture toughness data to assess the fatigue lives of functionally graded components containing cracks that possess arbitrary front profiles. Ó 2011 Elsevier B.V. All rights reserved.
1. Introduction Functionally graded materials (FGMs) are multiphase composites that possess spatial variations in the volume fractions of the constituents. These variations are intentionally introduced so as to take advantage of desirable properties of the constituents, such as the high thermal resistance of ceramics and superior fracture toughness of metals. FGMs have been used successfully in a number of technological applications including thermal barrier coatings, solid oxide fuel cells, high performance cutting tools, and biomedical materials. The spatial variations in the volume fractions of the constituents of FGMs require that certain physical properties be represented by continuous functions of spatial coordinates for analytical and computational purposes. Fracture mechanics of FGMs has been studied extensively in order to develop an understanding of the influence of the continuous gradation on the cracking behaviour. Although great bulk of this work is based on the two-dimensional fracture mechanics theory, in recent years there have been attempts to build computational models for three dimensional fracture mechanics analysis of functionally graded materials. Among the computational methods ⇑ Corresponding author. Present address: Interuniversity Microelectronics Center (IMEC), Kapeldreef 75, B-3001 Leuven, Belgium. Tel.: +44 (0) 1509 227520. E-mail addresses: [email protected] (B. Sabuncuoglu), [email protected] (S. Dag). 0927-0256/$ - see front matter Ó 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2011.06.010
proposed for three dimensional fracture analysis of FGMs, we can mention the domain integral method [1], the displacement correlation technique (DCT) [2,3], the method of enriched finite elements [4], the interaction integral approach [5], and the boundary domain integral method [6]. Note that in all these studies on three dimensional computational fracture mechanics of FGMs, the focus is on stationary cracks, and there seems to be no work on fatigue behaviour of FGMs. Present study is directed towards developing a three dimensional computational analysis technique to study the phenomenon of fatigue crack propagation in functionally graded structures. Fatigue crack growth computations are carried out by considering an embedded crack in the examined functionally graded medium. The initial front profile of the embedded crack is assumed to be in the shape of an ellipse. Elliptical cracks or voids are known to be the source of various structural failures especially under repeated loading conditions. Fatigue crack growth calculations involving elliptical cracks require the computation of the stress intensity factors (SIFs) around the crack border. We use the displacement correlation technique [2,3] in conjunction with the finite element method to evaluate the stress intensity factors. The subsequent front profiles of the initially-elliptical cracks subjected to cyclic loading are generated by adopting the fatigue crack growth analysis methods proposed by Paris and Erdogan [7], and Joseph and Erdogan [8]. The main outcome of this study is the finite elements based three dimensional computational analysis tool, which can be
B. Sabuncuoglu et al. / Computational Materials Science 52 (2012) 246–252
247
used to predict the fatigue lives of functionally graded components. Presented numerical results illustrate the influence of material property gradation on the fatigue propagation characteristics of cracks located in cyclically loaded ceramic–metal FGMs.
2. The computational procedure In this section, we provide the details of the three dimensional computational analysis methodology developed to examine fatigue crack growth in functionally graded materials. The geometry of the particular problem considered is depicted in Fig. 1. An elliptical crack resides in an elastic functionally graded medium, whose material properties are functions of the x-coordinate. The major and minor radii of the elliptical crack are denoted by c and a, respectively. The crack is loaded through the normal stress r applied over the remote bounding planes. Since, the material is assumed to possess property gradation in x-direction, one-quarter of the medium is modelled by considering symmetries about xyand xz-planes. The finite element model for the problem depicted in Fig. 1 is constructed by using the general purpose finite element analysis software ANSYS [9]. Fig. 2 shows the quarter model generated and close-up view of the region around the crack front. The dimensions were taken such that the bounding planes would have no effect on the behaviour of the elliptical crack. In order to simulate the square-root singularity of the strain field near the crack border, the volume around this border is modelled by utilizing collapsed 20-noded isoparametric brick elements, whereas the domain excluding this volume is discretized by means of brick and tetrahedral elements. Smooth spatial variations of the material properties of the graded medium are taken into account in the finite element model by computing the properties of each finite element at its centroid. The properties are uniform over each element; however they vary from element to element according to the position of the centroid. This is the so called homogeneous finite element approach and known to lead to highly accurate numerical results so long as there is appropriate mesh refinement in the finite element model [2,3]. If the displacement field for the cracked functionally graded medium is known, the mode I stress intensity factor distribution around the crack front can be evaluated by using one of the available computational methods. In the present study, we use the displacement correlation technique in the evaluation of the SIFs. In this technique, the expression of the stress intensity factor is derived by considering a point P on the crack front. Fig. 3 shows such
Fig. 2. (a) The quarter model and (b) the crack front region.
a point P lying on the border of the crack. The origin of a Cartesian coordinate system is also located at point P. This coordinate system comprises tangential (t), normal (n), and binormal (b) directions. Fig. 4 depicts the deformed shape of the crack in the normal plane passing through point P. The normal plane is formed by n and b
Fig. 1. (a) An elliptical crack in an infinite volume and (b) a point P on the elliptical crack front.
248
B. Sabuncuoglu et al. / Computational Materials Science 52 (2012) 246–252
crack front is specified as Dq0. The number of cycles corresponding to this increment is computed using
DN 0 ffi
Dq0 : C 0 DK nI00
ð5Þ
The increments at other points around the crack front are then calculated by substituting Eq. (5) into Eq. (4): n
Dqi ffi Fig. 3. Crack front and the local coordinate system at point P.
C i DK Iii Dq0 ; C 0 DK nI00
ð6Þ
where the subscript i stands for a generic point on the crack front. A sufficiently small value needs to be used for Dq0 to be able to obtain accurate results through Eqs. (5) and (6). After calculating the crack growth increments for a specified number of points on the crack front, the new crack front profile is generated by connecting the new locations by splines. The procedure is then repeated to form the crack front profile corresponding to the next specified increment.
3. Verification of the developed procedure Fig. 4. Deformed shape of the crack surface in the normal plane.
axes. The nodes 1–3 in Fig. 4 are considered to be at the edge of a collapsed 20-noded brick element. Under mode I loading conditions, the asymptotic displacements of these nodes in b-direction can be found from [3]
ub ðnÞ ¼
4ð1 ðmP Þ2 Þ EP
rffiffiffiffiffiffiffi jnj KI; 2p
ð1Þ
where KI is the mode I stress intensity factor; EP and mP are the modulus of elasticity and the Poisson’s ratio computed at point P, respectively; and the superscript P implies that the corresponding material property should be calculated at point P. The mode I SIF at P is defined as [2,3]
) ub ðjnjÞ KI ¼ lim pffiffiffiffiffiffi : jnj 4ð1 ðmP Þ2 Þ jnj!0 pffiffiffiffiffiffiffi P 2pE
(
ð2Þ
Then, using Eqs. (1) and (2) and the displacements of the nodes 2 and 3, it is possible to express the mode I stress intensity factor in the following form:
# 3=2 3=2 R3 ub2 R2 ub3 pffiffiffiffiffipffiffiffiffiffi : KI ¼ R2 R3 ðR3 R2 Þ 4ð1 ðmP Þ2 Þ pffiffiffiffiffiffiffi P 2pE
"
ð3Þ
R2 and R3 in Eq. (3) are the distances of nodes 2 and 3 with respect to the origin. The fatigue crack growth calculations for the configuration shown in Fig. 1 are based on the mathematical models proposed by Paris and Erdogan [7], and Joseph and Erdogan [8]. Suppose that the stress applied to the graded medium r varies cyclically in time with a constant stress range Dr and a constant mean stress. Then, the growth rate at a point on the crack front can be calculated using Paris–Erdogan law [7,8], which is expressed as:
dq ¼ C DK nI ; dN
ð4Þ
where q is the radius of curvature at the considered point on the crack front, N is the number of cycles, DKI is the stress intensity factor range, and C and n are material parameters. In the finite element implementation, first, the increment for a particular point on the
We validated the methods described in the previous section on problems involving stationary and propagating cracks. For the stationary cracks, the application of the displacement correlation technique is verified by considering two different problems, the first of which is the problem of an elliptical crack in a homogeneous medium, and the second is that of a surface crack in a functionally graded plate. As for the propagating cracks, the validation study is conducted by examining fatigue growth in a homogeneous medium. In what follows below, we present the results generated for these validation problems. By considering the elliptical crack problem shown in Fig. 1, Irwin [10] derived a simple closed form expression for the mode I stress intensity factor. In this derivation, the material was assumed to be homogeneous. We first focus on this problem solved by Irwin [10] and provide comparisons of the mode I SIFs computed by the DCT developed in the present study to those obtained through Irwin’s formula. These comparisons are given in Table 1. In the calculation of these results, the medium shown in Fig. 1 is assumed to be homogeneous and made of Ti–6Al–4V. The material properties are taken from Refs. [3,11] and provided in Table 2. The minor and major radii of the crack are respectively set as 0.5 cm and 1.5 cm. The normalized mode I stress intensity factor provided in Table 1 at various locations around the crack front is defined as
KI K In ¼ pffiffiffiffiffiffi : r pa
ð7Þ
Examining the results given in Table 1, it is seen that the normalized mode I stress intensity factors computed by the DCT are in excellent agreement with those obtained through Irwin’s for-
Table 1 Comparisons of the normalized stress intensity factors computed by the displacement correlation technique to those obtained by the formula given by Irwin [10]. 2Ø/p
Irwin [10]
DCT
% Difference
1 (minor axis) 0.875 0.75 0.5 0.375 0.25 0.125 0 (major axis)
0.8979 0.8902 0.8671 0.7752 0.7075 0.6293 0.5540 0.5184
0.8977 0.8898 0.8675 0.7755 0.7068 0.6270 0.5523 0.5180
0.02 0.05 0.04 0.03 0.10 0.37 0.31 0.09
249
B. Sabuncuoglu et al. / Computational Materials Science 52 (2012) 246–252 Table 2 Properties of Ti–6Al–4V [3,11]. Modulus of elasticity (GPa)
Poisson’s ratio
Paris–Erdogan law coefficient (C)
Paris–Erdogan law exponent (n)
105.8
0.298
5.2 1012
3.17
Table 3 Comparisons of the normalized mode I stress intensity factors computed by the displacement correlation technique to those given by Walters et al. [1]. a/h = 0.5, a/ c = 2, E(h)/E1 = 5.
mula. Hence, it is possible to compute the mode I stress intensity factors to within a high degree of accuracy by employing the displacement correlation technique developed in the present study. To be able to demonstrate the applicability of the developed method for functionally graded materials, we provide additional comparisons by considering the static surface crack problem studied by Walters et al. [1]. The problem geometry is shown in Fig. 5. The figure depicts a semi-elliptical surface crack located in an isotropic FGM plate, that is under the effect of remote normal stress r at the ends y = ±L. The modulus of elasticity of the plate is represented by an exponential function, and given by
EðxÞ ¼ E1 expðbxÞ;
ð8Þ
where b is a nonhomogeneity constant. Poisson’s ratio is assumed to be constant and equal to 0.25. L and B are taken as sufficiently large by Walters et al. [1]; and as a result the stress field around the crack front is not affected by the boundaries at y = ±L and those at z = ±B. Table 3 tabulates the normalized mode I stress intensity factors computed by using the displacement correlation technique developed in the present study and those given by Walters et al. [1], which are evaluated by means of a three dimensional domain integral method. The normalized mode I SIF is defined by
Fig. 5. (a) Surface cracked FGM plate problem considered by Walters et al. [1] and (b) plan view of the semi-elliptical surface crack.
2Ø/p
Walters et al. [1]
DCT
% Difference
1 (major axis) 0.875 0.75 0.5 0.375 0.25 0.125 0 (minor axis)
0.580 0.595 0.629 0.679 0.685 0.677 0.656 0.596
0.587 0.600 0.631 0.686 0.693 0.684 0.664 0.613
1.21 0.84 0.32 1.03 1.17 1.03 1.22 2.85
KI K In ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi :
r
pa 1þ1:464ðc=aÞ1:65
ð9Þ
The results obtained by the developed displacement correlation technique are in very good agreement with those provided by Walters et al. [1], which is indicative of the high level of accuracy achieved through the application of this technique. Percent difference at the point where the crack front intersects the free surface, i.e. at / = 0, is relatively larger compared to the percent differences computed at the other points. The degree of singularity at this point is different from 0.5 and depends on the Poisson’s ratio [12]. Both in the article by Walters et al. [1] and in the present study, the special singular behaviour at the free surface is not considered, and the results are obtained by assuming that there is a square-root singularity at all points around the crack front including the point of intersection at the free surface. The change in the singularity at this single point could be the main reason behind the difference between the results computed for / = 0. However, although the percent difference is relatively larger at this point, it is less than 3% and the agreement is still acceptable. The finite element mesh used in the calculation of the results tabulated in Table 3 is shown in Fig. 6. Next, we consider fatigue crack propagation in a homogeneous medium subjected to cyclic loading for further verification. The initial geometry of the crack is as depicted in Fig. 1; and the medium is assumed to be subjected to a cyclically varying stress with a constant stress range given as Dr = 100 MPa and a constant mean stress of 50 MPa. The material is homogeneous Ti–6Al–4V. The initial minor and major radii of the elliptical crack are taken as 0.5 cm and 1.5 cm, respectively. Fatigue growth behaviour of this crack is examined by using two different methods, the first of which is the method elucidated in the previous section. The second method is based on the assumption that the initially-elliptical crack retains
Fig. 6. Finite element model used in the solution of the semi-elliptical crack problem considered by Walters et al. [1].
250
B. Sabuncuoglu et al. / Computational Materials Science 52 (2012) 246–252
its elliptical form during propagation. Hence, in the second method it is possible to generate the crack front profile at each step by considering only two points on the crack border. In our calculations pertaining to this second method, we used the points on the minor and major axes to generate the crack front data. In fact, it is previously shown that, for homogeneous materials it is sufficient to use just two points along the crack front for accurate prediction of the crack shape and the fatigue life [13]. The crack front profiles generated by these two different approaches are shown in Fig. 7. The circular symbols and the full lines in this figure represent the crack front profiles formed by using the fatigue crack growth analysis procedure described in Section 2, whereas the dashed lines are for the results obtained through the numerical procedure based on the invariant elliptical crack front assumption. The number of cycles corresponding to each crack front profile is provided in Table 4. As can be seen from Fig. 7, the results obtained using the two different approaches are in very good agreement, which is a validation of the proposed finite elements based crack growth analysis algorithm.
Table 4 Cycle numbers corresponding to the crack front profiles. Profile
Number of cycles (N)
1 2 3 4
178,943 307,326 409,383 494,869
4. Fatigue crack propagation in functionally graded materials The numerical results regarding fatigue crack propagation in FGMs are generated for the initially-elliptical crack depicted in Fig. 1. Two different types of gradation profiles are considered in the analyses. In the first case, material properties are assumed to vary in the direction of the minor axis as shown in Fig. 8a, whereas in the second case the variation is along the major axis as illustrated in Fig. 8b. In both examples, the total length H is taken as 15 cm, and the initial major and minor radii are set as 1.5 cm and 0.5 cm, respectively. The applied stress range is 100 MPa and the mean stress is 50 MPa. The medium is 100% metal at x = H/3 and 100% ceramic at x = 2H/3. The particular metallic and ceramic components used are Ti–6Al–4V and zirconia (ZrO2). Material properties of Ti–6Al–4V are given Table 1 and the properties of zirconia are tabulated in Table 5. The spatial variation of each of the material properties is represented by the following relation:
AðxÞ ¼ Acr exp
2 x Am ln 3 H Acr
;
Fig. 8. (a) Material property gradation in the direction of the minor axis and (b) material property gradation in the direction of the major axis.
ð10Þ
where A symbolizes any of the material parameters E, m, C, and n; and the subscripts cr and m refer to the 100% ceramic and 100% metallic material properties, respectively. Notice that A(H/3) = Am and A(2H/3) = Acr. Three different types of finite elements are used to simulate fatigue crack propagation in functionally graded materials. The volume around the crack border is modelled by utilizing collapsed 20-noded isoparametric brick elements, whereas the domain excluding this volume is discretized by means of brick and tetrahe-
Fig. 7. Crack front profiles for an elliptical crack in a homogeneous medium.
dral elements. Note that remodelling is required to compute stress intensity factors, for the crack front morphology changes under repeated loading. The number of finite elements used at a given step in the finite element computations is different from those used at the other steps. Extensively refined finite element meshes are used throughout the computations so as to assure the generation of convergent results. For example, in the initial configuration of the problem in which the direction of elastic gradation coincides with the minor axis, a total of 71,514 finite elements are used; while for FGMs with elastic gradation in the direction of the major axis total initial element number is 75,635. For each case, 960 collapsed singular elements are utilized in the modelling of the crack front region. The element number is determined through a mesh refinement study. After successive mesh refinements, the mesh resulting in convergent results is chosen to compute the stress intensity factors. Since the examined medium is functionally graded, the material properties vary from element to element. In our analyses, we followed the homogeneous finite element method, and computed the material properties of a given finite element at its centroid. The approach of utilizing homogeneous finite elements leads to the generation of highly accurate results so long as there is sufficient mesh refinement in the model. The validation study we describe in Section 3 indicates that, highly accurate fracture analyses can be conducted through the use of the homogeneous finite element method employed in the present study.
251
B. Sabuncuoglu et al. / Computational Materials Science 52 (2012) 246–252 Table 5 Properties of ZrO2 [3,14]. Modulus of elasticity (GPa)
Poisson’s ratio
Paris–Erdogan law coefficient (C)
Paris–Erdogan law exponent (n)
116.4
0.333
2.7 1016
19
Results illustrating the crack front profiles in FGMs subjected to cycling loading are presented in Figs. 9 and 10. The numbers of cycles corresponding to the crack front profiles are provided in Tables 6 and 7. From Figs. 9 and 10, it can be observed that the elliptical shape of the crack is distorted in functionally graded materials due to the existence of the material property gradation. In both cases, the crack grows more rapidly towards the 100% ceramic surface. This is the expected result since the exponent of the Paris–Erdogan law for ceramic materials is much larger compared to that for metallic materials. Due to the brittle nature of ceramics, crack growth in these materials occurs at much larger rates. Moreover, crack growth in the direction orthogonal to the direction of gradation is seen to be significant when the orthogonal direction coincides with the minor axis and suppressed when the orthogonal direction coincides with the major axis. The results given in Figs. 9 and 10 and Tables 6 and 7 indicate that the method presented in this article could prove useful in keeping the track of crack front morphology and computing the fatigue lives of functionally graded structures.
Fig. 10. Crack front profiles generated by considering property variations along the major axis.
Table 6 Numbers of cycles corresponding to the crack front profiles shown in Fig. 9. Number of cycles
1 2 3 4
5359 8529 10,472 11,863
Table 7 Numbers of cycles corresponding to the crack front profiles shown in Fig. 10.
5. Concluding remarks The main objective in the present study is to develop a computational method capable of generating useful data regarding fatigue crack propagation behaviour of functionally graded materials. For this purpose, we considered the general configuration of an initially-elliptical crack in an elastic functionally graded medium. The crack and the graded medium are modelled by means of a three dimensional finite element approach. The square-root singular behaviour of the strain field around the crack front is integrated into the finite element model by employing singular finite elements. The smooth spatial variations in the material properties of FGMs are taken into account by specifying material properties at the centroid of each finite element. Mode I stress intensity factors around the crack front are computed through the application of the displacement correlation technique. Fatigue crack propagation algorithm is based on the Paris–Erdogan law [7] and can be considered as an extension of the technique proposed by Joseph and Erdogan [8]. The comparisons of SIFs given in Section 3 display the high level of accuracy achieved in the applications of the three dimensional displacement correlation technique for both homoge-
Crack profile
Crack profile
Number of cycles
1 2 3 4
5349 8512 10,442 11,768
neous and functionally graded materials. Moreover, results on fatigue crack growth provided in Section 3 illustrate that, the elliptical crack front form, that occurs during propagation in homogeneous materials, can be captured quite accurately by means of the proposed computational algorithm. In Section 4, we present results pertaining to fatigue crack growth in functionally graded materials, which demonstrate the distortion of the elliptical shape of the crack due to the influence of the material property gradation. Under cyclic loading, the crack is seen to grow much more rapidly towards the 100% ceramic plane; and crack growth increments are generally relatively small towards the 100% metallic plane. Furthermore, crack growth rates in the direction orthogonal to the direction of gradation are larger when the orthogonal direction coincides with the minor axis than when the orthogonal direction coincides with the major axis. It is concluded that the computational method proposed in this study is an effective way of generating crack front profiles and computing fatigue lives of functionally graded structures. This method is deemed to be among the first in the literature dealing with the phenomenon of fatigue in functionally graded materials. It can further be extended to study the propagation of mixed-mode cracks in FGMs. This requires the computation of the modes I–III stress intensity factors as well as the utilization of a relation that correlates crack growth rates with mixed-mode stress intensity factor ranges. References
Fig. 9. Crack front profiles generated by considering property variations along the minor axis.
[1] M.C. Walters, G.H. Paulino, R.H. Dodds Jr., International Journal of Solids and Structures 41 (2004) 1081–1118. [2] O. Inan, S. Dag, F. Erdogan, Materials Science Forum 492–493 (2005) 373–378.
252
B. Sabuncuoglu et al. / Computational Materials Science 52 (2012) 246–252
[3] B. Yildirim, S. Dag, F. Erdogan, International Journal of Fracture 132 (2005) 369–395. [4] A.O. Ayhan, International Journal of Solids and Structures 46 (2009) 796– 810. [5] H. Yu, L. Wu, L. Guo, H. Wu, S. Du, International Journal of Solids and Structures 47 (2010) 2178–2189. [6] Ch. Zhang, M. Cui, J. Wang, X.W. Gao, J. Sladek, V. Sladek, Engineering Fracture Mechanics 78 (2011) 585–604. [7] P. Paris, F. Erdogan, ASME Journal of Basic Engineering 85 (1963) 528–534. [8] P.F. Joseph, F. Erdogan, International Journal of Fracture 41 (1989) 105–131.
[9] [10] [11] [12] [13] [14]
ANSYS, ANSYS Basic Analysis Procedures Guide, Release 5.4, ANSYS Inc., Canonsburg PA, USA, 1997. G.R. Irwin, ASME Journal of Applied Mechanics 29 (1962) 651–654. R.O. Ritchie, D.L. Davidson, B.L. Boyce, J.P. Campbell, O. Roder, Fatigue and Fracture of Engineering Materials and Structures 22 (1999) 621–631. A.O. Ayhan, H.F. Nied, International Journal for Numerical Methods in Engineering 54 (2002) 899–921. J.M. Martinez-Esnaola, A. Martin-Meizoso, International Journal of Fracture 109 (2001) L17–L22. J. Acala, M. Anglada, Materials Science and Engineering A 232 (1997) 103–109.