# Stable crack propagation in steel at 1173K: Experimental investigation and simulation using 3D cohesive elements in large-displacements

Nicolas Tardif, Alain Combescure, Michel Coret, Philippe Matheron

## ► To cite this version:

Nicolas Tardif, Alain Combescure, Michel Coret, Philippe Matheron. Stable crack propagation in steel at 1173K: Experimental investigation and simulation using 3D cohesive elements in large-displacements. Engineering Fracture Mechanics, 2010, 77 (5), pp.776-792. 10.1016/j.engfracmech.2009.12.015. hal-01007318

HAL Id: hal-01007318

https://hal.science/hal-01007318

Submitted on 3 Jun 2017

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Stable crack propagation in steel at 1173 K: Experimental investigation and simulation using 3D cohesive elements in large-displacements

N. Tardif<sup>a,b</sup>, A. Combescure<sup>b</sup>, M. Coret<sup>b</sup>, P. Matheron<sup>c</sup>

<sup>a</sup> ISRN, BP17 92262 Fontenay-aux-Roses Cedex, France

<sup>b</sup> Université de Lyon, CNRS INSA-Lyon, LaMCoS, UMR5259, 20 Avenue Albert Einstein, F69621 Villeurbanne Cedex, France

<sup>c</sup> CEA/DEN/DM2S, Service d'Etudes Mécaniques et Thermiques, 91191 Gif-sur-Yvette Cedex, France

This paper deals with the numerical simulation of a stable crack propagation experiment at 1173 K in a 16MND5 steel. At this temperature, the material is viscoplastic. A cohesive zone model is formulated in order to simulate the rupture of a CT specimen. A large displacement 3D cohesive element with eight nodes is implemented in the finite element code ABAQUS. The associated traction-separation law is of Tvergaard and Hutchinson type, in which an hardening slope has been added. This hardening simulates the material strengthening associated to the increasing strain rate in front of the crack tip when crack tip starts to propagate.

We show that in this case the form of the cohesive law has great impact on the simulated propagation velocity.

Keywords:
Quasi-static crack growth
High temperature
Ductile fracture
3D cohesive zone model

## 1. Introduction

In the event of a serious accident on a pressurized water reactor (PWR) involving fusion of the fuel core and internal structural elements, the bottom head is subjected to significant thermal and mechanical solicitations. The TMI-II accident of March 28th, 1979 at Three Mile Island, USA is a reference in terms of serious accident scenarios due to the data collected during the event. During this accident, 20 tons of highly calorific corium, combining molten metal and oxide, were relocated in the bottom head. This corium bath led to a 1300 K temperature increase of the bottom head in 30 min, with a 1 m-diameter hot spot. During the accident, the pressure on the reactor vessel was 10 MPa. However, in spite of numerical predictions to the contrary, no cracking occurred in the vessel.

In the hypothetical case where the reactor vessel would suffer cracking, it is necessary to predict the time of occurrence of a crack, its size and its location in order to determine an out-of-vessel accident management strategy within the confinement system. Experimental programs [1-3] have been carried out in order to characterize these parameters.

In this paper, our objective is to study more particularly the initiation and quasi-static (or stable) Mode I crack propagation in pressure vessel steel (16MND5) at 1173 K [4]. In the course of this work, cracking tests were carried out on side-grooved CT specimens. We propose a numerical simulation of these tests using a 3D finite element model and the ABAQUS 6.6-1 analysis code [5]. The crack's initiation and propagation are modeled using cohesive surface elements under large-displacement and large-strain assumptions.

The cohesive zone concept was proposed independently by Barenblatt [6] and Dugdale [7] in order to regularize the crack-tip singularity predicted by linear elastic fracture mechanics [8]. Barenblatt reckoned that near the crack's tip the
Nomenclature

t time

T Cauchy stress tensor

D strain rate tensor

u displacement vector

Δu displacement vector increment

Δu' virtual displacement vector increment

F force density on the bulk

t traction in the cohesive zone

δ(δn, δs, δt) separation vector in the cohesive zone, with its normal, shear and tangential components

R global to local coordinate system matrix rotation

U(U, V, W) nodal displacement vector of components U, V and W (global coordinate system)

U(ξ, η, ζ) nodal displacement vector of components ξ, η and ζ (local coordinate system)

N<sub>i</sub> finite element shape function expressed in the local coordinate system of the element

K<sub>f</sub> tangent stiffness matrix

J Jacobian

w<sub>g</sub> integration weight of Gauss point g

S element area

R toughness of a ductile material

B thickness of the specimen

U internal energy

E<sub>e</sub> elastic energy

E<sub>pl</sub> inelastic energy

a crack progression

Γ<sub>0</sub> cohesive surface energy

C traction-separation law

E Young modulus

U<sub>UN</sub> inelastic energy dissipated in volume

U<sub>CZM</sub> energy dissipated in the cohesive zone

σ<sub>h</sub> average stress

σ<sub>eq</sub> Von Mises stress

ε<sub>eq</sub> accumulated plastic strain

interatomic forces which tend to bring the two surfaces being pulled apart back together cannot be ignored in the analysis. Dugdale introduced an expanded crack-tip zone in order to limit the development of stresses to a few times the yield stress.

Thus, the role of cohesive zone models is to reproduce the mechanical effects of the fracture process zone at the crack's tip. The zone is modeled using two surfaces which, initially, coincide. During the loading, the two surfaces come apart until complete decohesion. Then, they simulate the lips of the crack. A traction-separation law governs the traction between the two surfaces as a function of their separation and, thus, models the fracture phenomena associated with the material being studied.

Cohesive models have been used in the fracture simulation of numerous materials, including concrete [9–11], polymers [12–14], ceramics [15–17], and the interfaces between different materials or material regions [18–20].

