Computers and Structures 87 (2009) 1015–1021
Contents lists available at ScienceDirect
Computers and Structures journal homepage: www.elsevier.com/locate/compstruc
Strict bounds for computed stress intensity factors J. Panetier a,*, P. Ladevèze a,b, F. Louf a a
Laboratoire de Mécanique et Technologie, ENS Cachan, CNRS, Université Paris 6, PRES Universud Paris, 61, Avenue du Président Wilson, 94230 Cachan, France EADS Foundation Chair ‘Advanced computational structural mechanics’, Laboratoire de Mécanique et Technologie, ENS Cachan, CNRS, Université Paris 6, PRES Universud Paris, 61, Avenue du Président Wilson, 94230 Cachan, France b
a r t i c l e
i n f o
Article history: Received 16 November 2006 Accepted 27 November 2008 Available online 12 January 2009 Keywords: Stress intensity factors Error estimation Strict bounds Extended finite element method
a b s t r a c t This paper addresses the subject of local error estimators for finite element methods involving local enrichment. More precisely, we focus on the extended finite element method, which was developed in the last decade for fracture mechanics. We present a method for the calculation of accurate strict bounds of mixed-mode stress intensity factors associated with an extraction operator and the theory of error in constitutive relation. This method has already been shown to be effective for standard finite elements. The objective of this paper is to show that it can easily be adapted to the extended finite element method. Ó 2008 Civil-Comp Ltd and Elsevier Ltd. All rights reserved.
1. Introduction The three main concepts developed in the last few decades for the estimation of the global quality of a numerical solution obtained by a finite element method (FEM) are the constitutive relation error [1], the error based on the residuals of the equilibrium equations [2] and the error based on smoothing techniques [3]. A review of these concepts can be found in [1]. They have proved to be effective. However, they lead only to a global error in the energy norm, which is often insufficient. Today, research is focusing on the control of quantities of interest (see Refs. [4–13]). In this paper, we deal with the estimation of the quality of calculated stress intensity factors (SIFs). The first important paper dealing with error estimation for fracture mechanics is [25]. It addresses the problem of the error on the path taken by a propagating crack, but it does not involve local error estimations. An a posteriori error estimation for SIFs was proposed in [9], but this technique does not lead to guaranteed bounds. In [10] and [11] the authors were interested in the J-integral. The J-integral is a quadratic functional of the displacement: thus, the bounds are based on linearization and are not strict. Guaranteed bounds for the J-integral were established in [12]. Since the extraction operator for the SIF is a linear functional of the displacement, it is possible to apply the goal-oriented error estimation technique developed in [4–8]. This technique is based on the use of the FEM to solve the adjoint problem. Thus, the error in the energy norm is calculated both for the original problem and for the adjoint problem. Using these errors, one can deduce strict
* Corresponding author. Tel.: +33 1 4740 5301; fax: +33 1 4740 2785. E-mail address:
[email protected] (J. Panetier).
upper and lower bounds. In [13] we demonstrated the accuracy of this method for a standard finite element analysis. The objective of the present study is to adapt the goal-oriented error estimation method to the extended finite element method (XFEM) presented in [14,15]. The XFEM, based on the Partition of Unity Method (PUM) introduced in [16], involves a local enrichment of the finite element basis functions in order to take into account both the discontinuity across the crack’s lips and the singularity in the vicinity of the crack’s tip. In this early work, only local enrichments involving the asymptotic basis functions were studied. We show that for this type of enrichment a specific stress recovery technique enables one to calculate strict and accurate bounds. The present paper is organized as follows: in Section 2, we present the reference problem to be solved. In Section 3, we briefly recall the basics of the constitutive relation error. Section 4 introduces the goal-oriented error estimation and establishes the formulation of the bounds to be calculated. The extraction operator for the SIF is presented briefly in Section 5. In Section 6, dedicated to the XFEM, we detail the specific method for the recovery of the equilibrated stress fields. Finally, Section 7 presents numerical examples which demonstrate the accuracy of the method.
2. The continuous reference problem Let us consider a 2D linear elastic structure which occupies a domain X bounded by @ X. The external loading is represented by a prescribed displacement ud over @ 1 X and a prescribed surface density force F d over @ 2 X ¼ @ X @ 1 X. The Hooke’s operator is denoted K. Thus, the problem to be solved can be formulated as follows: find a solution ðu; rÞ defined in X which verifies:
0045-7949/$ - see front matter Ó 2008 Civil-Comp Ltd and Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2008.11.014
1016
J. Panetier et al. / Computers and Structures 87 (2009) 1015–1021
with:
the kinematic constraints:
u2U
uj@ 1 X ¼ ud
ð1Þ
the equilibrium equations:
r2S 8u 2 U0 : Z Z Tr½reðu Þ dX þ X
F d u dS þ
@2 X
Z
fd u dX ¼ 0
ð2Þ
¼ ð3Þ
1 2
eðuÞij ¼ ðui;j þ uj;i Þ
ð4Þ
U is the space in which the displacement field is sought (U ¼ ½H1 ðXÞ3 ). S is the space of the equilibrated stress fields (S ¼ ½L2 ðXÞ6 ). U0 is the subspace of U with the prescribed displacement set to 0. In the following sections, the solution of the reference problem will be denoted ðuex ; rex Þ.
Let us consider the solution pair ðuh ; rh Þ calculated by a finite element method using a classical displacement approach. The approximate solution ðuh ; rh Þ is such that uh belongs to subspace Uh of U and rh to subspace Sh of S. The subspace Uh is associated with a finite element mesh of characteristic size h and with a set of finite element basis function. Sh is defined by the derivative of the finite element basis functions of Uh due to the constitutive relation given by Eq. (3). The principle of the constitutive relation error consists in associating a new displacement stress pair with the approximate pair ^h ; r ^h Þ which verifies the kinematic constraints (1) and the equiðu librium Eq. (2). In the classical displacement approach, uh verifies ^ h is taken as uh . However, the finite ele(1); therefore, usually, u ment stress field rh is equilibrated in the finite element sense, but does not verify (2). The general technique for building an equilibrated stress field from rh has been developed for different types of problems in the last few decades and is described in Refs. [17,18]. The key of the method is to recover an admissible stress ^h 2 S (see Eq. (2)).The principle confield denoted r ^h such that: r sists in writing a prolongation condition (5) between the finite element stress field and the equilibrated stress field. Thus, r ^h can be calculated explicitly at a very small cost
^h Þeðui Þ dE ¼ 0 Tr½ðrh r ð5Þ
^ h ÞÞK1 ð^ ^ h ÞÞ dE Tr½ð^ rh Keðu rh Keðu
ð6Þ
E
The global constitutive relation error over X is:
e2cre ¼
X
e2cre;E
ð7Þ
E
Therefore, the relative error is:
2cre ¼ R X
e2cre Tr½rh;m K1 rh;m dX
Tr½ðrex rh ÞK1 ðrex rh Þ dX
ð10Þ
In order to assess the quality of the error estimation, we use the effectivity index:
c¼
ecre eh
ð11Þ
In practice, the effectivity index must be close to 1 to be relevant and greater than 1 to be reliable. The Prager–Synge theorem [19] leads to the property:
^h k2r;X þ kuex uh k2u;X e2cre ¼ krex r
ð12Þ
with the energy norm:
Z ZX
k k2u;X ¼
Tr½K1 dX and Tr½eðÞKeðÞ dX
ð13Þ
X
Thus, this theorem guarantees that:
eh 6 ecre
ð14Þ
This means that for the error estimated using the constitutive relation error the effectivity index is always greater than 1. 4. Goal-oriented error estimation We will now consider the assessment of a linear quantity of interest I. Let us assume that L is a linear functional of the displacement such that:
I ¼ LðuÞ
ð15Þ
The calculated quantity of interest is:
Ih ¼ Lðuh Þ
ð16Þ
Because of the linearity assumption, we can express the error on I as:
jIex Ih j ¼ jLðuex uh Þj ¼ jLðeh Þj
ð17Þ
Classically, the functional L can be expressed in terms of a pree and a body force ~f : stress R
LðuÞ ¼
Z
X
Then, the quality of the solution pair is measured by the residual related to the verification of the constitutive relation. We define the constitutive relation error for each element E:
Z
Tr½eðuex uh ÞKeðuex uh Þ dX
X
Z
k k2r;X ¼
3. The constitutive relation error in elasticity
for each element E and for each basis function ui
Z
X
where eðuÞ denotes the strain associated with the small displacement
e2cre;E ¼
e2h ¼
X
r ¼ KeðuÞ
E
ð9Þ
The discretization error is defined as:
the constitutive relation:
Z
1 ^ h ÞÞ ð^ rh þ Keðu 2
^h;m ¼ r
e eðuÞ dX þ Tr½ R
Z
~f u dX
ð18Þ
X
The technique for calculating the goal-oriented error estimate is detailed in [4–8]. It involves the resolution of a new finite element problem, called the adjoint problem. This is a new elasticity problem defined in the same domain X, in which the prescribed displacement ud is set to zero and the loading consists of the e and the body force ~f . prestress R The reference adjoint problem is defined as:
~ 2 U0 : find u Z ~ Þeðu Þ dX ¼ Lðu Þ Tr½Keðu
8u 2 U0
ð19Þ
X
ð8Þ
Thus, we obtain an approximate solution of the adjoint problem de~ h Þ. noted ð~ rh ; u
1017
J. Panetier et al. / Computers and Structures 87 (2009) 1015–1021
From the solution of the original problem, we deduce the admis^h ; r ^h Þ and the associated constitutive relation error sible pair ðu Ecre ; From the solution of the adjoint problem, we deduce the admis^ ;r ^ ~ sible pair ðu h ~ h Þ and the associated constitutive relation error ~cre . E In our previous works, we used the upper bound described in [1]:
~cre jI Ih j 6 Ecre E
ð20Þ
An improvement to this bound, introduced in [20], leads to the new upper bound:
jI Ih Ihh j 6
1 ~cre Ecre E 2
ð21Þ
Proof
Iex Ih ¼
Z
e K1 ðrex rh Þ dX þ Tr½ R
X
Iex Ih ¼
Z
X
¼
Z
~f ðuex u ^ h Þ dX þ
X
Z
Z
e K1 ð^ Tr½ R rh rh Þ dX
X
Z
~f ðu ^ h uh Þ dX
X
e K1 ðrex Keðu ^ h ÞÞ dX þ Tr½ R
X
~f ðuex u Þ dX h
X
e K1 ðrex r Tr½ R ^h Þ dX þ
þ
Z
Z X
Z
e K1 ðrh Keððu ^ h Þ dX þ Tr½ R
X
~f ðuex u ^ h Þ dX Z
~f ðu u ^ h Þ dX h
X
Introducing the adjoint problem, we get:
Iex Ih ¼
Z
^~ ÞÞ dX ^~h Keðu ^ h ÞÞK1 ðr Tr½ðrex Keðu h X Z ^~ Þ dX ^ h ÞÞeðu þ Tr½ð^ rh Keðu h X Z Z e K1 ðrh Keðu ^ h Þ dX þ ~f ðuh u ^ h ÞdX Tr½ R X
X
Z
e K1 ðrh Keðu ^ h Þ dX þ Tr½ R
X
Z
Z K a ¼ La ðuÞ ¼ Tr½ðKeð/v a Þ /ra ÞeðuÞ dx x2 Z ðra gradð/ÞÞ u dx a ¼ I; II
X
^h;m k ¼ E2cre 2krex r Thus, we get:
where ra and v a are the singular analytical solutions of the reference problem (for further details, see [21,24,13]) and / is the following linear function of the polar coordinates, defined in x2 :
R2 r R2 R1
ð25Þ
1 ~cre Ecre E 2
In practice, our estimate of the error is:
^eh ¼ jI Ih j 6
X1 ~cre;E þ jIhh j Ecre;E E 2 E
ð22Þ
lo and we get an upper bound Iup h and a lower bound Ih :
^ Iup h ¼ Ih þ e h
and
^ Ilo h ¼ Ih e h
ð24Þ
x2
R1 and R2 denote respectively the radii of C 1 and C 2 in Fig. 2. The approximate values of K I;h and K II;h are calculated from the approximate finite element solution uh . The integration is carried out over a patch of elements E such that E \ x2 –ø. Then, following the technique described in Section 4, we can calculate bounds of each K a;h . The accuracy of the method for the standard finite element method was demonstrated in [13].
~f ðu u ^ h Þ dX h
and introduce Eq. (9), the Prager–Synge theorem leads to:
jI Ih Ihh j 6
the lips are traction-free and that there are no body forces f d . The problem to be solved is the same as in Section 1: find a pair ðu; rÞ which verifies Eqs. (1)–(3). The extraction operator for the SIFs K I and K II is based on a modification of the method proposed in [21]. It involves singular functions which are non-finite energy solutions of the reference problem. This technique is developed for several kinds of singularity in [24]. The calculation of the SIFs consists in an integral over an arbitrary crown x2 around the tip of the crack (see Fig. 2) and leads to K I and K II separately without introducing any errors except for those due to the calculation of u. The expression of the extractor operator is:
/ðr; hÞ ¼
If we denote:
Ihh ¼
Fig. 1. The reference problem.
ð23Þ
It is important to emphasize that this technique leads to strict bounds of the calculated quantity of interest. In practice, Ihh is very small, so the new error ^eh can be divided by two compared with the previous error. h
5. Application to the control of stress intensity factors in mixedmode problems Let us consider a two-dimensional structure X with an existing crack tip at O (see Fig. 1). For the sake of simplicity, we assume that
Fig. 2. Arbitrary crown.
1018
J. Panetier et al. / Computers and Structures 87 (2009) 1015–1021
Fig. 3. Subdivision of the problem on X.
Fig. 4. The models for the numerical examples.
1019
J. Panetier et al. / Computers and Structures 87 (2009) 1015–1021 Table 1 Problem 1: calculated K Ih and their bounds.
Table 3 Problem 2: calculated K IIh and their bounds.
Number of elements
100
400
1600
6400
Number of elements
100
400
1600
6400
K Ih ^ eh
5.1904 0.9126
4.9724 0.1774
5.0620 0.0538
5.0774 0.0197
K IIh ^ eh
1.8117 1.2067
1.8369 0.2342
1.8602 0.0626
1.8652 0.0175
K up Ih K lo Ih
6.1030
5.1498
5.1158
5.0971
K up IIh
3.0184
2.0711
1.9228
1.8827
4.2778 5.0841
4.7950 5.0841
5.0082 5.0841
5.0577 5.0841
K lo IIh Reference
0.605 1.8653
1.6027 1.8653
1.7976 1.8653
1.8477 1.8653
6. Error estimation for the extended finite element method 6.1. The extended finite element method The principle of the XFEM, presented in [14,15], is to introduce some a priori knowledge of the exact solution near the crack’s tip to avoid having to remesh at each step of the propagation. This knowledge takes the form of a local enrichment based on the Partition of Unity Method (PUM) [16]. The position of the crack is described, as in [22], by segments with level sets. Two types of enrichment functions are usually introduced: For each node i 2 Ncrack , a Heavyside function, denoted H, which takes into account the discontinuity of the displacement along the crack. Only the nodes of the elements which are cut by the crack are enriched by function H. For each node i 2 N tip , a set of functions, denoted F j , which span the space of the asymptotic solutions in the vicinity of the crack’s tip. These functions F j enrich the nodes of the elements enclosing the tip. In our approach, only the second type is used. The displacement (in polar coordinates) belongs to the space spanned by the following functions:
F j ðxÞ ¼
pffiffiffi pffiffiffi h pffiffiffi h pffiffiffi h h ; r sin ; r sin sinðhÞ; rcos sinðhÞ r cos 2 2 2 2 ð26Þ
The approximate displacement field uh within Domain X takes into account the enrichment with the PUM:
uh ðxÞ ¼
X
ui ðxÞui þ
i2N
X
ui ðxÞðF j ðxÞbij Þ þ
i2N tip
X
ui ðxÞHðxÞai
ð27Þ
i2Ncrack
where function ui is the classical finite element shape function P1 in our study. It was proved in [15] that unless the convergence rates for the standard and enriched problems are the same (see [23]) enrichment has a decisive influence on the calculation of K I and K II . The quality of these quantities improves to a great extent when they are calculated with the XFEM method, but it still remains important to assess this quality seriously. A few works are devoted to error estimation for the XFEM. In [26], the problem of a guaranteed error estimator is addressed in the general case of the GFEM (generalized finite element method). The authors give an analysis of the quality of several error estimators based on patch residuals
indicators. The principle and the details of our method are the same as in the previous sections. Since field uh still verifies the ^ h ¼ uh . The problem consists kinematic constraints (1), we choose u ^h , which requires a special treatment. in recovering the field r 6.2. Stress field recovery in the XFEM For the sake of simplicity, these preliminary studies were carried out with a meshed crack, so functions F j alone were kept in the enrichment. Only the nodes whose support contains the crack’s tip are enriched. The specific recovery of r ^h is achieved by dividing the XFEM problem into two problems. The structure X is divided into: a sub-structure X1 which contains all the enriched elements; a sub-structure X2 which contains all the elements with classical approximation functions. Let Cint denote the internal segments of X which constitute the common border between X1 and X2 . Fig. 3 shows how structure X is separated into sub-structures X1 and X2 along with mesh examples. On the submesh of X2 , square marks designate nodes enriched by functions F j and cross marks designate the integration points. Sixteen-point Gauss integration was used for the enriched elements. 6.2.1. First step: recovery over X1 X1 contains the enriched elements near the crack’s tip. The ^1h which verifies Eq. (2). The solution conproblem is: find a field r sists in using the well-known asymptotic functions solutions of the reference problem. Usually, only the first-order development is taken into account (b1 ¼ 12). Here, we will develop the asymptotic solution to a higher order. The general finite energy solution is, in polar coordinates:
30 problem 1 problem 2
25
Effectivity index
Reference
20
15
10
Table 2 Problem 2: calculated K Ih and their bounds. Number of elements
100
400
1600
6400
K Ih ^ eh
8.1120 1.8887
7.8538 0.3585
7.9978 0.1019
8.0346 0.0352
K up Ih
10.0007
8.2123
8.0997
8.0698
6.2233 8.0455
7.4953 8.0455
7.8959 8.0455
7.9994 8.0455
K lo Ih Reference
5
0 5
10
15
20
25
1/h Fig. 5. Effectivity indexes for K I .
30
35
40
1020
rrr ¼
J. Panetier et al. / Computers and Structures 87 (2009) 1015–1021
1 1 U;hh þ U;r r2 r
ð28Þ
rhh ¼ U;rr 1 rrh ¼ U;h r ;r
ð29Þ ð30Þ
with the Airy function:
Uðr; hÞ ¼
n X
rbi þ2 Wi ðhÞ
ð31Þ
i¼1
and:
Wi ðhÞ ¼ Ai sinðbi hÞ þ Bi cosðbi hÞ þ C i sin½ðbi þ 2Þh þ Di cos½ðbi þ 2Þh: ð32Þ
where bi 2 12 ; 0; 12 ; 1; 32 ; 2; . . . and:
for bi ¼ n;
n2N
8. Conclusion
B i þ Di ¼ 0 Ai n þ C i ðn þ 2Þ ¼ 0 1 for bi ¼ n þ ; n 2 N 2 Ai þ C i ¼ 0 1 5 Bi n þ þ Di n þ ¼0 2 2 The coefficients (A; B; C and D) are calculated in order to minimize the error:
^1h k ¼ Eh ðA; B; C; DÞ ¼ krh r
Z X1
^1h ÞK1 ðrh r ^1h ÞdX Tr½ðrh r
ð33Þ
A fourth-order development is found to be sufficient to calculate an accurate stress field r ^1h . 6.2.2. Second step: recovery over X2 In order to ensure continuity of the normal stress between X1 and X2 over Cint , one must solve a new standard finite element problem in which the external load is the same as that defined on X and the load on Cint is given by:
^1h nj@X1 ¼ 0 on Cint r nj@ X2 þ r
ð34Þ
^2h is Then, using the classical recovery method (see [1,17,18]), Field r calculated from the new solution on X2 . 6.2.3. Third step: recovery over X ^1h and r ^2h have been calculated, Field r ^h is given simply by: Once r
^ 1h r ^2h ^h ¼ r r
The reference stress intensity factors were obtained with a refined mesh of 20,000 elements. These tables illustrate the accuracy of the calculated bounds for K I as well as K II . Fig. 5 shows the effectivity indexes for K I (Eq. (11)). One can see that the effectivity index c is close to 1 and does not change much with mesh refinement. However, one can note that for the coarse mesh c is really high: when the mesh is too coarse, the finite element basis functions represent the loading due the asymptotic ^1h on Cint rather poorly. This problem could easily be fixed stress r by refining Contour Cint . However, we never use so much coarse meshes for that kind of problems. Plotting the effectivity indexes for K II in the case of Problem 2 would be meaningless. The error between K IIh and K II for the 400-element mesh is still about one percent. It is interesting to compare these results with those obtained with no local enrichment in [13]. The bounds and the effectivity indexes are much better with the singularity functions of the XFEM.
ð35Þ
7. Numerical examples Two examples of a single-edge crack in a plate will be presented. In the first problem (Fig. 4a), the loading consists in a uniform unit tension. In the second problem (Fig. 4b), the loading consists in a unit end shear. The geometry is the same for both problems: a ¼ 1; L ¼ 2; w ¼ 2. The Young’s modulus is E ¼ 210; 000 MPa, and the Poisson’s ratio is m ¼ 0:3. All the calculations were carried out with four-node rectangular elements and included local enrichment with the F j functions described in Section 6. We studied the evolutions of the calculated stress intensity factors and their bounds as functions of the mesh size. The coarse mesh is a 100 elements mesh which is a very coarse mesh for a plate containing a crack. The fine mesh is a 6400 elements mesh.Table 1 gives the results of the first problem for K I . Tables 2 and 3 give the results of Problem 2 for K I and K II .
We presented a technique for calculating strict bounds of the stress intensity factors in an extended finite element analysis. We showed that we can easily adapt the previous goal-oriented error estimation methods to the XFEM. The main difficulty resides in the construction of equilibrated stress fields from the enriched solution. The technique we propose consists in separating the enriched elements from the standard elements. We used asymptotic solutions near the crack’s tip to get equilibrated fields. The numerical examples demonstrate the accuracy of the method for mixedmode problems. A companion article will present the admissible stress field recovery method for the discontinuous enrichment. References [1] Ladevèze P, Pelle JP. Mastering calculations in linear and non linear mechanics. Mechanical engineering series. Springer; 2005. [2] Babuska I, Strouboulis T. The finite element method and its reliability. Oxford: Oxford University Press; 2001. [3] Zienckiewicz OC, Zhu JZ. A simple error estimator and adaptative procedure for practical engineering analysis. Int J Numer Meth Eng 1987;24:337–57. [4] Rannacher R, Stuttmeier FT. A feedback approach to error control in finite element methods: application to linear elasticity. Comput Mech 1997;19:434–46. [5] Rannacher R, Stuttmeier FT. A posteriori error control and mesh adaptation for f. e. models in elasticity and elastoplasticity. In: Ladevèze P, Oden JT, editors. Advances in computational methods in mechanics. Amsterdam: Elsevier; 1998. [6] Péraire J, Patera A. Bounds for linear-functional outputs of coercive partial differential equations local error indicators and adaptive refinement. In: Ladevèze P, Oden JT, editors. Advances in adaptive computational methods. Elsevier; 1998. [7] Ladevèze P, Rougeot Ph, Blanchard P, Moreau JP. Local error estimators for finite element linear analysis. Comput Meth Appl Mech Eng 1999;176:231–46. [8] Prudhomme S, Oden JT. On goal-oriented error estimation for elliptic problems: application to the control of pointwise error. Comput Meth Appl Mech Eng 1999;176:313–31. [9] Strouboulis T, Babuska I, Datta DK, Copps K, Gangaraj SK. A posteriori estimation and adaptive control of the error in the quantity of interest. Part I: A posteriori estimation of the error in the Von Mises stress and the stress intensity factors. Comput Meth Appl Mech Eng 2000;181:261–94. [10] Heintz P, Samuelson K. On adaptative strategies and error control in fracture mechanics. Comput Struct 2004;82:485–97. [11] Rüter M, Stein E. Goal oriented a posteriori error estimates in linear elastic fracture mechanics. Comput Meth Appl Mech Eng 2006;195:251–78. [12] Xuan ZC, Parès N, Peraire J. Computing upper and lower bounds for the Jintegral in two dimensional linear elasticity. Comput Meth Appl Mech Eng 2006;195:430–43. [13] Gallimard L, Panetier J. Error estimation of stress intensity factors for mixedmode crack. Int J Numer Meth Eng 2006;68:299–316. [14] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. Int J Numer Meth Eng 1999;45:601–20. [15] Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. Int J Numer Meth 1999;46:131–50. [16] Melenk JM, Babuska I. The partition of unity finite element method: basic theory and applications. Comput Meth Appl Mech Eng 1996;139:289–314. [17] Ladevèze P, Pelle JP, Rougeot Ph. Error estimation and mesh optimization for classical finite elements. Eng Comput 1991;8:69–80.
J. Panetier et al. / Computers and Structures 87 (2009) 1015–1021 [18] Ladevèze P, Rougeot Ph. New advances on a posteriori error on constitutive relation in f.e. analysis. Comput Meth Appl Mech Eng 1997;150:239–49. [19] Prager W, Synge JL. Approximation in elasticity based on the concept of functions space. Quart Appl Math 1947;5:261–9. [20] Ladevèze P. Upper error bounds of calculated outputs of interest for linear and non linear structural problems. CR Méc 2006:334. [21] Stern M, Becker EB, Dunham RS. A contour integral computation of mixedmode stress intensity factors. Int J Fract 1976;12:359–68. [22] Osher S, Sethian SA. Fronts propagating with curvature dependent speed: algorithms based on Hamilton–Jacobi formulations. J Comput Phys 1988;79: 12–49.
1021
[23] Laborde P, Pommier J, Renard Y, Salaün M. High-order extended finite element method for crack domains. Int J Numer Meth Eng 2005;24:354–81. [24] Babuska I, Miler A. The post-processing approach in the finite element method. Part 2. The calculation of stress intensity factors. Int J Numer Meth Eng 1984;20:1111–29. [25] Stone TJ, Babuska I. A numerical method with a posteriori error estimation for determining the path taken bay a propagating crack. Comput Meth Appl Mech Eng 1998;160:245–71. [26] Strouboulis T, Lin Z, Delin W, Babuska I. A posteriori error estimation for generalized finite element methods. Comput Meth Appl Mech Eng 2006;195: 852–79.