A Phantom-Node Method with Edge-Based Strain Smoothing for Linear Elastic Fracture Mechanics

Abstract

This paper presents a novel numerical procedure based on the combination of an edge-based smoothed finite element (ES-FEM) with a phantom-node method for 2D linear elastic fracture mechanics. In the standard phantom-node method, the cracks are formulated by adding phantom nodes, and the cracked element is replaced by two new superimposed elements. This approach is quite simple to implement into existing explicit finite element programs. The shape functions associated with discontinuous elements are similar to those of the standard finite elements, which leads to certain simplification with implementing in the existing codes. The phantom-node method allows modeling discontinuities at an arbitrary location in the mesh. The ES-FEM model owns a close-to-exact stiffness that is much softer than lower-order finite element methods (FEM). Taking advantage of both the ES-FEM and the phantom-node method, we introduce an edge-based strain smoothing technique for the phantom-node method. Numerical results show that the proposed method achieves high accuracy compared with the extended finite element method (XFEM) and other reference solutions.

(ProQuest: ... denotes non-US-ASCII text omitted.)

N. Vu-Bac 1,2 and H. Nguyen-Xuan 3 and L. Chen 4 and C. K. Lee 5 and G. Zi 6 and X. Zhuang 1 and G. R. Liu 7 and T. Rabczuk 2,8

Recommended by Song Cen

1, Department of Geotechnical Engineering, Tongji University, Shanghai 200092, China
2, Institute of Structural Mechanics, Bauhaus-University Weimar, Marienstraße 15, D-99423 Weimar, Germany
3, Division of Computational Mechanics, Ton Duc Thang University, Ho Chi Minh City 70000, Vietnam
4, School of Chemistry, Physics and Mechanical Engineering, Queensland University of Technology, Brisbane, QLD 4001, Australia
5, School of Engineering, Institute of Mechanics and Advanced Materials, Theoretical and Computational Mechanics, Cardiff University, Wales CF24 3AA, UK
6, Department of Civil, Environmental & Architectural Engineering, Korea University, 5 Ga 1, An-Am Dong, Sung-Buk Gu, Seoul 136-701, Republic of Korea
7, Aerospace Systems Ohio Eminent Scholar, University of Cincinnati, Cincinnati, OH 45221-0070, USA
8, School of Civil, Environmental and Architectural Engineering, Korea University, 5 Ga 1, Anam-dong, Seongbuk-gu, Seoul 136-701, Republic of Korea

Received 15 January 2013; Accepted 1 May 2013

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

The extended finite element method (XFEM) [1] has become a standard tool to model arbitrary crack growth. However, the implementation of XFEM in an existing finite element code requires severe modifications. An alternative method to model arbitrary crack growth was proposed by [2], subsequently implemented by [3] in a static setting and by [4] in a dynamic setting. It was shown by [4, 5] that the method proposed by [2] is identical to a step-enriched XFEM; [4] refer to this method as phantom node method. The main difference to the original XFEM is that the discontinuity jump is not obtained by introducing additional unknowns but by the so-called overlapping paired elements. In other words, a new "overlapped" element is introduced to handle the crack kinematics when an underlying element is cracked. It is accomplished by integrating these overlapped elements up to the crack. Though it was shown by [4] that the crack kinematics obtained with the phantom node method are identical to the step-enriched XFEM, and it has some advantages over step-enriched XFEM.

(1) As no additional degrees of freedom are introduced, the implementation of the phantom node method in an existing finite element code is simpler. For example, arbitrary crack growths for nonlinear materials and cohesive zone models even for multiple cracks in two and three dimensions have already been implemented in ABAQUS [6] while an additional plug-in [7] is required to model crack growth using XFEM.

(2) No mixed terms ( ^{K u a} and ^{K a u} ) occur improving conditioning.

(3) Standard mass lumping schemes can be used due to the absence of an enrichment. There are several contributions to develop diagonalized mass matrices in standard XFEM [8, 9], but they are based on certain assumptions.

(4) The development of complex FE-formulations is much easier due to the lack of an enrichment. For example: when techniques such as EAS (enhanced assumed strain) or ANS (assumed natural strain) are used, special attention is required in a standard XFEM-formulation, particularly for problems with constraints. Those difficulties do not occur in the phantom node method [10].