The application of the cohesive model (CM) to the ductile fracture of metals was initiated primarily by Tvergaard and Hutchinson [21], Siegmund and Brocks [22] and Tvergaard [23]. The traction-separation laws they proposed are appropriate to the phenomena of nucleation, growth and coalescence of voids [24] which are encountered in ductile fracture.

Only few studies have been carried out using cohesive elements in a 3D formulation [25,26] and under large-strain assumptions [27–29]. Chen et al. [27], however, showed in an analysis of the tearing of CT specimens that a plane strain simulation is unable to reproduce the triaxiality state at the crack's tip which can be obtained near the middle section of the 3D model using an identical traction-separation law.

Consequently, the two formulations lead to different crack propagation predictions. Therefore, a 3D analysis is necessary. Besides, the high ductility of our material at high temperatures requires the use of a large-displacement, large-strain formulation.

It is commonly admitted that the energy and the cohesive stress are the only relevant parameters which characterize crack propagation in a ductile material. A few authors have cast doubt on this assumption [30,31]. We will show that in our case of a viscoplastic, highly ductile material at 1173 K the form of the traction-separation law affects the propagation significantly.

Section 2 presents the tests on CT specimens and describes the material, the geometry of the specimens, the testing apparatus and the results.

Section 3 introduces the cohesive zone model. The large-displacement 3D formulation of the cohesive elements is presented in detail. We also discuss the modification of the form of Tvergaard and Hutchinson's traction-separation law with respect to the viscous behavior of the material.
<table><caption>**Table 1**
Chemical analysis of "Krakatoa"-grade 16MND5 pressure vessel steel in thousandths of % mass.</caption><thead><tr><th>C</th><th>Si</th><th>Mn</th><th>P</th><th>S</th><th>Ni</th><th>Cr</th><th>Mo</th><th>V</th><th>Cu</th><th>Co</th><th>Fe</th></tr></thead><tbody><tr><td>140</td><td>160</td><td>1310</td><td>10</td><td>7</td><td>790</td><td>190</td><td>490</td><td>7</td><td>50</td><td>10</td><td>Balance</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Geometry of the specimen.</figcaption></figure>

<table><caption>**Table 2**
Precrack profile in relation to the initial groove.</caption><tr><td>Thickness (mm)</td><td>0.0</td><td>2.5</td><td>5.0</td><td>7.5</td><td>10.0</td><td>12.5</td><td>15</td><td>17.5</td><td>20</td><td>22.5</td><td>25</td></tr><tr><td>Precrack (mm)</td><td>3.39</td><td>2.45</td><td>2.23</td><td>2.16</td><td>2.12</td><td>2.15</td><td>2.11</td><td>2.13</td><td>2.20</td><td>2.47</td><td>3.57</td></tr></table>

The modeling of the tests is described in Section 4.

Section 5 presents the optimized fitting results of the traction-separation law in order to reproduce the tests. Based on these results, different assumptions concerning the simulation options and the traction-separation law (the need to introduce some viscosity, the dependence on the degree of triaxiality) are validated.

Section 6 proves that in the case of tearing tests of a CT20 specimen at high temperature (1173 K) the form of the traction-separation law has significant impact, not only locally in the fracture process zone, but also on the global response of the tearing tests.

## 2. The material and test setup

The material being studied is a "Krakatoa"-grade 16MND5 steel whose chemical analysis is presented in Table 1.

The specimen had an unconventional geometry based on a CT20 (see Fig. 1): some dimensions were increased in order to reduce necking at high temperature. Lateral grooves were made in order to approximate a plane strain state. A precrack whose profile is described in Table 2 was obtained by fatigue.

The tests were carried out under primary vacuum in a furnace consisting of a water-cooled dual-wall chamber. The temperature was controlled through a B-type thermocouple and its homogeneity measured with three K-type thermocouples. The temperature loading was first applied at a rate of $1200 \text{ K h}^{-1}$ up to 1173 K. Throughout that phase, a constant force was maintained to compensate for depressurization.

The actual tearing test was performed under prescribed displacement at a rate of $5 \text{ mm min}^{-1}$ in order to minimize load relaxation. The displacement set point was ended at 15 mm of the actuator's range. The load and displacement were measured using the sensors of the testing device.

The crack's propagation was determined using an ANS crack monitoring device. This system works by injecting a current into the specimen through welded wires in order to measure the potential difference between the two sides of the defect. A crack's propagation leads to a decrease in section characterized by an increased resistance and, therefore, an increased potential difference [32]. This method was calibrated during preliminary tests performed at different final extensions of the actuator. Besides, the signal thus obtained is noisy and must be post-processed using linear smoothing. In our case, this smoothing was achieved over a 2 mm crack increment. Consequently, the results at both the initiation and the beginning of the crack's path were extrapolated.

The tests were performed twice. The two results were similar in terms of load-displacement and crack's progression curves. Fig. 2 illustrates the results of one of the tests. Fig. 3 shows the deformed shape and the crack's profile at the end of the test.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Load and crack's progression curves.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. The deformed specimen at the end of the test.</figcaption></figure>

### 3. The cohesive zone model

The cohesive element and its traction-separation law were implemented as a User Element Fortran subroutine into the ABAQUS/Standard analysis code.

#### 3.1. Principle of virtual work under large-displacements

A solid subjected at time t to surface loads F over $\partial B_{text}^1$ and prescribed displacements $\mathbf{u}_{imp}$ over $\partial B_{text}^2$ is traversed by a cohesive surface $S_{tcoh}$ (Fig. 4). Neglecting the terms relative to volume forces and inertia, and using the configuration of the solid at time t as the reference, one can express the principle of virtual work as:

$$ \int_{B_t^+ \cup B_t^-} \mathbf{T} : \mathbf{D}(\Delta \mathbf{u}^*) dv + \int_{S_{tcoh}} \mathbf{t} \cdot \delta(\Delta \mathbf{u}^*) ds = \int_{\partial B_{text}^1} \mathbf{F} \cdot \Delta \mathbf{u}^* ds \quad \text{with} \quad \Delta \mathbf{u}^* \in \mathcal{U}_{ad} = \{u \in \mathcal{U} : u|_{\partial B_{text}^2} = \mathbf{u}_{imp}\}, \quad (1) $$

where $\delta$ denotes the separation in the cohesive zone, D the strain rate tensor and T the Cauchy stress tensor.
<figure><img src="image_5.png" /><figcaption>Fig. 4. The reference problem.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. An interface element in its initial configuration (a) and its deformed configuration (b), and characterization of the separation of two nodes in the deformed configuration (c).</figcaption></figure>

## 3.2. Formulation of the surface cohesive elements under large-displacements and large strains

The isoparametric bilinear interface element, initially with zero thickness, has eight nodes and four integration points. The nodes are divided into two 4-node surfaces representing a portion of the lower and upper lips of the crack undergoing decohesion.

The movement of the upper and lower parts of the bulk induces a displacement of the nodes of the interface element. Only the relative displacements which represent the three fracture modes are taken into account in the formulation. Here, we are focusing on Mode I alone (opening mode), i.e. the mode which is solicited in the tests. These displacements generate a traction force which prevents their development. This traction is governed by the traction-separation law described in the following section.

In a large-displacement, large-strain formulation, the choice of the reference surface is paramount in characterizing the separation mode of the element [33]. In order to define the local coordinate system of the separation, Ortiz and Pandolfi [29] proposed to use the midsurface $S_t$ calculated from the nodal positions of the upper and lower surfaces. Besides, the large-strain formulation enables one to take into account the area reduction due to necking of the bulk during the crack's progression.

Fig. 5 shows the geometry of the element. Nodes 1–4 belong to the lower side of the element while nodes 5–8 belong to the upper side. Initially, the two surfaces coincide. Each node has three translational degrees of freedom.
The displacements of the element's nodes in the global coordinate system (U, V, W) can be described by a vector U of dimension 24, in which the index represents the node being considered:

$$ \mathbf{U}_{24 \times 1} = (U_1, V_1, W_1, \dots, U_8, V_8, W_8)^T. \qquad (2) $$

The calculation of the separation at any point of the element is associated with a local coordinate system at the center of the element's midplane in its deformed configuration (Fig. 5). This reference coordinate system $(\xi, \eta, \zeta)$ is deduced from the global system by a rotation $\mathbf{R}_{3 \times 3}^T$. $\zeta$ is the normal to this midplane, and $\xi$ and $\eta$ are the two tangent vectors:

$$ (\xi, \eta, \zeta)^T = \mathbf{R}_{3 \times 3}^T (U, V, W)^T. \qquad (3) $$

Now we can express the value of the relative separation $\delta_{12 \times 1}$ of the two surfaces at each of the four Gauss points G, of coordinates $(\xi_g, \eta_g)$ in the element's local coordinate system:

$$ \delta_{12 \times 1}(\xi_g, \eta_g) = \underbrace{\begin{array}{c} -\mathcal{N}_1(\xi_g, \eta_g)\mathbf{R}_{3 \times 3}^T, \quad -\mathcal{N}_2(\xi_g, \eta_g)\mathbf{R}_{3 \times 3}^T, \quad -\mathcal{N}_3(\xi_g, \eta_g)\mathbf{R}_{3 \times 3}^T, \quad -\mathcal{N}_4(\xi_g, \eta_g)\mathbf{R}_{3 \times 3}^T \\ \mathcal{N}_1(\xi_g, \eta_g)\mathbf{R}_{3 \times 3}^T, \quad \mathcal{N}_2(\xi_g, \eta_g)\mathbf{R}_{3 \times 3}^T, \quad \mathcal{N}_3(\xi_g, \eta_g)\mathbf{R}_{3 \times 3}^T, \quad \mathcal{N}_4(\xi_g, \eta_g)\mathbf{R}_{3 \times 3}^T \end{array}}_{\mathbf{B}_{12 \times 24}(\xi_g, \eta_g)} \mathbf{U}_{24 \times 1}. \qquad (4) $$

Functions $\mathcal{N}_i$ correspond to the four usual bilinear polynomial interpolation functions expressed in the local coordinate system of the element.

Then, the traction vector $\mathbf{t}_{12 \times 1}$ at the Gauss point can be calculated using the traction-separation law C presented in the next section:

$$ \mathbf{t}_{12 \times 1} = \mathbf{C} \cdot \delta_{12 \times 1}. \qquad (5) $$

After linearization of the principle of virtual work (Eq. (1)) and finite element discretization, the calculation of the internal forces is expressed on the $n+1$ configuration at the end of the step. Let $J$ denote the Jacobian which transforms the element's local coordinates (Domain $\square$) into global Cartesian coordinates.

$$ \mathbf{F} = \int_{S_{\text{coh}}} \mathbf{B}_{n+1}^T \mathbf{t}_{n+1} ds = \int_{\square} \mathbf{B}_{n+1}^T \mathbf{t}_{n+1} J_{n+1} d\square. \qquad (6) $$

Using a Gauss method for numerical integration, one gets the nodal forces:

$$ \mathbf{F}_{24 \times 1} = S_{n+1} \sum_g \omega_g \mathbf{B}_{n+1}^T(\xi_g, \eta_g) \mathbf{t}_{n+1}(\xi_g, \eta_g), \qquad (7) $$

where $S$ represents the element's area, and $\omega_g$ the integration weight of Gauss point $g$.

The tangent stiffness matrix $\mathbf{K}_T$ is:

$$ \mathbf{K}_T = \int_{S_{\text{coh}}} \mathbf{B}_{n+1}^T \mathbf{D}\mathbf{B}_{n+1} ds = \int_{\square} \mathbf{B}_{n+1}^T \mathbf{D}\mathbf{B}_{n+1} J_{n+1} d\square, \qquad (8) $$

