# Finite Element Model of Crack Growth under Mixed Mode Loading

## Transcript Of Finite Element Model of Crack Growth under Mixed Mode Loading

International Journal of M aterials Engineering 2012, 2(5): 67-74 DOI: 10.5923/j.ijme.20120205.02

Finite Element Model of Crack Growth under Mixed Mode Loading

Souiyah Miloud1,*, A. Muchtar2, A. K. Ariffin2, Malek Ali1, M. I. Fadhe l1, Base m Abu Zne id3

1Faculty of Engineering & Technology (FET), M elaka, 75450, M ultimedia University (MMU) 2Faculty of Engineering & Built Environment, Bangi, 43000, University Kebangsaan M alaysia (UKM )

3Universiti Teknologi M alaysia Skudai, 81310, Johor Bahru

Abstract In this paper, in order to predict the crack growth trajectory and to evaluate the SIF under mixed modes (I & II),

one proposes a new finite element program for crack growth using the source code written in FORTRAN. The fin ite element mesh is generated using an advancing front method, where the generation of the background mesh and the construction of singular elements are also added to this developed programme to facilitate the crack process and the fracture analysis. Displacement Ext rapolation Technique (DET) was employed to evaluate the SIFs under mixed mode loading conditions. Therefore, the accuracy of both SIF`s values and the crack path predictions results are compared and validated with other relevant published research work. Ho wever, the assessment indicated that this developed fin ite element programme is reliable and robust to evaluate the SIFs and predicts the crack trajectories successfully based on the applied loading conditions.

Keywords Linear Elastic Fracture Mechanics, Adaptive Refinement Mesh, Stress Intensity Factors, Finite Elements

Method, Crack Propagation

1. Introduction

Structural failu re can be generally associated with one or more fracture of the materials making that structure. In fact, fractu re mechanics deals with the irrevers ible process of rupture due to nucleation or sudden crack and crack gro wth. The influence of cracking on the structural strength has been wid ely ap p reciated s in ce th e end o f th e 19th centu ry. Ho wever, so me aspects o f its natu re and influen ce st ill remain un kno wn[1]. Th e use o f crack p ropagat ion laws based on stress intensity factor range is the most successful engineering applicat ion of fracture mechanics. In the elastic fracture analysis, the stress intensity factors sufficient ly define the stress field close to the crack tip and prov ide fundamental in fo rmat io n o f ho w th e crack is go ing to p rop ag at e. Bas ically , th e est imat io n meth ods can be categorized into two groups , thos e bas ed on field extrapolation near the crack t ip and those which make use of the energy release when the crack propagates. The latter group includes the J-contour integration, the virtual crack extension and the strain energy release rate method. The main d isadv antag e o f these methods is th at th e st ress intens ity factor co mponents, KI and KII in mi xed mode prob lems are eit her imposs ib le o r very d ifficu lt t o be

* Corresponding author: [email protected] (Souiyah Miloud) Published online at http://journal.sapub.org/ijme Copyright © 2012 Scientific & Academic Publishing. All Rights Reserved

separated. Nevertheless, the first groups which are based on near-tip field fitting procedures require finer meshes to produce a good numerical representation of crack-tip fields generated to facilitate the calculat ion.

2. Finite Element Model

Most established research was on computer base by applying numerical methods which were mostly uses commercial software such as;[2] an enriched partit ion of unity finite element method (PUFEM ), which is known as one of the meshless method to calculate the SIF in LEFM under plane stress and plane strain condition. Another, numerical examp les were performed to evaluate the generated stress intensity factors directly fro m the scaled boundary ﬁnite element solution for the singular stress ﬁeld[3]. More recent work was developed variat ional mesh-free method, to evaluate the stress intensity factors of mixed mode crack problems, using the element free Galerkin method[4]. Accordingly, another new method was proposed to determine the SIFs for the indentation problem based on the conservation integral[5]. In contrast, Numerical experiments were conducted to evaluate the effectiveness of two proposed techniques on near crack tip singularity[6]. In fact, an observation was stated that a general theoretical model for automatically evaluating the increments of crack growth during a loading process does not exist yet[7]. Expe riments work on a central crack specimen with holes, to predict the crack path and to evaluate the mixed mode stress

68

Souiyah M iloud et al.: Finite Element M odel of Crack Growth under M ixed M ode Loading

intensity factors (KI & KII)[8]. Somehow, this limit nowadays still characterizes the

computer codes for crack growth simulat ion in elastic and elasto-plastic materials. Generally, the software should address three different aspects of the problem, specifically are; (i) SIFs evaluation, (ii) the crack propagation direction, and (iii) mesh modification to accommodate the crack advancement. The first and second aspects can be solved with no particular d ifficult ies and many comme rcial FE codes have built-in capabilit ies to evaluate the SIF. The last aspect is much mo re co mplex and is rarely solved by FE codes, e.g. FRANC2D, FRANC3D and ZENCRA CK[9]. Eventually, the FE model is modified in order to accommodate the obtained evaluation on short and straight crack propagation directions. The whole process is repeated many times and both is time-consuming and a source of errors if it is performed manually. This is part icularly t rue for the mesh modification. However, software on the other hand, is capable in handling the solution process in an automat ic way.

However, in the current study, we developed a new fin ite element programme for crack growth using the source code written in FORTRAN. In order to predict the crack growth trajectory and to evaluate the SIFs under mixed modes (I & II) loading in the frame of LEFM. The displacement extrapolation techniques with the adaptive refinement mesh method are used, to determine the stress intensity factors on two different geometries. Specifically are; four points bending plate and, centre double cracked plate with two-holes. Additionally, the obtained results of the current study are evaluated and validated with other relevant experimental and numerical results selected from the literature.

2.1. Generation of Singular Element

An unstructured triangular mesh is automatically generated by employing the advancing front method. The singular elements have to be constructed correctly to get a proper field of singularity around the crack tip as shown in Fig. 1.

Figure 1. The cut and patch procedure of generating singular elements around a crack tip

The number o f elements depends on the distributed nodes around the crack tip, wh ich can be set by the user as shown in Fig. 2.

In fact, the natural triangular quarter point elements are used[10] instead of the collapsed quadrilateral element, which is proposed by Barsoum[11]. The success of the adaptivity in general depends to a large extent on the efficient coupling between the error estimator, refinement scheme and automatic mesh generator. The importance of

these adaptive techniques in practical applications has led to a considerable research on fully automatic mesh generators that require only the specificat ion of the boundary and mesh size d istribution over the domain. Moreover, several comprehensive numerical tools have been developed to enhance the accuracy at the crack tip.

c b d

e

v,y

u,x

L

3L

4

4

Figure 2. A typical arrangement of the natural triangle quarter-point element s around a crack t ip

In this work, the unstructured triangle mesh is automatically generated by employing the advancing front method. Many researchers applied the FEM with remeshing to study the fracture propagation and its SIFs analysis ([12]; and[13]), though it is not an easy task. In contrast, it is almost impossible to automatically remesh finite elements of an arbitrarily growing crack[14]. One of the advantages of the advancing front method is that the new triangle element formation is coinstantaneous with new node generation, and this advantage makes it possible to control the shape and size of the element through the adjustment of the location of the node[15]. However, a lot o f intersection checking between the new generated triangular elements and the existing elements must be computed in order to ensure that the triangular elements are valid. Accordingly, an observation was stated that the algorithm for the advancing front method has been shown to be robust in two-dimensional mesh generation of triangles and could be extended easily to generate quadrilaterals[16].

2.2. Mesh Generati on and Adapti ve Refinement Mesh

In general, the s maller mesh size g ives more accurate fin ite element appro ximate solution. Ho wever, reduction in the mesh size leads to greater computational effort. The adaptive mesh refinement is emp loyed as the optimization scheme. This scheme bases on a posteriori error estimator which is obtained from the solution fro m the previous mesh. Basically, the success of the adaptivity in overall is depends to a large extent on the efficient coupling between the error estimator, refinement scheme and automatic mesh generator.

International Journal of M aterials Engineering 2012, 2(5): 67-74

69

The importance of these adaptive techniques in practical applications has led to a considerable research on fully automatic mesh generators that require only the specification of the boundary and mesh size distribution over the do main. An examp le of geo metry shows the whole process of generating the mesh is illustrated in Fig.3.

(a)

(b)

(c)

(d)

(ii) Determine the average error norm over the whole domain

∑ ∫ =eˆ 1 m σ Tσd Ω

(2)

m e =1 Ωe

Where m is the total number of elements in the whole do

ain. (iii)Determine a variable, εe for each element as

e 1/ 2

( ) 1 e

( ) εe = η eˆ 1/ 2

(3)

Where η is a percentage that measures the permissible

error for each element. If εe > 1 the size of the element is reduced and vice versa. (iv) The new element size is

determined as

hˆ = he

e

(4)

(εe )1/ p

Where he is the old element size and p is the order of the interpolation shape function.

2.3. The Dis placement Extrapolation Techni que

(g)

(e)

(f)

Figure 3. The mesh generation stages

The geometry of a p late with six holes and two notches is

illustrated Fig. 3a. Six connector lines as shown in Fig. 3b,

are fo rcing the internal boundaries to be the continuous part

of the external boundary. Fig. 3c shows the cutting out of the

rosette templates around each crack tip. The background

mesh for this domain is then set up automatically using

dichotomy technique as shown in Fig.3d and Fig. 3e shows

the conventional mesh being generated by the advancing

front method. The first generation produces mesh with initial

size set by user. Later, during adaptive refinement, this first

generated mesh will be taken as the background mesh. In Fig.

3f, for each rosette template, quarter-point elements are then

constructed. Fig.3g, shows the enlargement of the

quarter-point element at each crack t ip. In general, the

smaller mesh size gives more accurate fin ite element

approximate solution. However, reduction in the mesh size

leads to greater computational effort.

The strategy used to refine the mesh during analysis

process as follo ws[17]:

(i) Determine the error norm for each element

∫ ( ) ( ) e e = σ −σ * T σ −σ * d Ω

(1)

Ωe

Where σ is the stress field obtained fro m the fin ite element

calculation and σ* is the smoothed stress field inclusion of

quarter point elements.

One of the simplest and most frequently used methods is

displacement ext rapolation technique[18]. It consists

typically in the effective SIF concept by which, the fracture

evaluation can be easily carried out[19]. The crack length

was evaluated by linear-elastic analysis fro m the compliance

of single-edge-notched specimen in three-points bending

test[20]. A new ﬁn ite element of quasi-static brittle fracture

was developed, based on computational framework for

three-dimensional cracks propagation bodies. They uses

consistent thermodynamical framework for crack

propagation in elastic solids[3].

The asymptotic expression for the displacement norma l to

crack p lane, v under mode I loading[21]:

=v KI 14+Eν 2πr (2κ +1) sin θ2 − sin 32θ + A1(1E+ν )r (κ − 3) sinθ (5)

+

A2 (1+ν )r3

2

(2κ

−1)

sin

3θ

− sin θ

+ ...

E 3

2

2

Where KI is the stress intensity factor for mode I, E is the

modulus of elasticity, v is the Poisson’s ratio, K an elastic

parameter defined in equation (6), Ai are parameters

depending on the geometry and load on the specimen, and

r and θ are the polar coordinates defined at the crack t ip.

(3 − 4ν )

plane stress

κ = (3 − 4ν ) / (1+ν ) planestrain (6)

The near tip nodal displacements at nodes b, c, d and e

shown below in Fig. 3, are of interest. The displacements are

e xtrapolated by evaluating Equation (5) a long the crack faces

𝜃𝜃 = ± 𝜋𝜋. Particu larizing for nodes b and c on the singular element at the upper face of the crack yields:

= v K 2 (1+ν )(κ +1) L − A2 (1+ν )(κ +1) L3 2 + O(L5 2 ) (7)

b

Iπ

4E

12E

70

Souiyah M iloud et al.: Finite Element M odel of Crack Growth under M ixed M ode Loading

v= K 2 (1+ν )(κ +1) L − 2A2 (1+ν )(κ +1) L3 2 + O(L5 2 ) (8)

c

Iπ

2E

3E

Where L, is the length of the element side which is

connected to the crack tip. Equations (7) and (8) can be

solved to obtain the value of K I as:

= K E

2π (8v −v )

(9)

I 3(1+ν )(1+ κ ) L

b

c

The nodal displacements at the other two nodes can be evaluated by the similar way. If one considers the all four nodal displacements, the stress intensity factor expression b eco mes :

KI

E

2π

4

(v

b

−v ) − (v c d

−v e )

(10)

3(1+ν )(1+ κ ) L

2

Under pure mode II loading, the displacement tangential

to the crack plane, u is given by:

=u K 1+ν

2r

(2κ +1) sin θ

−

sin

3θ

+

A1 (1 +ν

)r

(κ

−

3) sinθ

(11)

II 4E π

2

2

E

+

A2 (1+ν )r 3

2

(2κ

−1)

sin

3θ

− sin θ

+ ...

E 3

2

2

Similarly, the stress intensity factor for mode II, using the

nodal displacements of the four nodes can be estimated by:

K II

E

2π

4

(u

b

−ud

)−

(uc

−ue )

(12)

3(1+ν )(1+ κ ) L

2

In order to simulate crack propagation under linear elastic

condition, the crack path direction must be determined.

There are several methods use to predict the d irection of

crack trajectory such as the maximu m circu mferential stress

theory, the maximu m energy release rate theory and the

minimu m strain energy density theory. The maximu m

circu mferential stress theory asserts that, for isotropic

materials under mixed-mode loading, the crack will

propagate in a direction normal to maximu m tangential

tensile stress. In polar coordinates, the tangential stress is

given by

σ

1 cos θ K cos2 θ − 3 K sin θ (13)

θ 2π r 2 I

2 2 II

The direct ion normal to the maximu m tangential stress can

be obtained by solving d σθ / d θ = 0 for θ . The nontrivial solution is given by

K I sinθ + K II (3cosθ −1) =0

(14)

Which can be solved as:

−1

3K

2 II

+KI

K

2 I

+

8K

2 II

θ0 = ± cos

2

2

(15)

K I + 9K II

y'

y'

In order to ensure that the opening stress associated with the crack direction of the crack extension is maximu m, the

sign of θ0 should be opposite to the sign of K II [19]. The

two possibilities are illustrated in Fig. 4. The criterion for crack to propagate fro m crack t ip is based

on the material toughness, KC . If the calculated stress

intensity factor, K I ≥ KC then the crack will propagate to the direction θ0 expressed by Equation (11). The crack increment length ∆a is taken 10%-20% of the initial crack

length a , inversely proportional to the ratio of K II K I . The ratio represents the mixed mode p roportionality, therefo re shorter increment length should be taken to carefully justify the crack path curvature when 𝐾𝐾𝐼𝐼𝐼𝐼 is relat ively large compare to 𝐾𝐾𝐼𝐼 [13]. Thus the crack length increment is approximated by the Lagrange interpolat ion as:

∆a=

1 −

K II

(20%) +

K II

(10% ) a

(16)

KI

KI

3. Numerical Analysis and Validation

A developed finite element model of crack growth for brittle materials has been employed on four points bending geometry and, double centre cracked plate with two holes under mixed mode loading conditions.

3.1. Four Poi nts Bending Geometry

having the dimensions of B = 2 mm, W = 5 mm and

length > 43 mm, made of Al2O3-Ceramics. An enlargement

of mesh around crack t ip represented the elements around the

crack tip is illustrated in Fig.3. Four points bending geometry

with final adaptive mesh are shown in Fig. 5.

The stresses in this loading case were also computed. The

geometric functions of FI and FII as follo ws[22]:

K

=

P

1

−

d

πaF

I WB L

I

K

=

P

1

−

d

πaF

(17)

II WB L

II

(a)

θ0 x'

θ0

x '

KII Positive

KII Negative

Figure 4. Sign of the propagation angle

(b)

Figure 5. (a) Four-point bending geometry and (b) the final adaptive mesh

In this study, a comparison and analysis were applied between the present results with those empirical results selected from the literature ([3] and[22]). The stress intensity

International Journal of M aterials Engineering 2012, 2(5): 67-74

71

factors evaluation ( K I & K II ), are tabulated in Table 1 and Table 2 with d/W = 0.5. The dimensionless form of the estimated stress intensity factor was obtained by using Eq. (17), presented for the range of 0.1 < a/ W < 0.4.

Table 1. Dimensionless stress intensity factor for the central cracked plate for Mode I

(a/W)

0.1 0.2 0.3 0.4 0.6 0.8

K

=

P

1

−

d

πaF

I WB L

I

Fett, et al.(1995)

0.3450 0.6633 0.9399 1.1702 1.5507 2.0684

present study

0.3461 0.6660 0.9354 1.1721 1.5512 2.0697

In order to show the trends of the dimensionless SIFs with the ratio value of initial crack per width under Mode I conditions. Fig. 6 shows, the comparison between the current study results of the dimensionless SIFs and the emp irical res u lts [22].

Table 2. Dimensionless stress intensity factor for the central cracked plate for Mode II

(a/W)

0.1 0.2 0.3 0.4 0.6 0.8

K

=

P

1

−

d

πaF

II WB L

II

Fett et al. (1995)

0.3841 0.2448 0.1580 0.1098 0.0566 0.0106

present study

0.386 0.248 0.1583 0.1079 0.0570 0.013

The relat ion showed that the values of SIFs and a/W increased from 0.31 with a/W = 0.1 to the maximu m dimensionless SIF value of 2.1 accordingly with a/W = 0.8, as shown in Fig. 6.

Fig. 7. Fro m these results, it indicated that under Mode I loading condition, high stresses were generated when a/W = 0.8, which increased the dimensionless SIFs and also exceeded the maximu m stiffness of this material. Whereas, under Mode II loading condition shows a low effect which influenced the values of the dimensionless SIFs. In fact, this was due to the brittleness of this material.

Fi gure 7. Dimensionless SIF vs. (a/W) for mode II

At the beginning, the stress was higher under Mode I and after that, a stress exceeded the criterion value it’s decreased gradually where produced a low effect under Mode II.

This section presents comparison of the current results with these results of Gürse and Miehe[3] in terms of crack propagation trajectories. Fig. 8 shows, the comparison results of predicted crack gro wth of the current study with those results of[3]. In this case, the crack direction was dominated by Mode I stress intensity factor, which became larger by the crack growth non-linear to the top right of the geometry. Due the shear stress which is generated by of Mode II loading condition.

Figure 8. Comparison results of the crack trajectory; (a) current study (b) Gürse & Miehe[3]

Fi gure 6. Dimensionless SIF vs. (a/W) for mode I

Therefore, under Mode II a condition, the values of the dimensionless SIFs was decreased gradually, where a/W was increased respectively the maximu m value of 0.8 as shown in

Fi gure 9. Three st eps of crack propagat ion direct ion of four-point bending geo met ry

72

Souiyah M iloud et al.: Finite Element M odel of Crack Growth under M ixed M ode Loading

Further, the effectiveness of the shear stress was significant on crack trajectory with different direct ions, which is mostly reflected on the materia l properties as shown in Fig. 9. However, the representation of the deformation was controlled by the user, to show the deformat ion step by step with direct ions through the initial crack propagation until the final crack as illustrated in fig. 9. Obviously, the current results of this developed fin ite element programme has good agreement as shown in fig. 8b.

K I and K II tended to decrease. Therefore, when the crack direction declined fro m both holes with the curvature, both KI and K II decreased. However, when the crack growth linearly K I tended to increased again whereas, K II decreased. Furthermore, based on this phenomenon, it could be understood that Mode I significant effect than Mode II, and this was due to the brittleness of this material.

3.2. Center Double cracked Pl ate with Two-Holes under Mi xed Mode Loading

A comp licated mixed mode fracture problem was studied to demonstrate the performance o f the developed programme on double centre cracked plate with two holes. The geometry and its final adaptive mesh are shown in Fig. 10. Meanwh ile, the enlargement of the rosette around both crack tips is shown in Fig. 11.

Figure 10. The geometry and its final adaptive mesh

Figure 12. The relationship of KI & KII with crack extension

A comparison study between the current study and the experimental results[8] was performed, to validate the accuracy and reliability of this developed FE model. Notably, all the crack lengths were measured along the cracked path. The crack trajectory of the nu merical results for this study is shown in Fig. 13. It appeared that, in the v icin ity of the holes, the crack direct ion was curved as a consequence of the mixed mode (I and II) loading. The crack was propagated non-linearly towards the hole by the attraction of the holes, which are generated high stresses from both sides of the crack t ip.

Figure 11. The enlargement of the rosette around both crack tips

The geometry was made of high ca rbon steel and had the following parameters: Young’s modulus E = 2.1 ×105 MPa and Poisson’s ratio ν = 0.3. The trend showed the relationship between the stress intensity factors of KI and K II as illustrated in Fig. 12.

In the beginning, the crack exh ibited a pure Mode I state, which increased the KI values, whereas in Mode II state, KII values became lower so that the crack curved closely to the left side of the hole and then curved downwards where

Figure 13. Numerical result of the crack path for the CCT specimen with two holes

These generated stresses obviously influenced the crack propagation direction as well as the values of SIFs. The

International Journal of M aterials Engineering 2012, 2(5): 67-74

73

experiment results[8] revealed the same behaviour of the crack growth trajectory as illustrated in Fig. 14. Furthermore, Fig. 15 shows the symmetric results of the maximu m principal stress with the crack propagation trajectory.

background mesh. Overall, all these advantages were successfully revealed the reliability and the capability of this new developed FE model in dealing with plane crack behaviours. In fact, in term of flexib ility the user could also control the problem type, which is either a plane strain or a plane stress to clarify the deformation. The comparisons and analysis shown that, this developed program was truly evaluated and validated with relevant research works under different methods selected from the literature.

Fi gure 14. Current work result s of an enlargement of a crack traject ory for the CCT specimen with two holes

Figure 15. The triangle points correspond to experimentally obtained crack path, where the full line holds for one parameter and the dash line for two parameter fracture mechanics[8]

The predicted crack trajectory fro m the numerical results of the current results was obviously shows, a similarity and a good agreement to the experimental results[8]. In fact, the crack direction is influenced whenever there is hole in the geometry and this is due to the huge of stresses were generated around the holes.

4. Conclusions

In this paper, a co mprehensive Finite Element model was developed using the source code written in FORTRAN© with an advancing front method for crack growth analysis. Displacement Extrapolat ion Technique (DET) was emp loyed to evaluate the SIFs under mixed mode loading conditions. The automatic crack propagation was characterised by the successive propagation steps performed without user interaction. The developed FE programme allo wed the user to control the process by changing the size of the crack gro wth increment, maximu m and minimu m element size in the mesh and init ial mesh size for the

REFERENCES

[1] Sukla, A. Practical Fracture Mechanics in Design. 2nd Edition. New York: M arcel Dekker (2005).

[2] Fan, S.C., Liu, X., Lee, C.K. Enriched partition-of-unity finite element method for stress intensity factors at crack tips. Comp. & Struct. 82, 445-461(2004).

[3] Gürses, E., M iehe, C. A computational framework of three-dimensional conﬁgurational-force-driven brittle crack propagation. Comp. M eth. Appl. M ech. Eng., 198, 1413–1428(2009).

[4] Wen, P.H., Aliabadi, M .H. A variational approach for evaluation of stress intensity factors using the element free Galerkin method. Int. J. of Sol. & Stru., 8, 1171–1179(2011).

[5] Xie, Y.J., Lee, K.Y., Hu, X.Z., Cai, Y.M . Applications of conservation integral to indentation with a rigid punch. Eng. Fract. M ech. 76, 949–957(2009).

[6] Abdelaziz ,Y., Benkheira, S., Rikioui, T., M ekkaoui, A. A double degenerated ﬁnite element for modeling the crack tip singularity, App. M ath. M odl., 34, 4031–4039(2010).

[7] Fortino, S., Bilotta, A. Evaluation of the amount of crack growth in 2D LEFM problem. Eng. Fract. M ech. 71, 1403–19(2004).

[8] Stanislav, S., Zdenek, K. Two parameter fracture mechanics: Fatigue crack behavior under mixed mode conditions. Eng. Fract. M ech. 75, 857(2008).

[9] Colombo, D., Giglio, M . A methodology for automatic crack propagation modelling in planar and shell FE models. Eng. Fract. M ech. 73, 490-504(2006).

[10] Freese, C.E., Tracey, D.M . The natural triangle versus collapsed quadrilateral for elastic crack analysis. Int. J. of Fract. 12, 767-770 (1976).

[11] Barsoum, R. S. On the use of isoparametric finite elements in linear fracture mechanics. Int. J. for Num. M eth. in Eng., 10, 25-37(1976).

[12] Xie, M ., Gerstle, W.H., Rahulkumar, P. Energy-based automatic mixed-mode crack-propagation modeling. J. of Eng. M ech. ASCE. 121, 914–923(1995).

[13] T Bittencourt, N., Wawrzynek, P.A., Ingraffea, A.R., Sousa, J.L. Quasi-automatic simulationof crack propagation for 2D LEFM problems, Eng. Fract. M ech., 55, 321–334(1996).

[14] Chiou, Y.J., Lee, Y.M ., Jowtsay, R. M ixed mode fracture propagation by manifold method. Int. J. of Fract. 114,

74

Souiyah M iloud et al.: Finite Element M odel of Crack Growth under M ixed M ode Loading

327–347(2002).

[15] Guan, Z.Q., Song, C., Gu, Y.X. Recent advances of research on finite element mesh generation methods. J. of Comp. Aided Design and Comp. Graph. 15 (1), 1–14(2003).

[16] Zienkiewicz, O., Tay lor, R., Zhu, J. The finite elemen t method: its basis and fundamenta. 2005.

[17] Ariffin, A.K. : PhD Thesis, University of Wales Swansea. 1995.

[18] Phongthanapanich, S., Dechaumphai, P. Adaptive Delaunay triangulation with object oriented programming for crack propagation analysis. Fin. Elem. Anal. Des. 40, 1753–1771(2004).

[19] Araújo, T., Bittencourt, T., Roehl, D., M artha, L. Numerical

estimation of fracture parameters in elastic and elastic-plastic analysis, European Congress on Computational M ethods in Applied Sciences and Engineering, 11-14 September. Barcelona (2000).

[20] Anlas, G., Santare, M ., Lambros, J. Numerical calculation of stress intensity factors in functionally graded materials. Int. J. of Fract. 104, 131–143(2000).

[21] Guinea, G.V., Planan, J., Elices, M . KI evaluation by the displacement extrapolation technique. Eng. Fract. M ech. 66, 243-255(2000).

[22] Fett,T., Gerteisen, G., Hahnenberger, S., M artin, G., M unz, D. Fracture tests for ceramics under mode-I, mode-II and mixed-mode loading. J. of the Eur.Cer. Soc., 5(4), 307-312(1995)

Finite Element Model of Crack Growth under Mixed Mode Loading

Souiyah Miloud1,*, A. Muchtar2, A. K. Ariffin2, Malek Ali1, M. I. Fadhe l1, Base m Abu Zne id3

1Faculty of Engineering & Technology (FET), M elaka, 75450, M ultimedia University (MMU) 2Faculty of Engineering & Built Environment, Bangi, 43000, University Kebangsaan M alaysia (UKM )

3Universiti Teknologi M alaysia Skudai, 81310, Johor Bahru

Abstract In this paper, in order to predict the crack growth trajectory and to evaluate the SIF under mixed modes (I & II),

one proposes a new finite element program for crack growth using the source code written in FORTRAN. The fin ite element mesh is generated using an advancing front method, where the generation of the background mesh and the construction of singular elements are also added to this developed programme to facilitate the crack process and the fracture analysis. Displacement Ext rapolation Technique (DET) was employed to evaluate the SIFs under mixed mode loading conditions. Therefore, the accuracy of both SIF`s values and the crack path predictions results are compared and validated with other relevant published research work. Ho wever, the assessment indicated that this developed fin ite element programme is reliable and robust to evaluate the SIFs and predicts the crack trajectories successfully based on the applied loading conditions.