The key drawback of the phantom node method compared to standard XFEM is its lower flexibility. It was developed for problems involving crack growth "only." However, avoiding a crack tip enrichment significantly facilitates the enrichment strategy and the crack tracking algorithm.

(i) A crack tip enrichment introduces more additional unknowns. It is well known that a topological enrichment is needed for accuracy reasons [11] leading to a substantial increase of additional unknowns (compared to "pure" step-enriched formulations) and increasing difficulties due to increasing the condition number.

(ii) The nonpolynomial (and singular) crack-tip enrichment complicates integration [12-15] and requires special attention (blending).

(iii): The enrichment strategy and the crack growth algorithms are complicated, in particular in 3D.

Modeling crack growth with the phantom node method on the other hand is quite simple. Commonly, plane crack segments are introduced through the entire element though crack tip elements were developed [16] that allow cracks to close inside an element.

Recently, Liu et al. constructed a new class of finite element methods based on strain smoothing. Among those methods, the so-called ES-FEM edge-based smoothed finite element method (ES-FEM) has been proven to be the most efficient and accurate one. In numerous application [17-19], it was shown that particularly low-order SFEM-formulations are superior in terms of efficiency and accuracy over "standard" low-order finite element formulations. In particular, it was shown for many applications [15, 20] that results obtained by triangular ES-FEM elements are of the same accuracy as standard Q4-elements.

Therefore, we propose to couple the ES-FEM with the phantom-node method. We name the new element edge-based phantom node method (ES-Phantom node). In this paper, we focus on two-dimensional problems in linear elastic fracture mechanics (LEFM). However, our long term goal is to model fracture in nonlinear materials in 3D. Numerical results show high reliability of the present method for analysis of fracture problems.

This paper is organized as follows. In Section 2, we briefly summarize the basic theory of phantom-node method. A brief description of ES-FEM is called back in Section 3. The combination between the phantom-node method and the ES-FEM is elaborated in Section 4. Section 5 presents the integration technique. Benchmark numerical problems taken from linear elastic fracture mechanics are studied in Section 6. Finally, we give some concluding remarks.

2. A Brief Description of Phantom-Node Method

Consider a deformable body occupying domain Ω in motion, subjected to body forces b , external applied traction t on boundary _{Γ t} , and displacement boundary conditions u = u - on _{Γ u} containing a crack as shown in Figure 1 with the corresponding finite element discretization. In the phantom-node method, a completely cracked element is replaced by two partially active superimposed elements 1 and 2 whose nodes consist of real nodes and phantom nodes marked by solid and empty circles, respectively. The active part of element 1 ( _{Ω 1} ) , ^{u 1} ( x ) , which holds for f ( x ) < 0 and the other active part ( _{Ω 2} ) , ^{u 2} ( x ) , which holds for f ( x ) > 0 . The two parts of the model do not share nodes, and therefore they displace (deform) independently. Areias and Belytschko [5] demonstrated that the Hansbo and Hansbo [2] formulation is equivalent to the XFEM formulation relying on discontinuous enrichment with the Heaviside function.

Figure 1: A two-dimensional body containing a crack and boundary conditions.

[figure omitted; refer to PDF]

The displacement field within an element _{Ω e} in Figure 2 is rewritten as [4] [figure omitted; refer to PDF] where _{S 1} and _{S 2} are the nodes of superimposed elements 1 and 2, respectively. As illustrated in Figure 2, each element contains real nodes and phantom nodes marked by solid and empty circles, respectively; _{N I} is the finite element shape function associated with node I , while ^{u I 1} and ^{u I 2} are nodal displacements of original nodes in superimposed element 1 and 2, respectively.

Figure 2: The decomposition of a cracked element into two superimposed elements.

[figure omitted; refer to PDF]

H is the Heaviside function given in [1, 21-24] and defined by [figure omitted; refer to PDF]

Here, we choose the physical domain up to the crack line. Note that the crack line is a boundary in phantom node method. It is like the elements near the external boundary. So, we avoid singularity in phantom node method. The corresponding strain terms are written the same.

