Microsoft Word - abd-ed-f.doc ETASR - Engineering, Technology & Applied Science Research Vol. 3, �o. 5, 2013, 506-510 506 www.etasr.com Boulenouar et al.: Two-dimensional �umerical Estimation of Stress Intensity Factors and … Two-dimensional Numerical Estimation of Stress Intensity Factors and Crack Propagation in Linear Elastic Analysis A. Boulenouar Mechanical Engineering Department, Djilali Liabes University of Sidi-Bel- Abbes, Algeria aek_boulenouar@yahoo.fr N. Benseddiq Mechanics Laboratory of Lille Ecole Polytech’Lille University of Lille1, noureddine.benseddiq@univ-lille1.fr M. Mazari Mechanical Engineering Department, Djilali Liabes University of Sidi-Bel- Abbes, Algeria Mazari_m@yahoo.fr Abstract—When the loading or the geometry of a structure is not symmetrical about the crack axis, rupture occurs in mixed mode loading and the crack does not propagate in a straight line. It is then necessary to use kinking criteria to determine the new direction of crack propagation. The aim of this work is to present a numerical modeling of crack propagation under mixed mode loading conditions. This work is based on the implementation of the displacement extrapolation method in a FE code and the strain energy density theory in a finite element code. At each crack increment length, the kinking angle is evaluated as a function of stress intensity factors. In this paper, we analyzed the mechanical behavior of inclined cracks by evaluating the stress intensity factors. Then, we presented the examples of crack propagation in structures containing inclusions and cavities. Keywords-stress intensity factor; crack propagation; mixed mode; inclusion I. INTRODUCTION The use of crack propagation laws based on Stress Intensity Factors (SIFs) is the most successful engineering application of fracture mechanics. This characterizes the SIFs as the most important parameters in fracture analysis. In elastic fracture analysis, SIFs sufficiently define the stress field close to the crack tip and provide fundamental information on how the crack is going to propagate. Basically, the estimation methods can be categorized into two groups, those based on field extrapolation near the crack tip and those which use the energy release when the crack propagates. The latter group includes J- contour integration, virtual crack extension and the strain energy release rate method. The main disadvantage of these methods is that the SIF components, KI and KII in mixed mode application are either impossible or very difficult to be separated. Nevertheless, the first group which is based on near- tip field fitting procedures requires finer meshes to produce a good numerical representation of crack-tip fields. Usually, the singular point elements are generated to facilitate the calculation [1]. One of the simplest and most frequently used methods is the displacement extrapolation method. It functions typically in obtaining displacement jumps along the crack faces and then applying the elasticity relations to compute a set of estimated SIF values [2]. In order to predict the fracture direction and loading based on the concept of Maximum Potential Energy Release Rate (MPERR), a new general mixed-mode brittle facture criterion was applied by Chang et al. [3]. The calculation and comparison of SIFs for a cracked Element Free Galerkin Method (EFGM) plate have been obtained by using several different numerical techniques in [4]. This paper aims to determine the SIFs for the crack propagation problem under linear-elastic fracture analysis by means of the displacement extrapolation method implemented in the Ansys finite element program. The method developed in this paper was applied to compute the stress intensity factors in elastic-plastic crack growth in plane stress problems. II. STRESS INTENSITY FACTOR AND CRACK PROPAGATION In linear elastic fracture mechanics the important parameters used are the stress intensity factors in various modes. Several methods have been proposed to determine the stress intensity factors, such as the displacement extrapolation near the crack tip [5], the J-integral [6] and the energy domain integral [7]. In this paper, the displacement extrapolation method [8] is used to calculate the stress intensity factors as follows: ( )( ) ( ) ( ) , 2 4 2 113     − −− ++ = ec dbI vv vv Lk E K π ν (1) ( )( ) ( ) ( ) , 2 4 2 113     − −− ++ = ec dbII uu uu Lk E K π ν (2) where E is the modulus of elasticity, ν is the Poisson’s ratio, L is the element length, u and v are the displacement components in the x and y directions, respectively and k is an elastic parameter defined by: ETASR - Engineering, Technology & Applied Science Research Vol. 3, �o. 5, 2013, 506-510 507 www.etasr.com Boulenouar et al.: Two-dimensional �umerical Estimation of Stress Intensity Factors and … 3 stress plane 1 . 3-4 strain plane k ν ν ν −  +  =     The near tip nodal displacements at nodes b, c, d and e, shown in Figure 1 are of a great interest. The tangential and normal displacements to crack plane are denoted as u and v respectively. In order to simulate crack propagation under linear elastic condition, the crack path direction must be determined. There are several methods used to predict the direction of crack trajectory such as the maximum circumferential stress theory, the maximum energy release rate theory and the minimum strain energy density theory. The maximum circumferential stress theory asserts that, for isotropic materials under mixed- mode loading, the crack will propagate in a direction normal to the maximum tangential tensile stress. In polar coordinates, the tangential stress is given by:       −= θ θθ π σθ sin 2 3 2 cos 2 cos 2 1 2 III KK L (3) The direction normal to the maximum tangential stress can be obtained by solving 0=θσθ dd for θ. The nontrivial solution is given by: 0)1cos3(sin =−+ θθ III KK (4) which can be solved as:         +      ±= 8 4 1 4 1 arctan2 2 II I II I K K K K θ (5) III. NUMERICAL ANALYSIS AND VALIDATION A. Single edge cracked plate under plane stress condition The geometry of the single edge cracked plate under plane stress condition and its final mesh in the first step before crack propagation are shown in Figure 2. The plate has an initial crack length of a=0.4 units, plate length of L=2 units, plate width of W=1 unit and thickness of t=1 unit. The modulus of elasticity E is 1, and Poisson’s ratio ν is 0.3. The far-field tensile stress is σ=1 unit. Fig. 1. Quarter-point triangular elements around the crack tip Figure 3 shows the final configuration corresponding to the last evaluated crack length obtained from the finite element software Ansys in the present study, which are found to be almost the same as the results obtained from the FRANC2D/L program and the results given by Alshoaibi et al. using the adaptive mesh strategy [9]. Figure 4 shows the calculated values of stress intensity factors KI and KII during crack propagation steps for a cracked plate with an initial crack length of a=1 unit. The SIF KI is increase nonlinearly as crack increment length increases and the values of factor KII are negligible i.e. the crack propagates towards the horizontal path under mode I loading condition. (a) (b) Fig. 2. (a) Problem statement and (b) the final mesh of the initial crack before crack propagates Fig. 3. Final configuration corresponding to the last evaluated crack length Fig. 4. Stress intensity factors KI and KII as functions of crack extension in a cracked plate ETASR - Engineering, Technology & Applied Science Research Vol. 3, �o. 5, 2013, 506-510 508 www.etasr.com Boulenouar et al.: Two-dimensional �umerical Estimation of Stress Intensity Factors and … B. Single edge cracked plate with an off-center hole The following example concerns a cracked part with a hole of radius r=0.2 units. In this section, we show the effect of the existence of a hole on the path of crack propagation under mode I loading conditions. This plate is meshed by an 8 nodes quadratic element and by special element meshing at the crack tip (Figure 5). The computation of the stress intensity factor is determined using the displacement extrapolation technique. At each increment of crack propagation the fork angle is evaluated using the maximum circumferential stress theory. The crack increment length ∆a is 8% of the initial crack length a. Figure 6 shows the crack propagation trajectories for our geometric model. This figure shows that the crack is moving towards the hole since it creates a stress drop. Once the crack tip has moved beyond the hole, the crack reorients horizontally in mode I loading. We obtain a similar crack path as in [10-14]. Figure 7 illustrates the predicted values of the SIFs for mode I loading during the crack propagation steps for our geometric model. Fig. 5. Problem statement and the final mesh for of the single edge cracked plate with an off-center hole Fig. 6. Crack trajectory in a cracked plate with a hole Fig. 7. Stress intensity factors KI and KII as functions of crack extension in a cracked plate with an off-center hole C. Single edge cracked plate with an off-center inclusion In the present example, we replace the hole by an inclusion of a circular form and we study the influence of this inclusion on the crack path direction. A rectangular part is pre-cracked and submitted to a tensile test. This part contains an inclusion which may be more rigid or less rigid than the matrix. If Ematrix is the Young modulus associated with the matrix, and Eincl the one associated with the inclusion, we define R as the ratio: R=Ematrix/Eincl. This example shows the ability of the discrete approach, implemented in the finite element software, to deal with multimaterial applications. Figure 8 shows the geometrical of the cracked plate and the typical model mesh at the inclusion and at the crack. In Figure 9, the crack propagation simulation is presented with a ratio R unit (Ematrix/Eincl=1), which means that we have an inclusion, but is made of the same material as the matrix. This figure shows five steps of crack propagation trajectories. As expected, the crack propagates horizontally, as in a homogeneous single-material. Figure 10 shows five steps of crack propagation trajectories for the structure with a soft inclusion (R=10 and Eincl=0.1), the inclusion is less rigid than the matrix; the crack is still attracted by the inclusion. The crack reorientation is however less pronounced than the one obtained with a hole. Conversely, if the inclusion is more rigid than the matrix (R=0.1 and Eincl=10), the crack is moving away from the inclusion (Figure 11). Good agreement is obtained for the crack path with the results presented in [15]. Fig. 8. Geometrical model and typical mesh model near the inclusion and the crack Fig. 9. Propagation of a crack in a part with an off-center inclusion (R=1) Fig. 10. Propagation of a crack in a part containing a soft inclusion (R=10) ETASR - Engineering, Technology & Applied Science Research Vol. 3, �o. 5, 2013, 506-510 509 www.etasr.com Boulenouar et al.: Two-dimensional �umerical Estimation of Stress Intensity Factors and … Fig. 11. Propagation of a crack in a part containing a hard inclusion R=0.1 The crack growth and its direction are determined by the calculated SIFs. Figure 12 shows the predicted values of the SIFs for mixed mode during crack propagation steps. As shown, the increase in crack length causes an increase of KI. This factor has the same curves for the three cases (R=1, 0.1 and 10), i.e. it does not depend on the crack path direction. In the case of R=1, the values of KII are negligible i.e. the crack propagates towards the horizontal path under mode I loading condition. For a negative angle (hard inclusion case), KII takes negative values and for the soft inclusion case takes positive values. Fig. 12. Stress intensity factors KI and KII as functions of crack extension in cracked plate with : (a) R=1 (b) R=0.1 (c) R=10 Figure 13 illustrates a comparison between the paths of propagation, for the various applications presented in this study. Figure 14 shows the normal stress distribution for the final step chosen of the crack propagation, for four cases presented in this study. The maximum value of the normal stress is localized around the perturbated zone by the existence of a crack. Fig. 13. Crack trajectories comparison Fig. 14. Contour plots of the normal stress for cracked plate with (a) an off- center hole (b) R=10 (c) R=1 (d) R=0.1 IV. CONCLUSION The predicted values of stress intensity factors for the crack propagation problem under linear elastic fracture analysis using the finite element method (FEM) and the displacement extrapolation technique was presented. In order to obtain a better approximation of the stress field near the crack tip, special quarter point finite elements as proposed by Barsoum are used. The estimated SIFs are used to further approximate the crack increment length and to predict the crack path ETASR - Engineering, Technology & Applied Science Research Vol. 3, �o. 5, 2013, 506-510 510 www.etasr.com Boulenouar et al.: Two-dimensional �umerical Estimation of Stress Intensity Factors and … direction. The prediction of the latter is based on the maximum circumferential stress theory. Finally, the numerical finite element analysis for 2D with displacement extrapolation method and maximum circumferential stress theory, have been successfully employed to determine the effect of the inclusion on the crack path direction under linear-elastic analysis. REFERENCES [1] T. Denyse De Araújo, T. N. Bittencourt, D. Roehl, L. F. Martha, “Numerical Estimation of Fracture Parameters in Elastic and Elastic- plastic Analysis”, ECCOMAS 2000, European Congress on Computational Methods in Applied Sciences and Engineering, Barcelona, September, 2000. [2] A. B. de Morais, “Calculation of stress intensity factors by the force method”, Engineering Fracture Mechanics, Vol. 74, No. 5, pp. 739-750, 2007 [3] J. Chang, J. Xu, Y. Mutoh, “A general mixed-mode brittle fracture criterion for cracked materials”, Engineering Fracture Mechanics, Vol. 73, No. 9, pp. 1249-1263, 2006 [4] G. Anlas, M. H. Santare, J. Lambros, “Numerical calculation of stress intensity factors in functionally graded materials”, International Journal of Fracture, Vol. 104, No. 2, pp. 131-143, 2000 [5] S. K. Chan, I. S. Tuba, W. K. Wilson, “On the finite element method in linear fracture mechanics”, Engineering Fracture Mechanic, Vol. 2, No. 1, pp.1-17, 1970 [6] D. M. Parks, “A stiffness derivative finite element technique for determination of crack tip stress intensity factors”, International Journal of Fracture, Vol. 10, No. 4, pp. 487-502, 1974 [7] B. Moran, C. F. Shih, “A general treatment of crack tip contour integrals”, International Journal of Fracture, Vol. 35, No. 4, pp. 295-310, 1987 [8] S. Phongthanapanich, P. Dechaumphai, “Adaptive Delaunay triangulation with object-oriented programming for crack propagation analysis”, Finite Element in Analysis and Design, Vol. 40, No. 13-14, pp.1753-1771, 2004 [9] A. M. Alshoaibi, A. K. Ariffin, “Finite element simulation of stress intensity factors in elastic-plastic crack growth”, Journal of Zhejiang University SCIENCE A, Vol. 7, No. 8, pp. 1336-1342, 2006 [10] P. O. Bouchard, F. Bay, Y. Chastel, I. Tovena, , “Crack propagation modelling using an advanced remeshing technique”, Computer Methods in Applied Mechanics and Engineering, Vol. 189, No. 3, pp. 723–742, 2000 [11] M. Souiyah, A. Alshoaibi, A. Muchtar, A. K. Ariffin, “Finite element model for linear-elastic mixed mode loading using adaptive mesh strategy”, Journal of Zhejiang University SCIENCE A, Vol. 9, No. 1, pp. 32-37, 2008 [12] T. N. Bittencourt, P. A. Wawrzynek, A. R.Ingraffea, J. L. Sousa, “Quasi- automatic simulation of crack propagation for 2D LEFM problems”, Engineering Fracture Mechanics, Vol. 55, No. 2, pp. 321-334, 1996 [13] M. M. Rashid, “The arbitrary local mesh replacement method: an alternative to remeshing for crack propagation analysis”, Computer Methods in Applied Mechanics and Engineering, Vol. 154, No. 1-2, pp. 133-150, 1998 [14] D. Azócar, M. Elgueta, M. C. Rivara, “Automatic LEFM crack propagation method based on local Lepp–Delaunay mesh refinement”, Advances in Engineering Software, Vol. 41, No. 2, pp. 111–119, 2010 [15] P. O. Bouchard, F. Bay, Y. Chastel, “Numerical modelling of crack propagation: automatic remeshing and comparison of different criteria”, Computer Methods in Applied Mechanics and Engineering, Vol. 192, No. 35-36, pp. 3887–3908, 2003