Keywords Linear Elastic Fracture Mechanics, Adaptive Refinement Mesh, Stress Intensity Factors, Finite Elements

Method, Crack Propagation

1. Introduction

Structural failu re can be generally associated with one or more fracture of the materials making that structure. In fact, fractu re mechanics deals with the irrevers ible process of rupture due to nucleation or sudden crack and crack gro wth. The influence of cracking on the structural strength has been wid ely ap p reciated s in ce th e end o f th e 19th centu ry. Ho wever, so me aspects o f its natu re and influen ce st ill remain un kno wn[1]. Th e use o f crack p ropagat ion laws based on stress intensity factor range is the most successful engineering applicat ion of fracture mechanics. In the elastic fracture analysis, the stress intensity factors sufficient ly define the stress field close to the crack tip and prov ide fundamental in fo rmat io n o f ho w th e crack is go ing to p rop ag at e. Bas ically , th e est imat io n meth ods can be categorized into two groups , thos e bas ed on field extrapolation near the crack t ip and those which make use of the energy release when the crack propagates. The latter group includes the J-contour integration, the virtual crack extension and the strain energy release rate method. The main d isadv antag e o f these methods is th at th e st ress intens ity factor co mponents, KI and KII in mi xed mode prob lems are eit her imposs ib le o r very d ifficu lt t o be

