An adaptive singular ES-FEM for mechanics problems with singular field of arbitrary order

Comput. Methods Appl. Mech. Engrg. 253 (2013) 252–273

Contents lists available at SciVerse ScienceDirect

Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma

An adaptive singular ES-FEM for mechanics problems with singular field of arbitrary order H. Nguyen-Xuan a,b,⇑, G.R. Liu c, S. Bordas d, S. Natarajan e, T. Rabczuk f a

Department of Mechanics, Faculty of Mathematics & Computer Science, University of Science, VNU-HCMC, Nguyen Van Cu Street, District 5, Ho Chi Minh City 700000, Viet Nam Division of Computational Mechanics, Ton Duc Thang University, Nguyen Huu Tho Street, District 7, Ho Chi Minh City 700000, Viet Nam c School of Aerospace Systems, University of Cincinnati, 2851 Woodside Dr, Cincinnati, OH 45221, USA d IMAM, Theoretical and Computational Mechanics, Cardiff University, UK e Department of Aerospace Engineering, Indian Institute of Science, Bangalore, India f Department of Civil Engineering, Bauhaus-Universität Weimar, Germany b

a r t i c l e

i n f o

Article history: Received 23 August 2011 Received in revised form 24 July 2012 Accepted 31 July 2012 Available online 14 August 2012 Keywords: Adaptive finite elements Smoothed finite element method Singular ES-FEM Singularity Crack propagation

a b s t r a c t This paper presents a singular edge-based smoothed finite element method (sES-FEM) for mechanics problems with singular stress fields of arbitrary order. The sES-FEM uses a basic mesh of three-noded linear triangular (T3) elements and a special layer of five-noded singular triangular elements (sT5) connected to the singular-point of the stress field. The sT5 element has an additional node on each of the two edges connected to the singular-point. It allows us to represent simple and efficient enrichment with desired terms for the displacement field near the singular-point with the satisfaction of partition-of-unity property. The stiffness matrix of the discretized system is then obtained using the assumed displacement values (not the derivatives) over smoothing domains associated with the edges of elements. An adaptive procedure for the sES-FEM is proposed to enhance the quality of the solution with minimized number of nodes. Several numerical examples are provided to validate the reliability of the present sES-FEM method. Ó 2012 Elsevier B.V. All rights reserved.

1. Introduction Smoothed finite element methods (S-FEMs) have been studied widely in the past few years [1]. The essential idea in the S-FEMs is to reconstruct the compatible strain field in finite element settings using a strain smoothing technique [2]. In the S-FEMs, the reconstructed strain field is obtained over various sub-domains called smoothing domains created on cells, nodes, edges or faces of the background mesh. The numerical operations used in S-FEMs can bring information from the neighboring elements depends on the requirements of the analyst, and hence are distinct in many ways from the standard FEM. The S-FEM can also be viewed as a combination of the treatments in both FEM and meshfree methods [3–5]. Since smoothing domains in S-FEM usually cover part of the adjacent elements, therefore the number of supporting nodes associated with the smoothing domain is larger than the set of nodes of the element. As a result, except a cell-based smoothing (CS-FEM), the bandwidth of the stiffness matrix in S-FEM is increased, and the computational cost is hence higher than the standard FEM with ⇑ Corresponding author at: Department of Mechanics, Faculty of Mathematics & Computer Science, University of Science, VNU-HCMC, Nguyen Van Cu Street, District 5, Ho Chi Minh City 700000, Viet Nam. E-mail address: [email protected] (H. Nguyen-Xuan). 0045-7825/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cma.2012.07.017

the same sets of nodes. On the other hand, thanks to the propagation of non-local information brought by the adjacent elements, the S-FEM often produces more accurate solutions than the standard FEM. For a given computational cost, all smoothed finite element methods have been shown to lead to a lower error than the displacement-based finite element method [6,7]. In addition, SFEM models can provide upper bound solutions, superconvergent stresses, and are volumetric free-locking [8] and shear locking [19–24], depending on types of the smoothing domains employed. Similar to the standard FEM, the S-FEM employs a finite element mesh. However, S-FEM evaluates the weak form based on smoothing domains created from the elements such as cells, or nodes, or edges or faces. As a result, a family of the S-FEMs was proposed as a cell-based smoothing (CS-FEM) [9], a node-based smoothing (NS-FEM) [10,11], an edge-based smoothing (ES-FEM) [12] and a face-based smoothing (FS-FEM) [13]. Applying the strain smoothing technique on smoothing domains helps to soften the over-stiffness of the standard FEM and hence can improve significantly the accuracy of both primal and dual quantities. Several theoretical aspects of the S-FEM were provided in [14–17]. The S-FEM models were also applied to a wide range of practical mechanics problems such as dynamics [12,18], plates and shells [19–24], piezoelectricity [25], fracture mechanics [26–34], visco-elastoplasticity [35–37], limit and shakedown analysis [38], fluid–structure

H. Nguyen-Xuan et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 252–273

interaction [39], etc. It can be seen that the S-FEMs have been becoming a simple and effective tool for analysis of various types of practical problems. Among these S-FEMs, the ES-FEM was found to be the most computationally efficient [1]. A singular five-noded triangular (sES-FEM-T3-5) (‘s’ for singular) element [30] that can produce a proper order of stress singularity near the crack tip was recently proposed based upon triangular simplex (T3) elements. This element requires only an additional node on each of the two edges of each element connected directly to the crack tip while linear T3 elements remain unchanged. Although the sES-FEM-T3-5 was found to perform well for various fracture problems [31–34], it is computationally expensive as uniform mesh was used in [31–34]. This is because no error control technique is available for the sES-FEM for adaptive mesh refinement near the crack tip. As already investigated in [33] for crack propagation problem, the mesh refinement at the crack tip is performed simply using the standard Delaunay triangulation procedure associated with the Laplacian smoothing technique without the use of any error indicator. It therefore is a purely geometric (crack location) based refinement. In addition, a study on other types of singular problems such as the singularity at the corner of L-shape re-entrants has not been conducted yet. In this paper, we present a formulation of the singular ES-FEM able to reproduce the stress singularities of arbitrary order using a base mesh of T3 elements with a layer of sT5 elements in the vicinity of the singularity. The stiffness matrix of the discretized system is computed using the assumed displacement values (not the derivatives) over smoothing domains associated with the edges of elements. As a result, the need to integrate singular terms r k1 ðk < 1Þ is eliminated. The present method will be denoted as sES-FEM-T3-5. An adaptive sES-FEM-T3-5 procedure is then proposed to enhance the computational efficiency of the original sES-FEM-T3-5 form. We show that the present approach performs well for various problems with weak and strong singularities. Also, the formulation is used to model crack growth. The effectiveness and reliability of the present method are confirmed by several numerical examples. The paper is organized as follows: In the next section, a brief on the problem at hand and the ES-FEM formulation is given. Section 3 introduces the general form of singular five-node formulation for singularity problems in linear elasticity. An adaptive procedure is presented Section 4. Solution procedure is summarized in Section 5. Several numerical examples are provided in Section 6. Section 7 closes some concluding remarks and opens future work.

Fig. 1. A two-dimensional domain with a reentrant corner resulting in a singular stress field of arbitrary order depending on the vertex angle.

 i on CD ; ui ¼ u

ð2Þ

where oi = o/oxi are the first partial derivatives corresponding to xi e {x, y} and o denotes a matrix of differential operators:

" @¼

@ @x

0

@ @y

0

@ @y

@ @x

#T :

