- Email: [email protected]
Determination of stress intensity factors of 3D curved non-planar cracks in FGMs subjected to thermal loading
Accepted Manuscript Determination of stress intensity factors of 3D curved non-planar cracks in FGMs subjected to thermal loading Ali Shaghaghi Moghaddam, Marco Alfano PII: DOI: Reference:
S0013-7944(15)00416-6 http://dx.doi.org/10.1016/j.engfracmech.2015.07.040 EFM 4786
To appear in:
Engineering Fracture Mechanics
Received Date: Revised Date: Accepted Date:
6 February 2015 30 June 2015 3 July 2015
Please cite this article as: Moghaddam, A.S., Alfano, M., Determination of stress intensity factors of 3D curved non-planar cracks in FGMs subjected to thermal loading, Engineering Fracture Mechanics (2015), doi: http:// dx.doi.org/10.1016/j.engfracmech.2015.07.040
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
2
Determination of stress intensity factors of 3D curved non-planar cracks in FGMs subjected to thermal loading
3
Ali Shaghaghi Moghaddama , Marco Alfanob,∗
1
4 5 6 7
8
a
Young Researchers and Elite Club, Takestan Branch, Islamic Azad University, Takestan, Iran b Department of Mechanical, Energy and Management Engineering, University of Calabria, P. Bucci 44C, 87036 Rende (CS), Italy
Abstract A finite element analysis is presented about 3D curved non-planar cracks in functionally graded materials (FGMs) subjected to steady state temperature gradients. The interaction energy integral is employed in conjunction with a computational strategy that does not require a priori information about crack front and crack surface curvatures. The influence of graded thermal and mechanical properties on mixed mode stress intensity factors (SIFs) is explored in detail. Results highlighted the robustness of the proposed computational framework, which therefore can be effectively deployed in the analysis of thermal fracture in FGMs containing curved non-planar cracks.
10
Keywords: FGMs, non-planar cracks, thermal SIFs, interaction integral, micromechanics models
11
1. Introduction
9
12 13 14 15 16 17 18 19 20 21 22 23
Functionally Graded Materials (FGMs) are nonhomogenous materials whose microstructure and composition are engineered so that to achieve a continuos spatial variation of mechanical and physical properties. A significant body of literature has explored multiple aspects of FGMs, including manufacturing and analysis, as recently reviewed in [1]. Several fabrication techniques have been developed (e.g. spark plasma sintering, plasma spray technique and physical vapor deposition) and analytical and numerical efforts have been carried out in order to sustain experimental investigations and the design of FGMs systems [2, 3, 4, 5, 6, 7, 8]. One of the most important design aspect that has been emphasized in earlier works relates to the integrity of cracked components made up of FGMs subjected to thermomechanical loading. The use of graded materials for thermal protection sys∗ Corresponding author. Tel.: +39 0984 494156 Preprint submitted to Engineering Fracture Mechanics Email address: [email protected] (Marco Alfano)
July 25, 2015
24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
tems, such as first-wall composite materials within nuclear fusion reactors, is very attractive and for this reason the analysis of cracks in FGMs under thermo-mechanical loading has been vigorously pursued in recent times. KC and Kim [9] employed the interaction energy integral to study 2D cracks in FGMs subjected to thermo-mechanical loading. The 3D finite domain integral expression and the two-dimensional and axisymmetric specializations for the energy release rate are given in [10]. Walters et al. [11] studied 3D planar cracks in FGMs under thermo-mechanical loading under mode I using the domain form of the J-integral. Yildirim et al. [12] focused on 3D planar cracks in FGMs subjected to mode I mechanical and transient thermal loading and employed the displacement correlation technique (DCT). Later on, Yildirim [13] studied 2D crack in FGMs under mode I steady-state and transient thermal loading using domain integrals. Ayhan [14, 15] proposed a different approach, based on the use of special enriched finite elements, to compute the SIFs of 3D planar cracks in FGMs under mode I and mixed mode loading. Previous works mostly addressed planar cracks, while for 3D nonplanar cracks, such as conical cracks arising from impact or indentation loading [16, 17], only few efforts can be found in the published literature. For instance Gosz and Moran [18] employed the interaction integral approach for 3D non-planar cracks, but their analysis addressed homogenous materials; similarly, Chang and Wu [19] computed mixed mode SIFs for 3D non-planar cracks in homogenous materials by using the concept of the Jk and GIII integrals. Shaghaghi et al. and Ghajar et al. [20, 21] extended the work carried out in [18] to FGMs. However, to the best of the authors knowledge none of the previous works addressed 3D non-planar cracks under thermal loading. The objective of the present study is to complement existing works by analyzing the mixed mode SIFs of 3D curved nonplanar crack in FGMs under thermal loading. An un-coupled thermo-mechanical finite element analysis is solved where the temperature field is firstly retrieved from a thermal analysis and then is employed as initial condition in a subsequent general static analysis step. In both cases gradation of material properties is accounted for1 . The determination of the effective material properties is herein carried out by resorting to micromechanics models. The self-consistent model is used to determine Young’s modulus and Poisson’s ratio, while the rule of 1
Contrary to fracture in homogeneous materials, in FGMs the SIFs and the J-integral not only depend on external loading and crack geometry, but may be heavily affected by the spatial variation of material properties.
2
87
mixture is employed to quantify the coefficient of thermal expansion and thermal conductivity. Upon completion of the finite element analysis, a post processing procedure enables the evaluation of the J-integral and interaction energy integrals. However, using the interaction energy integral, auxiliary fields need to be determined. Gosz and Moran [18] employed an orthogonal curvilinear coordinate system to set-up auxiliary fields for a 3D non-planar crack in homogenous material. A similar procedure was adopted by Shaghaghi et al. [20] to analyze 3D non planar cracks in FGMs. In both works, analytical descriptions of crack front and crack surface curvatures were required to compute the interaction integral and, in turn, the mixed mode SIFs. The approach employed herein does not require any a priori information about crack front and crack surface curvatures to locate integration stations within the finite element framework. Therefore it allows to analyze arbitrary non-planar cracks with relative ease. This simplified strategy, originally developed by the authors in [21], is herein extended to account for the effect of thermal loading on mixed mode SIFs of 3D nonplanar cracks on FGMs. To demonstrate the capabilities of the proposed computational scheme, the mixed mode SIFs for two selected boundary value problems are determined, i.e. a conical crack in an infinite half space and a lens shaped crack in a cylindrical domain. In order to assess the effectiveness of the computational scheme, results are compared with those obtained using the classical displacement correlation technique (DCT) [12]. The paper is organized as follows. In section 2 the computational procedure employed to extract the SIFs is firstly outlined and the determination of effective material properties of FGMs using micromechanics models is described. In section 3 details about the finite element implementation are provided. In turn, section 4 illustrates the results of two selected boundary value problems involving nonplanar crack configurations, i.e. lens shaped and conical cracks subjected to thermal loading.
88
2. Computational procedures
59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86
89 90 91 92 93 94 95 96
The mixed mode stress intensity factors of 3D non-planar cracks are herein obtained by solving an un-coupled thermo-mechanical finite element analysis. Specifically, the heat transfer analysis is performed considering gradation of thermal conductivity. The temperature field is then employed as initial condition in a general static analysis where gradation of mechanical properties is accounted for. Upon completion of the stress analysis, a post processing stage enables the determination of the J-integral, interaction energy integral and, in turn, the mixed mode SIFs. 3
120
In particular, the 3-D domain form of the J-integral found in [10] is employed to extract energy released (JA ) by the body per unit of crack advance over a finite segment of the crack front. The interaction energy integral is then employed to determine the mixed mode SIFs. Upon evaluation of mixed mode SIFs, the energy release rate is obtained (JK ) and a crosscheck of the results is made by computing the ratio of JK /JA . Since analytical solutions are not available for the problem geometries and boundary conditions analyzed herein, additional comparisons and judgements on the accuracy of the results is also carried out by resorting to the classical displacement correlation technique (DCT) [12]. The DCT has a very simple formulation since it essentially correlates numerical results for displacement at crack front location to evaluate SIFs. The size of the singular elements used around the crack front is an important factor to be considered to obtain sufficiently accurate results in the finite element analysis. Indeed, the accuracy is acceptable for most engineering applications provided quarter-point elements are used in the first layer of elements around the crack tip, and the mesh in the near-tip region is refined enough. On the other hand the domain integral formulation is extremely versatile and relatively simple to implement numerically. It removes the need to precisely capture he details of the singular fields near the crack tip, as a consequence there is no need to use singular elements and/or a very refined mesh. As it will be shown later, the details of the finite element mesh have been properly selected to accommodate the requirements of both methods, to provide consistent results and to ensure an objective comparative analysis.
121
2.1. The interaction energy integral for FGMs
97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119
122 123 124 125 126 127 128 129 130 131 132
Mixed mode stress intensity factors are calculated independently by means of the interaction energy integral. This is done by superimposing suitable auxiliary fields to the actual fields. A finite element procedure which enables the determination of SIFs using the interaction energy integral was originally proposed in [20, 21], and it is herein extended to include the effect of thermal loading on the mixed mode SIFs of 3D non-planar cracks in FGMs. Let consider a schematic depiction of a 3D nonplanar crack within a FGM as shown in Fig.1. Considering a finite domain V surrounding the crack front at location s, the J-integral for graded materials under general thermo-mechanical loading can be written in domain form as follows: Z 1 J(s) = − (Hlj ql,j + Hlj,j ql ) dV , (1) γ V 4
133 134 135 136 137 138
139 140 141 142
where ql is a sufficiently smooth vector field in V. In order to account for thermal loading the Eshelby’s energy-momentum tensor (H) needs to be properly modified. In particular, the mechanical strain, εmech = εtotal mn mn − α (X) T (X), the gradient of the thermal expansion coefficient, α (X),l , as well as the temperature gradient, T (X),l , need to be included in H. It follows that the tensor Hlj and its gradient, Hlj,j , are given by: Hlj
= 12 σmn εmech mn δlj − σij ui,l ,
(2)
Hlj,j
= Cmnpq,l εmech εmech pq mn − σpp {α,l T + αT,l }.
(3)
Notice that Cmnpq (X) is the spatially varying constitutive tensor for graded materials; α = α(X) is the coefficient of thermal expansion; T = T (X) is the temperature within the body; X is the spatial coordinate vector, while γ is given by the following equation: Z γ=
Lc 2
∆a(ξ)dξ,
(4)
− L2c
145
where ∆a(ξ) is the magnitude of the crack advance at crack front segment Lc , as shown in Fig.1. On the other hand the interaction energy integral is given by: Z 1 (1,2) (1,2) (1,2) J (s) = − (5) Hlj ql,j + Hlj,j ql dV , γ V
146
where
143 144
(1,2)
(1) (2)
(1,2)
Hlj,j 147
149 150 151 152 153
(2) (1)
(2) mech(1)
= Cmnpq,l εpq εmn
(2)
(6) (2)
(1)
− σpp {α,l T + α T,l } − σmn,n um,l .
(7)
Computation of SIFs can be effectively carried out by the following equation: 2 (1 + ν ) 2 1 − νs2 (1) (2) (1) (2) (1) (2) s (1,2) J (s) = KI KI + KII KII + KIII KIII , (8) Es Es (1)
148
(1) (2)
= σmn εmn δlj − {σij ui,l + σij ui,l },
Hlj
(2)
where Ki and Ki (i=I,II,III) are the SIFs for actual and auxiliary fields, respectively; Es and νs represent Young’s modulus and Poisson’s (1) ratio at crack front location s. In turn, SIFs for actual state, e.g. KI , can be obtained by setting the counterpart SIF for auxiliary state equal to (2) (2) (2) KI = 1, KII = 0 and KIII = 0. It follows that Eq. (8) yields the required SIF: Es (1) KI = J (1,2) (s) . (9) 2 (1 − νs2 ) 5
154 155 156 157
2.2. Auxiliary fields The Williams solution for the 2D asymptotic crack tip displacement fields [22], is used to set-up the auxiliary displacement and strain fields in the crack tip local coordinate system: (2) u ¯j
158
(2) r
K = I 2µ
(2) r
K r I gj (θ, ν)+ II 2π 2µ
and (2) ε¯ij
159 160 161 162
163 164 165 166
(2)
K r II gj (θ, ν)+ III 2π 2µ
1 = 2
r
r III g (θ, ν) , (10) 2π j
(2) !
(2) ∂u ¯j ∂u ¯i + ∂xj ∂xi
,
(11)
where gj (θ, ν) is the standard angular function. In the case of FGMs the effect of material gradation needs to be taken into account. In this work the auxiliary stress field is related to the auxiliary strain field through the constitutive tensor as: (2) (2) σ ¯ij = Cijkl (X) ε¯kl . (12) To be consistent with actual fields, the auxiliary fields are transformed to cartesian coordinate system. Considering the transformation tensor Q which maps the local basis into the cartesian one, the auxiliary displacement, strain and stress fields can be written as: ¯, u = QT u
(13)
T
¯ Q, ε=Q ε ¯ Q. σ = QT σ 167 168 169 170 171 172 173 174 175 176 177 178 179
Since the auxiliary fields are now referred to a fixed local coordinate system, it is observed that strain-displacement compatibility conditions is satisfied and auxiliary stresses are equilibrated. Note that due to the spatially varying constitutive tensor required to represent graded materials, the auxiliary stress field is not in equilibrium and this has to be considered. Gosz and Moran [18] presented an integral method for the analysis of mixed-mode SIFs along 3D non-planar crack fronts in homogeneous materials. They highlighted that for a proper definition of the auxiliary fields a curvilinear coordinate system should be employed, and crack surface and crack front curvatures need to be accounted for. In particular, the gradient of auxiliary displacement field requires a knowledge of the divergence of the auxiliary stress field as well as the location of integration points with respect to crack front. For this reason, analytical equations describing crack front
6
192
and crack surface were deployed. It should be noted the curvilinear coordinates obtained using analytical equations can yield incorrect locations of integration points, thereby leading to auxiliary fields inconsistent with actual fields generated by the finite element mesh. For this reason, Ghajar et al. [21] recently extended the work by Gosz and Moran [18] to FGMs and proposed an improved computational approach based on the use of a local coordinate system located at the crack front. The proposed method does not require a priori knowledge of the geometry of the crack in order to locate integration stations with respect to the crack front. Its accuracy and relatively simple implementation were demonstrated through the analysis of 3D curved non planar cracks in FGMs. This approach is herein employed in order to include the capabilities to analyze non-planar cracks in FGMs under thermo-mechanical loading.
193
2.3. Determination of material properties using micromechanics models
180 181 182 183 184 185 186 187 188 189 190 191
194 195 196 197 198 199 200
201 202 203 204 205 206
The determination of the effective material properties is carried out by resorting to micromechanics models, whose inputs are the volume fractions of the constituent materials. Micromechanics models, such as the MoriTanaka [23], the self-consistent [24] and rule of mixture, allow to determine the effective material properties by considering the volume fraction of constituents. In this work the rule of mixture is employed to obtain the coefficient of thermal expansion, α, and thermal conductivity, k, as follows: α = V1 α1 + V2 α2 ,
(14)
k = V1 k1 + V2 k2 ,
(15)
where Vi = Vi (X), ki = ki (X) and αi = αi (X) are volume fraction, thermal conductivity and thermal expansion coefficient of the material constituent i, which are all functions of the spatial coordinates. For estimation of Young’s modulus and Poisson’s ratio, the self-consistent model is employed. Specifically, the shear and bulk moduli (µ, κ) are given as a function of the constituents [24]: 1 V1 2V1 = + , κ + 4µ/3 κ1 + 4µ/3 κ2 + 4µ/3 V1 κ1 V 2 κ2 V 1 µ2 V2 µ 1 + +5 + + 2 = 0. κ1 + 4µ/3 κ2 + 4µ/3 µ − µ2 µ − µ1
207 208
(16)
The above equations are firstly solved at a specific spatial location for µ and κ, then Young’s modulus and Poisson’s ratio of the FGMs are obtained 7
209
using the following equations: E=
210 211 212 213
9µκ 3κ − 2µ , ν= . µ + 3κ 2 (µ + 3κ)
(17)
For graded materials the gradient of material properties and derivatives of material constitutive tensor, Cmnpq (X),l , are incorporated in the J-integral and interaction integral. The constitutive tensor for isotropic elastic FGMs is given by: Cmnpq (X) = λ (X) δmn δpq + µ (X) (δmp δnq + δmq δnp ) ,
214 215
where δmn is the Kronecker delta, and µ (X) and λ (X)are spatially varying Lame0 constants defined as follows: λ (X) =
216 217 218 219
E (X) ν (X) , (1 + ν (X)) (1 − 2ν (X))
221 222 223 224 225 226 227 228 229 230 231 232 233 234
µ (X) =
E (X) , 2 (1 + ν (X))
(19)
E(X) and ν(X) are modulus of elasticity and Poisson’s ratio which are determined using the self-consistent model. The derivative of constitutive tensor in Cartesian coordinate system is given in [21] and is recalled here for the sake of completeness: Cmnpq (X),l =
220
(18)
∂Cmnpq (X) ∂E (X) ∂Cmnpq (X) ∂ν (X) + (l = 1, 2, 3). (20) ∂E (X) ∂Xl ∂ν (X) ∂Xl
3. Numerical aspects 3.1. Finite element discretization of J-integral and interaction energy integral The actual fields are extracted from finite element analysis carried out using ABAQUS/Standard. A heat transfer analysis is firstly performed to obtain temperature distribution in the body. This step has been performed using the subroutine UMATHT available in ABAQUS/Standard which allows to account for the gradation of k. In turn, output model temperatures are set as initial conditions in the subsequent mechanical analysis, where gradation of mechanical properties is incorporated through subroutine UMAT. Upon completion of the thermo-mechanical analysis a post processing step is carried out to extract mixed mode SIFs. In the finite element context the domain V embeds a few elements surrounding the crack tip at crack front location s, which is in turn represented by a node in the finite element mesh (see Fig. 2). Having defined the finite 8
235 236
domain V the discretized versions of J-integral and interaction integrals are obtained as: elems gps 1 X X J(s) = [Hlj ql,j ]p + [Hlj,j ql ]p detJ wp , (21) γ p J
237 238 239 240 241 242 243 244 245
(1,2)
1 (s) = γ
V gps elems X X h V
(1,2) Hlj ql,j
p
i p
+
h
(1,2) Hlj,j ql
i p
detJ wp ,
where summation is done over the integration points (gps) for each element included in the domain V . Parameters wp and detJ denote the weight factors and the Jacobian determinant, respectively. The effect of element type and number of Gaussian points on the accuracy of finite element results have been studied in [25]. It was shown that 20-noded brick elements with 2×2×2 Gauss quadrature rule are well suited for modeling 3D FGMs. Accordingly, a similar approach was undertaken in this work. Material properties and their derivatives at integration points are computed using an isoparametric approach, such that for a given generic field variable (Φ): Φ (X) =
20 X
Φ i hi ,
(22)
∂hi , ∂Xi
(23)
i=1
∂Φ (X) = ∂Xi 246 247
248 249 250 251 252 253 254
20 X i=1
Φi
where Φ (X) = {E (X) , ν (X) , α (X) , k (X) , T (X)}, Φi is the nodal values and hi is the shape function at node i. 3.2. Evaluating the q function In order to evaluate the q function the finite domain V drawn in Fig. ¯ on crack front segment Lc is ∆a(ξ) in the 2 is considered. The values of q ¯1 while it is zero on the boundary of the considered domain. As direction e ¯ is set to zero for the nodes indicated by ◦ shown in Fig. 2, the value of q and unity for those indicated by •. It follows that for an element inside the ¯ function at node i is given by: domain V , the q 0 for nodes on boundary of domain V i q¯l = (24) ¯1 e for inner nodes (i)
255 256
¯ at node i in local coordinate system. where q¯l are the components of q ¯ is then mapped to the cartesian coordinate system as follows: The vector q ¯, q = QT q 9
(25)
257 258
and consistent with isoparametric finite element formulation, ql and ql,j at integration points are interpolated using standard shape functions: ql =
20 X
(i)
hi ql ,
(26)
i=1
ql,j =
20 X
(i)
hi,j ql .
i=1
259
260 261 262 263 264 265 266 267 268 269 270
4. Numerical examples 4.1. Boundary value problems selected for procedure verification In this section the developed numerical procedures are applied to the analysis of selected boundary value problems involving a conical crack and a lens shaped crack in three-dimensional bodies made up of FGMs. A conical crack in an infinite half space is shown in Fig. 3(a). The gradation of mechanical and thermal properties occurs in the Y -direction and extends from zero to Y2 . The material properties at top surface of the domain (Y =0) are E1 = 5, ν1 = 0.2, α1 = 1, k1 = 1, while at the bottom of the graded region (Y = Y2 ) E2 = 1, ν2 = 0.4, α2 = 5, k2 = 5. Then the volume fraction of constituent materials is considered to vary in the Y -direction according to a power function with exponent n as follows: ( n |y| if y < y2 , y2 V2 = (27) 1 otherwise, V1 = 1 − V2 .
271 272 273 274 275 276 277 278 279 280 281 282 283
The volume fraction V2 , Young’s modulus, Poisson’s ratio, coefficient of thermal expansion and thermal conductivity of the FGMs are drawn in Fig. 4 for varying values of the power exponent n. It is observed that for increasing values of the exponent n, the volume fraction of the stiffer material (i.e. material #1) is mostly dominant in the graded region. For increasing values of n the softer material (i.e. material #2) is introduced toward the end of the graded region. Notice that outside the graded region, that is Y > Y2 , the material is assumed to be homogenous and isotropic with properties equal to those of material #2. Therefore for higher n the graded region exhibits mostly the behavior of material #1 and the rate of change of V2 across the boundary, Y = Y2 , is quite sharp. On the other hand for lower n, for instance n = 0.2, the volume fraction of material #2 quickly increases and the transition across the boundary is rather smooth. It is 10
284 285 286 287 288 289 290 291 292 293 294
295 296 297 298 299 300 301 302
303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318
then concluded that lower values of the index n lead to a smooth variation of material properties. As it will be shown in the next section, the conical crack is studied considering both uniform and non-uniform temperature loading and free boundary conditions (free grips). Fig. 3(b) illustrates a cylinder of radius r and height H with an embedded lens-shaped crack whose surface geometry is characterized by radius R and angle β. In particular the geometric parameters of the model are given as follows: 2r /H=1, 2r /R=10, R=1 and β = π/4. The material properties at r=0 are considered to be E1 = 1, ν1 = 0.4, α1 = 5, k1 = 5, while at the outer surface of the cylinder are E2 = 5, ν2 = 0.2, α2 = 1, k2 = 1. The volume fraction varies as a power function with exponent n: r n V2 = , (28) 5 V1 = 1 − V2 . (29) In Fig. 5 the volume fraction V1 , Young’s modulus (E) and coefficient of thermal expansion (α) are also shown. In addition, the product (E · α) is displayed. This parameter is considered because the thermal strain in the cylinder, εth is proportional to α and, in turn, the associated normal stress scales with E · εth . As a consequence it provides an approximate measure of the normal stress acting on the cylinder. Finally, notice that boundary conditions include uniform temperature variation while both fixed and free grips conditions will be analyzed. 4.2. Finite element models The finite element models of the conical and lens shaped cracks are given in Fig. 6. In particular, a conical crack is embedded in a cubic domain whose dimensions are 600 × 600 × 300 units length and which is shown in Fig. 6(a). In order to decrease the computational costs, only one quarter of the domain has been included in the model. The crack circumscribes a circle with a radius r=10 units length on the free surface (cfr. Fig. 3(a)), and intersects this last at an angle of 45◦ . In addition it extends 15 units length into the solid domain. The mesh shown in the figure comprises around 9500 brick elements (C3D20) and 42000 nodes. Notice that a sensitivity analysis to the number of elements along the crack front was already performed in [21]. The results obtained with finite element meshes comprising six, eight and ten elements along the 90-deg sector were compared with an analytical solution. The results indicated that fairly accurate results could be obtained using only six elements along crack front. In the present study an even higher number of elements was 11
330
deployed, i.e. 18 elements are inserted along crack front and 16 elements are surrounding the crack tip as shown in Fig. 6(b). The size of the elements incident at crack front, Le , is such that Le /a = 1/50. On the other hand, the finite element discretization of the lens shaped crack is given in Fig. 6(c). Once again, symmetry conditions were exploited and only one quarter of the cylinder was considered for modeling purposes. The mesh comprised around 2500 brick elements (C3D20) and the size of the smallest element incident at crack front was equal to Le /R = 0.01. Notice that the finite element discretization around the crack front is similar to that employed for the conical crack. However, the work by Gosz and Moran indicated that for a lens shaped crack [18] ten elements along a 90-deg sector may suffice, a similar choice was also made herein.
331
4.3. Results and discussion
319 320 321 322 323 324 325 326 327 328 329
332 333 334 335 336 337 338 339 340 341 342
343 344 345 346 347 348 349 350 351 352 353
Since reference solutions are not available for the boundary value problems analyzed in the present work, the results were compared with those obtained through an objective and consistent use of the DCT. From this standpoint, a proper selection of the mesh details around the crack tip region is very important, since the size of singular elements around the crack front is an important factor in order to obtain accurate results. A in-depth discussion on this point was already provided in previous related works [3, 26, 27], and it was shown that element size ranging from 1/20a to 1/32a, where a is the crack length, can provide accurate results. Similarly, preliminary convergence analyses carried out in the present work indicated that element size equal to 1/50a was able to provide satisfactory results. 4.3.1. Procedure verification considering a conical crack under uniform temperature change The results pertaining to the conical crack are illustrated first. Notice that the body is initially assumed to be at a uniform temperature T0 = 20 and subsequently it is heated up to a uniform temperature T = 100 (cfr. Fig. 3(a)). For the sake of verification, the normalized2 mixed mode SIFs obtained using the domain integral for different domain sizes, rd /Le , and the DCT approach were compared as shown in Fig. 7. Computations were made assuming n = 1. The deviations among the two sets of results is within 2% thereby testifying the robustness of the developed computational procedures. An additional cross-check was made by analyzing the ratio 2
Mixed mode SIFs were normalized against K0 =
12
√ E1 α T πa. 1−ν1 1 0
354 355 356 357 358 359 360 361 362 363 364 365 366
367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392
between the energy release rate obtained from stress intensity factors, Jk , and that stemming from the domain integral, JA , which is also given in Fig. 7. Again, the deviation from unity is less than 1% which implies that the results are consistent and domain independent. The sensitivity of the SIFs to the power exponent n was also analyzed and the results are shown in Fig. 8. The normalized SIFs, KIN and KIIN , display an increase as the power exponent n increases. This effect is also explained by the above discussion concerning the trend of V2 with respect to power exponent n. In particular, the induced material gradation, see Fig. 4, enables an increase in the crack face opening displacements as the power exponent n increases. However, both SIFs consistently increases and, as a consequence, the mode mixity is almost unchanged in the investigated range of n. 4.3.2. Steady state temperature gradients To study the influence of thermal conductivity on temperature profile in the infinite half-space containing a conical crack, it is assumed that the body is exposed to T = 100 on the upper surface, Y = 0, and to T = 20 on the lower surface, Y = 300. In this circumstance, a heat transfer analysis is first run and temperature distribution is extracted. Then temperature distribution is imported to a static analysis considering initial temperature T0 = 20 and mixed mode SIFs are in turn obtained. The temperature variation which is a function of Y is shown in Fig. 9. It is observed that as n increases, the temperature profile deviates from linearity, especially around the transition boundary located at Y2 . The temperature distribution is then imported to the general static analysis and mixed mode SIFs are calculated. The results for KIN and KIIN are shown in Fig. 10. The trend is quite similar to that of uniform expansion as Fig. 8. As the power exponent n increases, KIN and KIIN increase. It is important to highlight that negative values of KIN were also observed for n ≤0.3, which are induced by the occurrence of negative opening displacement, i.e. crack closure. In general terms, crack closure can be associated to a variety of mechanisms, for instance it can be induced by the applied loading whenever compression stresses near a crack are generated. On the other hand, crack closure might also be geometry induced because of the specific shape of the crack, such as kinked cracks. It might be also generated by material inhomogeneity. In particular, it is speculated that in the present work the observed crack closure is caused by the combined effects of thermal loading, crack geometry and the non-homogeneity brought by FGMs. Indeed, the material gradation profile strongly affects the stress and strain 13
413
fields around a crack. This point was not further analyzed herein because in doing so the effect of crack surface contact and friction should be taken into account. In particular, the treatment of crack closure would require to heavily modify the proposed formulation. For instance, the solution of frictional contact problems is strongly history dependent and tailored incremental solution procedures would be needed. Moreover, from the analysis of Fig. 10 it is possible to infer that crack opening displacements are almost always positive for the considered model material system and boundary conditions. For the lens shaped crack the initial temperature of the body is assumed to be T0 = 100 and then it is uniformly decreased to T = 20. At first it has been considered that the body is gripped at the upper and lower surfaces, and the corresponding normalized mixed mode SIFs3 are shown in Fig. 11. It is inferred that KI and KII follow similar trends, i.e. they consistently increase until n = 0.5 is reached and then decrease. A possibile explanation for such a trend can be given if one consider that near the crack front the parameter (E · α) increases as the power exponent decreases, see Fig. 5. Free boundary conditions have been also explored, and the body could freely contract in this case. The mixed mode SIFs are given in Fig. 12 and quite similar trend to that of fixed grip was observed. Although the trend with the power exponent is still the same, the SIFs are now smaller than that obtained on fixed grip boundary conditions.
414
5. Conclusion
393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412
415 416 417 418 419 420 421 422 423 424 425 426 427
In the present work 3D non-planar cracks in FGMs under thermal load were analyzed using the domain form of the J-integral and interaction energy integrals. Gradation of both thermal and mechanical properties was considered simultaneously. Micromechanics models were deployed in order to determine the effective properties of the FGM, and the volume fraction of constituent materials was assumed to display a power function form. The overall computational approach consisted of a few interrelated steps. At first, an heat transfer analysis was performed considering gradation of thermal conductivity. The temperature field was then employed as initial condition in a general static analysis where gradation of mechanical properties was accounted for. Upon completion of the stress analysis, a post processing stage enabled the determination of the J-integral, interaction energy integral and, in turn, the mixed mode SIFs. Two selected boundary value problems 3
Normalization was made against K0 =
√ 2 E1 α T πa. π 1−ν1 1 0
14
445
are analyzed, i.e. a conical crack in an infinite half space and a lens shaped crack in a cylindrical domain. In order to assess the effectiveness of the computational scheme, a comparative analysis was carried out by using the classical displacement correlation technique (DCT). The results indicated that the proposed procedure is domain independent and that the mixed mode SIFs are in a good agreement with those obtained using the DCT. However, the domain integral formulation proposed herein is extremely versatile and relatively simple to implement numerically. Moreover, if compared to the DCT, it removes the need to precisely capture the details of the singular fields near the crack tip. Therefore, singular elements and/or very refined meshes, which are computationally expensive when conducting three-dimensional analyses, are not needed. In addition, a parametric study was carried out to reveal the effect of material gradation profile on the SIFs. It was concluded that the mixed mode SIFs of curved non planar cracks are considerably influenced by material gradation. The obtained results testify that the proposed procedure is an accurate and useful tool for the analysis of 3D non planar cracks in FGMs which can assist the design of FGMs against fracture.
446
References
428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444
447 448 449
450 451 452
453 454 455
456 457 458
459 460 461
[1] Birman, V., Byrd, L.W., 2007. Modeling and analysis of functionally graded materials and structures. Applied Mechanics Reviews, 60(1-6), 195-216. [2] Dolbow, J.E., Gosz, M., 2002. On the computation of mixed mode stress intensity factors in functionally graded materials. International Journal of Solids and Structures, 39(9), 2557-2574. [3] Kim, J.-H, Paulino, G.H., 2002. Finite element evaluation of mixed mode stress intensity factors in functionally graded materials. International Journal for Numerical Methods in Engineering, 53(8), 1903-1935. [4] Ghajar, R., Shaghaghi, A., 2011. Numerical investigation of the mode III stress intensity factors in FGMs considering the effect of graded Poisson’s ratio. Engineering Fracture Mechanics, 78, 1478-1486. [5] Shaghaghi A, Ghajar R, Alfano M, 2013. Determining the mixed mode stress intensity factors of surface cracks in functionally graded hollow cylinders. Materials and Design, 43, 475-484
15
462 463 464 465
466 467 468 469
470 471 472 473
474 475 476
477 478 479
480 481 482 483
484 485 486
487 488 489
490 491 492
493 494 495
[6] Nikishkov, G.P., Atluri, S.N., 1987. Calculation of fracture mechanics parameters for an arbitrary three-dimensional crack, by the equivalent domain integral method. International Journal for Numerical Methods in Engineering, 24(9), 1801-1821. [7] Yu, H., Wu, L., Guo, L., Wu, H., Du, S., 2010. An interaction integral method for 3D curved cracks in nonhomogeneous materials with complex interfaces. International Journal of Solids and Structures, 47(16), 2178-2189 [8] Walters, M.C., Paulino, G.H., Dodds JR., R.H., 2006. Computation of mixed-mode stress intensity factors for cracks in three-dimensional functionally graded solids. Journal of Engineering Mechanics, 132(1), 1-15. [9] KC A, Kim J-H, 2008. Interaction integrals for thermal fracture of functionally graded materials. Engineering Fracture Mechanics 75, 25422565. [10] Shih, C.F., Moran, B., Nakamura, T., 1986. Energy release rate along a three-dimensional crack front in a thermally stressed body. International Journal of Fracture 30, 79-102. [11] Walters, M.C., Paulino, G.H., Dodds JR., R.H., 2004. Stress-intensity factors for surface cracks in functionally graded materials under mode-I thermo-mechanical loading. International Journal of Solids and Structures, 41(3-4), 1081-1118. [12] Yildirim, B., Dag, S., Erdogan, F., 2005. Three dimensional fracture analysis of FGM coatings under thermo-mechanical loading. International Journal of Fracture, 132(4), 369-395. [13] Yildirim, B. An equivalent domain integral method for fracture analysis of functionally graded materials under thermal stresses. Journal of Thermal Stresses, 29: 371-397, 2006. [14] Ayhan, A.O., 2009. Three-dimensional mixed-mode stress intensity factors for cracks in functionally graded materials using enriched finite elements. International Journal of Solids and Structures, 46(3-4), 796-810. [15] Ayhan, A.O., 2007. Stress intensity factors for three-dimensional cracks in functionally graded materials using enriched finite elements. International Journal of Solids and Structures, 44(25-26), 8579-8599. 16
496 497
498 499 500
501 502 503 504
505 506 507
508 509 510
511 512 513 514
515 516
517 518 519
520 521
522 523 524
525 526 527
528 529
[16] Gogotsi, G.A., 2009. Fracture behaviour of Mg-PSZ ceramics: Comparative estimates. Ceramics International, 35(7), 2735-2740. [17] Mencik, J.,1996. Mechanics of Components with Treated or Coated Surfaces. Kluwer Academic Publishers, Solid Mechanics and its Applications, Vol. 42, London, pp. 360. [18] Gosz, M, Moran, B., 2002. An interaction energy integral method for computation of mixed-mode stress intensity factors along nonplanar crack fronts in three dimensions. Engineering Fracture Mechanics, 69(3), 299-319. [19] Chang, J.H., Wu, D.J., 2007. Stress intensity factor computation along a non-planar curved crack in three dimensions. International Journal of Solids and Structures, 44(2), 371-386. [20] Shaghaghi A, Ghajar R, Alfano M, 2011. Finite element evaluation of stress intensity factors in curved non-planar cracks in FGMs. Mechanics Research Communications, 38, 17-23. [21] Ghajar R., Shaghaghi A., Alfano M., 2011. An improved numerical method for computation of stress intensity factors along 3D curved nonplanar cracks in FGMs. International Journal of Solids and Structure, 48, 208-216. [22] Williams, M.L.,1957. On the stress distribution at the base of a stationary crack. ASME Journal of Applied Mechanics. 24, 109-114. [23] Mori T, Tanaka K, 1973. Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta Materialia, 21(5),571-574. [24] Hill R, 1965. A self-consistent mechanics of composite materials. Journal of the Mechanics and Physics of Solids,13(4), 213-222. [25] Shaghaghi A., Ghajar R., Alfano M., 2010.An ABAQUS implementation of the isoparametric graded finite element. Technical Report. University of Calabria (Available upon request). [26] Yehia A. B.; Shephard M. S., 1985. On the effect of quarter-point element size on fracture criteria. International Journal for Numerical Methods in Engineering, 21(10),1911-1924. [27] Li C., Zou Z., Duan Z., 1999. Stress intensity factors for functionally graded solidcylinders. Engineering Fracture Mechanics, 63(6), 735-749. 17
List of symbols X α(X) k(X) T (X) Vi (X) mech Cmnpq (X) KI , KII , KIII µ(X), κ(X) λ(X), µ(X) J(s) 1,2 J (s) Q J H wp
position vector coefficient of thermal expansion thermal conductivity temperature volume fraction of constituent material i mechanical strain constitutive tensor mixed mode SIFs shear and bulk moduli Lame0 constants J-integral at crack front location s interaction energy integral at crack front location s trasformation tensor Jacobian operator Eshelby’s energy-momentum tensor weight factors for Gauss numerical integration
530
18
Figure 1: Finite domain V along a curved non-planar crack in a functionally graded material. The tubular domain employed for computation of the SIFs and a schematic variation of the q-function along the crack from segment Lc are also displayed.
19
Figure 2: Finite element discretization of the domain V employed in the subsequent finite element computation of the stress intensity factors.
20
Figure 3: Boundary value problems addressed in the present work. (a) Conical crack in an infinite half-space subjected to non-uniform temperature variation. (b) Lens shaped crack subjected to uniform temperature variation. Notice that the dashed line represents the initial temperature set in the model.
21
Figure 4: Evolution of volume fraction (material #2) and material properties as a function of the Y -coordinate and the power exponent n within the graded region of the infinite half space embedding the conical crack. The vertical line corresponds to Y = Y2 .
22
Figure 5: Evolution of volume fraction (material #2) and material properties as a function of the radius r and the power exponent n within the cylindrical domain embedding the lens shaped crack.
23
Figure 6: Finite element discretization of the boundary value problems analyzed in the present work. (a) Finite element mesh of the conical crack. (b) Details of the mesh around the crack surface for the conical crack, illustrating the smallest size of the element incident at a crack front (Le ) as well as the domain size (rd ) employed for numerical computations. (c) Finite element mesh employed for the analysis of the lens shaped crack. A mesh similar to that of the conical crack was employed in the near tip region.
24
Figure 7: Stress intensity factors as a function of the domain size (rd ) for the conical crack embedded in an infinite half-space. The dashed curves pertain to the results obtained using the displacement correlation technique (DCT). The ratio JK /JA is also reported to cross-check the accuracy of the obtained mixed mode SIFs.
25
Figure 8: Evolution of the SIFs with respect to the power exponent n for a conical crack in an infinite half-space subjected to uniform heating from T0 =20 to T =100.
26
Figure 9: Temperature distribution within the infinite half space and for varying values of the power exponent n. The insert shows details of the temperature distribution within the graded region.
27
Figure 10: Evolution of the SIFs with respect to the power exponent n for a conical crack in an infinite half-space subjected to non-uniform heating as displayed in Fig. 9.
28
Figure 11: Evolution of the SIFs with respect to the power exponent n for a lens shaped crack in a cylindrical domain subjected to uniform cooling from T0 =100 to T =20 and under fixed grip conditions.
29
Figure 12: Evolution of the SIFs with respect to the power exponent n for a lens shaped crack in a cylindrical domain subjected to uniform cooling from T0 =100 to T =20 and under free grip conditions.
30
List of symbols X α(X) k(X) T (X) Vi (X) mech Cmnpq (X) KI , KII , KIII µ(X), κ(X) λ(X), µ(X) J(s) 1,2 J (s) Q J H wp
position vector coefficient of thermal expansion thermal conductivity temperature volume fraction of constituent material i mechanical strain constitutive tensor mixed mode SIFs shear and bulk moduli Lame0 constants J-integral at crack front location s interaction energy integral at crack front location s trasformation tensor Jacobian operator Eshelby’s energy-momentum tensor weight factors for Gauss numerical integration
1
Highlights Manuscript Number: Title: Authors:
EFM-D-15-00069 Determination of stress intensity factors for 3D curved non-planar cracks in FGMs subjected to thermal loading A. M. Shaghaghi, M. Alfano
Highlights
3D curved non-planar cracks in FGMs under thermal loading are analyzed using the FE method The procedure does not require a priori knowledge of crack front and crack surface curvatures Effective material properties of the FGM are determined by resorting to micromechanics models SIFs of curved non-planar cracks are significantly affected by gradation of material properties The computational framework is shown to be robust and able to provide consistent results
S0013-7944(15)00416-6 http://dx.doi.org/10.1016/j.engfracmech.2015.07.040 EFM 4786
To appear in:
Engineering Fracture Mechanics
Received Date: Revised Date: Accepted Date:
6 February 2015 30 June 2015 3 July 2015
Please cite this article as: Moghaddam, A.S., Alfano, M., Determination of stress intensity factors of 3D curved non-planar cracks in FGMs subjected to thermal loading, Engineering Fracture Mechanics (2015), doi: http:// dx.doi.org/10.1016/j.engfracmech.2015.07.040
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
2
Determination of stress intensity factors of 3D curved non-planar cracks in FGMs subjected to thermal loading
3
Ali Shaghaghi Moghaddama , Marco Alfanob,∗
1
4 5 6 7
8
a
Young Researchers and Elite Club, Takestan Branch, Islamic Azad University, Takestan, Iran b Department of Mechanical, Energy and Management Engineering, University of Calabria, P. Bucci 44C, 87036 Rende (CS), Italy
Abstract A finite element analysis is presented about 3D curved non-planar cracks in functionally graded materials (FGMs) subjected to steady state temperature gradients. The interaction energy integral is employed in conjunction with a computational strategy that does not require a priori information about crack front and crack surface curvatures. The influence of graded thermal and mechanical properties on mixed mode stress intensity factors (SIFs) is explored in detail. Results highlighted the robustness of the proposed computational framework, which therefore can be effectively deployed in the analysis of thermal fracture in FGMs containing curved non-planar cracks.
10
Keywords: FGMs, non-planar cracks, thermal SIFs, interaction integral, micromechanics models
11
1. Introduction
9
12 13 14 15 16 17 18 19 20 21 22 23
Functionally Graded Materials (FGMs) are nonhomogenous materials whose microstructure and composition are engineered so that to achieve a continuos spatial variation of mechanical and physical properties. A significant body of literature has explored multiple aspects of FGMs, including manufacturing and analysis, as recently reviewed in [1]. Several fabrication techniques have been developed (e.g. spark plasma sintering, plasma spray technique and physical vapor deposition) and analytical and numerical efforts have been carried out in order to sustain experimental investigations and the design of FGMs systems [2, 3, 4, 5, 6, 7, 8]. One of the most important design aspect that has been emphasized in earlier works relates to the integrity of cracked components made up of FGMs subjected to thermomechanical loading. The use of graded materials for thermal protection sys∗ Corresponding author. Tel.: +39 0984 494156 Preprint submitted to Engineering Fracture Mechanics Email address: [email protected] (Marco Alfano)
July 25, 2015
24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
tems, such as first-wall composite materials within nuclear fusion reactors, is very attractive and for this reason the analysis of cracks in FGMs under thermo-mechanical loading has been vigorously pursued in recent times. KC and Kim [9] employed the interaction energy integral to study 2D cracks in FGMs subjected to thermo-mechanical loading. The 3D finite domain integral expression and the two-dimensional and axisymmetric specializations for the energy release rate are given in [10]. Walters et al. [11] studied 3D planar cracks in FGMs under thermo-mechanical loading under mode I using the domain form of the J-integral. Yildirim et al. [12] focused on 3D planar cracks in FGMs subjected to mode I mechanical and transient thermal loading and employed the displacement correlation technique (DCT). Later on, Yildirim [13] studied 2D crack in FGMs under mode I steady-state and transient thermal loading using domain integrals. Ayhan [14, 15] proposed a different approach, based on the use of special enriched finite elements, to compute the SIFs of 3D planar cracks in FGMs under mode I and mixed mode loading. Previous works mostly addressed planar cracks, while for 3D nonplanar cracks, such as conical cracks arising from impact or indentation loading [16, 17], only few efforts can be found in the published literature. For instance Gosz and Moran [18] employed the interaction integral approach for 3D non-planar cracks, but their analysis addressed homogenous materials; similarly, Chang and Wu [19] computed mixed mode SIFs for 3D non-planar cracks in homogenous materials by using the concept of the Jk and GIII integrals. Shaghaghi et al. and Ghajar et al. [20, 21] extended the work carried out in [18] to FGMs. However, to the best of the authors knowledge none of the previous works addressed 3D non-planar cracks under thermal loading. The objective of the present study is to complement existing works by analyzing the mixed mode SIFs of 3D curved nonplanar crack in FGMs under thermal loading. An un-coupled thermo-mechanical finite element analysis is solved where the temperature field is firstly retrieved from a thermal analysis and then is employed as initial condition in a subsequent general static analysis step. In both cases gradation of material properties is accounted for1 . The determination of the effective material properties is herein carried out by resorting to micromechanics models. The self-consistent model is used to determine Young’s modulus and Poisson’s ratio, while the rule of 1
Contrary to fracture in homogeneous materials, in FGMs the SIFs and the J-integral not only depend on external loading and crack geometry, but may be heavily affected by the spatial variation of material properties.
2
87
mixture is employed to quantify the coefficient of thermal expansion and thermal conductivity. Upon completion of the finite element analysis, a post processing procedure enables the evaluation of the J-integral and interaction energy integrals. However, using the interaction energy integral, auxiliary fields need to be determined. Gosz and Moran [18] employed an orthogonal curvilinear coordinate system to set-up auxiliary fields for a 3D non-planar crack in homogenous material. A similar procedure was adopted by Shaghaghi et al. [20] to analyze 3D non planar cracks in FGMs. In both works, analytical descriptions of crack front and crack surface curvatures were required to compute the interaction integral and, in turn, the mixed mode SIFs. The approach employed herein does not require any a priori information about crack front and crack surface curvatures to locate integration stations within the finite element framework. Therefore it allows to analyze arbitrary non-planar cracks with relative ease. This simplified strategy, originally developed by the authors in [21], is herein extended to account for the effect of thermal loading on mixed mode SIFs of 3D nonplanar cracks on FGMs. To demonstrate the capabilities of the proposed computational scheme, the mixed mode SIFs for two selected boundary value problems are determined, i.e. a conical crack in an infinite half space and a lens shaped crack in a cylindrical domain. In order to assess the effectiveness of the computational scheme, results are compared with those obtained using the classical displacement correlation technique (DCT) [12]. The paper is organized as follows. In section 2 the computational procedure employed to extract the SIFs is firstly outlined and the determination of effective material properties of FGMs using micromechanics models is described. In section 3 details about the finite element implementation are provided. In turn, section 4 illustrates the results of two selected boundary value problems involving nonplanar crack configurations, i.e. lens shaped and conical cracks subjected to thermal loading.
88
2. Computational procedures
59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86
89 90 91 92 93 94 95 96
The mixed mode stress intensity factors of 3D non-planar cracks are herein obtained by solving an un-coupled thermo-mechanical finite element analysis. Specifically, the heat transfer analysis is performed considering gradation of thermal conductivity. The temperature field is then employed as initial condition in a general static analysis where gradation of mechanical properties is accounted for. Upon completion of the stress analysis, a post processing stage enables the determination of the J-integral, interaction energy integral and, in turn, the mixed mode SIFs. 3
120
In particular, the 3-D domain form of the J-integral found in [10] is employed to extract energy released (JA ) by the body per unit of crack advance over a finite segment of the crack front. The interaction energy integral is then employed to determine the mixed mode SIFs. Upon evaluation of mixed mode SIFs, the energy release rate is obtained (JK ) and a crosscheck of the results is made by computing the ratio of JK /JA . Since analytical solutions are not available for the problem geometries and boundary conditions analyzed herein, additional comparisons and judgements on the accuracy of the results is also carried out by resorting to the classical displacement correlation technique (DCT) [12]. The DCT has a very simple formulation since it essentially correlates numerical results for displacement at crack front location to evaluate SIFs. The size of the singular elements used around the crack front is an important factor to be considered to obtain sufficiently accurate results in the finite element analysis. Indeed, the accuracy is acceptable for most engineering applications provided quarter-point elements are used in the first layer of elements around the crack tip, and the mesh in the near-tip region is refined enough. On the other hand the domain integral formulation is extremely versatile and relatively simple to implement numerically. It removes the need to precisely capture he details of the singular fields near the crack tip, as a consequence there is no need to use singular elements and/or a very refined mesh. As it will be shown later, the details of the finite element mesh have been properly selected to accommodate the requirements of both methods, to provide consistent results and to ensure an objective comparative analysis.
121
2.1. The interaction energy integral for FGMs
97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119
122 123 124 125 126 127 128 129 130 131 132
Mixed mode stress intensity factors are calculated independently by means of the interaction energy integral. This is done by superimposing suitable auxiliary fields to the actual fields. A finite element procedure which enables the determination of SIFs using the interaction energy integral was originally proposed in [20, 21], and it is herein extended to include the effect of thermal loading on the mixed mode SIFs of 3D non-planar cracks in FGMs. Let consider a schematic depiction of a 3D nonplanar crack within a FGM as shown in Fig.1. Considering a finite domain V surrounding the crack front at location s, the J-integral for graded materials under general thermo-mechanical loading can be written in domain form as follows: Z 1 J(s) = − (Hlj ql,j + Hlj,j ql ) dV , (1) γ V 4
133 134 135 136 137 138
139 140 141 142
where ql is a sufficiently smooth vector field in V. In order to account for thermal loading the Eshelby’s energy-momentum tensor (H) needs to be properly modified. In particular, the mechanical strain, εmech = εtotal mn mn − α (X) T (X), the gradient of the thermal expansion coefficient, α (X),l , as well as the temperature gradient, T (X),l , need to be included in H. It follows that the tensor Hlj and its gradient, Hlj,j , are given by: Hlj
= 12 σmn εmech mn δlj − σij ui,l ,
(2)
Hlj,j
= Cmnpq,l εmech εmech pq mn − σpp {α,l T + αT,l }.
(3)
Notice that Cmnpq (X) is the spatially varying constitutive tensor for graded materials; α = α(X) is the coefficient of thermal expansion; T = T (X) is the temperature within the body; X is the spatial coordinate vector, while γ is given by the following equation: Z γ=
Lc 2
∆a(ξ)dξ,
(4)
− L2c
145
where ∆a(ξ) is the magnitude of the crack advance at crack front segment Lc , as shown in Fig.1. On the other hand the interaction energy integral is given by: Z 1 (1,2) (1,2) (1,2) J (s) = − (5) Hlj ql,j + Hlj,j ql dV , γ V
146
where
143 144
(1,2)
(1) (2)
(1,2)
Hlj,j 147
149 150 151 152 153
(2) (1)
(2) mech(1)
= Cmnpq,l εpq εmn
(2)
(6) (2)
(1)
− σpp {α,l T + α T,l } − σmn,n um,l .
(7)
Computation of SIFs can be effectively carried out by the following equation: 2 (1 + ν ) 2 1 − νs2 (1) (2) (1) (2) (1) (2) s (1,2) J (s) = KI KI + KII KII + KIII KIII , (8) Es Es (1)
148
(1) (2)
= σmn εmn δlj − {σij ui,l + σij ui,l },
Hlj
(2)
where Ki and Ki (i=I,II,III) are the SIFs for actual and auxiliary fields, respectively; Es and νs represent Young’s modulus and Poisson’s (1) ratio at crack front location s. In turn, SIFs for actual state, e.g. KI , can be obtained by setting the counterpart SIF for auxiliary state equal to (2) (2) (2) KI = 1, KII = 0 and KIII = 0. It follows that Eq. (8) yields the required SIF: Es (1) KI = J (1,2) (s) . (9) 2 (1 − νs2 ) 5
154 155 156 157
2.2. Auxiliary fields The Williams solution for the 2D asymptotic crack tip displacement fields [22], is used to set-up the auxiliary displacement and strain fields in the crack tip local coordinate system: (2) u ¯j
158
(2) r
K = I 2µ
(2) r
K r I gj (θ, ν)+ II 2π 2µ
and (2) ε¯ij
159 160 161 162
163 164 165 166
(2)
K r II gj (θ, ν)+ III 2π 2µ
1 = 2
r
r III g (θ, ν) , (10) 2π j
(2) !
(2) ∂u ¯j ∂u ¯i + ∂xj ∂xi
,
(11)
where gj (θ, ν) is the standard angular function. In the case of FGMs the effect of material gradation needs to be taken into account. In this work the auxiliary stress field is related to the auxiliary strain field through the constitutive tensor as: (2) (2) σ ¯ij = Cijkl (X) ε¯kl . (12) To be consistent with actual fields, the auxiliary fields are transformed to cartesian coordinate system. Considering the transformation tensor Q which maps the local basis into the cartesian one, the auxiliary displacement, strain and stress fields can be written as: ¯, u = QT u
(13)
T
¯ Q, ε=Q ε ¯ Q. σ = QT σ 167 168 169 170 171 172 173 174 175 176 177 178 179
Since the auxiliary fields are now referred to a fixed local coordinate system, it is observed that strain-displacement compatibility conditions is satisfied and auxiliary stresses are equilibrated. Note that due to the spatially varying constitutive tensor required to represent graded materials, the auxiliary stress field is not in equilibrium and this has to be considered. Gosz and Moran [18] presented an integral method for the analysis of mixed-mode SIFs along 3D non-planar crack fronts in homogeneous materials. They highlighted that for a proper definition of the auxiliary fields a curvilinear coordinate system should be employed, and crack surface and crack front curvatures need to be accounted for. In particular, the gradient of auxiliary displacement field requires a knowledge of the divergence of the auxiliary stress field as well as the location of integration points with respect to crack front. For this reason, analytical equations describing crack front
6
192
and crack surface were deployed. It should be noted the curvilinear coordinates obtained using analytical equations can yield incorrect locations of integration points, thereby leading to auxiliary fields inconsistent with actual fields generated by the finite element mesh. For this reason, Ghajar et al. [21] recently extended the work by Gosz and Moran [18] to FGMs and proposed an improved computational approach based on the use of a local coordinate system located at the crack front. The proposed method does not require a priori knowledge of the geometry of the crack in order to locate integration stations with respect to the crack front. Its accuracy and relatively simple implementation were demonstrated through the analysis of 3D curved non planar cracks in FGMs. This approach is herein employed in order to include the capabilities to analyze non-planar cracks in FGMs under thermo-mechanical loading.
193
2.3. Determination of material properties using micromechanics models
180 181 182 183 184 185 186 187 188 189 190 191
194 195 196 197 198 199 200
201 202 203 204 205 206
The determination of the effective material properties is carried out by resorting to micromechanics models, whose inputs are the volume fractions of the constituent materials. Micromechanics models, such as the MoriTanaka [23], the self-consistent [24] and rule of mixture, allow to determine the effective material properties by considering the volume fraction of constituents. In this work the rule of mixture is employed to obtain the coefficient of thermal expansion, α, and thermal conductivity, k, as follows: α = V1 α1 + V2 α2 ,
(14)
k = V1 k1 + V2 k2 ,
(15)
where Vi = Vi (X), ki = ki (X) and αi = αi (X) are volume fraction, thermal conductivity and thermal expansion coefficient of the material constituent i, which are all functions of the spatial coordinates. For estimation of Young’s modulus and Poisson’s ratio, the self-consistent model is employed. Specifically, the shear and bulk moduli (µ, κ) are given as a function of the constituents [24]: 1 V1 2V1 = + , κ + 4µ/3 κ1 + 4µ/3 κ2 + 4µ/3 V1 κ1 V 2 κ2 V 1 µ2 V2 µ 1 + +5 + + 2 = 0. κ1 + 4µ/3 κ2 + 4µ/3 µ − µ2 µ − µ1
207 208
(16)
The above equations are firstly solved at a specific spatial location for µ and κ, then Young’s modulus and Poisson’s ratio of the FGMs are obtained 7
209
using the following equations: E=
210 211 212 213
9µκ 3κ − 2µ , ν= . µ + 3κ 2 (µ + 3κ)
(17)
For graded materials the gradient of material properties and derivatives of material constitutive tensor, Cmnpq (X),l , are incorporated in the J-integral and interaction integral. The constitutive tensor for isotropic elastic FGMs is given by: Cmnpq (X) = λ (X) δmn δpq + µ (X) (δmp δnq + δmq δnp ) ,
214 215
where δmn is the Kronecker delta, and µ (X) and λ (X)are spatially varying Lame0 constants defined as follows: λ (X) =
216 217 218 219
E (X) ν (X) , (1 + ν (X)) (1 − 2ν (X))
221 222 223 224 225 226 227 228 229 230 231 232 233 234
µ (X) =
E (X) , 2 (1 + ν (X))
(19)
E(X) and ν(X) are modulus of elasticity and Poisson’s ratio which are determined using the self-consistent model. The derivative of constitutive tensor in Cartesian coordinate system is given in [21] and is recalled here for the sake of completeness: Cmnpq (X),l =
220
(18)
∂Cmnpq (X) ∂E (X) ∂Cmnpq (X) ∂ν (X) + (l = 1, 2, 3). (20) ∂E (X) ∂Xl ∂ν (X) ∂Xl
3. Numerical aspects 3.1. Finite element discretization of J-integral and interaction energy integral The actual fields are extracted from finite element analysis carried out using ABAQUS/Standard. A heat transfer analysis is firstly performed to obtain temperature distribution in the body. This step has been performed using the subroutine UMATHT available in ABAQUS/Standard which allows to account for the gradation of k. In turn, output model temperatures are set as initial conditions in the subsequent mechanical analysis, where gradation of mechanical properties is incorporated through subroutine UMAT. Upon completion of the thermo-mechanical analysis a post processing step is carried out to extract mixed mode SIFs. In the finite element context the domain V embeds a few elements surrounding the crack tip at crack front location s, which is in turn represented by a node in the finite element mesh (see Fig. 2). Having defined the finite 8
235 236
domain V the discretized versions of J-integral and interaction integrals are obtained as: elems gps 1 X X J(s) = [Hlj ql,j ]p + [Hlj,j ql ]p detJ wp , (21) γ p J
237 238 239 240 241 242 243 244 245
(1,2)
1 (s) = γ
V gps elems X X h V
(1,2) Hlj ql,j
p
i p
+
h
(1,2) Hlj,j ql
i p
detJ wp ,
where summation is done over the integration points (gps) for each element included in the domain V . Parameters wp and detJ denote the weight factors and the Jacobian determinant, respectively. The effect of element type and number of Gaussian points on the accuracy of finite element results have been studied in [25]. It was shown that 20-noded brick elements with 2×2×2 Gauss quadrature rule are well suited for modeling 3D FGMs. Accordingly, a similar approach was undertaken in this work. Material properties and their derivatives at integration points are computed using an isoparametric approach, such that for a given generic field variable (Φ): Φ (X) =
20 X
Φ i hi ,
(22)
∂hi , ∂Xi
(23)
i=1
∂Φ (X) = ∂Xi 246 247
248 249 250 251 252 253 254
20 X i=1
Φi
where Φ (X) = {E (X) , ν (X) , α (X) , k (X) , T (X)}, Φi is the nodal values and hi is the shape function at node i. 3.2. Evaluating the q function In order to evaluate the q function the finite domain V drawn in Fig. ¯ on crack front segment Lc is ∆a(ξ) in the 2 is considered. The values of q ¯1 while it is zero on the boundary of the considered domain. As direction e ¯ is set to zero for the nodes indicated by ◦ shown in Fig. 2, the value of q and unity for those indicated by •. It follows that for an element inside the ¯ function at node i is given by: domain V , the q 0 for nodes on boundary of domain V i q¯l = (24) ¯1 e for inner nodes (i)
255 256
¯ at node i in local coordinate system. where q¯l are the components of q ¯ is then mapped to the cartesian coordinate system as follows: The vector q ¯, q = QT q 9
(25)
257 258
and consistent with isoparametric finite element formulation, ql and ql,j at integration points are interpolated using standard shape functions: ql =
20 X
(i)
hi ql ,
(26)
i=1
ql,j =
20 X
(i)
hi,j ql .
i=1
259
260 261 262 263 264 265 266 267 268 269 270
4. Numerical examples 4.1. Boundary value problems selected for procedure verification In this section the developed numerical procedures are applied to the analysis of selected boundary value problems involving a conical crack and a lens shaped crack in three-dimensional bodies made up of FGMs. A conical crack in an infinite half space is shown in Fig. 3(a). The gradation of mechanical and thermal properties occurs in the Y -direction and extends from zero to Y2 . The material properties at top surface of the domain (Y =0) are E1 = 5, ν1 = 0.2, α1 = 1, k1 = 1, while at the bottom of the graded region (Y = Y2 ) E2 = 1, ν2 = 0.4, α2 = 5, k2 = 5. Then the volume fraction of constituent materials is considered to vary in the Y -direction according to a power function with exponent n as follows: ( n |y| if y < y2 , y2 V2 = (27) 1 otherwise, V1 = 1 − V2 .
271 272 273 274 275 276 277 278 279 280 281 282 283
The volume fraction V2 , Young’s modulus, Poisson’s ratio, coefficient of thermal expansion and thermal conductivity of the FGMs are drawn in Fig. 4 for varying values of the power exponent n. It is observed that for increasing values of the exponent n, the volume fraction of the stiffer material (i.e. material #1) is mostly dominant in the graded region. For increasing values of n the softer material (i.e. material #2) is introduced toward the end of the graded region. Notice that outside the graded region, that is Y > Y2 , the material is assumed to be homogenous and isotropic with properties equal to those of material #2. Therefore for higher n the graded region exhibits mostly the behavior of material #1 and the rate of change of V2 across the boundary, Y = Y2 , is quite sharp. On the other hand for lower n, for instance n = 0.2, the volume fraction of material #2 quickly increases and the transition across the boundary is rather smooth. It is 10
284 285 286 287 288 289 290 291 292 293 294
295 296 297 298 299 300 301 302
303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318
then concluded that lower values of the index n lead to a smooth variation of material properties. As it will be shown in the next section, the conical crack is studied considering both uniform and non-uniform temperature loading and free boundary conditions (free grips). Fig. 3(b) illustrates a cylinder of radius r and height H with an embedded lens-shaped crack whose surface geometry is characterized by radius R and angle β. In particular the geometric parameters of the model are given as follows: 2r /H=1, 2r /R=10, R=1 and β = π/4. The material properties at r=0 are considered to be E1 = 1, ν1 = 0.4, α1 = 5, k1 = 5, while at the outer surface of the cylinder are E2 = 5, ν2 = 0.2, α2 = 1, k2 = 1. The volume fraction varies as a power function with exponent n: r n V2 = , (28) 5 V1 = 1 − V2 . (29) In Fig. 5 the volume fraction V1 , Young’s modulus (E) and coefficient of thermal expansion (α) are also shown. In addition, the product (E · α) is displayed. This parameter is considered because the thermal strain in the cylinder, εth is proportional to α and, in turn, the associated normal stress scales with E · εth . As a consequence it provides an approximate measure of the normal stress acting on the cylinder. Finally, notice that boundary conditions include uniform temperature variation while both fixed and free grips conditions will be analyzed. 4.2. Finite element models The finite element models of the conical and lens shaped cracks are given in Fig. 6. In particular, a conical crack is embedded in a cubic domain whose dimensions are 600 × 600 × 300 units length and which is shown in Fig. 6(a). In order to decrease the computational costs, only one quarter of the domain has been included in the model. The crack circumscribes a circle with a radius r=10 units length on the free surface (cfr. Fig. 3(a)), and intersects this last at an angle of 45◦ . In addition it extends 15 units length into the solid domain. The mesh shown in the figure comprises around 9500 brick elements (C3D20) and 42000 nodes. Notice that a sensitivity analysis to the number of elements along the crack front was already performed in [21]. The results obtained with finite element meshes comprising six, eight and ten elements along the 90-deg sector were compared with an analytical solution. The results indicated that fairly accurate results could be obtained using only six elements along crack front. In the present study an even higher number of elements was 11
330
deployed, i.e. 18 elements are inserted along crack front and 16 elements are surrounding the crack tip as shown in Fig. 6(b). The size of the elements incident at crack front, Le , is such that Le /a = 1/50. On the other hand, the finite element discretization of the lens shaped crack is given in Fig. 6(c). Once again, symmetry conditions were exploited and only one quarter of the cylinder was considered for modeling purposes. The mesh comprised around 2500 brick elements (C3D20) and the size of the smallest element incident at crack front was equal to Le /R = 0.01. Notice that the finite element discretization around the crack front is similar to that employed for the conical crack. However, the work by Gosz and Moran indicated that for a lens shaped crack [18] ten elements along a 90-deg sector may suffice, a similar choice was also made herein.
331
4.3. Results and discussion
319 320 321 322 323 324 325 326 327 328 329
332 333 334 335 336 337 338 339 340 341 342
343 344 345 346 347 348 349 350 351 352 353
Since reference solutions are not available for the boundary value problems analyzed in the present work, the results were compared with those obtained through an objective and consistent use of the DCT. From this standpoint, a proper selection of the mesh details around the crack tip region is very important, since the size of singular elements around the crack front is an important factor in order to obtain accurate results. A in-depth discussion on this point was already provided in previous related works [3, 26, 27], and it was shown that element size ranging from 1/20a to 1/32a, where a is the crack length, can provide accurate results. Similarly, preliminary convergence analyses carried out in the present work indicated that element size equal to 1/50a was able to provide satisfactory results. 4.3.1. Procedure verification considering a conical crack under uniform temperature change The results pertaining to the conical crack are illustrated first. Notice that the body is initially assumed to be at a uniform temperature T0 = 20 and subsequently it is heated up to a uniform temperature T = 100 (cfr. Fig. 3(a)). For the sake of verification, the normalized2 mixed mode SIFs obtained using the domain integral for different domain sizes, rd /Le , and the DCT approach were compared as shown in Fig. 7. Computations were made assuming n = 1. The deviations among the two sets of results is within 2% thereby testifying the robustness of the developed computational procedures. An additional cross-check was made by analyzing the ratio 2
Mixed mode SIFs were normalized against K0 =
12
√ E1 α T πa. 1−ν1 1 0
354 355 356 357 358 359 360 361 362 363 364 365 366
367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392
between the energy release rate obtained from stress intensity factors, Jk , and that stemming from the domain integral, JA , which is also given in Fig. 7. Again, the deviation from unity is less than 1% which implies that the results are consistent and domain independent. The sensitivity of the SIFs to the power exponent n was also analyzed and the results are shown in Fig. 8. The normalized SIFs, KIN and KIIN , display an increase as the power exponent n increases. This effect is also explained by the above discussion concerning the trend of V2 with respect to power exponent n. In particular, the induced material gradation, see Fig. 4, enables an increase in the crack face opening displacements as the power exponent n increases. However, both SIFs consistently increases and, as a consequence, the mode mixity is almost unchanged in the investigated range of n. 4.3.2. Steady state temperature gradients To study the influence of thermal conductivity on temperature profile in the infinite half-space containing a conical crack, it is assumed that the body is exposed to T = 100 on the upper surface, Y = 0, and to T = 20 on the lower surface, Y = 300. In this circumstance, a heat transfer analysis is first run and temperature distribution is extracted. Then temperature distribution is imported to a static analysis considering initial temperature T0 = 20 and mixed mode SIFs are in turn obtained. The temperature variation which is a function of Y is shown in Fig. 9. It is observed that as n increases, the temperature profile deviates from linearity, especially around the transition boundary located at Y2 . The temperature distribution is then imported to the general static analysis and mixed mode SIFs are calculated. The results for KIN and KIIN are shown in Fig. 10. The trend is quite similar to that of uniform expansion as Fig. 8. As the power exponent n increases, KIN and KIIN increase. It is important to highlight that negative values of KIN were also observed for n ≤0.3, which are induced by the occurrence of negative opening displacement, i.e. crack closure. In general terms, crack closure can be associated to a variety of mechanisms, for instance it can be induced by the applied loading whenever compression stresses near a crack are generated. On the other hand, crack closure might also be geometry induced because of the specific shape of the crack, such as kinked cracks. It might be also generated by material inhomogeneity. In particular, it is speculated that in the present work the observed crack closure is caused by the combined effects of thermal loading, crack geometry and the non-homogeneity brought by FGMs. Indeed, the material gradation profile strongly affects the stress and strain 13
413
fields around a crack. This point was not further analyzed herein because in doing so the effect of crack surface contact and friction should be taken into account. In particular, the treatment of crack closure would require to heavily modify the proposed formulation. For instance, the solution of frictional contact problems is strongly history dependent and tailored incremental solution procedures would be needed. Moreover, from the analysis of Fig. 10 it is possible to infer that crack opening displacements are almost always positive for the considered model material system and boundary conditions. For the lens shaped crack the initial temperature of the body is assumed to be T0 = 100 and then it is uniformly decreased to T = 20. At first it has been considered that the body is gripped at the upper and lower surfaces, and the corresponding normalized mixed mode SIFs3 are shown in Fig. 11. It is inferred that KI and KII follow similar trends, i.e. they consistently increase until n = 0.5 is reached and then decrease. A possibile explanation for such a trend can be given if one consider that near the crack front the parameter (E · α) increases as the power exponent decreases, see Fig. 5. Free boundary conditions have been also explored, and the body could freely contract in this case. The mixed mode SIFs are given in Fig. 12 and quite similar trend to that of fixed grip was observed. Although the trend with the power exponent is still the same, the SIFs are now smaller than that obtained on fixed grip boundary conditions.
414
5. Conclusion
393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412
415 416 417 418 419 420 421 422 423 424 425 426 427
In the present work 3D non-planar cracks in FGMs under thermal load were analyzed using the domain form of the J-integral and interaction energy integrals. Gradation of both thermal and mechanical properties was considered simultaneously. Micromechanics models were deployed in order to determine the effective properties of the FGM, and the volume fraction of constituent materials was assumed to display a power function form. The overall computational approach consisted of a few interrelated steps. At first, an heat transfer analysis was performed considering gradation of thermal conductivity. The temperature field was then employed as initial condition in a general static analysis where gradation of mechanical properties was accounted for. Upon completion of the stress analysis, a post processing stage enabled the determination of the J-integral, interaction energy integral and, in turn, the mixed mode SIFs. Two selected boundary value problems 3
Normalization was made against K0 =
√ 2 E1 α T πa. π 1−ν1 1 0
14
445
are analyzed, i.e. a conical crack in an infinite half space and a lens shaped crack in a cylindrical domain. In order to assess the effectiveness of the computational scheme, a comparative analysis was carried out by using the classical displacement correlation technique (DCT). The results indicated that the proposed procedure is domain independent and that the mixed mode SIFs are in a good agreement with those obtained using the DCT. However, the domain integral formulation proposed herein is extremely versatile and relatively simple to implement numerically. Moreover, if compared to the DCT, it removes the need to precisely capture the details of the singular fields near the crack tip. Therefore, singular elements and/or very refined meshes, which are computationally expensive when conducting three-dimensional analyses, are not needed. In addition, a parametric study was carried out to reveal the effect of material gradation profile on the SIFs. It was concluded that the mixed mode SIFs of curved non planar cracks are considerably influenced by material gradation. The obtained results testify that the proposed procedure is an accurate and useful tool for the analysis of 3D non planar cracks in FGMs which can assist the design of FGMs against fracture.
446
References
428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444
447 448 449
450 451 452
453 454 455
456 457 458
459 460 461
[1] Birman, V., Byrd, L.W., 2007. Modeling and analysis of functionally graded materials and structures. Applied Mechanics Reviews, 60(1-6), 195-216. [2] Dolbow, J.E., Gosz, M., 2002. On the computation of mixed mode stress intensity factors in functionally graded materials. International Journal of Solids and Structures, 39(9), 2557-2574. [3] Kim, J.-H, Paulino, G.H., 2002. Finite element evaluation of mixed mode stress intensity factors in functionally graded materials. International Journal for Numerical Methods in Engineering, 53(8), 1903-1935. [4] Ghajar, R., Shaghaghi, A., 2011. Numerical investigation of the mode III stress intensity factors in FGMs considering the effect of graded Poisson’s ratio. Engineering Fracture Mechanics, 78, 1478-1486. [5] Shaghaghi A, Ghajar R, Alfano M, 2013. Determining the mixed mode stress intensity factors of surface cracks in functionally graded hollow cylinders. Materials and Design, 43, 475-484
15
462 463 464 465
466 467 468 469
470 471 472 473
474 475 476
477 478 479
480 481 482 483
484 485 486
487 488 489
490 491 492
493 494 495
[6] Nikishkov, G.P., Atluri, S.N., 1987. Calculation of fracture mechanics parameters for an arbitrary three-dimensional crack, by the equivalent domain integral method. International Journal for Numerical Methods in Engineering, 24(9), 1801-1821. [7] Yu, H., Wu, L., Guo, L., Wu, H., Du, S., 2010. An interaction integral method for 3D curved cracks in nonhomogeneous materials with complex interfaces. International Journal of Solids and Structures, 47(16), 2178-2189 [8] Walters, M.C., Paulino, G.H., Dodds JR., R.H., 2006. Computation of mixed-mode stress intensity factors for cracks in three-dimensional functionally graded solids. Journal of Engineering Mechanics, 132(1), 1-15. [9] KC A, Kim J-H, 2008. Interaction integrals for thermal fracture of functionally graded materials. Engineering Fracture Mechanics 75, 25422565. [10] Shih, C.F., Moran, B., Nakamura, T., 1986. Energy release rate along a three-dimensional crack front in a thermally stressed body. International Journal of Fracture 30, 79-102. [11] Walters, M.C., Paulino, G.H., Dodds JR., R.H., 2004. Stress-intensity factors for surface cracks in functionally graded materials under mode-I thermo-mechanical loading. International Journal of Solids and Structures, 41(3-4), 1081-1118. [12] Yildirim, B., Dag, S., Erdogan, F., 2005. Three dimensional fracture analysis of FGM coatings under thermo-mechanical loading. International Journal of Fracture, 132(4), 369-395. [13] Yildirim, B. An equivalent domain integral method for fracture analysis of functionally graded materials under thermal stresses. Journal of Thermal Stresses, 29: 371-397, 2006. [14] Ayhan, A.O., 2009. Three-dimensional mixed-mode stress intensity factors for cracks in functionally graded materials using enriched finite elements. International Journal of Solids and Structures, 46(3-4), 796-810. [15] Ayhan, A.O., 2007. Stress intensity factors for three-dimensional cracks in functionally graded materials using enriched finite elements. International Journal of Solids and Structures, 44(25-26), 8579-8599. 16
496 497
498 499 500
501 502 503 504
505 506 507
508 509 510
511 512 513 514
515 516
517 518 519
520 521
522 523 524
525 526 527
528 529
[16] Gogotsi, G.A., 2009. Fracture behaviour of Mg-PSZ ceramics: Comparative estimates. Ceramics International, 35(7), 2735-2740. [17] Mencik, J.,1996. Mechanics of Components with Treated or Coated Surfaces. Kluwer Academic Publishers, Solid Mechanics and its Applications, Vol. 42, London, pp. 360. [18] Gosz, M, Moran, B., 2002. An interaction energy integral method for computation of mixed-mode stress intensity factors along nonplanar crack fronts in three dimensions. Engineering Fracture Mechanics, 69(3), 299-319. [19] Chang, J.H., Wu, D.J., 2007. Stress intensity factor computation along a non-planar curved crack in three dimensions. International Journal of Solids and Structures, 44(2), 371-386. [20] Shaghaghi A, Ghajar R, Alfano M, 2011. Finite element evaluation of stress intensity factors in curved non-planar cracks in FGMs. Mechanics Research Communications, 38, 17-23. [21] Ghajar R., Shaghaghi A., Alfano M., 2011. An improved numerical method for computation of stress intensity factors along 3D curved nonplanar cracks in FGMs. International Journal of Solids and Structure, 48, 208-216. [22] Williams, M.L.,1957. On the stress distribution at the base of a stationary crack. ASME Journal of Applied Mechanics. 24, 109-114. [23] Mori T, Tanaka K, 1973. Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta Materialia, 21(5),571-574. [24] Hill R, 1965. A self-consistent mechanics of composite materials. Journal of the Mechanics and Physics of Solids,13(4), 213-222. [25] Shaghaghi A., Ghajar R., Alfano M., 2010.An ABAQUS implementation of the isoparametric graded finite element. Technical Report. University of Calabria (Available upon request). [26] Yehia A. B.; Shephard M. S., 1985. On the effect of quarter-point element size on fracture criteria. International Journal for Numerical Methods in Engineering, 21(10),1911-1924. [27] Li C., Zou Z., Duan Z., 1999. Stress intensity factors for functionally graded solidcylinders. Engineering Fracture Mechanics, 63(6), 735-749. 17
List of symbols X α(X) k(X) T (X) Vi (X) mech Cmnpq (X) KI , KII , KIII µ(X), κ(X) λ(X), µ(X) J(s) 1,2 J (s) Q J H wp
position vector coefficient of thermal expansion thermal conductivity temperature volume fraction of constituent material i mechanical strain constitutive tensor mixed mode SIFs shear and bulk moduli Lame0 constants J-integral at crack front location s interaction energy integral at crack front location s trasformation tensor Jacobian operator Eshelby’s energy-momentum tensor weight factors for Gauss numerical integration
530
18
Figure 1: Finite domain V along a curved non-planar crack in a functionally graded material. The tubular domain employed for computation of the SIFs and a schematic variation of the q-function along the crack from segment Lc are also displayed.
19
Figure 2: Finite element discretization of the domain V employed in the subsequent finite element computation of the stress intensity factors.
20
Figure 3: Boundary value problems addressed in the present work. (a) Conical crack in an infinite half-space subjected to non-uniform temperature variation. (b) Lens shaped crack subjected to uniform temperature variation. Notice that the dashed line represents the initial temperature set in the model.
21
Figure 4: Evolution of volume fraction (material #2) and material properties as a function of the Y -coordinate and the power exponent n within the graded region of the infinite half space embedding the conical crack. The vertical line corresponds to Y = Y2 .
22
Figure 5: Evolution of volume fraction (material #2) and material properties as a function of the radius r and the power exponent n within the cylindrical domain embedding the lens shaped crack.
23
Figure 6: Finite element discretization of the boundary value problems analyzed in the present work. (a) Finite element mesh of the conical crack. (b) Details of the mesh around the crack surface for the conical crack, illustrating the smallest size of the element incident at a crack front (Le ) as well as the domain size (rd ) employed for numerical computations. (c) Finite element mesh employed for the analysis of the lens shaped crack. A mesh similar to that of the conical crack was employed in the near tip region.
24
Figure 7: Stress intensity factors as a function of the domain size (rd ) for the conical crack embedded in an infinite half-space. The dashed curves pertain to the results obtained using the displacement correlation technique (DCT). The ratio JK /JA is also reported to cross-check the accuracy of the obtained mixed mode SIFs.
25
Figure 8: Evolution of the SIFs with respect to the power exponent n for a conical crack in an infinite half-space subjected to uniform heating from T0 =20 to T =100.
26
Figure 9: Temperature distribution within the infinite half space and for varying values of the power exponent n. The insert shows details of the temperature distribution within the graded region.
27
Figure 10: Evolution of the SIFs with respect to the power exponent n for a conical crack in an infinite half-space subjected to non-uniform heating as displayed in Fig. 9.
28
Figure 11: Evolution of the SIFs with respect to the power exponent n for a lens shaped crack in a cylindrical domain subjected to uniform cooling from T0 =100 to T =20 and under fixed grip conditions.
29
Figure 12: Evolution of the SIFs with respect to the power exponent n for a lens shaped crack in a cylindrical domain subjected to uniform cooling from T0 =100 to T =20 and under free grip conditions.
30
List of symbols X α(X) k(X) T (X) Vi (X) mech Cmnpq (X) KI , KII , KIII µ(X), κ(X) λ(X), µ(X) J(s) 1,2 J (s) Q J H wp
position vector coefficient of thermal expansion thermal conductivity temperature volume fraction of constituent material i mechanical strain constitutive tensor mixed mode SIFs shear and bulk moduli Lame0 constants J-integral at crack front location s interaction energy integral at crack front location s trasformation tensor Jacobian operator Eshelby’s energy-momentum tensor weight factors for Gauss numerical integration
1
Highlights Manuscript Number: Title: Authors:
EFM-D-15-00069 Determination of stress intensity factors for 3D curved non-planar cracks in FGMs subjected to thermal loading A. M. Shaghaghi, M. Alfano
Highlights
3D curved non-planar cracks in FGMs under thermal loading are analyzed using the FE method The procedure does not require a priori knowledge of crack front and crack surface curvatures Effective material properties of the FGM are determined by resorting to micromechanics models SIFs of curved non-planar cracks are significantly affected by gradation of material properties The computational framework is shown to be robust and able to provide consistent results