Investigation of stress intensity factor for internal cracks in FG cylinders under static and dynamic loading

This paper investigates the variations of mode I stress intensity factor (KI) for inner penny-shaped and circumferential cracks in functionally graded solid and hollow thick walled cylinders, respectively with the changes of crack geometry, material gradation and loading conditions. The functionally graded material of cylinders consists of epoxy and glass. It is assumed that the mechanical properties vary with a power law in the radial direction of cylinders. Micromechanical models for conventional composites are used to estimate the material properties of functionally graded cylinders. The equations of motion obtained from the extended finite element discretization are solved by the Newmark method in the time domain. The interaction integral method is employed to calculate the mode I stress intensity factor (KI). The MATLAB programming environment was implemented to solve the problem.


INTRODUCTION
unctionally graded materials (FGMs) are a new class of composite materials in which their composition and structure gradually change resulting in a corresponding change in their properties. These materials are used extensively in applications such as thermal barrier coatings are exposed to severe stress gradients induced by F thermal and/or mechanical loading [1]. FGMs can be used in the construction of tanks and cylindrical furnace such as cement kiln. Circumferential cracks are occasionally developed in cylindrical structures during service or production. These cracks are threats to the safety and reliability of these structures. Subsequent fracture and fatigue analysis of such cracks is of great interest, and requires the determination of stress intensity factors. The stress intensity factor (SIF) is an important parameter to determine the safety of a cracked part. Although several stress intensity factor handbooks have been published, the available solutions of stress intensity factors for pressure vessels are not always adequate for particular engineering applications. [2] Solutions to the problem concerning a circumferential crack in a circular cylinder made of homogeneous or composite materials are relatively few. Closed form stress intensity factor of an arbitrarily located inner-surface circumferential crack in an edge-restraint homogeneous cylinder under linear radial temperature distribution derived by Meshii and Watanabe [3]. Chen [4] evaluated stress intensity factors in a cylinder with a circumferential crack using the computing compliance and the finite difference methods to solve the boundary value problem. Lee [5] analyzed the stress distribution in a long circular cylinder of isotropic elastic material with a circumferential edge crack under uniform shearing stress and determined the stress intensity factor. Jones [6] applied the impulse response method to the analysis of the thermally striped internal surface of a hollow cylinder containing a circumferential crack on this surface. He calculated stress intensity factor and strain energy density factor ranges as functions of crack depth for various sinusoidal striping frequencies. Moulick and Sahu [2] derived weight functions for the surface and the deepest point of an internal semi elliptical crack in a thick-wall cylinder. Tran and Geniaut [7] developed an extended finite element method (X-FEM) axisymmetric model and employed it to compute stress intensity factors for cracked industrial specimens and components. They used the X-FEM model to assess the integrity and durability of a cracked rotor coil retaining ring during the power plant operation. Wu et al. [8] described a three-dimensional domain integral method for extracting mixed-mode stress intensity factors. Predan et al. [9] calculated the stress-intensity factor for the circumferential semielliptical surface cracks in a hollow cylinder cross section under torsion using a finite element technique. Sharma [10] et al. used the extended finite element method to evaluate the stress intensity factors of a semi-elliptical part through thickness axial/circumferential crack. The pipe or pipe-bend having a crack on the outer surface was subjected to internal pressure or opening bending moment in their research. Seifi [11] determined the stress intensity factors for internal surface cracks in autofrettaged functionally graded cylinders. Eshraghi and Soltani [12] obtained stress intensity factors for functionally graded cylinders with internal circumferential cracks using the weight function method for different combinations of cylinder geometry, crack depth, and material gradation. In this paper, two problems are considered; the first is the determination of the stress intensity factor (KI) of an inner penny-shaped crack in a FG rod loaded by uniform axial tension and the second problem is the calculation of the SIF for an inner circumferential crack in a FG thick walled cylinder under uniform axial tension. Both problems are solved under static and dynamic loading conditions. The extended finite element method is used in this study, to compute the stress and displacement fields necessary for determining the stress intensity factor. Stress intensity factors are obtained using the interaction integral method. The Newmark time integration scheme is used to solve the dynamical system of matrix equations obtained from the spatial discretization of equations of motion. The effects of crack radius, rotational speed of cylinders and material gradation on SIFs are studied. The programming is done in MATLAB.

