A new M-integral parameter for mixed-mode crack growth in orthotropic viscoelastic material

Engineering Fracture Mechanics 75 (2008) 4450–4465

Contents lists available at ScienceDirect

Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech

A new M-integral parameter for mixed-mode crack growth in orthotropic viscoelastic material R. Moutou Pitti, F. Dubois *, C. Petit, N. Sauvat, O. Pop Laboratory of Mechanics and Modeling of Materials and Structures in Civil Engineering (3MSCE), University of Limoges, Centre Universitaire de Génie Civil, Boulevard Jacques Derche, 19300 Egletons, France

a r t i c l e

i n f o

Article history: Received 7 September 2007 Received in revised form 23 April 2008 Accepted 29 April 2008 Available online 16 May 2008 Keywords: Creep crack growth Finite element analysis Mixed-mode fracture Viscoelastic media Orthotropic materials

a b s t r a c t This paper deals with a new independent path integral which provides the mixed-mode during a creep crack growth process in viscoelastic orthotropic media. The developments are based on an energetic approach using conservative laws. The mixed-mode fracture separation is introduced according to the generalization of the virtual work principle. The fracture algorithm is implemented in a finite element software and coupled with an incremental viscoelastic formulation and an automatic crack growth simulation. This Mintegral provides the computation of stress intensity factors and energy release rate for each fracture mode. A numerical validation, in terms of energy release rate and stress intensity factors, is carried out on a CTS specimen under mixed-mode loading for different crack growth speeds. Ó 2008 Elsevier Ltd. All rights reserved.

1. Introduction The widespread use of viscoelastic and anisotropic polymers for structural designs implies having a better understanding of their fracture mechanical behavior. Many defects, resulting from the fabrication process or other natural fact such as corrosion, are often the cause of initial cracks in materials. Also, crack propagation in practical engineering problems is often of mixed-mode type, and may have been caused by fatigue [1], overload or creep loading [2]. A way to understand this is to study the roles and the combinations of various modes in the fracture process. Numerous authors have proposed new criterions for mixed-mode fracture in anisotropic elastic [3] and in ductile solids [4,5]. But, for the polymers characterized by a viscoelastic behavior and isotropic or anisotropic symmetries [6], current methods do not allow mixed-mode separation during a time dependant crack growth process. The crack growth initiation phase is generally studied when operating an energetic balance during a unit crack tip advance. The critical time which induces the crack growth process is based on Griffith’s criterion [7]. According to a finite element approach, the main difficulty is to characterize the mechanical fields in the crack tip vicinity. In this context, numerous authors [8,9], have proposed non linear viscoelastic behavior in the process zone. With this theory, Dugdale [10] has performed a model in which the process zone size a is given by the constant tensile stress assumption. In the same way, Zhang [11] has shown the non-singularity of stress due to the existence of the cohesive zone. Dubois [12] has used this approach for the simulation of the crack growth process in wood. The high non-linearity of the time dependent behavior, characterizing this process zone, invites authors to develop non dependent path integrals. The main difficulties consist in the separation of dissipative energies caused by viscous properties and the crack tip advance. Dubois [12] has proposed the new integral

* Corresponding author. E-mail address: [email protected] (F. Dubois). 0013-7944/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.engfracmech.2008.04.021

R. Moutou Pitti et al. / Engineering Fracture Mechanics 75 (2008) 4450–4465

Nomenclature a, Da crack length and crack increment A1A2,B1B2 crack lips A(C1),A(C2) areas enclosed by C1 and C2 Aa,Ba(a 2 {1 . . . 7}) solicitation points Jacobian and transposed Jacobian matrix B,BT b width specimen C p1 ; C p2 elastic compliances in opening and shear mode in pth spring Ca(a 2 {1,2}) reduced viscoelastic compliances in a mode unit elastic compliance tensor C0 E(t) tangent modulus EX,EY,GXY longitudinal, transverse and shear elastic modulus F Helmholtz strain energy density F bilinear form of strain energy density F(tn) external load {Fp}(tn1) supplementary viscous load vector Gpv ; 2 Gpv Gv, 2Gv I1,I2 Jijkl J u p v p KI ; KI 1

1

total energy release rate in opening and shear mode storage in the pth spring energy release rate in opening and shear mode unvariable notations four order creep compliance tensor Rice’s integral real and virtual stress intensity factors in opening mode induced by pth spring

K pII ; v K pII real and virtual stress intensity factors in shear mode induced by pth spring K ra ða 2 f1; 2gÞ stress intensity factors in a mode u e K a ða 2 f1; 2gÞ opening intensity factors in amode K pT tangent matrix u u

p

kijkl

spring rigidity components

M Mh Mhv ~ n ~ n2 n1 ; ~ N spj

new M-integral modeling form of M-integral viscoelastic form of Mh-integral normal vector of components nj normal vectors to contours C1 and C2 total number of Kelvin Voigt cells roots of characteristic equation

Sp11 ; Sp22 ; Sp12 ; Sp33 ; compliance tensor components S surface domain t,s,b time variables u,v real and virtual displacement fields of components ui and vi up,vp real and virtual displacement fields in the pth spring elastic strain energy Ue V volume domain p spring of pth Kelvin Voigt cell Wvis viscous energy fracture energy WS Wi geometry of crack length ai x1,x2 axis a(x1,x2,t) spatial and temporal function b solicitation angles dL Lagrangian variation ~ ; dv ~ du real and virtual Euleurian representation du*,dv* real and virtual Lagrangian representation f~ep gðt n1 Þ strain vector e e ðt n1 Þ influence of global and local mechanical past history ~eij ðt n1 Þ; K a

gm ijkl

dash-pot viscositie components

~ h

vector field elastic stress and strain tensor component total real and virtual stress components

rij,eij

ruij ; rvij rpij ðuÞ; rpij ðvÞ real and virtual stress components in the pth spring

4451

4452

R. Moutou Pitti et al. / Engineering Fracture Mechanics 75 (2008) 4450–4465

t Poisson’s ratio oV contour around the crack tip Dt time increment Deij(tn),Drkl (tn) strain and stress increments {Dup}(tn) displacement vector increment fDF pext gðt n Þ increment nodal load vector

DK ea ðtn Þ; Drra ðtn Þ opening and stress intensity factor increments in a mode C1,C2 surface contours Pijkl strain components X integration domain Wijkl,Wa equivalent and reduced compliance matrix

