Reformulation of XFEM based on PUFEM for solving problem caused by blending elements

Finite Elements in Analysis and Design 45 (2009) 806 -- 816

Contents lists available at ScienceDirect

Finite Elements in Analysis and Design journal homepage: w w w . e l s e v i e r . c o m / l o c a t e / f i n e l

Reformulation of XFEM based on PUFEM for solving problem caused by blending elements Kazuki Shibanuma ∗ , Tomoaki Utsunomiya Department of Civil and Earth Resources Engineering, Kyoto University, Japan

A R T I C L E

I N F O

Article history: Received 22 September 2008 Received in revised form 20 May 2009 Accepted 12 June 2009 Available online 10 July 2009 Keywords: XFEM PUFEM The corrected XFEM Blending elements Linear fracture mechanics

A B S T R A C T

A lack of numerical accuracy in the standard extended finite element method (XFEM) is caused by `blending elements', whose nodes are partially enriched. `The corrected XFEM' proposed by Fries showed the effective improvement of this problem with a lot of numerical results. The theoretical approach of this proposal was however not sufficiently described. In the present paper, an approximation of the XFEM is reformulated based on the concept of the partition of unity finite element method (PUFEM) approximation, which assures the numerical accuracy in the entire domain, for solving the problem of blending elements. The form of the reformulated XFEM results in the coincidence with that of the corrected XFEM. It is therefore found out that the theoretical validation of the corrected XFEM is based on the PUFEM approximation. It is also found out that the problem of the blending elements in the application to two dimensional linear fracture mechanics has been sufficiently solved for actual use by the XFEM based on the PUFEM. © 2009 Elsevier B.V. All rights reserved.

1. Introduction The partition of unity finite element method (PUFEM) [1] includes local approximations reflecting a priori knowledge about the solution in the framework of FEM by using `partition of unity (PU)'. The extended finite element method (XFEM) [2] is a meshfree method also uses PU and employs the local enrichment function, which enables the approximation allowing the reproduction of singularity or discontinuity such as crack in the local parts of the domain. In the XFEM, there are therefore many investigations in application to the crack analysis, including additional modifications and proposals such as the crack representation [3] and simplification of post processing tasks [4,5]. The features of the XFEM in two dimensional linear fracture mechanics are summarized below [6]: (a) Crack is defined as arbitrary `line' independently from the meshes without any double nodes representing the discontinuity of crack. (b) High accuracy analysis is possible by using relatively simple and coarse meshes without any sub-divided fine meshes near the crack tip.

∗ Corresponding author. E-mail address: [email protected] (K. Shibanuma). 0168-874X/$ - see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.finel.2009.06.007

(c) Crack propaga tion is modeled by only elongating the `line' defined as crack without any complicated remeshing procedures. In the XFEM, there are elements with some of their nodes being enriched and the remaining nodes being standard elements in the numerical model. These partially enriched elements are called `blending elements', by which the problems in numerical accuracy were reported [7]. Therefore, the following studies were performed in the past for improvement of the numerical accuracy of the XFEM. In crack analyses, the range of the enrichment based on asymptotic solution of a singularity near a crack tip was defined as the interior of the circle with a certain radius from crack tip, instead of only the support of nodal interpolation function containing the crack tip. As a result, it was shown that numerical accuracy was consequentially improved by enlarging the range of the enrichment area [6]. However, relatively large range of the enrichment around the crack tip causes the limitation in modeling of the enrichment close to the boundary of the model, in particular in the crack propagation process. Recently, Gracie tried to avoid the problem of the accuracy reduction caused by the blending elements by applying the assumed strain (AS) method and the discontinuous Galerkin (DG) method to the XFEM, providing a certain improvement [8]. However, in the application to the fracture mechanics, a sufficient numerical accuracy is not achieved when coarse meshes are used near the crack tip. This

K. Shibanuma, T. Utsunomiya / Finite Elements in Analysis and Design 45 (2009) 806 -- 816

is because the application of the coarse meshes to the crack analyses is one of the advantages of the XFEM. In addition, complicated processing is required for application of the DG method or the AS method to the XFEM. This method may therefore not be a realistic solution for the problem of the blending element. On the other hand, Fries proposed `the corrected XFEM' employing a new definition of enrichment [9]. The corrected XFEM showed the effective improvement of the accuracy reduction caused by the blending elements with a lot of numerical results. In addition, this method can be implemented by simple revision of the classical FEM or the standard XFEM programming codes. Therefore, this proposal is an effective method for actual use to improve the problem caused by the blending elements. This proposal is however based on a simple redefinition of enrichment function satisfying to be zero on the boundary between standard elements and blending elements by using a rump function in the framework of the approximation in the standard XFEM. That is, the usefulness of this proposal is unfortunately based only on the results of the inductive evaluations because the theoretical approach was not sufficiently described, though the problem of blending elements was clarified by the numerical evaluations on the interpolation error [7]. As a theoretical approach to solve the problem caused by the blending elements, an application of the PUFEM to the XFEM will be useful, because it is known that the approximate accuracy based on PUFEM is assured in the entire numerical domain [1]. However, the compositions of the approximation of the PUFEM are different from that of the XFEM, and the relationship between the XFEM and the PUFEM is differently defined by the respective researches [10–12]. From the above background, a reformulation of the XFEM is required based on the useful theory to solve the accuracy reduction caused by the blending elements. Therefore, we try to reformulate the XFEM based on the PUFEM in this paper in order to solve the problem caused by the blending elements, showing the relationship between the XFEM and the PUFEM. In addition, we also try to find out the reason why the corrected XFEM can improve the accuracy reduction caused by the blending elements. An outline of this paper is as follows: In Section 2, we describe the definition of the existing standard XFEM and its approximation accuracy including the problem of blending elements. In Section 3, we propose a reformulation of the XFEM based on concept in the approximation of the PUFEM, which assures the approximation accuracy of the entire numerical domain in order to solve the problem of blending elements. Then, the approximate accuracy of this reformulated XFEM for one dimensional problem is evaluated to clarify the potential ability based on a typical approach using convergence study, including an assessment of the theoretical background of the corrected XFEM. In Section 4, an approximation formula of displacement field in two dimensional linear fracture mechanics using the reformulated XFEM is described. In Section 5, some numerical results by using the reformulated XFEM in two dimensional linear fracture mechanics are presented. Finally, in Section 6, the conclusions of this study are summarized.

807

2.1. Approximation of displacement field in standard XFEM The standard XFEM approximation uap (x) has the form uap (x) = ustd (x) + uenr (x)

(1)

where ustd (x) is the standard FEM approximation and uenr (x) is the approximation using the enrichment. The standard part of the approximation is  I (x)uI (2) ustd (x) = I∈N

where N is the set of all nodes, I (x) are the interpolation functions and uI are the nodal degrees of freedom associated for standard finite element approximation. The enriched part of the approximation is  I (x)(x)aI (3) uenr (x) = I∈Nenr

where (x) is called enrichment function, which is generally a basic function including a priori knowledge of the solution. Nenr is the set of nodes enriched by (x) and aI are the nodal degrees of freedom associated for enrichment. In this paper, we consider that I (x) are first order interpolation functions for simplification. 2.2. Blending elements In the approximation of the XFEM, three sets of elements are defined. The first type is a set of elements Estd , in which none of the nodes is enriched. The second type is a set of elements Eenr , in which all of the nodes are enriched. The third type is a set of elements Eblnd , in which only some of the nodes are enriched. The sub-domain composed of each set of elements is defined as

std = ∪ e e∈Estd

enr = ∪enr e e∈E

blnd



=

∪

e∈Eblnd

e

(4)

where e is the domain composed of the element e. The elements blnd are called `blending elements' [7–9]. They blend composing  enr with that in the the approximation in the enriched elements in  std standard elements in  . In the component i of the approximation uenr (x) in (3), if aIi = 1 then ⎧ std ∀x∈ ⎪ ⎪ I (x)(x) = 0,  ⎨ I (x)(x) = (x), ∀ x ∈ enr (5) uenr i (x) = ⎪ ⎪ I∈Nenr ⎩ blnd I (x)(x)
(x), ∀ x ∈  The influence of the blending elements on the approximation accuracy in the standard XFEM will be evaluated in Section 2.3 using the convergence study for a one dimensional problem.