MODELING OF FUNCTIONALLY GRADED CYLINDER
n the present study, we assume that the material properties change along the radius of the cylinder and the volume fraction of inclusion i V follows a simple power function, where R and W are inner radius and thickness of cylinder, respectively. p is the power exponent determining the volume fraction profile. The volume fraction of matrix m V is obtained as below.
We assume that the functionally graded cylinder is made of epoxy-phase and glass-phase. In this study, micromechanical models for conventional composites [13,14] are employed to calculate the properties of functionally graded cylinder. According to these models, the shear modulus of rigidity ( ) and bulk modulus ( K ) of FGM cylinder are obtained by the following equations.
In these equations subscript i and m refer to inclusion and matrix, respectively. Young's modulus ( E ), Poisson's ratio ( ) and density of FGM cylinder are calculated as follows.
To incorporate these formulas into the extended finite element model, the value of each material property is computed at each individual node based on micromechanical models. Utilizing the generalized isoparametric graded finite elements, introduced by Kim and Paulino [15], material properties gradation is considered in an element. In this method, material properties such as elastic modulus ( E ), Poisson's ratio ( ), and mass density ( ) at Gauss points can be interpolated using shape functions from their nodal point values as follows [15]:

GENERAL PROBLEM FORMULATION
ith regard to the geometry of cracked cylinders and loading conditions, considered problems are axisymmetric. In this kind of problems, material properties and forces are constant in circumferential direction; also, displacement along -axis (u ) vanishes. The general governing equations are the equations of motion in cylindrical coordinates, which are simplified as follows for axisymmetric problems [16].
and are radius, density, displacement, body force, time, normal and shear stress, respectively. Subscripts r and z refer to radial and longitudinal directions of cylinders, respectively. Also, in axisymmetric problems we have [17] r z To solve equations of motion, discontinuous-Galerkin-based extended finite element method is employed. The extended finite element model of the problem is obtained by discretizing the solution domain into a number of arbitrary elements. The formulation of the XFEM for displacement components can be written as [18] n n n n n n all nodes n N cr where cr N is the set of nodes that the discontinuity has in its influence domain, while tip N is the set of nodes inside a predefined area around the crack tip ( Fig. 1). H r z , is Heaviside eichment function and m F represents crack tip enrichment functions [19]. r and are the usual crack-tip polar coordinates. Also, n a , n b and nm c are vectors of the nodal unknowns. Displacement components for the element (e) with n nodal points are approximated by compact forms as follows h = 1, 2, …, ne and m = 1, 2, 3, 4, where ne is the number of nodes in an element. Also and exhibit the enriched parts of displacements related to crack path and crack tip enrichment, respectively [20]. Applying the weighted residual integral to the equations of motion with respect to the weighting functions l S r z , , the formal Galerkin approximations reduce to where ns is the number of shape functions of the element e and l S is the component of the vector S .
For axisymmetric problems in cylindrical coordinates relations between the stresses and displacements can be expressed in the form below [21].
In this equation M and K are the mass and stiffness matrices, respectively. Also and F are the nodal displacements and force vectors, respectively. Generally, for the fictional element e which is enriched with both Heaviside and crack tip enrichment functions, these matrices and vectors can be written as follows: Matrices S 1 and S 2 are derived as follows: The material stiffness matrix D for axisymmetric problems is defined as below [22].
We use the Newmark method to solve equations of motion. The Newmark family is the most widely used family of direct methods for solving the equation of motion which consists of the following equations [23]. The Newmark family includes many widely used methods. The average acceleration method is one of them for structural dynamics applications, which is unconditionally stable. In this method, and are equal to 0.5 and 0.25, respectively. We choose the mean acceleration method in this study.

