- About this Journal
- Abstracting and Indexing
- Aims and Scope
- Annual Issues
- Article Processing Charges
- Articles in Press
- Author Guidelines
- Bibliographic Information
- Citations to this Journal
- Contact Information
- Editorial Board
- Editorial Workflow
- Free eTOC Alerts
- Publication Ethics
- Reviewers Acknowledgment
- Submit a Manuscript
- Subscription Information
- Table of Contents
Journal of Applied Mathematics
Volume 2013 (2013), Article ID 978026, 12 pages
A Phantom-Node Method with Edge-Based Strain Smoothing for Linear Elastic Fracture Mechanics
1Department of Geotechnical Engineering, Tongji University, Shanghai 200092, China
2Institute of Structural Mechanics, Bauhaus-University Weimar, Marienstraße 15, D-99423 Weimar, Germany
3Division of Computational Mechanics, Ton Duc Thang University, Ho Chi Minh City 70000, Vietnam
4School of Chemistry, Physics and Mechanical Engineering, Queensland University of Technology, Brisbane, QLD 4001, Australia
5School of Engineering, Institute of Mechanics and Advanced Materials, Theoretical and Computational Mechanics, Cardiff University, Wales CF24 3AA, UK
6Department of Civil, Environmental & Architectural Engineering, Korea University, 5 Ga 1, An-Am Dong, Sung-Buk Gu, Seoul 136-701, Republic of Korea
7Aerospace Systems Ohio Eminent Scholar, University of Cincinnati, Cincinnati, OH 45221-0070, USA
8School 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
Academic Editor: Song Cen
Copyright © 2013 N. Vu-Bac et al. 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.
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.
The extended finite element method (XFEM)  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 , subsequently implemented by  in a static setting and by  in a dynamic setting. It was shown by [4, 5] that the method proposed by  is identical to a step-enriched XFEM;  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  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  while an additional plug-in  is required to model crack growth using XFEM. (2)No mixed terms ( and ) 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 .
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  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  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 , external applied traction on boundary , and displacement boundary conditions on 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 , , which holds for and the other active part , , which holds for . The two parts of the model do not share nodes, and therefore they displace (deform) independently. Areias and Belytschko  demonstrated that the Hansbo and Hansbo  formulation is equivalent to the XFEM formulation relying on discontinuous enrichment with the Heaviside function.
The displacement field within an element in Figure 2 is rewritten as  where and 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; is the finite element shape function associated with node , while and are nodal displacements of original nodes in superimposed element 1 and 2, respectively.
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: where is the standard strain-displacement matrix.
The jump in the displacement field across the crack is calculated by
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 , the domain is partitioned into a set of nonoverlapping no-gap smoothing domains constructed using element edges of the triangular elements. satisfies the conditions and , for all , in which 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 , and the smoothing domain for a boundary edge are illustrated.
Introducing the edge-based smoothing operation, the compatible strain is smoothed over cell associated with edge as follows: where is a given normalized smoothing function that satisfies
Using the following constant smoothing function where is the area of the smoothing domain , is the boundary of the smoothing domain , and is the outward unit normal matrix which can be expressed as
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): where and 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 (DEFG) and boundary smoothing domain (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 () = ; nodes of smoothing domain 2 () = .
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 () = ; nodes of smoothing domain 2 () = .
Using the strain smoothing operation, the smoothed strain associated with edge created from the displacement approximation in (9) can be rewritten as where is the smoothed strain gradient matrix for the standard ES-FEM part. Those matrices write as follows
In (11), , is computed by
Using Gauss-Legendre integration along the segments of boundary , we have where is the number of segments of the boundary , is the number of Gauss points used along each segment, are the corresponding Gauss weights, is the th Gaussian point on the th segment of the boundary , whose outward unit normal is denoted by , the subscript “” is either 1 or 2 as shown in Figure 4, and the superscript “” indicates a domain restriction to element .
The stiffness matrix associated with a smoothing domain is assembled by a similar process as in the FEM:
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  and used by [27, 28], the assumed displacement and the smoothed strains satisfy the “smoothed” Galerkin weak form: find , for all such that with , , discontinuous on and , , discontinuous on .
Substituting the trial and test functions into (16), we finally obtain the familiar equation: where is the nodal force vector that is identical to that in the standard phantomnode. The edge-based smoothed stiffness matrix for all subcells follows (15).
The smoothed stress is obtained in the same way from the as in FEM, which is constant over a smoothing cell. In particular, for linear elastic problems, is calculated on the level of the smoothing cell.
4.3. Crack Growth and Stress Intensity Factor
Fracture parameters such as mode and mode stress intensity factors (SIFs) are determined using the domain form [29, 30] of the interaction integral . All the finite elements within a radius of from the crack-tip are used. Herein, is the crack tip element size, and 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 , where the circumferential (hoop) stress is maximum. The angle of crack propagation satisfies the following equation: Solving this equation, we have Once and 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 norm and 1 in the norm ( and , respectively, where 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 Pa, Poisson’s ratio . The plane strain condition is assumed. The reference mode SIF is given by where is the crack length, is the sheet width and is given by
The strain energy and the error in the energy norm are defined as 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 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.
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 versus computation time () 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 .
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 Pa and Poisson’s ratio . The exact stress intensity factors for this load case are given by 
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 . Figures 19, 21, and 23 again are shown as a reference for readers.
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 mm is considered. Plane stress conditions are assumed with Young’s modulus, MPa as well as the Poisson ratio, , and the load is taken to be unity. By symmetry, a crack on the mid-plane of the cantilever beam is dominated by pure mode 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.
The crack growth increment, , 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 × 104 used to enable a clear description and the evolution of the crack path. The result shows that the crack path for an initial angle agrees with the published results .
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 . This will be studied in the future.
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).
- T. Belytschko and 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.
- A. Hansbo and 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.
- J. Mergheim, E. Kuhl, and 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.
- J.-H. Song, P. M. A. Areias, and 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.
- P. M. A. Areias and 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.
- J. Shi, D. Chopp, J. Lua, N. Sukumar, and 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.
- 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.
- T. Menouillard, J. Réthoré, A. Combescure, and 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.
- T. Menouillard, J. Réthoré, N. Moës, A. Combescure, and 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.
- T. Chau-Dinh, G. Zi, P.-S. Lee, T. Rabczuk, and J.-H. Song, “Phantom-node method for shell models with arbitrary cracks,” Computers and Structures, vol. 92-93, pp. 242–246, 2012.
- P. Laborde, J. Pommier, Y. Renard, and 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.
- E. Béchet, H. Minnebo, N. Moës, and 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.
- G. Ventura, R. Gracie, and 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.
- R. Gracie, H. Wang, and 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.
- S. P. A. Bordas, T. Rabczuk, N.-X. Hung et al., “Strain smoothing in FEM and XFEM,” Computers and Structures, vol. 88, no. 23-24, pp. 1419–1443, 2010.
- T. Rabczuk, G. Zi, A. Gerstenberger, and 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.
- G. R. Liu, K. Y. Dai, and T. T. Nguyen, “A smoothed finite element method for mechanics problems,” Computational Mechanics, vol. 39, no. 6, pp. 859–877, 2007.
- G. R. Liu, T. T. Nguyen, K. Y. Dai, and 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.
- N. Nguyen-Thanh, T. Rabczuk, H. Nguyen-Xuan, and 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.
- S. P. A. Bordas and 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.
- T. Rabczuk, P. M. A. Areias, and 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.
- T. Rabczuk and 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.
- T. Rabczuk and 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.
- S. Bordas, P. V. Nguyen, C. Dunant, A. Guidoum, and H. Nguyen-Dang, “An extended finite element library,” International Journal for Numerical Methods in Engineering, vol. 71, no. 6, pp. 703–732, 2007.
- G. R. Liu, T. Nguyen-Thoi, and 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.
- G. R. Liu, “A G space theory and a weakened weak (W2) 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.
- L. Chen, T. Rabczuk, S. P. A. Bordas, G. R. Liu, K. Y. Zeng, and 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.
- N. Vu-Bac, H. Nguyen-Xuan, L. Chen et al., “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.
- F. Z. Li, C. F. Shih, and A. Needleman, “A comparison of methods for calculating energy release rates,” Engineering Fracture Mechanics, vol. 21, no. 2, pp. 405–421, 1985.
- B. Moran and 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.
- S. S. Wang, J. F. Yau, and 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.
- F. Erdogan and 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.
- A. Menk and 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.
- D. F. Li, C. F. Li, S. Q. Shu, Z. X. Wang, and 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.