$$ \mathbf{K}_{24 \times 24} = S_{n+1} \sum_g \omega_g \mathbf{B}_{n+1}^T(\xi_g, \eta_g) \frac{d\mathbf{t}}{d\delta_{n+1}}(\xi_g, \eta_g) \mathbf{B}_{n+1}(\xi_g, \eta_g). \qquad (9) $$

In this expression, $\mathbf{D}$ is the tangent constitutive relation which links the stress to the strain. The expression of $\mathbf{K}_T$ given here is a simplified version valid only for small strains and small rotations.

This definition is sufficient for the case being considered in this work, where the cohesive zone remains in the plane of symmetry and does not rotate.

### 3.3. The traction-separation law

The traction-separation law represents, on the macroscale, the dissipative phenomena leading to the complete decohesion of the crack's lips in the fracture process zone at the crack's tip.

In the case of a metal at low temperature, the traction-separation law represents the local processes of void nucleation, growth and coalescence. The form of the traction-separation law is unknown a priori.

Tvergaard and Hutchinson [20,21] proposed the law described in Fig. 6. This law is relatively simple, with four parameters, and consistent with an irreversible deformation of the process zone and perfect plasticity of the surrounding bulk.

The same authors showed that in the case of a nonviscous elastic-plastic material the global form of their law has little influence on the crack's propagation. Only two parameters (the maximum traction $t_{max}$ and the rate of energy dissipated during the decohesion of the two surfaces $\Gamma_0$) have an influence.

A good candidate for studying this influence, called the material's ductile toughness ($R$), was proposed by Turner and Kolednik [34]. By extension of the concept of energy recovery rate [35], they proposed the following expression for the toughness of a ductile material:

$$ R = \frac{1}{B} \frac{d(U - E_e)}{da} = \frac{1}{B} \frac{dE_{pl}}{da} + \Gamma_0, \qquad (10) $$
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Traction-separation laws.</figcaption></figure>

where *U* denotes the internal energy, *E*<sub>e</sub> the elastic energy, *B* the thickness of the specimen, *E*<sub>pl</sub> the inelastic energy and *a* the crack's progression.

*R* represents the dissipated energy required for the extension of the cracked surface by an increment *Bda*. This energy is divided into two contributions: an energy $\frac{dE_{pl}}{da}$ dissipated inelastically in the bulk and a surface energy $\Gamma_0$ dissipated during the creation of the two lips of the crack. The latter contribution corresponds to the parameter $\Gamma_0$ of the traction-separation law. While these two contributions can be easily separated by modeling a crack with cohesive elements, their experimental characterization has scarcely been investigated [36,37].

Let us consider a constitutive relation of the elastic-plastic bulk with any type of strain hardening. For a constant $\Gamma_0$, the larger the maximum traction $t_{max}$ relative to the yield strength of the bulk, the larger the parameter *R* because the rate of energy dissipated in the bulk through plasticity increases [30]. For a constant $t_{max}$, the larger the surface energy $\Gamma_0$, the larger *R*, both intrinsically and because the size of the fracture process zone increases and causes an increase in the plastic energy dissipated.

In our investigation, the temperature is very high (1173 K). Viscous effects in the bulk material become important. Hence, the localization process leading to fracture is time-dependant.

If the loading rate is slow, the damage process leading to fracture includes diffusive mechanisms (grain boundary cavi-
tation and grain boundary diffusion). The traction-separation law is used in this paper to simplify the localization process
which leads to fracture. Consequently, the traction-separation law should be described by a viscous model.

If it is fast, as in this experiment (the experimental loading pin velocity is set constant at 5 mm min⁻¹), the fracture is of
plastic ductile type depending on the strain rate effect [38]. Hence, the traction-separation law should be simply rate depen-
dant elastoplastic curves.

However, the only available tests are constant strain rate (loading pin velocity is constant). Therefore, the cohesive zone model will be elastic-plastic and independent of the global opening rate, which is supposed to be constant during the whole test.

Near the crack tip, in the first stage of the propagation, the strain rates increase and results in the strengthening of the
properties of the bulk due to the singularity. However, this strengthening is attenuated by the damage in the zone surround-
ing the crack's tip. In order to account for this phenomenon simply, Tvergaard and Hutchinson's constitutive law was mod-
ified by adding a hardening slope, as shown in Fig. 6b.

Besides, the cracking mode of the CT specimen is Mode I. Therefore, mixed modes were not implemented. Spurious dis-
placements in mixed mode are controlled through a simple elastic stiffness and do not affect the response in Mode I.

The interface is considered to be ruptured once all cohesion has been lost. Then, the separation is greater than $\delta_c$.

**4. The finite element model**

4.1. *Mesh and boundary conditions*

The mesh of the CT specimen, which consists of about 46,000 elements, is shown in Fig. 7.

The elements which constitute the bulk are brick elements with reduced integration (eight nodes and one integration point). The mesh size transition near the specimen's fillets is provided through wedge elements (six nodes and two integration points).

The elements of the cohesive zone are laid out appropriately in order to follow the profile of the precrack. The surface is meshed regularly with 43 elements in the direction of the propagation by 20 elements in the direction of the front (860 elements). Two opposite nodes near the central position of the precrack are fixed in the direction of the propagation (U₁).

The pins are modeled with rigid shell elements with their reference points located along the axis of the pin. Rotations of the pins other than about their axis (U₂) and translation about this axis are not allowed. The displacement *u*(*t*) of the axis of each pin in the traction direction (U₃) is set at 2.5 mm min⁻¹. Thus, the midplane of the cohesive zone remains orthogonal to (U₃) and undergoes no rotation (see Section 3.2).
<figure><img src="image_8.png" /><figcaption>Fig. 7. Mesh and boundary conditions.</figcaption></figure>