parameter entitled Ghv which computes the energy release rate in the time domain for a viscoelastic and orthotropic media. Actually, his algorithm is limited to crack growth initiation and to opening mode configurations. In order to complete fracture tools, this paper deals with a new integral parameter allowing the fracture mode separation during the creep crack growth process. Firstly, the development of the M-integral concept, based on a conservative law approach, is presented. A pseudo potential, combining the real and auxiliary displacement fields, is employed in order to isolate different fracture modes by employing a generalization of the virtual work principle. The Noether’s theorem is used to define the general form of the M-integral. The generalization for a viscoelastic behavior is proposed in the second section. The global resolution is based on a complex algorithm which computes step by step the viscoelastic response, the separation of the energy release rate and the crack growth process by remeshing. According to finite element method, the resolution in the time domain is given by this algorithm. A validation, based on a CTS specimen, is presented in the last section. Also, during the crack growth process, the non-dependence domain is proven. Results are presented in terms of a parametric analysis of the energy release rate and stress intensity factor evaluations versus time, during the crack tip advance, for different mixed-mode configurations. 2. Integral parameters In a quasi-static case, if we consider a fractured viscoelastic media, the energy W, given by the external efforts, is the distribution of the elastic strain energy Ue, the viscous energy Wvis, and the energy Ws dissipated during a crack process such as

W ¼ U e þ W vis þ W s

ð1Þ

The elastic strain energy is the integration of the Helmoltz strain energy density noted F. For an elastic case, this free energy density is employed in the J-integral form [13]:

J¼

Z C1



F  n1  rij  nj 

 oui dC ox1

ð2Þ

C1 is an arbitrary curve, oriented by the normal vector ~ n characterized by it components nj, surrounding the crack tip in the external failure zone. ui and rij ((i,j) 2 {1,2}2) designate the displacement and stress components, respectively. The integral form supposes that the crack is oriented in the x1 direction as shown in Fig. 1.

x2

Failure zone

n

x1

a Γ1

Fig. 1. Path integration.

A(Γ1)

4453

R. Moutou Pitti et al. / Engineering Fracture Mechanics 75 (2008) 4450–4465

For a viscoelastic behavior, the Helmoltz strain energy density is defined by Stavermann [14] as follow:

F¼

1  2

Z 0

t

Z



t

0



½2  J ijkl ðt  sÞ  J ijkl ð2  t  s  bÞ

orij orkl dsdb os ob

ð3Þ

Jijkl are the components of the four order creep compliance tensor. t, s and b are three integration time variables. In order to develop an integral which enables us to separate fracture modes during a crack growth process in viscoelastic and orthotropic media, Noether’s theorem [15] is employed. This approach has been used by Chen [16] and Attigui [17] for elastic and dynamic problems, respectively. A generalization of this method to study a creep crack growth in viscoelastic material is proposed. Let us consider a perturbation of the virtual cinematically admissible displacement du. The transformed Noether’s theorem [15], based on a Lagrangian stationary, can be written as

dL ¼

Z Z t

dFdVdt ¼ 0

ð4Þ

V

In order to separate mixed-modes, Chen [16] has proposed a bilinear form of the elastic strain energy:

Fðu; vÞ ¼ F 

ð5Þ

In which u and v are real and virtual singular displacement fields, respectively. In an Arbitrary Lagrangian Eulerian config~i , are given by ˜ i and dv uration (ALE) [2–17], the variations of displacement fields, denoted d u

~i ¼ dvi ; dv ~i ¼ 0; dv ~i ¼ dvi þ dvi dv ~ i ¼ dui ; du ~ i ¼ 0; du ~ i ¼ dui þ dui du

ð6Þ

Considering the relation (6), the expression (4) of Lagrangian can be rewritten as

dL ¼

Z Z  V

t

 oF  oF  oF  oF   dui;a þ dui;a þ dvi;a þ dvi;a dVdt oui;a oui;a ovi;a ovi;a

ð7Þ

in which a(x1,x2,t) is a spatial and temporal function. If the virtual displacement is considered along the x1 axis, the following equations can be established for a crack tip advanced da

dui;a ¼

oui;a ovi;a da and dvi;a ¼ da ox1 ox1

ð8Þ

Using relation (8), Eq. (7) can be written as

dL ¼

 Z Z   oF oF  oF   da þ dui;a þ dvi;a dVdt ox1 oui;a ovi;a V t

ð9Þ

Considering (5) and (6), the last terms of Eq. (9) become:

     oF  oF oF  dui;a ¼  dui   dui oui;a oui;a oui;a ;a ;a      oF  oF oF  dvi;a ¼  dvi   dvi ovi;a ovi;a ovi;a ;a ;a

ð10Þ

In this case, the Lagrangian expression (9) can also take the following form:

! !         Z Z   oF  oF oF oF oF     dL ¼ da þ du   dui dVdt þ dv   dvi dVdt ox1 oui;j i ;j oui;j ;j ovi;j i ;j ovi;j ;j V t V t ! !       Z Z   Z Z   oF oF oF oF þ dui   dui dVdt þ dvi   dvi dVdt ou ou ov ov i;t i;j ;t i;t i;j ;t V t V t ;t ;t Z Z

ð11Þ

