Cast3m implementation of the extended finite element method for cohesive crack

Vietnam Journal of Mechanics, VAST, Vol. 33, No. 1 (2011), pp. 55 – 64 CAST3M IMPLEMENTATION OF THE EXTENDED FINITE ELEMENT METHOD FOR COHESIVE CRACK Nguyen Truong Giang, Ngo Huong Nhu Institute of Mechanics, VAST Abstract. In this paper, the finite element for cohesive crack for quasi-brittle materials is constructed by XFEM, which allows incorporation of the displacement discontinuities in the element. The algorithm of construction and procedures for involving this finite element into code Cast3M are presented. The numerical calculations in fracture mechan- ics are presented to demonstrate the benefits of the proposed implementation. Key words: Fracture, cohesive crack, quasi-brittle material, XFEM, Cast3M. 1. INTRODUCTION The extended finite element method (XFEM ) was first developed for two-dimensional linear elastic fracture mechanics [1, 2, 3], then it was used for crack growth, intersecting cracks and fluid mechanics. In the recent years, XFEM has been also applied to the cohe- sive crack model [4], [6]-[9]. The model which Hillerborg introduced in [5] has been used extensively in non-linear fracture mechanics of quasi-brittle materials and is applied in this paper. There are other developments of XFEM to existing codes. An implementation of the extended finite element method for linear elastic fracture mechanics problem within the finite element software Abaqus is introduced in [11]. In the same interest, an object enriched finite element library, named OpenXFEM++ which has been built and used to solve various 2D linear elastic fracture problems with great success, is presented in [12]. The finite element for the brittle crack by XFEM is constructed and involved into code Cast3M [3]. In order to produce extendable code, the objective of this paper is to present the algorithm of construction of new extended cohesive crack element 4 nodes and the way to add it into an existing finite element code Cast3M. 2. COHESIVE CRACK FOR QUASI-BRITTLE MATERIALS The fracture process in quasi-brittle materials like concrete, rocks, ceramics, fiber reinforced composites etc., is characterized by the contribution of the resistance from aggregate interlocking, and another inelastic effect. At the crack tip region there exists a narrow strip, where occur non-linear phenomena, such as plasticity yielding, micro- cracks, etc., however the rest of the body presents an elastic behaviour (Fig. 1a). These inelastic located processes can be represented by cohesive models, which simulate the damage region close to the crack tip. These models use to represent those non-linear phenomena such 56 Nguyen Truong Giang, Ngo Huong Nhu as the "cohesive forces", which are applied along the "cohesive crack". These forces are defined by relations between stresses and displacements, the material resistance decreases proportionally to the increase of the displacements [4]. Fig. 1. Sketch showing the process zone of a cohesive law with smooth crack closure Cohesive stress t(w) used by Hillerborg [5] The stress closing or cohesive stress t(w) is applied to the crack surface. These cohesive stresses depend on the crack opening w, and vary from zero at a characteristic crack opening wc, to the tensile strength of the material ft at the tip of the crack (Figure. 1b). It is assumed that the cohesive stresses close the crack smoothly, thus even when the un-cracked material is considered as linear elastic, and often used in the modelling of cementitious materials, stresses are finite in the un-cracked material at the crack tip. In other words, during Mode I crack propagation, the stress intensity factor KI is zero and the condition for crack propagation is that the stress at the crack tip has reached the tensile strength ft. The zone in which the cohesive stresses are located was originally called the micro-cracked zone [5], and later the process zone or fracture zone [6]. Mathematical formulation for cohesive crack problem Consider a domain Ω as shown in Fig. 2, containing a crack Γd. The part of the crack where a cohesive law acted is denoted by Γcoh. Prescribed tractions F are imposed on the boundary ΓF whereas prescribed displacements u¯ (assumed to be zero for simplicity) are imposed on Γu. The stress field inside the domain, σ is related to the external loading F and the tractions t+, t− in the cohesive zone through the equilibrium equations. The equilibrium conditions for the quasi-brittle materials having cohesive crack without body force are [7]: ∇ · σ = 0, σ · n = F on ΓF , σ ·m = −t+ on Γ+ coh ,σ ·m = t− on Γ− coh , (1) m,n are the normal vectors for ΓF and Γcoh, respectively. The equilibrium condition across surface Γcoh is : t ≡ t− = −t+ Cast3M implementation of the extended finite element method for cohesive crack 57 G d n F Fig. 2. Geometry of modelled domain and notation The displacement discontinuity w across Γd can be expressed in terms of the dis- placement vector u computed on two sides of the discontinuity: w = u ∣∣∣Γ+ d − u ∣∣∣Γ− d Let U be the space of admissible displacements u in Ω; i.e. such that u = u¯ on Γu, u possibly discontinuous on Γd and u ∈ C0 (C0 is the space of Sobolev) in Ω \ Γd. By introducing the test functions v ∈ U0 (the displacements are equal zero on Γu), the weak form of the equilibrium equations reads [7]:∫ Ω σ (u) : ε (v) dΩ = ∫ ΓF F · vdΓ + ∫ Γ + d t+ · vdΓ+ ∫ Γ − d t− · vdΓ ∀v ∈ U0∫ Ω σ (u) : ε (v) dΩ = ∫ ΓF F · vdΓ− ∫ Γd t · wdΓ∫ Ω σ (u) : ε (v) dΩ+ ∫ Γd t · wdΓ = ∫ ΓF F · vdΓ (2) The equation (2) is weak form of (1) and it includes cohesive stresses t on the two crack faces. 3. THE EXTENDED FINITE ELEMENT METHOD Partition of unity A collection of functions ϕi(x), each belonging to a node, defined over a body Ω(x ∈ Ω) forms a partition of unity [1] if: n∑ i=1 ϕi (x) = 1 (3) 58 Nguyen Truong Giang, Ngo Huong Nhu n is the number of nodal points. The displacement approximation using partition of unity can be formed by: uh (x) = n∑ i=1 ϕi (x) ( m∑ α=1 ψα (x) a α i ) (4) where ψα are enrichment functions, m is the number of enrichment functions. Extended finite element method Finite element shape functions Ni are also a partition of unity because n∑ i=1 Ni (x) = 1 Therefore, finite element shape functions are chosen as the partition of unity func- tions. From equation (4), we note that the finite element space (ψ1 = 1, ψα = 0(α 6= 1)) is a subspace of the enriched space uh. The enriched displacement approximation is written as [2] uh (x) = ∑ I∈N NI (x)  uI +H (x) aI︸ ︷︷ ︸ I∈NΓ + 4∑ α=1 Fα (x) b α I︸ ︷︷ ︸ I∈NΛ   (5) where uI is the nodal displacement vector associated with the continuous part of the finite element, aI is the nodal enriched degree of freedom vector associated with the Heaviside function, and bαI is the nodal enriched degree vector associated with the elastic asymptotic crack tip function. N is the set of all nodes in the mesh; NΓ is the set of nodes whose shape function support is cut by the crack interior; and NΛ is the set of nodes whose shape function support is cut by the crack tip (NΓ ∩NΛ = ∅) (see [3]). The interior of crack is modeled by the generalized Heaviside enrichment function H , where H takes the value +1 above the crack and -1 below the crack [2] H (x, y) = { 1 if y > 0 −1 if y < 0 (6) The functions Fα are defined as: {Fα (r, θ)} = {√ r sin ( θ 2 ) , √ r cos ( θ 2 ) , √ r sin ( θ 2 ) cos θ, √ r cos ( θ 2 ) cos θ } (7) where (r, θ) are the local polar coordinates at the crack tip. Note that the first function in (7), √ r sin (θ/2), is the discontinuous across the crack faces whereas the last three functions are continuous. 4. DISCRETISED EQUATIONS AND FINITE ELEMENT FOR THE COHESIVE CRACK The numerical simulation of the cohesive crack is hereby performed within the ex- tended finite element method which is presented in the previous section. The displacement Cast3M implementation of the extended finite element method for cohesive crack 59 field u of the body can be additively decomposed into a continuous part ucont and a dis- continuous part udisc: u(x) = ucont(x) + udisc(x). (8) Inserting the Heaviside function into the discontinuous part, the displacement field in element where the extra degrees of freedom are acted can be written as u(x) = N (x)a+HN (x)b, vectors a has the standard degrees and b has the extra degrees. The displacement discontinuity across Γd can be given w = 2N (x)b (9) The strain field and stress field in element can be expressed ε = Ba+HBb σ =D(Ba+HBb) (10) where D relates the stress and strain. If T is called the relation between the traction and crack displacement, the traction can be written as t = T .w or { tn ts } = [T ] { wn ws } (11) Here, we used the cohesive law T which proposed by Hillerborg [5], the normal traction force tn transmitted across a discontinuity is made by an exponentially decaying function of the history parameter: tn = ft exp ( − ft Gf k ) , where ft is the tensile strength of the material and Gf is the fracture energy and k (internal variable) is the largest value of the normal separation wn. The crack shear stiffness is also made by a function of the history parameter. The shear traction ts acting on the discontinuity surface is calculated from ts = dint exp (hsk)ws, where dint is the initial crack shear stiffness (k= 0), ws is the crack sliding displacement and hs is equal to: hs = ln (dk=1.0/dint) where dk=1.0 is crack shear stiffness when k = 1.0. Replace the displacement field (9), (10), (11) into the equation (2), the discrete equations can be obtained as follow:[ ∫ Ω BTDBdΩ ∫ Ω HBTDBdΩ∫ Ω HBTBdΩ ∫ Ω HBTDBdΩ ]{ a b } + [ 0 0 0 4 ∫ Γd NTTNdΓ ]{ a b } = ∫ ΓF NTFdΓ 60 Nguyen Truong Giang, Ngo Huong Nhu The discretised form of virtual work for incremental displacements (da, db) is written as:[ ∫ Ω BTDBdΩ ∫ Ω HBTDBdΩ∫ Ω HBTBdΩ ∫ Ω HBTDBdΩ + 4 ∫ Γd NTTNdΓ ]{ da db } = { fext 0 } − { f inta f intb } (12) So, the element matrix can be written: Ke = [ ∫ Ω BTDBdΩ ∫ Ω HBTDBdΩ∫ Ω HBTBdΩ ∫ Ω HBTDBdΩ+ 4 ∫ Γd NTTNdΓ ] , (13) where the matrices of shape functions derivatives B for node i are: B=i   Ni,x 00 Ni,y Ni,y Ni,x   (14) f ext = ∫ ΓF NTFdΓ (15) f inta = ∫ Ω BTσdΩ, f intb = ∫ Ω HBTσdΩ+ 2 ∫ Γd NT tdΓ (16) Numerical integration for the element with the cohesive crack As in the standard FEM, it is necessary to perform numerical integration over the element domain for computing the element stiffness matrix (13). The elements which contain the crack include a displacement discontinuity due to the X-FEM formulation. These elements must be subdivided into sub-domains in which the crack is one of the sub-domain boundaries to carry out the numerical integrations. Therefore, the elements non-crossed by a discontinuity are integrated by standard four-point Gauss quadrature. When an element is crossed by a discontinuity, the domains Ω+ and Ω− on either side of the discontinuity are triangulated into sub-domains [3]. Within each triangular sub- domain, three-point Gauss quadrature is applied. For the cohesive crack sub-domain, two integration points are positioned on the discontinuity in order to integrate the traction forces. 5. CAST3M IMPLEMENTATION FOR THE COHESIVE ELEMENT In this part, a procedure for the implementation of the XFEM within the FE code Cast3M [10] for two-dimensional fracture problems is presented. The quadrangular element 4 nodes CRCOH is constructed and involved into Cast3M. The implementation was based on the use of this element and performed in source code of Cast3M by language FORTRAN and enables in modeling of different crack locations and orientations. In order to develop XFEM in numerical calculation, a number of subroutines contain XFEM formulation and some operators represented below such as RIGI for matrix stiffness K of cohesive element, EPSI for the displacement ε, COMP for stress σ, BSIG for nodal force,. . . are inserted into software Cast3M . Cast3M implementation of the extended finite element method for cohesive crack 61 5.1. Algorithm - Step 0: Initially, a mesh for quadrilateral element, the location of the crack and the parameters of concrete are inputted. - Step 1: Using level set to choose crack tip element and crack element. For the element without crack, the standard FEM is applied. - Step 2: Calculating of the matrix stiffness K for cohesive elements by using (13) and then assembling in the matrix stiffness standard. - Step 3: Applying a given loading step to solve for the displacement field { da db } using Eq. (12). - Step 4: Calculating the incremental deformation dε, and then determinating the stress σ. Note that the cohesive stresses are local stresses, leading to the co-ordinate system changing. - Step 5: Computing the nodal forces through the program BSIG using Eq. (15) and (16). - Step 6: Verifying the condition of convergence. If it isn’t satisfied, go to step 3. On the contrast, if the condition is satisfied, the numerical results as the stress, deformation, and displacement for each the loading step can be obtained. 5.2. Numerical example and discussion Example 1: Consider a plate of width w = 100mm and height L = 200mm with an edge crack of length a = 20mm. The tensile stress σ =10MPa. The concretes parameters Ec = 35000MPa, ν = 0.3, ft = 3MPa, Gf = 0.1N/mm. Fig. 3. Comparison of the stresses 62 Nguyen Truong Giang, Ngo Huong Nhu Fig. 4. Comparison of the deformations Remark: Separation of crack in the linear case (=3.677x10−2 mm) is larger than that of cohesive crack (=2.73086x10−2 mm). The stress at a crack-tip (= 45MPa) in cohesive crack is also represents a contribution of the cohesive stress in comparison to stress intensity (= 57 MPa) of the linear elastic crack. Example 2: We use the same dimension of specimen for the incline crack α = 45˚, length a = 20mm. Fig. 5. Comparison of the deformations and the σyy for incline crack Cast3M implementation of the extended finite element method for cohesive crack 63 Remark: In the incline crack case, the influence of the quasi-brittle material on the crack side displacements is presented (see Fig. 5). Clearly, the stress at a crack-tip of cohe- sive crack (σmax = 26 MPa) are smaller than that of linear elastic crack (σmax = 47 MPa). 6. CONCLUSIONS In this paper, a procedure for the implementation of the XFEM within the FE code Cast3M for two-dimensional fracture problems is presented. The quadrangular finite element CRCOH for cohesive crack are constructed and involved in to Cast3M for quasi- brittle materials. This finite element is enriched by the cohesive stress which is related with a separation of the crack sides through the law of Hillerborg. The obtained results are compared with the results of the linear elastic crack showing that the existence of cohe- sive stress in quasi-brittle materials reduces the separation of crack and stress at crack-tip. So, using the above mentioned model with cohesive crack and given finite element CRCOH is necessary. In the near future, relied upon this framework, the growth of cohesive crack will be continued to develop. ACKNOWLEDGMENTS The authors are grateful to Professor Alain Millard from French Atomic energy Com- mission (CEA) for his assistance. This paper is completed with the partly financial support from the NAFOSTED (National Foundation for Science and TechnologyDevelopment). REFERENCES [1] Melenk JM, Babuska I, The partition of unity finite element method: Basis theory and appli- cation, Computer Methods in Applied Mechanics and Engineering, 45(5) (1999), 601-620. [2] Moes N, Dolbow J, Belytschko T, A finite element method for crack growth without remeshing, International Journal for Numerical Methods in Engineering, 46 (1999), 131-150. [3] Nguyen Truong Giang, Numerical simulation the crack 2D by the finite element incorporated the discontinuity, Vietnam Journal of Mechanics, VAST, 30(2) (2008). [4] Estela M.R. Bueno, Tulio N. Bittencourt, Luiz F. Martha, Carlos V.A. Carvalho, A graphics and object-oriented system for modeling 2D cohesive crack problems, European Congress on Computational Methods in applied Sciences and Engineering, (2000). [5] Hillerborg, A., Modéer, M. and Petersson, P.E., Analysis of crack formation and crack growth concrete by means of fracture mechanics and finite elements, Cements and Concrete Research, 6 (1976), 773-782. [6] Henrik Stang, John Forbes Olesen, Peter Noe Poulsen, Lars Dick Nielsen, On the application of cohesive crack modeling in cementious materials, Journal of Materials and Structures, 40 (2007), 365-374. [7] Stefano Mariani, Umberto Perego, Extended finite element method for quasi-brittle fracture, International Journal for Numerical Methods in Engineering, 58 (2003), 103-126. [8] Belytschko T, Black T, Elastic crack growth in finite elements with minimal remeshing, In- ternational Journal for Numerical Methods in Engineering, 45 (1999), 601-620. 64 Nguyen Truong Giang, Ngo Huong Nhu [9] Sukumar N, Prevost JH, Modeling quasi-static crack growth with the extended finite element method, Part I: Computer implementation, International Journal of Solids and Structure, 40 (2003), 7513-7537. [10] [11] Giner E, Sukumar N, Tarancon JE, Fuenmayor FJ, An Abaqus implementation of the extended finite element method, Engineering fracture mechanics, 76 (2009), 347-368. [12] Bordas S, Phu Vinh Nguyen, Dunant C, Hung Nguyen-Dang, Guidoum A, An extended finite element library, International Journal for Numerical Methods in Engineering, 2 (2006), 1-33. Received March 29, 2010