* Corresponding author: [email protected] (Souiyah Miloud) Published online at http://journal.sapub.org/ijme Copyright © 2012 Scientific & Academic Publishing. All Rights Reserved

separated. Nevertheless, the first groups which are based on near-tip field fitting procedures require finer meshes to produce a good numerical representation of crack-tip fields generated to facilitate the calculat ion.

2. Finite Element Model

Most established research was on computer base by applying numerical methods which were mostly uses commercial software such as;[2] an enriched partit ion of unity finite element method (PUFEM ), which is known as one of the meshless method to calculate the SIF in LEFM under plane stress and plane strain condition. Another, numerical examp les were performed to evaluate the generated stress intensity factors directly fro m the scaled boundary ﬁnite element solution for the singular stress ﬁeld[3]. More recent work was developed variat ional mesh-free method, to evaluate the stress intensity factors of mixed mode crack problems, using the element free Galerkin method[4]. Accordingly, another new method was proposed to determine the SIFs for the indentation problem based on the conservation integral[5]. In contrast, Numerical experiments were conducted to evaluate the effectiveness of two proposed techniques on near crack tip singularity[6]. In fact, an observation was stated that a general theoretical model for automatically evaluating the increments of crack growth during a loading process does not exist yet[7]. Expe riments work on a central crack specimen with holes, to predict the crack path and to evaluate the mixed mode stress

