Dynamic crack growth: Analytical and numerical cohesive zone models approaches from basic tests to industrial structures

Engineering Fracture Mechanics 90 (2012) 1–29

Contents lists available at SciVerse ScienceDirect

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

Dynamic crack growth: Analytical and numerical cohesive zone models approaches from basic tests to industrial structures G. Debruyne a,⇑, J. Laverne a,1, P.E. Dumouchel b a b

LaMSID EDF CNRS CEA 8193, 1 av du Général de Gaulle, 92141 Clamart, France LMM-UPMC, Paris VI University, 4 place Jussieu, 75252 Paris Cedex 5, France

a r t i c l e

i n f o

Article history: Received 22 June 2006 Received in revised form 21 February 2012 Accepted 2 April 2012

Keywords: Dynamic fracture Crack arrest Thin films Material interfaces Cohesive zone models Pressurized components

a b s t r a c t Some dynamic crack growth problems are investigated, from basic academic tests to actual industrial situations, using analytical or numerical methods. The purpose of this paper is to describe the successive fast crack growth and arrest, driven by a discontinuity in fracture toughness. The main purpose of this article is the description of both crack growth and arrest, with the same governing equations for a wide range of examples. The initial problem is the peel test of a thin film bonded to a flat rigid surface. The film is divided in two zones of different bonding properties. This entails a fast debonding process, followed by an arrest. The problem is analytically solved, with the combined use of the characteristics method and the Griffith criterion. Then, a bimaterial Double Cantilever Beam (DCB) is considered, the materials concerned having different surface energies. This test involves a dynamic crack growth, which is numerically handled with Cohesive Zone Models (CZM). These models are derived from general energy concepts of the fracture process. Comparative predictions with dynamic and static analyses are discussed for these two problems. Finally, a real survey of a Pressure Water Reactor vessel shell, affected by an edge crack and submitted to an inner pressure loading, is carried out with CZM. Two situations are investigated. First, the initial flaw is assumed to propagate in a homogeneous base steel of constant toughness. Secondly, a small elastic zone of low toughness is embedded in the base metal along the crack path. We will focus on the possible crack jump and arrest in these two configurations, depending on whether the base metal exhibits elastic or plastic behavior, and on the relative toughness of the small zone with respect to the surrounding material. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction Predicting dynamic fracture propagation and arrest constitutes a major challenge. Even if the loadings of a structure are slowly prescribed, the presence of a material interface between two zones with contrasting fracture toughness, or an inclusion of low toughness, may lead to local inertia effects and to a possible crack jump (or rather a high speed propagation), possibly followed by an arrest. This article proposes some physical understanding, analytical predictions and numerical simulations of dynamic crack propagation and arrest for a variety of structures, presenting such material interfaces. From here onwards, the term ‘‘material interface’’ will be related to a toughness jump between two materials (which may have the same elastic properties) while the related term ‘‘interface’’ will denote a cohesive separation zone.

⇑ Corresponding author. Tel.: +33 147654575; fax: +33 147654118. 1

E-mail addresses: [email protected] (G. Debruyne), [email protected] (J. Laverne), [email protected] (P.E. Dumouchel). Tel.: +33 147653201; fax: +33 147654118.

0013-7944/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.engfracmech.2012.04.002

2

G. Debruyne et al. / Engineering Fracture Mechanics 90 (2012) 1–29

Nomenclature a cantilever beam crack length {Ci} characteristic lines set CR Rayleigh wave speed c one dimensional wave speed d horizontal displacement in (9) Es internal energy in (10) Ec kinetic energy in (11) Eelas elastic energy in (32) total energy in (32) Etot E Young modulus ET hardening modulus G energy release rate Gic surface energy for the material i H Heaviside function h beam height I moment of inertia KIC fracture toughness for mode one ‘D = ‘ totally debonded length ‘c partially debonded length debonding arrest length ‘A ‘_ debonding velocity L Double Cantilever Beam length Lcoh cohesive zone length N tensile force in the film Nw number of wave round trips P⁄ critical inner pressure in the vessel, at the crack onset t time variable u film deflection _ U € U; U; Double Cantilever Beam deflection, velocity and acceleration V0 initial film deflection velocity x space variable v film deflection velocity a ¼ G2c =G1c surface energy ratio b, c Newmark scheme parameters d cohesive displacement jump j cohesive displacement jump separating linear and dissipative mode j0 initial yield displacement jump in the cohesive model m Poisson ratio q specific volumic mass ric maximum cohesive stress for the material i rY plastic yield stress x film rotation

The first example is the peel test consisting in the debonding of a thin flexible film from a rigid flat surface. This is one of the few cases for which an analytical solution is available. This problem has been analyzed in a dynamic way by Freund [1], who considered an initial blunted crack, rapidly running along an homogeneous bonding surface, and whose propagation is governed by the Griffith criterion [2]. This problem is equivalent to the shear beam model in a Double Cantilever Beam (DCB) specimen [3], so that this simple 1D film model allows to account for the observed crack propagation in the DCB specimen. It is revisited here under some slightly different conditions, such as the existence of a material interface between two zones of the glued surface [4], or a narrow weak zone with a low fracture energy. This generates strong dynamic effects, due to the fracture energy discontinuity across these material interfaces. The closed-form dynamic solution of this problem is entirely determined in the first section of the paper, giving some clarifications for the analysis of more complex structures. A specific paragraph is dedicated to the presentation of a cohesive zone model. The notion of cohesive forces was introduced by Barenblatt [5]. An early model with a governing cohesive law was proposed by Dugdale [6]. CZM has been extensively developed by many authors since then, among which Hillerborg et al. [7] who associated a surface energy and a cohesive law dedicated to concrete, Needleman [8] who defined a cohesive law derived from a potential, and Tvergaard [9] who introduced an irreversibility notion in the cohesive law and a post decohesion behavior. Dynamic investigations have been also performed [10]. The model used here is based on an original framework [11], derived from the Francfort– Marigo theory [12,13], and is available with two types of cohesive laws.

3

G. Debruyne et al. / Engineering Fracture Mechanics 90 (2012) 1–29

The second example is the deflection of a Double Cantilever Beam with an initial edge crack, which grows across material interfaces (separating zones with contrasting fracture energy). The initial crack is growths slowly into the material with a constant surface energy, under quasi-static loading. Then, the crack tip reaches the material interface and the loading is frozen. But, owing to the surface energy jump, the crack is pushed ahead and runs through the lowest tough material, before a definitive arrest. A number of works have been already performed about this problem, but with simplifications, avoiding the consequences of waves reflections occurring during the dynamic stages [14,15]. Because there exists no exact analytical solution to this problem (except for the special case of the antiplane shear loading [24], similar to the film peeling problem) numerical simulations, using a Finite Element Model in plane strain, with dynamic conditions, are carried out using CZM described before. The effect of the fracture energy ratio on crack growth and arrest is studied, along with the effect of the critical cohesive stress, and the cohesive law shape (two cohesive laws, respectively based on a polynomial and an exponential energy shape, have been used). A comparison with a predictive node release method is performed [16,17], and an alternate static method is introduced to predict the crack arrest. The case of an embedded brittle zone in the specimen is also investigated, which makes a link with the upcoming industrial survey. The last section deals with the dynamic analysis of a Pressure Water Reactor vessel shell affected by an undercladding flaw. Service life extension of Pressure Water Reactor (PWR) vessels is an important issue for numbers of nuclear operators [18]. This goal is partially subjected to improvement of safety margins. These safety margins may be increased if it can be demonstrated that under the most severe accidental loadings, the crack growth resulting from the presence of initial flaws is limited [19]. For some accidental transient loadings, the possible suddenness of the propagation and arrest, due to variations in toughness in part of the vessel, caused by neutron irradiation and a temperature collapse may involve important dynamic effects [20]. Furthermore, some flaws characterized by a low toughness, and located along the crack path may enhance these phenomena [21]. The geometry is reduced here to a torus designed by an axisymmetric geometry, submitted to a pressure loading on the inner cladding surface. Contrary to previous examples, the geometry and loading involve unstable crack growth even in a homogeneous medium (at least within the elastic frame), but the arrest may occur because of changes in toughness through the vessel wall. Cohesive elements are used here for the analysis of the onset and the propagation of a 6 mm deep initial subclad flaw. Due to the simplified axisymmetric geometry, the latter is considered as a circumferential crack all around the shell. A first analysis is performed, assuming that the base metal toughness is uniform, and the behavior is either elastic or elasto-plastic. A second analysis is carried out, considering that a small inclusion with a low toughness, embedded in a zone along the crack path, involves a fracture energy jump which will speed up the crack growth because of favorable dynamic conditions (hence, this configuration presents some similarities with the preceding examples). Considerations about the metal base behavior (elastic or elasto-plastic) on crack growth stability and arrest are emphasized. 2. Analytical aspect of the dynamic fracture: film peeling from a bonding surface with two surface energy zones This section is devoted to analytical investigations of the peel test of a thin film glued on a rigid flat surface, with one or two material interfaces involving surface energy jumps. This specific academic one-dimension problem serves as basic test for the following investigations of more complex 2D and 3D structures using Finite Element Method associated to cohesive zone models. The film is considered as inextensible but perfectly flexible, with an homogeneous lineic mass density q. Small displacements assumption and Griffith criterion are taken into account for the debonding process. Two situations are considered: the first one is related to a bonding surface, divided into two zones of respective surface energy G1c and G2c 1 with G2c < G1c (Fig. 1a). The second one deals with a small zone of a surface energy G2c , embedded in a homogeneous bonding surface with a surface energy G1c with G2c < G1c (Fig. 1b). In both cases, due to slowly prescribed deflection at the end of the film, the peeling process is quasi-static in the toughest zone. Then, due to the toughness jump at the material interface, the debonding runs into or throughout the low tough zone before it stops. Therefore, the conditions of debonding arrest and its location need to be considered. In both situations, the loading continues to grow with an infinitesimal rate, even during the fast peeling phase,

u ( 0, t ) = V0t

u0 ( t ) = u ( 0, t ) = V0t N

N

Z

d ( x, t )

0

Gc1

u ( x, t )

(t )

Gc2

x1

u ( x, t )

x

0

Gc1

(t )

Gc2

x1

Gc1

x2

Fig. 1. Peeling of a thin film from a surface: (a) One material interface. (b) Two interfaces (weak inclusion).

x

4

G. Debruyne et al. / Engineering Fracture Mechanics 90 (2012) 1–29