The strain field is obtained as follows: [figure omitted; refer to PDF] where _{B I} is the standard strain-displacement matrix.

The jump in the displacement field across the crack is calculated by [figure omitted; refer to PDF]

In this paper, the crack tip is forced to be located on the element's boundary.

3. Brief on Edge-Based Strain Smoothing Method in Finite Elements

3.1. Displacement and Strain Field

In the ES-FEM [25], the domain Ω is partitioned into a set of nonoverlapping no-gap smoothing domains constructed using element edges of the triangular elements. ^{Ω ( k )} satisfies the conditions Ω = ^{... k = 1 _{N e}} ... ^{Ω ( k )} and ^{Ω ( i )} ∩ ^{Ω ( j )} = ∅ , for all i ...0; j , in which _{N e} is the total number of edges of elements in the problem domain. In Figure 3, the smoothing domain, the smoothing domain corresponding to an inner edge k , and the smoothing domain for a boundary edge m are illustrated.

Figure 3: Construction of edge-based strain smoothing domains.

[figure omitted; refer to PDF]

Numerical integration is implemented on chosen Gauss points as illustrated in Figures 6 and 7 corresponding with split smoothing domain in Figure 4 and tip smoothing domain in Figure 5, respectively.

Figure 4: The decomposition of a completely cracked smoothing domain into two superimposed smoothing domains.

[figure omitted; refer to PDF]

Figure 5: The decomposition of a cracked smoothing domain containing crack tip into two superimposed smoothing domains.

[figure omitted; refer to PDF]

Figure 6: The decomposition of a completely cracked smoothing domain into two superimposed smoothing domains.

[figure omitted; refer to PDF]

Figure 7: The decomposition of a cracked smoothing domain containing crack tip into two superimposed smoothing domains.

[figure omitted; refer to PDF]

Distribution of the stress components _{σ x x} and _{σ y y} for and unstructured mesh are shown in Figures 27 and 28, respectively.

Introducing the edge-based smoothing operation, the compatible strain ... = _{∇ s} ^{u k h} is smoothed over cell ^{Ω ( k )} associated with edge k as follows: [figure omitted; refer to PDF] where _{Φ k} is a given normalized smoothing function that satisfies [figure omitted; refer to PDF]

Using the following constant smoothing function [figure omitted; refer to PDF] where ^{A ( k )} = _{∫ ^{Ω ( k )}} d Ω is the area of the smoothing domain ^{Ω ( k )} , ^{Γ ( k )} is the boundary of the smoothing domain ^{Ω ( k )} , and _{L n} is the outward unit normal matrix which can be expressed as [figure omitted; refer to PDF]

4. Edge-Based Strain Smoothing Phantom Node Method

4.1. Displacement and Strain Field

The approximation of the displacement field is written similarly to (1): [figure omitted; refer to PDF] where ^{S 1 es-pht} and ^{S 2 es-pht} are nodes associated with smoothing domains 1 and 2, respectively, consisting of real nodes and phantom nodes illustrated in Figures 4 and 5. The associated nodes of the inner smoothing domain ^{Ω ( k )} (DEFG) and boundary smoothing domain ^{Ω ( m )} (ABC) are shown in Figure 3.

The connectivities of these superimposed smoothing domains which are cracked completely and the corresponding active parts are shown in Figure 4:

: nodes of smoothing domain 1 ( ^{Ω 1 ( k )} ) = [ 1,2 , 3 ~ , 4 ~ ] ;

: nodes of smoothing domain 2 ( ^{Ω 2 ( k )} ) = [ 1 ~ , 2 ~ , 3,4 ] .

The connectivity of a superimposed smoothing domain containing the crack tip and the corresponding active parts is shown in Figure 5 so that crack tip is guaranteed to locate on the element's edge:

: nodes of smoothing domain 1 ( ^{Ω 1 ( k )} ) = [ 1,2 , 4 ] ;

: nodes of smoothing domain 2 ( ^{Ω 2 ( k )} ) = [ 1 ~ , 2,4 ] .

Using the strain smoothing operation, the smoothed strain associated with edge k created from the displacement approximation in (9) can be rewritten as [figure omitted; refer to PDF] where _{B - I} ( _{x k} ) is the smoothed strain gradient matrix for the standard ES-FEM part. Those matrices write as follows [figure omitted; refer to PDF]

