疲劳裂纹扩展

Extended finite element methods independent of the mesh structure but determined by the mixed-mode loading cases. d mo rial fa tive f ward, based on the partition of unity approach, Belytschko and his co-workers [10,11] presented the extended finite element concept in which the crack can propagate through elements in an arbitrary manner. Wells and Sluys [12] combined this method with the cohesive zone model to investigate fracture of concrete materials and obtained excellent mixed-mode crack prediction compared with experiments. In addition, another method of incorporating a discontinuity in the FEM has been reported by Hansbo and Hansbo [13]. Mergheim et al. [14] applied this method with cohesive zone model to study crack 0013-7944/$ - see front matter � 2008 Elsevier Ltd. All rights reserved. * Corresponding author. Tel.: +49 202 439 2124; fax: +49 202 439 2027. E-mail address: h.yuan@uni-wuppertal.de (H. Yuan). Engineering Fracture Mechanics 76 (2009) 165–181 Contents lists available at ScienceDirect Engineering Fracture Mechanics doi:10.1016/j.engfracmech.2008.08.011 have been further put forward and validated to be suitable for simulating crack nucleation as well as crack propagation in both brittle and ductile materials. However, so far the development of CCZM still stands at the beginning. As known, for engi- neering problems, the components are mostly subjected to mixed-mode loading, in which the crack usually advances along an unknown curved path. Actually, the curved crack propagation is still difficult to be simulated by CZM [6], and even impracticable for 3D mixed-mode loading [7]. The primary reason is that the method requires pre-assumption of the crack path in the FE models beforehand, and crack only can propagate along the boundaries of elements. Therefore, application of the cohesive zone models are mainly limited to mode I loading cases. Recently, several promising methods of modelling solid discontinuities induced by, e.g., cracking have been developed since Melenk and Babüska [8] proposed the approach of partition of unity. The idea was exploited firstly by Oden and his co-workers [9] for solving the internal boundary problem which is referred to as generalized finite element (GFEM). After- Mixed-mode loading Fatigue crack growth Cyclic cohesive model 1. Introduction Since the damage mechanics base has been proposed to deal with mate recently by considering the accumula Numerical simulations illustrate the ability of this method to simulate fracture with unstructured meshes. The computational results agree with known fracture experiment data. Known fatigue observations can be predicted using the present model. � 2008 Elsevier Ltd. All rights reserved. delling confronts difficulties of strain localizations, the cohesive zone model (CZM) ilure [1,2], which provides an alternative way to describe crack propagation. More atigue damage in the cohesive zone, the cyclic cohesive zone models (CCZM) [3–5] Computational analysis of mixed-mode fatigue crack growth in quasi-brittle materials using extended finite element methods Yangjian Xu a,b, Huang Yuan a,* aUniversity of Wuppertal, Department of Mechanical Engineering, Gaussstr. 20, 42119 Wuppertal, Germany b Zhejiang University of Technology, College of Mechanical and Electrical Engineering, Hangzhou, China a r t i c l e i n f o Article history: Received 29 November 2007 Received in revised form 13 August 2008 Accepted 25 August 2008 Available online 6 September 2008 Keywords: a b s t r a c t In the present paper an extended finite element method (XFEM) containing strong discon- tinuity within elements is introduced and implemented in the commercial general purpose software ABAQUS. The algorithm allows introducing a new crack surface at arbitrary loca- tions and directions in elements. To consider fatigue crack nucleation and propagation in quasi-brittle materials the XFEM is combined with a cyclic cohesive model. Accumulative material damage is described by separate evolution equations. The crack path is completely journal homepage: www.elsevier .com/locate /engfracmech Nomenclature ai conventional degrees of freedom on node i a conventional nodal degrees of freedom a^ conventional degrees of freedom Da crack growth length Ai area related to Gauss point i bij enhanced degrees of freedom associated with node i and the jth term of enriched user function b enriched nodal degrees of freedom ~b enhanced degrees of freedom B matrix containing derivatives of shape function N D constitutive tangent matrix dn normal separation at discontinuity dt tangent separation at discontinuity d0 cohesive length dn,max maximum value of normal separation before unloading du accumulated cohesive length E Young’s modulus f0 cohesive endurance limit ft0 initial maximum traction of material fi,ext equivalent nodal exterior force related to degree of freedom i fi,int equivalent nodal interior force related to degree of freedom i ft material tensile strength in cohesive law ui standard basis function for node i G fracture energy release rate Gt shear modulus of the cohesive zone Cd discontinuous border Ct prescribed traction border Cu prescribed displacement border H height of specimen H0 Heaviside step function for damage accumulation HCd Heaviside step function e strain tensor j separation related to the normal separation dn Kij stiffness matrix related to i j degrees of freedom l characteristic length L length of specimen m0 inward unit normal vector which is perpendicular to discontinuity Cd m term number of enriched user function. n0 outward vector normal to body n number of nodes ngp Gauss points N shape function r gradient operator X+/� positive and negative part of FEM domain x damage indicator of material x0 threshold of damage indicator P force Pmax maximal force wj enriched user function of the jth term (with j = 0,1,. . ..,m) for node i r distance from Gauss point to cohesive zone tip R stress ratio t inner traction at discontinuity �t prescribed traction vector at boundary Ct T tangent matrix for the discontinuity Tn normal traction at discontinuity Tn,max normal traction corresponding to dn,max Tt shear traction at discontinuity u vector of displacement [u] displacement jump across the discontinuity �u prescribed displacement at boundary Cu u^ standard term of displacement u ~u enhanced term of displacement u 166 Y. Xu, H. Yuan / Engineering Fracture Mechanics 76 (2009) 165–181 propa ation w Y. Xu, H. Yuan / Engineering Fracture Mechanics 76 (2009) 165–181 167 ments [24,27]. Moreover, it allows multi-cracks nucleation, growth and coalescence without remeshing. In this class of methods, the discontinuity is embedded by means of introducing additional degrees of freedom on those elements where the discontinuity crosses. The special appealing feature is that it inherits the finite element framework and its property such as sparsity and symmetry, which provides a robust and efficient method to model the strong discontinuity. A comprehensive review and comparisons among various discontinuity approaches have been provided by Jirasek et al. [16]. Additionally, it is worthy of note that the cohesive zone model can be naturally integrated into XFEM, by which one may predict crack nucle- ation and propagation under static and dynamical conditions [11,15]. In quasi-brittle materials failure is characterized by damage within narrow zone caused by high stress concentrations [17]. The CZM was first applied to predict crack propagation in the 1970s [18]. Most known works are, however, limited in simple loading configurations or for monotonic loading conditions [23,25]. Under mixed-mode loading conditions failure prediction of the specimen is getting more complicated. From computational point of view this problem is coupled with effi- cient numerical algorithms as well as modelling of strong discontinuities in continuum solid. In the present work the CCZM has been implemented into the extended finite element formulation in conjunction with the general commercial FEM package ABAQUS [19]. After description of computational algorithms and XFEM implementa- tion we have verified the programming using several known numerical examples under monotonic as well as cyclic loading conditions. Since this method is incorporated within the package of ABAQUS the robust solution ability and other features such as contact algorithms in ABAQUS can be retained. Hereby, the treatments of forming stiffness matrix, enriching degrees of freedom, introducing discontinuity into the intersected element and visualizing the results have been illustrated. 2. The method of partition of unity 2.1. Discontinuity in finite elements A comprehensive overview of development of the discontinuity method in conjunction with partition of unity property has been given by numerous publications, e.g., in [16,20]. The XFEM is an efficient way to reduce mesh dependencies for crack growth analysis. The key lies in the fact that the discontinuity field may be directly approximated by the standard degrees of freedom together with the enriched ones, meanwhile the specific user function wj can be incorporated with the standard basis function ui to describe the character of discontinuity field. The product, uiwj, not only inherits the good approximate property from standard basis function ui, but also the specific discontinuous property from user function wj. The approximation is only locally enriched in some region of interest such as in the damage zone. According to the approach of partition of unity, if the standard basis function ui fulfils the partition of unity (i.e., Pn i¼1ui ¼ 1, the enriched field approximation can be expressed as follows: u ¼ Xn ui ai þ Xm wjbij ! ; ð1Þ where terms hance the m In riched the fin where is defi gation under mixed-mode conditions. Using extended finite element methods (XFEM), one can simulate crack nucle- at an arbitrary material point and crack propagation in an arbitrary direction without introducing extra nodes and ele- w^ standard term of test function w ~w enhanced term of test function w Wi weight function related to Gauss point i r Cauchy stress tensor r1 first principal stress ru cohesive strength r^ nonlocal stress test function Du load amplitude m Poisson’s ratio i¼1 j¼1 n is the number of the nodes associated with the standard basis function uiwj is the enriched user function with m for the given node i. ai is the vector of the conventional degrees of freedom on node i, and bij is the vector of the en- d degrees of freedom associated with node i. Since degrees of freedom are enriched on nodes of the existing element, odification of topology can be avoided. the finite element notation, the element shape function N fulfils the requirement of the partition of unity. If the en- user function takes the Heaviside step function HCd , the approximation form of discontinuity displacement field in ite element frame can be written as u ¼ Nðaþ HCdbÞ ¼ Naþ HCdNb; ð2Þ a and b represent the standard and enriched nodal degrees of freedom, respectively. The Heaviside step function HCd , ned as tion, o jump Co ven as to bod placem Th eleme equat Herein space The n ment That i 168 Y. Xu, H. Yuan / Engineering Fracture Mechanics 76 (2009) 165–181 0 d 0 y. Correspondingly, t is the inner traction and �t is the prescribed traction vector at boundary Ct. �u is a prescribed dis- ent at boundary Cu. e foregoing strong form of governing equation will be recast into the weak form, since the solution obtained by finite r � r ¼ 0 ð5Þ where r is the gradient operator and r is the Cauchy stress tensor. The natural and essential boundary conditions are given as r �m0 ¼ t on Cd; r � n0 ¼ �t on Ct; u ¼ �u on Cu: ð6Þ Herein m is the inward unit normal vector, which is perpendicular to discontinuity C and n is the outward vector normal nsidering the domain occupied by a body with a discontinuity Cd as shown in Fig. 1, the equilibrium equations are gi- in Eq. (5), without taking body forces into account. ½u� ¼ Nbjx¼Cd : ð4Þ In nonlinear analysis of the finite element framework combined with discontinuity, the nodal displacement is computed iteratively from the global equilibrium equations. The XFEM can be applied to solve nonlinear problems directly. 2.2. XFEM governing equations across the discontinuity which is expressed as (2) will be considered in these nodes. Otherwise, the standard shape function is retrieved. By the enrichment of degrees of freedom, the discontinuity is introduced straight forward, without modifying finite element mesh. Effectively, the standard nodal degrees of freedom, a, are used to represent the continuum part, while the enriched DOFs, b, stand for displacement d nly the nodes related to the elements containing the discontinuity are enriched. That is, the extra enriched item in Eq. HCd ¼ 1 x 2 Xþ 0 x 2 X� ( : ð3Þ The whole FEM domain consists of part X+ and X� with the discontinuous border C as shown in Fig. 1. In the formula- Fig. 1. The FEM domain crossed by the discontinuity Cd and the applied boundary condition. Before the material failure starts Cd does not exist. nt is the approximate solution to the strong form. The weak form can be obtained by taking product of the equilibrium ion expressed in the form of trial function and the corresponding test function in accordance with the Galerkin method. , the displacement spatial trial function is expressed as Eq. (2). At the same time, the test function is taken the same as displacement field u, which is written as w ¼ w^þ HCd ~w: ð7Þ otations (^) and (�) denote the standard and enhanced term in approximations, respectively. The extended finite ele- equations are obtained using the Galerkin method from the equilibrium Eq. (5) asZ X r � rð Þ �w dV ¼ 0: sZ X ðr � rÞ � w^þ HCd ~w � � dV ¼ 0: ð8Þ Th rizing Th 2.3. D In o where per, th can be hand Y. Xu, H. Yuan / Engineering Fracture Mechanics 76 (2009) 165–181 169 KTab Kbb Db ¼ fb;ext � fb;int ð15Þ side of (13) cannot be performed without linearization. Finally the finite element equation can be written as [27] Kaa Kab � � Da � � fa;ext � � fa;int � � by this expression, the cohesive law can be incorporated for describing connections of solids X� and X+ along Cd. The cohe- sive law is nonlinear and will be discussed separately. Due to the nonlinear behaviour at Cd, the second integral in the left generally written as t ¼ T½u� ¼ TðNbÞ at Cd; ð14Þ Xþ B r dV þ Cd N t dC ¼ Cþt N �t dC: Eq. (13) is the discretized XFEM governing equations. The relation between the traction and the displacement jump in Cd X BTr dV ¼ Ct NT�t dC;Z T Z T Z T ð13Þ Cd Following the Bubnov–Galerkin procedure, Eq. (10) can be discretized asZ Z rSw^ ¼ Ba^ rS ~w ¼ B~b rS means the symmetric term of the gradient operator numerically calculated through the matrix B [28]. In this pa- e elastic materials are considered. r ¼ De ¼ DrSu ¼ DðBaþ H BbÞ linearization of Eq. (10) should be carried out. The nonlinearity is induced due to cohesive model. The resulting nonlinear discrete weak form is consistently solved by Newton–Raphson scheme. We use the conventional shape function, N, in the isoparametric element formulation [28] as u^ ¼ Na ~u ¼ Nb w^ ¼ Na^ ~w ¼ N~b ð11Þ where the displacement u is decomposed into two parts as the weight function w, u ¼ u^þ HCd ~u. Using the discretized inter- polation (11) the gradients can be expressed as rSu^ ¼ Ba rS~u ¼ Bb ð12Þ iscretization and linearization rder to apply the weak form of governing Eq. (10) in the framework of finite element, the discretized interpolation and tional function to describe behaviour of Cd. Generally speaking, Cd is determined during the computation and behaves non- linearly so that the whole finite element equations are to be solved incrementally. Xþ Cd C þ t e Eq. (10) forms the governing equation for the XFEM. To connect the region X+ and X� we have to introduce an addi- X ðrw^Þ : r dV ¼ Ct w^ � �t dC;Z ðr~wÞ : r dV þ Z ~w � t dC ¼ Z ~w � �t dC ð10Þ Z Z In accordance with the decomposition of the displacement field, the admissible variationw can be taken the same fashion of decomposition into an admissible variation to w^ and ~w. It follows that the set of coupled equations yields [27] d the results above the Eq. (8) turns intoZ X ðrw^Þ : r dV þ Z Xþ ðr ~wÞ : r dV þ Z Cd ~w � t dC ¼ Z Ct w^ � �t dCþ Z Cþt ~w � �t dC ð9Þ X ðr � rÞ � HCd ~w dV ¼ Cþt ~w � �t dC� Cd ~w � t dC� Xþ ðr ~wÞ : r dV : In the volume integration of the above equation, HC is eliminated by transferring integral region from X to X+. Summa- Z Z Z Z X r � rð Þ � w^ dV ¼ Ct w^ � t dC� X ðrw^Þ : r dV : The enhanced term in Eq. (8) makes the XFEM equation differ from the conventional FEM. Attention has to be paid to the discontinuity curve Cd, e integral with w^ in Eq. (8) leads to the governing equation in conventional FEM which reads [12]Z Z � Z with Kaa ¼ Z X BTDB dV Kab ¼ Z Xþ BTDB dV Kbb ¼ Z Xþ BTDB dV þ Z Cd NTTN dV fa;ext ¼ Z Ct NT�t dC f b;ext ¼ Z Cþt NT�t dC fa;int ¼ Z X BTr dV f b;int ¼ Z Xþ BTr dV þ Z Cd NTt dC In Eq. (15)Da and Db denote increments of nodal displacements in an incremental step. In the infinitesimal displacement formulation the stiffness matrices Kaa and Kab are independent of loading intensity. The nonlinearity of Kbb is caused by the nonlinear cohesive law which will be discussed in the following section. 3. Cohesive zone models (CZM) 3.1. The CZM for monotonic loading in quasi-brittle materials The CZM was introduced by Dugdale as well as by Barrenblatt in the 1960s to analyze cracks. The cohesive zone should take into account the key features of inelastic material ahead of the crack tip [1,2]. The constitutive relation there is replaced by a so-called cohesive law which connects traction(s) and separation(s) in the cohesive zone [22]. Under mixed-mode load- ing conditions both normal and tangent separations (dn, dt) are related to normal traction Tn and shear stress Tt, respectively. behaves linearly elastic. That is, a cohesive zone is not yielded. As soon as f is exceeded, the material begins to fail and be- haves soften Fig. 2. 170 Y. Xu, H. Yuan / Engineering Fracture Mechanics 76 (2009) 165–181 reloading conditions the traction decreases/increases linearly with separation. Material damage is characterized by diminishing the maximum tensile traction ft or/and decreasing the cohesive zone stiffness. 0 0.2 0.0 1.0 2.0 3.0 4.0 5.0 6.0 δn 0/δ Cohesive law for CCZM. In monotonic loading case the cohesive law in dn > 0 is identical with the model of Xu–Needleman [22]. For unloading and Tn ¼ ft jd0 exp � j d0 � d 2 t d20 ; ð16Þ where ft is the tensile strength of the material, which is the threshold for cohesive zone propagation. j is directly related to the normal separation dn, which reads j = d0 + dn. It means a linear coordinate shifting in the traction–separation diagram. For the shear separation, the present model does not coincide with Xu–Needleman suggestion, since the normal separation dominant failure happen in quasi-brittle materials 0.4 0.6 0.8 1 T n /f t0 Monotonic Loading Cyclic Loading (one cycle) Degradation after one cycle t unstable. For this reason we neglect the hardening curve in the Xu–Needleman model and just consider the material ing regime in the normal separation, as in Fig. 2. The normal traction–separation relationship can be written as ! The material damage is modelled using a softening traction–separation law. In our analysis we use the popular model intro- duced by Xu and Needleman [22] with certain modifications that suits quasi-brittle and cyclic failure. Brittle material fails without significant inelastic deformations. For this reason all knownmodels for quasi-brittle materials do not consider mate- rial hardening in their cohesive laws [12–14,16]. Before the first principal stress reaches the tensile strength, f