68

Souiyah M iloud et al.: Finite Element M odel of Crack Growth under M ixed M ode Loading

intensity factors (KI & KII)[8]. Somehow, this limit nowadays still characterizes the

computer codes for crack growth simulat ion in elastic and elasto-plastic materials. Generally, the software should address three different aspects of the problem, specifically are; (i) SIFs evaluation, (ii) the crack propagation direction, and (iii) mesh modification to accommodate the crack advancement. The first and second aspects can be solved with no particular d ifficult ies and many comme rcial FE codes have built-in capabilit ies to evaluate the SIF. The last aspect is much mo re co mplex and is rarely solved by FE codes, e.g. FRANC2D, FRANC3D and ZENCRA CK[9]. Eventually, the FE model is modified in order to accommodate the obtained evaluation on short and straight crack propagation directions. The whole process is repeated many times and both is time-consuming and a source of errors if it is performed manually. This is part icularly t rue for the mesh modification. However, software on the other hand, is capable in handling the solution process in an automat ic way.

However, in the current study, we developed a new fin ite element programme for crack growth using the source code written in FORTRAN. In order to predict the crack growth trajectory and to evaluate the SIFs under mixed modes (I & II) loading in the frame of LEFM. The displacement extrapolation techniques with the adaptive refinement mesh method are used, to determine the stress intensity factors on two different geometries. Specifically are; four points bending plate and, centre double cracked plate with two-holes. Additionally, the obtained results of the current study are evaluated and validated with other relevant experimental and numerical results selected from the literature.