<table><caption>**Table 3**

Behavior of the bulk: plastic part.</caption><thead><tr><th colspan="2">1E-5</th><th colspan="2">5E-4</th><th colspan="2">1E-3</th><th colspan="2">5E-3</th></tr><tr><th>ε<sub>eq</sub><sup>b</sup></th><th>σ<sub>eq</sub><sup>a</sup></th><th>ε<sub>eq</sub></th><th>σ<sub>eq</sub></th><th>ε<sub>eq</sub></th><th>σ<sub>eq</sub></th><th>ε<sub>eq</sub></th><th>σ<sub>eq</sub></th></tr></thead><tbody><tr><td colspan="8">Logarithmic strain rates (s<sup>-1</sup>)</td></tr><tr><td>0.000</td><td>19.93</td><td>0.000</td><td>35.33</td><td>0.000</td><td>41.01</td><td>0.000</td><td>46.81</td></tr><tr><td>1.422E-3</td><td>22.73</td><td>2.050E-3</td><td>43.17</td><td>6.778E-3</td><td>50.63</td><td>9.320E-3</td><td>63.86</td></tr><tr><td>3.147E-3</td><td>26.13</td><td>2.232E-2</td><td>52.58</td><td>2.522E-2</td><td>59.64</td><td>1.647E-2</td><td>71.08</td></tr><tr><td>5.699E-2</td><td>29.21</td><td>4.655E-2</td><td>58.52</td><td>5.294E-2</td><td>66.12</td><td>3.218E-2</td><td>79.17</td></tr><tr><td>1.001E-1</td><td>29.93</td><td>7.119E-2</td><td>62.29</td><td>7.998E-2</td><td>71.08</td><td>4.686E-2</td><td>84.37</td></tr><tr><td>3.000E-1</td><td>29.93</td><td>9.334E-2</td><td>64.59</td><td>1.096E-1</td><td>75.6</td><td>7.521E-2</td><td>91.88</td></tr><tr><td></td><td></td><td>1.067E-1</td><td>66.32</td><td>1.438E-1</td><td>79.55</td><td>1.139E-1</td><td>98.53</td></tr><tr><td></td><td></td><td>3.000E-1</td><td>66.32</td><td>1.733E-1</td><td>82.21</td><td>1.738E-1</td><td>105.18</td></tr><tr><td></td><td></td><td></td><td></td><td>3.000E-1</td><td>85.17</td><td>1.904E-1</td><td>106.33</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td>2.500E-1</td><td>107.49</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td>3.000E-1</td><td>107.49</td></tr></tbody></table>

<footnote>a Von Mises' (Cauchy) stress in MPa.</footnote>

<footnote>b Accumulated (logarithmic) plastic strain.</footnote>
<figure><img src="image_9.png" /><figcaption>Fig. 8. Set of curves used to define the behavior of the bulk.</figcaption></figure>

## 4.2. The calculation procedure

A quasi-static calculation was performed, but the loading rate was respected in order to account properly for the viscous phenomena associated with the material.

The subintegration associated with the elements of the bulk resulted in hourglass displacement modes [39] around the elements adjacent to the cohesive zone. These displacement modes can perturb the propagation of the crack because of the formation of grooves whose amplitude increases with the loading. We chose to stabilize these modes using additional stiffnesses.

## 4.3. The constitutive relation of the bulk

Because of the temperature (1173 K), the material's behavior is highly dependent on the strain rate. This dependence was taken into account by introducing, at different logarithmic strain rates, a scattered pattern of elastic-plastic curves with strain hardening. The code used linear interpolation to define the behavior at intermediate strain rates.

The Young's modulus E of the material is independent of the strain rate and equal to 23,466 MPa.

The elastic-plastic characteristics used in the analysis are given in Table 3 and represented graphically in Fig. 8.

## 4.4. Constitutive law of the interface

For the purpose of the simulations, certain parameters characterizing the form of the traction-separation law, given in Fig. 6b, were assigned.

Parameter $\delta_1$ was chosen such that:

$$ \frac{t_1}{\delta_1} = \frac{E}{10^{-5}}. \qquad (11) $$

That was because some authors showed the limits of validity of a traction-separation law with zero initial traction [31,40]. The cracking behavior predicted by cohesive zone models with elastic initial response is mesh-dependent and, thus, different from that predicted by models with infinite initial stiffness [31]. We cannot avoid using this type of law because of the solver used by ABAQUS (Newton Raphson). However, we minimized its influence by choosing the highest possible stiffness with respect to the stiffness of the bulk (convergence limit). Besides, forcing the crack to travel through a single cohesive zone reduces the undesirable effects of this type of law.

The value of $\delta_2$ was chosen as proposed by Cornec et al. [25]:

$$ \delta_2 = 0.75\delta_c. \qquad (12) $$

Under these conditions, the value of $\delta_c$ can be deduced from the choice of the three other parameters $t_1$, $t_2$, $\Gamma_0$:

$$ \delta_c = \frac{2\Gamma_0 + \delta_1 t_2}{t_2 + 0.75t_1}. \qquad (13) $$

Now let us study the influence of parameters $t_1$ and $t_2$ on the simulation of the rupture.
<figure><img src="image_10.png" /><figcaption>Fig. 9. The traction-separation law.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Comparison with experimental results.</figcaption></figure>

# 5. Simulation of the tests

This section presents the best simulation of the tests, which corresponds to the following choice of the parameters: $Γ_0 = 47,000$ J/m$^2$, $t_1 = 164.52$ MPa, $t_2 = 197.43$ MPa. The corresponding traction-separation law is shown in Fig. 9.

## 5.1. Comparison with the test results

The global load-displacement and crack propagation-displacement curves are in very good agreement with the experimental results (Fig. 10). The difference between the simulation and the test at the beginning of the crack's propagation is due to the simplification of the experimental results (see Section 2).

