Analysis of stress intensity factors for functionally graded cylinders with multiple longitudinal cracks using finite element method

In this paper, the cracked thick-walled functionally graded cylinder has been analyzed using the finite element method. The internally pressurized cylinder contains fully longitudinal cracks. The analyses have been done with two, four, six and eight longitudinal cracks at the inner surface of the cylinder in four different conditions. For this purpose, using the USDFLD subroutine coding in ABAQUS software, the variations of properties of the functionally graded material are considered based on a power-law function model in the cracked cylinder. The inner and outer surfaces of the cylinder are made of aluminium oxide and titanium carbide respectively. The continuous composition of the two materials has been considered in the form of varying elastic modulus and Poisson’s ratio along the radial direction as the power-law function. The J-integral has been used to calculate the stress intensity factors, taking into account the variable properties at the crack tip. The effect of non-homogeneity of the properties, the number of the cracks and cracks’ relative depth on the stress intensity factors have been evaluated. The results have been compared with those of other available literature to verify the finite element modelling, in which very good agreement has been found. The results show that similar to the isotropic cylinder with multiple longitudinal cracks, in the cylinder of functionally graded material, the two-crack model has the highest amount of stress intensity factors and the variation of Poisson’s ratio has to be taken into account in the shallow cracks. c © 2019 University of West Bohemia. All rights reserved.


Introduction
Functionally graded materials (FGM) are component of advanced engineering material consisting of two different materials. In this material, properties change as a continuous function from one surface to the other one. In multilayer composites, varying the mechanical and thermal properties of the layers creates a stress concentration at the interfaces which leads to debonding between fibers and matrix or between laminas. The use of FGMs eliminates the creation of stress concentration sources along the interfaces and reduce the mechanical mismatch. The two materials' composition nature is so that the material properties change from one side of the thickness to the other one continuously. Typically, a combination of ceramic and metal is used in these materials to utilize the thermal resistance of ceramics on one side and the flexibility and toughness of the metal on the other side. Various industries such as aerospace and biomedicine focus on the FGM due to its special feature. Pressure vessels and pipes containing corrosive fluid are more prone to cracking. Consequently, examining crack's behaviour is highly critical to estimate its service life and measure its safety and replacing these tubes with an FGM type can prolong life. Amid this, risk analysis can significantly help define the requirements for selecting the material through the probability crack creation analysis approach. Accordingly, the analysis and determination of stress intensity factors (SIFs) is an effective step in this process.
A considerable amount of literature has been published on the isotropic material containing fully longitudinal cracks with diverse methods in terms of its use and loading, each of which has its own merits and demerits. Andrasic and Parker [6], Perl and Arone [17], Kirkhope et al. [13], Shu et al. [21] determined the SIFs for longitudinal cracks using the finite element method (FEM). All analyses have been modelled two-dimensionally in accordance with the crack type. Cylinder's geometric dimensions along with the number of cracks have been the criterion for doing these researches. It is worth noting that among the number of the existing cracks, in all papers, the cylinder with two cracks has achieved the highest SIF. In other words, the crack number increase does not mean the structure being weak. As a result, the analysis with two cracks is more likely to be found in other articles. The analysis of the isotropic cylinder with one and two longitudinal cracks was carried out using the weight function method, respectively, by Varfolomeyev and Hodulak [23] and Nabavi et al. [16]. Presenting the crack weight function (WF) allows using for different loads. Gonzalez and Martinez [11] analyzed one longitudinal crack using the boundary element method. Applying standard codes in the study by Miura et al. [15] has been reviewed in the cylinder containing different kinds of crack.
The last two decades have seen a growing trend towards analyzing SIFs in the cracked cylinders made of an FGM compared to the isotropic ones. Afsar et al. [2][3][4][5] determined the SIFs in an FGM cylinder containing two longitudinal cracks using the FEM. The effect of properties and geometries of the cylinder on the SIFs has investigated. Mirahmadi et al. [14] also obtained the WF of two longitudinal cracks in an FGM cylinder. The FGM model has considered as an exponential function. The unknown parameters in the WF meant that it was not possible to employ the weight function in a practical situation. By solving a singular integral equation, Ç etin and Kadıoglu [7] analyzed the circumferential embedded crack in an FG cylinder. Exponential function employed to describe the variation of properties in the study. Shariati and Rokhi [20] have used the extended FEM to analyse the internal crack in an FG cylinder by a power-law function model. In the continuation, Elahi and Rokhi [9] analyzed the thermal loading induced SIFs in this problem. Eshraqi and Soltani [10] determined the WF for an internal crack in an FG cylinder with a power-law model using reference loadings. Zhou [24] dealt with implementing the crack problem in an FGM using ABAQUS software. In this paper, an implementation method with virtual crack closure technique for the analysis of fracture problems in non-homogeneous materials like the FGM has been addressed using UMAT and UEL subroutine in the software. Seifi and Dejam [19] investigated the variations of the SIF and energy release rate based on the displacement correlation and modified crack closure integral methods for external semi-elliptical cracks in an autofrettaged FG cylinder. Seifi [18] also solved this problem using the WFM.
As seen, most of the activities performed on the longitudinal crack area in the FG cylinder are subject to two cracks. No previous study has investigated the multiple cracks in an FG cylinder. Besides, in defining the properties, it is limited to the elastic modulus and Poisson's ratio has been assumed constant. In this paper, the pressurized cylinder made of the FGM containing multiple longitudinal cracks is analyzed using the FEM. The number of the cracks has been considered in four various modes as 2, 4, 6 and 8 fully longitudinal cracks and has been analysed separately. The four-crack mode has presented in Fig. 1 schematically. Choosing these four modes has been achieved with the results evaluation approach in the FG cylinder mode and studying its behaviour compared with that of the isotropic cylinder mode. It has been assumed that the inner and outer surface of the cylinder is made of aluminium oxide and titanium carbide, respectively. The FGM's properties have been considered as the changing elastic modulus and also Poisson's ratio. The power-law model has included in the form of varying properties of aluminium oxide in the cylinder's inner surface to titanium carbide's properties in its external surface. USDFLD subroutine coding in ABAQUS software [1] has been employed to implement material nonhomogeneity. The J-integral has been used to determine the stress intensity factors in the cracked FG cylinder. In the isotropic material, the J-integral is path independent, while in the FGM, the J-integral related to the integral path. This independence has to be provided due to the continuous change of material properties and nonhomogeneity [8]. The dependency is suppressed by the use of the sufficiently fine mesh around the crack tip. The effect of nonhomogeneity ratio, the crack number and the crack geometries on SIFs have considered as an important factor in identifying the crack behaviour in the FGM. The results have been compared with those of other literature in special cases to validate the analysis. These results agree with the findings of other studies well and indicate the possibility to use the above-mentioned subroutine coding.