2.1. Generation of Singular Element

An unstructured triangular mesh is automatically generated by employing the advancing front method. The singular elements have to be constructed correctly to get a proper field of singularity around the crack tip as shown in Fig. 1.

Figure 1. The cut and patch procedure of generating singular elements around a crack tip

The number o f elements depends on the distributed nodes around the crack tip, wh ich can be set by the user as shown in Fig. 2.

In fact, the natural triangular quarter point elements are used[10] instead of the collapsed quadrilateral element, which is proposed by Barsoum[11]. The success of the adaptivity in general depends to a large extent on the efficient coupling between the error estimator, refinement scheme and automatic mesh generator. The importance of

these adaptive techniques in practical applications has led to a considerable research on fully automatic mesh generators that require only the specificat ion of the boundary and mesh size d istribution over the domain. Moreover, several comprehensive numerical tools have been developed to enhance the accuracy at the crack tip.

c b d

e

v,y

u,x

L

3L

4

4

Figure 2. A typical arrangement of the natural triangle quarter-point element s around a crack t ip

In this work, the unstructured triangle mesh is automatically generated by employing the advancing front method. Many researchers applied the FEM with remeshing to study the fracture propagation and its SIFs analysis ([12]; and[13]), though it is not an easy task. In contrast, it is almost impossible to automatically remesh finite elements of an arbitrarily growing crack[14]. One of the advantages of the advancing front method is that the new triangle element formation is coinstantaneous with new node generation, and this advantage makes it possible to control the shape and size of the element through the adjustment of the location of the node[15]. However, a lot o f intersection checking between the new generated triangular elements and the existing elements must be computed in order to ensure that the triangular elements are valid. Accordingly, an observation was stated that the algorithm for the advancing front method has been shown to be robust in two-dimensional mesh generation of triangles and could be extended easily to generate quadrilaterals[16].