contrary to the Freund analysis [1] where the loading is frozen during the fast debonding stage. For the case of one material interface, we may refer to Dumouchel et al. [4], where the authors, start from a revisiting quasi-static theory of fracture [12], and perform a complete dynamic analysis whose solution converges toward a specific quasi-static solution for the arrest point, when the loading rate tends to zero (for the other kinematics stage, only the dynamic solution holds). For the case of a brittle inclusion embedded in a safe material, we may refer to Dumouchel [23] who investigates the case of multi surface energy zones. In this section, some of these results are summarized and only features prone to extension to computation of actual structures, are presented. 2.1. The film peeling from a surface with two contrasting surface energies The main features of this model are depicted in Fig. 1a (the small deflection is magnified on the figure to make it more clear). The film is initially completely bonded on the surface, the debonding/bonding separation interface is described by the space-variable x, x P 0, with x = ‘(t). A horizontal dead load tensile force of constant magnitude N is prescribed to stretch thepffiffiffiffiffiffiffiffiffi flexible film. Then, at x = 0, a vertical displacement, with a constant rate, u0(t) = V0t is slowly applied ffi (V 0  c ¼ N=q, where c is the wave velocity and u0(t)  ‘(t)). Meanwhile, the film debonding activates a (cohesive) surface energy Gc, internally stored in the separation interface between the film and the supporting surface. u(x, t) denotes the film deflection (the horizontal displacement d(x, t) will be considered of a second order magnitude in this problem where the film is inextensible). For sake of simplicity, a Griffith criterion is adopted, so that the debonding energy does not depend on the separation deflection u (as is the case for usual cohesive energy which will also be considered further in this article), and takes the constant values G1c for 0 6 x < x1 and G2c for x P x1, with a ¼ G1c =G2c > 1. Some definitions and concepts are discussed here without demonstrations, further details can be found in [4,22,23]. 2.1.1. Equations of motion In addition to the classical boundary conditions described before, the existence, in the space–time domain XT = [0, +1]  [0, T], of a set of m smooth curves S ¼ fC i ; i ¼ 1; mg is considered. On these curves, velocity (v = u,t) and rotation (x = u,x) jumps generate shock waves. fC ‘ g is a particular curve defining the spatial–temporal path of the peeling point x = ‘(t) (which is similar to a crack tip). On the regular part XT =S, the film motion is governed by the classical wave equation:

Nu;xx  qu;tt ¼ 0;

x 6 ‘ðtÞ; ðx; tÞ R fC i g:

ð1Þ

On the set of points where rotations x and velocities v are discontinuous, the fields have to verify both the equation of motion and Hadamard compatibility relations, so that on {Ci}, the Rankine–Hugoniot conditions hold:

ckxk þ kv k ¼ 0;

ðx; tÞ 2 C þi –C ‘ ;

ð2Þ

ckxk  kv k ¼ 0;

ðx; tÞ 2 C i –C ‘ ;

ð3Þ

_ xk þ kv k ¼ 0; ‘k

ðx; tÞ 2 C ‘ ;

ð4Þ

Cþ i

C i

where and are respectively related to forward and backward shock waves in the (x, t) plane, and kxk, kvk denote the rotation and velocity jumps throughout {Ci}. 2.1.2. Griffith law of debonding The peeling process is assumed to be governed by the Griffith law, and formulated in terms of dynamic energy release rate G(t), which is defined, in a general 2D situation, as the limit of a path integral when the path tends to the crack tip or by an overall energy rate balance [1]. For our 1D case, the expression of G(t) reduces to (with x0 = x(‘, t)) [1]:

GðtÞ ¼

qc2 2

"

1

# ‘_2 x2 : c2 0

ð5Þ

As G(t) P 0, this last expression points out that ‘_ 6 c. The Griffith’s law is then described by the two following relations:

_ P 0; ‘ðtÞ

GðtÞ 6 Gc ð‘Þ

