An adaptive finite element procedure for crack propagation analysis

228Alshoaibi et al. 1 J Zhejiang Univ SciA 2007 8(2);228-236Journal of Zhejiang University SCIENCE AISSN 1009 3095 (Prit); ISSN 1862-1775 (Online)www. zju.edu.cn/jzus; www.springerlink.comJzusE-mail: jzus@zju.edu.cnAn adaptive finite element procedure for crack propagation analysisALSHOAIBI Abdulnaser M.*, HADI M.S.A., ARIFFIN A.K.(Department of Mechanical & Materials Engineering, Universiti Kebangsaan Malaysia, 43600 Bangi, Selangor Darul Ehsan, Malaysia)'E-mail: alshoaibi@gmail.com; alhager01@yahoo.comReceived May 14, 2006; revision accepted Oct. 10, 2006problems. The mesh is generated by the advancing front method and the norm stress error is taken as a posteriori error estimator forthe h-type adaptive refinement. The stress intensity factors are estimated by a displacement extrapolation technique. The nearcrack tip displacements used are obtained from specific nodes of natural six-noded quarter-point elements which are generatedaround the crack tip defined by the user. The crack growth and its direction are determined by the calculated stress intensity factors.The maximum circumference theory is used for the ltter. In evaluating the accuracy of the estimated stress intensity factors, fourcases are tested consisting of compact tension specimen, three-point bending specimen, central cracked plate and double edgenotched plate. These were carried out and compared to the results from other studies. The crack trajectories of these specimen testsare also ilustrated.Key words: Linear elastic fracture mechanics, Adaptive refinement, Stress intensity factors, Crack propagationdoi: 10.163 l/jzus.2007.A0228Document code: ACLC number: O346.1; TB303INTRODUCTIONments known as quarter-point elements are used toprovide the element polynomial representing the st-The finite element method (FEM) has proved to resses and strainsfield singularity near the crack tip.be very well suited for the study of fracture mechanics. Nodal relaxation is used to release nodes, one by one,Nevertheless, modelling the propagation of a crack in order to enable the crack tip to propagate throughthrough a finite element mesh turns out to be difficult the mesh. This method which is based on near-tipbecause of the modification of the mesh topology. field fitting procedures requires finer meshes to pro-Use of crack propagation laws based on stress inten- duce a good numerical representation of crack-tipsity factor range is the most successful engineering fields with the most accurate methods being thoseapplication of fracture mechanics. The stress intensity based on nodal displacements, which are a primaryfactors are a very important parameter in fracture output of the finite element program (Guinea et al,analysis. These factors define the stress field close to 2000). In adaptive mesh refinement, most analyststhe crack tip of a crack and provide fundamental in-favour either the Delaunay technique or the advanc-formation on how the crack is going to propagate. In ing front method over other techniques when gener-linear elastic fracture problem, the prediction of the ating meshes due to the quality of unstructuredcrack growth and the crack direction are determined meshes generated (EI-Hamalawi, 2004). The mainby the stress intensity factor. Several methods have advantage of the advancing front method is that itbeen proposed to numerically estimate stress intensity tends to produce nicely graded meshes and highfactors using FEM such as the displacement ex- quality triangles that are usually very close in shape totrapolation technique, the J-integral and the energy equilaterals. The boundary integrity is also preserved,domain integral. The first method is used when thesince the discretisation of the domain boundary con-singular elements are present at the crack tip. In the stitutes the initi中国煤化工to the De-finite element fracture analysis, these special ele- launay triangulagrity is not .HCNMH(.Alshoaibi et al. 1 J Zhejiang Univ SciA 2007 8(2);228-236229usually preserved for complicated domains, which is each internal boundary to the external boundary. Thisa key requirement for mesh generation procedures will force the internal boundaries to be part of a con-(EI-Hamalawi, 2004).tinuous line of the externals and therefore set theThe main objective of this paper is to determine computational domain to be a polygon. In order to dohe stress intensity factor for crack propagation this, the orientation direction of internal boundaries isproblem under linear elastic fracture analysis using set clockwise while for external boundary is set thethe displacement extrapolation technique with adap- other way round. The connector line is introduced bytive FEM. The computational code is written in finding the nearest distance between any internalFORTRAN. The mesh for finite elements is the boundary points to any of the external ones (Sezer andunstructuredype, generated using the advancing Zeid, 1991).front method. The global h-type adaptive mesh isIn the proposed dichotomy method, the divisionadopted based on the norm stress error estimator. The starts at any first found boundary point with large facequarter-point singular elements are uniformly gener- angle and setting an angle range for searching theated around the crack tip in the form of a rosette. The nearest nonadjacent point to be connected by a divi-displacement extrapolation technique used in the sion line. The angle range is set in such a way that thecalculation is explained. The algorithm is assessed by division can produce high quality polygon subsetconsidering four standard test specimen geometries，shape. If the search for the nearest nonadjacent pointi.e. compact tension specimen, three-point bending failed, then the division can be initiated at a boundaryspecimen, double edge notched plate, central cracked point with smaller face angle. The classification orderplate, and single edge cracked plate.of the face angle size 0 is set as π≤θ<2π, π/2<θ<π,0<0≤n/2.In order to properly represent the field singular-ADAPTIVE MESH GENERATION AND RE- ity around the crack tip, the singular elements have toFINEMENTbe constructed as well. Since the advancing frontmethod generates the triangle elements starting fromIn this work, the unstructured triangle mesh is the boundary faces, the area around the crack tip forautomatically generated by employing the advancing the construction of the singular elements is supposedfront method (Lohner, 1997). This technique however to be isolated. This area is isolated by first generatingrequires. generating background mesh in order to nodes around the crack tip in the rosette form and thenaccurately control the ditribution of the geometrical the crack tip node and the jointed boundary segmentscharacteristics such as the element size, element are removed. New boundary segments are then in-stretching and stretching directions for the new mesh. troduced linking all the new nodes to temporarily‘cutThe background does not have to be precisely repre- out'’ the template area from the original domain.senting the geometry; however the accuracy of the Subsequently the advancing front triangulation can beditribution depends on this excellence and it must executed. Finally singular elements are 'patched' intocompletely cover the computational domain (Zien- the rosette template to complete the process. Thiskiewicz et al, 2005). Here the strategy taken to gen- procedure is ilustrated by Fig.1. The numbers oferate the background mesh is to utilize all the initial elements depend on the distributed nodes around theboundary nodes of geometry and construct the crack tip, which can be set by user.boundary triangles as the background mesh by thedichotomy technique. In this technique the computa-tional domain must be a polygon since the boundarytriangulations are carried out by means of dividingand repeatedly dividing the polygon into two subsets>5505米until the simplest polygon subsets i.e. the boundarytriangles, are yielded. Therefore, if there are any in-ternal boundaries representing for example holes,Fig.1 The cut and patch procedure of generating sin-then connector lines must be introduced connectinggular elements:中国煤化工fYHCNMH G.230Alshoaibi et al. 1 J Zhejiang Univ SciA 2007 8(2);228-236Fig.2 shows example geometry where the wholeThe adaptive mesh refinement is based on aprocess of generating the mesh is ilustrated for better posteriori error estimator which is obtained from theunderstanding. Fig.2a ilustrates the geometry of a solution from the previous mesh. The error estimatorplate with six holes and two notches. Fig.2b shows six used in this paper is based on stress error norm. Theconnector lines forcing the internal boundaries to be strategy used to refine the mesh during analysisthe continuous part of the external boundary. Fig.2c process is adopted from (Ariffin, 1995) as follows:shows the cutting out of the rosette templates around(1) Determine the error norm for each elementeach crack tip. The background mesh for this domainis then set up automatically using dichotomy tech-|EIll= j(σ-σ°)°(σ_σ”)dS， (1)nique as shown in Fig.2d. Fig.2e shows the conven-tional mesh being generated by the advancing frontmethod. The first generation produces mesh with where E is the error estimator, e is the current element,initial size set by user. Later, during adaptive re- Le is the area of the current element, σ is the stressfinement, this first generated mesh will be taken as the field obtained from the finite element calculation andbackground mesh. In Fig.2f, for each rosette template,σ is the smoothed stress field.quarter-point elements are then constructed. Fig.2g(2) Determine the average error norm over theshows the enlargement of the quarter-point element at whole domainone of the crack tip.|1==之jσ°σd2,(2)m三n。[0° |[98 |[98where m is the total number of elements in the wholedomain.(3) Determine a variable ce for each element asεe=.1 (|EI|[)/2(3)| 0%18oη (I|II)4/2(a(b)(c)d)where η is a percentage that measures the permissibleerror for each element. If c2>1 the size of the elementis reduced and vice versa.(4) The new element size is determined ash。=he /(eo)",(4)where he is the old element size and p is the order ofthe interpolation shape function.e)(f(gFig.2 (a)~(g) showing the mesh generation stages withSTRESS INTENSITY FACTOR AND CRACKinclusion of quarter point elementsPROPAGATIONIn general, smaller mesh size gives more accu-In this paper, the displacement extrapolationrate finite element approximate solution. However, method (Phongthanapanich and Dechaumphai, 2004)reduction in the mesh size leads to greater computa- is used to calculate the stress intensity factors as fol-tional effort.lows:中国煤化工YHCNMH G ..Alshoaibi et al. 1 J Zhejiang Univ SciA 2007 8(2);228-236231E2π「40/;-1*)-(v.-")](5)The direction normal to the maximum tangential3(1+v)(1+x)V L2stress can be obtained by solving dσ'dθ=0 for 0. Thenontrivial solution is given byKu=2π44(u,"-u)-(u'-u')一，(6K, sinθ+ K.(3cosθ-1)=0,(8)where E is the modulus of elasticity, vis the Poisson'swhich can be solved as:ratio, K is the elastic parameter defined by[(3- 4v)plane stress,|3K}+K、K?+8K2K=l(3-4v)/(1+v)plane strain,θ=土cos.K? +9K}and L is the quarter-point element length. Theu' and In order to ensure that the opening stress; associatedv are the displacement components in thex' andy’with the crack direction of the crack extension isdirections, respectively; the subscripts indicate theirmaximum, the sign of 6 should be opposite to theposition as shown in Fig.3.sign of Kr (Andersen, 1998).The criterion for crack to propagate from crack.y▲tip is based on the material toughness K。. If the cal-culated stress intensity factor, K2Kc then the crackwill propagate in the direction 0o expressed by Eq.(9).+u,xu',x'The crack increment length△a is taken as 10%~ -20%of the initial crack length a, inversely proportional tothe ratio of Km/K. The ratio represents the mixedmode proportionality, therefore shorter incrementlength should be taken to carefully justify the crackpath curvature when K1 is relatively large comparedL/4 3L/4to K (Bittencourt et al., 1996).fig.3 The arrangement of the natural quarter-pointtriangular elements around the crack tipNUMERICAL ANALYSIS AND RESULTSIn order to simulate crack propagation underIn order to carry comprehensively evaluate thelinear elastic condition, the crack path direction muststress intensity factors approximated by the devel-be determined.There are several methods used to predict theoped program, four well-known plate geometries,compact tension specimen, three-points bendingdirection of crack trajectory such as the maximumspecimen, central cracked plate and double edgecircumferential stress theory, the maximum energynotched plate are considered.release rate theory and the minimum strain energydensity theory.The maximum circumferential stress theoryCompact tension specimenhe compact tension test specimen geometrystates that, for isotropic materials under mixed-modeand the final adaptive mesh are shown in Fig.4. Theloading, the crack will propagate in a direction normalspecimen has an initial crack length a=9 cm, widthto maximum tangential tensile stress. In polar coor-W=18.8 cm, and the thickness B, has various values asdinates, the tangential stress is given byexplained below .The analytical stress intensity factor for thisσ,=cos→| K; cos2θ_Kμ sin0. (7) geometry can be calculated from (Anderson, 1994) asV2πr”2|22follows:中国煤化工MHCNMH G.32Alshoaibl et al. 1J Zhejiang Univ Sc/A 2007 82:2232640●Present study, Theon studya-9cn! Experiment3022.6 cm20W= 18.8 cm .0t10 12 14(aLoad (kN)。Present study: ExpermentExperiment48c)(b)Fig.4 Compact tension geometry (a), the final adaptiveFig.5 Stress intensity fator values for steel specimen.mesh (b) and the zooming of mesh (c) around crack tip(a) h=8.3 mm; (b) h=13.6 mmthe stress intensity factors results.K=P| 2+号| 0.886+4.64-13.32|wThe present values of stress intensity factors arevery close to the theoretical solutions and are in good+14.72()-5.6(agreement with the experimental results.()][=[(-品门]Fig.6 shows four steps of crack propagation. Thepredicted crack propagation seems to follow the modewhere P is the applied load.I trijctory very well.The computed values of the stress intensity fac-tor under plane stress condition are compared with theThree-point bend specimenexperimental and numerical results obtained fromThe gcometry of the three point bend specimen(Parnas et al, 1996) as shown in Fig.s. In their study,and the final adaptive mesh are shown in Fig.7. Thesteel plates, with modulus of elasticity 210 GPa and analytical stress intensity factor for this problem canPoisson's ratio of 0.3, are used. Steel specimens withbe calculated from (Broek, 1986) as:thickness of 8.3 and 13.6 mm are considered. Theexperimental results are compared with finite elementPS( a\3/+21.8| ., >5/2results using ANSYS software (Fig.5a). The com-|20(别)(w)parison also comprises the analytical solution (Fig.5).The results of the steel specimens show obviously the( a9/27-37.6effect of the variation of the two thicknesses on中国煤化工TYHCNMH GAlshoaibi et al. 1 J Zhejiang Univ SciA 2007 8(2);228-236233those obtained from (Freese and Baratta, 2006; Orang,1988; Fett, 1999) are shown in Table 1 for S/W=4.The dimensionless form of the estimated stress in-tensity factor is obtained by Eq.(12). Freese andBaratta (2006) provided values for a range of0