INTERACTION INTEGRAL AND SIF COMPUTATIONS
he general form of domain integral for axisymmetric problems, introduced by Moran and Shih [24] is as follows: where, c r is the crack tip radial coordinate, r q r z e , q , lj lj ij i l P W u , is the energy-momentum tensor and W is summation of strain and kinetic energy. In this work, the interaction integral method is used to compute the mode I stress intensity factor (K I ). By superimposing the actual and auxiliary fields on the domain integral, in the absence of thermal strains the general axisymmetric form of interaction integral ( MI ) in local Cartesian coordinate systems on crack tip (Fig.  2) for FGMs is obtained as below.  where q is a weight function varying from unity at the crack tip to zero on boundary of domain A * . For a stationary crack in axisymmetric state, the relation between the M-integral and the SIF is identical to its relation in plane strain state [25].

VALIDATION OF THE XFEM AXISYMMETRIC MODEL
n order to validate the results of the XFEM code, four problems were solved and their results were compared with reported values in previous researches.
In the first problem, a homogeneous solid cylinder with a central penny-shaped crack under uniform tensile load of 10 MPa, was modeled using the XFE model described in this work. The radius and length of cylinder are 50 mm and 100 mm, respectively. The central crack radius is 5 mm. The obtained SIF solution for this crack geometry is compared with the analytical and numerical solutions given by Eshraghi and Soltani [12] in Tab. 1. In the second example, a penny-shaped crack embedded in a homogeneous steel cylinder was considered. The ratio of crack radius to cylinder radius and the ratio of crack radius to cylinder height were selected to 0.2 and 0.1, respectively. The cylinder was under a uniform tensile stress of 1 MPa. The calculated SIF (K I ) is compared with the computational and analytical stress intensity factors reported by Tran and Geniaut [7], in Tab. 1. Previous problem under the loading condition of a uniform tensile stress of 1 MPa and rotation of 150 RPM was considered as the third example. In the fourth example, a complete circumferential surface crack at the inner wall of a hollow cylinder under constant axial tension of 105 MPa was studied. The values of the inner and outer radius of the cylinder were chosen to 50 mm and 55 mm, respectively with a crack length of 2.5 mm. The material was assumed to be linear elastic with values of the Young's modulus and the Poisson's ratio of E = 200 MPa and = 0.3. The computed SIF (K I ) is compared with the stress intensity factor reported by Grebner and Ustrathmeier [26], in Tab. 1.

Problem
No. It seems that the obtained results based on the XFEM method agree very well with that reported in the literature.

Penny-shaped crack embedded in a FG cylinder under static loading
n this example, an epoxy/glass functionally graded cylinder with a central penny-shaped crack (Fig. 3) Table 2: The properties of the epoxy and glass [27] The obtained stress intensity factors KI for different values of crack radius are plotted in Fig. 4. As expected, it is shown that increasing the crack radius increases the stress intensity factor, since the stresses near the crack tip increase with increasing the crack radius (Fig. 5).
The deformed mesh and the von Mises stress contours for an axisymmetric cross section are illustrated in Fig. 5. Note that a 7000 times larger scale for displacements was used to plot Fig 5(a).

Penny-shaped crack embedded in a FG cylinder under dynamic loading
An epoxy/glass functionally graded cylinder with the same dimensions as the cylinder in the previous section (Fig. 3) including a central penny-shaped crack with various radiuses under uniform impulsive tensile stress of 1 MPa is analyzed in this section. The time step and the material gradient parameter are chosen to t 1 s and P 0.2 , respectively. The impulsive tensile stress at time t 1 s is applied to the upper surface of the cylinder and the stress intensity factors are obtained at given time steps. The curves of SIF versus time for different crack radius are plotted in Fig. 6. Before reaching the stress wave to the crack tip, the SIF is zero and as the wave approaches to the crack tip, the SIF dramatically increases. After passing the wave through the crack tip, the stress intensity factor reduces. Wave arrival time to the crack tip decreases with increasing crack radius.

