Two-dimensional Numerical Estimation of Stress Intensity Factors and Crack Propagation in Linear Elastic Analysis

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


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 Jcontour integration, virtual crack extension and the strain energy release rate method.The main disadvantage of these methods is that the SIF components, K I and K II in mixed mode application are either impossible or very difficult to be separated.Nevertheless, the first group which is based on neartip 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: 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:

www.etasr.com Boulenouar et al.: Two-dimensional umerical Estimation of Stress Intensity Factors and …
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 mixedmode loading, the crack will propagate in a direction normal to the maximum tangential tensile stress.In polar coordinates, the tangential stress is given by: The direction normal to the maximum tangential stress can be obtained by solving 0 = θ σ θ d d for θ.The nontrivial solution is given by: which can be solved as:

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.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].

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][11][12][13][14]. Figure 7 illustrates the predicted values of the SIFs for mode I loading during the crack propagation steps for our geometric model.

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 E matrix is the Young modulus associated with the matrix, and E incl the one associated with the inclusion, we define R as the ratio: 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 (E matrix /E incl =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 E incl =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 E incl =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].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 K I .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 K II are negligible i.e. the crack propagates towards the horizontal path under mode I loading condition.For a negative angle (hard inclusion case), K II takes negative values and for the soft inclusion case takes positive values.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.

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 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.

Figure 4 Fig. 2 .Fig. 3 .Fig. 4 .
Figure4shows the calculated values of stress intensity factors K I and K II during crack propagation steps for a cracked plate with an initial crack length of a=1 unit.The SIF K I is increase nonlinearly as crack increment length increases and the values of factor K II are negligible i.e. the crack propagates towards the horizontal path under mode I loading condition.

Fig. 5 .Fig. 6 .Fig. 7 .
Fig. 5.Problem statement and the final mesh for of the single edge cracked plate with an off-center hole