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.

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 (𝐊^𝑢𝑎 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 [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 𝐛, 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 𝑓(𝐱)<0 and the other active part (Ω₂), 𝐮²(𝐱), which holds for 𝑓(𝐱)>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.

The displacement field within an element Ω_𝑒 in Figure 2 is rewritten as [4] ∀𝐱∈Ω_𝑒,𝐮(𝐱)=_𝐼∈𝑆1𝐮¹_𝐼𝑁_𝐼(𝐱)_𝐮1(𝐱)𝐻−𝑓(𝐱)+_𝐼∈𝑆2𝐮²_𝐼𝑁_𝐼(𝐱)_𝐮2(𝐱)𝐻𝑓(𝐱), (1) 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.

Figure 2 

The decomposition of a cracked element into two superimposed elements.

𝐻 is the Heaviside function given in [1, 21–24] and defined by 𝐻(𝑥)=1,𝑥>0,0,𝑥⩽0. (2)

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: ∀𝐱∈Ω_𝑒,𝝐(𝐱)=_𝐼∈𝑆1𝐁_𝐼(𝐱)𝐮¹_{𝐼𝝐1(𝐱)}𝐻−𝑓(𝐱)+_𝐼∈𝑆2𝐁_𝐼(𝐱)𝐮²_{𝐼𝝐2(𝐱)}𝐻𝑓(𝐱), (3) where 𝐁_𝐼 is the standard strain-displacement matrix.

The jump in the displacement field across the crack is calculated by ⟦𝐮(𝐱)⟧=𝐮¹(𝐱)−𝐮²(𝐱)onΓ_𝑐;𝐱𝐼isaphantomnodeinelement1if𝑓_𝐼𝐱>0,element2if𝑓_𝐼<0. (4)

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. Ω^(𝑘) satisfies the conditions Ω=⋃^𝑁𝑒_𝑘=1Ω^(𝑘) 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.

Figure 3 

Construction of edge-based strain smoothing domains.

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 5 

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

Figure 6 

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

Figure 7 

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

Distribution of the stress components 𝜎_𝑥𝑥 and 𝜎_𝑦𝑦 for and unstructured mesh are shown in Figures 27 and 28, respectively.

Introducing the edge-based smoothing operation, the compatible strain 𝝐=∇_𝑠𝐮^ℎ_𝑘 is smoothed over cell Ω^(𝑘) associated with edge 𝑘 as follows: 𝝐_𝑘=_Ω(𝑘)𝝐(𝐱)Φ_𝑘(𝐱)dΩ=_Ω(𝑘)∇_𝑠𝐮^ℎ(𝐱)Φ_𝑘(𝐱)dΩ, (5) where Φ_𝑘 is a given normalized smoothing function that satisfies _Ω(𝑘)Φ_𝑘(𝐱)dΩ=1. (6)

Using the following constant smoothing function 1Φ=𝐴^(𝑘),𝐱∈Ω^(𝑘)0,𝐱∉Ω^(𝑘),𝝐_𝑘=1𝐴^{(𝑘)Ω(𝑘)}∇_𝑠𝐮^ℎ(𝐱)1dΩ=𝐴^{(𝑘)Γ(𝑘)}𝐋_𝑛𝐮^ℎ(𝐱)dΓ, (7) where 𝐴^(𝑘)=∫_Ω(𝑘)dΩ 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 𝐋_𝑛=𝑛_𝑥00𝑛_𝑦𝑛_𝑦𝑛_𝑥. (8)

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): 𝐮^ℎ(𝐱)=_{𝐼∈𝑆es-pht1}𝐮¹_𝐼𝑁_𝐼(𝐱)_𝐮1(𝐱)𝐻𝑓(𝐱)+_{𝐼∈𝑆es-pht2}𝐮²_𝐼𝑁_𝐼(𝐱)_𝐮2(𝐱)𝐻−𝑓(𝐱), (9) where 𝑆^es-pht₁ and 𝑆^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 Ω^(𝑘) (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 (Ω₁^(𝑘)) = [1,2,3,4]; nodes of smoothing domain 2 (Ω₂^(𝑘)) = [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,2,4]; nodes of smoothing domain 2 (Ω₂^(𝑘)) = [1,2,4].

Using the strain smoothing operation, the smoothed strain associated with edge 𝑘 created from the displacement approximation in (9) can be rewritten as𝝐_𝑘=_{𝐼∈𝑆es-pht1}𝐁_𝐼𝐱_𝑘𝐮¹_{𝐼𝝐1(𝐱)}𝐻−𝑓(𝐱)+_{𝐼∈𝑆es-pht2}𝐁_𝐼𝐱_𝑘𝐮²_{𝐼𝝐2(𝐱)}𝐻𝑓(𝐱), (10) where 𝐁_𝐼(𝐱_𝑘) is the smoothed strain gradient matrix for the standard ES-FEM part. Those matrices write as follows 𝐁_𝐼𝐱_𝑘=𝑏_𝐼𝑥𝐱_𝑘00𝑏_𝐼𝑦𝐱_𝑘𝑏_𝐼𝑦𝐱_𝑘𝑏_𝐼𝑥𝐱_𝑘. (11)

In (11), 𝑏_𝐼ℎ(𝐱_𝑘), ℎ∈𝑥,𝑦 is computed by 𝑏_𝐼ℎ𝐱_𝑘=1𝐴^𝑠_𝑘Γ^𝑠_𝑘𝑛_ℎ(𝐱)𝑁_𝑖(𝐱)𝐻(−1)^𝑒𝑓(𝐱)dΓ. (12)

Using Gauss-Legendre integration along the segments of boundary Γ^𝑠_𝑘, we have 𝑏_𝐼ℎ=1𝐴^𝑠_𝑘^𝑁seg_𝑚=1^𝑁gauss_𝑛=1𝑤_𝑚,𝑛𝑁_𝑖𝐱_𝑚,𝑛𝐻×(−1)^𝑒𝑓𝐱_𝑚,𝑛𝑛_ℎ𝐱_𝑚,𝑛, (13) where 𝑁_seg is the number of segments of the boundary Γ^𝑠_𝑘, 𝑁_gauss 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: 𝐊_𝐼𝐽=^𝑁𝑠_𝑘=1𝐊^𝑠_{𝐼𝐽,𝑘}=^𝑁𝑠_{𝑘=1Ω𝑠𝑘}𝐁^𝑢_𝐼^𝑇𝐃𝐁^𝑢_𝐽dΩ. (14)

All entries in matrix 𝐁_𝐼 in (11) with triangular meshes are constants over each smoothing domain; the stiffness matrix in (14) is therefore calculated by 𝐊_𝐼𝐽=^𝑁𝑠_𝑘=1𝐊^𝑠_{𝐼𝐽,𝑘}=^𝑁𝑠_𝑘=1𝐁^𝑢_𝐼^𝑇𝐃𝐁^𝑢_𝐽𝐴^𝑠_𝑘. (15)

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 𝐮^ℎ and the smoothed strains 𝝐 satisfy the “smoothed” Galerkin weak form: find 𝐮^ℎ∈𝑉, for all 𝛿𝐮^ℎ∈𝑉₀ such that _Ω𝛿𝝐𝐮^ℎ𝑇𝐃𝝐𝐮^ℎdΩ−_Ω𝛿𝐮^ℎ𝑇−𝐛dΩ_Γ𝛿𝐮^ℎ𝑇𝐭_ΓdΓ=0 (16) with 𝑉={𝐮∣𝐮∈𝐻¹(Ω⧵Γ_𝑐), 𝐮=𝐮onΓ_𝑢, 𝐮 discontinuous on Γ_𝑐} and 𝑉₀={𝛿𝐮∣𝛿𝐮∈𝐻¹(Ω⧵Γ_𝑐), 𝛿𝐮=0onΓ_𝑢, 𝛿𝐮 discontinuous on Γ_𝑐}.

Substituting the trial and test functions into (16), we finally obtain the familiar equation: 𝐊𝐝=𝐟, (17) 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 [31]. 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: 𝐾_𝐼sin𝜃_𝑐+𝐾_𝐼𝐼3cos𝜃_𝑐−1=0. (18) Solving this equation, we have 𝜃_𝑐𝐾=2arctan−2_𝐼𝐼/𝐾_𝐼1+1+8𝐾_𝐼𝐼/𝐾_𝐼². (19) 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 (𝒪(ℎ^1/2) 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 𝐸=3×10⁷ Pa, Poisson’s ratio 𝜈=0.3. The plane strain condition is assumed. The reference mode 𝐼 SIF is given by 𝐾^exact_𝐼𝑎=𝐹𝑏𝜎√𝑎𝜋=1.6118Pa√mm, (20) where 𝑎=0.5 is the crack length, 𝑏 is the sheet width and 𝐹(𝑎/𝑏) is given by 𝐹𝑎𝑏𝑎=1.12−0.231𝑏𝑎+10.55𝑏²𝑎−21.72𝑏³𝑎+30.39𝑏⁴. (21)

Figure 8 

Sheet with edge crack under tension.

The strain energy and the error in the energy norm are defined as 𝐸_(Ω)=12_Ω𝝐^𝑇𝑒𝐃𝝐dΩ,_𝑒=𝐸^num_(Ω)−𝐸^ref_(Ω)𝐸^ref_(Ω)^1/2,𝑒_𝑘=𝐾^num_sif−𝐾^ref_sif𝐾^ref_sif^1/2×100%,sif=𝐼,𝐼𝐼, (22) 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.

Figure 9 

Strain energy for the sheet with edge crack under tension.

Figure 10 

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

Figure 11 

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

Figure 12 

The convergence in the stress intensity factor 𝐾𝐼 versus ℎ (mesh size) for the sheet with edge crack under tension.

Figure 13 

The convergence in the stress intensity factor 𝐾𝐼 of XFEM versus ℎ (mesh size) for the sheet with edge crack under tension.

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 𝐾_𝐼.

Figure 14 

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

Figure 15 

Computational efficiency of mode 𝐼 SIF 𝐾𝐼 for the problem of a sheet with an edge crack under remote tension.

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 𝐸=3×10⁷ Pa and Poisson’s ratio 𝜈=0.25. The exact stress intensity factors for this load case are given by [31] 𝐾_𝐼=34.0Pa√mm;𝐾_𝐼𝐼=4.55Pa√mm. (23)

Figure 16 

Sheet with edge crack under shear.

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.

Figure 17 

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

Figure 18 

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

Figure 19 

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

Figure 20 

The convergence in the stress intensity factor 𝐾𝐼 versus ℎ (mesh size) for a sheet with an edge crack under shear.

Figure 21 

The convergence in the stress intensity factor 𝐾𝐼 of XFEM versus ℎ (mesh size) for a sheet with an edge crack under shear.

Figure 22 

The convergence in the stress intensity factor 𝐾𝐼𝐼 versus ℎ (mesh size) for sheet with edge crack under shear.

Figure 23 

The convergence in the stress intensity factor 𝐾𝐼𝐼 of XFEM versus ℎ (mesh size) for sheet with edge crack under shear.

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 𝑎=2 mm is considered. Plane stress conditions are assumed with Young’s modulus, 𝐸=100 MPa as well as the Poisson ratio, 𝜈=3, 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.

Figure 24 

Double cantilever beam with an edge crack.

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 × 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 𝜃_𝑐=0 agrees with the published results [1].

(a)

(b)

Figure 25 

(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)

(b)

Figure 26 

(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)

(b)

Figure 27 

Stress (a) 𝜎𝑥𝑥 and (b) 𝜎𝑦𝑦 contours in the sheet (structured mesh) after the crack propagates.

(a)

(b)

Figure 28 

Stress (a) 𝜎𝑥𝑥 and (b) 𝜎𝑦𝑦 contours in the sheet (unstructured mesh) after the crack propagates.

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).