In (11), _{b - I h} ( _{x k} ) , h ∈ x , y is computed by [figure omitted; refer to PDF]

Using Gauss-Legendre integration along the segments of boundary ^{Γ k s} , we have [figure omitted; refer to PDF] where _{N seg} is the number of segments of the boundary ^{Γ k s} , _{N gauss} is the number of Gauss points used along each segment, _{w m , n} are the corresponding Gauss weights, _{x m , n} is the n th Gaussian point on the m th segment of the boundary ^{Γ k s} , whose outward unit normal is denoted by _{n h} , the subscript " e " is either 1 or 2 as shown in Figure 4, and the superscript " e " indicates a domain restriction to element e .

The stiffness matrix K - associated with a smoothing domain is assembled by a similar process as in the FEM: [figure omitted; refer to PDF]

All entries in matrix _{B - I} in (11) with triangular meshes are constants over each smoothing domain; the stiffness matrix in (14) is therefore calculated by [figure omitted; refer to PDF]

4.2. Weak Formulation and Discretized Equation

We return to the two-dimensional body in Figure 1. Since the smoothed strain over smoothing domains is variationally consistent as proven in [26] and used by [27, 28], the assumed displacement ^{u h} and the smoothed strains ... - satisfy the "smoothed" Galerkin weak form: find ^{u h} ∈ V , for all δ ^{u h} ∈ _{V 0} such that [figure omitted; refer to PDF] with V = { u |" u ∈ ^{H 1} ( Ω \ _{Γ c} ) , u = u - on _{Γ u} , u discontinuous on _{Γ c} } and _{V 0} = { δ u |" δ u ∈ ^{H 1} ( Ω \ _{Γ c} ) , δ u = 0 on _{Γ u} , δ u discontinuous on _{Γ c} } .

Substituting the trial and test functions into (16), we finally obtain the familiar equation: [figure omitted; refer to PDF] where f is the nodal force vector that is identical to that in the standard phantomnode. The edge-based smoothed stiffness matrix K - for all subcells follows (15).

The smoothed stress ^{σ - h} is obtained in the same way from the as ^{... - h} in FEM, which is constant over a smoothing cell. In particular, for linear elastic problems, ^{σ - h} = D ^{... - h} is calculated on the level of the smoothing cell.

4.3. Crack Growth and Stress Intensity Factor

Fracture parameters such as mode I and mode I I stress intensity factors (SIFs) are determined using the domain form [29, 30] of the interaction integral [31]. All the finite elements within a radius of _{r d} = _{r k h e} from the crack-tip are used. Herein, _{h e} is the crack tip element size, and _{r k} is a scalar.

In this paper, crack growth is governed by the maximum hoop stress criterion [32, 33], which assumes that the crack will propagate from its tip in the direction _{θ c} , where the circumferential (hoop) stress _{σ θ θ} is maximum. The angle of crack propagation satisfies the following equation: [figure omitted; refer to PDF] Solving this equation, we have [figure omitted; refer to PDF] Once _{K I} and _{K I I} are known, (19) may be used to compute the direction of propagation.

5. Numerical Examples

In all numerical examples, we are not using near-tip enrichment; that is, only discontinuous enrichment is used. This means that the best convergence rate attainable is 1/2 in the _{H 1} norm and 1 in the _{L 2} norm ( ...AA; ( ^{h 1 / 2} ) and ...AA; ( h ) , respectively, where h is the mesh size).

5.1. Sheet with an Edge Crack under Uniaxial Tension

Consider a sheet under uniaxial tension as shown in Figure 8. The dimensions of the sheet are in unit of [mm]. The material parameters are Young's modulus E = 3 × ^{1 0 7} Pa, Poisson's ratio ν = 0.3 . The plane strain condition is assumed. The reference mode I SIF is given by [figure omitted; refer to PDF] where a = 0.5 is the crack length, b is the sheet width and F ( a / b ) is given by [figure omitted; refer to PDF]

Figure 8: Sheet with edge crack under tension.

[figure omitted; refer to PDF]

