Numerical implementation of the eXtended Finite Element Method for dynamic crack analysis

Available online at www.sciencedirect.com

Advances in Engineering Software 39 (2008) 573–587 www.elsevier.com/locate/advengsoft

Numerical implementation of the eXtended Finite Element Method for dynamic crack analysis Ionel Nistor, Olivier Pantale´ *, Serge Caperaa L.G.P C.M.A.O – E.N.I.T, 47 Av d’Azereix BP 1629, 65016 Tarbes Cedex, France Received 18 September 2006; received in revised form 25 May 2007; accepted 8 June 2007 Available online 14 September 2007

Abstract A numerical implementation of the eXtended Finite Element Method (X-FEM) to analyze crack propagation in a structure under dynamic loading is presented in this paper. The arbitrary crack is treated by the X-FEM method without re-meshing but using an enrichment of the classical displacement-based finite element approximation in the framework of the partition of unity method. Several algorithms have been implemented, within an oriented object framework in C++, in the home made explicit FEM code. The new module, called DynaCrack, included in the dynamic FEM code DynELA, evaluates the crack geometry, the propagation of the crack and allow the post-processing of the numerical results. The module solves the system of discrete equations using an explicit integration scheme. Some numerical examples illustrating the main features and the computational efficiency of the DynaCrack module for dynamic crack propagation are presented in the last section of the paper.  2007 Elsevier Ltd. All rights reserved. Keywords: Partition of unity; eXtended Finite Element Method; Finite element programming; Dynamic crack propagation; Dynamic energy release rate

1. Introduction The development of computational techniques for the analysis of dynamic fracture and their implementation in numerical codes are becoming more and more important in recent years. Such interest is motivated by the desire to predict both the initiation of a crack and its propagation through the structure under dynamic loading. This is a typical case concerning impact applications where severe dynamic loading induces damage and fracture of the material. Several numerical approaches have been proposed in the last decades for analyzing some discontinuous phenomena, such as cracks and shear bands, occurring in structures under quasi-static or dynamic loading. The first category concerns the re-meshing methods that are usually used for modeling cracks or other strong discontinuities in structures. Based on classical finite element *

Corresponding author. E-mail address: [email protected] (O. Pantale´). URL: http://www.enit.fr (O. Pantale´).

0965-9978/$ - see front matter  2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.advengsoft.2007.06.003

method (FEM), the geometry is usually re-meshed at each time step during the discontinuity propagation. In the most recent developments, the re-meshing area has been limited to the immediate vicinity of the discontinuity to save computational time. Because of its simplicity (a standard FEM program and a re-meshing algorithm are sufficient to evaluate crack initiation and propagation), different versions of this technique have been implemented in commercial codes, especially for quasi-static analysis. Nevertheless, several important drawbacks remain. The mesh dependence of the crack is one of the main. The user must have ‘‘a priori’’ knowledge of the response of the model in order to generate an accurate initial mesh in the crack-tip region; beside that, the direction of the crack propagation is usually very sensitive with nodes alignment. Another important difficulty is the remapping of the data attached to physical points situated around the crack between the old mesh and the new one. For dynamic fracture problems, this approach remains quite difficult to apply. Discontinuity methods appeared as an innovating technique to model crack growth using cohesive segments at

574

I. Nistor et al. / Advances in Engineering Software 39 (2008) 573–587

the element’s interface without the necessity to introduce supplementary nodes. This approach was used by Xu and Needleman [1] and Camacho and Ortiz [2]. Its numerical implementation for an explicit time integration scheme was presented by Remmers and co-workers [3]. Based on the use of the partition of unity method (PUM) developed by Babuska and Melenk [4], a crack is represented by a number of overlapping cohesive segments which are inserted as discontinuities in the displacement field of the elements cut by the crack. In this approach, the crack direction is limited to the element edges and thus the crack paths are limited to specified directions [5]. Even based on substantial theoretical foundations, this approach is not still completely well-contained. For instance, the link between the parameters of the cohesive law used in the models with measurable material properties is not well known. The embedded discontinuity methods represent another class of fracture methods that consider the crack at the element level as a band of high strain. This approach was proposed by Belytschko, Fish and Engelmann [6]. A good description of these methods is given by Dvorkin et al. [7], Simo et al. [8] and references quoted. Jirasek [9] published a comparative study of these methods and showed that there is three major classes of models with embedded discontinuities. One of these classes is known as the statically and kinematically optimal non-symmetric (SKON) formulation [9] which allow to represent both the kinematic and the static aspects properly and leads to an improved numerical performance. It is able to effectively represent the complete separation at later stages of the fracturing process without the transfer of spurious stress. Recently Oliver and Huespe [10] presented a numerical implementation of this approach based on the use of finite elements with embedded discontinuities where both nodal and elemental enrichments are taken into account. Finally we present in this literature survey the eXtended Finite Element Method (X-FEM). The general idea of this method is to enrich the displacement approximation space spanned by standard finite element shape functions with some specific discontinuous functions. It is about an approach based on PUM. To our knowledge, Belytschko and Black [11] have been the first one to model a crack using this approach. The name X-FEM was given by Moes et al. [12] and Dolbow et al. [13] after the introduction of a step function enrichment in the displacement field for the elements entirely cut by the crack. This method we chosen for implementation in the home made DynaCrack code in order to carry out dynamic analysis of structures containing discontinuities. Concerning the restricted application field relative to dynamic crack analysis, a few contributions based on the X-FEM formulation have been published. The major contribution has been proposed by Belytschko et al. [5] with an adaptation of the X-FEM approach for dynamic crack analysis and the development of a new discontinuous enrichment. They proposed a crack evolution model based on the loss of hyperbolicity criterion. Re´thore´ et al. [14]