Distribution of properties in functionally graded material
Since functionally graded materials have their own features, then it is imperative to consider a model for the properties' distribution. The nature of the FGM is defined by the change in the volume fraction of its two structural materials in which three models of properties distribution, namely, power-law function, sigmoid distribution and exponential function are employed more than the other ones. In each of these models, the properties have been considered in the general state of a non-linear combination in the volume fraction of its two materials. Although there is no restriction to simulate the properties for the FG cylinder in the USDFLD subroutine code, the power-law distribution, which is applicable to any volume fraction, is assumed to analyse the cracked FG cylinder. The power-law distribution can be used to construct a wide variety of distributions. The modulus of elasticity and the Poisson's ratio varies continuously along the wall-thickness of the cylinder as the following: where n stands for the volume fraction index of the FGM and characterizes the manner combining the volume fraction of the material used in the internal surface (aluminium oxide) and the external surface (titanium carbide) for the distribution of the elastic modulus and the Poisson's ratio according to Figs. 2 and 3 along the thickness, respectively. The magnitude of the volume fraction index determines the curvature of E(r) and ν(r). The curvature can be made concave upward for n > 1 and concave downward for 1 > n > 0. As the exponent becomes much bigger than one, the gradient tends to be smoother near aluminium oxide. Adjusting the volume fraction index will control the rate of transition between aluminium oxide and titanium carbide.
Depending on the type of production of the FGM, the proportion of aluminium oxide is greater than that of titanium carbide.