Using an Eulerian’s description, the real and virtual displacement variations (dui and dvi Þ are given by

dui ¼ 

oui ovi da and dvi ¼  da ox1 ox1

ð12Þ

Combining relations (11) and (12), the following expression is obtained:

      ! Z Z oF  oF oF dVdtda  dL ¼  dui;1 þ dvi;1 ox1 oui;j ovi;j V t V t ;j ;j !   Z Z   oF oF þ  duþi;1  dvi;1 dVdtda ¼ 0 oui;a ;a ovi;a ;a V t Z Z



oF  dui;1 oui;t





oF  þ dvi;1 ovi;t ;t

 !

dVdtda

;t

ð13Þ

4454

R. Moutou Pitti et al. / Engineering Fracture Mechanics 75 (2008) 4450–4465

The Eq. (10) enables us to write:



    oF  oF  ui;1 ¼  ui;1  F ;1 ðuÞ oui;a ;a oui;a ;a      oF oF  vi;1 ¼  vi;1  F ;1 ðvÞ ovi;a ;a ovi;a ;a

ð14Þ

Assuming Eq. (14), the last integral of Eq. (13) can be transformed into

      ! oF  oF oF dVdtda dL ¼   ui;1   vi;1 ox1 oui;j ovi;j V t ;j ;j !     Z Z   oF oF    dui;1  F ;1 ðuÞ þ dvi;1  F ;1 ðvÞ dVdtda ¼ 0 oui;a ovi;a V t ;a ;a Z Z

ð15Þ

Applying the Gauss–Ostrogradski theorem to the first term of Eq. (15), it becomes:

dL ¼

     Z Z  oF oF  F   n1   ui;1 þ  vi;1  nj dS dtda oui;j ovi;j t oV ! #     Z "Z     oF oF dui;1 þ dvi;1  F ;1 ðuÞ þ F ;1 ðvÞ dV dtda ¼ 0 þ oui;a ovi;a t V ;a ;a

ð16Þ

oV designates the boundary curve of V, see Fig. 2. These integration domains V and oV are supposed as constants versus time. According to the Lagrangian stationary for all variations of crack length da and time dt, relation (16) can be simplified as follow:

dL ¼

Z  oV



oF  oF  F  n1   ui;1 þ  vi;1 oui;j ovi;j 



 Z  nj dS þ



oF  dui;1 oui;a

V





oF  þ dvi;1 ovi;a ;a

 



F ;1 ðuÞ



!

þ

F ;1 ðvÞ

dV ¼ 0

;a

ð17Þ In order to simplify the expression (17), the following notations are introduced:

 oF  oF   ui;1 þ  vi;1 oui;j ovi;j !         oF oF dui;1 þ dvi;1  F ;1 ðuÞ þ F ;1 ðvÞ with I1 ðu; vÞ þ I2 ðu; vÞ ¼ 0 oui;a ovi;a ;a ;a

I1 ðu; vÞ ¼ F   n1  I2 ðu; vÞ ¼



ð18Þ

For plane problems, the integration volume, introduced by a crown around the crack tip with a constant width b, can be written as follow:

dV ¼ b  dS

ð19Þ

oV becomes a closed contour around the crack tip, composed by two contours (C1 and C2), and two segments (A1A2 and B2B1), Fig. 2. With these considerations, oV and V can be replaced by

x2

Γ1 A1 B1

n1

Γ2

A2

n2

B2

V Ω Fig. 2. Integration domain.

x1

4455

R. Moutou Pitti et al. / Engineering Fracture Mechanics 75 (2008) 4450–4465

oV ¼ C1 þ A1 A2 þ C2 þ B2 B1

ð20Þ

V ¼ AðC1 Þ  AðC2 Þ In which A(C1) and A(C2) are the areas enclosed by C1 and C2, respectively. Taking into account notations (20), the expression (18) can be rewritten as

Z

I1 ðu; vÞdC þ

Z

I2 ðu; vÞdV ¼

Z

V

oV

I1 ðu; vÞdC þ

Z

C1 þA1 A2 þC2 þB2 B1

I2 ðu; vÞ dV ¼ 0

ð21Þ

AðC1 ÞAðC2 Þ

Decomposing integration path, see expressions (20), the relation (21) can take the following form:

Z

I1 ðu; vÞdC þ

Z

I2 ðu; vÞdV þ

AðC1 Þ

C1

Z

I1 ðu; vÞdC þ

A1 A2

Z

I1 ðu; vÞdC ¼

Z

B2 B1

I1 ðu; vÞdC þ

Z

I2 ðu; vÞdV

ð22Þ

AðC2 Þ

C2

The non-dependent domain is justified by the relative orientation of normal vectors ~ n1 and ~ n2 , see Fig. 2. In these conditions, let us replace contour C2 by

C02 ¼ C2 with AðC2 Þ ¼ AðC02 Þ

ð23Þ

Eq. (22) can be rewritten as follow:

Z

I1 ðu; vÞdC þ

Z

I2 ðu; vÞdV þ

AðC1 Þ

C1

Z

I1 ðu; vÞdC ¼

Z C02

A1 A2 þB2 B1

I1 ðu; vÞdC þ

Z AðC02 Þ

I2 ðu; vÞdV

ð24Þ

If we suppose that the crack lip surfaces are not loaded, the contributions of A1A2 and B1B2 are zero. Introducing the M-notations, Eq. (24) shows us the path dependence property:

M¼

Z

I1 ðu; vÞdC þ

Z

I2 ðu; vÞdV ¼

Z

AðC1 Þ

C1

C02

I1 ðu; vÞdC þ

Z AðC02 Þ

I2 ðu; vÞdV

ð25Þ

Considering Eqs. (25) and (18), the M-integral form can take the following form:

M¼



Z

F   n1 

C1



oF  oF   ui;1 þ  vi;1 oui;j ovi;j



 Z  nj dC1 þ

 AðC1 Þ

oF  dui;1 oui;a



 þ

;a

oF  dvi;1 ovi;a

!    F ;1 ðuÞ þ F ;1 ðvÞ dV

 ;a

ð26Þ In a linear elastic configuration, the bilinear elastic strain energy density is defined such as

oF  1 v ¼ r oui;j 2 ij

and

oF  1 u ¼ r ovi;j 2 ij

ð27Þ

ðuÞ rðvÞ ij and rij are virtual and real stresses respectively. Eq. (27) allows transforming (26):

M¼



Z

C1

þ

Z

1 2



1 2

rvij  ui;j  n1  ð rvij  ui;j þ ruij  vi;j Þ  nj dC1

  1 v 1 1 rij  ðeuij Þ;1 þ ruij  ðevij Þ;1  ððrvij  euij Þ;1 þ ðruij  evij Þ;1 Þ dV 2 2 AðC1 Þ 2

Introducing the static admissible assumption for

M¼

1 2

Z



C1



rvij;1  ui  ruij  vi;j  nj dC1 þ

1 2

Z

ð28Þ

ðuÞ rðvÞ ij and rij , the M-integral form (28) can be simplified as

AðC1 Þ





rvij  ðui;j Þ;1 þ ruij  ðvi;j Þ;1 



rvij  ui;j

 ;1

þ ðruij  vi;j Þ;1

 dV

ð29Þ

The integral (29) is defined with a curvilinear integration domain C1. In order to implement this integral in a finite element software, it is easier to take into account a surface domain integral. In this context, the curvilinear domain must be transformed by introducing a vector field ~ h. This mapping function is continuously differentiable and takes these values: h1 = 1 and h2 = 0 inside the ring, ~ h¼~ 0 outside it, Fig. 3. Let us introduce the following definition:

p_ j;k ¼

 1 v rij;k  ui  ruij  vi;k 2

ð30Þ

Considering the Eq. (30) in (29), the following form ofMh is obtained:

Mh ¼

Z C1

p_ j;k  nj  hk dC1 þ

1 2

Z AðC1 Þ

ððrvij  ðui;j Þ;k þ ruij  ðvi;j Þ;k Þ  ððrvij  ui;j Þ;k þ ðruij  vi;j Þ;k ÞÞ  hk dV

ð31Þ

4456

R. Moutou Pitti et al. / Engineering Fracture Mechanics 75 (2008) 4450–4465

x2

⎛1⎞

θ = ⎜⎜ ⎟⎟ ⎝ 0⎠ Γ1

x1 V

⎛ 0⎞

θ = ⎜⎜ ⎟⎟ ⎝ 0⎠

Ω

Fig. 3. Integration domain.

The application of the Gauss–Ostrogradski theorem, to the first term of Eq. (31), enables us to obtain the following solution:

Mh ¼

Z

ðp_ j;kj  hk þ p_ j;k  nj  hk;j ÞdV þ

V

1 2

Z AðC1 Þ

ððrvij  ðui;j Þ;k þ ruij  ðvi;j Þ;k Þ  ððrvij  ui;j Þ;k þ ðruij  vi;j Þ;k ÞÞ  hk dV

ð32Þ

Considering the vector field ~ h, A(C1)can be replaced by the difference between the total surface X and the crown surface V. Then, the last term of Eq. (32) can be developed as follows:

1 2

Z

ððrvij  ðui;j Þ;k þ ruij  ðvi;j Þ;k Þ  ððrvij  ui;j Þ;k þ ðruij  vi;j Þ;k ÞÞ  hk dV Z 1 ¼ ððrv  ðui;j Þ;k þ ruij  ðvi;j Þ;k Þ  ððrvij  ui;j Þ;k þ ðruij  vi;j Þ;k ÞÞ  hk dV 2 X ij Z 1  ððrvij  ðui;j Þ;k þ ruij  ðvi;j Þ;k Þ  ððrvij  ui;j Þ;k þ ðruij  vi;j Þ;k ÞÞ  hk dV 2 V AðC1 Þ

ð33Þ

With relation (33), Eq. (32) can be rewritten as

Z

1 ðp_ j;kj þ ððrvij  ðui;j Þ;k þ ruij  ðvi;j Þ;k Þ  ððrvij  ui;j Þ;k þ ðruij  vi;j Þ;k ÞÞÞ  hk dV 2 ZV 1 þ p_ j;k  hk;j þ ððrvij  ðui;j Þ;k þ ruij  ðvi;j Þ;k Þ  ððrvij  ui;j Þ;k þ ðruij  vi;j Þ;k ÞÞ  hk dV 2 X

Mh ¼ 

ð34Þ

Finally, introducing Eqs. (33) and (30) into (34), and taking into account the conservation law (17), the Mh-integral form is given by:

Mh ¼

1 2

Z

X

ðruij  vi;k  rvij;k  ui Þ  hk;j dV þ

1 2

Z

X

ððrvij  ðui;j Þ;k þ ruij  ðvi;j Þ;k Þ  ððrvij  ui;j Þ;k þ ðruij  vi;j Þ;k ÞÞ  hk dV

ð35Þ

The first term of (35) provides us the mixed-mode separation for a stationary crack, see Destuynder and et al. [18]. The second term translated the dissipated energy induced by the crack growth process. 3. Generalization for a viscoelastic behavior and physical interpretations 3.1. Generalized Kelvin Voigt model The finite element model for the linear viscoelastic behavior is based on a generalized Kelvin Voigt model composed by N cells of Kelvin Voigt associated with a spring in series, Fig. 4. In this case, the Eq. (35) can be generalized for each elastic property as follow:

1 2

Z

ðrpij ðuÞ  vpi;k  rpij;k ðvÞ  upi Þ  hk;j dV Z 1 þ ððrpij ðvÞ  ðupi;j Þ;k þ rpij ðuÞ  ðvpi;j Þ;k Þ  ððrpij ðvÞ  upi;j Þ;k þ ðrpij ðuÞ  vpi;j Þ;k ÞÞ  hk dV 2 X with p ¼ ð0; 1; . . . ::NÞ

Mhpv ðu; vÞ ¼

X

ð36Þ

4457

R. Moutou Pitti et al. / Engineering Fracture Mechanics 75 (2008) 4450–4465 1 k ijkl

k ijklp

N k ijkl

0 k ijkl

σ kl

kl

η ε ij0

η

1 ijkl

ε ij1

p ijkl

ε ijp

η

N ijkl

ε ijN

ε ij Fig. 4. Generalized Kelvin Voigt model.

For the pth spring, rpij ðuÞ and rpij ðvÞ indicate the real and virtual stresses, respectively. In the same way, upi and vpi are real and virtual displacements for this spring, respectively. While referring to the definition of the energy release rate G, the superposition principle enables us to write [19]

Mhpv ðu; vÞ ¼ C p1 

u

u p v p K pI  v K pI K  K II þ C p2  II 8 8

ð37Þ

The real and virtual stress distribution is characterized by the real and virtual stress intensity factors for each fracture mode. They are noted (u K pI ; u K pII Þ and (v K pI ; v K pII Þ, respectively. C p1 and C p2 designate the reduced viscoelastic compliances in opening and shear modes, respectively [20]. Their form is given by

  qp  sp  qp1  sp2  C p1 ¼ 4  Re i  2 1p p s1  s2 2

2

with ppj ¼ Sp11  spj þ Sp12 ; qpj ¼

  pp  pp1 and C p2 ¼ 4  Re i  2p p s1  s2

ð38Þ

Sp22 þ Sp12  spj spj

and qpj ¼ cosðhÞ þ i  spj  sinðhÞ; with j 2 ð1; 2Þ

ð39Þ

Where (r,h) are polar coordinates, centered in the crack tip. spj are the roots of the characteristic equation: 4

2

Sp11  spj þ ð2  Sp12 þ Sp33 Þ  spj þ Sp22 ¼ 0 According to the plane behavior law,

ep11

¼

Sp11

r

p 11

þ

Sp12

r

p p 22 ; e22

¼

ð40Þ

Sp11 ; Sp22 ; Sp12 ;

Sp12

p 11

r

þ

Sp22

and p 22

r

Sp33

designate the compliance tensor components, such as

and

2  ep12 ¼ Sp33  rp12

ð41Þ

The expression (37) allows the mixed-mode separation by performing two distinct calculations, considering judicious values for the virtual stress intensity factors v K pI and v K pII : u

K pI ¼

8  Mhpv ðu; vÞðv K pI ¼ 1; v K pII ¼ 0Þ C p1

and

u

K pII ¼

8  Mhpv ðu; vÞðv K pI ¼ 0; v K pII ¼ 1Þ C p2

ð42Þ

The combination of expressions (37) and (42) gives the viscoelastic energy release rate Gpv as follow:

Gpv ¼ 1 Gpv þ 2 Gpv

with

1

Gpv ¼ C p1

ðu K pI Þ2 8

and

2

Gpv ¼ C p2

ðu K pII Þ2 8

ð43Þ

Gpv and 2 Gpv represent the specific energy release rate contribution induced by the free energy stored in the pth spring for each fracture modes, respectively. Finally, the total energy release rate for each mode is given by the following summation:

1

1

Gv ¼

X

1

Gpv

and

k

2

Gv ¼

X

2

Gpv

with p ¼ ð0; 1 . . . ::NÞ

ð44Þ

k

3.2. Physical interpretations For an orthotropic media, virtual field vp is given by the Sih’s singular form [21], defined in the crack tip vicinity by

rffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffiffi     0:5 r 1 r 1 p p p0:5 p p p0:5 p p0:5 v p ¼2   p1  s2  q1 Þ þ 2  K 2   pp1  qp1 Þ  Re p  Re p p  ðp2  s1  q2 p  ðp2  q2 2p 2p s1  s2 s1  s2 rffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffiffi     r 1 r 1 p p p0:5 p p p0:5 p p0:5 p p0:5 v p vp2 ¼ 2  v K p1   ðq  s  q  q  s  q Þ þ 2  K   ðq  q  q  q Þ  Re p  Re 2 1 2 1 2 1 1 2 2 2 1 2p 2p s1  sp2 sp1  sp2

vp1

v

K p1

ð45Þ

4458

R. Moutou Pitti et al. / Engineering Fracture Mechanics 75 (2008) 4450–4465

Dubois [20] and Chazal [22] have shown that mechanical fields, in the crack tip vicinity, can be defined by two viscoelastic stress intensity factors u K ra ða ¼ f1; 2g for plane problems) and two viscoelastic opening displacement intensity factors u K ea such as

½u1  ¼

rffiffiffiffiffiffiffiffiffiffi r  u K e2 2p

rffiffiffiffiffiffiffiffiffiffi r  u K e1 2p

and ½u2  ¼

ð46Þ

[u1] and [u2] are the components of the crack opening displacement. It designates relative displacement vector of crack lip nodes. Introducing a Boltzmann integral form into expressions (5) and (46), the relationship between stress and crack opening intensity factors takes the following form: u

K ea ¼

Z

t

C a ðt  sÞ 

0

ou K ra ds os

ð47Þ

Ca is the viscoelastic compliance function for a mode. Its form, which is analogue to a creep function, is given by

1

C a ðtÞ ¼

0

ka

þ

m N X 1 ka km t ia Þwith km m  ð1  exp a ¼ m ga m¼1 ka

ð48Þ

p

p

ka and gpa ðp 2 f0; 1 . . . ; NgÞ are the contribution of creep tensor components kijkl and gpijkl , respectively. u K pa is the distribution of the crack opening intensity factors such as u

K ea ¼ u K 0a þ

N X

u

Km a

ð49Þ

m¼1

with u

K 0a ¼

1

 u K ra 0

u

and

ka

K pa ¼

Z

t 0

p 1 ou K ra ka ðtsÞ Þ ds p  ð1  exp os ka

ð50Þ

By taking into account the Helmoltz strain energy density (3), the part of viscoelastic energy release rate for a a mode fracture using local creep properties, can be expressed by:

Gav ðtÞ ¼

1  8

Z

t

0

Z

t

0

½2  C a ðt  sÞ  C a ð2  t  s  bÞ

ou K ra ou K ra dsdb ðwithout summation onaÞ os ob

ð51Þ

The combination of expressions (43), (44) and (51) gives the partition of Gv such as

Gav ¼ a Gov þ

N X

a m

Gv

with

a p

Gv ¼

m¼1

1 1 u p2 a o 1 1 u o2   ½ K a  ; Gv ¼  0  ½ K a  8 kpa 8 ka

ð52Þ

4. Finite element formulation and algorithm crack growth 4.1. Finite element formulation for viscoelastic analysis According to a generalized Kelvin Voigt model, the relationship between stress components rkl and strain components eij is given by a Boltzmann’s integral including components of the four order creep tenor Jijkl(t) as follows [23,24]:

eij ðtÞ ¼

Z 0

t 

J ijkl ðt  sÞ 

orkl ds os

ð53Þ p

The evolution Jijkl(t) can be described using spring rigidities kijkl ðp 2 f0; 1; . . . ; NgÞ and dash-pot viscosities gm ijkl ðm 2 f1; . . . ; NgÞ such as, see Fig. 4:

J ijkl ðtÞ ¼

m N X kijkl 1 1 km t ijkl Þ with km o þ m  ð1  exp ijkl ¼ m gijkl kijkl m¼1 kijkl

ð54Þ

The spectral creep tensor decomposition (54) invites us to consider a distribution of the total strains such as [25]:

eij ¼

X k;l

Pijkl with Pijkl ¼ Pð0Þ ijkl þ

N X

ðmÞ Pijkl m ¼ ð1; . . . ::NÞ

ð55Þ

m¼1

Pijkl designates the part of strain component generated by the stress component rkl. The analogy between local and global behavior Eqs. (47) and (53), respectively and spectral decompositions (48) and (54), allows us to perform a finite difference integration and a linear approximation of stress and stress intensity factors for each time step Dtn = tn  tn1. With these considerations, incremental constitutive equations are defined by

4459

R. Moutou Pitti et al. / Engineering Fracture Mechanics 75 (2008) 4450–4465

Deij ðt n Þ ¼ Wijkl  Drkl ðt n Þ þ ~eij ðtn1 Þ e e ðtn1 Þ DK e ðt n Þ ¼ Wa  DK r ðt n Þ þ K a

a

ð56Þ ð57Þ

a

During the time increment Dtn, D eij(tn) and Drkl (tn) designate the strain and stress increments, respectively. In the same way, DK ea ðtn Þ and DK ra ðt n Þ are increments of crack opening and stress intensity factors for each a mode. ~eij ðtn1 Þ and e e ðtn1 Þ represent the global and local influences of the mechanical past history. Wijkl andWa are equivalent and reduced K a compliance functions, respectively. Their form is given by

! m m N X kijkl 1 1  expkijkl Dtn Wijkl ¼ 0 þ 1 with km m  m ijkl ¼ m gijkl kijkl  Dtn kijkl m¼1 kijkl   m m N 1 X 1 1  expka Dtn ka with km Wa ¼ 0 þ m  1 m a ¼ m g k  D t k n ka m¼1 a a a 1

ð58Þ ð59Þ

Introducing the creep function form (54) in the definition of the isothermal Helmoltz free energy F (3), it can be rewritten as [26]:

F¼

N X 1 0 1 m m m  kijkl  e0ij  e0kl þ k e e 2 2 ijkl ij kl m¼1

ð60Þ

with

e0ij ¼

1 0 kijkl

and em ij ¼

 rkl

Z 0

t

1



kijkl

0

m

 ð1  expkijkl ðtsÞ Þ 

orkl ds os

ð61Þ

The evaluation of F(tn), Eq. (60), for each time step, requires computation of the increments of e0ij and em ij . Coupling expressions (61) with the incremental formulation (56), we obtain:

De0ij ðt n Þ ¼

1 0 kijkl

 Drkl ðtn Þ and Dem ij ðt n Þ ¼

X

DPm ijkl ðt n Þ

ð62Þ

k;l

with m

m

kijkl Dt n DPm  1Þ  Pm ijkl ðt n Þ ¼ ðexp ijkl ðt n1 Þ þ

DK 0a ðt n Þ ¼

1 0

ka

1  expkijkl Dtn m kijkl

!

m

 rkl ðtn1 Þ þ

1 1  expkijkl Dtn 1 m  km kijkl ijkl  Dt n

!  Drkl ðt n Þ

 DK ra ðtn Þ

ð63Þ ð64Þ

and m

m

DðK ea Þm ðt n Þ ¼ ðexpka Dtn  1Þ  ðK ea Þm ðtn1 Þ þ

1  expka Dtn m ka

!  K ra ðt n1 Þ þ

  m 1 1  expka Dtn  DK ra ðt n Þ m  1 m ka  Dtn ka

ð65Þ

In order to solve the incremental formulation (56) with a finite element algorithm, the method proposed by Ghazlan [27] is employed. With the displacement vector increment noted {Dup}(tn), the balance equation can be written for all integration nodes, as

K pT  fDup gðtn Þ ¼ fDF pext gðt n Þ þ f e F p gðtn1 Þ

ð66Þ

K pT is the tangent matrix defined by using the Jacobean matrix B (connecting displacements and strains) and the equivalent compliance tensor defined by it components, Eq. (58):

K pT ¼

Z

BT  ½Wp 1  BdX

ð67Þ

X

fDF pext gðt n Þ denotes the increment of nodal force vector induced by the external increment loading. f e F p gðtn1 Þ is the supplementary viscous load vector which represents the complete mechanical history. It is given by

fe F p gðtn1 Þ ¼

Z

BT  Wp  fee p gðtn1 ÞdX

ð68Þ

X

where f~ep gðtn1 Þ is the strain vector defined for each integration point and derived from the strain history tensor f~ep gðt n1 Þ. 4.2. Finite element algorithm for crack growth process A geometry variation versus time (crack evolution) is considered in the crack growth model. For symmetric geometries and loadings, this crack growth is modeled by moving boundary elements in a semi-mesh. In a mixed mode configuration,

4460

R. Moutou Pitti et al. / Engineering Fracture Mechanics 75 (2008) 4450–4465

re-mesh Wi +1

Fi

•step 1

Fi

a + Δa

a

p

i

σ n

ect io

•step 0 a

Wi

σ i (t n )

( )

tn

eq

− Δσ i (t n )

•step 3 Δt n

•step 2

( )

Δσ i t n = σ i +1 t n − σ p i +1

ro j stre ss p

Fi

Δa σ i +1 (t n )

σ i +1 (t n ) instantaneous crack growth Δa

Wi +1

on rati igu f n t co alen u iv

(t n ) Fi

viscoelastic crack growth

Wi +1

σ i +1 (t n +1 ) Initial state

+ Δσ i (t n )

a + Δa

t n +1

n = n +1 Fig. 5. Viscoelastic crack growth algorithm.

this symmetric loading is lost. Then, it is necessary to operate a re-mesh driven by the crack tip advance. In our cases, a fixed crack orientation is supposed. To overcome this difficulty, the hereditary mechanical fields have been projected in the transformed mesh. This complex algorithm must separate the time and geometry variations, see Fig. 5. Firstly, all mechanical fields and the crack length a are supposed known at time tn. Then, stress and external loading, named ri(tn) and Fi(tn), respectively, are defined in the initial mesh noted Wi characterized by a crack length a (step 0). Firstly, the stress perturbation induced by a non-time dependent crack growth is evaluated: Step 1: Let us consider an instantaneous crack tip advance Da. New mesh noted Wi+1 is defined by re-meshing. With the same external loading, an elastic calculation provides the new stress state ri+1(tn). We note rpiþ1 ðt n Þ the geometric projection of ri(tn) on mesh Wi+1. Step 2: With a differentiation, the stress perturbation Dri+1 between the two mesh configurations is computed. Step 3: Dri+1 is applied as a cohesion stress (equivalent external loading) by using the superposition principle. This supplementary loading allows to close the new crack on Da. In this case, we obtain an equivalent configuration between steps 0 and 3 (same mechanical state and same geometry) but with two different meshes. Dri+1 can be interpreted as the stress cohesion in the process zone. Secondly, the crack length advance is fixed, and time is continuous: Employing the incremental viscoelastic procedure, the stress cohesion Dri+1 is employed as an external load vector during the time increment Dtn, see expression (68). In the time domain, the progressive un-cohesion of crack lips un-cohesion D a is obtained. The true mechanical state can then be released. Finally, fracture parameters can be computed and this algorithm is incremented. 5. Algorithm verification 5.1. CTS specimen and numerical mesh The global formulation, developed in this paper, is implemented in the finite element software Castem [28]. In order to validate M-integral expressions (36) and (44), a square, 80 mm wide CTS wood specimen (constant tension shear) is used. This specimen was developed by Richard [29] and used by Zhang [4,5]. The initial crack length chosen is 25 mm. The external load is an unitary loading applied to a perfect rigid steel arm (which presents a large crack growth zone), see Fig. 6. Points Aa and Ba with a 2 (1. . .7) are holes where forces can be applied with the angle b oriented according to the trigonometrically direction for different mixed-mode ratios. The simple opening mode is obtained by applying opposite forces in A1 and B1 with b = 0°. The loading b = 90°, in and, corresponds to the case of simple shear mode. Intermediary positions induce different mixed-mode configurations.

4461

R. Moutou Pitti et al. / Engineering Fracture Mechanics 75 (2008) 4450–4465

β = 90º°

a = 25mm

A7

β A7

B1 Wood specimen

β = 0º° A1

steel arms timber element

β B7

Fig. 6. CTS specimen.

C0 C2 C4 C6 C8

Crack tip Fig. 7. Crowns around the crack type.

1

G v (J / mm 2 )

2

(

Gv J / mm 2

4.0E-02

7.0E+00

C2

C4

C6

) C2

C8

C4

C6

C8

3.5E-02

6.0E+00

3.0E-02

5.0E+00

2.5E-02 4.0E+00

2.0E-02 3.0E+00

1.5E-02 2.0E+00

1.0E-02

1.0E+00

Δa / Δt (mm / h )

0.0E+00 0

1

2

3

4

5

6

7

8

9 10 11 12 13 14 15 ∞16

5.0E-03

Δa / Δt (mm / h )

0.0E+00 0

1

2

3

4

5

6

7

8

9 10 11 12 13 14 15 ∞16

Fig. 8. Evolution of energy release rate versus vs. time for each crown for pure mode (a: mode I, b: mode II).

4462

R. Moutou Pitti et al. / Engineering Fracture Mechanics 75 (2008) 4450–4465

1

G v (J / mm 2 )

2

3.5E+00

C2

C4

C6

(

Gv J / mm 2

2.5E-02

C8

) C2

C4

C6

C8

3.0E+00 2.0E-02

2.5E+00 1.5E-02

2.0E+00 1.5E+00

1.0E-02

1.0E+00 5.0E-03

5.0E-01 0.0E+00

Δa / Δt (mm / h )

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ∞16

Δa / Δt (mm / h )

0.0E+00 0

1

2

3

4

5

6

7

8

9 10 11 12 13 14 15 ∞16

Fig. 9. Evolution of energy release rate vs. time and crowns in mixed-mode (a: part of mode I; b: part of mode II).

Gv (J/mm² )⋅10 3

2

1

8.10

-5

5. 10

-6

4. 10

-6

3. 10

-6

2. 10

-6

1. 10

-6

Gv (J/mm² )⋅10 3

0˚ 45˚ 90˚ 6.10

4.10

2.10

-5

0˚ 45˚ 90˚

-5

-5

crowns

crowns

0

0 2

4

6

8

10

12

2

4

a: opening mode

6

8

10

12

b: Shear mode

Fig. 10. The non-dependence domain in a separate mode process (a: mode I; b: mode II).

In order to simplify analytic developments, a time proportionality for the creep tensor is chosen such as

JðtÞ ¼

1  C0 EðtÞ

ð69Þ

In which C0 is a constant and unit compliance tensor composed by a unity elastic modulus and a constant Poisson coefficient of 0.4. E(t) represents the tangent modulus for the longitudinal direction. In this context, the creep properties are given in terms of creep function by interpolating 1/E(t) with six Kelvin Voigt cells. 1/E(t) takes the following form:

" # 1 1 1 1 ð1  e22t Þ þ 74 ð1  e2:2t Þ þ 23 ð1  e2:210 t Þ 1 1 1 þ 74:3 ðtime given in hourÞ ¼ EðtÞ EX þ 1 ð1  e2:2102 t Þ þ 1 ð1  e2:2103 t Þ þ 1 ð1  e2:2104 t Þ 28 7:8 3:2

ð70Þ

in which EX(=15 GPa) is the elastic longitudinal Young modulus. C0 admits the definition for plane configurations:

2

1

6 t C0 ¼ 6 4 0

t EX EY

0

0

3

7 0 7 5

EX GXY

EY(=0.6 GPa) and GXY(=0.7 GPa) are the transverse and shear modulus, respectively.

ð71Þ

4463

R. Moutou Pitti et al. / Engineering Fracture Mechanics 75 (2008) 4450–4465

6. Numerical results In order to observe the non-dependence of the M-integral, the energy release rate is computed versus different integration crowns, see Fig. 7. The mixed-mode fracture separation is operated for different mixed mode ratios. The crack growth step is fixed to Da = 8 mm and the initial crack length is fixed to 25 mm. For different crack growth speeds, the simulation presents a crack growth process to a final crack length of 65 mm. For pure opening and shear modes, Fig. 8a (1Gv) and b (2Gv) illustrate the energy release rate versus different crack growth speed. In the same way, the part of opening mode (a) and the part of shear mode (b) in terms of energy release rate for a mixed-mode configuration (b = 45°) is presented in Fig. 9.

u

6

(

K1σ N ⋅ mm −3 / 2

5

)

u

elastic response (Δt = 0) Δt = 1h Δt = 10 days

(

K 2σ N ⋅ mm −3 / 2

0.5

)

elastic response (Δt = 0 ) Δt = 1h Δt = 10 days

0.4

4

0.3

3

2 0.2

1

a(mm) 0 20

30

40

50

60

a(mm)

0.1 20

70

30

40

pure opening mode

50

60

70

pure shear mode

Fig. 11. Stress intensity factor vs. crack length.

u

4

(

K1σ N ⋅ mm −3 / 2

)

u

21

elastic response (Δt = 0) Δt = 1h Δt = 10 days

18

3

(

K 2σ N ⋅ mm −3 / 2

)

elastic response (Δt = 0 ) Δt = 1h Δt = 10 days

15 12

2 9 6

1

3

a(mm)

0

a(mm) 0

20

30

40

50

part of opening mode

60

70

20

30

40

50

part of shear mode

Fig. 12. Stress intensity factor vs. crack length for a mixed mode configuration (b = 45°).

60

70

4464

R. Moutou Pitti et al. / Engineering Fracture Mechanics 75 (2008) 4450–4465

We observe, for all fracture modes, an increase in the energy release rate versus crack growth speed for each crown. This fact translates the impacts of viscoelastic effects and the relaxation phenomenon caused by the continuous un-cohesion, in the process zone. We also observe a stability of results which increases with the crown size. In this context, the non-dependence domain is proven with a specific simulation characterized by a crack growth step Da of 1 mm. The time increment is fixed to 240 h and the initial crack length is 25 mm. Fig. 10 shows results in terms of energy release rate distribution for a final crack length of 30 mm. We can observe the stability of energy release rate versus different crowns in mode I and mode II, respectively (0° = pure mode I, 45° = mixed-mode, 90° = pure mode II). These results show that the evolution of G, versus time, presents a non integration domain dependence defined by h field. For next results, the integration domain is based on crown C8. The crack growth step is fixed to 1 mm. The indicated time designates the corresponding time increment. According to the creep loading, the local stress field is characterized by a non time dependent function.

Gv (J/mm² )⋅10 3

1

4.5.10

-2

3.0.10

-2

1.5.10

-2

2

elastic response (Δt = 0) Δt = 1h Δt = 10 days

6. 10

-5

4. 10

-5

2. 10

-5

Gv (J/mm²)⋅103 elastic response (Δt = 0) Δt = 1h Δt = 10 days

a(mm) 0

a(mm) 0

25

35

45

55

pure opening mode

65

75

25

35

45

55

pure shear mode

65

75

Fig. 13. Energy release rate vs. crack tip advance for pure open and shear mode configurations.

Gv (J/mm² )⋅10 3

2

1

12.5.10

-2

-2

10.0.10

-2

1.5.10

-2

7.5.10

-2

1.0.10

-2

5.0.10 -2

5.0.10

-3

2.5.10

2.5.10

-2

2.0.10

elastic response (Δt = 0 ) Δt = 1h Δt = 10 days

a(mm) 0 25

35

45

55

part of opening mode

65

75

Gv (J/mm² )⋅10 3 elastic response (Δt = 0 ) Δt = 1h Δt = 10 days

-2

a(mm)

0 25

35

45

55

part of shear mode

Fig. 14. Energy release rate vs. crack tip advance for a mixed mode configuration (b = 45°).

65

75

R. Moutou Pitti et al. / Engineering Fracture Mechanics 75 (2008) 4450–4465

4465

The instantaneous response indicates reference values versus crack length in terms of stress intensity factors for each mode. This reference is obtained by performing elastic calculations. Different simulations are proposed for different constant crack tip speed characterized by the imposed time increment step. Fig. 11 shows results for pure opening and shear mode for u r K 1 and u K r2 , respectively. We can observe perfect time independence. For a mixed-mode configuration, the stress intensity factors distribution is plotted in Fig. 12. The time independence is also validated. However, the shear mode separation shows a law difference which can be explained by errors induced by stress field projection operated during the crack growth algorithm. Fig. 13 shows the energy release rate during the crack tip progression for pure opening and shear modes. In order to appreciate viscoelastic effects, different time increments have been applied. Results are compared to an elastic response corresponding to an instantaneous crack growth process. The same simulation is proposed in Fig. 14 for a mixed-mode configuration (b = 45°). These two last figures show clearly the continuous crack lip un-cohesion in the process zone caused by the relaxation time induced by viscoelastic effects. This fact is amplified versus the crack growth length and denotes also, the effects of viscoelastic material and the good behavior of the proposed model. 7. Conclusion and perspectives A new M-integral parameter is developed in the time domain in order to separate fracture modes during the creep crack growth in an orthotropic and viscoelastic media for mixed-mode configurations. This non path dependent integral permits us to compute the distribution of stress intensity factor and energy release rate in terms of opening and shear mode for mixed-mode configurations. The coupling between the M-integral and an incremental viscoelastic formulation is implemented in a finite element software and introduced in an automatic crack growth algorithm, considering the crack lip un-cohesion and the crack tip vicinity in a process zone. The performance and the validation of this model are clearly shown in terms of path-independence results. Relaxation effects can be observed on the energy release rate evolution versus crack growth speed. In order to validate this approach, we must develop an experimental strategy to obtain stable crack growth by defining a specific specimen for mixed-modes crack configurations. In this approach, the definition of crack growth criteria, based on an energy release rate balance for opening and shear modes, is necessary. References [1] Bian L, Fawaz Z, Behdiman K. A mixed mode crack growth model taking account of fracture surface contact and friction. Int J Fract 2006;139:39–58. [2] Dubois F, Chazal C, Petit C. A finite element analysis of creep-crack growth in viscoelastic media. Mech Time-Depend Mater 1999;2:269–86. [3] Chang JH, Wu DJ. Computation of mixed-mode stress intensity factors for curved cracks in anisotropic elastic solids. Engng Fract Mech 2007;74:1360–72. [4] Ma S, Zhangl XB, Recho N, Li J. The mixed-mode investigation of the fatigue crack in CTS metallic specimen. Int J Fatigue 2006;28:1780–90. [5] Zhang XB, Ma S, Recho N, Li J. Bifurcation and propagation of a mixed-mode crack in a ductile material. Engng Fract Mech 2006;3:1925–39. [6] Van der Put TACM. A new fracture mechanics theory for orthotropic materials like wood. Engng Fract Mech 2007;74:771–81. [7] Griffith AA. The phenomena of rupture and flow in solids. Philos Trans Roy Soc Lond 1921;221:163–97. [8] Nielsen LF. A lifetime analysis of cracked linear-viscoelastic material with special reference to wood. In: Proceedings of IUFRO S5.02 Sweden; 1982. [9] Shapery RA. Correspondence principles and a generalized J integral for large deformation and fracture analysis of viscoelastic media. Int J Fract 1984;25:195–223. [10] Dugdale DF. Yielding of steel sheets containing slits. J Mech Phys Solids 1960;8:100–4. [11] Zhang W, Deng X. Mixed-mode I/II fields around a crack tip with a cohesive zone ahead of the crack tip. Mech Res Commun 2007;34:172–80. [12] Dubois F, Chazal C, Petit C. Viscoelastic crack growth process in wood timbers: an approach by the finite element method for mode I fracture.Int J Fract 2002;113:367–88. [13] Rice JR. A path independent integral and the approximate analysis of strain concentrations by notches and cracks. J Appl Mech 1968;35:379–85. [14] Staverman AJ, Schwarzl P. Thermodynamics of viscoelastic behavior. Proc Acad Sci 1952;55:474–92. [15] Noether E. Invariant variations problem. Transport Theory Stat Phys 1918;1:183–207. [16] Chen FMK, Shield RT. Conservation laws in elasticity of the J-integral type. J Appl Mech Phys 1977;28:1–22. [17] Attigui M, Petit C. Mixed-mode separation in dynamic fracture mechanics: new path independent integrals. Int J Fract 1997;84:19–36. [18] Destuynder Ph, Djaoua M, Lescure S. Quelques remarques sur la mécanique de la rupture élastique. J Méc Théor Appl 1983;2:113–35. [19] Dubois F, Chazal C, Petit C. Modelling of crack growth initiation in a linear viscoelastic material. J Theor Appl Mech 1999;37:207–22. [20] Valentin G, Morlier P. A criterion of crack propagation in timber. Matér Construct 1982;88:291–8. [21] Sih GC. Strain energy density factor applied to mixed mode crack problems. Int J Fract ;10:10;305305–321. [22] Chazal C, Dubois F. A new incremental formulation in the time domain of crack initiation in an orthotropic linearly viscoelastic solid. Mech TimeDepend Mater 2001;5:3–21. [23] Chrestensen RM. Theory of viscoelasticity: an introduction. New York: Academic Press; 1982. [24] Salençon J. Viscoélasticité. Presse de l’école nationale des ponts et chaussées Paris; 1983. [25] Mandel J. Cours de mecanique des milieux continus. Paris; 1966. [26] Dubois F, Petit C. Modelling of the crack growth initiation in viscoelastic media by the Ghv-integral. Engng Fract Mech 2005;72:2821–36. [27] Ghazlan G, Caperaa S, Petit C. An Incremental formulation for the linear analysis of thin viscoelastic structures using generalized variables.Int J Numer Meth Engng 1995;38:3315–33. [28] Charvet-Quemin F, Combescure A, Ebersol L, Charras Th, Millard A. Méthode de calcul du taux de restitution de l’énergie en élastique et en non linéaire matériau. Report DEMT; 1986. p. 86/438. [29] Richard HA, Benitz K. A loading device for the creation of mixed mode in fracture mechanics. Int J Fract 1983;22:R55–8.