ð6Þ (

_ ‘ðtÞ½G c  GðtÞ ¼ 0;

Gc ¼ (

_ ‘ðtÞ½G c  GðtÞ ¼ 0; Gc ¼

G1c ; x 6 x0 G2c ; x > x0

;

for the case depicted in Fig: 1a or

G1c ; x 6 x0 ; x P x1 G2c ; x0 < x < x1

;

for the case depicted in Fig: 1b:

ð7Þ

ð8Þ

G. Debruyne et al. / Engineering Fracture Mechanics 90 (2012) 1–29

5

2.1.3. Energy contributions The energetic interpretation of the arrest suggested in [4,23], involves the consideration of the three energy contributions: the potential energy of the tensile force N, the internal energy and the kinetic energy. By assuming the inextensibility and the small magnitude of u, an estimated form of the potential energy generated by the tensile force N may be considered:

PðuÞ ¼ Ndð0; tÞ 

N 2

Z 0

þ1

u2;x ðx; tÞ dx

ð9Þ

The internal energy is expressed as:

es ðuÞ ¼

Z

þ1

Gc HðuÞ dx where HðuÞ ¼

0



0;

if

u60

1; if

u>0

:

ð10Þ

The kinetic energy, which is considered to depend only on the transverse velocity, is approximated by:

ec ðuÞ 

1 2

Z 0

þ1

qu2;t dx

ð11Þ

Note that in Eq. (2) H(f) is representative of a particular non-smooth cohesive law corresponding to a Griffith criterion. This function may be generalized to any cohesive law shape, keeping the same analytical solving process. But the investigation should then be much more complex, involving besides the totally debonded length evolution ‘_D ðtÞ, a partially debonded zone evolution ‘_C ðtÞ [22], (here, in the framework of the Griffith criterion ‘D(t) = ‘C(t) = ‘(t)). 2.1.4. General structure for the dynamic solution The curves {Ci}, i = 1, m, describing backward (slope 1/c) and forward (slope +1/c) shock waves, and the front of debonding fC ‘ g define characteristic lines in the (x, t) plane, circumscribing sectors {Si}, including a specific solution for the fields (x, v) (see Fig. 2 where for the sake of visibility the picture is not at the right scale). Inside these sectors, piecewise constant   _ tÞ ¼ ‘_i (the relation (1) holds solutions are sought, so that the fields be constant within each Si xðx; tÞ ¼ xi ; v ðx; tÞ ¼ v i ; ‘ðx; therefore automatically). To show the difference between the successive shock waves reflection times at x = 0 and x = ‘(t), they will respectively be denoted by t2i, and t2i1. The solutions are distinguished between two classes, depending on whether the sector design number is even or odd. A recurrent process links the even sectors solutions to the odd sectors fields [4]. Therefore, starting from S0, which is at rest (u = 0), a general piecewise constant solution for rotations and velocities fields, and for the debonding speed is estimated. The process is performed in a recurrent way, sector by sector, connecting S2i1 to S2i fields, using in particular the relations (2), (3), (4), (5), and (7) [4]. Therefore, debonding occurs in four stages. Up to ‘ = x1, the debonding point moves at a moderate speed, and its path converges to the classic quasi-static solution, when the loading rate tends to zero. Then, the material interface involves a fast ‘‘crack growth’’ while a backward shock wave arises and is reflected at x = 0 in a forward wave, which intersects the moving debonding front. Hence, the debonding speed is affected and the wave is transformed in a backward wave reflected again at x = 0 and so on. This process continues until a possible equilibrium state is established. Each stage is investigated with a dynamic formulation, so that the solutions are available for any V0 value. Note that only the particular situation where V0  c is investigated here.

Fig. 2. Shock waves and path of the peeling point in the (x, t) plane.

6

G. Debruyne et al. / Engineering Fracture Mechanics 90 (2012) 1–29

2.1.5. First stage: dynamic ‘‘slow’’ growth (Sector S1) At the beginning of the debonding process, before the peeling point (also referred to as ‘‘crack tip’’) reaches the position x1, this latter moves slowly, because V0  c (Sector S1 in Fig. 2). Following the Eqs. (2)–(7) the crack length evolution is:

‘ðtÞ ¼ V 0 t

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi N

qV 20 þ 2G1c

¼ V 0t

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qc2

qV 20 þ 2G1c

ð12Þ

;

and the fields in S1 are:

‘_ ¼ V 0

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qc2 ; 1 2Gc þ qV 20

V x1 ¼  0 c

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2G1c þ 1; qV 20

v1 ¼ V 0;

ð13Þ

Note that when V0 ? 0, the solution converges to the standard quasi-static solution, which is obtained from the equilibrium equation and the Griffith criterion:

‘ðtÞ ¼ V 0 t

sffiffiffiffiffiffiffiffiffi qc2 2G1c

;

sffiffiffiffiffiffiffiffiffi 2G1c ; x1 ¼  qc2

ð14Þ

This is not the case for the next kinematics stage, as we will see in the next paragraph. Furthermore, in (13), the ‘‘crack tip’’ velocity is limited to c when V0 ? 1 while it is unlimited in the quasi-static relation (14). 2.1.6. Second stage: dynamic ‘‘fast’’ growth (Sector S2) Due to the bonding strength discontinuity, a shock wave is generated at the material interface x = x1, while the crack tip _ 1 ; t 1 Þ. Thus, abruptly propagates into the second zone activating the surface energy G2c at the time t1, with a non-zero speed ‘ðx rotations and velocities jumps are generated at the point O(x1, t1) and propagate along the characteristic line until the time t2 when the backward wave is reflected at x = 0 to a forward wave up to (x1 + ‘A, t3) (see Fig. 2). During this back-and-forth wave propagation, the crack tip velocity is constant and not proportional to V0 as before, so that even when V0 ? 0, the crack velocity keeps a finite non-zero value depending on the energy ratio: 1 2 _ 1 Þ=c ¼ GC  GC ¼ a  1 : ‘ðt 1 GC þ G2C a þ 1

ð15Þ

2.1.7. Third stage: arrest (Sector S3) At the time t3 the shock wave reaches the moving peeling point and a ‘‘crack arrest’’ occurs during a fair number Nw of wave back-and-forth propagation between x = 0 and x = x1 + ‘A, depending of V 0 ; G1C ; G2C [4]. For V0 strictly equal to zero, Nw = 0, and the arrest is permanent, as found by Freund [1]. For V0 ? 0:

‘A ¼ ða  1Þx1 :

ð16Þ

This solution is in disagreement with the classic quasi-static one, obtained by expressing the Griffith criterion both at x = x1 ðG ¼ G1C Þ and x = x1 + ‘A ðG ¼ G2C Þ, and leading to the following underestimated arrest:

pffiffiffi ‘stat ¼ ð a  1Þx1 : A

ð17Þ

Let us note that the crack velocity (15) and arrest relation (16) only depend upon the surface energy ratio a. In particular, they do not depend on the energy values ðG1C ; G2C Þ themselves. 2.1.8. Fourth stage: slow growth around the homogeneous media solution (Sector S4) ffiffiffiffiffiffiffiffiffiffiffiffiffiffi resumes with a velocity oscillating around the solution in a homogeneous media with a At t = t2N+3, the peelingr process 2 surface energy G2C : ‘_ ¼ V 0 2G2qþc qV 2 . c

0

2.1.9. Numerical solutions achieved with the step by step procedure The solutions of this problem (represented in Fig. 3 with an adimensional scale) are analytically computed , using the characteristics network , and they have been numerically achieved with the recurrent process introduced in [4]. In Fig. 3 the shock wave paths and crack front evolution including arrest are illustrated with the ratio a = 2, for the dynamic and the standard quasi-static analysis. For the first stage (up to x1), both cases yield similar solutions (because V0  c). Then, the crack jump arising in the simple quasi-static analysis is underestimated as mentioned before ð‘stat < ‘A Þ. In the third A stage, dynamic arrest occurs with a dynamic analysis while quasi-static solution exhibits no phase of arrest, but it recovers the dynamic solution after the stage of arrest. At last, the quasi-static kinematics fits to the ‘‘dynamic one’’ after the dynamic arrest, with a dynamic solution oscillating around the quasi-static one with alternating short phases of acceleration and slowing down.

7

G. Debruyne et al. / Engineering Fracture Mechanics 90 (2012) 1–29

Fig. 3. Peeling of the film with one material interface, for dynamic and static analysis in the x, t plane.

2.1.10. Energy balance investigation If the investigation aims at predicting arrest length ‘A, and if the complete description of the debonding process is no subject of interest, some simple energy considerations may provide the arrest length. Starting from x = x1, and assuming a quasistatic solution, the deflection is: uðxÞ ¼ u‘0 ð‘  xÞ. According to relation (9), the potential energy is Pð‘Þ ¼ pated energy is

Nu20 2‘

and the dissi-

1

eS ð‘Þ ¼ G2c ‘ ¼ Gac ðx1 þ D‘Þ. Therefore, as Griffith criterion holds for ‘ = x1, we may write Pðx1 Þ ¼ G1c x. Thus

the conservation law may be expressed by the fact that the energy ratio is equal to one:

W TOT ð‘Þ Pðx1 Þ

eS ð‘Þ ¼ Pð‘Þþ ¼ x1xþ1D‘ þ aDx‘1 ¼ 1. Pðx1 Þ

Therefore, the solution of the above equation is x1 + D‘ = ax1, which is exactly the dynamic solution for ‘ = ‘A (see relation (16)). Indeed, the static energy balance only holds for ‘ = ‘A because the equilibrium equation then holds, but it does not during the dynamic stage. To illustrate this method, Fig. 4 represents the different energy quantities for x1 = 1, a = 2, along with the dynamic analysis. For the purpose of comparison, note that the minimum of

W TOT ð‘Þ Pðx1 Þ

corresponds to the Griffith criterion pffiffiffi a  1.

associated to equilibrium condition (that is not fulfilled here), and predicts an underestimated arrest length D‘stat ¼ A

1,2

1

Minimum of energy=Griffith criterion

Energy/W0

0,8

0,6

dynamic total energy (strain energy+kinetic energy+dissipation) static total energy (strain energy + dissipation Gc(L-x0)

0,4

dynamic strain energy static strain energy kinetic energy

0,2

Static arrest prediction

exact arrest solution

0 0

0,1

0,2

0,3

0,4

0,5

0,6

0,7

0,8

0,9

(L-x1)/x1 Fig. 4. Energy balance in the static and the dynamic analysis for a = 2.

1

8

G. Debruyne et al. / Engineering Fracture Mechanics 90 (2012) 1–29

Fig. 5. Peeling of the film with two material interfaces, for dynamic and static analysis in the x, t plane.

2.2. Film peeling from a surface with a low surface energy inclusion The particular case of a zone with a low surface energy G2c , with a length ‘Z ¼ x2  x1 , located in the space interval [x1, x2] and embedded in a tougher media with an uniform energy G1c (see Fig. 1b) is considered now. This model may be seen as an analogy with some industrial problem where some flaw may weaken the safe material. The boundary conditions are the same as in the previous paragraph, and we consider again a ¼ G1c =G2c > 1. This problem is formulated with the equations presented before, and solved in the same way, but with a further complexity involved by an additional wave reflecting from the second material interface at x = x2. As depicted in Fig. 5, the peeling path exhibits a running growth through the inclusion, involving shock waves first at (x = x1) and then at (x = x2), followed by a stop-and-go propagation in the tougher media until it finally gets around the quasi-static solution. Each short arrest and restart corresponds to the arrival of a shock wave from one of the inclusion interface (x = x1 or x = x2) to the peeling point. Considering the arrest in the case of a constant load from ‘ = x1, two situations are possible: – The inclusion is wide enough for the shock wave to hit the crack tip and stops it within the inclusion. The problem is then equivalent to the one material interface problem. – the inclusion is short enough for the peeling point to reaches its end (x = x2) before the wave returns from the prescribed loaded end (as illustrated in Fig. 5) and the arrest takes place at x = x2. In this case, the static approach presented earlier fails to predict the exact solution, thus making dynamic computations mandatory.

3. 2D interface cohesive zone model In this part a two dimensional cohesive zone model is introduced to predict the crack evolution in a given direction. The mechanical problem is defined by an energetic formulation. Each potential crack surface in the material is associated to a surface energy density depending on the displacement jumps. This energy takes into account the dissipative process during the crack growth. According to the Francfort–Marigo theory [12], the displacement field solution is obtained by minimizing the sum of the bulk and surface energies at each loading step. The first cohesive models were introduced by Barenblatt [5] and Dugdale [6] at the begining of the sixties. Observing that infinite stresses at the crack tip, predicted by the elastic model, do not have physical meaning, they postulate the existence of cohesive forces acting on both crack lips all along a ‘‘cohesive zone’’ ahead the crack tip. This cohesive zone may be related to experimental observations such as the breaking zone of atomic bonds in the cleavage range or the development of plastic zones ahead the crack tip (Dugdale found good agreement between measured lengths of plastic zones in steels and predictions based upon his strip model [6]). It corresponds to the process zone between the sound material and the crack surface. In the 1970s Hillerborg et al. [7] introduced the concept of fracture energy and proposed a cohesive force versus crack opening relationship for concrete. A large number of cohesive model were developed ever since, with a wide range of applications. Needleman [8] and Tvergaard [9] may be quoted for interface debonding modeling concerning composite materials with metal matrix. Our approach follows these authors, but the most original feature of our study is based on the least energy principle, derived from the latter developments of the Francfort–Marigo theory [12] (see also Charlotte et al. [13] and Laverne [11]). In the following paragraph two interface cohesive laws are developed. Then interface finite element properties are explained. Finally, the different steps of the minimization problem to be solved will be detailed.

G. Debruyne et al. / Engineering Fracture Mechanics 90 (2012) 1–29

9

3.1. Interface cohesive laws In this section two interface cohesive laws are suggested for mode I or mode II crack opening (in the next sections, only applications with mode I will be considered). The first law adopts an exponential expression for the force opening displacement relationship and the second one with a polynomial expression. The last one is relevant for a real free stress cracked surface (cohesive forces are exactly nil when the crack opening reaches a critical value). The exponential one only describes a ‘‘damaged zone’’ all along the crack lips because cohesive forces smoothly vanish, but they never reach zero. Simulations will be performed to evaluate the crack evolution sensitivity to the cohesive law. Let’s denote by C some surfaces in the body subjected to displacement jumps and W the surface energy activated during crack lips motion. The displacement jumps are d ¼ ðd  n; d  tÞ with ðn; tÞ the crack normal and tangential unit vectors (cf. Fig. 6). The dissipative process is governed by the threshold variable j which gives the maximum norm of the jumps ever reached during the crack history. The surface energy is then expressed as:

Wðd; jÞ ¼ Hðkdk  jÞWdis ðkdkÞ þ ½1  Hðkdk  jÞWlin ðkdk; jÞ þ IRþ ðd:nÞ with H the Heavyside function (HðxÞ ¼ 1 if x P 0 and HðxÞ ¼ 0 if x < 0), and I an indicator function, preventing the crack lips interpenetrating each other (IRþ ðd  nÞ ¼ þ1 if d  n < 0 and IRþ ðd  nÞ ¼ 0 otherwise). Thus, after the value of kdk  j, the surface energy exhibits either a dissipative mode or a non dissipative mode that will be denoted Wdis or Wlin respectively:

Wdis ðkdkÞ ¼

Z

wdis ðkdkÞ dC and Wlin ðkdk; jÞ ¼

C

Z

wlin ðkdk; jÞ dC;

C

with Wdis and Wlin the surface energy densities for each mode. The first mode kdk > j deals with the damage activation of the crack. The bulk energy is partially transformed into surface energy in an irreversible process, which will be designed by ‘‘dissipative mode’’. The second mode kdk < j corresponds to displacement jumps without energy dissipation (for instance crack closure or re-opening without cohesive zone growth). This reversible process, with a stationary threshold variable, is referred to as ‘‘linear mode’’. Let’s give the expression of surface energy density in each case. 3.1.1. Surface energy density for the dissipative mode In case of cohesive zone growth kdk P j, the surface energy density is defined as:

r2c

kdk2 þ rc kdk 4Gc    rc Exponential law wdis ðkdkÞ ¼ Gc 1  exp  kdk : Gc

Polynomial law wdis ðkdkÞ ¼ 

ð18Þ

Let us note that wdis is not differentiable with respect to the displacement jump d in zero. To overcome this difficulty, a new parameter j0 (a priori known) is introduced to regularize the energy. In the subset ½0; j0  the energy is defined as a quadratic function with a C 1 connection in j0 (see Fig. 7). Note that wdis is not exactly the dissipated energy. For instance, for kdk = j0, wdis – 0. The terms dissipative and linear modes should be understood as referring to the respective cases kdk P j and kdk < j, and the energy for kdk = j0 should be continuous and derivable. The indicator function which prevents the lips interpenetration is numerically estimated by a continuous penalty function denoted wpen and defined as:

IRþ ðd  nÞ  wpen ðd  nÞ ¼

1 Pðj0 Þhd  ni2 þ cst 2

with P defined by (21), h i ¼ minð; 0Þ, and cst a constant ensuring the continuity of Wpen for kdk = 0.

Fig. 6. Crack surface and cohesive zone.

ð19Þ

10

G. Debruyne et al. / Engineering Fracture Mechanics 90 (2012) 1–29

Fig. 7. Surface energy density vs. the norm of the jump displacement for the polynomial (left) and exponential (right) laws.

3.1.2. Surface energy density for the linear mode In the case of a cohesive zone evolution without dissipated energy (kdk < j), a quadratic surface energy density depending on the jump displacement norm is selected:

wlin ðkdk; jÞ ¼

1 PðjÞkdk2 þ cst; 2

ð20Þ

where the expression of P depends on the cohesive law:

 rc  j 2Gc   rc rc exp  j Exponential law PðjÞ ¼ j Gc

Polynomial law

PðjÞ ¼ rc



1

ð21Þ

the expression of P is chosen to ensure the continuity of the derivative of W with respect to j, and cst is the same constant as in (19), adjusted to ensure the continuity and the derivability of W at kdk = j. The material parameters are the critical stress rc and the fracture energy Gc, which are related to the toughness by the Irwin relation Gc ¼ K 2Ic =E0 or Gc ¼ K 2IIc =E0 , with K Ic ; K IIc , the mode I and mode II toughness, and E0 ¼ Eð1  m2 Þ for plane strain conditions and E0 ¼ E for plane stress conditions, where E is the material Young modulus. Note that the quantities KI, KII are irrelevant in a cohesive model, since the solution is not singular at the crack tip, but the toughness parameters KIc, KIIc are kept here, because they are widely used in the area of material science. 3.1.3. Cohesive stress vector The cohesive stress vector denoted ~ r corresponds to the energy derivative with respect to the displacement jump d. Its depends both on the cohesive law and the energy mode. The general expression of the stress vector is:

~ r ¼ Hðkdk  jÞ~ rdis þ ½1  Hðkdk  jÞ~ rlin þ ~ rpen ;

ð22Þ

with the following definitions. 3.1.4. Stress vector in dissipative mode

  @wdis 1 rc ;  ¼ rc d kdk 2Gc @d   @wdis d rc exp  kdk : ¼ ¼ rc kdk @d Gc

Polynomial law ~ rdis ¼ Exponential law ~ rdis

ð23Þ

3.1.5. Stress vector in linear mode

~ rlin ¼

@wlin ¼ PðjÞd @d

ð24Þ

with P defined by (21). 3.1.6. Penalization Stress vector

~ rpen ¼

  @wpen hd  ni ¼ Pðj0 Þ : @d 0

ð25Þ

G. Debruyne et al. / Engineering Fracture Mechanics 90 (2012) 1–29

11

Fig. 8. Norm of the stress vector versus norm of the displacement jump for the polynomial (left) and exponential (right) laws.

Fig. 8 illustrates the stress vector norm evolution versus the displacement jump for both cohesive laws. The arrows point at stress vector paths related to irreversible or reversible energy (dissipative or linear mode), which are associated with the current jump value kdk (regarding the maximum jump j). The initial slope, due to the regularization of the law, is governed by the parameter j0 . In other words, the initial value of j is equal to j0 . 3.2. Interface finite element properties The displacement jump is modeled by a four node interface finite element (a linear shape function, and two gauss points are considered). The facing nodes (1–4 and 2–3 in Fig. 9) distance may be equal to zero, since only the displacement jumps are considered degrees of freedom. This interface element will use a local system in the basis (n, t) defined by the normal and the tangent vector of the large side (nodes [1 2] and [3 4] in Fig. 9). Elsewhere the global system will be denoted by ðX; YÞ. The jump displacement at a gauss point dg is defined as a linear function of the node displacements U j :

dg ¼ M g U j ;

ð26Þ

where the matrix M g , defined at each gauss point, is a linear combination of the shape functions. 3.2.1. Total energy minimization The total energy of an elastic body X is defined as the sum of its bulk energy UðUÞ inside the sound material X n C and its surface energy Wðd; jÞ on a discontinuity surface C, the external work W ext ðUÞ being substracted:

ET ðU; d; jÞ ¼ UðU Þ þ Wðd; jÞ  W ext ðU Þ Due to some mathematical difficulties, the problem will be restricted to a local minimum search:

We seek ðU  ; d Þ as a local minimum of ET ðU; d; jÞ:

ð27Þ

3.2.2. Discrete FE problem To simplify the notations, the subscripts designing the node values will be omitted from here onward. Using the relationship (26): d = MU, the minimization problem (27) may be expressed again as:

The search for a local minimum U  of ET ðU; MU; jÞ:

ð28Þ

Let us detail the different steps for solving the problem of an elastic body X submitted to a volume force density f and to a surface force density F on a part CN of the external boundary @ X, while displacements U d are prescribed on the part CD of @ X. The total energy reads:

Fig. 9. Scheme of the interface finite element.

12

G. Debruyne et al. / Engineering Fracture Mechanics 90 (2012) 1–29

ET ðU; MU; jÞ ¼ UðUÞ þ WðMU; jÞ  W ext ðUÞ; with U, the bulk energy defined by:

UðUÞ ¼

1 2

Z

ðBUÞt ABUdX

XnC

where A is the elasticity tensor and B the finite element matrix corresponding to the symmetrical gradient of the shape functions. The stress tensor is defined by r ¼ Ae with the strain tensor e ¼ BU. The energy on a discontinuity surface C reads:

WðMU; jÞ ¼

Z

wðMU; jÞdC; C

with W the surface energy density defined earlier. The external work is given by:

W ext ðUÞ ¼

Z

f :UdX þ

XnC

Z

F:UdCN CN

The necessary equilibrium condition for U takes a local minimum reads:

@ET ðU; MU; jÞ  V P 0; @U

ð29Þ

for any test fields V. Knowing that:

@w @w @d rðUÞ; ¼ ¼ Mt ~ @U @d @U condition (29) becomes:

Z

Bt rðUÞ:VdX þ

XnC

Z C

Mt ~ rðUÞ  VdC 

Z XnC

f :VdX 

Z

F  VdCN P 0

CN

for each test field V. The previous inequality becomes an equality for quite selected V, thus we obtain the equilibrium between internal and external forces:

F int ðUÞ ¼ F ext ; with:

(

R R F int ðUÞ ¼ XnC Bt rðUÞdX þ C M t ~ rðUÞdC R R ext F ¼ XnC f dX þ CN FdCN

Let us denote by C the linear operator corresponding to the dirichlet boundary conditions. The dual formulation of the previous equations leads to the system:

(

F int ðUÞ þ C t k ¼ F ext CU ¼ U d

ð30Þ

The unknown of the problem is then, for each loading step, the couple ðU; kÞ, where k is the lagrange multiplier associated with dirichlet conditions (in fact the reaction forces). To solve the previous non linear system a Newton algorithm is used. Finally, the threshold variable j is updated at each loading step i, as follows:

ji ¼ maxfMU i ; ji1 g

ð31Þ

4. Dynamic crack growth in a Double Cantilever Beam with material interfaces In order to deal with realistic problems, rather than 1D academic tests, a bimaterial DCB of length L and height 2h , with an initial crack length a, will be considered. Although analytical investigations have been already performed for this specimen [14,15], an exact solution seems out of reach. The elastic properties are the same for both materials: E = 200000 MPa, m = 0.3, q = 7500 kg/m3 but they differ in their respective surface energy density G1c and G2c , with a ratio a ¼ G1c =G2c > 1. As already described for the film problem, two situations are considered: the first one where there is a single surface energy discontinuity located at x = x0, beyond which Gc ¼ G2c , and a second case where a small brittle flaw ðGc ¼ G2c Þ is embedded into the tough material ðGc ¼ G1c Þ between two material interfaces located at x = x0 and x = x1 (cf. Fig 10). An initial edge crack is considered on the mid plane. A displacement ±Uy is slowly prescribed at each end-arm of the specimen, so that the crack is quasi-statically growing along the mid plane, until it reaches the material interface at x = x0. At this stage, the displacement pffiffiffiffiffi is prescribed to a constant value so that G ¼ G1c (here G1c ¼ 16; 000 N=m to fit to a toughness K 1Ic  60 MPa m related to a brittle steel). The energy surface jump at the material interface x = x0 involves inertial effects and the crack should run throughout

G. Debruyne et al. / Engineering Fracture Mechanics 90 (2012) 1–29

13

Fig. 10. DCB geometry and loading.

the first material interface in the weakest material with possible stop-and-go growth or definitive crack arrest. In the situation where a small brittle flaw with a low toughness is embedded into the tough material, the further question which arises is whether the crack arrest will occur inside or beyond the flaw at x > x1. The case where the material interface itself may be debonded along the y axis is not considered here. This point has been widely investigated by Arata and Needleman [25]. We focus here on the conditions of crack growth and arrest along the x direction, as essentially a function of the ratio a. Moreover, concerning the simulations using CZM, some analysis has been carried out to evaluate the influence of the cohesive strength and the cohesive law type on the crack kinematics. Provided symmetry of the geometry and the loading, only the upper half beam has been meshed, with 87959 four node quadrilateral isoparametric elements for the bulk material and the analysis is performed with plane strain conditions. The DCB model is close to the film peeling problem. Hence the numerical results are expected to have the same qualitative properties than the film peeling analytical solution two numerical methods are used here: the main one uses the Cohesive Model described above, and, for comparison, a subsidiary one deals with a special version of the node release method, detailed in Section 4.1.2. Concerning the first one, the potential crack path along the mid-line is discretized by 4000 interface elements uniformly positioned ahead of the initial crack (each element length is fixed to 25 lm so that the estimated cohesive zone suggested by Rice for a linear cohesive law [26,27], Lcoh  932p

EG1c ð1m2 Þaðrc Þ2

 a1 2:15 mm is covered by at least four ele-

ments if a < 20, which is the maximum ratio tested here). The small initial crack opening parameter j0, governing the energy regularization and described in the second section, is written here as: j0 ¼ e rGcc , where e is an arbitrary small adimensional factor. This latter has been successfully tested in the range [105, 103] without unwanted side-effects results. For e < 105, the numerical convergence becomes tough, and for e > 103, some divergent results arise. Therefore, in all the following investigations, e =104. 4.1. The bimaterial Double Cantilever Beam with one material interface For this configuration, some parametric studies have been carried out, to test the sensitivity to different numerical crack growth procedures, various cohesive parameters such as variation in the cohesive strength or the cohesive law type. The effects of the stress waves and the bulk material plasticity have also been examined. 4.1.1. A simplified CZM process The investigation of the crack growth with CZM, step by step, from the quasi-static propagation in the first material to the dynamic running in the second one, is the most straightforward method. From now on, this latter will be called the natural process. Nevertheless, it entails the following drawback: the time orders of magnitude in the two different propagation phases are so different that it is difficult to catch with certainty the very moment where the crack tip reaches the material interface, and thus sharp running begins. So many time steps have to be used to catch the sharp crack onset. This is a very time consuming process, and it has been only applied to a very specific case detailed at the end of the section, to validate a simplified procedure in which the first stage is disregarded (this ‘‘simplified process’’ will be used in most of the present investigations). So as not to perform the quasi-static phase analysis, which is not a crucial element, the initial crack tip at time t = 0 will be considered at the exact position x = x0. The load level should involve G ¼ G1c , and it will be analytically 3 determined, considering a Bernoulli beam (G ¼ 108EI2 U 2Y =h a4 , where I is the moment of inertia). During this preliminary stage, the interface element displacements are blocked, so that no initial decohesion occurs ahead of the initial crack tip. At this point, the interface element displacements are suddenly released, and the cohesive zone growth is computed with a dynamic procedure within the interval [0, 15104 s] divided in uniform time steps (Dt = 2  107 s will be used most of the time) and using a classical non dissipative Newmark scheme for time integration:

) 2 _ € þ 2bUðt € þ DtÞ Uðt þ DtÞ ¼ UðtÞ þ DtUðtÞ þ D2t ½ð1  2bÞUðtÞ c ¼ 0:5; b ¼ 0:25: _ þ DtÞ ¼ UðtÞ € þ DtÞ _ € þ cUðt Uðt þ Dt½ð1  cÞUðtÞ

14

G. Debruyne et al. / Engineering Fracture Mechanics 90 (2012) 1–29

4.1.2. Comparison of the simplified CZM process and a node release method An alternative method to CZM is based on the node release techniques, which used to be very popular in the seventies and the eighties, even using forces relaxation [28,29], and may still be encountered in more recent works with more complex fracture behavior [30]. The basic discontinuous version is rather rough, but due to the lack of experimental tests for our problem, it still provides some elements of comparison with more sophisticated methods such as cohesive models. Furthermore, it is based on a mere Griffith criterion, so that it is interesting to check that it gives close results to those of a model with a sufficiently small cohesive zone. These two points explain why this method is used here. However, as the usual node release method is unable to predict crack tip velocity, it is associated with an iterative procedure in which, for each crack growth step, the dynamic energy release rate is estimated for several virtual crack velocities. The actual velocity is found by dichotomy, fulfilling the condition G ¼ G2c in the weak material, and the crack length is updated by releasing nodes to reach the next crack growth step. This procedure is relatively sensitive to mesh refinement and a large number of iterations is necessary to get accurate results. So, this is a very time consuming process, and it is only compared to the Cohesive Model for a relatively short time interval, with the aim of assessing CZM ability to perform dynamic crack growth analysis, and to evaluate the influence of the cohesive zone length (governed by the cohesive strength) on the crack behavior in comparison with a pure Griffith criterion. Figs. 11 and 12 display the comparative crack velocity, derived from a Griffith criterion associated to the releasing nodes iterative process, and from the polynomial CZM simplified process (with a cohesive strength r2c ¼ 1200 MPa). The results are in good accordance, except for the first point where the node release procedure exhibits an irrelevant velocity above the Rayleigh wave speed (CR = 3460 m/s), probably because the complete node release is instantaneous, without force relaxation. Moreover, the CZM accounts for crack tip stop-and-go growth, which the node release model does not (see Fig. 12). 4.1.3. Parametric analysis using CZM with different surface energy ratios This section is devoted to an overview of the crack propagation and arrest characteristics, for a large range of surface energy ratios from a = 1.33 to a = 20 simulated by CZM. Different configurations were tested (the cohesive law type, the cohesive stress values, the comparison with the node release method, the crack arrest report, and the considered procedure (natural or simplified process)). The results are summarized in Table 1. Fig. 13 shows an overview of the crack velocity for each selected surface energy ratio, obtained by a polynomial simplified CZM analysis. The same qualitative aspects may be observed for all ratios, that is to say a decreasing crack velocity followed by several ‘‘stop-and-go’’ stages. The higher the ratio, the faster the crack growth, however the initial crack tip speed (with respect to the Rayleigh wave velocity) is lower than the peeling point velocity observed for the film peeling problem in the first section ð‘_ ¼ aa1 cÞ. þ1 4.1.4. Parametric analysis using CZM with different cohesive parameters In this section, we focus on the influence of the cohesive strength and of the traction–separation law (polynomial or exponential as described in second section) on the crack kinematics. Crack and cohesive zone growth versus time, were plotted in Fig. 14, for a = 10. It can be seen in this figure that the cohesive length is keeping a constant value (of about 0.15 mm) during the propagation. This value obeys the Rice formula, so that the location of the crack tip and the cohesive tip coincides in Fig. 14. Comparing two different cohesive strength values for the weak material (rc = 1200 MPa, rc = 600 MPa), only a very slight difference can be observed on the resulting crack growth,

7000 Crack velocity α=10, σc=1200 MPa

6000

Velocity (m/s)

5000 CZM Nodes releasing

4000 3000 2000 1000 0 0

0,005

0,01

0,015

0,02

0,025

0,03

0,035

Δa (m) Fig. 11. Crack velocity versus crack growth using linear CZM or a node releasing procedure for a = 10, and

r2c ¼ 1200 MPa.

15

G. Debruyne et al. / Engineering Fracture Mechanics 90 (2012) 1–29

4500

Crack velocity α =3, σ c=1200 MPa

4000

Velocity (m/s)

3500

CZM Nodes releasing

3000 2500 2000 1500 1000 500 0 0

0,001

0,002

0,003

0,004

0,005

0,006

Δ a (m) Fig. 12. Crack velocity versus crack growth using linear CZM or a node releasing procedure for a = 3, and

r2c ¼ 1200 MPa.

Table 1 Outline of investigations with CZM.

a ¼ G1c =G2c

Polynomial CZM (simplified process)

Exponential CZM (simplified process)

Node release procedure

Crack arrest detection

Polynomial CZM (natural process)

a ¼ 1:33

rc = 1200 MPa

Not tested

Not tested

Not tested

a=2

rc = 1200 MPa

Tested rc = 1200 MPa

Tested

Crack arrest detected Crack arrest detected

a=3

rc = 1200 MPa

Not tested

Tested

a=5

rc = 1200 MPa

Not tested

Not tested

a = 10

rc = 1200 MPa, rc = 600 MPa

Not tested

Tested

a = 20

rc = 1200 MPa

Not tested

Not tested

r1c ¼ 2400 MPa r2c ¼ 1200 MPa

Tested

Crack arrest detected No arrest detected No arrest detected No arrest detected

Not tested Not tested Not tested Not tested

2500

VELOCITY (m/s)

2000

alpha=2 alpha=5 alpha=10 alpha=20

1500

1000

500

0 0

0,005

0,01

0,015

0,02

0,025

0,03

Δ a (m) Fig. 13. Crack velocity versus crack tip position for different a (polynomial separation law).

which means that, compared to the previous node release analysis, this range of values corresponds to a Griffith criterion. A key point concerns the acceleration and transient arrest at the exact time when the Rayleigh wave comes back to the crack tip (see Fig. 14). Unlike the film problem, the wave strike does not involve a definitive arrest, because the crack, after a short arrest, restarts, then stops and starts again (the analysis has not been pursued to detect a definitive arrest or a complete

16

G. Debruyne et al. / Engineering Fracture Mechanics 90 (2012) 1–29

0,035

First arrival of the reflected Rayleigh wave on the crack tip

0,03

Δ a (m)

0,025 0,02

crack growth: cohesive strength 1200 MPa cohesive zone growth: cohesive strength 1200 MPa crack growth: cohesive strength 600 MPa

0,015 0,01 0,005 0 0,00E+00

2,00E-05

4,00E-05

6,00E-05

8,00E-05

1,00E-04

1,20E-04

Δ t (s) Fig. 14. Crack growth versus time using linear CZM for a = 10.

fracture). Because of pressure and shear wave reflections on the longitudinal boundaries, the problem is much more complex than that of the one dimensional film case, though some similarities may be noticed. The next paragraph will be devoted to giving further details about the wave reflections effects. Henceforth, the analysis is focused on some lower ratios, according to the usual characteristics of industrial metal defects, as depicted in the last section (2 6 a 6 3).Fig. 15 is similar to Fig. 14 except that a = 2, and shows the comparative predictions for each cohesive model (it is not possible to describe a real crack surface with the exponential model because the cohesive energy never entirely vanished in that case, so that only the cohesive zone growth is plotted in Fig. 15). However, the two models are in good accordance. As for the other ratios, the crack stops and restarts many times, and it is difficult to know if the last stop is the definitive arrest or not, without carrying on computations for a very long time. As we will see in a next paragraph, an energy balance will be helpful to answer this question. 4.1.5. A special case for the ratio a = 2, using the natural process with different initial fracture energy values G1c The previous case with a = 2 is investigated again but with the new following initial surface energy values: G1c ¼ 32; 000 N=m, and G1c ¼ 64; 000 N=m. The cohesive stress considered in the second material is: r2c ¼ 1200 MPa. This

0,05 0,045 0,04

Crack growth (m)

0,035 0,03 0,025

Polynomial CZM (crack length)

0,02

Exponential CZM (crack+cohesive length)

0,015 0,01 0,005 0 0,0E+00 1,0E-04

2,0E-04

3,0E-04

4,0E-04

5,0E-04

6,0E-04

7,0E-04

8,0E-04

Time (s) Fig. 15. Crack growth versus time using polynomial and exponential CZM for a ¼ 2.

9,0E-04

17

G. Debruyne et al. / Engineering Fracture Mechanics 90 (2012) 1–29

analysis is pursued with the natural CZM process, including the quasi-static propagation in the first material. Contrary to the simplified process, a cohesive stress has to be introduced for the first material as well. It will be r1c ¼ 2400 MPa, and the separation law used here will be the polynomial one. The results for the two initial energy values are first compared with each other, and then compared with the previous simplified CZM process solutions. Fig. 16a shows the crack growth in the second material for both cases, starting from the crack onset. This one does not correspond to the same critical load for the different initial surface energy values. Therefore, it is difficult to adjust the two trials starting time. Thus, there is a gap of less than one millisecond between the two curves depicted in Fig. 16a. Nevertheless, the different stages of the crack kinematics are very similar, and the curves are almost parallel. The most outstanding feature is that the kinematics, and in particular the arrest, is not very sensitive to the initial fracture energy value, and mainly depends on the ratio a. This statement was already noticed for the 1D film problem. Three stages are spotted: the first one is about quasi-static, then a stage with a moderate crack velocity arises (the velocity increases when the initial surface energy decreases). Finally, after about a 2.5 mm crack extension, a fast growth stage appears with very similar aspects for both trials. Fig. 16b shows the dynamic third crack growth stage with the natural process, in addition to the crack growth computed with the simplified process. The surface energy of the first material is respectively G1c ¼ 16; 000 N=m for the simplified process and G1c ¼ 64; 000 N=m for the natural one, but with the same ratio a = 2. It is little bit tricky to fit the starting point of each curve, because in the simplified process, the cohesive forces are suddenly released, while they are progressively activated in the natural process. Therefore, the initial velocities corresponding to the slope of the curves in Fig. 16b cannot be the

Fig. 16a. Crack growth versus time using natural CZM process for a ¼ 2, with two initial surface energy values.

Natural process versus simplified process

0,05 0,045

crack growth (m)

0,04 0,035 0,03

Natural process Gc1=64000 N/m Simplified process Gc1=16000 N/m

0,025 0,02 0,015 0,01 0,005 0 0

0,0001

0,0002

0,0003

0,0004

0,0005

0,0006

0,0007

0,0008

0,0009

time (s) Fig. 16b. Crack growth versus time using natural or simplified CZM process for a ¼ 2, with two different initial surface energy values.

18

G. Debruyne et al. / Engineering Fracture Mechanics 90 (2012) 1–29

same for both processes. Furthermore, space and time are not exactly discretized in the same way for the finite element model related to each process. Indeed, even if curves are not perfectly matching, the major parameter of influence to predict the arrest remains the energy ratio a. The absolute value of the surface energy in the first material, and the initial crack velocity in the second one, may be regarded as phenomena of a second order. 4.1.6. Time discretization and wave travel effects A set of time steps in the range [108, 106]seconds has been tested to simulate the early propagation stage with the ratio a = 3 (until t = 6  106 s, so that the time scale of interest here is much shorter than the time interval investigated in the previous paragraph t = 6  106 s). A first observation (Fig. 17) is the short crack onset delay (about 3  107 s), caused by the cohesive zone process before complete fracture. Then, the graph (Fig. 17) showing crack growth versus time exhibits some oscillations with a period of about 2  106 s which seem to correspond to the back and forth propagation of pressure waves reflected from the beam horizontal edges. These waves enhance the crack growth, when hitting the crack tip. The crack slows down just before a new wave strikes it, and so on. A very short time step should be used to capture these oscillations (at least 2  107 s, which seems to be a good compromise between numerical accuracy and computational time). To accurately examine the influence of the stress waves on the crack growth, a modified DCB test, with absorbing boundary conditions on the beam edges is performed, to prevent any wave reflection. These specific boundary conditions are based on paraxial approximations of the wave equations [31]. A specific boundary condition resulting of the wave impedance expansion, and similar to some viscous dampers, cancels out the reflected wave, whatever the incidence angle and the frequency. In Fig. 18, crack growth is plotted versus time for both models. Until about t = 2  106 s, the kinematics for both models coincide, but just after this time, the models more and more differ from each other. The model with the absorbing conditions tends to slow down the crack kinematics. This result is in good agreement with the computational results obtained by Kanninen who compared the crack propagation inside a DCB specimen with a crack growth in an infinite medium [32,33]. In Fig. 19 this phenomenon is stressed out from an energetic point of view, by showing the kinetic energy of the whole body during the crack growth for both models. The divergence on energy between the two different models appears at time t  1  106 s, which is about the time when the pressure wave hits the horizontal edge of the DCB specimen (at this time the crack propagation is not perturbed by wave reflection yet, as shown in Fig. 17). The double conclusion which may be drawn from this paragraph is that it is vital to take wave propagation in account is essential to avoid the risk of underestimating crack arrest, but this simulation needs a tight time discretization involving a significant increase in computational time. This raises the following issue: would a simpler and cheaper method, limited to the quasi-static frame, be able to predict a first approximation of the crack arrest? (even with a reasonable overestimation of the cracked length). The next paragraph suggests some tracks based on energy considerations (which are close to those that has been used for the peeling test [4]). 4.1.7. Energy balance considerations: a tool to predict the crack arrest It has been proved for the film peeling test [4], that the dynamic solution for the debonding arrest converges to a quasistatic one, which is governed by a principle of conservation of the total quasi-static energy, and by the method presented in

4,0E-3

Pressure wave Round-trips

3,5E-3

Crack growth (m)

3,0E-3 2,5E-3 2,0E-3 0,01 microsecond step(12h CPU time)

1,5E-3

0,02 microsecond step (9h CPU time) 0,1 microsecond step (3h CPU time)

1,0E-3

0,2 microsecond step (2h CPU time) 1 microsecond step(1h CPU time)

500,0E-6 000,0E+0 000,0E+0

1,0E-6

2,0E-6

3,0E-6

4,0E-6

5,0E-6

6,0E-6

7,0E-6

time (s) Fig. 17. Various time discretizations and their ability to catch wave reflection effects.

8,0E-6

19

G. Debruyne et al. / Engineering Fracture Mechanics 90 (2012) 1–29

3,5E-3

crack growth (m)

3,0E-3 2,5E-3 2,0E-3 1,5E-3 0.1 microsecond step, 3h CPU time

1,0E-3

0.1 microsecond step absorbing boundaries 8h CPU time

500,0E-6 000,0E+0 000,0E+0

1,0E-6

2,0E-6

3,0E-6

4,0E-6

5,0E-6

6,0E-6

7,0E-6

time (s) Fig. 18. Comparative crack growth for the actual DCB and the modified model with absorbing boundary conditions.

3,5 3

Energy (J)

2,5 2 Kinetic energy actual boudaries

1,5

Kinetic energy absorbing boundaries

1 0,5 0 0,E+00

1,E-06

2,E-06

3,E-06

4,E-06

5,E-06

6,E-06

7,E-06

time (s) Fig. 19. Comparative kinetic energy evolution for the actual DCB and the modified model with absorbing boundary conditions.

Section 2.1.10. It is suggested to transpose this principle to the DCB specimen, assuming that the kinetic energy may be neglected after the crack arrest. Strain, kinetic and total energy (including the fracture energy dissipated) are plotted against propagation time in Fig. 20, for a DCB with a surface energy ratio a = 2. Though these curves may seem quite similar to those of the film peeling at first, with common kinetic energy growth and strain energy decrease, followed by stages with a kinetic energy decrease and a reduced elastic energy decline, several differences may be noticed. The first Raleigh wave hit on the crack does not correspond to a cancellation of the kinetic energy, and unlike the one dimensional film problem, many other wave reflections disturb the mechanical state, as already described in the previous paragraph. At the end of the computation, the kinetic energy is about to vanish, while the crack growth seems to be steady around Da = 43 mm (see Fig. 15). In fact, it does not totally cancel out, but bounces a little after the crack arrest (this bouncing part of kinetic energy corresponds to the beam free vibrations). Another analysis with a = 3 shows (Fig. 21) that the kinetic energy also vanishes, while the crack seems to stop at about Da = 78 mm, but then bounces a little bit more than in the previous case with a small crack restart until arrest at Da = 79.5 mm. 4.1.8. Quasi-static prediction Since there is no inertial force at the beginning of the process and only little kinetic energy at the end, it may be suggested that a static energy balance is valid to detect the crack arrest. This energy balance states that the elastic energy plus the fracture energy at the crack arrest time is equal to the elastic energy at crack growth onset time. This is the same process as described in Section 2.1.10, concerning the film peeling (see Fig. 4), and some similar to an energy balance, applied to laboratory tests with a small crack, for predicting crack arrest [35]. Using the Bernoulli beam theory, the energies involved in this static method can be assessed in an analytical way.

20

G. Debruyne et al. / Engineering Fracture Mechanics 90 (2012) 1–29

3,00E+02

2,50E+02

Energy (Joule)

2,00E+02 Total energy Kinetic energy Elastic energy

1,50E+02

1,00E+02

5,00E+01

0,00E+00 0,00E+00 1,00E-04 2,00E-04 3,00E-04 4,00E-04 5,00E-04 6,00E-04 7,00E-04 8,00E-04 9,00E-04

time (s)

70

0,09

60

0,08

Kinetic energy (J)

0,07 50

0,06

kinetic energy (J) crack growth (m)

40

0,05 0,04

30

0,03

20

Crack growth (m)

Fig. 20. Energy balance versus time using polynomial CZM for a = 2.

0,02 10 0 0,00E+00

0,01 0 5,00E-04

1,00E-03

1,50E-03

2,00E-03

Time (s) Fig. 21. Kinetic energy and crack growth versus time using polynomial CZM for a = 3.

The elastic energy is Eelas ¼

Etot ðaÞ ¼

108EI 3

2

U 2Y

3h a3

108EI2 U 2Y 3h3 a3

and thus the total energy relation holds:

þ G2C Da

ð32Þ

The analytical relation (32) is plotted in Fig. 22 along with the kinetic energy predicted by the numerical dynamic CZM analysis, with a narrow range of ratios from a = 1.33 to a = 3, which represents usual ratios related to brittle flaws in industrial components. It is first emphasized that the total energy minimum (at which the Griffith criterion G ¼ G2C holds) is located at respectively Da = 6.5 mm, Da = 18 mm and Da = 31 mm for a = 1.33, a = 2, and a = 3, that is far below the dynamic crack arrest (respectively Da = 16 mm, Da = 43 mm and Da = 79.5 mm), which confirms the lack of safety margin of usual static arrest predictions using the Griffith criterion associated to equilibrium conditions (which was already noticed for the film peeling). In return, the balance energy gives a total crack growth of Da = 14 mm, Da = 43 mm and Da = 82 mm which is pretty close to the dynamic predictions. 4.2. Double Cantilever Beam featuring a brittle flaw embedded in the host material This section is devoted to the case of a 3 mm wide brittle inclusion embedded into the host metal (see Fig. 10). This inclusion has a surface energy Gc ¼ G2c , with a ¼ G1c =G2c ¼ 3 (which is a typical order of magnitude for brittle flaws in industrial

21

G. Debruyne et al. / Engineering Fracture Mechanics 90 (2012) 1–29

300

250

« Quasi-static » arrest Prediction by energy balance

200

Energy (j)

total static energy alpha=1.33 total static energy alpha=2 total static energy alpha=3

150

Quasi-static Griffith criterion (G=Gc2)

kinetic energy a=1.33 kinetic energy alpha=2 kinetic energy alpha=3

100

50

0 -0,01

0

0,01

0,02

0,03

0,04

0,05

0,06

0,07

0,08

0,09

Δ a(m) Fig. 22. Kinetic energy versus crack growth using polynomial CZM for a ¼ 1:33; a ¼ 2; a ¼ 3, compared to an analytical total energy balance for a quasistatic process.

components). The same geometry, loading, and discretization with a polynomial interface law as in the previous section are considered. Furthermore, the safe metal is considered either elastic or elasto-plastic (the small inclusion always remains elastic). In the safe metal, the yield stress is set to rY = 600 MPa and the hardening modulus to ET = 20,000 MPa. 4.2.1. Comment The minimization problem (see relation (27)), still holds for the deformation theory of plasticity (the bulk energy U is then considered as a non linear elastic potential). But when some unloading occurs in the material surrounding the cohesive elements (which is likely for some areas during the propagation), this theory is no longer valid. Thus, it is not guaranteed that a minimization problem is the appropriate process to seek the solution. Nevertheless, our numerical model still works, but without the same theoretical convergence guarantee as in elasticity. One of the question which is raised is whether the crack abruptly stops or keeps growing at the second flaw material interface, after the high velocity growing stage through the brittle flaw. We would also like to assess the influence of the plasticity on crack kinematics. Fig. 23 shows the crack tip kinematics throughout the flaw for an elastic or elasto-plastic safe material. For both behaviors, the crack runs through the flaw and the arrest occurs just on the second flaw material interface (x = x1), the only difference which emerges between the two behaviors is the more stable crack kinematics in the flaw, when the host metal (x < x0 designed as the starter material and x > x1 designed as the target material) is elasto-plastic. Hence, the crack kinematics within a given elastic media is damped, when the starting and the target materials are elasto-plastic. This damping is due to the plastic energy dissipated in the host metal. This point is emphasized in Fig. 24, by comparing the kinetic energy developed in the beam, for both host metal behaviors. It is worth pointing out that the amount of energy is smaller for the plastic behavior than it is for the elastic one, and that there is a drop of periodic energy enhancement. Kinetic energy still increases once the crack has stopped, especially in the elastic case, where this energy soars after the arrest, due to the arrest shock. Thus, the quasi-static prediction, based on the principle of total energy conservation, that neglects inertial effects after the crack arrest, is no longer available here, as it has already been mentioned in Section 2.2 in the case of the peeling test with an inclusion. This inability is illustrated in Fig. 25, where the total energy, analytically computed with a quasi-static assumption, is plotted against the virtual crack growth. The energy balance is verified for Da  10 mm which is far for the dynamic prediction Da = 3 mm. 5. Subclad flaw growth and arrest in a Pressure Water Reactor vessel shell submitted to a pressure loading This section is devoted to some investigations about the crack growth and arrest in a Pressure Water Reactor vessel under some particular geometry and loading conditions. Many authors have addressed this issue with thermal-shock transients,

22

G. Debruyne et al. / Engineering Fracture Mechanics 90 (2012) 1–29

Crack growth through a brittle inclusion 0,0035

Elastic or elastoplastic host metal Elastic inclusion

0,0025 0,002

Elastic host metal Elastoplastic host metal

0,0015 0,001 0,0005 0 0,00E+00

4,00E-06

8,00E-06

1,20E-05

1,60E-05

2,00E-05

time (s) Fig. 23. Crack tip growth through the brittle flaw.

7 kinetic energy evolution

kinetic energy (J)

6 5 4 Elastic host metal

3

Elastoplastic host metal

2 1 0 0

0,0005

0,001

0,0015

0,002

0,0025

0,003

0,0035

crack growth (m) Fig. 24. Comparative kinetic energies for elastic or elasto-plastic material.

Total energy (quasi-static evaluation)

266 264 262

Total energy (J)

crack growth (m)

0,003

Inclusion

260

Target

258 256 254 252 250 248 0

0,002

0,004

0,006

0,008

0,01

crack length Fig. 25. Total energy evaluated with a quasi-static assumption, for the elastic material.

0,012

G. Debruyne et al. / Engineering Fracture Mechanics 90 (2012) 1–29

23

invoking the conservatism or the unconservatism of the quasi-static approach [29,34], and the role of crack tip plasticity [36]. This is a wide subject, and we will only focus on a simplified geometry and loading here, regardless of thermal loading or 3D effects. The main goal is to enlighten the part of dynamics, material interfaces and plasticity about the crack arrest, following the previous investigations concerning simple specimens. The vessel 2D geometry is only submitted to an inner pressure loading P increasing up to P at the crack onset. The vessel base metal is a ferritic steel and the inner cladding is an austenitic steel, known to be more ductile. A 6 mm deep radial flaw preexists in the base metal, under the cladding. This loading and this geometry involve more unstable crack growth than the DCB case explored in the previous section. Therefore, the crack growth stability will be the major point examined here. In actual reactor vessels, the base metal gets tougher as we go deeper within the ring wall, because of the decrease of irradiation (and also the rise in temperature inside the wall, which is neglected here). Thus, a crack arrest is expected inside the wall. Two situations are considered here: (1) an homogeneous base metal with an uniform fracture energy is considered, so that no arrest is expected, but the influence of the medium behavior (elastic or elasto-plastic) on stability is analyzed and (2) a thin elastic and brittle inclusion is embedded in the base metal, ahead of the initial flaw, which could enhance the crack velocity. Beyond this zone, the metal exhibits a high toughness again, and the crack arrest conditions are analyzed as function of the toughness value. To take into account the close end effect, the vessel shell ring geometry is modelized by an axisymmetric tore with a large curvature radius (cf. Fig. 26), which involves that the initial flaw is assumed to be axisymmetric too. This latter assumption should lead to a crack growth overestimation (the real flaws are rather of elliptical shape with a finite circumferential width). Each material behavior is considered either elastic or elasto-plastic (except for the inclusion which is always assumed to be elastic). Mechanical computations are carried out with the Finite Element Software Code_Aster [37]. Isoparametric linear triangular and rectangular elements are used to mesh the vessel (cf. Fig. 27). Due to loading and geometry symmetry, only the upper half part of the shell is considered, so that only mode I Fracture is activated. The elastic–plastic and fracture parameters for both materials are shown in Table 2a and ab. Interface elements with the polynomial cohesive law, introduced in Section 2, are disposed ahead of both tips of the initial crack, along the expected rectilinear crack path. To lighten the Finite Element model, the cohesive elements are only set along half of the wall width, and the length of each element is in the range [0.05 mm, 0.2 mm] which appears to be small enough, where compared to the previous cantilever beam experiments carried out with similar materials. The fracture energy associated to the cohesive elements is linked to the toughness, through the Irwin relationship:

Gc ¼ g

K 2Ic ; E

ð33Þ

with g = 1 for plane stress conditions, and g = (1  m2) for plane strain conditions. Note that the cladding is much more ductile than the base metal (cf. Tables 2a and b where K Base < K Cladding ). Thus, it is expected that the crack should run through the Ic Ic ferritic steel and not through the cladding. Nevertheless, cohesive elements have been displayed on the side of the cladding because cohesive energy should be dissipated in this zone. 5.1. Flaw onset predictions Our first study is devoted to the crack onset, which can be predicted with a mere quasi-static analysis since the pressure loading is slowly prescribed. Both behaviors (elastic and plastic) are considered for the base metal and the cladding, and

Fig. 26. Schematic geometry for the pressure vessel shell.

24

G. Debruyne et al. / Engineering Fracture Mechanics 90 (2012) 1–29

7.5mm Cladding

Base metal

Brittle inclusion initial 6mm flaw

Cohesive elements Fig. 27. Sketch of the mesh in the vicinity of the flaw.

Table 2a Base metal properties. Fracture, Gc ðkJ=m2 Þ

Toughness, pffiffiffiffiffi K Ic ðMPa mÞ

Critical stress, rc (MPa)

Young modulus, E (GPa)

Plastic yield stress, rY (MPa)

Hardening modulus, ET (MPa)

7

40

1500

205

583

3577

Table 2b Cladding properties. Fracture energy, ½Gc ðkJ=m2 Þ

Toughness, pffiffiffiffiffi K Ic ðMPa mÞ

Critical stress, rc (MPa)

Young modulus, E (GPa)

Plastic yield stress, rY (MPa)

Hardening modulus, ET (GPa)

87

140

1000

205

370

1525

Table 3a Crack onset comparative predictions with a classical Griffith criterion and with cohesive elements for a linear elastic material. Elasticity

6 mm subclad flaw

16 mm subclad flaw

Critical pressure P at the crack onset with CZM pffiffiffiffiffi Apparent toughness for a Griffith crack with P (actual K Ic ¼ 40 MPa m)

P=38.2 MPa pffiffiffiffiffi K Ic ¼ 37:4 MPa m

P=23.4 MPa pffiffiffiffiffi K Ic ¼ 40:8 MPa m

Table 3b Crack onset comparative predictions with a classical Griffith criterion and with CZM for an elasto-plastic material. Plasticity

6 mm subclad flaw

16 mm subclad flaw

Critical pressure P at the crack onset with CZM pffiffiffiffiffi Apparent toughness for a Griffith crack with P (actual K Ic ¼ 40 MPa m)

P=37.4 MPa pffiffiffiffiffi K Ic ¼ 45:4 MPa m

P=24.7 MPa pffiffiffiffiffi K Ic ¼ 44:5 MPa m

compared with a classical analysis based on Griffith criterion. In addition to the actual short flaw of 6 mm, another deep crack with a 16 mm length has been considered. Tables 3a and b show the comparative results for both methods and both cracks, for elastic and elasto-plastic materials. For the Cohesive Model investigation, the critical inner pressure P is associated to the rupture of the first Gauss point, located ahead of the crack tip in the base metal. This pressure value is then used to assess the energy release rate, using a classical Griffith crack model (the cohesive elements around the initial flaw are then desactivated and replaced by a crack ligament). This energy release rate is estimated by a Gh method [38], so that Gh (P), for the critical pressure, is expected to be close to Gc. For the elastic behavior, the critical pressure obtained with the CZM, leads pffiffiffiffiffi to an apparent toughness close to the genuine base metal toughness ðK Ic ¼ 40 MPa mÞ and which gets closer as the crack length increases. This point is well explained by the small cohesive zone length which makes the cohesive model close to the Griffith one.

25

G. Debruyne et al. / Engineering Fracture Mechanics 90 (2012) 1–29

Pressure (Pa)

Crack growth analysis 5,0E+07 4,5E+07 4,0E+07 3,5E+07 3,0E+07 2,5E+07 2,0E+07 1,5E+07 1,0E+07 5,0E+06 0,0E+00 0,0E+00 2,0E-03

Plasticity Yield stress : 583Mpa Elasticity Plasticity Yield stress : 450Mpa

4,0E-03

6,0E-03

8,0E-03

1,0E-02

1,2E-02

Crack length (m) Fig. 28. Quasi-static crack growth analysis, according three constitutive laws.

The elasto-plastic behavior leads to a poor agreement between the two methods as shown in Table 3b, but it is well known that the crack resistance for an elasto-plastic media may be underestimated, when using the concept of energy release rate.

5.2. Crack growth into the base metal, with an homogeneous toughness, using the Cohesive Model First of all, a quasi-static analysis with a loading path following algorithm [39] to match with unstable propagation stages, is performed, in the frame of elasticity and plasticity. A first result is that the propagation always occurs towards the base metal owing to the cladding ductility. For an elastic behavior, the crack kinematics is unstable all along the crack path as expected, which is shown by the uniformly decreasing loading path following the green2 curve in Fig. 28. On the other hand, the blue curve, derived from the elasto-plastic analysis considering usual plasticity characteristics (Table 3a), exhibits a first stable phase followed by a strong unstable propagation. If the base metal ductility is increased, by decreasing the yield stress, a stronger stable phase is first observed and then a quasi-steady growth is obtained, as described by the red line. Beyond these preliminary conclusions, quasi-static investigations are inadequate to predict a potential crack arrest. To confirm the previous conclusions about crack growth stability as a function of material behavior, further dynamic investigations are mandatory. Using the same methodology as in the previous section for the beam analysis, a quasi-static simulation is achieved, up to the critical loading P corresponding to the crack onset (which has different values for an elastic or an elasto-plastic material, as mentioned in Tables 3a and b). At this stage, the dynamic algorithm is triggered, to observe the running crack phase, while the pressure P is frozen. As expected, accordingly to the quasi-static observations summarized in Fig. 28, no propagation is detected in the case of an homogeneous plastic material (so that the pressure value P should exceed P and keep growing afterward, to involve further growth, which characterizes a stable crack growth). On the other hand, for an elastic material, a fast crack growth without arrest and with an increasing velocity is observed, as shown by the red curve of Fig. 29, and so the crack growth is fully unstable. To go even further and to promote a full propagation even with plasticity with the constant pressure P, an elastic inclusion of 2 mm width and 20 mm high, embedded into the base metal ahead of the initial crack tip, is now considered (see Fig. 27). This inclusion is elastic but keeps the same fracture properties as the base metal (Table 2a), and the surrounding material keeps its plastic properties. Under these conditions, the crack quickly runs into the inclusion and stops abruptly in the middle of inclusion (Da = 1.2 mm, see blue curve of Fig. 29). This configuration shows the importance of the bulk material properties on the crack growth stability, and more precisely the stabilizing properties of plasticity (as already been mentioned in Section 4.2 related to the DCB with material interfaces).

5.3. Crack growth into the heterogeneous base metal 1 This elastic inclusion is now embrittled. Its toughness is half the one of the regular base metal ðK incl Ic ¼ K Ic =2 ¼ pffiffiffiffiffi 20 MPa mÞ, which enhances crack growth instability. Beyond this brittle defect, the media is considered now tougher than the initial base metal to take into account the lack of neutron embrittlement inside the wall (so that along the crack propagation, three sections are considered: the regular base metal, with a surface energy G1c , the inclusion with Gincl ¼ 0:25G1c , and c 2 1 a tougher section with Gc > Gc ).

2

For interpretation of color in Figs. 1–8 and 11–32, the reader is referred to the web version of this article.

26

G. Debruyne et al. / Engineering Fracture Mechanics 90 (2012) 1–29

Crack velocity for homogeneous fracture properties in the base metal 2500

crack velocity (m/s)

2000

1500 Elastoplastic material with an elastic 2mm width inclusion Elastic material

1000

500

0 0

0,002

0,004

0,006

0,008

0,01

0,012

crack growth Fig. 29. Dynamic crack growth in the base metal with homogeneous fracture properties.

2500

Crack velocity (m/s)

2000

Gc2=1.59Gc1 1500

Gc2=1.6Gc1

1000

500

0 0

0,001

0,002

0,003

0,004

0,005

0,006

0,007

0,008

0,009

-500

Crack growth (m) Fig. 30. Crack tip velocity for an elastic media with multi fracture properties.

To check arrest conditions, and after some fitting tests, an amount of surface energy around G2c ¼ 1:6G1c has been considered for the tough material section. The base metal, including both sections, is first considered elastic. The crack front velocity through and beyond the inclusion is plotted in Fig. 30. For G2c ¼ 1:6G1c , a crack arrest occurs on the material interface. For a slightly smaller energy surface (G2c ¼ 1:59G1c ), there is no arrest and the crack accelerates. For an elasto-plastic behavior of each section of the base metal (except the inclusion which remains elastic), the same parametric analysis is performed. For a tough third section ðG2c P G1c Þ, an arrest always occurs at the joined section between the inclusion and the tough metal. To enhance instability, a new parametric analysis is performed with G2c < G1c . The crack runs beyond the inclusion (without arrest), only under the condition that G2c 6 0:86G1c (Fig. 31), which confirms the strong stabilizing part of plasticity, and the importance to take into account that feature for crack arrest investigation. An important remark is that there is whether an arrest on the material interface or no arrest at all (no arrest has been observed beyond the material interface inside the third metal section, for both behaviors).

27

G. Debruyne et al. / Engineering Fracture Mechanics 90 (2012) 1–29

2500

Crack velocity (m/s)

2000

1500

1000

Gc2/Gc1=0.86 500

Gc2/Gc1=0.87

0 0

0,002

0,004

0,006

0,008

0,01

0,012

-500

Crack growth (m) Fig. 31. Crack tip velocity for an elasto-plastic media with multi fracture properties.

6. Conclusions and future outlooks In this paper, some dynamic phenomena due to material interfaces separating zones with different fracture toughness, are studied, with the help of analytical or numerical tools. One of the main results concerns the crack arrest. Dynamic and static solutions concerning this aspect are very different from each other (the predicted crack jump in the beam specimen or the predicted debonding length in the peel test are much larger in the dynamic analysis than in the quasi-static one). Numerical results, derived from the Cohesive Model developed in the second section, are used for the cantilever beam and the pressure vessel investigations. The main results connected to each section may be summarized as follows: – Regarding to the peel test of a thin film, debonded from a rigid substrate with two zones of different interfacial properties, the main conclusion is that the toughness ratio is the guide line for the kinematics, regardless of the absolute value of each zone’s toughness. The peel arrest point is underestimated by a conventional quasi-static method (that is equilibrium associated to a Griffith criterion). On the other hand, from a simple energy balance, a static prediction is achieved, which gives the same result as the exact dynamic solution (obtained with the characteristics method) for the arrest length (and only for that point). – When an inclusion, with weak interfacial properties, is inserted into the film, the peel kinematics depends on the inclusion length. For a sufficiently long inclusion, it amounts to the previous problem. For a short inclusion, the peel arrest occurs at the ending inclusion interface (with a transient stage converging to the quasi-static solution if the loading is not frozen). In this case, the static prediction relying on energy balance is inappropriate for predicting the arrest. – Concerning the fast crack growth in a cantilever beam with a material interface, finite elements computations associated to a cohesive zone model have been performed, considering a large number of toughness ratios. To assess the reliability of the method, some comparisons have been done with a node release technique associated to the energy release rate concept. The results are in good accordance (which is natural because the cohesive zone is small compared to the crack length, so that the cohesive solution is close to a pure Griffith solution). Besides, no significant differences have to be noticed between the different cohesive laws (linear or exponential) during the propagation stage. Once again, the toughness ratio is a key point for the crack arrest length (other features as the initial toughness in the first material and the initial crack velocity have much less influence on the crack kinematics). Compared to the peel test, for a given ratio, the arrest length is shorter than the debonding length in the peel test). An important result is that the static prediction of the crack arrest, relying energy balance, still works for this problem (but the prediction is not as accurate as in the peel test, it’s only a conservative estimate). – In the case where a small elastic inclusion, with a weak toughness, is inserted into the beam, the cohesive model predicts an arrest precisely at the ending interface (for a moderate ratio a = 3), at that whatever the elastic or plastic behavior of the host material. In the case of an elastic host material, the analytical static arrest prediction, mentioned above, does not hold because much kinetic energy remains at the arrest time. In the case of a plastic material, the kinetic energy is much smaller, but this method could not be tested, due to the lack of plastic analytical solution. – Concerning the pressure vessel affected by a subclad crack, and submitted to an inner pressure loading, the main feature investigated is the role of plasticity on the crack kinematics. The study shows two situations.

28

G. Debruyne et al. / Engineering Fracture Mechanics 90 (2012) 1–29

The first deals with crack growth inside a vessel wall with an uniform toughness, and in this case two extreme situations occur according to the material behavior: (1) If the material is elastic, the kinematics is fully unstable, without crack arrest. (2) If the material is elasto-plastic, the kinematics is stable. Even if the plastic behavior is disturbed by putting a small elastic part ahead the initial crack tip, which keeps the same fracture properties, only a short unstable growth is observed and the crack stops inside the elastic part. The second situation deals with a small brittle elastic defect embedded in a tougher base metal, and located on the crack path. Once again, two situations may occur according to the material behavior: (1) If the base material is elastic, the crack is pushed forward throughout the defect, and the arrest takes place at the very end of the brittle defect, and that only if the material beyond this material interface is sufficiently tough with respect to the initial base metal. Otherwise, no arrest is detected. (2) If the base metal is elasto-plastic (the defect keeping its elastic properties), the arrest occurs at the same point as in the previous case, even if the metal beyond the defect is slightly softer than the initial metal before the inclusion. Naturally, if this material is brittle enough, no arrest occurs. At any case, no arrest has been observed within the brittle defect or beyond the material interface, but always on the brittle inclusion interface, which confirms the conclusions drawn from the academic peeling test. A next step in further investigations should be the analysis of 3D and viscosity effects. An important point would be the search for a simple method, predicting the crack arrest, avoiding time-consuming dynamic computations. A first step has been done here in the right direction, using an energy balance between the crack onset and the crack arrest time. Acknowledgments The authors acknowledge Prof. J.J. Marigo from Paris VI university for advices about fundamental dynamic features, and appreciatively thank Prof. S. Andrieux and Dr. P. Massin from LaMSID-EDF for helpful discussions. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13]

[14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28]

Freund LB. Dynamic fracture mechanics. Cambridge monographs on mechanics and applied mathematics; 1990. Griffith AA. The phenomena of rupture and flow in solids. Philos Trans Roy Soc Lond A 1921;221:163–98. Freund LB. A simple model of the double cantilever beam crack propagation specimen. J Mech Phys Solids 1977;25:69–79. Dumouchel P-E, Marigo J-J, Charlotte M. Dynamic fracture: an example of convergence towards a discontinuous quasistatic solution. Continuum Mech Thermodyn 2008;20:1–19. Barenblatt GI. The formation of equilibrium cracks during brittle fracture: general ideas and hypothesis, axially symmetric cracks. Appl Math Mech 1959;23:622–36. Dugdale DS. Yielding of steel sheets containing slits. J Mech Phys Solids 1960;8:100–4. Hillerborg A, Modeer M, Petersson PE. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cem Concr Res 1976;6(6):163–8. Needleman A. A continuum model for void nucleation by inclusion debonding. J Appl Mech 1987;54:525–31. Tvergaard V. Effect of fibre debounding in a whisker-reinforced metal. Mater Sci Engng 1990;125:203–13. Costanzo F, Walton JR. A study in dynamic crack growth using a cohesive zone model. Int J Engng Sci 1997;35(12–17):1085–114. Laverne J. Formulation énergétique de la rupture par des modèles de forces cohésives. PhD thesis. Paris 13 University; 2004. Francfort GA, Marigo J-J. Revisiting brittle fracture as an energy minimization problem. J Mech Phys Solids 1998;46(8):1319–42. Charlotte M, Francfort GA, Marigo J-J, Truskinovsky L. Revisiting brittle fracture as an energy minimization problem: comparison of Griffith and Barenblatt surface energy models. In: Benallal A, editor. Proceedings of the symposium on ‘‘Continuous Damage and Fracture’’. Paris: The Data Science Library, Elsevier; 2000. p. 7–18. Kanninen MF, Popelar C, Gehlen PC. Dynamic analysis of crack propagation and crack arrest in the double-cantilever-beam specimen. In: Hahn GT, Kanninen MF, editors. Fast fracture and crack arrest. ASTM STP 627. Philadelphia: American Society for Testing and Materials; 1977. p. 19–38. Hellan K. An alternative one-dimensional study of dynamic crack growth in DCB specimens. Int J Fracture 1981;17. Lo CY, Nakamura T, Kushner A. Computational analysis of dynamic crack propagation along bimaterial interface. Int J Solids Struct 1994;31:145–68. Nakamura T, Kushner A, Lo CY. Interlaminar dynamic crack propagation. Int J Solids Struct 1995;32:2657–75. Masson R, Nicolas L, Moinereau D. RPV structural integrity assessment during a PTS event: application of an extended Beremin model consistent with WPS test results. In: ASME PVP conference, Vancouver, August 2002. Cheverton RD, Ball DG. The role of crack arrest in the evaluation of PWR pressure vessel integrity during PTS transients. Engng Fract Mech 1986;23(1):71–80. Jung J, Kanninen MF. An analysis of dynamic crack propagation and arrest in a nuclear pressure vessel under thermal shock conditions. J Press Vessel Technol 1983;105(2):111–6. Burdekin FM, Knott JF, Sumpter JDG, Sherry AH. TAGSI views on aspects of crack arrest philosophies for pressure vessels with thicknesses up to 100 mm. Press Vessels Pip 1999;76:879–83. Charlotte M, Debruyne G, Marigo J-J. Static and dynamic peeling of a thin film and crack arrest. EDF R&D internal report H-T64-2007-63259-EN; 2008. Dumouchel P-E. Propagation brutale de fissures et effects dynamiques: application industrielle à la longueur d’arrêt. PhD thesis. Paris 6 University; 2008. Freund LB. A simple model of the double cantilever beam propagation specimen. J Mech Phys Solids 1977;25:69–79. Arata JJM, Needleman A. The effect of plasticity on dynamic crack growth across an interface. Int J Fracture 1998;94:383–99. Rice JR. The mechanics of earthquake rupture. In: Proceeding of International School of Physics Enrico Fermi, North Holland; 1980. p. 555–649. Falk ML, Needleman A, Rice JM. A critical evaluation of cohesive zone models of dynamic fracture. J Phys IV France 2001;11:43–50. Andersson H. A finite element representation of stable crack-growth. J Mech Phys Solids 1973;21(5):337–56.

G. Debruyne et al. / Engineering Fracture Mechanics 90 (2012) 1–29

29

[29] Brickstad B, Nilsson F. Dynamic analysis of crack growth and arrest in a pressure vessel subjected to thermal and pressure loading. Engng Fract Mech 1986;23(1):61–70. [30] Batra RC, Love BM. Crack propagation due to brittle and ductile failures in microporous thermoelastoviscoplastic functionally graded materials. Engng Fract Mech 2005;72:1954–79. [31] Engquist B, Majda A. Absorbing boundary conditions for the numerical simulation of waves. Math Comput 1977. [32] Kanninen MF. An analysis of dynamic crack propagation and arrest for a material having a crack speed dependent fracture toughness. In: Sih GC et al., editors. Prospects of fracture mechanics. Leyden: Noordhoff; 1974. p. 251–66. [33] Kanninen MF, Popelar CH. advanced fracture mechanics. Oxford Engng Sci Ser 1985;15:216–7. [34] Jung J, Kanninen MF. An analysis of dynamic crack propagation and arrest in a nuclear pressure vessel under thermal shock conditions. J Press Vessel Technol 1983;105:111–6. [35] Priest AH. An energy balance in crack propagation and arrest. Engng Fract Mech 1998;61:231–51. [36] Ahmad J, Jung J, Barnes CR, Kanninen MF. Elastic–plastic finite element analysis of dynamic fracture. Engng Fract Mech 1983;17(3):235–46. [37] Code_Aster website. . [38] Destuynder Ph, Djaoua M. Sur une interpretation Mathématique de l’Intégrale de Rice en Théorie de la Rupture Fragile. Math Meth Appl Sci 1981;3:70–87 [AMS subject classification: 73 M05]. [39] Lorentz E, Badel P. A new path-following constraint for strain-softening finite element simulations. Int J Numer Meth Engng 2004;60(2):499–526.