The strain energy and the error in the energy norm are defined as [figure omitted; refer to PDF] where the superscript "ref" denotes the exact or reference solution, and "num" denotes the numerical solution.

The results of the ES-Phantom node are compared with those of the standard Phantom-node using triangular meshes and the XFEM-T3(0t) (the "standard" XFEM formulation without tip enrichment that only employs the Heaviside enrichment of (2)). Figure 9 shows that the strain energy of the ES-Phantom node method is more accurate than both the original Phantom-node and the XFEM-T3(0t). The convergence rates in terms of the strain energy and the stress intensity factor _{K I} are depicted in Figures 10 and 12, respectively. Furthermore, Phantom-node using triangular meshes (Phantom-T3) is superior to XFEM-T3(0t) although they are equivalent to each other. We also have included two more figures comparing the phantom-node method to the tip-enriched XFEM in Figures 11 and 13 as a reference, although it would not be fair to compare a method that includes the asymptotic crack tip enrichment to a method that models the crack in a much simpler way.

Figure 9: Strain energy for the sheet with edge crack under tension.

[figure omitted; refer to PDF]

Figure 10: The convergence in the energy norm versus h (mesh size) for the sheet with an edge crack under tension.

[figure omitted; refer to PDF]

Figure 11: The convergence in the energy norm of XFEM versus h (mesh size) for the sheet with an edge crack under tension.

[figure omitted; refer to PDF]

Figure 12: The convergence in the stress intensity factor _{K I} versus h (mesh size) for the sheet with edge crack under tension.

[figure omitted; refer to PDF]

Figure 13: The convergence in the stress intensity factor _{K I} of XFEM versus h (mesh size) for the sheet with edge crack under tension.

[figure omitted; refer to PDF]

Note that the proposed method leads to a similar convergence rate to the standard XFEM and standard phantom-node, which is close to optimal (1/2) given the lack of tip enrichment. Also note that the error level of the proposed method is a fifth of an order of magnitude lower than the method compared with.

The computational efficiency in terms of the error in the energy norm and the relative error of _{K I} versus computation time ( s ) is compared for the ES-Phantom, the standard Phantom and the XFEM-T3(0t). The results are plotted in Figures 14 and 15, respectively. It is clear that the present method always produces higher computational efficiency, that is, accuracy to computational time ratio, compared to the other methods. The accuracy of the present method is approximate (1) (ES-Phantom10^-1.88 )/(Phantom10^-1.93 ) = 1.12 times as much as that of the standard Phantom, (ES-Phantom10^-1.88 )/(XFEM-T3(0t)10^-2.01 ) = 1.34 times of the XFEM-T3(0t) in term of error in energy norm; (2) (ES-Phantom10^-0.48 )/(Phantom10^-0.58 ) = 1.26 times as much as that of the standard Phantom and (ES-Phantom10^-0.48 )/(XFEM-T3(0t)10^-0.58 ) = 1.62 times of the XFEM-T3(0t) in term of relative error for _{K I} .

Figure 14: Computational efficiency of energy norm for the problem of a sheet with an edge crack under remote tension.

[figure omitted; refer to PDF]

Figure 15: Computational efficiency of mode I SIF _{K I} for the problem of a sheet with an edge crack under remote tension.

[figure omitted; refer to PDF]

5.2. Sheet with Edge-Crack under Shear

In this example, we consider the edge crack geometry subjected to a shear load as shown in Figure 16. The material parameters are Young's modulus E = 3 × ^{1 0 7} Pa and Poisson's ratio ν = 0.25 . The exact stress intensity factors for this load case are given by [31] [figure omitted; refer to PDF]

Figure 16: Sheet with edge crack under shear.

[figure omitted; refer to PDF]

The results shown in Figures 16, 17, 18, 19, 20, 21, and 22 show that ES-Phantom node results are more accurate than both those of the standard Phantom-node and the XFEM-T3(0t). ES-Phantom node maintains slight superconvergent solutions in the strain energy. Furthermore, Phantom-node using triangular meshes (Phantom-T3) is superior to XFEM-T3(0t) although they are equivalent to each other with respect to the convergence in energy norm and the stress intensity factor _{K I} . Figures 19, 21, and 23 again are shown as a reference for readers.