2.2. Mesh Generati on and Adapti ve Refinement Mesh

In general, the s maller mesh size g ives more accurate fin ite element appro ximate solution. Ho wever, reduction in the mesh size leads to greater computational effort. The adaptive mesh refinement is emp loyed as the optimization scheme. This scheme bases on a posteriori error estimator which is obtained from the solution fro m the previous mesh. Basically, the success of the adaptivity in overall is depends to a large extent on the efficient coupling between the error estimator, refinement scheme and automatic mesh generator.

International Journal of M aterials Engineering 2012, 2(5): 67-74

69

The importance of these adaptive techniques in practical applications has led to a considerable research on fully automatic mesh generators that require only the specification of the boundary and mesh size distribution over the do main. An examp le of geo metry shows the whole process of generating the mesh is illustrated in Fig.3.

(a)

(b)

(c)

(d)

(ii) Determine the average error norm over the whole domain

∑ ∫ =eˆ 1 m σ Tσd Ω

(2)

m e =1 Ωe

Where m is the total number of elements in the whole do

ain. (iii)Determine a variable, εe for each element as

e 1/ 2

( ) 1 e

( ) εe = η eˆ 1/ 2

(3)

Where η is a percentage that measures the permissible

error for each element. If εe > 1 the size of the element is reduced and vice versa. (iv) The new element size is

determined as

hˆ = he

e

(4)

(εe )1/ p

Where he is the old element size and p is the order of the interpolation shape function.

2.3. The Dis placement Extrapolation Techni que

(g)

(e)

(f)

Figure 3. The mesh generation stages

The geometry of a p late with six holes and two notches is

illustrated Fig. 3a. Six connector lines as shown in Fig. 3b,

are fo rcing the internal boundaries to be the continuous part

of the external boundary. Fig. 3c shows the cutting out of the

rosette templates around each crack tip. The background

mesh for this domain is then set up automatically using

dichotomy technique as shown in Fig.3d and Fig. 3e shows

the conventional mesh being generated by the advancing

front method. The first generation produces mesh with initial

size set by user. Later, during adaptive refinement, this first

generated mesh will be taken as the background mesh. In Fig.

3f, for each rosette template, quarter-point elements are then

constructed. Fig.3g, shows the enlargement of the

quarter-point element at each crack t ip. In general, the

smaller mesh size gives more accurate fin ite element

approximate solution. However, reduction in the mesh size

leads to greater computational effort.

The strategy used to refine the mesh during analysis

process as follo ws[17]:

(i) Determine the error norm for each element

∫ ( ) ( ) e e = σ −σ * T σ −σ * d Ω

(1)

Ωe

Where σ is the stress field obtained fro m the fin ite element

calculation and σ* is the smoothed stress field inclusion of

quarter point elements.

One of the simplest and most frequently used methods is

displacement ext rapolation technique[18]. It consists

typically in the effective SIF concept by which, the fracture

evaluation can be easily carried out[19]. The crack length

was evaluated by linear-elastic analysis fro m the compliance

of single-edge-notched specimen in three-points bending

test[20]. A new ﬁn ite element of quasi-static brittle fracture

was developed, based on computational framework for

three-dimensional cracks propagation bodies. They uses

consistent thermodynamical framework for crack

propagation in elastic solids[3].

The asymptotic expression for the displacement norma l to

crack p lane, v under mode I loading[21]:

=v KI 14+Eν 2πr (2κ +1) sin θ2 − sin 32θ + A1(1E+ν )r (κ − 3) sinθ (5)

+

A2 (1+ν )r3

2

(2κ

−1)

sin

3θ

− sin θ

+ ...

E 3

2

2

Where KI is the stress intensity factor for mode I, E is the

modulus of elasticity, v is the Poisson’s ratio, K an elastic

parameter defined in equation (6), Ai are parameters

depending on the geometry and load on the specimen, and

r and θ are the polar coordinates defined at the crack t ip.

(3 − 4ν )

plane stress

κ = (3 − 4ν ) / (1+ν ) planestrain (6)

The near tip nodal displacements at nodes b, c, d and e

shown below in Fig. 3, are of interest. The displacements are

e xtrapolated by evaluating Equation (5) a long the crack faces

𝜃𝜃 = ± 𝜋𝜋. Particu larizing for nodes b and c on the singular element at the upper face of the crack yields:

= v K 2 (1+ν )(κ +1) L − A2 (1+ν )(κ +1) L3 2 + O(L5 2 ) (7)

b

Iπ

4E

12E

70

Souiyah M iloud et al.: Finite Element M odel of Crack Growth under M ixed M ode Loading

v= K 2 (1+ν )(κ +1) L − 2A2 (1+ν )(κ +1) L3 2 + O(L5 2 ) (8)

c

Iπ

2E

3E

Where L, is the length of the element side which is

connected to the crack tip. Equations (7) and (8) can be

solved to obtain the value of K I as:

= K E

2π (8v −v )

(9)

I 3(1+ν )(1+ κ ) L

b

c

The nodal displacements at the other two nodes can be evaluated by the similar way. If one considers the all four nodal displacements, the stress intensity factor expression b eco mes :

KI

E

2π

4

(v

b

−v ) − (v c d

−v e )

(10)

3(1+ν )(1+ κ ) L

2

Under pure mode II loading, the displacement tangential

to the crack plane, u is given by:

=u K 1+ν

2r

(2κ +1) sin θ

−

sin

3θ

+

A1 (1 +ν

)r

(κ

−

3) sinθ

(11)

II 4E π

2

2

E

+

A2 (1+ν )r 3

2

(2κ

−1)

sin

3θ

− sin θ

+ ...

E 3

2

2

Similarly, the stress intensity factor for mode II, using the

nodal displacements of the four nodes can be estimated by:

K II

E

2π

4

(u

b

−ud

)−

(uc

−ue )

(12)

3(1+ν )(1+ κ ) L

2

In order to simulate crack propagation under linear elastic

condition, the crack path direction must be determined.

There are several methods use to predict the d irection of

crack trajectory such as the maximu m circu mferential stress

theory, the maximu m energy release rate theory and the

minimu m strain energy density theory. The maximu m

circu mferential stress theory asserts that, for isotropic

materials under mixed-mode loading, the crack will

propagate in a direction normal to maximu m tangential