investigated some instability problems occurring in XFEM dynamic crack analysis. They proposed a technique called ‘‘balance recovery method’’, that provide the ability to evaluate both numerical stability and accuracy for any type of projection used with varying meshes. Recently, Menouillard et al. [15] proposed a new technique to evaluate the lumped mass matrix for the enriched elements in XFEM. This is a non-trivial result because of the additional degrees of freedom linked to the enriched nodes. This later allows to use an explicit integration scheme where the critical time step does not tend to zero when the crack is in the immediate vicinity of a node. The present paper is organized as follows. A description of the X-FEM with a focus on the explicit time integration scheme is presented in Section 2. The main features of the crack evolution model and its numerical implementation are presented in Section 3. The entire procedure used for the numerical integration in the DynELA FEM code [16] is presented in Section 4. Numerical simulations illustrating the robustness and the effectiveness of the implemented algorithms are presented in Section 5. A brief review of the problems encountered during this work and some conclusions and future works are reported in the last section of this paper. 2. eXtended Finite Element Method for dynamic crack analysis The crack representation in X-FEM is based on the enrichment of the classical displacement-based finite element approximation through the framework of the partition of unity method. Therefore, a crack is modeled by introducing a set of additional degrees of freedom to the nodes whose nodal shape function space intersects this crack. Within the approach proposed by Moes et al. [12], two types of nodal enrichments are considered in order to model the crack. When an element is splitted in two parts by the crack, all the nodes of this element are enriched by the Heaviside step function, while, the Westergaard asymptotic function is used to enrich the nodes of the elements containing the crack tips. Therefore, in a 2-D analysis for example, the nodes of fully-cut elements have two classical degrees of freedom and two enriched degrees of freedom modeling the strong displacement jump. The nodes of the elements containing the crack-tip have two classical degrees of freedom and four enriched degrees of freedom based on the radial and angular behavior of the asymptotic displacement field. As the crack propagates and the crack-tip crosses the edges of the elements, the enrichment status of some nodes of the structure changes. This type of enrichment was implemented by Sukumar et al. [17] for a quasi-static crack. Fig. 1a shows a finite element mesh where the circled nodes are enriched with the Heaviside step function, while the squared nodes are enriched with the Westergaard functions. As discussed by Belytschko et al. [5], this mixed enrichment is not easy to incorporate in methods with timedependent solutions. Therefore, they proposed recently

I. Nistor et al. / Advances in Engineering Software 39 (2008) 573–587

a

b

y

y

x

575

x

modeling with crack-tip functions

modeling with Heaviside functions

Fig. 1. Crack modeling using X-FEM approach.

another approach for the enrichment in dynamic applications. The basic idea is to avoid the near-tip enrichment by imposing the crack-tip to cross one element at a time. Essentially the crack-tip goes from edge to edge and the enrichment for such a situation can be treated using only Heaviside step functions. In Fig. 1b, we present for the same mesh and arbitrary crack the new enrichment and the new path for the crack resulting from this approach. The major consequence of limiting the crack-tip position at the element’s edge is that the modification of the direction will only occur at the element edges. The comparison of Fig. 1a and b shows that the mesh size plays a more important role, since a finer mesh helps to minimize the errors due to the crack path approximation and the fact that the crack tips are not taken into account. In the same time, being a simpler procedure for implementation in dynamic crack analysis, we have chosen to use it. This involves a special choice in the crack evolution models (as presented further in Section 3), in order to avoid that the propagation models require quantities not accurately computed within this framework.