Figure 17: Strain energy for a sheet with an edge crack under shear.

[figure omitted; refer to PDF]

Figure 18: The convergence in the energy norm versus h (mesh size) for a sheet with an edge crack under shear.

[figure omitted; refer to PDF]

Figure 19: The convergence in the energy norm of XFEM versus h (mesh size) for a sheet with an edge crack under shear.

[figure omitted; refer to PDF]

Figure 20: The convergence in the stress intensity factor _{K I} versus h (mesh size) for a sheet with an edge crack under shear.

[figure omitted; refer to PDF]

Figure 21: The convergence in the stress intensity factor _{K I} of XFEM versus h (mesh size) for a sheet with an edge crack under shear.

[figure omitted; refer to PDF]

Figure 22: The convergence in the stress intensity factor _{K I I} versus h (mesh size) for sheet with edge crack under shear.

[figure omitted; refer to PDF]

Figure 23: The convergence in the stress intensity factor _{K I I} of XFEM versus h (mesh size) for sheet with edge crack under shear.

[figure omitted; refer to PDF]

5.3. Crack Growth Simulation in a Double Cantilever Beam

In this section, the ES-Phantom node with structured and unstructured meshes is used for crack grow simulation. The dimensions of the double cantilever beam Figure 24 are 6 mm × 2 mm and an initial pre-crack with length a = 2 mm is considered. Plane stress conditions are assumed with Young's modulus, E = 100 MPa as well as the Poisson ratio, ν = 3 , and the load P is taken to be unity. By symmetry, a crack on the mid-plane of the cantilever beam is dominated by pure mode I and the crack would propagate in a straight line. We also have included simulations with a "distorted" mesh and showed that the crack path does not change.

Figure 24: Double cantilever beam with an edge crack.

[figure omitted; refer to PDF]

The crack growth increment, Δ a , is taken so that the tip is always located at an element's edge and the crack growth is simulated for 10 steps. The domain is discretized with a structured and unstructured mesh of 2730 nodes. The crack path is simulated using both the proposed ES-Phantom node method and XFEM-T3(0t), and Figures 25 and 26 show the deformed shape of the double cantilever beam with the magnification factor of 25 × 10⁴ used to enable a clear description and the evolution of the crack path. The result shows that the crack path for an initial angle _{θ c} = 0 agrees with the published results [1].

(a) Deformed shape of the double cantilever beam (structured mesh) and (b) crack path simulated by ES-Phantom node method (structured mesh) after ten-step growing in which the filled circles are the new crack tip after each step.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(a) Deformed shape of the double cantilever beam (unstructured mesh) and (b) Crack path simulated by ES-Phantom node method (unstructured mesh) after ten-step growing in which the filled circles are the new crack tip after each step.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

Stress (a) _{σ x x} and (b) _{σ y y} contours in the sheet (structured mesh) after the crack propagates.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

Stress (a) _{σ x x} and (b) _{σ y y} contours in the sheet (unstructured mesh) after the crack propagates.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

6. Conclusions

A numerical Phantom-node method for analysis of two-linear elastic fracture problems was developed in framework of the ES-FEM to create the novel ES-Phantom node method. In this method, a cracked element is replaced by two superimposed elements and a set of additional phantom nodes. The two first examples were performed to investigate convergence rate in terms of strain energy and stress intensity factors. The results have shown that the ES-Phantom node is able to produce superconvergent solutions. Meanwhile, the last example has demonstrated the capability of the method to deal with the growing crack.

Future applications of this method may deal with the interactions among a large number of cracks in linear elastic solids with the purpose of obtaining the higher accuracy and efficiency in solving complicated crack interactions as shown in [34]. This will be studied in the future.

Acknowledgments

The authors gratefully acknowledge the support by the Deutscher Akademischer Austausch Dienst (DAAD). X. Zhuang acknowledges the supports from the NSFC (51109162), the National Basic Research Program of China (973 Program: 2011CB013800), and the Shanghai Pujiang Program (12PJ1409100).

References

[1] T. Belytschko, T. Black, "Elastic crack growth in finite elements with minimal remeshing," International Journal for Numerical Methods in Engineering , vol. 45, no. 5, pp. 601-620, 1999.