## 5.2. Results in terms of energy: validation of the control of the displacement hourglass modes

Fig. 11 shows the internal energy distribution during the propagation. Two scales are used in order to facilitate the interpretation: the lines use the left-hand scale and the dots use the right-hand scale. The main dissipation source is the energy dissipated inelastically in the bulk (88.9% of the internal energy at the end of the calculation), whereas the energy dissipated
<figure><img src="image_12.png" /><figcaption>Fig. 11. Internal energy distribution.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. Energy dissipation rate.</figcaption></figure>

in the creation of the crack's lips in the cohesive zone represents only 9.3% of the internal energy at the end of the calculation. Besides, one can observe that the energy used to control the hourglass modes is negligible compared to the energy dissipated in the cohesive zone (≈0.4% of the internal energy at the end of the calculation), which justifies the use of this mode of control in the case of our simulation.

Fig. 12 shows the distribution of the energy dissipation rate R during the propagation. The numerical calculation was carried out by fitting the curves of the inelastic energy dissipated in volume $U_{UN}$ and the energy dissipated in the cohesive zone $U_{CZM}$ as functions of the average progression of the crack $\bar{a}$. The width of the crack's front used was approximated by the initial thickness of the specimen $B_i$ (25 mm in the useful zone):

$$R \equiv \frac{1}{B_i} \left( \frac{dU_{UN}}{d\bar{a}} + \frac{dU_{CZM}}{d\bar{a}} \right). \qquad (14)$$

One can verify that the rate of energy dissipated in the cohesive zone corresponds to $\Gamma_0$ in the stationary part of the propagation.
<figure><img src="image_14.png" /><figcaption>Fig. 13. Separation velocity in the median zone of the specimen's width.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 14. Maximum triaxiality along the front.</figcaption></figure>

### 5.3. Validation of the constant opening velocity assumption in the simulation of the test

The separation velocity of the crack's lips was calculated as the average over the Gauss points located along a line in the propagation direction (≡$U_1$), near the median position of the specimen's width. Only Gauss points including a separation between $\delta_c/4 < \delta < \delta_c$ were taken into account in the calculation. Fig. 13 shows that the separation velocity remained constant throughout the test ($2 \times 10^{-2} < V < 4 \times 10^{-2}$). Therefore, it was unnecessary to use a viscous traction-separation to model this test.

### 5.4. Dependence of the traction-separation law on the triaxiality rate

The local process of nucleation, growth and coalescence of cavities which is specific to ductile rupture has been shown to be dependent on the stress state [41]. An appropriate way to detect this phenomenon is to measure the triaxiality rate $H$, which can be calculated from the average stress $\sigma_h$ and Von Mises' equivalent stress $\sigma_{eq}$:

$$H = \frac{\sigma_h}{\sigma_{eq}}. \qquad (15)$$

In Siegmund and Brocks' works [22], the traction-separation law of the cohesive zone depends on the triaxiality rate, which is assumed to be equal to that of the adjacent element layer in the bulk. Then, the parameters (surface energy dissipation
<figure><img src="image_16.png" /><figcaption>Fig. 15. Location of the Gauss points at maximum triaxiality.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 16. Study of the influence of the second slope of the traction-separation law.</figcaption></figure>

pation rate and maximum stress) are modified. Based on a 2D plane strain analysis, Siegmund and Brocks showed that in the case of a CT specimen (high-constraint specimen) it is not necessary to take this dependence into account. This result is no longer valid for MT specimens (low-constraint specimens).

However, Chen et al. [27] showed on a 3D analysis of CT10 specimens that the variation of the triaxiality rate along the crack's front throughout the thickness must be taken into account in order for the crack's propagation, and particularly the tunneling effect, to be represented correctly.

We calculated the triaxiality rates for our specimen's geometry (side-grooved CT 25). Fig. 14 shows the maximum triaxiality rate calculated near the Gauss points of the upper layer of bulk elements adjacent to the cohesive zone. Fig. 15 illustrates the position of the corresponding Gauss points and the position of the crack's front in the (U1, U2) plane (Fig. 7). The triaxiality rate is constant over 80% of the width of our specimen, compared to 50% in the case of the CT10 specimen calculated by Chen et al. Therefore, we chose not to make the traction-separation law dependent on the triaxiality rate.

**6. Dependence on the form of the traction-separation law**

As already mentioned in Section 3.3, it is usually agreed that for a ductile material only the surface energy recovery rate $\Gamma_0$ and the maximum stress $t_2$ have an influence on the crack's propagation.
<figure><img src="image_18.png" /><figcaption>Fig. 17. Influence on the toughness ($\Gamma_0$ = 50 kJ m⁻², $t_2$ = 183 MPa).</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 18. Influence on the global load-displacement curve and on the propagation ($\Gamma_0$ = 50 kJ m⁻², $t_2$ = 183 MPa).</figcaption></figure>

<table><caption>**Table 4**
Global variations.</caption><thead><tr><th rowspan="2">* ≡ 100 · (x<sub>1</sub> - x<sub>2</sub>) / x<sub>1</sub></th><th rowspan="2">Case 1</th><th colspan="2">Case 2</th><th colspan="2">Case 3</th></tr><tr><th>Value</th><th>Gain* (%)</th><th>Value</th><th>Gain</th></tr></thead><tbody><tr><td>Propagation velocity (10<sup>-2</sup> mm s<sup>-1</sup>)</td><td>5.3</td><td>4.8</td><td>-9.4</td><td>3.8</td><td>-28.3</td></tr><tr><td>10 &lt; load line displacement &lt; 12 mm</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Maximum load (N)</td><td>5966</td><td>6374</td><td>6.8</td><td>6830</td><td>14.4</td></tr><tr><td>Corresponding displacement (mm)</td><td>3.26</td><td>4.58</td><td>40.5</td><td>7.37</td><td>126.1</td></tr></tbody></table>