The J-integral method for determining stress intensity factors in functionally graded material
The nature of the J-integral in the linear elastic fracture mechanics is the energy release rate defines as the difference in the potential energy value of a cracked structure in two different crack lengths as follows: in which Γ is the arbitrary closed contour that has to start from a point on the lower crack face and ends on its upper one. Additionally, W is the strain energy density, n i is the unit vector component perpendicular to the closed path, σ ij and u j are the stress tensor and displacement field, respectively. In practice, the values within this integral are calculated by solving the FEM and then, J is calculated through numerical integration. The advantage of using this method is related to elastic and plastic analysis. The SIFs for homogeneous material are obtained by determination of the J-integral amount as follows: The J-integral can be determined at high accuracy by selecting suitably the closed contours away from the crack tip in homogeneous material. In contrast to isotropic material, the J-integral is path dependent in the FGM generally. There are two approaches to overcome this problem. In the first approach, the modified J-integral is used as the following [8]: where A I is the enclosing the area within the arbitrary closed contour Γ and q is the weight function. In this state, by adding the second term to the integral, the non-homogeneous effect of the FGM is taken into account. In the second approach, based on Gu et al. [12], by choosing very fine mesh around the crack tip, the second term of the modified J-integral containing the energy derivative can be neglected. Moreover, this can make it possible to assume constant properties inside elements near the crack tip in FGMs. The SIF can be determined based on this approach as the following: where E tip and ν tip are the modulus of elasticity and Poisson's ratio of the FGM in the position of the crack tip. In the subroutine used in ABAQUS software, it is possible to define different properties with this method in the FGM and there is no obstacle in determining the SIFs.

The finite element modeling
The ABAQUS finite element software [1] has been used to model an FGM cylinder. It is possible to define the variation of the physical properties of the FGM using FORTRAN coding in the USDFLD subroutine of the software. Besides, regarding the study cracks being of the fully longitudinal type, 2-D modelling has been used. In term of the existing crack number, it is viable to divide the cylinder into equal parts by applying the boundary conditions of symmetry over the equivalent border. In the two and four crack modes, the two planes of symmetry have been employed as shown in Fig. 4 for the latter one. The FG cylinder was meshed by 2D plane strain elements with reduced integration. At the crack tip, six-node modified quadratic singular elements with hourglass control (CPE6M) have employed and in other parts, eight-node biquadratic elements (CPE8R) have used for meshing the models. The standard J-integral is path-dependent when applied to the FGM directly. Instead of using the modified J-integral, in order to determine the SIFs, the approximated model of Gu et al. [12] has employed by applying the fine mesh around the crack tip. A sufficient fine mesh near the crack tip has been applied to achieve the convergence of the results. The size of elements controls the error in neglecting the nonhomogeneity effect and the the J-integral is path-independent and accrued. The size of elements around the crack tip is one percent of the crack depth in all FE models. In other words, the first ring of elements is arranged within the radius a/100 from the crack tip to suppress the J-integral path dependency. The difference among the J-integral values for 10 contours is less than 1% by applying very fine mesh around the crack tip and are accessible by determining SIF from (4).
To verify the functionally graded material's modelling, at first the simulation and analysis of the thick-walled cylinder made of the FGM subjected to internal pressure have been dealt with and the results have been evaluated and examined with the available data. Secondly, after getting assured of the modelling material simulation accuracy in the software, the crack modelling analysis is addressed. For this purpose, the pressurized homogeneous cylinder containing longitudinal cracks have analyzed and the results have compared with those of other papers. Good agreement is observed among the results. And finally, the FG cylinder with two cracks has been analysed with the combination of the two preceding models, which indicates high precision and accuracy of the modelling in the subroutine coding.

Validation of functionally graded material's modelling
The FG cylinder subjected to the internal pressure P has been considered according to the study of Tutuncu and Ozturk [22] in order to validate the FGM modelling in the subroutine coding of the software. In this reference, the Poisson's ratio is constant and the variation of modulus of elasticity have assumed as E(r) = E 0 r (for n = 1 representing the nonhomogeneous ratio of 1). The geometric dimensions and properties have been listed in Table 1.
The coding derived results have been compared with those of the analytical solution of the paper by Tutuncu and Ozturk [22], considering the nonhomogeneous ratio 1 (see Table 2). The maximum error achieved between the numerical and accurate analytical solution is less than 0.2%, which indicates very high accuracy of the modelling in the coding.

Validation of fully longitudinal crack modeling in the homogeneous cylinder modelling
In order to validate the fully longitudinal crack modelling, three modes of 1, 2 and 4 crack(s) in the internally pressurized homogeneous cylinder have been modelled and compared with the results obtained by Shu et al. [21]. The geometric dimensions and properties according to this reference, the ones listed in Table 3, have been taken into account to compare suitably. The normalized SIFs (K N = K I /P √ πa) for the homogeneous ratio of zero, indicating the material used in the code is isotropic, have been given in Table 4 and compared with the results of Shu et al. [21] in all three crack modes and the maximum difference has been gained less than 1%. Table 3. Geometric parameters and properties of the homogeneous cylinder containing longitudinal crack [21]