FG hollow cylinder with internal circumferential crack under static loading
In this section, first an epoxy/glass functionally graded hollow cylinder with an inner circumferential crack with various radiuses (Fig. 7)  The influence of the change in the material gradient parameter P on the value of stress intensity factor KI is considered. The values of SIF for different crack length are plotted (in Fig. 8) versus P . For P 0 the constituent material of the cylinder is pure epoxy. Fig. 8 shows that for large cracks the stress intensity factor decreases with increasing P , while for cracks which their lengths are less than half of the wall thickness of the cylinder, the SIF curves are downward for P 0.5 approximately, and for larger values of P the SIF curves are upward. To study the effect of the rotational speed of cylinder on the SIF of circumferential crack, cylinders with different circumferential crack length under various rotational speed are considered. For this purpose, the stress intensity factors are calculated for different values of rotational speed of cylinders, which vary from 100 to 1900 rpm. In these analysis the value of material gradient parameter is selected to P 0.2 . The SIF vs. rotational speed of cylinder curves for different values of crack length are illustrated in Fig. 9. As expected, the stress intensity factor increases with increasing rotational speed because the radial component of body force is proportional to the square of rotational speed of cylinder.  For investigating the effect of the inner radius of the cylinder on the SIF, hollow cylinders with three different circumferential crack lengths under static uniform tensile load of 1 MPa are analyzed. In this analysis, the values of the height and thickness of the cylinders are chosen to 0.2 m and 0.1 m, respectively. In addition, the material gradient parameter is considered to be 0.2. The results are presented in Fig. 11. It is observed that K I increases with increasing the inner radius of hollow cylinder. For small crack sizes, the curve rises gradually. These behaviors have been reported in reference [12].

FG hollow cylinder with internal circumferential crack under dynamic loading
According to the previous section, an epoxy/glass functionally graded hollow cylinder with an inner circumferential crack (Fig. 7) under impulsive tensile load is analyzed. The values of the height, inner and outer radius of the cylinder are chosen to H 0.2 m, i R 0.1 m and o R 0.2 m, respectively. In addition, the material gradient parameter is chosen to P 0.2 for this example. A uniform impulsive tension 1 MPa at time t 1 s is applied to the upper and lower surfaces of cylinder and the stress intensity factors are obtained at specified time steps ( t s 1 ). The SIF vs. time curves for different crack sizes are plotted in Fig.12. The stress wave arrival time to the circumferential crack tip decreases with increasing crack radius similar to a pennyshaped crack. Fig. 12 shows that the maximum value of SIF in dynamic loading condition is so much larger than its value at static loading condition. The von Mises stress and Longitudinal displacement contours for an axisymmetric cross section of a hollow cylinder with a circumferential crack of length 0.05 m at time t s 26 is illustrated in Fig. 13  To investigate the effect of the material gradient parameter P on the SIF of circumferential crack in dynamic loading condition, some functionally graded hollow cylinders composed of epoxy/glass with different values of parameter P under impulsive tensile load are analyzed. The impulsive tension 1 MPa at time t 1 s is applied to the upper surface of cylinders and the stress intensity factors are obtained at specified time steps ( t s 1 ). SIF curves in term of time are plotted in Fig. 14. It is concluded from Fig. 14 that increasing the material gradient parameter P decreases the maximum SIF value so long as P 1 while for P 1 the maximum SIF value increases with increasing P . In addition, the speed of stress wave decreases with increasing P .

CONCLUSIONS
n this study, equations of motion in cylindrical coordinate system were solved using the XFE and Newmark's methods in functionally graded materials. The stress intensity factors for penny-shaped and circumferential cracks in the epoxy/glass circular cylinders under static and dynamic loading conditions were computed using the interaction integral method. Finally, the following results were obtained from this research: 1-Increasing the crack radius increases the stress intensity factor. 2-Wave arrival time to the crack tip decreases with increasing crack radius. 3-For hollow cylinders with large circumferential crack sizes ( o i a R R / 0 . 5 ), the stress intensity factor decreases with increasing the material gradient parameter P (increasing glass inclusion). 4-For hollow cylinders with circumferential cracks which their lengths are less than half of the wall thickness of the cylinder, the SIF curves are downward for P 0.5 approximately, and for larger values of P the SIF curves are upward. 5-Increasing the rotational speed of the cylinder increases the stress intensity factor. 6-The SIF increases with increasing the inner radius of hollow cylinder. However, for small crack sizes, the effect of inner radius on SIF value is small. 7. For epoxy/glass hollow cylinders with circumferential cracks under impulsive tension, increasing the material gradient parameter P decreases the maximum SIF value so long as P 1 while for P 1 the maximum SIF value increases with increasing P . In addition, the speed of stress wave decreases with increasing P .