( ~Þ ¼ H ðX

~  X~ Þ  ~ þ1; if ðX nP0 1; otherwise

ð1Þ

~ is the considered point in the initial configuration, where X  ~ ~ onto the crack and ~ n is the unit X is the projection of X outward normal to the crack at X~ . Each node I of the ~Þ. The new dismesh is associated to a shape function /I ðX ~ h ~ continuous displacement field u ðX Þ for a N nodes mesh, including NC enriched nodes, is therefore approximated by ~Þ ¼ u~h ðX

X

~Þ~ / I ðX uI þ

I2N

X

~ÞH ðX ~Þ~ / I ðX aI

ð2Þ

I2N C

where u~I denotes the classical degrees of freedom and a~I are the enriched degrees of freedom relative to node I. This approach is quite the same as the one proposed by Moes et al. [12] except that the Westergaard contribution has been removed. The numerical implementation of X-FEM presented in this paper was achieved for a 2-D analysis and a single type of element: quadrilateral four nodes element and an example illustrating the enrichment used is presented in Fig. 1b.

2.1. Crack modeling in 2-D 2.2. Discrete equations To model a given crack geometry within the X-FEM approach, a criterion for the selection of the enriched nodes is necessary. In the mainly used approach in the literature [12,17], the support of nodal shape functions is defined as the union of the elements connected to the node I. If this support of nodal shape functions is intersected by the crack, then the node I is enriched with a discontinuous ~Þ. The function based on the Heaviside step function H ðX ~Þ allows use of the generalized Heaviside step function H ðX to represent the ‘‘jump’’ in the displacement field across the ~Þ takes the value +1 crack. In the proposed approach, H ðX above the crack and 1 below the crack for a given direction of the crack:

The proposed approach is based on an elastodynamics behavior for the X-FEM analysis of a cracked homogeneous domain X, in the current configuration, as presented in Fig. 2. The crack is represented by the boundary Cc used t is applied to represent the two lips. A traction force vector ~ u is the applied displaceon the Neumann boundary Ct and ~ ment vector on the Dirichlet boundary Cu. It can be noted that Cu [ Ct = C and Cu \ Ct = ;. Crack lips are considered traction-free. Thus, we can write the strong form of the momentum conservation law in terms of the Cauchy stress tensor, for the current configuration described in Fig. 2, as follows:

576

I. Nistor et al. / Advances in Engineering Software 39 (2008) 573–587

Γu

n

M aa IJ ¼

Γc

Ω

Γt Fig. 2. Notations used for a 2-D domain.

orij þ qbi  q€ ui ¼ 0 2 X oxj rij nj ¼ ti 2 Ct ui ¼ ui 2 Cu

ð3Þ ð4Þ ð5Þ

where q is the current density, ~ b is the body force vector per unit mass, r is the Cauchy stress tensor, ~ n is the external unit vector to C and (ÆÆ) is the second time derivative of (). The weak form of the momentum equation in the current configuration is then given by Z Z Z dui q€ui dX ¼ dui qbi dX þ duiti dCt X X Ct Z oðdui Þ  rij dX ð6Þ oxj XnCc where ui is the trial displacement field (see Eq. (2) for the definition of ui) and dui is the test displacement field. The equilibrium discrete system of equations for dynamic analysis with X-FEM is obtained from Eq. (6) using the standard Bubnov–Galerkin procedure. Substituting trial and test displacement fields and their derivatives yield to the following system:  uu   " int #  ext  M IJ M ua € F iI uJ F iI IJ  þ ¼0 ð7Þ au aa int M IJ M IJ € Qext aJ QiI iI where: Z oð / I Þ ¼ rij dX F int iI XnCc oxj Z Z  F ext ¼ / dC þ /I qbi dX t I i iI Ct X Z oð/I H Þ Qint ¼ rij dX iI oxj XnCc Z Z  Qext ¼ ð/ H Þ t dC þ ð/I H Þqbi dX i I iI Ct X Z M uu ¼ q/I /J dX IJ ZX M ua q/I ð/J H ÞdX IJ ¼ X

qð/I H Þð/J H ÞdX

ð14Þ

X

u

t

Z

The details concerning the development of the terms in Eqs. 8–14 for the quadrilateral four nodes finite element are presented further in Section 4.3 and the assembling procedure in Section 4.4. In this work the consistent mass matrix is used because the enriched degrees of freedom obstruct its direct lumping as reported by Belytschko et al. [5]. In fact, neglecting these terms leads to suppress one of the essential information concerning the coupling of the regular and enriched degrees of freedom as de Borst et al. [18] observed. Even the lumping technique proposed by Menouillard et al. [15], for the case of Heaviside step function, does not takes into account the coupling terms. The use of the entire mass matrix increases the CPU time and requires a more powerful processor, but the analyzed models presented here usually contains a quite small number of degrees of freedom, therefore, this choice has been adopted in order to preserve the informations concerning the enrichment. 2.3. Explicit integration scheme The explicit integration procedure used in the X-FEM module DynaCrack uses the Chung–Hulbert [19] explicit time integration scheme already implemented in the DynELA code [20]. The time integration scheme is given by   int u€n M 1 F ext  aM ~ n  Fn ~ u€nþ1 ¼ ð15Þ 1  aM h i ~ u_nþ1 ¼ ~ u_n þ Dtnþ1 ð1  cÞ~ u€n þ c~ u€nþ1 ð16Þ    1 ~ u€n þ b~ b ~ u_n þ Dt2nþ1 un þ Dtnþ1~ u€nþ1 ð17Þ unþ1 ¼ ~ 2 The main feature of this algorithm is the presence of a numerical dissipation through its characteristic parameters aM, b and c. The values of these parameters are given by the following relations [20]: aM ¼

2qb  1 ; 1 þ qb

b¼

5  3qb 2

ð 1 þ qb Þ ð 2  q b Þ

;

c¼

3  aM 2 ð18Þ

ð8Þ ð9Þ ð10Þ ð11Þ ð12Þ ð13Þ

where qb 2 [0, 1] defines the numerical dissipative character of the algorithm. Setting qb = 1.0 leads to a conservative algorithm while qb < 1.0 introduces numerical dissipation in the scheme. In this work, the conservative algorithm is considered. The integration time step is computed using the following relation: Dt ¼ f

le wd

ð19Þ

where f is a safety factor (the value of f = 0.82 has been used here) that accounts for numerical instabilities, le is the characteristic length of the smallest element of the structure and wd is the dilatational elastic wave speed of the material. Note that the elements intersected by the

I. Nistor et al. / Advances in Engineering Software 39 (2008) 573–587

crack are not taken into account for this computation in order to avoid that the time step tend to zero when the crack path is very close to an enriched node.

577

sub-triangulation associated to a three or six-point integration rule. The entire geometrical procedure for partitioning the cut elements and the numerical integration algorithm are presented in Section 4.

2.4. Cut elements partition 3. Crack evolution model For the elements cut by the crack, a special procedure is applied in X-FEM in order to integrate the discrete system of equations. The main idea is not to perform a re-meshing of these elements by adding supplementary nodes on the intersecting points with the crack because this would be contrary to the main principle of X-FEM (modeling of the discontinuities without any dependence with the mesh size and orientation). The main accepted concept in XFEM is the partitioning of these elements. A clear analysis of ‘‘element partitioning versus re-meshing’’ is provided by Sukumar and Pre´vost [17]. The partitioning procedure in X-FEM is done for numerical integration purpose only, and no additional degrees of freedom are introduced into the discrete space during this operation. Subdividing these elements into triangles in 2-D was proposed by many authors [12,17,21], and several numerical integration options were also presented. On the other way, an original method to perform this integration without any subdivision of the cut elements has been recently proposed by Ventura [22]. Special (higher-order) quadrature rules are usually used for the numerical integration of the elements that are partitioned in this way [12] like a six-point integration rule for triangular elements. In this work, we subdivide the two zones on both sides of the crack into sub-quadrilaterals, as shown in Fig. 3. The main reason for this partitioning solution is related to the numerical integration accuracy. The numerical computation of the requested quantities such as stiffness matrix, on the sub-quadrilaterals is achieved using the same integration scheme as for all other element of the mesh. In the same time, bilinear shape functions are used to interpolate the fields in order to integrate them. This approach has given more accurate results than the classic

For the complete characterization of the dynamic crack propagation, beside the strong form given by Eqs. 3–5, a crack evolution model providing the crack advancing criteria (its direction and its velocity) is necessary. As mentioned earlier in the opening part in Section 2, the choice of the crack evolution model implemented in DynaCrack was strongly influenced by the enrichment. Since only the Heaviside function is taken into account for the crack modeling, the propagation is restricted from edge to edge of the elements. The limitations related to this approach are a quite inaccurate stress–strain evaluation in the immediate vicinity of the crack tips and a loss of smoothness for some fields because the crack propagates in fits and starts. On the other hand, in the context of adapting classical models of crack propagation to this enrichment type, this approach is very attractive for numerical integration and good results can be obtained as presented further in Section 5. Several propagation criteria were studied and implemented for dynamic crack propagation. The maximum circumferential stress criterion, also called the maximum principal stress criterion by Erdogan et al. [23] was implemented by Moes and Sukumar [12,17] in X-FEM for quasi-static propagation and by Belytschko et al. [5] for dynamic propagation. It sets that a crack will propagate from its actual tip in the direction hc where the circumferential stress rhh is maximum. Physical models used for computing crack propagation based on the energy release rate calculation represent an other important class of criteria for quasi-static and dynamic cracks. Freund [24] has developed a criterion giving an analytical solution for the DSIF (dynamic stress intensity factor) and the connection to the energy release rate for a dynamic stationary crack.

Fig. 3. Partitioning of a four nodes element.

578

I. Nistor et al. / Advances in Engineering Software 39 (2008) 573–587

Those criteria have been used in other numerical methods than the finite element method. Krysl and Belytschko [25] and Duarte et al. [26] implemented these criteria in EFGM (element free Galerkin method) and GFEM (generalized finite element method), respectively, for a 3D dynamic crack propagation. In our X-FEM code, we have chosen to introduce a physical crack evolution model based on Nishioka [27] and Freund [24] approaches and have adapted it to our considered enrichment. The definition of a physical crack evolution model enhances some particular problems since the modeling of the field in the immediate vicinity of the tip is not very accurate in our approach. Therefore it was difficult to obtain an accurate numerical solution around the crack-tip. One of the possible solutions is then to use the path-independent dynamic J-integral characterized by the following features [28]: • it has the physical meaning of a dynamic energy release rate, • it gives a unique value for an arbitrary path surrounding the crack-tip, • it can be related to the dynamic stress intensity factors (DSIF). The second feature mentioned above is the most interesting one since it allows to use a contour-path quite far from the crack-tip to evaluate the dynamic J-integral. This allow to avoid the inaccurately computed asymptotic field zone around the crack-tip. The analytical form of the J 0 integral developed by Nishioka and Alturi [29] for moving dynamic cracks is considered here: Z  J 0k ¼ ðW þ U Þnk  rij ui;k nj dC CþCc Z ð20Þ þ ðqu€i ui;k  qu_ i u_ i;k ÞdS S

where, as presented in Fig. 4, W and U are the strain and kinetic energy densities, respectively, S is an area inside  of C and Cc ¼ Cþ c þ Cc represents the crack edge inside of the considered contour C. Considering a plane strain approach, the J 0 -integral components denoted by k in Eq. (20) can be related to the DSIF by  1

AI ðcs ÞK 2I þ AII ðcs ÞK 2II ; 2l 1 ¼  AIV ðcs ÞK I K II 2l

J 00 1 ¼

J 00 2 ð21Þ

where AI,II,IV(cs) are coefficients depending on the propagation crack speed cs (see [24] for more details) and l is the shear modulus of the material. The numerical value of the dynamic energy release rate is given by G ¼ J 01 cos h0 þ J 02 sin h0

ð22Þ

DSIF numerical values are extracted using the separation components method proposed by Nishioka and Atluri

V–Vε

+

y0

Γ

θ0

Γc Γc

x0

– Vε

Γ ε

n

Fig. 4. Evaluation path scheme for the J 0 -integral.

[29]. This method has been retained because it avoids the need to compute the J 02 value since this one is known to be sensitive to the near crack-tip singular stress solution. In the separation components method, only the J 01 (non-affected by the near crack-tip singular solution) and the normal and tangential crack-tip opening displacement components (dn and dt) contributes to the DSIF evaluation: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2lJ 00 as  2 1 ; K I ¼ dn AI ðcs Þ dn as þ d2t ad sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2lJ 00 as  21  ð23Þ K II ¼ dt AII ðcs Þ dn as þ d2t ad where as and ad are functions depending on the crack propagation speed cs, the shear wave speed ws and dilatational wave speed wd by sffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffi c2s c2 as ¼ 1  2 ; ad ¼ 1  s2 ð24Þ ws wd The crack propagation model uses the DSIF to answer for the essential questions: is the crack propagates, and if it’s true, in what direction and how quickly? Hence, the crack will propagate if the value of the energy release rate G computed from Eq. (22) is greater or equal to a critical limit Gcrit given by Gcrit ¼ K ID ðcs Þ

1  m2 E

ð25Þ

where KID(cs) is the dynamic fracture toughness, assumed b ID a constant value considered as an intrinhere equal to K sic material property, E and m are the Young coefficient and the Poisson ratio, respectively. The direction of the crack propagation is then given by 8 0 s ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 19 2 <1 K = KI @ I  signðK II Þ hc ¼ 2 arctan þ 8A ; :4 K II ; K II if

K II 6¼ 0

ð26Þ

and by hc = 0 if KII = 0. The term sign(KII) inside of the previous equation leads to a positive stress intensity factor along the direction given by hc. The crack speed is provided by the numerical propagation algorithm, since the crack-tip advances one edge at a time. As mentioned in Ref. [5], this edge-to-edge approach achieves closure by the discretiza-

I. Nistor et al. / Advances in Engineering Software 39 (2008) 573–587

tion but the coarse meshes must be avoid. The numerical algorithms implemented for computing the parameters of this crack evolution model will be presented in Section 4.

579

N

4. Numerical implementation of the DynaCrack module

R The DynELA explicit finite element program [16] has been used for the implementation of the X-FEM crack propagation module. This FEM program is written in C++ and developed within an object-oriented programming (OOP) approach. This feature presents a very well defined mechanism for modular design and re-use of code: in our case this has allowed to develop the new module DynaCrack, re-using several objects already implemented. Some new classes have been implemented and other have been specialized for this application using the inheritance mechanism. In the following sub-sections, we present the major steps of this implementation, pointing out the characteristic keys of X-FEM. Concerning the notation, all keywords related to the DynaCrack code such as the class names, method names, object names . . . are written using a typewriter font. 4.1. Crack representation and implementation The crack is represented by the class CrackFunction containing a list of p + 1 points Xc = {xi}i2[0,p]. Within this representation, the two points x0 and xp are the crack tips. The initial crack at time t = 0 is defined through the input data file according to the initial geometry of the crack. The CrackFunction class contains a dynamic list updated when adding new points as the crack propagates. One of the most useful method of this class returns the sign of the Heaviside function for any given point. The so called method IsOnPositiveSide() returns a Boolean value: true for positive side (i.e. the Heaviside function H = 1) and false for negative side (i.e. the Heaviside function H = 1) based on geometric predicate evaluation. This is done by computing the projection of this point onto the closest crack segment. In the presented X-FEM approach, the mesh contains both classical and enriched nodes, therefore, the class XNode, inherited from the class Node is used for handling the supplementary degrees of freedom (dof) corresponding to the enriched nodes. A specific flag XType is then used to identify the status of the nodes. During the computation, the status of some nodes changes from normal to enriched because of the crack propagation. The method getXRType() is called at each integration time step and is used to update the status of the elements. A method based on an efficient BoundingBox algorithm is used to find the nodes around the crack. For all the elements of the structure three flags are considered according to theirs status: 8 > < TypeX  enriched elements ðfrom 1 to 3 enriched nodesÞ x elements TypeR  completely cut elements ð4 enriched nodesÞ > : TypeN  normal elements ðno enriched nodesÞ

X

y

x Fig. 5. Enriched element types definition.

Fig. 5 illustrates this aspect for an X-FEM mesh; the enriched nodes are the encircled ones. The two nodes at both end of the edge where the crack-tip is situated are enriched with the Heaviside step function. In the approach proposed by Moes et al. [12] those two nodes are enriched using the Westergaard functions leading to the step function for h = p. 4.2. Partition algorithm Once the elements status have been updated for the current time step, the partitioning of the so-called TypeR elements is done using the methods of the XPartition class. An arbitrary crack geometry in a structure discretized with 4-nodes quadrilateral finite elements leads to one of the two situations for cut elements: the crack intersects two opposite or two adjoining edges (see Fig. 3 for more details). The algorithms implemented in the XPartition class compute the partition of the elements. As an illustration, the algorithm for the partition of the element presented in Fig. 6 is reported in Box 1.

Box 1 Partitioning algorithm 1. Compute the intersection points, xc1 and xc2; 2. Build the {pi}i=1. . .4 xpoints using nodes location; 3. Compute the positive and negative centroids, C+ and C; 4. Get the median points, mj, for all sub-domain segments; 5. Build pointsPositifs and pointsNegatifs lists; 6. Build Surf2D surfaces.

580

I. Nistor et al. / Advances in Engineering Software 39 (2008) 573–587

(p + 2) n3

(p + 3) m5

(p + 4) n4

+

+ S2

(p+5 ) m6

+

–

S4

m1 (p 1– )

7

(p – ) m7 7

(+1,+1)

xc2 (p + ) 8 – (p )

g4

g3

ξ

6

–

S3

–

S2 m2 (p 3– )

x

g1

m3 (p 5– )

– C

S1 n1 (p 2– )

S4

(p + )

–

y

η (–1,+1)

+

+

(p 8– )

(p + 1) m4

C

S3

(p + 6 ) xc1

S1

g2 (+1,–1)

(–1,–1)

n2 (p 4– )

Fig. 6. Integration of a partitioned element.

Two more classes are used to hold the informations for points and surfaces created during the partition, XPoint and Surf2D, respectively. As the element is sub-divided into two sub-domains, a ‘‘positive’’ one and a ‘‘negative’’ one (according to enriched nodes sign), two lists of XPoint objects are created, pointsPositifs fpþ i g and pointsNegatifs fp g, as reported in Fig. 6. The j XPoint list contains the intersection points between the crack and the element edges, xc1 and xc2, the node locations and all middle points ml of the edges arranged in counterclockwise order. From the computation of the sub-domain centroids (C+ and C), two lists of Surf2D quadrilateral  sub-surfaces are generated, fS þ i g and fS i g. Those two sub-surfaces are used for the numerical integration of the conservative laws and their respective contributions are summed. 4.3. Setting-up of the matrices One of the main consequence of the additionals dof from a numerical point of view is the variable size of the elementary matrices, according to the element status: cut, enriched or normal. In this work, we adopted a block formalism for the setting-up of the matrices as described here after. The numerical implementation of the DynaCrack module has been done for a 4-node quadrilateral element. For a standard element (TypeN), the elementary dof vector contains eight terms (two dof for each node) and the displacement field is approximated by standard shape functions. For the enriched or cut elements, the size of the elementary dof vector is larger (10, 12 or 14 dof for TypeX and 16 dof for TypeR). A choice concerning the dof placement is necessary: do we put all enriched dof in the second part of the vector after the classical ones as proposed by Sukumar et al. [17] or do we keep in order all dof (both classical

and enriched) for each node? As this choice has no real impact on the final solution, we adopted the later and called this one the ‘‘block approach’’ as each node is stored with his ‘‘block’’ of dof. For a TypeR element, the elementary dof vector is therefore given by T

ue ¼ f ux1 ;

uy1 ;

ax1 ; ay1 ; . . . ;

ux4 ; uy4 ; ax4 ; ay4 g ð27Þ

where uxi, uyi are the standard dof and axi, ayi are the enriched dof for the node i. Starting from the displacement approximation Eq. (2), the shape functions matrix N and the derivatives of the shape functions matrix B, for the TypeR element, are given by N ¼ ½ Ns1

Nh1

Ns2

B ¼ ½ Bs1

Bh1

Bs2

Nh2 Bh2

Ns3 Bs3

Nh3 Bh3

Ns4 Bs4

Nh4  Bh4 

where the block matrices Ns, Nh, Bs and Bh are 2 3 /i;x 0 " # /i 0 6 0 / 7 Bsi ¼ 4 Nsi ¼ i;y 5; 0 /i /i;y /i;x 2 3 H /i;x 0 " # H /i 0 6 0 7 H /i;y 5; Nhi ¼ Bhi ¼ 4 0 H /i H /i;y H /i;x

ð28Þ ð29Þ

ð30Þ

ð31Þ

For TypeX elements, both N and B matrices contains Nh and Bh blocks only for the enriched nodes. A detailed flowchart of the algorithm for the computing of the elementary stiffness matrix for a TypeR element is reported in Fig. 7. This one returns the stiffness matrix resulting from the numerical integration over both positive and negative domains Surf2D by a Gaussian quadrature (see Fig. 6). The algorithm for the mass matrix computation is obtained straightforward.

I. Nistor et al. / Advances in Engineering Software 39 (2008) 573–587

581

manner, looping over all elements of the mesh and taking into account the status of each node. 4.5. Implementation of the crack evolution model Numerical algorithms implemented for the crack evolution model refer mainly to the energy release rate computation and the DSIF extraction presented earlier. The components of the J 0 -integral are numerically evaluated in DynaCrack considering a symmetrical rectangular path centered at the crack-tip as presented in Fig. 9. As the edges of the path are set parallel to the system axis and the terms inside Eq. (20) are computed with respect to the normal unit vector orientation, the first term in Eq. (20), noted J Ck , is given by DE EF FA B J Ck ¼ J AD k þ Jk þ Jk þ Jk þ Jk

C

þ

þ J CB k

ð32Þ

where the six components in the right hand side are computed along the corresponding segments with respect to notations reported in Fig. 9. The first component of the J 0 -integral (for k = 1) is given by  Z D ou ov AD J1 ¼ ðW þ U Þ þ r11 þ r12 ð33Þ dx2 ox1 ox1 A   Z E ov ou J DE ¼ r22 þ r12 ð34Þ dx1 1 ox1 ox1 D  Z F ou ov J EF ¼ ð W þ U Þ  r  r ð35Þ dx2 11 12 1 ox1 ox1 E  Z A ov ou J FA ¼ r  r ð36Þ dx1 22 12 1 ox1 ox1 F  Z C ov ou  J 1B C ¼ r22 þ r12 ð37Þ dx1 ox1 ox1 B  Z Bþ  ov ou CBþ r22  r12 ð38Þ J1 ¼ dx1 ox1 ox1 C

Fig. 7. Flowchart for stiffness matrix computation of a partitioned element.

4.4. Assembly procedure The global mass matrix, stiffness matrix and external force vector assembling procedure is specific since the elementary corresponding matrices sizes differ from element to element depending on the status. To get round this, the XAllocation class has been specifically developed. The main feature of this class is to allow a dynamic mapping between both local and global positions of a dof as illustrated in Fig. 8. This bi-directional link is performed by the methods loc2glob and glob2loc. The former returns the global position of a dof depending on the node Id and the local dof position while the later gives the reversed mapping. The assembly of the mass matrix, stiffness matrix and external force vector are achieved in a like

The second component for k = 2 are obtained in a straightforward manner. In the above equations, the strain and the kinetic energies are given by 1 W ¼ ½r11 e11 þ r22 e22 þ 2r12 e12  2 1 U ¼ q½u_ 1 u_ 1 þ u_ 2 u_ 2  2 The second term, noted J S1 , of Eq. (20) is given by    Z  ou ov ou_ o_v q €u þ €v þ v_ J S1 ¼  u_ dS ox1 ox1 ox1 ox1 S

ð39Þ ð40Þ

ð41Þ

The computation of this integral is numerically done by evaluating the integrand terms for all Gauss points inside of the considered path. Concerning the elements intersected by the integral path, in light gray in Fig. 9, only theirs integration points inside of the path are considered for the above computation. This is done by the method computeJ of the XExplicitSolver class dedicated to the explicit time integration and the crack propagation processing. A partitioning algorithm is used to split the path’s

582

I. Nistor et al. / Advances in Engineering Software 39 (2008) 573–587

Fig. 8. Allocation procedure illustration.

d¼2

4 X

ð42Þ

N j aj

j¼1

Some numerical results proving the effectiveness of the numerical evaluation of the J 0 -integral are presented in Section 5. 5. Numerical examples In this section, we present some numerical results concerning the propagation of a crack in Mode I for a rectangular plate subjected to an impact load and the analysis of a finite plate with an inclined crack subjected to a mixed mode fracture. The first problem, illustrated in Fig. 10, was proposed by Lu et al. [31], Krysl et al. [25] and Belytschko et al.

y

Fig. 9. J-integral rectangular path.

segments into a set of equal sub-segments Dx, and to generate a list of geometric points corresponding to the middle of those sub-segments. The numerical values for stresses, ou ov strains, ox and ox are interpolated to the above mentioned i i points. The numerical integration is based on a classic Gauss quadrature. As observed in Fig. 9, this evaluation is done over a 8L · 8L quadrilateral domain (L being the largest edge of the element containing the crack-tip). Once the J 0 -integral is computed, the dynamic energy release rate is determined by Eq. (22) and the DSIF components are extracted using Eq. (23). Chessa et al. [30] have shown that the numerical evaluation of the crack-tip opening displacements (dn and dt), contributing to KI and KII, depends only on the enrichment dof by

σ (t)

h

a h/2

x

σ (t) l Fig. 10. Model used for the crack propagation.

I. Nistor et al. / Advances in Engineering Software 39 (2008) 573–587

[32] using EFGM, by Duarte et al. [26] using GFEM, and by Belytschko et al. [5] using X-FEM among many others. Two uniform quadrilateral meshes are presented in this paper to illustrate the independence of the crack propaga-

250

7e+06

Dynamic stress intensity factor KI

Energy release rate G

150

100

50

0 0

tion with the mesh size: a coarse mesh (51 · 21 elements), and a finer mesh (65 · 25 elements). The crack-tip is located at the center of the plate, as shown in Fig. 10, and the crack is horizontal.

fine mesh coarse mesh theory

200

0.0005

0.001

0.0015

0.002

583

fine mesh coarse mesh theory

6e+06 5e+06 4e+06 3e+06 2e+06 1e+06 0 0

0.0005

Time (s)

0.001

0.0015

0.002

Time (s)

Fig. 11. Energy release rate time-history plot for the stationary crack.

Fig. 12. Mode I DSIF time-history plot for the stationary crack.

Fig. 13. von Mises stress field for the stationary crack.

584

I. Nistor et al. / Advances in Engineering Software 39 (2008) 573–587

A linear elastic material behavior with E = 200 GPa, m = 0.3 and q = 7833 kg/m3 is considered and a critical pffiffiffiffi value for fracture toughness K IDcrit ¼ 1:5 MPa m is used in the dynamic crack evolution model. Only a stationary crack is considered in the first part and the dimensions are l = 10 m and h = 4m. Only a boundary traction force with a magnitude of 63750 Pa is applied in accordance with the model proposed by Duarte et al. [26]. The computed and the analytical mode-I DSIF and the dynamic energy release rate reported in Figs. 11 and 12. Both figures show a quite good agreement and the non-dependence of the solution with the mesh size. On both figures, one can observe, after t = 7.5 · 104 s, quite strong oscillations disturbing the numerical solutions. Actually, as mentioned in [26], for the same test, it occurs because of the finite dimensions of the plate whose boundary reflect the elastic waves. This elastic waves perturbing the displacements, strains and stresses fields used for the J 0 -integral evaluation. The explicit integration algorithm favour and maintain also these oscillations in absence of plasticity behavior. Fig. 13 is a contour-plot of the von Mises stresses field at the end of the computation for the coarse mesh. In this later the expected form of the stress field gradient near

the crack-tip can be observed. Fig. 14 reports the vertical displacement field contour-plot and shows the jump in the vertical displacement on both sides of the crack. It must be mentioned that because of the finite dimensions of the plate, some elastic waves are reflected by the boundaries when the simulation time increases. Those reflected elastic waves are perturbing the numerical solution and this one is no longer in accordance with the analytic solution. For the analysis of the crack propagation, the dimensions of the rectangular plate are l = 0.1 m and h = 0.04 m and the applied traction force magnitude is 1 MPa. The same conditions were adopted by Belytschko et al. [5]. Fig. 15 reports the evolution of the mode I DSIF for both fine and coarse meshes. The crack propagation initiates quite at the same instant for the two considered meshes, t ’ 14 ls for the fine mesh and t ’ 15 ls for the coarse mesh. The crack-tip advances through 8 elements for the fine mesh and through seven elements for the coarse mesh. The exact final crack-tip propagations are dx = 6.23 · 102 m and dx = 6.37 · 102 m for the fine and coarse meshes, respectively. The crack propagation speed is in the range [1000, 2200] m/s. This is below the Rayleigh wave speed. Those results agree quite well with

Fig. 14. Vertical displacement component for the stationary crack.

I. Nistor et al. / Advances in Engineering Software 39 (2008) 573–587

Fig. 15. Mode I DSIF time-history plot for the moving crack.

the one obtained by Belytschko et al. [5] and allow to validate the proposed approach. The second problem presented in Fig. 16a is the analysis of an inclined crack in a finite plate. The purpose here is to compare the numerical results obtained by DynaCrack for the DSIF in a mixed mode fracture with the results obtained by the Abaqus FEM code for the same simulation. The reference length of the plate is L = 1 m and a boundary traction force with a magnitude r = 1 · 103 Pa is applied via a step function. The material considered is the same as the one used in the previous example and three different uniform quadrilateral meshes are considered: (i) 20 · 35, (ii) 25 · 44 and (iii) 30 · 53 elements. The equiva-

585

lent model simulated with the Abaqus software is composed of 911 non-structured elements and the mesh conforms with the crack geometry. Stress distribution fields given by the X-FEM analysis at the end of the computation time (t = 0.1 s) are illustrated in Fig. 16b. Here again, the stress field has the expected shape around the crack-tip. Table 1 reports the numerical values for the fracture parameters (DSIF and dynamic energy release rate) obtained at the end of the analysis by both FEM codes. The results obtained by the path-independent integral technique with DynaCrack agree quite well with the ones obtained using Abaqus. The differences observed for the KII values are related to the inaccurate field representation around the crack-tip because we used a quite simplistic X-FEM enrichment used in this work. More explicitly, these differences are related to the different manner used in Abaqus for computing KII since the second equation for extracting KI and KII (the first one being the one of the energy release rate) is related to the stress field around crack-tip, obviously more accurately represented in this case.

Table 1 Numerical results obtained for the inclined crack problem Parameter

Fine mesh

Middle mesh

Coarse mesh

Abaqus

Elements number G [J/m2] pffiffiffiffi K I ½Pa m pffiffiffiffi K II ½Pa m

1590 1.18 · 103 1521.7 404.4

1100 1.26 · 103 1617.9 460.9

700 1.36 · 103 1624.2 493.8

911 1.23 · 103 1480.4 722.4

Fig. 16. Inclined crack problem in a finite plate.

586

I. Nistor et al. / Advances in Engineering Software 39 (2008) 573–587

6. Conclusions The development of a dynamic crack propagation module based on the X-FEM approach and its implementation in an explicit finite element code has been presented in this paper. Some specific algorithms have been developed and programmed using the C++ language in order to implement the main features of the X-FEM. The main challenges of this work concerns the treatment of the dynamic nodes enrichments with additional degrees of freedom for the representation of the crack while this one propagates through the structure. The approach used for the numerical integration of the mass and stiffness matrices and the vector force have also been presented with the retained approach to manage of a variable size for the global degree of freedom vector. The implementation of a crack propagation criterion adapted to our specific crack-tip enrichment has also been presented. Finally two numerical examples have also been presented to show the effectiveness of the proposed algorithms. Concerning the results obtained for these numerical examples, it was clearly pointed out that some supplementary works on our code are needed in order to improve them. First of all a new type of finite element (3-node triangle) will be tested for searching a better conditioning of the elementary matrix. The impact of the Westergaard cracktip enrichment must be studied also to evaluate the potential gain in accuracy. Several contributions on the dynamic analysis using X-FEM were published since our works began, especially concerning the numerical integration schemes, and next we will try to apply some of such concepts recently developed in order to improve the time integration algorithm. Future works and further developments of this code concerns the initiation of the crack, i.e. the evolution from a continuous structure to a fissured one, the computation of the crack propagation within a plastic deformable body and the introduction of the contact between the lips of the crack. The most recent developments concern a new formalism for the treatment of the crack based on a mixed extended element free Galerkin (X-EFG) and a FEM model. The main advantage of this approach is to propose an implicit enrichment, i.e. no additional dof are used. This may allow to get round of some problems encountered in this work, for example the mass matrix diagonalization. References [1] Xu XP, Needleman A. Numerical simulation of fast crack growth in brittle solids. J Mech Phys Solids 1994;42:1397–434. [2] Camacho GT, Ortiz M. Computational modeling of impact damage in brittle materials. Int J Solids Struct 1996;33:2899–938. [3] Remmers, JC, Borst, R, Needleman, A. Simulation of fast crack growth using cohesive segments. In: VII international conference on computational plasticity, COMPLAS; 2003. [4] Babuska I, Melenk JM. The partition of unity method. Int J Numer Meth Eng 1997;40:727–58.

[5] Belytschko T, Chen H, Xu J, Zi G. Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment. Int J Numer Meth Eng 2003;58:1873–905. [6] Belytschko T, Fish J, Engelmann BE. A finite element with embedded localisation zones. Comput Meth Appl Mech Eng 1988;70:59–89. [7] Dvorkin EN. Finite elements with displacements interpolated embedded localization lines insensitive to mesh size and distorsions. Int J Numer Meth Eng 1990;30:541–64. [8] Simo JC, Oliver J, Armero F. An analysis of strong discontinuities induced by strain softening in rate-independent inelastic solids. Comput Mech 1993;12:277–96. [9] Jirasek M. Comparative study on finite elements with embedded discontinuities. Comput Meth Appl Mech Eng 2000;188:307–30. [10] Oliver J, Huespe AE. Theoretical and computational issues in modelling material failure in strong discontinuity scenarios. Comput Meth Appl Mech Eng 2004;193:3195–220. [11] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. Int J Numer Meth Eng 1999;45:601–20. [12] Moes N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. Int J Numer Meth Eng 1999;46: 131–50. [13] Dolbow J, Moes N, Belytschko T. Discontinuous enrichment in finite elements with a partition of unity method. Finite Elem Anal Des 2000;36:235–60. [14] Re´thore´ J, Gravouil A, Combescure A. A stable numerical scheme for the finite element simulation of dynamic crack propagation with remeshing. Comput Meth Appl Mech Eng 2004;193:4493–510. [15] Menouillard T, Re´thore´ J, Combescure A, Bung H. Efficient explicit time stepping of the extended finite element method (X-FEM). Int J Numer Meth Eng 2006;68(9):911–39. [16] Pantale´ O. An object-oriented programming of an explicit dynamics code: application to impact simulation. Adv Eng Softw 2002;33(5):275–84. [17] Sukumar N, Pre´vost JH. Modelling quasi-static crack growth with the extended finite element method. Part I: Computer implementation. Int J Solids Struct 2003;40:7513–37. [18] Borst R, Remmers JJC, Needleman A. Mesh-independent discrete numerical representations of cohesive-zone models. Eng Fract Mech 2006;73:160–77. [19] Hulbert GM, Chung J. Explicit time integration for structural dynamics with optimal numerical dissipation. Comput Meth Appl Mech Eng 1996;137:175–88. [20] Pantale´ O. Parallelization of an object-oriented FEM dynamics code: influence of the strategies on the speedup. Adv Eng Softw 2005;36(6):361–73. [21] Dolbow JE, Devan A. Enrichment of enhanced assumed strain approximations for representing strong discontinuities: addressing volumetric incompressibility and the discontinuous patch test. Int J Numer Meth Eng 2004;59:47–67. [22] Ventura G. On the elimination of quadrature subcells for discontinuous functions in the extended finite-element method. Int J Numer Meth Eng 2006;66:761–95. [23] Erdogan F, Sih G. On the crack extention in plates under plane loading and traverse shear. J Basic Eng 1963;85:519–27. [24] Freund LB. Dynamic fracture mechanics. Cambridge University Press; 1998. [25] Krysl P, Belytschko T. The element free Galerkin method for dynamic propagation of arbitrary 3-D cracks. Int J Nume Meth Eng 1999;44:767–800. [26] Duarte CA, Hamzeh ON, Liszka TJ, Tworzydlo WW. A generalized finite element method for the simulation of three-dimensional dynamic crack propagation. Comput Meth Appl Mech Eng 2001;190:2227–62. [27] Nishioka T. Computational aspects of dynamic fracture in comprehensive structural integrity. Elsevier; 2003. [28] Nishioka T. Computational dynamic fracture mechanics. Int J Frac 1997;86:127–59.

I. Nistor et al. / Advances in Engineering Software 39 (2008) 573–587 [29] Nishioka T, Atluri SN. On the computation of mixed-mode k-factors for a dynamically propagating crack, using path-independent integrals j’. Eng Frac Mech 1984;20:193–208. [30] Chessa C, Belytschko T. A local space-time discontinuous finite element method. Comput Meth Appl Mech Eng 2006;195:1325–43.

587

[31] Lu YY, Belytschko T, Tabbara M. Element-free Galerkin method for wave propagation and dynamic fracture. Comput Meth Appl Mech Eng 1995;126:131–53. [32] Belytschko T, Tabbara M. Dynamic fracture using element-free Galerkin methods. Int J Numer Meth Eng 1996;39:923–38.