tensile stress. In polar coordinates, the tangential stress is

given by

σ

1 cos θ K cos2 θ − 3 K sin θ (13)

θ 2π r 2 I

2 2 II

The direct ion normal to the maximu m tangential stress can

be obtained by solving d σθ / d θ = 0 for θ . The nontrivial solution is given by

K I sinθ + K II (3cosθ −1) =0

(14)

Which can be solved as:

−1

3K

2 II

+KI

K

2 I

+

8K

2 II

θ0 = ± cos

2

2

(15)

K I + 9K II

y'

y'

In order to ensure that the opening stress associated with the crack direction of the crack extension is maximu m, the

sign of θ0 should be opposite to the sign of K II [19]. The

two possibilities are illustrated in Fig. 4. The criterion for crack to propagate fro m crack t ip is based

on the material toughness, KC . If the calculated stress

intensity factor, K I ≥ KC then the crack will propagate to the direction θ0 expressed by Equation (11). The crack increment length ∆a is taken 10%-20% of the initial crack

length a , inversely proportional to the ratio of K II K I . The ratio represents the mixed mode p roportionality, therefo re shorter increment length should be taken to carefully justify the crack path curvature when 𝐾𝐾𝐼𝐼𝐼𝐼 is relat ively large compare to 𝐾𝐾𝐼𝐼 [13]. Thus the crack length increment is approximated by the Lagrange interpolat ion as:

∆a=

1 −

K II

(20%) +

K II

(10% ) a

(16)

KI

KI

3. Numerical Analysis and Validation

A developed finite element model of crack growth for brittle materials has been employed on four points bending geometry and, double centre cracked plate with two holes under mixed mode loading conditions.

3.1. Four Poi nts Bending Geometry

having the dimensions of B = 2 mm, W = 5 mm and

length > 43 mm, made of Al2O3-Ceramics. An enlargement

of mesh around crack t ip represented the elements around the

crack tip is illustrated in Fig.3. Four points bending geometry

with final adaptive mesh are shown in Fig. 5.

The stresses in this loading case were also computed. The

geometric functions of FI and FII as follo ws[22]:

K

=

P

1

−

d

πaF

I WB L

I

K

=

P

1

−

d

πaF

(17)

II WB L

II

(a)

θ0 x'

θ0

x '

KII Positive

KII Negative

Figure 4. Sign of the propagation angle

(b)

Figure 5. (a) Four-point bending geometry and (b) the final adaptive mesh

In this study, a comparison and analysis were applied between the present results with those empirical results selected from the literature ([3] and[22]). The stress intensity

International Journal of M aterials Engineering 2012, 2(5): 67-74

71

factors evaluation ( K I & K II ), are tabulated in Table 1 and Table 2 with d/W = 0.5. The dimensionless form of the estimated stress intensity factor was obtained by using Eq. (17), presented for the range of 0.1 < a/ W < 0.4.

Table 1. Dimensionless stress intensity factor for the central cracked plate for Mode I

(a/W)

0.1 0.2 0.3 0.4 0.6 0.8

K

=

P

1

−

d

πaF

I WB L

I

Fett, et al.(1995)

0.3450 0.6633 0.9399 1.1702 1.5507 2.0684

present study

0.3461 0.6660 0.9354 1.1721 1.5512 2.0697

In order to show the trends of the dimensionless SIFs with the ratio value of initial crack per width under Mode I conditions. Fig. 6 shows, the comparison between the current study results of the dimensionless SIFs and the emp irical res u lts [22].

Table 2. Dimensionless stress intensity factor for the central cracked plate for Mode II

(a/W)

0.1 0.2 0.3 0.4 0.6 0.8

K

=

P

1

−

d

πaF

II WB L

II

Fett et al. (1995)

0.3841 0.2448 0.1580 0.1098 0.0566 0.0106

present study

0.386 0.248 0.1583 0.1079 0.0570 0.013

The relat ion showed that the values of SIFs and a/W increased from 0.31 with a/W = 0.1 to the maximu m dimensionless SIF value of 2.1 accordingly with a/W = 0.8, as shown in Fig. 6.

Fig. 7. Fro m these results, it indicated that under Mode I loading condition, high stresses were generated when a/W = 0.8, which increased the dimensionless SIFs and also exceeded the maximu m stiffness of this material. Whereas, under Mode II loading condition shows a low effect which influenced the values of the dimensionless SIFs. In fact, this was due to the brittleness of this material.

Fi gure 7. Dimensionless SIF vs. (a/W) for mode II

At the beginning, the stress was higher under Mode I and after that, a stress exceeded the criterion value it’s decreased gradually where produced a low effect under Mode II.

This section presents comparison of the current results with these results of Gürse and Miehe[3] in terms of crack propagation trajectories. Fig. 8 shows, the comparison results of predicted crack gro wth of the current study with those results of[3]. In this case, the crack direction was dominated by Mode I stress intensity factor, which became larger by the crack growth non-linear to the top right of the geometry. Due the shear stress which is generated by of Mode II loading condition.

Figure 8. Comparison results of the crack trajectory; (a) current study (b) Gürse & Miehe[3]

Fi gure 6. Dimensionless SIF vs. (a/W) for mode I

Therefore, under Mode II a condition, the values of the dimensionless SIFs was decreased gradually, where a/W was increased respectively the maximu m value of 0.8 as shown in

Fi gure 9. Three st eps of crack propagat ion direct ion of four-point bending geo met ry

72

Souiyah M iloud et al.: Finite Element M odel of Crack Growth under M ixed M ode Loading

Further, the effectiveness of the shear stress was significant on crack trajectory with different direct ions, which is mostly reflected on the materia l properties as shown in Fig. 9. However, the representation of the deformation was controlled by the user, to show the deformat ion step by step with direct ions through the initial crack propagation until the final crack as illustrated in fig. 9. Obviously, the current results of this developed fin ite element programme has good agreement as shown in fig. 8b.

K I and K II tended to decrease. Therefore, when the crack direction declined fro m both holes with the curvature, both KI and K II decreased. However, when the crack growth linearly K I tended to increased again whereas, K II decreased. Furthermore, based on this phenomenon, it could be understood that Mode I significant effect than Mode II, and this was due to the brittleness of this material.

3.2. Center Double cracked Pl ate with Two-Holes under Mi xed Mode Loading

A comp licated mixed mode fracture problem was studied to demonstrate the performance o f the developed programme on double centre cracked plate with two holes. The geometry and its final adaptive mesh are shown in Fig. 10. Meanwh ile, the enlargement of the rosette around both crack tips is shown in Fig. 11.

Figure 10. The geometry and its final adaptive mesh

Figure 12. The relationship of KI & KII with crack extension