[2] A. Hansbo, P. Hansbo, "A finite element method for the simulation of strong and weak discontinuities in solid mechanics," Computer Methods in Applied Mechanics and Engineering , vol. 193, no. 33-35, pp. 3523-3540, 2004.

[3] J. Mergheim, E. Kuhl, P. Steinmann, "A finite element method for the computational modelling of cohesive cracks," International Journal for Numerical Methods in Engineering , vol. 63, no. 2, pp. 276-289, 2005.

[4] J.-H. Song, P. M. A. Areias, T. Belytschko, "A method for dynamic crack and shear band propagation with phantom nodes," International Journal for Numerical Methods in Engineering , vol. 67, no. 6, pp. 868-893, 2006.

[5] P. M. A. Areias, T. Belytschko, "A comment on the article "A finite element method for simulation of strong and weak discontinuities in solid mechanics"," Computer Methods in Applied Mechanics and Engineering , vol. 195, no. 9-12, pp. 1275-1276, 2006.

[6] J. Shi, D. Chopp, J. Lua, N. Sukumar, T. Belytschko, "Abaqus implementation of extended finite element method using a level set representation for three-dimensional fatigue crack growth and life predictions," Engineering Fracture Mechanics , vol. 77, no. 14, pp. 2840-2863, 2010.

[7] M. Duflot, "Industrial applications of XFEM for 3D crack propagation with Morfeo/Crack and Abaqus," in ECCOMAS Thematic Conference on XFEM, Cardiff, UK, June 2011.

[8] T. Menouillard, J. Réthoré, A. Combescure, H. Bung, "Efficient explicit time stepping for the eXtended Finite Element Method (X-FEM)," International Journal for Numerical Methods in Engineering , vol. 68, no. 9, pp. 911-939, 2006.

[9] T. Menouillard, J. Réthoré, N. Moës, A. Combescure, H. Bung, "Mass lumping strategies for X-FEM explicit dynamics: application to crack propagation," International Journal for Numerical Methods in Engineering , vol. 74, no. 3, pp. 447-474, 2008.

[10] T. Chau-Dinh, G. Zi, P.-S. Lee, T. Rabczuk, J.-H. Song, "Phantom-node method for shell models with arbitrary cracks," Computers and Structures , vol. 92-93, pp. 242-246, 2012.

[11] P. Laborde, J. Pommier, Y. Renard, M. Salaün, "High-order extended finite element method for cracked domains," International Journal for Numerical Methods in Engineering , vol. 64, no. 3, pp. 354-381, 2005.

[12] E. Béchet, H. Minnebo, N. Moës, B. Burgardt, "Improved implementation and robustness study of the X-FEM for stress analysis around cracks," International Journal for Numerical Methods in Engineering , vol. 64, no. 8, pp. 1033-1056, 2005.

[13] G. Ventura, R. Gracie, T. Belytschko, "Fast integration and weight function blending in the extended finite element method," International Journal for Numerical Methods in Engineering , vol. 77, no. 1, pp. 1-29, 2009.

[14] R. Gracie, H. Wang, T. Belytschko, "Blending in the extended finite element method by discontinuous Galerkin and assumed strain methods," International Journal for Numerical Methods in Engineering , vol. 74, no. 11, pp. 1645-1669, 2008.

[15] S. P. A. Bordas, T. Rabczuk, N.-X. Hung, V. P. Nguyen, S. Natarajan, T. Bog, D. M. Quan, N. V. Hiep, "Strain smoothing in FEM and XFEM," Computers and Structures , vol. 88, no. 23-24, pp. 1419-1443, 2010.

[16] T. Rabczuk, G. Zi, A. Gerstenberger, W. A. Wall, "A new crack tip element for the phantom-node method with arbitrary cohesive cracks," International Journal for Numerical Methods in Engineering , vol. 75, no. 5, pp. 577-599, 2008.

[17] G. R. Liu, K. Y. Dai, T. T. Nguyen, "A smoothed finite element method for mechanics problems," Computational Mechanics , vol. 39, no. 6, pp. 859-877, 2007.