2. Standard XFEM 2.3. Approximation accuracy of standard XFEM In this section, we describe the definition of the approximation using the XFEM based on the previous study. This section contains as follows: In Section 2.1, the definition of the approximation in the standard XFEM is described. In Section 2.2, the definition of blending elements, which affects on the accuracy in the standard XFEM, is presented. In Section 2.3, the problem in the approximate accuracy of the standard XFEM is evaluated in the one dimensional problem.

The evaluation of the approximation accuracy in the XFEM for a one dimensional problem is described below using the L2 norm of the interpolation error. The entire domain is assumed as  ={x|0 ⱕ x ⱕ 1}, then each subenr domain defined in (4) is also assumed as  = {x|0 ⱕ x ⱕ (M − 1)h}, blnd std  = {x|(M − 1)h ⱕ x ⱕ Mh} and  = {x|Mh ⱕ x ⱕ 1} as shown in Fig. 1, where h is the element size related to the number of elements

808

K. Shibanuma, T. Utsunomiya / Finite Elements in Analysis and Design 45 (2009) 806 -- 816

Ωenr

h

0

Ωblnd

(M−1)h

Ψ(x)

Ωstd

1−h

Mh

1

: enriched nodes

: enriched elements

: standard nodes

: blending elements

x Ψ(x)M(x)

: standard elements

M+1(x)

M(x)

Fig. 1. Discretized domain for one dimensional problem in the standard XFEM.

N, defined as h = 1/N. The blending element is composed of nodes M and (M+1). If there are singularities in this problem, their locations enr are supposed to be in  . The interpolation functions in which are used in the element n, composed of nodes n and (n+1), are

x − (n − 1)h n+1 (x) = h

(6)

The interpolation error is defined as (7)

The approximation accuracy in the element composing each subdomain is evaluated as follows. The approximation of the displacement field in the element n enr is (< M), which compose the sub-domain  n+1 

{I (x)uI + I (x)(x)aI }

enr

on 

(8)

I=n

nh

(n−1)h

|(x)|2 dx ⱕ Chp

{I (x)uI }+M (x)aM {(y)+ (y˜ )(x−y)}

M+1 

enr

on 

(9)

{I (x)uI } + M (x)(x)aM

on 

blnd

(12)

I=M

uap (x) = −

2aM   (y˜ ) on blnd h

on 

blnd

(10)

(13)

In the approximation, the maximum interpolation error occurs at the point x¯ where  (x¯ ) = 0. Then a Taylor series expansion of the interpolation error (x) about x = x¯ gives

(x) =

∞ 

(m) (x¯ )

(x − x¯ )m m!

= (x¯ ) +  (x¯ )(x − x¯ ) + = (x¯ ) +

1   (x¯ )(x − x¯ )2 + O(h3 ) 2

blnd

on 

1   (x˜ )(x − x¯ )2 2

(14)

where (M − 1)h ⱕ x¯ ⱕ Mh. If we let x = (M − 1)h, then ((M − 1)h) = 0, since uap (x) is the finite element interpolation of u(x). Therefore, we obtain 1 2

(x¯ ) = {(M − 1)h − x¯ }2  ((M − 1)h)

where C is a generic constant independently from h and p. The approximation of the displacement field in the blending element blnd , is n ( = M), which composes the sub-domain  uap (x) =

M+1 

m=0

If each nodal degree of freedom in (8) is assumed as un = un+1 = 0, an = an+1 = 1, the approximation uap (x) obtained as an exact solution (x). However, this can only be true in special cases. It is difficult enr in actual to reproduce the exact solution in the sub-domain  analyses due to the restriction on the consistency of the interpolation enr blnd . Thus, the on the boundary between the sub-domains  and  following inequality is defined using a constant p ( > 0) as: 

uap (x)=

Therefore, the second order derivative of the approximation uap (x) is

(x) ≡ uap (x) − u(x)

uap (x) =

Fig. 2. Approximation on a blending element for one dimensional problem in the standard XFEM.

where (M − 1)h ⱕ y˜ ⱕ Mh. Substituting (11) in (10), we obtain