Validation of fully longitudinal cracks in the functionally graded cylinder modelling
Integrating the two preceding parts, it is possible to cope with modelling the cracks in the FGMs. For this reason, the functionally graded materials containing two fully longitudinal cracks have been modelled and analyzed according to the properties and dimensions of Table 5 based on the study by Mirahmadi et al. [14]. The elasticity modulus vary has been considered as an exponential function Table 6, the values of the dimensionless SIFs have been presented for the two longitudinal symmetric cracks resulted from the coding and this reference. In Table 6, the values of the normalized SIFs have been presented for the two longitudinal symmetric cracks resulted from the coding and this reference. The comparison with each other implies very high congruence and the highest difference has obtained 0.1%. Considering the used code validation in diverse above modes and comparing it with the available results in the literature, it can certainly be stated that this coding can be applied for modelling and analyzing the FG cylinder containing fully longitudinal crack.

Results and discussion
Non-dimensional SIFs in the thick-walled cylinder made of the FGM with the properties powerlaw distribution containing two, four, six and eight fully longitudinal cracks have been determined using the written code. The cylinder radius ratio, R o /R i , is 1.5 and the cylinder is subjected to internal pressure. The cylinder is long enough to neglect the edge effects on the SIFs. The inner surface of the FGM cylinder is made of Al 2 O 3 with the modulus of elasticity and Poisson's ratio 380 GPa and 0.26, respectively and the outer one is made of TiC with the ratios as 462 GPa and 0.19, respectively. All the material properties of the FGM are the function of radius only. The properties power-law distribution has been studied as a change of elasticity modulus and Poisson's ratio for the nonhomogeneity ratios as 1, 2, 4, 8 and 10. In the power-law model, as the nonhomogeneity ratio increases, the rate of properties change is slow initially and then increases. The relative crack depth, a/t, varies from 0.2 to 0.8 in the analysis.
The dimensionless values of SIFs ( Tables 7-10 have been presented for four types of cracks in the FGM cylinder. As can be seen from the tables, in a constant nonhomogeneity ratio in all four models, as the relative crack depth increases the values of the dimensionless SIFs increase. The rate of this increase declines with increment nonhomogeneity of material. However, if the material change rate from the inner to the outer surface of the cylinder increases, the SIFs in the shallow cracks increase and in the deep cracks decline. It is worth mentioning that as a result of the slight variation in the properties of the inner portion of the cylinder, due to the change in the nonhomogeneity ratio, changes in shallow cracks are also small. Therefore, in low-depth cracks, with increasing the nonhomogeneity ratio to more than eight, the SIFs approaches a constant value.     Another important practical implication is that as the crack number rises, the SIFs reduce in all different relative crack depths. This means that regarding the life estimation of the cylinders containing several cracks, the two-crack mode has the shortest working life. There are similarities between the attitudes expressed by multiple cracks in the FGM cylinder in this study and those described in the isotropic cylinder by various references [6,13,17,21]. In other words, the crack number increase inherently will cause about decreasing the SIFs in the FGM. This phenomenon is not related to material properties and is preserved as the nonhomogeneity ratio changes. As stated in other studies, for the isotropic cylinders, the two-crack mode is the critical state of multicrack mode. This phenomenon similarly reappears in the non-homogeneous FG cylinder and the two-crack mode has greater SIFs relative to four, six and eight crack modes. The effect of Poisson's ratio changes on the SIFs in the FG cylinder with linear variation of properties (n = 1) is considered and presented in Table 11. It can be seen from the data in

Conclusion
The purpose of the current study was to determine the SIFs in the FG cylinder containing multiple longitudinal crack using the finite element method. This was the first study to undertake a multiple longitudinal cracks analysis in the FG cylinder. The analysis of the FGM has been accomplished using the subroutine coding in ABAQUS software. The nonhomogeneity effect of the material composition has been considered in the form of changing elastic modulus and Poisson's ratio of the cylinder with the power-law model. Four modes have been investigated and assessed in terms of the crack number.
The results indicate that the FG cylinder with two cracks is the most critical type of multiple cracks in the assessment of the cylinder's life with the criterion of linear elastic fracture mechanics. It is worth noting that applying variation of Poisson's ratio in modeling the FGM in shallow cracks is significant.