[18] G. R. Liu, T. T. Nguyen, K. Y. Dai, K. Y. Lam, "Theoretical aspects of the smoothed finite element method (SFEM)," International Journal for Numerical Methods in Engineering , vol. 71, no. 8, pp. 902-930, 2007.

[19] N. Nguyen-Thanh, T. Rabczuk, H. Nguyen-Xuan, S. P. A. Bordas, "A smoothed finite element method for shell analysis," Computer Methods in Applied Mechanics and Engineering , vol. 198, no. 2, pp. 165-177, 2008.

[20] S. P. A. Bordas, S. Natarajan, "On the approximation in the smoothed finite element method (SFEM)," International Journal for Numerical Methods in Engineering , vol. 81, no. 5, pp. 660-670, 2010.

[21] T. Rabczuk, P. M. A. Areias, T. Belytschko, "A meshfree thin shell method for non-linear dynamic fracture," International Journal for Numerical Methods in Engineering , vol. 72, no. 5, pp. 524-548, 2007.

[22] T. Rabczuk, T. Belytschko, "Application of particle methods to static fracture of reinforced concrete structures," International Journal of Fracture , vol. 137, no. 1-4, pp. 19-49, 2006.

[23] T. Rabczuk, T. Belytschko, "A three-dimensional large deformation meshfree method for arbitrary evolving cracks," Computer Methods in Applied Mechanics and Engineering , vol. 196, no. 29-30, pp. 2777-2799, 2007.

[24] S. Bordas, P. V. Nguyen, C. Dunant, A. Guidoum, H. Nguyen-Dang, "An extended finite element library," International Journal for Numerical Methods in Engineering , vol. 71, no. 6, pp. 703-732, 2007.

[25] G. R. Liu, T. Nguyen-Thoi, K. Y. Lam, "An edge-based smoothed finite element method (ES-FEM) for static, free and forced vibration analyses of solids," Journal of Sound and Vibration , vol. 320, no. 4-5, pp. 1100-1130, 2009.

[26] G. R. Liu, "A G space theory and a weakened weak (W² ) form for a unified formulation of compatible and incompatible methods: part I theory," International Journal for Numerical Methods in Engineering , vol. 81, no. 9, pp. 1093-1126, 2010.

[27] L. Chen, T. Rabczuk, S. P. A. Bordas, G. R. Liu, K. Y. Zeng, P. Kerfriden, "Extended finite element method with edge-based strain smoothing (ESm-XFEM) for linear elastic crack growth," Computer Methods in Applied Mechanics and Engineering , vol. 209-212, pp. 250-265, 2012.

[28] N. Vu-Bac, H. Nguyen-Xuan, L. Chen, S. Bordas, P. Kerfriden, R. N. Simpson, G. R. Liu, T. Rabczuk, "A node-based smoothed extended finite element method (NS-XFEM) for fracture analysis," CMES: Computer Modeling in Engineering and Sciences , vol. 73, no. 4, pp. 331-355, 2011.

[29] F. Z. Li, C. F. Shih, A. Needleman, "A comparison of methods for calculating energy release rates," Engineering Fracture Mechanics , vol. 21, no. 2, pp. 405-421, 1985.

[30] B. Moran, C. F. Shih, "Crack tip and associated domain integrals from momentum and energy balance," Engineering Fracture Mechanics , vol. 27, no. 6, pp. 615-642, 1987.

[31] S. S. Wang, J. F. Yau, H. T. Corten, "A mixed-mode crack analysis of rectilinear anisotropic solids using conservation laws of elasticity," International Journal of Fracture , vol. 16, no. 3, pp. 247-259, 1980.

[32] F. Erdogan, G. Sih, "On the crack extension in sheets under plane loading and transverse shear," Journal Basic Engineering , vol. 85, no. 6, pp. 519-527, 1963.

[33] A. Menk, S. Bordas, "Crack growth calculations in solder joints based on microstructural phenomena with X-FEM," Computational Materials Science , vol. 50, no. 3, pp. 1145-1156, 2011.

[34] D. F. Li, C. F. Li, S. Q. Shu, Z. X. Wang, J. Lu, "A fast and accurate analysis of the interacting cracks in linear elastic solids," International Journal of Fracture , vol. 151, pp. 169-185, 2008.