x − (n − 1)h h

n (x) = 1 −

Mh

(M−1)h

1 blnd {(M − 1)h − x¯ }2 | ((M − 1)h)| on  2     1 2a ⱕ {((M − 1)h) − x¯ }2 max u (x) + M  (y) 2 h

ⱕ

(15)

and since

I=M

An approximation of the displacement field in the blending element in (10) is shown in Fig. 2. blnd can be The evaluation for the approximation accuracy in  expressed using as a reference the approach presented by Chessa [7] as follows: We let y be a point of (M − 1)h ⱕ y ⱕ Mh. The exact solution (x) is expanded in a Taylor series about x = y to obtain

(x) =

∞ 

(m) (y)

m=0

(x − y) m!

m



2

− 1)h) − x¯ }2 ⱕ 18 h2

1 8

 

(x) ⱕ (x¯ ) ⱕ h2 max u (x) +

 2aM    (y) h

on 

blnd

(17)

blnd

has the order

Thus, the interpolation error in the sub-domain  h [7]. This leads to Mh (M−1)h

(11)

(16)

it follows that:



= (y) +  (y)(x − y) + O(h ) = (y) +  (y˜ )(x − y)

1 2 {((M

|(x)|2 dx ⱕ Ch3

blnd

on 

where C is a generic constant independently from h.

(18)

K. Shibanuma, T. Utsunomiya / Finite Elements in Analysis and Design 45 (2009) 806 -- 816

The error of the approximation of the displacement field in the std element n ( > M), which composes the sub-domain  corresponds to one of the standard finite element approximation as uap (x) =

n+1 

I (x)uI on std

(19)

I=n

If the approximation uap (x) in (19) corresponds to the exact soblnd

lution (x) at the nodes as the evaluation in the sub-domain  , the nodal degrees of freedom are determined as un = ((n − 1)/h), and un+1 = (nh), respectively. Thus, the approximation uap (x) on

std in (19) can be rewritten as

uap (x) = n (x)((n − 1)h) + n+1 (x)(nh) on 

std

(20)

std

The interpolation error on  can be given by using the expansions of (nh) and (x) in a Taylor series expansion around x = (n − 1)h as

809

3. Reformulation of XFEM based on PUFEM In Section 2, the limitation in the approximate accuracy of the standard XFEM was described. It was therefore necessary to reconsider the approximation of the XFEM to solve the problem of the standard XFEM. It is known that the approximate accuracy based on PUFEM is assured in the entire numerical domain [1]. Therefore, we propose a reformulation of the approximation of the XFEM based on the PUFEM so as to solve the lack of accuracy in the blending elements. The contents of this section are as follows: In Section 3.1, the approximation of the displacement field in the PUFEM is described, based on the work of Melenk and Babuska [1]. In Section 3.2, the approximation of displacement field in the XFEM including enrichment is reformulated based on the concept of the approximation in the PUFEM. In addition, we try to clarify the reason why the corrected XFEM proposed by Fries [9] can improve the accuracy reduction caused by the blending elements.

(x) = {n (x)((n − 1)h) + n+1 (x)(nh)} − (x)

=

3.1. Approximation of displacement field in the PUFEM

n (x)((n − 1)h)

+ n+1 (x) −

∞ 

(m) ((n − 1)h)

m=0 ∞ 

m=0



(m)

m

(nh − (n − 1)h) m!

(x − (n − 1)h)m ((n − 1)h) m!

 std

on 



j

1 =  ((n − 1)h)(x − (n − 1)h)(nh − x) + O(h3 ) 2

(21)

(see Ref. [13] for details the procedure in (21)). std Thus, the interpolation error in the sub-domain  has the order h2 . This leads to  nh std |(x)|2 dx ⱕ Ch5 on  (22) (n−1)h

where C is a generic constant, which is independent from h. Due to (9), (18) and (22), the L2 norm on the entire domain  by using the standard XFEM is evaluated as

 uap (x) − u(x) L2 () =

|uap (x) − u(x)|2 dx



1/2

enr

+

|uap (x) − u(x)| dx

blnd

 +

std

j L∞ , ⱕ C0

(27)

where j are the support of j and C0 is a constant. Note that j is PU but their definitions are different from the interpolation functions in the FEM. ap Using local approximation spaces Vj on j , the PUFEM space Vap is defined as ⎧ ⎫ ⎨  ⎬  ap ap  ap ap (28) j Vj = j vj  vj ∈ Vj Vap ≡ ⎩ ⎭ j

j

ap

Assume that the local approximation function Vj on j have the following approximation properties for the exact solution u:

1/2

j

ⱕ C{(M − 1) · hp + h3 + (N − M)h5 }1/2

(23)

In the case of p ⱖ 3, the convergence rate of the norm in (23) is determined as 32 due to the approximation accuracy in the subL2

Due to (29), the L2 norm on the entire domain  by using the PUFEM is evaluated as ⎧ ⎫1/2 2 ⎪  ⎪ ⎨ ⎬   j (vap − u) uap − u L2 () =    j ⎪ ⎪ ⎩ ⎭  j

blnd

domain  . On the other hand, in the case of p < 3, the convergence rate is determined as p/2 ( < 3/2) due to the approximation enr accuracy in the sub-domain  . Therefore, there is obviously no case that the convergence rate is more than 32 , and we obtain uap (x) − u(x) L2 () ⱕ Ch3/2

(29)

Then, the approximation of displacement field uap in the PUFEM space Vap on  is defined as  j vap ∈ Vap (30) uap = j

|uap (x) − u(x)|2 dx |uap (x) − u(x)|2 dx

(26)

ap

2



card{j|x ∈ j } ⱕ N0

vj − u L2 (j ∩) ⱕ j

 =

In this section, the approximation in the PUFEM is described, based on the work of Melenk and Babuska [1]. Assume the functions j , satisfying  j = 1 on  (25)

(24)

ⱕ

ⱕ

⎧ ⎨

N0

⎩

ⱕ



N0

⎩



 j

N0 C 0

⎫1/2 ⎬

ap j (vj

− u) 2L2 () ⎭

ap j (vj

− u) 2L2 ( ∩) j ⎭

j

⎧ ⎨

L2

The convergence rate on the norm of the interpolation error in the FEM is 2. However, the convergence rate on the L2 norm of the interpolation error in the standard XFEM as shown in (24) is lower than that of the FEM by 12 at least.

L2 ()

⎧ ⎨ ⎩

j

⎫1/2 ⎬

2j ⎭

⎫1/2 ⎬

(31)

810

K. Shibanuma, T. Utsunomiya / Finite Elements in Analysis and Design 45 (2009) 806 -- 816

Ω1

Ω0

Ωenr

Ωblnd

Ωstd 0

0

h

(M−1)h Mh

1−h

1

: enriched nodes

: enriched elements

: standard nodes

: blending elements

x

1=1

0=1

0=0

1=0

h

(M−1)h

Mh

1−h

1

: enriched nodes

: enriched elements

: standard nodes

: blending elements

x

: standard elements

: standard elements

Fig. 4. Schematic description of the functions 0 (x) and 1 (x).

Fig. 3. Discretized domain for one dimensional problem in the reformulated XFEM.

Due to (31) in the approximation based on the PUFEM, it is shown that the approximate accuracy is assured in the entire domain  according to the local approximate accuracy on j . See Ref. [1] for details of the concept of the PUFEM approximation. 3.2. Formulation of XFEM based on PUFEM and its evaluation As shown in Section 3.1, the approximate accuracy in the PUFEM, which is different from the standard XFEM, is assured in the entire numerical domain. Then, in this section, an approximation of the XFEM with enrichment is reformulated based on the concept of the approximation of the PUFEM. Let us consider the one dimensional problem that corresponds to the case of the standard XFEM in Section 2.3. Thus, the entire domain is assumed as  = {x|0 ⱕ x ⱕ 1}, then each sub-domain is enr blnd = {x|(M − 1)h ⱕ x ⱕ Mh} defined as  = {x|0 ⱕ x ⱕ (M − 1)h},  std and  = {x|Mh ⱕ x ⱕ 1}. If there are singularities in this problem, enr their locations are supposed to be in  . Let us redefine two new domains 0 and 1 shown in Fig. 3 by enr blnd std and  as using sub-domain  , 

0 = std ∪ blnd

(32)

1 = enr ∪ blnd ap

Corresponding to the local approximation function vj in the PUFEM described in the Section 3.1, the approximations of displaceap ap ment field v0 (x) on 0 and v1 (x) on 1 , are, respectively, defined as  ap I (x)uI (33) v0 (x) = I ap v1 (x) =



I (x)(uI + (x)aI )

(34)

I

When the approximation in the domain 0 consists of only one ap function v0 (x) in (33), the interpolation error in the domain 0 is the same as that of the classical finite element approximation, as shown in the following expression corresponding to (19): ap

v0 − u L2 (0 ∩) ⱕ Ch2

(=0 )

(35)

When the approximation in the domain 1 consists of only one ap function v1 (x) in (34), the interpolation error in the domain 1 can be also defined using a constant q ( > 0) as ap

v1 − u L2 (1 ∩) ⱕ Chq

(=1 ) ap

(36)

In the approximation v1 (x) of (34), it is possible to obtain the interpolation to the exact solution using the enrichment function (x). Therefore, the convergence order of q in (36) can be made more than the one in the classical finite element approximation in (35). Thus, we can assume that q > 2.

Based on the above discussion, in the following statements, the reformulation of the XFEM based on the PUFEM is defined in the entire numerical domain  and its approximation accuracy is evaluated using the L2 norm of the interpolation error. The functions 0 (x) and 1 (x) are defined below when they are active on 0 and 1 , respectively ⎧ enr 0 on  ⎪ ⎪ ⎪ ⎨ x − Mh blnd on  0 (x) = h ⎪ ⎪ ⎪ ⎩ std 1 on  ⎧ enr 1 on  ⎪ ⎪ ⎪ ⎨ x − Mh blnd on  1 (x) = 1 − (37) h ⎪ ⎪ ⎪ ⎩ std 0 on  where the functions 0 (x) and 1 (x) are PU, i.e. the following is satisfied:

0 (x) + 1 (x) = 1

(38)

These functions 0 (x) and 1 (x) are illustrated in Fig. 4. The cardinal number in (26) is therefore N0 = 2 in this case. By using (33), (34) and (37), a reformulated approximation of uap (x) on  in the XFEM based on the PUFEM is proposed as ap

ap

uap (x) = 0 (x)v0 (x) + 1 (x)v1 (x)

(39) ap

ap

In this approximation, local approximations v0 (x) and v1 (x) are blnd

. That is, the approximation of the displacement both active on  field of blending elements in the reformulated XFEM is defined in the common part of the domains 0 and 1 using PU based on the PUFEM, as shown in Fig. 5. It is therefore noted that the approximation in the reformulated XFEM is different from the one in the standard XFEM, as shown in Fig. 2. Due to (31), (35) and (36), the L2 norm on the entire domain  by using the reformulated XFEM in (39) is evaluated as ⎧ ⎫1/2 ⎨ ⎬  2 j uap − u L2 () ⱕ N0 C0 ⎩ ⎭ j

√ = 2C0 {(Ch2 )2 + (Chq )2 }1/2

(40)

The convergence rate q can be assumed to be q > 2, as described in (36). Thus, we obtain uap − u L2 () ⱕ Ch2

(41)

Due to (41), it is shown that the insufficient convergence rate on the L2 norm of the interpolation error in the standard XFEM is overcome in the reformulated XFEM. It is therefore found that the reformulated approximation of the XFEM in (39) provides a potential ability to solve the problem of the lack of accuracy for blending elements in the standard XFEM.

K. Shibanuma, T. Utsunomiya / Finite Elements in Analysis and Design 45 (2009) 806 -- 816

v0ap(x)

0(x)

811

v1ap(x)

1(x)

Ψ(x) Ψ(x)M(x) Ψ(x)M+1(x) 1(x)

0(x)

(M−1)h

Mh

M(x)

M+1(x)

(M−1)h

Mh

M(x)

(M−1)h

M+1(x)

(M − 1) h

Mh

Mh

Fig. 5. Approximation on a blending element for one dimensional problem in the reformulated XFEM.

In this section, the reformulation of the XFEM has been proposed using the one dimensional problem. The reformulated XFEM can be extended easily to higher dimensional problems by using the same concept of approximation. According to above discussion, the XFEM has been reformulated on the basis of the PUFEM, showing the relationship between the XFEM and the PUFEM. It is therefore found out that the approximation of the standard XFEM is not exactly based on that of the PUFEM. Fig. 6. Domain and boundaries for two dimensional linear fracture mechanics.

3.3. Theoretical background of the corrected XFEM As described in the introduction in this paper, it was recently reported that the problem of blending elements is effectively improved by employing `the corrected XFEM' proposed by Fries. In this section, the theoretical background of the corrected XFEM is considered. Substituting (33), (34) and (37) into (39), which is the approximation of the reformulated XFEM, we obtain   I (x)uI + 1 (x) I (x)(x)aI (42) uap (x) = I

This form of approximation (42) results in the coincidence with that of the corrected XFEM [9]. However, the approach of the formulation of the corrected XFEM was different from the reformulated one presented in this study. That is, the corrected XFEM was based on a simple redefinition of enrichment function satisfying zero in the blnd std and  by using a rump function in the boundary between  framework of the approximation in the standard XFEM. According to the above discussion, the theoretical validation of the corrected XFEM is found out based on the concept of the PUFEM approximation. Therefore, the corrected XFEM and the reformulated XFEM in this paper may be unified as a `PUFEM-based XFEM'. 4. Application to linear fracture mechanics In this section, the approximation of displacement field in two dimensional linear fracture mechanics by the PUFEM-based XFEM is described. 4.1. Governing equations

on u

r · n = t¯

on t

(43)

⑀ = 12 (∇u + (∇u) ) in 

(44)

r = c : ⑀ in 

(45)

(46)

where n is the unit normal vector, u¯ the specified displacement vector and t¯ the specified traction vector. According to these equations, the weak form of governing equation of a linear elastic problem is expressed as  

⑀ : r d = u · t¯ d (47) 

t

4.2. Approximation of displacement field In modeling of crack analyses of the XFEM, we use two types of enrichment functions as 1 on + (48) H(x) = −1 on − √

1 (x) = r cos √

A two dimensional linear elastic body with domain  and boundary  is considered as shown in Fig. 6. Under the assumption of absence of body forces, the equation of equilibrium, the relation of strain–displacement and the relation of stress–strain for a linear elastic body are

T

u = u¯

r · n = 0 on c

I

∇ · r = 0 in 

where is the stress tensor,  the strain tensor, u the displacement vector and c the elastic constant tensor. The boundary  is decomposed into the complementary sets u , t and c as

4 (x) = r sin

2

2

,

√

2 (x) = r sin

sin

2

,

√

3 (x) = r cos

2

sin , (49)

where H(x) is the Heaviside function used for the enrichment in the presence of strong discontinuities. + and − are the upper and lower sides of the domain defined on the supports I of node I containing the Heaviside enrichment [14]. It is noted that this Heaviside enrichment dose not lead to the problem caused by blending elements even in the standard XFEM. This is because the Heaviside enrichment function is a constant in the blending elements and is not assumed as an exact solution described in the Section 4. Detailed discussion of this topic is found in Ref. [9].

812

K. Shibanuma, T. Utsunomiya / Finite Elements in Analysis and Design 45 (2009) 806 -- 816

J

J

J

J

J

J

J

J

J

J

J

J

J

J

J

J

J

J

J

J

J

C

C

C

C

C

C

C

D

J

J

J

J

J

J

J

J

J

C

C

C

C

C

C

C

RC C

Crack Tip

1

C

Crack Tip 0 Fig. 8. The contour of values of C (x).

Fig. 7. Modeling of the sets of nodes J and C.

Therefore, the form of the Heaviside enrichment in the standard XFEM is useful for the reduction of the degrees of freedom. k (x) (k = 1, . . . ,4) are the set of singular enrichment functions, developed by Fleming et al. [15], based on the near-field asymptotic solution and (r, ) are the local polar co-ordinates at the crack tip. By using the above enrichment functions, the approximation ap ap functions v0 (x) and vC (x) of the displacement field are defined as   ap I (x)uI + I (x)H(x)bI v0 (x) = I ap

vC (x) =



⎛

I∈J

I (x) ⎝uI +

I

4 

⎞

k (x)ckI ⎠ +

Nodal degrees of freedom ( : 12, : 10, : 4, no marking: 2)

RC 

I (x)H(x)bI

ΩC

(50)

I∈J

k=1

Crack Tip

where bI and ckI (k = 1, . . . ,4) are the nodal degree of freedoms related to the enrichment function H(x) and k (x) (k = 1, . . . ,4), respectively. In the application of the PUFEM-based XFEM to a two dimensional problem, an extension of the functions j (x) of (37) in one dimensional problem is required. Then, the extension to two dimensional problem j (x) is defined as

C (x) =



4.3. Modeling of nodal degrees of freedom and numerical integration

I (x)

I∈C

0 (x) = 1 − C (x)

(51)

The functions 0 (x) and C (x) are PU, satisfying

0 (x) + C (x) = 1

In order to consider the number of nodal degrees of freedom, the form of (42) in one dimensional problem is useful. Thus, substituting (50) and (51) into (55), which is the approximation in the crack analysis using the PUFEM-based XFEM, we obtain

(52) ap v0 (x)

ap vC (x)

and in (50) Note that the approximation functions are defined on the domains 0 as the support of 0 (x) and on the domain C as the support of C (x), respectively. J and C in (50) and (51) are the sets of nodes adding the enrichments and are defined as J = {I ∈ N|I ∩ D
, I ∈/ C}

(53)

C = {I ∈ N||xI − xtip | ⱕ RC }

(54)

where D is the geometry of the crack and RC the radius using the definition of the nodal set C. D and C are treated as input data in numerical analyses [6,14]. An example of the sets of nodes J and C near the crack is shown in Fig. 7. The contours of values of C (x) that correspond to the Fig. 7 are shown in Fig. 8. From (48)–(54), the approximation of displacement field uap (x) on  in two dimensional linear fracture mechanics by using the PUFEM-based XFEM is defined as ap

Fig. 9. The supports of C and the nodal degrees of freedom.

ap

uap (x) = 0 (x)v0 (x) + C (x)vC (x)

(55)

The first and second terms of the right hand side in (55) are active on the domain 0 and on the domain C , respectively.

uap (x) =



I (x)uI +

I

+ C (x)



I (x)H(x)bI

I∈J

 I

I (x)

4 

k (x)ckI

(56)

k=1

The third term of the right hand side in (56) is active on the support C of C (x). This requires that the number of nodal degrees of freedom in the PUFEM-based XFEM increases than that of the standard XFEM if their enrichment areas are kept in the same ranges [9]. That is, the nodal degrees of freedom at the nodes which compose the blending elements but not belong to the set C are increased, as compared with that of the standard XFEM. The support C and the number of nodal degree of freedoms of each node are shown in Fig. 9, and their corresponding enrichments in Figs. 7 and 8. The numerical integration procedure in each element is carried out by using the standard Gauss quadrature in the FEM except elements that are cut by the crack discontinuity. In the elements cut by the crack, the sub-division of the domain is required. This procedure corresponds to that of the standard XFEM, and its detail can be seen in Ref. [2,6].

K. Shibanuma, T. Utsunomiya / Finite Elements in Analysis and Design 45 (2009) 806 -- 816

5. Numerical evaluation for crack analyses

5.1. Mode I crack model in the infinite plate (Model 1) The numerical analyses in a crack problem are evaluated by using a model which assumes a finite domain near the crack tip in the infinite plate. This model is generally used in the evaluation for the problem of the blending elements in the past studies [8,9], because the asymptotic solution described in (49) in Section 4.2, which is generally used as enrichment in the crack analyses, is exact solution in this model and therefore we can evaluate the influence of only the blending elements on the numerical accuracy. An infinite plate including a straight crack with a unit fracture mode I is considered and depicted in Fig. 10(a). Then, let us define a local finite square domain A including the crack tip in the center. The domain A is smaller than the crack length 2a in infinite plate. The size of this analytical domain A is 2×2 as shown in Fig. 10(b). The domain A is divided into the elements of N×N. The mesh size is therefore h = 2/N. The Gauss quadrature is used in the numerical integration for stiffness matrix or J-integral. Except the evaluations on the convergence study using L2 norm, the integration order is 4 in the support of enrichments C , based on our preliminary research, and the order is 2 in other areas because the standard finite element approximation for the bi-linear elements is used. The boundary conditions on the boundary A of the domain A assume the asymptotic solution of displacement field near the crack

A y p2 p1

x

2

A

∞

1.0E-03

Standard XFEM (RC = 1.0h), slope=0.91 Standard XFEM (RC = 2.0h), slope=0.92 PUFEM-based XFEM (RC = 1.0h), slope=1.08

1.0E-04 1.0E-02

1.0E-01 Mesh size : h

Fig. 11. Convergence of relative L2 norm of Model 1.

tip in the fracture mode I (KI = 1 and KII = 0). The nodal displacements are therefore directly specified on A as  ⎫ ⎧  cos 2 1 − 2 + sin2 2 ⎬ ⎨ ux r 1 uI = = (57)   ⎭ G 2 ⎩ uy sin 2 2 − 2 − cos2 2 where G corresponds to the modules of elasticity in shear. The nodes belonging to the set J on the boundary A are defined as the nodes p1 and p2 in Fig. 10(b), in which there are not only nodal degrees of freedom uI for the FEM approximation but also bI for the Heaviside enrichment. In such cases, the nodal degrees of freedom bI is specified by the relative displacement between the crack surfaces, and uI is specified to satisfy the asymptotic solution, respectively.  bx 0 r 1 bI = = (58a) G 2 2 − 2 by  ⎫ ⎧  cos 2 1 − 2 + sin2 2 ⎬ ux r ⎨ 1 (58b) = uI =   − H(x)bI ⎭ G 2 ⎩ uy sin 2 2 − 2 − cos2 2 First of the evaluations for the numerical results in this section, we investigate the convergence study using L2 norm. A convergence study of approximation accuracy is made by using following relative L2 norm e2 : e2 =

uap − u L2 (A) u L2 (A)

 =

A

|uap − u|2 dA  2 A |u| dA

1/2 (59)

Let us assume the following asymptotic solution u as the exact solution:  ⎫ ⎧  cos 2 1 − 2 + sin2 2 ⎬ ux r ⎨ 1 u= = (60)   ⎭ G 2 ⎩ uy sin 2 2 − 2 − cos2 2

ΓA

2a

Relative L2 norm: e2

1.0E-02

In this section, we present the evaluations of numerical results in the linear fracture mechanics to reconfirm the usefulness of the PUFEM-based XFEM, including comparisons with that by using the standard XFEM. Following two types of models including cracks will be used; a Mode I crack model in an infinite plate (Model 1), and an edge crack model under shear in a finite plate (Model 2). Using the Mode I crack model in the infinite plate, only the influence of the problem of the blending elements can be evaluated. We investigate the convergence study using L2 norm of the interpolation error as a basic evaluation of numerical approximation. Then, we evaluate the stress field around the crack tip and the effect of the enrichment area on the accuracy of linear fracture mechanics parameters, which are not evaluated sufficiently in [9]. Using the edge crack model under shear in the finite plate, the application of the PUFEM-based XFEM to the edge crack model under shear in the finite plate is also investigated for the evaluation of the accuracy of linear fracture mechanics parameters. In both examples, the plane stain conditions are assumed and their materials are defined by Young's modules (E) of 100 and Poisson ratio ( ) of 0.3.

∞

813

2

Fig. 10. An Mode I crack model in an infinite plate (Model 1): (a) infinite crack body with fracture mode I and (b) numerical model for the domain A.

In the evaluation of relative L2 norm, the numerical integration order is 5 on the entire domain, taking into account the asymptotic solution √ including the singularity of r. The conditions of ranges of the enrichment related to the nodes of C are defined as RC = 1.0h in the PUFEM-based XFEM and RC = 1.0h, 2.0h in the standard XFEM, respectively. The numerical results in the convergence study of L2 norm are shown in Fig. 11. In the standard XFEM, it is shown that enlarging the ranges of the enrichment RC reduces the errors but it is less effective in

814

K. Shibanuma, T. Utsunomiya / Finite Elements in Analysis and Design 45 (2009) 806 -- 816

Standard XFEM (RC = 1.0h)

Standard XFEM (RC = 2.0h)

PUFEM-based XFEM (RC = 1.0h)

2.8 2.4 2.0 1.6 1.2 0.8 0.4 0.0

Standard XFEM (RC = 1.0h)

Standard XFEM (RC = 2.0h)

PUFEM-based XFEM (RC = 1.0h) 0.7 0.5 0.3 0.1 -0.1 -0.3 -0.5 -0.7

Fig. 12. Numerical results of stress fields of Model 1: (a) y-axial stress field y and (b) errors from the asymptotic solution.

the convergence rate. On the other hand, the PUFEM-based XFEM with RC = 1.0h provides better results than the standard XFEM with RC = 2.0h in the convergence rate as well as improved numerical accuracy, although both of them are equivalent in nodal degrees of freedom. Then, the stress field calculated by using the standard XFEM and the PUFEM-based XFEM is evaluated by comparing for the asymptotic solution, respectively. The numerical models in this evaluation have elements of 11×11. The radius RC using definition of the set of nodes C, which is related to crack tip enrichment, is the same as in Section 5.2 specified as RC = 1.0h in the PUFEM-based XFEM and RC = 1.0h, 2.0h in the standard XFEM. These conditions of the standard XFEM provides similar conditions to those of the PUFEM-based XFEM, i.e. the range of enrichment and the number of nodal degrees of freedom, respectively. Under the fracture mode I (KI = 1 and KII = 0), the y-axial stress field of the asymptotic solution y is

y = √

1

2 r

cos

2

 1 + sin

2

sin

3 2

 (61)

The numerical result of y-axial stress field y and their errors from the asymptotic solution are shown in Fig. 12(a) and (b), respectively. The area of blending elements is shown by red border lines. We can see large errors in the blending elements when the standard XFEM is used. On the other hand, it is shown that the PUFEMbased XFEM enables to improve the problem and to reproduce the stress field near the crack tip more exactly. Then, the numerical accuracy of the fracture mechanics parameter is evaluated. The parameter used for the evaluation is the energy release rate  G, which are the basic linear fracture mechanics parameter. For the evaluation of the energy release rate  G, the J-integral method is

employed. These methods are generally used in crack analyses by using the XFEM. The procedure of the J-integral method in the PUFEMbased XFEM is the same as that of the standard XFEM. Therefore, the detailed formulation of the method can be seen in Refs. [14,16]. The numerical models for the evaluation have the elements of 15×15. The conditions of ranges of the enrichment related to the nodes of C correspond to five cases of RC = 1.0h, 2.0h, 3.0h, 4.0h and 5.0h. In addition, the integral radiuses for the J-integral correspond to two cases of RF = 1.0h, 2.0h for the respective ranges RC . Numerical results for the evaluation of the energy release rate G˜ are shown in Fig. 13. In the cases of the standard XFEM, it is found that the accuracy of the energy release rate  G is lower in the condition where the two radiuses of RC and RF correspond to the same value, i.e. RC = RF = 1.0h or RC = RF = 2.0h. It means that this lack of numerical accuracies is caused when the integration path of J-integral is overlapped with the domain of the blending elements. It is therefore demonstrated that the influence of the blending elements, which is assumed to be small enough in many past studies, is too large to be neglected in the evaluation of the energy release rate  G by using standard XFEM. On the other hand, in the cases of the PUFEM-based XFEM, the energy release rate  G can be evaluated with significantly higher accuracy than those of the standard XFEM in all evaluations, even if the blending elements are overlapped with the integration path of J-integral. 5.2. Edge crack model under shear in a finite plate (Model 2) In this section, we evaluate the numerical accuracy of the linear fracture mechanics parameter for a finite plate model using the PUFEM-based XFEM. The parameter used for the evaluation is the energy release rate  G, as evaluated for the Mode I crack model in the infinite plate in the Section 5.1.

K. Shibanuma, T. Utsunomiya / Finite Elements in Analysis and Design 45 (2009) 806 -- 816

815

14.0 Standard XFEM (RF= 1.0h)

12.0

Standard XFEM (RF = 2.0h)

10.0

PUFEM-based XFEM (RF = 1.0 h)

Error (%)

8.0

PUFEM-based XFEM (RF = 2.0h)

6.0 4.0 2.0 0.0 -2.0 -4.0 0.0

1.0

2.0 3.0 4.0 Radius for Definition of C-nodes : RC (×h)

5.0

6.0

Fig. 13. Numerical results for the energy release rate of Model 1.

= 1

L = 16 a = 3.5

According to the above results, the numerical accuracy can be improved by enlarging the radius RC in the standard XFEM in order to avoid the overlap with the radius RF . However, such a large range of the enrichment around the crack tip causes a limitation in modeling the enrichment close to the boundary of the structural model, in particular for the crack propagation process. On the other hand, in the cases of the PUFEM-based XFEM, the error of the energy release rate  G can be significantly reduced, compared with those of the standard XFEM. In addition, the PUFEM-based XFEM shows good accuracy even if the radiuses RC and RF are similarly minimal as RC = RF = 1.0h. Therefore, the PUFEM-based XFEM is effective for the actual crack analyses. 6. Conclusion

W=7 Fig. 14. An edge crack model under shear in a finite plate (Model 2).

The numerical model in this section is an edge crack model, as shown in Fig. 14, it is fixed at its bottom and subjected to a uniform shear stress  = 1 applied at its top. The model has a length L = 16, width W = 7, and the crack length is a = 3.5. The reference energy release rate in this evaluation is  G = 10.71 [17]. The digitization of the numerical models is 47×107 elements. The conditions of ranges of the enrichment related to the nodes of C correspond to five cases of RC = 1.0h, 2.0h, 3.0h, 4.0h and 5.0h. In addition, the integral radiuses for the J-integral correspond to two cases of RF = 1.0h, 2.0h for the respective ranges RC . The numerical results are shown in Fig. 15. The results show the similar tendency to those obtained from the Mode I crack model in the infinite plate presented in Section 5.1, as follows. In the cases of the standard XFEM, it is found that the accuracy of the energy release rate  G is lower in the condition where the radius RF for J-integral is overlapped with the radius RC . That is under the conditions of RC = RF .

In the past studies, it is pointed out that the approximation of the standard XFEM has the problems caused by the blending elements lowering the convergence rate. On the other hand, the corrected XFEM proposed by Fries [9] showed the significant improvement of the accuracy reduction caused by the blending elements. The usefulness of this method is however based on the results of the inductive evaluations, because the theoretical background is not sufficiently described in the proposal. In the present paper, an approximation of the XFEM has been reformulated based on the concept of the PUFEM approximation to solve the problem of the blending elements. In the evaluation of the approximation accuracies for a one dimensional problem, it is shown that the insufficient convergence rate on the L2 norm in the standard XFEM has been overcome in the reformulated XFEM based on the PUFEM. It is therefore found out that this reformulated XFEM provides a potential ability to solve the problem of blending elements in the standard XFEM, by showing the relationship between the XFEM and the PUFEM. The form of approximation in the reformulated XFEM results in the coincidence with that of `the corrected XFEM' proposed by Fries [9]. The theoretical validation of the corrected XFEM is therefore found out based on the concept of the PUFEM approximation. That is, `the corrected XFEM' and `the reformulated XFEM in this paper' may be unified as a `PUFEM-based XFEM'. The PUFEM-based XFEM is applied to the two dimensional linear fracture mechanics to reconfirm the validation for actual use. As the

816

K. Shibanuma, T. Utsunomiya / Finite Elements in Analysis and Design 45 (2009) 806 -- 816

8.0 Standard XFEM (RF = 1.0)h

6.0

Standard XFEM (RF = 2.0h) PUFEM-based XFEM (RF = 1.0h)

4.0 Error (%)

PUFEM-based XFEM (RF = 2.0h)

2.0 0.0 -2.0 -4.0 -6.0 0.0

1.0

2.0 3.0 4.0 Radius for Definition of C-nodes : RC (×h)

5.0

6.0

Fig. 15. Numerical results for the energy release rate of Model 2.

result, the following features of the PUFEM-based XFEM are developed and compared with the standard XFEM. (1) In the convergence study using L2 norm of the interpolation error, we found the improvement on the convergence rate as well as the accuracy. (2) The stress field around the crack tip is exactly reproduced even in the blending elements. (3) In the evaluation of the linear fracture mechanics parameter, the good accuracy is provided even if both of the ranges of the enrichments and the integration path of J-integral are minimal. This leads the effectiveness for the actual crack analyses. Based on the aforementioned facts, it is possible to conclude that the problem of the blending elements has been sufficiently solved by the PUFEM-based XFEM proposed in this paper. References [1] J.M. Melenk, I. Babuska, The partition of unity finite element method: basic theory and applications, Computer Methods in Applied Mechanics and Engineering 39 (1996) 289–314. [2] T. Belytschko, T. Black, Elastic crack growth in finite elements with minimal remeshing, International Journal for Numerical Methods in Engineering 45 (1999) 602–620. [3] M. Duflot, A study of the representation of cracks with level sets, International Journal for Numerical Methods in Engineering 70 (2007) 1261–1302. [4] Q.Z. Xiao, B.L. Karihaloo, Direct evaluation of accurate coefficients of the linear elastic crack tip asymptotic field, Fatigue & Fracture of Engineering Materials & Structures 26 (2003) 719–729. [5] X.Y. Liu, Q.Z. Xiao, B.L. Karihaloo, XFEM for direct evaluation of mixed mode SIFs in homogeneous and bi-materials, International Journal for Numerical Methods in Engineering 59 (2004) 1103–1118. [6] K. Shibanuma, T. Utsunomiya, Curved-crack modeling by X-FEM and its application for simulation of crack growth. Doboku Gakkai Ronbunshuu A 63 (2007) 108–121, DOI:10.2208/jsceja.63.108.

[7] J. Chessa, H. Wang, T. Belytschko, On the construction of blending elements for local partition of unity enriched finite elements, International Journal for Numerical Methods in Engineering 57 (1999) 1015–1038. [8] R. Gracie, H. Wang, T. Belytschko, Blending in the extended finite element method by discontinuous Galerkin and assumed strain methods, International Journal for Numerical Methods in Engineering 74 (2008) 1645–1669. [9] T.P. Fries, A corrected XFEM approximation without problems in blending elements, International Journal for Numerical Methods in Engineering 75 (2008) 503–532. [10] S. Bordas, B. Moran, Enriched finite elements and level sets for damage tolerance assessment of complex structures, Engineering Fracture Mechanics 73 (2006) 1176–1201. [11] E. Wyart, D. Coulon, M. Duflot, T. Pardoen, J.F. Remacle, F. Lani, A substructured FE-shell/XFE-3D method for crack analysis in thin-walled structures, International Journal for Numerical Methods in Engineering 72 (2007) 757–779. [12] G. Ventura, R. Gracie, T. Belytschko, Fast integration and weight function blending in the extended finite element method. International Journal for Numerical Methods in Engineering 2007, DOI:10.1002/nme.2387 (early view). [13] J.N. Reddy, An Introduction to the Finite Element Method, second ed., McGrawHill, 1993. [14] N. Mo¨es, J. Dolbow, T. Belytschko, A finite element method for crack growth without remeshing, International Journal for Numerical Methods in Engineering 46 (1999) 131–150. [15] M. Fleming, Y.A. Chu, B. Moran, T. Belytschko, Enriched element-free Galerkin methods for crack tip fields, International Journal for Numerical Methods in Engineering 40 (1997) 1483–1504. [16] K. Shibanuma, T. Utsunomiya, Proposal on approximation of path-independent M-integral by mapping and analyses of kinked or curved crack using X-FEM, Doboku Gakkai Ronbunshuu A 64 (2008) 303–316 DOI:10.2208/jsceja.64.303. [17] J. Yau, S. Wang, H. Corten, A mixed-mode crack analysis of isotropic solids using conservation laws of elasticity, Journal of Applied Mechanics 47 (1980) 335–341.