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

Authors: N. Vu-Bac; H. Nguyen-Xuan; L. Chen; C. K. Lee; G. Zi; X. Zhuang; G. R. Liu; T. Rabczuk
Journal: Journal of Applied Mathematics
 (2013)

---

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


---

## Body

## 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 (Kua and Kau) 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 Section2, 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), u1(x), which holds for f(x)<0 and the other active part (Ω2), u2(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.The displacement field within an elementΩe in Figure 2 is rewritten as [4]
(1)∀x∈Ωe,u(x)=∑I∈S1‍uI1NI(x)︸u1(x)H(-f(x))+∑I∈S2‍uI2NI(x)︸u2(x)H(f(x)),
where S1 and S2 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; NI is the finite element shape function associated with node I, while uI1 and uI2 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.H is the Heaviside function given in [1, 21–24] and defined by
(2)H(x)={1,x>0,0,x⩽0.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:(3)∀x∈Ωe,ϵ(x)=∑I∈S1‍BI(x)uI1︸ϵ1(x)H(-f(x))+∑I∈S2‍BI(x)uI2︸ϵ2(x)H(f(x)),
where BI is the standard strain-displacement matrix.The jump in the displacement field across the crack is calculated by(4)⟦u(x)⟧=u1(x)-u2(x)onΓc;Iisaphantomnodein{element1iff(xI)>0,element2iff(xI)<0.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=1Ne‍Ω(k) and Ω(i)∩Ω(j)=∅, for  all i≠j, in which Ne 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.Numerical integration is implemented on chosen Gauss points as illustrated in Figures6 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σxx and σyy for and unstructured mesh are shown in Figures 27 and 28, respectively.Introducing the edge-based smoothing operation, the compatible strainϵ=∇sukh is smoothed over cell Ω(k) associated with edge k as follows:
(5)ϵ-k=∫Ω(k)ϵ(x)Φk(x)dΩ=∫Ω(k)∇suh(x)Φk(x)dΩ,
where Φk is a given normalized smoothing function that satisfies
(6)∫Ω(k)Φk(x)dΩ=1.Using the following constant smoothing function(7)Φ={1A(k),x∈Ω(k)0,x∉Ω(k),ϵ-k=1A(k)∫Ω(k)∇suh(x)dΩ=1A(k)∫Γ(k)Lnuh(x)dΓ,
where A(k)=∫Ω(k)dΩ is the area of the smoothing domain Ω(k),  Γ(k) is the boundary of the smoothing domain Ω(k), and Ln is the outward unit normal matrix which can be expressed as
(8)Ln=[nx00nynynx].


## 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=1Ne‍Ω(k) and Ω(i)∩Ω(j)=∅, for  all i≠j, in which Ne 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.Numerical integration is implemented on chosen Gauss points as illustrated in Figures6 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σxx and σyy for and unstructured mesh are shown in Figures 27 and 28, respectively.Introducing the edge-based smoothing operation, the compatible strainϵ=∇sukh is smoothed over cell Ω(k) associated with edge k as follows:
(5)ϵ-k=∫Ω(k)ϵ(x)Φk(x)dΩ=∫Ω(k)∇suh(x)Φk(x)dΩ,
where Φk is a given normalized smoothing function that satisfies
(6)∫Ω(k)Φk(x)dΩ=1.Using the following constant smoothing function(7)Φ={1A(k),x∈Ω(k)0,x∉Ω(k),ϵ-k=1A(k)∫Ω(k)∇suh(x)dΩ=1A(k)∫Γ(k)Lnuh(x)dΓ,
where A(k)=∫Ω(k)dΩ is the area of the smoothing domain Ω(k),  Γ(k) is the boundary of the smoothing domain Ω(k), and Ln is the outward unit normal matrix which can be expressed as
(8)Ln=[nx00nynynx].


## 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):
(9)uh(x)=∑I∈S1es-pht‍uI1NI(x)︸u1(x)H(f(x))+∑I∈S2es-pht‍uI2NI(x)︸u2(x)H(-f(x)),
where S1es-pht and S2es-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 Figure4: 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 Figure5 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 edgek created from the displacement approximation in (9) can be rewritten as(10)ϵ-k=∑I∈S1es-pht‍B-I(xk)uI1︸ϵ-1(x)H(-f(x))+∑I∈S2es-pht‍B-I(xk)uI2︸ϵ-2(x)H(f(x)),
where B-I(xk) is the smoothed strain gradient matrix for the standard ES-FEM part. Those matrices write as follows
(11)B-I(xk)=[b-Ix(xk)00b-Iy(xk)b-Iy(xk)b-Ix(xk)].In (11), b-Ih(xk), h∈x,y is computed by
(12)b-Ih(xk)=1Aks∫Γksnh(x)Ni(x)H((-1)ef(x))dΓ.Using Gauss-Legendre integration along the segments of boundaryΓks, we have
(13)b-Ih=1Aks∑m=1Nseg‍[∑n=1Ngauss‍wm,nNi(xm,n)H×((-1)ef(xm,n))nh(xm,n)∑n=1Ngauss‍],
where Nseg is the number of segments of the boundary Γks, Ngauss is the number of Gauss points used along each segment, wm,n are the corresponding Gauss weights, xm,n is the nth Gaussian point on the mth segment of the boundary Γks, whose outward unit normal is denoted by nh, 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 matrixK- associated with a smoothing domain is assembled by a similar process as in the FEM:
(14)K-IJ=∑k=1Ns‍K-IJ,ks=∑k=1Ns‍∫Ωks(B-Iu)TDB-JudΩ.All entries in matrixB-I in (11) with triangular meshes are constants over each smoothing domain; the stiffness matrix in (14) is therefore calculated by
(15)K-IJ=∑k=1Ns‍K-IJ,ks=∑k=1Ns‍(B-Iu)TDB-JuAks.


### 4.2. Weak Formulation and Discretized Equation

We return to the two-dimensional body in Figure1. Since the smoothed strain over smoothing domains is variationally consistent as proven in [26] and used by [27, 28], the assumed displacement uh and the smoothed strains ϵ- satisfy the “smoothed” Galerkin weak form: find uh∈V, for all δuh∈V0 such that
(16)∫Ωδ(ϵ-(uh))TD(ϵ-(uh))dΩ-∫Ω(δuh)TbdΩ-∫Γ(δuh)TtΓdΓ=0
with V={u∣u∈H1(Ω∖Γc), u=u-onΓu, u discontinuous on Γc} and V0={δu∣δu∈H1(Ω∖Γc), δu=0onΓu, δu discontinuous on Γc}.Substituting the trial and test functions into (16), we finally obtain the familiar equation:
(17)K-d-=f,
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 modeI and mode II 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 rd=rkhe from the crack-tip are used. Herein, he is the crack tip element size, and rk 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:
(18)KIsin(θc)+KII(3cos(θc)-1)=0.
Solving this equation, we have
(19)θc=2arctan[-2(KII/KI)1+1+8(KII/KI)2].
Once KI and KII are known, (19) may be used to compute the direction of propagation.


## 4.1. Displacement and Strain Field

The approximation of the displacement field is written similarly to (1):
(9)uh(x)=∑I∈S1es-pht‍uI1NI(x)︸u1(x)H(f(x))+∑I∈S2es-pht‍uI2NI(x)︸u2(x)H(-f(x)),
where S1es-pht and S2es-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 Figure4: 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 Figure5 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 edgek created from the displacement approximation in (9) can be rewritten as(10)ϵ-k=∑I∈S1es-pht‍B-I(xk)uI1︸ϵ-1(x)H(-f(x))+∑I∈S2es-pht‍B-I(xk)uI2︸ϵ-2(x)H(f(x)),
where B-I(xk) is the smoothed strain gradient matrix for the standard ES-FEM part. Those matrices write as follows
(11)B-I(xk)=[b-Ix(xk)00b-Iy(xk)b-Iy(xk)b-Ix(xk)].In (11), b-Ih(xk), h∈x,y is computed by
(12)b-Ih(xk)=1Aks∫Γksnh(x)Ni(x)H((-1)ef(x))dΓ.Using Gauss-Legendre integration along the segments of boundaryΓks, we have
(13)b-Ih=1Aks∑m=1Nseg‍[∑n=1Ngauss‍wm,nNi(xm,n)H×((-1)ef(xm,n))nh(xm,n)∑n=1Ngauss‍],
where Nseg is the number of segments of the boundary Γks, Ngauss is the number of Gauss points used along each segment, wm,n are the corresponding Gauss weights, xm,n is the nth Gaussian point on the mth segment of the boundary Γks, whose outward unit normal is denoted by nh, 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 matrixK- associated with a smoothing domain is assembled by a similar process as in the FEM:
(14)K-IJ=∑k=1Ns‍K-IJ,ks=∑k=1Ns‍∫Ωks(B-Iu)TDB-JudΩ.All entries in matrixB-I in (11) with triangular meshes are constants over each smoothing domain; the stiffness matrix in (14) is therefore calculated by
(15)K-IJ=∑k=1Ns‍K-IJ,ks=∑k=1Ns‍(B-Iu)TDB-JuAks.


## 4.2. Weak Formulation and Discretized Equation

We return to the two-dimensional body in Figure1. Since the smoothed strain over smoothing domains is variationally consistent as proven in [26] and used by [27, 28], the assumed displacement uh and the smoothed strains ϵ- satisfy the “smoothed” Galerkin weak form: find uh∈V, for all δuh∈V0 such that
(16)∫Ωδ(ϵ-(uh))TD(ϵ-(uh))dΩ-∫Ω(δuh)TbdΩ-∫Γ(δuh)TtΓdΓ=0
with V={u∣u∈H1(Ω∖Γc), u=u-onΓu, u discontinuous on Γc} and V0={δu∣δu∈H1(Ω∖Γc), δu=0onΓu, δu discontinuous on Γc}.Substituting the trial and test functions into (16), we finally obtain the familiar equation:
(17)K-d-=f,
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 modeI and mode II 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 rd=rkhe from the crack-tip are used. Herein, he is the crack tip element size, and rk 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:
(18)KIsin(θc)+KII(3cos(θc)-1)=0.
Solving this equation, we have
(19)θc=2arctan[-2(KII/KI)1+1+8(KII/KI)2].
Once KI and KII 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 theH1 norm and 1 in the L2 norm (𝒪(h1/2) and 𝒪(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 Figure8. The dimensions of the sheet are in unit of [mm]. The material parameters are Young’s modulus E=3×107 Pa, Poisson’s ratio ν=0.3. The plane strain condition is assumed. The reference mode I SIF is given by
(20)KIexact=F(ab)σaπ=1.6118Pamm,
where a=0.5 is the crack length, b is the sheet width and F(a/b) is given by
(21)F(ab)=1.12-0.231(ab)+10.55(ab)2-21.72(ab)3+30.39(ab)4.Figure 8
Sheet with edge crack under tension.The strain energy and the error in the energy norm are defined as(22)E(Ω)=12∫ΩϵTDϵdΩ,ee=|E(Ω)num-E(Ω)refE(Ω)ref|1/2,ek=|Ksifnum-KsifrefKsifref|1/2×100%,sif=I,II,
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 KI 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 versush (mesh size) for the sheet with an edge crack under tension.Figure 11
The convergence in the energy norm of XFEM versush (mesh size) for the sheet with an edge crack under tension.Figure 12
The convergence in the stress intensity factorKI versus h (mesh size) for the sheet with edge crack under tension.Figure 13
The convergence in the stress intensity factorKI of XFEM versus h (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 ofKI 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 KI.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 modeI SIF KI 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 Figure16. The material parameters are Young’s modulus E=3×107 Pa and Poisson’s ratio ν=0.25. The exact stress intensity factors for this load case are given by [31]
(23)KI=34.0Pamm;KII=4.55Pamm.Figure 16
Sheet with edge crack under shear.The results shown in Figures16, 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 KI. 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 versush (mesh size) for a sheet with an edge crack under shear.Figure 19
The convergence in the energy norm of XFEM versush (mesh size) for a sheet with an edge crack under shear.Figure 20
The convergence in the stress intensity factorKI versus h (mesh size) for a sheet with an edge crack under shear.Figure 21
The convergence in the stress intensity factorKI of XFEM versus h (mesh size) for a sheet with an edge crack under shear.Figure 22
The convergence in the stress intensity factorKII versus h (mesh size) for sheet with edge crack under shear.Figure 23
The convergence in the stress intensity factorKII of XFEM versus h (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 Figure24 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.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 × 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 θ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)
(b)(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)Stress (a)σxx and (b) σyy contours in the sheet (structured mesh) after the crack propagates.
(a)
(b)Stress (a)σxx and (b) σyy contours in the sheet (unstructured mesh) after the crack propagates.
(a)
(b)


## 5.1. Sheet with an Edge Crack under Uniaxial Tension

Consider a sheet under uniaxial tension as shown in Figure8. The dimensions of the sheet are in unit of [mm]. The material parameters are Young’s modulus E=3×107 Pa, Poisson’s ratio ν=0.3. The plane strain condition is assumed. The reference mode I SIF is given by
(20)KIexact=F(ab)σaπ=1.6118Pamm,
where a=0.5 is the crack length, b is the sheet width and F(a/b) is given by
(21)F(ab)=1.12-0.231(ab)+10.55(ab)2-21.72(ab)3+30.39(ab)4.Figure 8
Sheet with edge crack under tension.The strain energy and the error in the energy norm are defined as(22)E(Ω)=12∫ΩϵTDϵdΩ,ee=|E(Ω)num-E(Ω)refE(Ω)ref|1/2,ek=|Ksifnum-KsifrefKsifref|1/2×100%,sif=I,II,
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 KI 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 versush (mesh size) for the sheet with an edge crack under tension.Figure 11
The convergence in the energy norm of XFEM versush (mesh size) for the sheet with an edge crack under tension.Figure 12
The convergence in the stress intensity factorKI versus h (mesh size) for the sheet with edge crack under tension.Figure 13
The convergence in the stress intensity factorKI of XFEM versus h (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 ofKI 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 KI.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 modeI SIF KI 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 Figure16. The material parameters are Young’s modulus E=3×107 Pa and Poisson’s ratio ν=0.25. The exact stress intensity factors for this load case are given by [31]
(23)KI=34.0Pamm;KII=4.55Pamm.Figure 16
Sheet with edge crack under shear.The results shown in Figures16, 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 KI. 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 versush (mesh size) for a sheet with an edge crack under shear.Figure 19
The convergence in the energy norm of XFEM versush (mesh size) for a sheet with an edge crack under shear.Figure 20
The convergence in the stress intensity factorKI versus h (mesh size) for a sheet with an edge crack under shear.Figure 21
The convergence in the stress intensity factorKI of XFEM versus h (mesh size) for a sheet with an edge crack under shear.Figure 22
The convergence in the stress intensity factorKII versus h (mesh size) for sheet with edge crack under shear.Figure 23
The convergence in the stress intensity factorKII of XFEM versus h (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 Figure24 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.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 × 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 θ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)
(b)(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)Stress (a)σxx and (b) σyy contours in the sheet (structured mesh) after the crack propagates.
(a)
(b)Stress (a)σxx and (b) σyy contours in the sheet (unstructured mesh) after the crack propagates.
(a)
(b)


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



---

*Source: 978026-2013-07-25.xml*