However, for the simulation of the tests, it was necessary to add an additional slope to the traction-separation law developed by Tvergaard and Hutchinson. Therefore, we studied the influence of this slope on the propagation results for an identical constitutive relation of the bulk, keeping the values of $\Gamma_0$ and $t_2$ the same. The results are given for 3 values of $t_1$. The corresponding traction-separation laws are shown in Fig. 16.
<figure><img src="image_20.png" /><figcaption>Fig. 19. Distribution of the traction load (about 4 mm propagation, $\Gamma_0$ = 50 kJ m⁻², $t_2$ = 183 MPa).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 20. Distribution of the triaxiality rate (about 4 mm propagation, $\Gamma_0$ = 50 kJ m⁻², $t_2$ = 183 MPa).</figcaption></figure>

## 6.1. Influence on the global behavior of the specimen

Fig. 18 and Table 4 show that the influence of the second slope is anything but negligible. The crack's initiation is not affected, but the propagation velocity is influenced by this parameter. The larger the slope, the higher the propagation velocity. This leads to a decrease in the global stiffness of the specimen. Besides, Fig. 17 illustrates the influence of the slope on the evolution of the toughness ($R$) during the propagation. For a fixed surface energy dissipation rate, the larger the slope of the traction-separation law, the lesser the energy rate dissipated in the bulk. This results in a decrease in the toughness. Besides, one can observe a decrease in the toughness during the crack's propagation. The smaller the slope of the traction-separation law, the sharper this decrease.

## 6.2. Local influence

The process zone is defined as the zone in which the separation is between $\delta_1$ and $\delta_c$. Fig. 19 shows the evolution of the traction load in the median position of the width of the cohesive zone. This evolution is plotted for a similar crack's progression of about 4 mm. Locally, one can observe a great variation in the distribution of the traction stress. Besides, the size of the fracture process zone also varies. The smaller the traction load $t_1$, the larger the process zone. One can see that the size of the
<figure><img src="image_22.png" /><figcaption>**Fig. 21.** Distribution of the average and equivalent stresses (about 4 mm propagation, $\Gamma_0 = 50 \text{ kJ m}^{-2}$, $t_2 = 183 \text{ MPa}$).</figcaption></figure>

fracture process zone is not negligible compared to the size of the crack. This can explain the great influence of the form of the traction-separation law on the crack's propagation and on the load-displacement curve.

Fig. 20 illustrates the influence of the slope on the triaxiality rate (calculated as in Section 5.4). One can see that the form of the law has no influence on the maximum value of the triaxiality rate. The triaxiality rate is small at the crack's tip, rises sharply upstream, then remains constant over the remaining length of the fracture process zone. Therefore, the constant part is longer for a small traction load $t_1$.

Fig. 21 shows the results in terms of average stress and equivalent stress. The average stress is highly dependent on the value of $t_1$ and follows the same evolution as the traction load calculated near the cohesive zone (Fig. 19). Indeed, it is the maximum principal stress at the bulk level which governs the traction load near the cohesive zone. The equivalent stress is maximum near the crack's front $\delta = \delta_c$, which explains that the triaxiality rate is small at that location. The maximum value of the equivalent stress is rather dependent on the form of the traction-separation law. However, in the rest of the fracture process zone, the influence of the slope is negligible.

## 7. Conclusion

This paper proposes a simulation of the tearing of a CT specimen at 1173 K. The crack is modeled using cohesive elements. The specificity of such an analysis lies in the viscoplastic nature of the material and the large strains induced by the temperature of the test. These specificities are taken into account in the cohesive elements at different stages.

The formulation of the elements is three-dimensional. We defined a local coordinate system based on the midplane calculated in the deformed configuration of the element [29]. Indeed, the choice of the coordinate system for the decomposition of the separation is very important. It defines the rupture mode of the element. Thus, due to the large strains, the coordinate system must follow the deformed configuration of the element and take into account the physical reality of the fracture phenomenon. Finally, the formulation takes into account the change in geometry induced by the large strains by calculating the residual and the stiffness matrix in the deformed configuration of the element and by taking the variation of its area into account.

The traction-separation law represents, in a macroscopic way, the local dissipative phenomena leading to the creation of the crack's lips. We used a traction-separation law of the Tvergaard and Hutchinson type [20], modified by introducing a hardening slope in place of the plateau. Indeed, the mechanical properties of the material strengthen as one gets closer to the crack's front because the strain rate increases due to the singularity. However, this phenomenon is lessened by the local damage near the singularity. This approach is validated by the fact that the separation velocity of the crack's lips, on average within the process zone, remained constant throughout the test.

The choice of the form of the traction-separation law is often arbitrary. Indeed, in the case of ductile fracture, it is commonly accepted that only the surface energy recovery rate and maximum traction parameters have an influence on the crack's propagation. In the case of our work, that assumption was not verified. We showed that the form of the traction law plays a significant role.

For given values of the surface energy recovery rate and maximum traction load, the stiffening slope influences the crack's propagation velocity. If one measures the toughness in the Turner and Kolednik sense [34], the larger the slope of the traction-separation law, the less tough the specimen. Therefore, the larger the slope, the higher the propagation velocity and the lower the overall stiffness of the specimen.
References

[1] OLIS, editor. OLHF seminar 2002, OECD, Madrid; 2002.

[2] Sehgal B et al. Assessment of reactor vessel integrity (ARVI). Nucl Engng Des 2005;235:213-32.

[3] Sehgal B et al. Assessment of reactor vessel integrity (ARVI). Nucl Engng Des 2003;221:23-53.

[4] Koundy V, Caroli C, Nicolas L, Matheron P, Gentzbittel JM, Coret M. Study of tearing behaviour of a PWR reactor pressure vessel lower head under severe accident loadings. Nucl Engng Des 2008;238:2411-9.