A comparison study between the current study and the experimental results[8] was performed, to validate the accuracy and reliability of this developed FE model. Notably, all the crack lengths were measured along the cracked path. The crack trajectory of the nu merical results for this study is shown in Fig. 13. It appeared that, in the v icin ity of the holes, the crack direct ion was curved as a consequence of the mixed mode (I and II) loading. The crack was propagated non-linearly towards the hole by the attraction of the holes, which are generated high stresses from both sides of the crack t ip.

Figure 11. The enlargement of the rosette around both crack tips

The geometry was made of high ca rbon steel and had the following parameters: Young’s modulus E = 2.1 ×105 MPa and Poisson’s ratio ν = 0.3. The trend showed the relationship between the stress intensity factors of KI and K II as illustrated in Fig. 12.

In the beginning, the crack exh ibited a pure Mode I state, which increased the KI values, whereas in Mode II state, KII values became lower so that the crack curved closely to the left side of the hole and then curved downwards where

Figure 13. Numerical result of the crack path for the CCT specimen with two holes

These generated stresses obviously influenced the crack propagation direction as well as the values of SIFs. The

International Journal of M aterials Engineering 2012, 2(5): 67-74

73

experiment results[8] revealed the same behaviour of the crack growth trajectory as illustrated in Fig. 14. Furthermore, Fig. 15 shows the symmetric results of the maximu m principal stress with the crack propagation trajectory.

background mesh. Overall, all these advantages were successfully revealed the reliability and the capability of this new developed FE model in dealing with plane crack behaviours. In fact, in term of flexib ility the user could also control the problem type, which is either a plane strain or a plane stress to clarify the deformation. The comparisons and analysis shown that, this developed program was truly evaluated and validated with relevant research works under different methods selected from the literature.

Fi gure 14. Current work result s of an enlargement of a crack traject ory for the CCT specimen with two holes

Figure 15. The triangle points correspond to experimentally obtained crack path, where the full line holds for one parameter and the dash line for two parameter fracture mechanics[8]

The predicted crack trajectory fro m the numerical results of the current results was obviously shows, a similarity and a good agreement to the experimental results[8]. In fact, the crack direction is influenced whenever there is hole in the geometry and this is due to the huge of stresses were generated around the holes.

4. Conclusions

In this paper, a co mprehensive Finite Element model was developed using the source code written in FORTRAN© with an advancing front method for crack growth analysis. Displacement Extrapolat ion Technique (DET) was emp loyed to evaluate the SIFs under mixed mode loading conditions. The automatic crack propagation was characterised by the successive propagation steps performed without user interaction. The developed FE programme allo wed the user to control the process by changing the size of the crack gro wth increment, maximu m and minimu m element size in the mesh and init ial mesh size for the

REFERENCES

[1] Sukla, A. Practical Fracture Mechanics in Design. 2nd Edition. New York: M arcel Dekker (2005).

[2] Fan, S.C., Liu, X., Lee, C.K. Enriched partition-of-unity finite element method for stress intensity factors at crack tips. Comp. & Struct. 82, 445-461(2004).

[3] Gürses, E., M iehe, C. A computational framework of three-dimensional conﬁgurational-force-driven brittle crack propagation. Comp. M eth. Appl. M ech. Eng., 198, 1413–1428(2009).

[4] Wen, P.H., Aliabadi, M .H. A variational approach for evaluation of stress intensity factors using the element free Galerkin method. Int. J. of Sol. & Stru., 8, 1171–1179(2011).

[5] Xie, Y.J., Lee, K.Y., Hu, X.Z., Cai, Y.M . Applications of conservation integral to indentation with a rigid punch. Eng. Fract. M ech. 76, 949–957(2009).

[6] Abdelaziz ,Y., Benkheira, S., Rikioui, T., M ekkaoui, A. A double degenerated ﬁnite element for modeling the crack tip singularity, App. M ath. M odl., 34, 4031–4039(2010).

[7] Fortino, S., Bilotta, A. Evaluation of the amount of crack growth in 2D LEFM problem. Eng. Fract. M ech. 71, 1403–19(2004).

[8] Stanislav, S., Zdenek, K. Two parameter fracture mechanics: Fatigue crack behavior under mixed mode conditions. Eng. Fract. M ech. 75, 857(2008).

[9] Colombo, D., Giglio, M . A methodology for automatic crack propagation modelling in planar and shell FE models. Eng. Fract. M ech. 73, 490-504(2006).

[10] Freese, C.E., Tracey, D.M . The natural triangle versus collapsed quadrilateral for elastic crack analysis. Int. J. of Fract. 12, 767-770 (1976).

[11] Barsoum, R. S. On the use of isoparametric finite elements in linear fracture mechanics. Int. J. for Num. M eth. in Eng., 10, 25-37(1976).

[12] Xie, M ., Gerstle, W.H., Rahulkumar, P. Energy-based automatic mixed-mode crack-propagation modeling. J. of Eng. M ech. ASCE. 121, 914–923(1995).

[13] T Bittencourt, N., Wawrzynek, P.A., Ingraffea, A.R., Sousa, J.L. Quasi-automatic simulationof crack propagation for 2D LEFM problems, Eng. Fract. M ech., 55, 321–334(1996).

[14] Chiou, Y.J., Lee, Y.M ., Jowtsay, R. M ixed mode fracture propagation by manifold method. Int. J. of Fract. 114,

74

Souiyah M iloud et al.: Finite Element M odel of Crack Growth under M ixed M ode Loading

327–347(2002).

[15] Guan, Z.Q., Song, C., Gu, Y.X. Recent advances of research on finite element mesh generation methods. J. of Comp. Aided Design and Comp. Graph. 15 (1), 1–14(2003).

[16] Zienkiewicz, O., Tay lor, R., Zhu, J. The finite elemen t method: its basis and fundamenta. 2005.

[17] Ariffin, A.K. : PhD Thesis, University of Wales Swansea. 1995.

[18] Phongthanapanich, S., Dechaumphai, P. Adaptive Delaunay triangulation with object oriented programming for crack propagation analysis. Fin. Elem. Anal. Des. 40, 1753–1771(2004).

[19] Araújo, T., Bittencourt, T., Roehl, D., M artha, L. Numerical

estimation of fracture parameters in elastic and elastic-plastic analysis, European Congress on Computational M ethods in Applied Sciences and Engineering, 11-14 September. Barcelona (2000).

[20] Anlas, G., Santare, M ., Lambros, J. Numerical calculation of stress intensity factors in functionally graded materials. Int. J. of Fract. 104, 131–143(2000).

[21] Guinea, G.V., Planan, J., Elices, M . KI evaluation by the displacement extrapolation technique. Eng. Fract. M ech. 66, 243-255(2000).

[22] Fett,T., Gerteisen, G., Hahnenberger, S., M artin, G., M unz, D. Fracture tests for ceramics under mode-I, mode-II and mixed-mode loading. J. of the Eur.Cer. Soc., 5(4), 307-312(1995)