2.1. Governing equations Consider a two dimensional (2D) linear-elastic solid defined in a domain X bounded by C that can contain reentrant corners such that C = CD [ CN, CD \ CN = £ as shown in Fig. 1. A distributed body force b acts within the domain. Boundary C is split into two parts, namely CD where displacements u are prescribed (Dirichlet conditions) and it is assumed to be composed of smooth open parts CDi , and CN where tractions t are prescribed (Neumann conditions). The relations between the displacement field u, the strain field e and the stress field r are: (1) The compatibility displacements

1 2

relation

between

the

8i; j ¼ 1; 2 : eij ¼ ð@ j ui þ @ i uj Þ or e ¼ @u in X;

strains

and

ð1Þ

ð3Þ

(2) The constitutive relations for solids is given as

rij ¼ Dijkl ekl in X or r ¼ De in X;

ð4Þ

where D is the Hooke matrix of elastic constants. (3) The equilibrium equations

@ j rij þ bi ¼ 0 in X rij nj ¼ ti on CN :

ð5Þ

It is well known that the displacements of such a singular problem in the vicinity of the singular point has the form of ui  rk , where k depends only on the vertex angles of p < / 6 2p and satisfies the real simple root of the following characteristic equation [40]:

sin kð1Þ u þ kð1Þ sin u ¼ 0; sin kð2Þ u  kð2Þ sin u ¼ 0:

2. Problem description and ES-FEM formulation

253

ð6Þ

For example, as / = 2p, we obtain the value k ¼ 1=2 resulting in the case of a fracture problem. Existing numerical methods for solving singular mechanics problems can be found in the literature such as the h-and p-versions of the finite element method [40], boundary element methods [41], meshfree methods [42–44], extended finite elements [45,46,26–28]. In these methods, the space of kinematically admissible displacements with continuous/discontinuous functions is commonly constructed such that the regularity of the solution is ensured. More detailed investigation of this space construction has already been well known in the literature [47]. Recently, a G space [17] that can accommodate certain types of continuous/discontinuous functions has also been introduced to construct simple yet effective shape functions for numerical methods. This is achieved by constructing the weakened weak (W2) formulation using, for example, various types of smoothing domains. In this work, we use the weakened weak (W2) form with edgebased smoothing domains and then present a singular edge-based smoothed finite element method (sES-FEM) for modeling singular mechanics problems.

254

H. Nguyen-Xuan et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 252–273

2

2.2. Brief on edge-based strain smoothing In the ES-FEM-T3 [12], the strains are ‘‘smoothed’’ over local smoothing domains, and naturally the computation for the stiffness matrix is no longer based on elements, but on these smoothing domains. These smoothing domains are Nconstructed based on Sed s the edges of the elements such that X  Xk and Xsi \ Xsj ¼ ; for i – j, in which Ned is the total number ofk¼1 edges of all elements in the entire domain. For triangular elements, the smoothing domain Xsk associated with edge k is created by connecting the two end-nodes of the edge to centroids of the adjacent elements as depicted in Fig. 2. Using the strain smoothing operator over the edge-based smoothing domain Xsk , we write

ek ¼

Z

Xsk

eh ðxÞUðxÞdX;

ð7Þ

where U(x) is a distribution function or a smoothing function that satisfies the following properties [2]:

UðxÞ P 0 and

Z Xsk

UðxÞdX ¼ 1:

ð8Þ

For simplicity, the smoothing function U is assumed to be a piecewise constant function and is given by

(

UðxÞ ¼

1=Ask 0

x 2 Xsk ; x R Xsk :

ð9Þ

Substituting Eq. (9) into Eq. (7) and applying the divergence theorem, the smoothed strains become

ek ¼ ¼

1 Ask 1 Ask

Z Xsk

Z

Csk

h

e

1 ðxÞdX ¼ s Ak

Z Xsk

nðkÞ ðxÞuh ðxÞdC;

h

ðkÞ

3

0

nx

6 nðkÞ ðxÞ ¼ 6 40

ðkÞ 7 ny 7 5 ðkÞ nx

ðkÞ

ny

ð11Þ

and Ask is the area of the smoothing cell Xsk computed by

Ask ¼

Z Xsk

Nk

dX ¼

e 1X As ; 3 i¼1 i

ð12Þ

where Nke is the number of elements containing edge k (Nke ¼ 1 for the boundary edges and N ke ¼ 2 for inner edges as shown in Fig. 2) and Ai is the area of the ith element around the edge k. 2.3. Weakened weak statement for ES-FEM  h 2 Vh  V such The discrete solution of the problem is to find u that the following ‘‘smoothed’’ Galerkin weak form is satisfied

8v h 2 Vh0  V0 ;

 h ; v h Þ ¼ f ðv h Þ;  ðu a

ð13Þ

(.,.) is a ‘‘smoothed’’ bilinear form of smoothed derivatives where a of functions and has the following form:

 h ; vhÞ ¼ ðu a

Ns Z X k¼1

Xsk

 h ÞDek ðv h ÞdX ¼ eTk ðu

Ns X  h ÞDek ðv h Þ; Ak eTk ðu

ð14Þ

k¼1

where V and V0 are the two spaces of kinematically admissible displacements u, respectively, defined by

V ¼ fu 2 ðH1 ðXÞÞ2 ; u ¼ uC on CD g;

ð15Þ

V0 ¼ fv 2 V; v ¼ 0 on CD g

ð16Þ

1

@u ðxÞdX ð10Þ

where Csk is the boundary of the smoothing domain Xsk , and n(k) is the outward normal matrix on the boundary Csk defined by

and H ðXÞ is a Hilbert space as defined in [17].It was proved that Eq. (13) has a unique solution for stable solids. Eq. (13) is known as a weakened weak (W2) statement, because in evaluating the smoothed bilinear form, we use only the assumed displacements themselves, and no derivatives of the assumed displacements are needed as in the usual weak form [17].

3. A singular ES-FEM formulation of with singularities of arbitrary order 3.1. Interpolation of displacements along the element edge

Fig. 2. Division of a cracked solid domain into three-noded triangular (T3) elements and smoothing cells Xsk connected to edge k of based on the T3 mesh.

One of the earliest studies on singularities focused on a fracture problem, where the stress field is singular near the crack tip and hence the polynomial basis functions used in the standard finite elements (FE) cannot reproduce directly the stress singularity near the crack tip. A simple way to reproduce such a singular field in the standard FE formulation is to construct the so-called six/eightnode quarter-point singular elements, in which the mid-side nodes are moved to the quarter points toward the crack tip [48,49]. The approximate displacement field can produce a r1/2 stress (derivative of displacement) singularity through an isoparametric mapping procedure. Such singular elements have been applied widely in practice and are implemented in most FEM packages. Recently, an alternative approach called the singular S-FEM [30] has been proposed, in which the displacement filed is also constructed based on the required features directly as a function of the coordinates. The formulation uses a simple point interpolation method (PIM). There is no mapping involved and only the shape function values (not the derivatives) are used to compute the stiffness matrix. Fig. 2 illustrates a problem domain with a horizontal opening crack, where nodes are added only on the edges of triangular elements connected to the crack tip. In the immediate layer of elements around the crack tip, each element has five nodes as

H. Nguyen-Xuan et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 252–273

255

Fig. 3. An illustration of 5-node triangular elements with an additional node on each of the triangular elements connected to the singular point: (a) a L-shape domain (/ = 3p/ 2); (b) a crack domain (/ = 2p); and (c) local coordinate originated at the singular node and the position of an additional node.

Fig. 4. Two-layer smoothing domains near the singular point: (a) edge-based smoothing domain in the layer of crack tip elements and (b) subdivision of smoothing domains into 2 sub-smoothing domains Xs;a k ða ¼ 1; 2Þ.

256

H. Nguyen-Xuan et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 252–273

Fig. 5. A subdivision of a mother triangle into two new children triangles: (a) the midpoint of the refinement edge is generated and (b) two sub-triangles share a new common peak.

(b) η = 23.15% (η : relative error percentage)

(a)

Fig. 6. A L-shape domain and a uniform mesh of 65 nodes. −4

2.45

x 10

Exact FEM−T3 NS−FEM−T3 ES−FEM−T3 sES−FEM−T3−5

2.4

2.35

−2.2

FEM−T3(r = 0.59) NS−FEM−T3 (r = 51) ES−FEM−T3 (r = 0.58) sES−FEM−T3−5 (r = 0.61)

−2.3 −2.4 −2.5

2.25

−2.6

log10(e )

2.2

e

Strain energy

2.3

2.15

−2.7 −2.8 −2.9

2.1

−3

2.05 −3.1

2

0

500

1000

1500

2000

2500

3000

Number of nodes

−3.2 −1.3

−1.2

−1.1

−1

−0.9

−0.8

−0.7

−0.6

−0.5

−0.4

−0.3

log10(h) Fig. 7. Convergence of strain energy for a L-shape domain using ‘‘uniform’’ T3 meshes.

Fig. 8. Convergence rate of a L-shape domain using uniform T3 meshes.

257

H. Nguyen-Xuan et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 252–273

shown in Fig. 3b. It is worth noting that the location of additional intermediate nodes can be in general at any position on the edge, and not necessarily at the quarter positions. Consider the component uh of the displacement field  2 Vh  V at an arbitrary point of interest on an edge of an eleu ment connected to the singular point assumed as

where d ¼ f u1 u2 u3 gT is the vector of nodal displacements and C is a 3  3 matrix depending on the nodal r-coordinates:

uh ðxÞ ¼ PT ðxÞa ¼ a1 þ a2 r þ a3 r k ;

Solving Eqs. (18) for a and substituting them back to Eq. (17), one obtains

ð17Þ

where 0 6 r 6 l is the radial coordinate originated at the singular node, 0 < v < 1 is used to tune the location of additional node 2 as shown in Fig. 3, ai (i = 0, 1, 2, 3) are the coefficients to be determined by interpolation of u at nodes, and 1=2  k < 1 is a singularitydependent parameter chosen to reproduce a desired displacement field around the singular point. Substituting the radial coordinates at nodes into Eq. (17), it yields

u1 ¼ a1 k

u2 ¼ a1 þ vla2 þ ðvlÞ a3

or d ¼ Ca;

ð18Þ

2

3 1 0 0 6 1 vl ðvlÞk 7 C¼4 5: 1 l

l

ð19Þ

k

uh ðxÞ ¼ PT ðxÞC1 d ¼ /ðxÞ d;

ð20Þ

where / is a matrix of shape functions defined by

/ ¼ ½ /1

/2

/3 :

ð21Þ

For example, when a crack is considered, the shape functions shall be with v ¼ 14 and k ¼ 12:

rffiffiffi r r /1 ¼ 1 þ 2  3 ; l l

rffiffiffi r r /2 ¼ 4 þ 4 ; l l

r /3 ¼ 2  l

rffiffiffi r : l

ð22Þ

k

u3 ¼ a1 þ la2 þ l a3

(a) η = 12.5%

(b) η = 7.16%

(c) η = 4.69%

(d) η = 3.0%

Fig. 9. An illustration of adaptive meshes automatically generated based on the error indicators for the L-shape plate: (a) 119 nodes; (b) 255 nodes; (c) 708 nodes; (d) 1369 nodes.

258

H. Nguyen-Xuan et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 252–273

In this case, the present sES-FEM-T3-5 formulation is identical to the singular ES-FEM presented in [30]. It is seen that /i satisfy all the basic properties of the shape functions: partition of unity, linear reproducibility and Kronecker delta function properties [8]. Most importantly, these shape functions are enriched with r k that can produce a stress field with singularity of order of r k1 in the vicinity of the singular point. Note that only the shape function values (not the derivatives) (see e.g. Eq. (10)) will be used to compute the stiffness matrix. Hence, integrations with the rk1 terms are completely avoided, but we will still need more integration points to integrate r k , 1=2  k < 1, as this is not a polynomial. 3.2. Displacement interpolation within a five-node singular element In the present formulation, only the elements connected to the singular point are with additional nodes on their edges. Consider a 5-node element 1–4–2–3–5 as shown in Fig. 4. We assume now that the displacement field for any point of interest forms the enriched form of Eq. (20). We now add two points (not nodes) 6 and 7 to the midpoints of lines 4–5 and 2–3, respectively. The displacement at these points can be assumed to be computed using the linear PIM:

u6 ¼

1 ðu4 þ u5 Þ; 2

u7 ¼

1 ðu2 þ u3 Þ: 2

ð23Þ

This assumption means that the displacement at a point #6 can be ‘‘linearly’’ evaluated (using only the linear PIM) from information of two additional nodes on two element edges connected to the singular point. Note that the linear interpolation u6 depends on u4 and u5 (or depends on the displacement field uh(x) (see in Eq. (20)) evaluated at u4 and u5) on two edges of the element connected to the singular point, respectively. Hence, Eq. (23) is valid for any linear interpolation of the present 5-node element at any point. Therefore, it is not required any justification in Eq. (23) for the 5-node element, even close to the singular point. In the present sES-FEM-T3-5, multi-layer smoothing domains are used in order to better capture the singularity. Each 5-node

element is comprised of five edge-based smoothing domains: four are the singular-point smoothing domains and one is a standard smoothing domain (see Fig. 4). It is worth noting that in the formulation of the sES-FEM-T3-5, only the shape function values at the Gauss points along the boundary segments are required to calculate the smoothed strain gradient matrix. The interpolation schemes for the boundary segments are dependant (slave) on what is constructed in Section 3.1 for element edges. In fact the strain displacement matrix relies on simply computing using the linear PIM. In the ES-FEM-T3, only normal smoothing domains are used, and the shape function is linearly compatible along boundary segments. Therefore, only one Gauss point is needed on each boundary segment. In the sES-FEM-T3-5, however, approximate number of Gauss points are required on the boundary segments of the singular-point smoothing domains, due to the complexity of the assumed displacement field (or shape function) along these boundary segments. Because the displacement in the vicinity of the singular point possesses a r k behavior, more Gauss points are required. As shown in Fig. 4, for the singular-point smoothing domains with boundary segments of A12–2B1–B1B2–B24–4A2–A2A1, there are two kinds of segments: (1) one is along the quasi-tangential direction (A12, 2B1, B24, and 4A2); (2) another is along the radial direction (A2A1 and B1B2). We now construct the shape function for Gauss points on two kinds of boundary segments. – For the Gauss point Gx on a segment along the quasi-tangential direction (e.g. 2B1, this point is also on the line 1–C–M as shown in Fig. 4), the displacement is interpolated as:

uh ¼ /1 u1 þ /2 uc þ /3 uM

ð24Þ

in which

uc ¼

  lc4 lc4 u4 þ 1 u5 ; l45 l45

 uM ¼

1

 lM2 lM2 u2 þ u3 l23 l23 ð25Þ

and li–j is the distance between two points i and j. It is proved that lc4 ¼ llM2 ¼ a and hence the displacement interpolation within a l45 23 crack can be rewritten as

−2.4 −2.5 −2.6

25

sES−FEM−T3−5 Adaptive sES−FEM−T3−5

−2.7

20

e

log10(e )

−2.8

−3 −3.1 −3.2

sES−FEM−T3−5 (r = 0.61) Adaptive sES−FEM−T3−5 (r = 1.08)

−3.3 −3.4 −1.3

−1.2

−1.1

−1

−0.9

−0.8 −0.7 log10(h)

−0.6

−0.5

−0.4

−0.3

Fig. 10. The convergence rate of L-shape domain using uniform and adaptive meshes. It is observed that the optimal convergence rate of 1.0 is achieved with adaptive refinement for linear weakform models: superconvergence. Note that the characteristic length h canpbe regarded as the average length of the sides of the ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi elements defined by h ¼ 2AX =N e where AX is the area of the whole problem domain and Ne is the total number of triangular elements in the whole problem domain.

Relative error (%)

−2.9

15

10

5

0

0

500

1000

1500

2000

2500

3000

3500

Number of nodes Fig. 11. Relative error percentage (g) in the energy norm of L-shape domain using uniform and adaptive meshes.

259

H. Nguyen-Xuan et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 252–273

uh ¼ /1 u1 þ ð1  aÞ/3 u2 þ a/3 u3 þ ð1  aÞ/2 u4 þ a/2 u5 : |{z} |{z} |{z} |fflfflfflfflfflffl{zfflfflfflfflfflffl} |fflfflfflfflfflffl{zfflfflfflfflfflffl} N1

N2

N3

N4

N5

ð26Þ – For the Gauss point Gy on a segment along the radial direction (e.g. B1–B2), the displacement is interpolated as:

uh ¼ /1 u1 þ /2 u6 þ /3 u7 1 1 1 1 ¼ /1 u1 þ /3 u2 þ /3 u3 þ /2 u4 þ /2 u5 : 2 2 2 2

ð27Þ

It is seen that the displacement formulation given in (27) is a specific case of (26) with a = 0.5. The displacement interpolation is simply obtained by setting (26) with a = 0 and a = 1 for the displacement at any point along edges 1-4-2 and 1-5-3, respectively. Thus, it can be said that Eq. (26) can be considered as the general form of the displacement interpolation for the five-node singular-point triangular element. In general, we can construct an n-noded singular element following the similar displacement interpolation procedure as the sES-FEM-T3-5 element. A 7-node crack-tip element for fracture mechanics problems has already been formulated in [50]. However, the mesh of this element for crack propagation simulation is refined only around the crack tip following the standard Delaunay triangulation procedure, which is a fixed level refinement without the use of any error indicator and the refinement is based only the geometry and this process is computationally expensive. 3.3. Smoothed stiffness matrix of the present element The stiffness matrix of the present element can be evaluated based on two types of smoothing domains: (1) smoothing domains which are not directly connected to the singular point (scheme #1) and (2) smoothing domains directly connected to the singular point (scheme #2). For scheme #1, the computation is identical to the ES-FEM-T3 [12]. The smoothed strain in (10) for the smoothing domain Xsk is obtained using the following matrix form

(a)

ek ¼

X

I ;  I ðxk Þd B

ð28Þ

ðkÞ I2Wn

 I is the unknown displacement at node I; WðkÞ is the set of where d n  I ðxk Þ is the nodes supporting to the smoothing domain Xsk and B smoothed strain-displacement matrix on the smoothing domain Xsk and is calculated numerically by

2 3 bIx ðxk Þ 0 Iy ðx Þ 7  I ðxk Þ ¼ 6 B b 4 0 k 5;   bIy ðxk Þ bIx ðxk Þ

ð29Þ

where

 ðx Þ ¼ 1 b Ih k Ask

Z Csk

ðkÞ

NI ðxÞnh ðxÞdC;

ðh ¼ x; yÞ:

ð30Þ

From the above equation, it is observed that we only need the values of NI (not the derivatives) on the boundary of the smoothing domain Csk , which can be also obtained by a simple linear interpolation technique [8]. In such cases, one Gauss point is sufficient for line integration along each segment of the boundary Cskb 2 Csk . Eq. (30) can be further simplified to its algebraic form n

eg X   ðkÞ ðkÞ  ðx Þ ¼ 1 b NI xGP Ih k s b nbh lb ; Ak b¼1

ðh ¼ x; yÞ;

ð31Þ

where neg is the total number of the segments Cskb 2 Csk ; xGP b is the ðkÞ midpoint (Gauss point) of the boundary segment of Cskb ; lb and ðkÞ s nbh are the length and the outward unit normal of Ckb , respectively. Scheme #2 is used for the sES-FEM-T3-5 element to improve the performance of the ES-FEM for elements containing a singularity. Based on the basic T3 mesh, only a small set of nodes is added to fit the 5-node element. Hence, the total number of nodes of the entire domain is significantly reduced when compared with the traditional FEM-T6, while additional nodes around the singular point only rely on the background T3 mesh. After having placed additional nodes, the displacement field is approximated by Eq. (26). To ensure the stability of the present element, each smoothing domain containing a singular point should be further sub-divided into several sub-smoothing domains similarly to what is

(b) η = 16.0% (η : relative error percentage)

Fig. 12. (a) Plate with an edge crack under a uniform tension; and (b) a uniform triangular mesh of 234 nodes and the relative error percentage.

260

H. Nguyen-Xuan et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 252–273

done in the extended finite element method [51] and in its smoothed version [26–28]. The number of subcells required depends on the singular element being formulated. In the 5-node singular element, only two sub-smoothing domains can be sufficient to ensure both the stability and accuracy of the element. The smoothed strain-displacement matrix on each sub-smoothing domain Xs;a k (a = 1, 2) can be expressed as

2

6  a ðxk Þ ¼ 6 B I 4

 a ðx Þ b k Ix

3

0

a ðx Þ 7 7 b k 5 Iy a a   bIy ðxk Þ bIx ðxk Þ 0

ð32Þ

in which

 a ðx Þ ¼ 1 b k Ih As;a k

Z Cs;a k

ðk;aÞ

NI ðxÞnh

ðxÞdC;

ðh ¼ x; yÞ;

ð33Þ

where Aks;a , Cs;a k is the area and the boundary of the sub-smoothing s cell Xs;a and C k associated with the edge k, respectively. k It is noted from (33) that the shape functions NI in (26) are now employed. Again, applying Gauss integration along the segments of boundary Cs;a k , one obtains

" # neg nX Gauss 1 X ðk;aÞ a GP GP  bIh ¼ s;a wb;c NI ðxb;c Þnbh ðxb;c Þ ðh ¼ x; yÞ; Ak b¼1 c¼1

ð34Þ

where nGauss is the number of Gauss points used in each segment, ðk;aÞ wb,c is the corresponding weight of Gauss points, nbh is the outward unit normal corresponding to each segment of the boundary s;a th th GP Cs;a kb 2 Ck and xb;c is the c Gauss point of the b boundary segment of Cs;a . The numerical integration in (34) requires more Gauss points k than the linear case given in (31). A detailed study of using approximate number of Gauss points for line integration along the boundary of the singular-point smoothing domains has already been addressed in [34]. In this work, three Gauss points on each boundary segment of sub-smoothing domains (nGauss = 3) are sufficient to provide reasonable accuracy. By substituting Eqs. (26) and (28) into Eq. (13), we obtain a linear algebraic system of equations

 ¼ f; d K

ð35Þ

where f is the nodal force vector that can be obtained by:

f¼

Z

NT bdX þ

Z

NT tC dC:

ð36Þ

Ct

X

−3

x 10

1.22

(a) 1.2

Strain energy

1.18

1.16

1.14

Ref FEM−T3 NS−FEM−T3 XFEM−Q4 ES−FEM−T3 sES−FEM−T3−5 sFEM−T6

1.12

1.1

(a) η = 11.1%

(b) η = 7.8%

0

1000

2000

3000

4000

5000

6000

7000

8000

Number of nodes 1.164

x 10

−3

(b) 1.162

Strain energy

1.16

1.158

1.156

Ref sES−FEM−T3−5 Adaptive sES−FEM−T3−5

1.154

(c) η = 5.9%

(d)η = 4.4%

Fig. 13. An illustration of adaptive meshes for a plate with an edge crack under remote tension and the relative error percentage (g): (a) 267 nodes; (b) 373 nodes; (c) 585 nodes; (d) 1038 nodes.

0

1000

2000

3000 4000 5000 Number of nodes

6000

7000

8000

Fig. 14. Convergence of the strain energy of a plate with an edge crack under remote tension: (a) uniform meshes and (b) adaptive meshes.

261

H. Nguyen-Xuan et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 252–273

 of the system is then assembled by The smoothed stiffness matrix K a similar process as in the FEM

1.7

(a)

1.68 Ned X ¼  ðkÞ ; K K

ð37Þ

1.66

k¼1

1.64

 ðkÞ is the stiffness matrix associated with smoothing domain where K Xsk and is calculated by

I

J

k

1.62

SIF KI

 ðkÞ K IJ

8 T  s s > < BI DBJ Ak for Xk without containing the singular point 3 ¼ X  aT DB  a As;a for Xs containing the singular point: > B :

1.6 1.58

k

a¼1

Ref FEM−T3 NS−FEM−T3 XFEM−Q4 ES−FEM−T3 sES−FEM−T3−5 sFEM−T6

1.56

ð38Þ

1.54

4. An adaptive procedure

1.52 1.5 0

Adaptive techniques are important for problems with singular stress fields, because of the drastic variation in displacements and stresses over the problem domain especially near the singular points. Naturally one needs a mesh that should be finer near the singular points and gradually coarser away from the singular

1000

2000

3000

4000 5000 Number of nodes

6000

7000

8000

1.612

(b)

−2

1.61

(a) −2.1

1.608

SIF KI

−2.2 −2.3

1.604

−2.4

log10(ee)

1.606

−2.5

1.602

−2.6

Ref sES−FEM−T3−5 Adaptive sES−FEM−T3−5

1.6 −2.7

FEM−T3 (r = 0.44) NS−FEM−T3 (r = 0.46) XFEM−Q4 (r = 0.49) ES−FEM−T3 ( r =0.49) sES−FEM−T3−5 (r = 0.57)

−2.8 −2.9 −3 −1.8

−1.7

−1.6

−1.5

−1.4

−1.3

−1.2

−1.1

−1

1.598 0

1000

2000

3000

4000 5000 Number of nodes

6000

7000

8000

Fig. 16. Convergence of stress intensity factor KI of a plate with an edge crack under remote tension: (a) uniform meshes and (b) adaptive meshes.

log10(h) −2.4

−2.5

sES−FEM−T3−5 (r = 0.57) Adaptive sES−FEM−T3−5 (r = 1.21)

(b)

−2.6

log10(ee)

−2.7

−2.8

−2.9

−3

−3.1

−3.2 −1.8

−1.7

−1.6

−1.5

−1.4

−1.3

−1.2

−1.1

−1

log10(h)

Fig. 15. Convergence rate in energy norm of a plate with an edge crack under remote tension: (a) uniform meshes; (b) adaptive meshes. It is observed that we have exceeded the optimal convergence rate of 1.0 for linear weakform models: superconvergence in the adaptive case.

(a)

(b) η = 21.3%

Fig. 17. Plate with an edge crack under shear: (a) a problem setup with dimensions in millimeters and (b) a uniform mesh of 331 nodes.

262

H. Nguyen-Xuan et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 252–273

points. Ideally, this should be done automatically or adaptively, and has been attempted in the past using the FEM and the extended finite element method [55,56]. In such an automated adaptive analysis procedure, we need an effective error estimator/ indicator to assess the error caused the numerical methods. We propose here an error indicator that is suitable for the sES-FEMT3-5 for singular problems. The 2nd important ingredient is automatic local refinement. Because the sES-FEM-T3-5 uses mesh of triangles, such local refinements can be done easily and automatically, by essentially adapting the existing techniques used in the FEM with T3 elements. For singular problems, Ródenas et al. [57,58] have proposed remarkably accurate recovery techniques and applied them to the extended finite element method with singular solutions (linear elastic fracture mechanics).

the ‘‘exact’’ and the FEM solutions. The ‘‘exact’’ one is, however, not known in general. Therefore one can produce an approximated ‘‘exact’’ solution by smoothing the FE stresses across the whole mesh with various recovery techniques [59]. Let rh be the FE approximation of the exact solutions r for stresses. The error in energy norm can be written as

kek2 ¼

Z

ðr  rh ÞT D1 ðr  rh ÞdX:

ð39Þ

X

^ be the recovery (smoothed) stress field which is expected Let r more accurate than the primitive solution rh, and hence can be used in lieu of the exact field r. This leads to a so-called posterior error estimate. A posteriori error estimate for FEM solutions can be obtained based on Zienkiewicz and Zhu (Z2) approach [59]:

Z

4.1. Recovery stresses

^k2 ¼ ke

^  rh ÞT D1 ðr ^  rh ÞdX: ðr

In the standard FEM for elasticity, the error analysis is performed using a posteriori error estimation process. The error in the energy norm should be measured by the difference between

It is obvious that Eq. (40) is different from Eq. (39). However, as the mesh is refined, the value yielded by (40) will tend towards that given by (39). Hence Eq. (40) is often called an error indicator.

ð40Þ

X

−5

9.5

x 10

−2.4

(a)

(b)

FEM−T3 (r = 0.47) NS−FEM−T3 (r = 0.56) XFEM−Q4 (r = 0.56) ES−FEM−T3 (r = 0.48) sES−FEM−T3−5 (r = 0.51)

−2.5

9 −2.6

log10(ee)

Strain energy

−2.7

8.5

8

Ref FEM−T3 NS−FEM−T3 XFEM−Q4 ES−FEM−T3 sES−FEM−T3−5

7.5

7

0

1000

2000

3000

4000

5000

6000

7000

8000

Number of nodes

−2.8

−2.9

−3

−3.1

−3.2

−3.3 −1

−0.7

−0.6

−0.5

−0.4

−0.3

−0.2

6

Ref FEM−T3 NS−FEM−T3 XFEM−Q4 ES−FEM−T3 sES−FEM−T3−5

38

−0.8

log10(h)

42

40

−0.9

Exact FEM−T3 NS−FEM−T3 XFEM−Q4 ES−FEM−T3 sES−FEM−T3−5

5.8 5.6 5.4

36

SIF KII

SIF KI

5.2

34

5 4.8

32

4.6

30 4.4

28

26

(c)

4

0

1000

2000

3000

4000

5000

Number of nodes

6000

7000

(d)

4.2

8000

0

1000

2000

3000

4000

5000

6000

7000

8000

Number of nodes

Fig. 18. Comparison of different methods for a plate with an edge crack under shear: (a) the strain energy; (b) the convergence rate; (c) SIF KI; (d) SIF KII.

263

H. Nguyen-Xuan et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 252–273

In this work, our error indicator uses a continuous recovery ^ is constructed by combining the smoothed nodal stress field. r stresses derived from the sES-FEM-T3-5 and the linear shape functions of the original triangular elements. Such an interpolation guarantees the stress field to be continuous on the whole domain. The recovery stresses can be defined as

r^ Rrh ¼

3 X ^I; NI r

ð41Þ

I¼1

where NI is the linear shape function at node I of the ^ I are the nodal stresses derived from the triangular element, r sES-FEM-T3-5 model via the simple nodal averaging technique. The nodal averaging technique is also known as first-order recovery stress [60] which has also been implemented in several commercial FEM codes. The difference in the sES-FEM-T3-5 model ^ I is computed using the stresses computed for is that the stress r each of the edge-based smoothing domains, rather than the elements. Note that the indicator used here is the simplest possible, as the goal of the paper is not to study in depth error estimation for singular problems. This was done recently in [55–58].

5. Solution procedure In this section, we describe the numerical implementation of the present method for solids with singular stress fields of arbitrary order. The algorithm is summarized as follows: 1. Define the problem domain containing the singular point (Step #1). 2. Discretize the problem domain with coarse mesh and: – add one node into each edge of triangular elements directly connected to the singular point; – create element connectivity including three-node triangular elements (FEM-T3) and a special layer of singular five-node triangular elements (FEM-T3-5), see e.g. Fig. 2 (Step #2). 3. Create edge-based smoothing domains as given in Section 3.2 (Step #3). 4. Loop over edge-based smoothing domains (Step #4) {

4.2. Error indicator and refinement strategy Once an error indicator is available, it is used for guiding the mesh refinements, leading to an automatic adaptive scheme. In this paper, we use Eq. (40) together with (46). Assume that the bounded domain X is discretized into a set L of nel elements such P ^ that X  Xh ¼ nel e¼1 Xe . The global error indicator, e, can be rewritten as the sum of the local refinement indicator for all the individual elements as

^ k2 ¼ ke

nel X ^e k2 ; ke

ð42Þ

e¼1

^e k is the estimated error norm of the element: where ke

^e k2 ¼ ke

Z

Xe

^  rh ÞT Dðr ^  rh ÞdX: ðr

ð43Þ

^e , we mark the elements Through the local refinement indicator, e Xe 2 L for refinement using the Dorfler criterion [61]. A minimal set M # L is determined such that

X

2

2

^k P #ke ^k ; for some # 2 ð0; 1Þ: ke

(a) η = 15.7%

(b) η = 12.3%

(c) η = 10.2%

(d) η = 6.8%

ð44Þ

X;in M

Then, a new mesh L0 is generated from L by refining (at least) the marked elements Xe 2 M. Now we briefly describe the algorithm of newest vertex bisection [62,63]. For each (mother) triangle Xe 2 L of the initial triangular mesh of X, we find the longest edge called base or refinement edge and the opposite vertex of the base is labeled peak or newest vertex, cf. Fig. 5. This process is known as a labeling of L. The rule of newest vertex bisection is drawn as follows: (1) Two new children triangles are generated from a mother triangle by connecting the peak to the midpoint of the base; (2) A midpoint (new node) of a base will be designated to create a peak of the children triangles. Once an initial triangulation is labeled, the proper triangulations inherit the label by the rule (2) such that the bisection process can continue. This refinement scheme using the newest vertex bisection is fast and is very simple to implement because the division is always on the longest edge, and mesh distortion is very well controlled.

Fig. 19. An illustration of adaptive meshes for a plate with an edge crack under shear and the relative error percentage: (a) 375 nodes; (b) 466 nodes; (c) 607 nodes; (d) 1288 nodes.

264

H. Nguyen-Xuan et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 252–273

{// call adaptive procedures in Section 5

Refinement

Recovery } // for crack problem If ProblemType = ‘‘crackType’’ – Compute the J –integral [51] and hence the stress intensity factors of KI and KII based on the smoothing domains [1]; – Calculate the crack propagation angle at the current tip; - Compute the equivalent stress intensity factor; if KIeq < Kc // Kc is the fracture toughness – Given an incremental crack Da; – update the crack path; – update the geometry of problem; – call Step #2; – call Step #4 else break out Loops exited early. end end }// End the loop over final propagation step

if Scheme #1 (or smoothing domains are not connected to singular point) (see Section 3.3)  I ðxk Þ in (29); – Evaluate the smoothed gradient matrix B  using Eq.(38)1; else // – Compute the stiffness matrix K Scheme #2 (or smoothing domains are not connected to the singular point)  I ðxk Þ in (32); – Evaluate the smoothed gradient matrix B  using Eq. (38)2; end // End – Compute the stiffness matrix K Schemes # 1&2 – Assemble the stiffness matrix and the load vector of the current smoothing domain into the global stiffness matrix and load vector; } // End the loop over the smoothing domains – Apply boundary conditions; – Solve system equations to determinate nodal displacements, strains and stresses; 5. Loop over each propagation step // step = 1 for all week and strong singularity static cases { – Compute recovery stresses based on Eq. (41) – Compute local and global error and global relative error (g) While global relative error (g) > given target error (gt) −5

8.46

x 10

−2.9

(a)

(b)

−3

8.42

−3.05

8.4

e

log10(e )

Strain energy

sES−FEM−T3−5 (r = 0.51) Adaptive sES−FEM−T3−5 (r = 0.98)

−2.95

8.44

8.38

−3.1

−3.15

8.36 −3.2

Ref sES−FEM−T3−5 Adaptive sES−FEM−T3−5

8.34

0

1000

2000

3000

4000

5000

6000

7000

−3.25

−3.3 −1

8000

−0.9

−0.8

−0.7

Number of nodes

−0.6

−0.5

−0.4

−0.3

−0.2

log10(h) 4.56

(d)

34

(c)

4.55

33.9 4.54

33.8

II

33.7

SIF K

SIF K

I

4.53

4.52

33.6 4.51

33.5

4.5

Exact sES−FEM−T3−5 Adaptive sES−FEM−T3−5

33.4

33.3

0

1000

2000

3000

4000

5000

Number of nodes

6000

7000

8000

Ref sES−FEM−T3−5 Adaptive sES−FEM−T3−5

4.49

0

1000

2000

3000

4000

5000

6000

7000

8000

Number of nodes

Fig. 20. Comparison of numerical results for a plate with an edge crack under shear using uniform and adaptive meshes: (a) strain energy; (b) convergence rate; (c) SIF KI; (d) SIF KII.

H. Nguyen-Xuan et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 252–273

(a) Initial coarse mesh of 98 DOFs (step 0)

(b) step 1: 1154 DOFs

(c) step 4: 1280 DOFs

(d) step 5: 1542 DOFs

(e) step 7: 1456 DOFs

(f) step 10: 1838 DOFs

(g) A zoom of deformation Fig. 21. Several mesh-refinements of a plate with an edge crack under shear after 10 steps, crack propagation trajectory and an illustration of deformation.

265

266

H. Nguyen-Xuan et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 252–273

6. Numerical examples

900

XFEM−Q4 Adaptive sES−FEM−T3−5

800

(a) In this section, some benchmark problems are studied. For comparison, the elements used in this paper are denoted as follows:

Equivalent SIF K

eq

700

XFEM-Q4 – the extended finite element method using standard 4-node quadrilateral elements that are locally enriched by Heaviside function and asymptotic functions to locate the crack [26].

FEM-T3 – the standard FEM using three-node element with shape linear function (T3).

NS-FEM-T3 – the node-based SFEM [10] using T3 elements.

ES-FEM-T3– the standard edge-based SFEM [12] using T3 elements with no enrichment.

sES-FEM-T3-5– the present singular five-node ES-FEM with enrichment of the singular point.

600 500 400 300 200 100 0 3.5

4

4.5

5

5.5

6

6.5

7

Crack length

This example aims to demonstrate that the present sES-FEMT3-5 can perform well for some week singularity problems using both uniform and adaptive meshes. Consider a L-shape plate with singularity at a re-entrant corner. The geometry is shown in Fig. 6a. The material parameters are Young’s modulus E = 105 Pa and Poisson’s ratio m = 0.3. Plane strain conditions are assumed. The analytical solution is given as [60]:

16

XFEM−Q4 (90986 DOFs) Adaptive sES−FEM−T3−5

14

6.1. L-shape with a re-entrant corner singularity

(b) 12

10

1 a r ½ða þ1Þcosðða þ1ÞhÞþðC 2 ða þ1ÞÞC 1 cosðða 1ÞhÞ; 2l 1 uh ðr;hÞ ¼ r a ½ða þ1Þsinðða þ1ÞhÞþðC 2 þ a 1ÞC 1 sinðða 1ÞhÞ; 2l ð45Þ

H=16

ur ðr;hÞ ¼ 8

6

4

where a  0.544483737 is the critical exponent that satisfies the equation of a sin 2x + sin 2xa = 0 for x = 3p/4, C1 = cos((a + 1)x)/cos((a  1)x), C 2 ¼ 2ðk þ 2lÞ=ðk þ lÞ and the Lamé constants k ¼ Em=ðð1 þ mÞð1  2mÞÞ; l ¼ E=ð2ð1 þ mÞÞ. The above solution is derived from the Navier-Lamé equation with zero prescribed forces and only Dirichlet conditions. This is considered

2

0

0

1

2

3

4

5

6

7

W=7 8

(c) 7.9

7.8

H/2

7.7

7.6

7.5

7.4

XFEM−Q4 Adaptive sES−FEM−T3−5 7.3 0

1

2

3

4

5

6

7

W/2 Fig. 22. Comparison of numerical results for a plate with an edge crack under shear load: (a) strain energy; (a) equivalent SIF; (b) crack growth trajectories; (c) a zoom of crack growth path.

Fig. 23. The PMMA beam with three holes subjected to a concentrated load (dimensions in inches).

H. Nguyen-Xuan et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 252–273

267

(a) 642 DOFs (initial step 0)

(b) 2354 DOFs (step 5)

(c) 2470 DOFs (step 9)

A zoom of deformation at step 7 Fig. 24. Crack growth trajectory for the PMMA beam (case I) with the crack increment at 0.5.

as a ‘‘displacement-driven’’ problem where external forces are zero and nonzero displacements are enforced on the boundary. In finite element analysis, the solution in equation (45) is used on the whole boundary of the problem. As a result, the energy bound of the dis-

placement-driven problem has the reverse sign the force-driven problem [64]. Fig. 6b illustrates an initial uniform mesh of 65 nodes. The convergence curves of strain energy are plotted in Fig. 7. It is seen that

268

H. Nguyen-Xuan et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 252–273

the FEM-T3, ES-FEM-T3 and sES-FEM-T3-5 produce an upper bound while the NS-FEM-T3 provides a lower bound of the exact strain energy. The solutions of the sES-FEM-T3-5 are more accurate than those of the ES-FEM-T3. The convergence of the numerical solution is proved by the theoretical analysis in [64]. Note that the NS-FEM-T3 possesses the essential properties of an equilibrium FEM model. More intuitional explanation related to the features of the equilibrium approach of the NS-FEM-T3 was given in [1]. Fig. 8 shows the convergence rate in the energy norm using a sequence of quasi-uniform meshes of various densities. It is seen that all methods can produce a convergence rate that is close to the theoretically expected rate a  0.544483737. In terms of accuracy levels, the sES-FEM-T3-5 is superior to the FEM-T3 and ESFEM-T3. However, it seems that the NS-FEM-T3 is slightly more accurate than the sES-FEM-T3-5 for this problem. This may be due to the super-accurate stress solution inherited from the equilibrium FEM model. To enhance the computational effect of the present sES-FEMT3-5, we analyze the problem with a sequence of adaptive meshes. The reliability of the present element compared with other elements has already been verified for uniform meshes. Therefore, we only focus on the accuracy of the sES-FEM-T3-5 itself with both uniform and adaptive meshes. For adaptive analysis, the local refinement indicator is derived from Eq. (43) and the criterion for refinement is provided in Eq. (44). For the problems tested below, the parameter # can be fixed at 0.25. Also, the algorithm of newest vertex bisection [62,63] is adopted to subdivide triangles refinement. The refinement process is started from a uniformly coarse mesh as given in Fig. 6b. Fig. 9 illustrates four adaptive meshes with respect to the relative error percentage (g) obtained. It is observed that the mesh refinement happens primarily near the re-entrant corner singularity as expected. Fig. 10 shows the convergence rate of the energy norm using adaptive meshes in comparison with uniform meshes. It is clear that using adaptive mesh-refinement results in the expected convergence order to 1.0 and improves very significant the computational cost. Also as seen from Fig. 11, the relative error percentage (g) in the energy norm shows again very good performance of the proposed method used in tandem with the adaptive algorithm. Now, in the following sections, we will prove the high performance of the present formulation for analyzing static crack and crack propagation. In order to illustrate the method, we restrict our regards to linear elastic fracture mechanics and material regime around the crack tip remains elastic. An extension of the present formulation to cohesive fracture model will be of our interest in future research.

Fig. 14a plots the convergence of the strain energy. As a result, the NS-FEM-T3 produces an upper bound of the strain energy. This is because the problem belongs to the force-driven problem class [64]. The convergence rates in terms of the strain energy norm and stress intensity factor KI for different numerical methods are displayed in Figs. 15a, and 16a, respectively. It is seen that the sES-FEM-T3-5 stands out clearly for this problem. The sES-FEMT3-5 solutions are closer to the reference values compared to the FEM-T3, NS-FEM-T3, and ES-FEM-T3, and slightly more accurate than the standard XFEM-Q4 results. Additionally, as illustrated in Figs. 14a and 16a, it is shown that the sES-FEM-T3-5 is a strong competitor to the six-node quarter-point singular element (sFEMT6) [49], which has been widely used in computational fracture mechanics. It is worth emphasizing that in the sES-FEM-T3-5, a mesh generation is much simpler than the sFEM-T6 because it only relies on a basic mesh of three-node linear triangular (T3) elements with a set of nodes added only on the edges of T3 elements connected to the crack tip. Some advantages of the singular ES-FEM compared to higher-order sFEM-T6 element can be found in [31– 34]. Next we show the effect of using an adaptive procedure. The mesh-refinement is started from a uniform mesh as shown in Fig. 12b. This adaptive analysis process is performed until the relative error in energy norm is less than the user specified relative error of 5%. Fig. 13 illustrates several steps of refinement. It is evident that the high density of elements is concentrated in the vicinity of the crack tip. As can be seen from Figs. 14b, 15b, and 16b, the computational effect is decreased significantly. As expected, the experimental convergence rate reached the optimal value of 1.0 and is slightly higher using the adaptive sES-FEM-T3-5. 6.3. Crack propagation of an edge-crack under mixed-mode loading In this example, we consider the edge crack geometry subjected to a shear load s = 1 N/mm2 as shown in Fig. 17a. The material parameters are Young’s modulus E = 3  107 N/mm2 and Poisson’s ratio m = 0.25. The exact stress intensity factors are given as

K I ¼ 34 mm3=2 ;

K II ¼ 4:55 mm3=2 :

ð48Þ

6.2. Plate with an single edge-crack under uniaxial tension Consider a plate under uniaxial tension r = 1 N/m2 as shown in Fig. 12a. The plate’s dimension (a = 0.3, H = 2, b = 1) is measured in meters. The material parameters are Young’s modulus E = 103 N/ m2 and Poisson’s ratio m = 0.3. Plane strain conditions are assumed. The empirical mode I SIF (a/b 6 0.6) is

pffiffiffiffiffiffi K I ¼ C f r ap;

ð46Þ

where a is the crack length, b is the plate width and the correction factor Cf is given as:

a a 2 a 3 a 4 C f ¼ 1:12  0:231  21:72 þ 30:39 : þ 10:55 b b b b ð47Þ First we study the performance of the present element with uniform meshes. The results of the sES-FEM-T3-5 are compared with those of the FEM-T3, NS-FEM-T3, ES-FEM-T3 and XFEM-Q4.

(a)

(b)

Fig. 25. Comparison between presented and experimental crack growth trajectories for case I: (a) prediction and (b) experience [52,53].

H. Nguyen-Xuan et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 252–273

269

(a) 624 DOFs (initial step 0)

(b) 2422 DOFs (step 3)

(c) 2396 DOFs (step 14)

(d) 2420 DOFs (step 21) Fig. 26. Crack growth trajectory for the PMMA beam (case II) with the crack increment at 0.15.

Fig. 18 shows the convergence of the strain energy and the stress intensity factors. It is again seen that the sES-FEM-T3-5 results are more accurate than those of the other methods. The convergence rate of the sES-FEM-T3-5 is asymptotically suboptimal to 0.5 which matches the theoretically exponential value of the strong singularity (namely, the r1/2 displacement).

For adaptive analysis, several steps of refinement meshes are depicted in Fig. 19. It can be seen that the refinement is performed around the crack tip. In addition, the refinement is also performed partly at two corners of the clamped boundary. The strain energy, the convergence rate and the stress intensity factors are shown in Fig. 20. The numerical results confirm again the reliability of the

270

H. Nguyen-Xuan et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 252–273

(e) 2530 DOFs (step 24)

(f) A zoom of deformation at step 21 Fig. 26. (continued)

present approach. As expected, the adaptive mesh-refinement leads to an improved convergence rate which is asymptotically optimal (tends to 1.0). In the following problems, we will use this element to predict the crack propagation trajectory for LEFM problems. Let us consider the crack propagation under mixed-mode loading. For crack propagation simulation, the maximum circumferential stress criterion is adopted [51]. The adaptive meshing technique is used such that the mesh is refined based on the error indicator and refinement strategy presented in Section 5. This helps to enhance the accuracy of the crack propagation trajectory whilst maintaining a reasonable computational cost. In each propagation step, we assume that the adaptive analysis process is performed until the relative error in energy norm is less than the target error (about 11– 12%). The crack increment at each step is assumed to be 0.35 mm (about 50% of the initial crack size). Fig. 21 shows meshes at the initial step (step 0), several intermediate steps, step 10, the crack propagation trajectory and the deformed shape. It is observed that the present element produces very similar results to the crack growth path when compared with published methods in [54]. Also, the present results are compared with the standard XFEM-Q4 ones using 90986 DOFs. It is evident from Fig. 22 that the proposed method is a good competitor to the standard XFEM-Q4 for the crack propagation simulation. 6.4. Crack propagation of a panel with rivet holes (PMMA beam) This example simulates the crack propagation of a polymethyl methacrylate (PMMA) beam. The beam has length L = 20, width W = 8. Three holes and their radius are given in Fig. 23. The PMMA beam is subjected to a concentrated load at the center of the top edge. The parameters are given as: Young’s modulus

E = 29  106 psi, Poisson’s ratio m = 0.3, load P = 1 lb. We consider two values for the initial crack length a and its distance b from the left side of the beam as given in the table in Fig. 23. This benchmark problem has been numerically and experimentally investigated by many authors [52–54]. For Case I, Fig. 24 illustrates various mesh steps for predicting the crack growth trajectory. As seen, a mesh refinement is performed not only in the vicinity of the crack tip, but also around holes and even at the load point and the constraint points. In numerical experiments, the crack increment being between 20% and 50% of the initial crack length [52,53] can be taken at 0.5 for seven first steps and 0.3 for the two latter steps. It is found that the crack path eventually grows into the second hole after 9 steps. A comparison between presented and actual crack paths is given in Fig. 25. The predicted crack path shows excellent agreement with several numerical and experimental results published in [52,53]. For Case II, it is more challenging than Case I. As shown in [52,53], the crack increment size should be chosen in the small range of 0.05–0.3 in., which can yield accurately the crack path in the region where the crack intersects the hole. In our work, a priori crack increment length can be chosen between 10% and 20% of the initial crack length, depending on the KII/KI ratio. As illustrated below, the crack increment can be fixed at 0.15 (about 10% of the initial crack size). Fig. 26 shows the crack propagation and the deformed shape. It is found from numerical experiments for the present method that the crack eventually grows into the second hole after 25 steps. In this case, the crack increment is taken to be fixed at 0.15. A comparison between the predicted, actual and XFEM-Q4 crack paths is displayed in Fig. 27. It is seen that the present method can produce very accurate results of the observed crack trajectory in [52,53]. It is worth emphasizing that our

H. Nguyen-Xuan et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 252–273

271

(b)

(a)

Fig. 28. Model of spalling of brittle plate problem and experimental result [65].

can be encountered in engineering problems such as offshore, flint knapping, and thin film structures [65]. Several meshes for simulating the crack propagation using adaptive procedure are illustrated in Fig. 29 along with crack trajectory. A comparison between the predicted and XFEM-Q4 crack paths is displayed in Fig. 30. It is again noted that the adaptive ES-FEM-T3-5 using about 2500 DOFs can reach quite accurate results compared with XFEM-Q4 using a structured Q4 mesh of 77726 DOFs. 7. Conclusions

(c) Fig. 27. Comparison between presented, experimental, and XFEM-Q4 crack growth trajectories for case II: (a) prediction; (b) experience [52,53]; (c) XFEM-Q4 (with 143554 DOFs).

method only uses meshes with about 2600 DOFs and the present result is well compared with that of XFEM-Q4 using a structured Q4 mesh of 143554 DOFs, and that reported in [54], in which their method required a final mesh of 8992 DOFs. 6.5. Edge cracking and spalling of brittle plate Finally, we consider a brittle PMMA plate (40 mm  60 mm) with an edge angled crack is subjected to a uniform traction r = 1 N/mm2 over its upper edge region of length 2h (h = 2 mm) as shown in Fig. 28. Other parameters are given as the initial crack a0 = 5 mm, d0 = 15.6 mm. This example aims to simulate spalling phenomena of a brittle plate under a load applied along a region of the edge. The model

In this paper, we have developed successfully an adaptive singular ES-FEM approach to model the singularity in solids, using a base mesh of triangles that can be generated automatically for complicated geometries. The adaptive meshing process is based on the algorithm of latest vertex bisection and an error indicator derived from smoothed (recovered) stresses of the present method. The singular ES-FEM formulation relies on adding a node on each edge of linear triangular elements directly connected to the singular point. The displacement interpolation was then performed on the sT5 elements. The present formulation uses the shape functions (not the derivatives) to compute the stiffness matrix and thus the complexity of the computations is significantly reduced. Numerical examples were provided to validate the high accuracy and the effectiveness of the present method. Some concluding remarks are pointed out as follows: – The present numerical procedure is simple to implement into the existing finite element packages. – The proposed adaptive algorithm gained high effectiveness for problems tested. – Both weak and strong singular problems were solved efficiently.

272

H. Nguyen-Xuan et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 252–273 60

50

40

30

20

10

(a) 118 DOFs (initial step 0)

XFEM−Q4 Adaptive sES−FEM−T3−5

(b) 1928 DOFs (step 1)

0

0

10

20

30

40

Fig. 30. Model of spalling of a brittle plate problem: comparison between presented and XFEM-Q4 approaches.

field, similar to Ródenas et al. [57,58]. It will also be of our great interest in studying a coupling technique of the basic functions derived from the isogeometric approach [66,67] into the present method to simulate curved crack paths. Acknowledgements The support of the Vietnam National Foundation for Science and Technology Development (NAFOSTED; Grant No. 107.02-2012.17) and of the Research Scholar Program from the University of Cincinnati (USA) is gratefully acknowledged.

(c) 2130 DOFs (step 3)

(d) 2420 DOFs (step 7) References

(e) 2280 DOFs (step 9) Fig. 29. Crack growth trajectory for the edge crack and spalling brittle plate with the crack increment at 5 mm.

– The method can simulate very well the crack propagation under mixed mode behaviors and the geometry with holes. – The present method showed very good agreement compared with several other published methods and experiments. Future research will be necessary to further improve the recovery process by using enriched moving least squares recovery as in [55,56] or splitting the smooth and the singular parts of the stress

[1] G.R. Liu, T. Nguyen-Thoi, Smoothed Finite Element Methods, CRC Press, Taylor and Francis Group, New York, 2010. [2] J.S. Chen, C.T. Wu, S. Yoon, Y. You, A stabilized conforming nodal integration for Galerkin mesh-free methods, Int. J. Numer. Methods Engrg. 50 (2001) 435– 466. [3] G.R. Liu, On G space theory, Int. J. Comput. Methods 6 (2009) 257–289. [4] G.R. Liu, G.Y. Zhang, A normed G space and weakened weak (W-2) formulation of a cell-based smoothed point interpolation method, Int. J. Comput. Methods 6 (2009) 147–179. [5] M.B. Liu, G.R. Liu, Z. Zong, An overview on smoothed particle hydrodynamics, Int. J. Comput. Methods 5 (2008) 135–188. [6] T.J.R. Hughes, The Finite Element Method, Prentice-Hall, Englewood Cliffs, NJ, 1987. [7] G.R. Liu, S.S. Quek, The Finite Element Method: A Practical Course, Butterworth Heinemann, Oxford, 2003. [8] G.R. Liu, Mesh-free Methods: Moving Beyond the Finite Element Method, second ed., CRC Press, Boca Raton, 2009. [9] G.R. Liu, K.Y. Dai, T. Nguyen-Thoi, A smoothed finite element for mechanics problems, Comput. Mech. 39 (2007) 859–877. [10] G.R. Liu, T. Nguyen-Thoi, H. Nguyen-Xuan, K.Y. Lam, A node based smoothed finite element method (NS-FEM) for upper bound solution to solid mechanics problems, Comput. Struct. 87 (2009) 14–26. [11] T. Nguyen-Thoi, G.R. Liu, H. Nguyen-Xuan, Additional properties of the nodebased smoothed finite element method (NS-FEM) for solid mechanics problems, Int. J. Comput. Methods 6 (2009) 633–666. [12] 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, J. Sound Vib. 320 (2009) 1100–1130. [13] T. Nguyen-Thoi, G.R. Liu, K.Y. Lam, G.Y. Zhang, A Face-based Smoothed Finite Element Method (FS-FEM) for 3D linear and nonlinear solid mechanics problems using 4-node tetrahedral elements, Int. J. Numer. Methods Engrg. 78 (2009) 324–353. [14] G.R. Liu, T. Nguyen-Thoi, K.Y. Dai, K.Y. Lam, Theoretical aspects of the smoothed finite element method (SFEM), Int. J. Numer. Methods Engrg. 71 (2007) 902–930. [15] H. Nguyen-Xuan, S. Bordas, H. Nguyen-Dang, Smooth finite element methods: convergence, accuracy and properties, Int. J. Numer. Methods Engrg. 74 (2008) 175–208.

H. Nguyen-Xuan et al. / Comput. Methods Appl. Mech. Engrg. 253 (2013) 252–273 [16] G.R. Liu, H. Nguyen-Xuan, T. Nguyen-Thoi, A theoretical study of S-FEM models: properties, accuracy and convergence rates, Int. J. Numer. Methods Engrg. 84 (2010) 1222–1256. [17] G.R. Liu, A G space theory and weakened weak (W2) form for a unified formulation of compatible and incompatible methods: Part I Theory, Int. J. Numer. Methods Engrg. 81 (2009) 1093–1126. [18] K.Y. Dai, G.R. Liu, Free and forced vibration analysis using the smoothed finite element method (SFEM), J. Sound Vib. 301 (2007) 803–820. [19] H. Nguyen-Xuan, T. Rabczuk, S. Bordas, J.F. Debongnie, A smoothed finite element method for plate analysis, Comput. Methods Appl. Mech. Engrg. 197 (2008) 1184–1203. [20] N. Nguyen-Thanh, T. Rabczuk, H. Nguyen-Xuan, S. Bordas, A smoothed finite element method for shell analysis, Comput. Methods Appl. Mech. Engrg. 198 (2008) 165–177. [21] X.Y. Cui, G.R. Liu, G.Y. Li, X. Zhao, G.Y. Zhang, G. Zheng, Analysis of plates and shells using edge-based smoothed finite element method, Comput. Mech. 45 (2009) 141–156. [22] H. Nguyen-Xuan, T. Rabczuk, N. Nguyen-Thanh, T. Nguyen-Thoi, S. Bordas, A node-based smoothed finite element method with stabilized discrete shear gap technique for analysis of Reissner-Mindlin plates, Comput. Mech. 46 (2010) 679–701. [23] H. Nguyen-Xuan, Loc V. Tran, T. Nguyen-Thoi, H.C. Vu Do, Analysis of functionally graded plates using an edge-based smoothed finite element method, Comput. Struct. 93 (2011) 3019–3039. [24] H. Nguyen-Xuan, G.R. Liu, C. Thai-Hoang, T. Nguyen-Thoi, An edge-based smoothed finite element method with stabilized discrete shear gap technique for analysis of Reissner-Mindlin plates, Comput. Methods Appl. Mech. Engrg. 199 (9–12) (2010) 471–489. [25] H. Nguyen-Xuan, G.R. Liu, T. Nguyen-Thoi, C. Nguyen Tran, An edge-based smoothed finite element method (ES-FEM) for analysis of two-dimensional piezoelectric structures, Smart. Mater. Struct. 18 (2009) 065015. [26] S. Bordas, T. Rabczuk, H. Nguyen-Xuan, P. Nguyen Vinh, S. Natarajan, T. Bog, Q. Do Minh, H. Nguyen Vinh, Strain smoothing in FEM and XFEM, Comput. Struct. 88 (23–24) (2010) 1419–1443. [27] S. Bordas, S. Natarajan, P. Kerfriden, C.E. Augarde, D.R. Mahapatra, T. Rabczuk, S.D. Pont, On the performance of strain smoothing for quadratic and enriched finite element approximations (XFEM/GFEM/PUFEM), Int. J. Numer. Methods Engrg. 86 (4–5) (2011) 637–666. [28] N. Vu-Bac, H. Nguyen-Xuan, L. Chen, S. Bordas, P. Kerfriden, R.N. Simpson, G.R. Liu, T. Rabczuk, A noded-based smoothed XFEM for fracture mechanics, CMES 73 (2011) 331–356. [29] G.R. Liu, L. Chen, T. Nguyen-Thoi, K. Zeng, G.Y. Zhang, A novel singular nodebased smoothed finite element method (NS-FEM) for upper bound solutions of fracture problems, Int. J. Numer. Methods Engrg. 83 (2010) 1466–1497. [30] G.R. Liu, N. Nourbakhshnia, Y.W. Zhang, A novel singular ES-FEM method for simulating singular stress fields near the crack tips for linear fracture problems, Engrg. Fract. Mech. 78 (2011) 863–876. [31] L. Chen, G.R. Liu, N. Nourbakhashnia, K. Zeng, A singular edge-based smoothed finite element method (ES-FEM) for bimaterial interface cracks, Comput. Mech. 45 (2010) 109–125. [32] G.R. Liu, N. Nourbakhshnia, L. Chen, A novel general formulation for singular stress field using the ES-FEM method for the analysis of mixed-mode crack, Int. J. Comput. Methods 7 (2010) 191–214. [33] N. Nourbakhshnia, G.R. Liu, A quasi-static crack growth simulation based on the singular ES-FEM, Int. J. Numer. Methods Engrg. 88 (2011) 473–492. [34] L. Chen, G.R. Liu, Y. Jiang, K. Zeng, J. Zhang, A singular edge-based smoothed finite element method (ES-FEM) for crack analyses in anisotropic media, Engrg. Fract. Mech. 78 (2011) 85–109. [35] T. Nguyen-Thoi, H.C. Vu-Do, T. Rabczuk, H. Nguyen-Xuan, A node-based smoothed finite element method (NS-FEM) for upper bound solution to viscoelastoplastic analyses of solids using triangular and tetrahedral meshes, Comput. Methods Appl. Mech. Engrg. 199 (45–48) (2010) 3005–3027. [36] T. Nguyen-Thoi, G.R. Liu, H.C. Vu-Do, H. Nguyen-Xuan, An edge-based smoothed finite element method (ES-FEM) for visco-elastoplastic analyses of 2D solids using triangular mesh, Comput. Mech. 45 (2009) 23–44. [37] T. Nguyen-Thoi, G.R. Liu, H.C. Vu-Do, H. Nguyen-Xuan, A face-based smoothed finite element method (FS-FEM) for visco-elastoplastic analyses of 3D solids using tetrahedral mesh, Comput. Methods Appl. Mech. Engrg. 198 (2009) 3479–3498. [38] T.N. Tran, G.R. Liu, H. Nguyen-Xuan, T. Nguyen-Thoi, An edge-based smoothed finite element method for primal-dual shakedown analysis of structures, Int. J. Numer. Methods Engrg. 82 (2010) 917–938. [39] Z.C. He, G.R. Liu, Z.H. Zhong, G.Y. Zhang, A.G. Chen, A coupled ES-FEM/BEM method for fluid–structure interaction problems, Engrg. Anal. Bound. Elem. 35 (2011) 140–147.

273

[40] B.A. Szabó, I. Babuška, Finite Element Analysis, John Wiley & Sons, New York, 1991. [41] R.K. Banerjee, The Boundary Element Methods in Engineering, Mcgraw-Hill College, 1994. Rev Sub Edition. [42] T. Belytschko, Y.Y. Lu, L. Gu, Crack propagation by element-free Galerkin methods, Engrg. Fract. Mech. 51 (1995) 295–315. [43] M. Duflot, H. Nguyen-Dang, A meshless method with enriched weight functions for fatigue crack growth, Int. J. Numer. Methods Engrg. 59 (2004) 1945–1961. [44] T. Rabczuk, S. Bordas, G. Zi, A three-dimensional meshfree method for continuous multiple crack initiation, propagation and junction in statics and dynamics, Comput. Mech. 40 (2007) 473–495. [45] N. Moes, J. Dolbow, T. Belytschko, A finite element method for crack growth without remeshing, Int. J. Numer. Methods Engrg. 46 (1999) 131–150. [46] G. Ventura, R. Gracie, T. Belytschko, Fast integration and weight function blending in the extended finite element method, Int. J. Numer. Methods Engrg. 77 (2009) 1–29. [47] P. Rrisvard, Elliptic Problems in Nonsmooth Domains, Pitman Publishing, Inc., Boston, 1985. [48] R.S. Barsoum, On the use of isoparametric finite elements in linear fracture mechanics, Int. J. Numer. Methods Engrg. 10 (1976) 25–37. [49] R.S. Barsoum, Triangular quarter-point elements as elastic and perfectly plastic crack tip elements, Int. J. Numer. Methods Engrg. 11 (1977) 85–98. [50] H. Nguyen-Xuan, G.R. Liu, N. Nourbakhshnia, L. Chen, A novel singular ES-FEM for crack growth simulation, Engrg. Fract. Mech. 84 (2012) 41–66. [51] N. Moes, J. Dolbow, T. Belytschko, A finite element method for crack growth without remeshing, Int. J. Numer. Methods Engrg. 46 (1999) 131–150. [52] A.R. Ingraffea, M. Grigoriu, Probabilistic Fracture Mechanics: A Validation of Predictive Capability, Report 90-8, Department of Structural Engineering, Cornell University, 1990. [53] T.N. Bittencourt, P.A. Wawrzynek, A.R. Ingraffea, J.L. Sousa, Quasi-automatic simulation of crack propagation for 2D LEFM problems, Engrg. Fract. Mech. (522) (1996) 321–334. [54] D. Azocar, M. Elgueta, M.C. Rivara, Automatic LEFM crack propagation method based on local Lepp-Delaunay mesh refinement, Adv. Engrg. Softw. 41 (2010) 111–119. [55] S. Bordas, M. Duflot, P. Le, A simple error estimator for extended finite elements, Commun. Numer. Methods Engrg. 24 (11) (2008) 961–971. [56] M. Duflot, S. Bordas, A posteriori error estimation for extended finite elements by an extended global recovery, Int. J. Numer. Methods Engrg. 76 (8) (2008) 1123–1138. [57] J.J. Ródenas, O.A. González-Estrada, J.E. Tarancón, F.J. Fuenmayor, A recoverytype error estimator for the extended finite element method based on singular + smooth stress field splitting, Int. J. Numer. Methods Engrg. 76 (2008) 545–571. [58] J.J. Ródenas, O.A. González-Estrada, P. Díez, F.J. Fuenmayor, Accurate recoverybased error upper bounds for the extended finite element framework, Comput. Methods Appl. Mech. Engrg. 199 (2010) 2607–2621. [59] O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method, fifth ed., vol. 1, Butterworth Heinemann, Oxford, 2000. [60] J. Alberty (Kiel), C. Carstensen (Vienna), S.A. Funken (Kiel), R. Klose (Kiel), Matlab implementation of the finite element method in elasticity, Computing 69 (2002) 239–263. [61] W. Dorfler, A convergent adaptive algorithm for Poisson’s equation, SIAM J. Numer. Anal. 33 (1996) 1106–1124. [62] L. Chen, Short implementation of bisection in MATLAB, in: P. Jorgensen, X. Shen, C.-W. Shu, N. Yan, (Eds.), Recent Advances in Computational Sciences – Selected Papers from the International Workshop on Computational Sciences and its Education, 2008, pp. 318–332. [63] S. Funken, D. Praetorius, P. Wissgott, Efficient implementation of adaptive p1FEM in Matlab, Preprintm, 2008. . [64] J.F. Debongnie, H.G. Zhong, P. Beckers, Dual analysis with general boundary conditions, Comput. Methods Appl. Mech. Engrg. 122 (1995) 183–192. [65] M.D. Thouless, A.G. Evans, M.F. Ashby, J.W. Hutchinson, The edge cracking and spalling of brittle plates, Acta Metall. 35 (1987) 1333–1341. [66] T.J.R. Hughes, J.A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Comput. Methods Appl. Mech. Engrg. 194 (2005) 4135–4195. [67] N. Nguyen-Thanh, H. Nguyen-Xuan, S. Bordas, T. Rabczuk, Isogeometric analysis using polynomial splines over hierarchical T-meshes for twodimensional elastic solids, Comput. Methods Appl. Mech. Engrg. 200 (21–22) (2011) 1892–1908.