[5] ABAQUS. Version 6.6. Hibbit, Karlson & Sorensen Inc, Pawtucket, RI (USA).

[6] Barenblatt G. The mathematical theory of equilibrium cracks in brittle fracture. Adv Appl Mech 1962;7:749-58.

[7] Dugdale DS. Yielding of steel sheets containing slits. J Mech Phys Solids 1960;8:100-4.

[8] Irwin GR. Analysis of stresses and strains near the end of a crack traversing a plate. J Appl Mech 1957;24:361-4.

[9] Hillerborg A, Modeer M, Petersson P. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cem Concr Res 1976;6:773-82.

[10] Planas J, Guinea GV, Elices M. Size effect and inverse analysis in concrete fracture. Int J Fract 1999;95:367-78.

[11] Bazant Z. Concrete fracture models: testing and practice. Engng Fract Mech 2002;69:165-205.

[12] Rahul-Kumar P. Polymer interfacial fracture simulations using cohesive elements. Acta Mater 1999;47:4161-9.

[13] Allix O, LTvOque D, Perret L. Identification and forecast of delamination in composite laminates by an interlaminar interface model. Compos Sci Technol 1998;58:671-8.

[14] Tijssens M, van der Giessen E, Sluys L. Simulation of mode I crack growth in polymers by crazing. Int J Solids Struct 2000;37:7307-27.

[15] Foote RML, Mai Y-W, Cotterell B. Crack growth resistance curves in strain-softening materials. J Mech Phys Solids 1986;34:593-607.

[16] Camacho G, Ortiz M. Computational modelling of impact damage in brittle materials. Int J Solids Struct 1996;33:2899-938.

[17] Swanson PL, Fairbanks CJ, Lawn BR, Mai Y-W, Hockey BJ. Crack-interface grain bridging as a fracture resistance mechanism in ceramics: I. Experimental study on alumina. J Am Ceram Soc 2005;70:279-89.

[18] Needleman A. A continuum model for void nucleation by inclusion debonding. J Appl Mech 1987;54:525-31.

[19] Needleman A. An analysis of tensile decohesion along an interface. J Mech Phys Solids 1990;38:298-324.

[20] Tvergaard V, Hutchinson JW. The relation between crack growth resistance and fracture process parameters in elastic-plastic solids. J Mech Phys Solids 1992;6:1377-97.

[21] Tvergaard V, Hutchinson JW. Effect of strain dependent cohesive zone model on predictions of crack growth resistance. Int J Solids Struct 1996;33:3297-308.

[22] Siegmund T, Brocks W. A numerical study on the correlation between the work of separation and the dissipation rate in ductile fracture. Engng Fract Mech 2000;67:139-54.

[23] Tvergaard V. Predictions of mixed mode interface crack growth using a cohesive zone model for ductile fracture. J Mech Phys Solids 2004;52:925-40.

[24] Tvergaard V, Hutchinson JW. Two mechanisms of ductile fracture: void by void growth versus multiple void interaction. Int J Solids Struct 2002;39:3581-97.

[25] Correc A, Scheider I, Schwalbe KH. On the practical application of the cohesive model. Engng Fract Mech 2003;70:1963-87.

[26] Roy YA, Dodds Jr RH. Simulation of ductile crack growth in thin aluminium panels using 3D surface cohesive elements. Int J Fract 2001;110:21-45.

[27] Chen CR, Kolednik O, Heerens J, Fischer FD. Three dimensional modeling of ductile crack growth: cohesive zone parameters and crack tip triaxiality. Engng Fract Mech 2005;72:2072-94.

[28] Roychowdhury S, Roy YDA, Dodds Jr RH. Ductile tearing in thin aluminium panels: experiments and analyses using large displacement 3D surface cohesive elements. Engng Fract Mech 2002;69:983-1002.

[29] Ortiz M, Pandolfi A. Finite deformation irreversible cohesive elements for three dimensional crack propagation analysis. Int J Numer Meth Engng 1999;44:1267-82.

[30] Li H, Chandra N. Analysis of crack growth and crack tip plasticity in ductile materials using cohesive zone models. Int J Plast 2003;19:849-82.

[31] Jin Z-H, Sun CT. Cohesive fracture model based on necking. Int J Fract 2005;134:91-108.

[32] Schwalbe KH, Hellmann D. Application of the electrical potential method to crack length measurements using Johnson's formula. J Test Evaluat 1981;9:218-20.

[33] den Bosch MV, Schreurs P, Geers M. On the development of a 3D cohesive zone element in the presence of large deformations. Comput Mech 2008;42:171-80.

[34] Turner CE, Kolednik O. A micro and macro approach to the energy dissipation rate model of stable ductile crack growth. Fatigue Fract Engng Mater Struct 1994;17:1089-107.

[35] Griffith A. The phenomena of rupture and flow in solids. Philos Trans Roy Soc Lond 1920;221A:163-98.

[36] Stampfl J, Kolednik O. The separation of the fracture energy in metallic materials. Int J Fract 2000;101:321-45.

[37] Marie S. Approche TnergTtique de la dTchirure ductile. Ph.D. Thesis, UniversitT de PoitiT; 1992.

[38] Kassner M, Hayes T. Creep cavitations in metals. Int J Plast 2003;19:1715-48.

[39] Belytschko T, Liu WK, Moran B. Nonlinear finite elements for continua and structures. Wiley; 2000. p. 492-8 [chapter Element technology].

[40] Charlotte M, Laverne J, Marigo JJ. Initiation of cracks with cohesive force models: a variational approach. Eur J Mech Solids 2006;25:649-69.

[41] Gurson AL. Continuum theory of ductile rupture by void nucleation and growth - I. ASME J Engng Mater Technol 1977;99(2):2-15.