# Mode-I crack propagation using a cohesive zone with rate-sensitivity

Nunziante Valoroso, Gilles Debruyne, Jérôme Laverne

► To cite this version:

Nunziante Valoroso, Gilles Debruyne, Jérôme Laverne. Mode-I crack propagation using a cohesive zone with rate-sensitivity. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01717092

HAL Id: hal-01717092

https://hal.science/hal-01717092

Submitted on 25 Feb 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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.

Public Domain
# Mode-I crack propagation using a cohesive zone with rate-sensitivity

Nunziante VALOROSO <fnref n="1" />,<fnref n="2" />, Gilles DEBRUYNE <fnref n="2" />, Jérôme LAVERNE <fnref n="2" />

**Abstract** — We study mode-I dynamic crack propagation using the cohesive-zone methodology. Rate-dependency is included in the model in a way that is directly related to the speed of the crack front that is typically measured during experiments. A representative numerical example is presented that demonstrates the effectiveness of the proposed approach.

**Keywords** — Dynamic fracture, cohesive-zone models, rate-dependency

## 1 Introduction

The subject of dynamic fracture has received increasing attention in recent years owing to its relevance in a variety of industrial applications where crack initiation cannot be precluded [1]. In this context design and verification of safe crack arrest in structures is the primary risk management strategy against unwanted and possibly catastrophic events.

Tackling the problem of dynamic fracture via the classical cohesive-zone approach requires a special care. This is mainly due to the fact that in dynamic fracture additional dissipative mechanisms can manifest that, if not properly accounted for, prevent from obtaining accurate numerical results. One of the most common undesired effects encountered when simulating rapidly propagating cracks is the very high speed of the crack front that is likely to be estimated in numerical computations, where velocities up to the Rayleigh wave speed of the material can be attained. This last one has to be understood as the upper bound speed for cracks propagating in an ideal elastic continuum [2] since in real situations the experimentally recorded crack velocities for mode-I fracture usually do not exceed 60% of the theoretical limit even in most brittle materials, see for instance [6] and references therein.

In the cohesive zone methodology, stemming from ideas initially proposed by Dugdale and Barenblatt, the description of the whole fracture process is provided by the cohesive law, i.e. a relationship between surface tractions and displacements discontinuities. Therefore, the cohesive law embodies all complex dissipative processes leading to the progressive decay of cohesive forces preceding the formation of traction-free surfaces and crack propagation. For situations where potential fracture trajectories are a priori known, cohesive models are used in conjunction with zero-thickness interface elements whose size has to be chosen in a way to adequately resolve the process zone.

Recent contributions, see e.g. [3, 4, 11] among others, have shown that use of classical rate-independent cohesive models to simulate dynamic fracture can produce unrealistic answers. Conceptually this is not surprising since dynamic fracture phenomena typically occur over a quite short time scale, whereby some form of rate-dependency at the crack tip is reasonably expectable. Basically, two different approaches have been used in the literature to account for this rate-dependency, that is either by using the classical cohesive model in conjunction with a rate-dependent constitution for the bulk material, either introducing rate-sensitivity directly into the cohesive law.

The objective of the present work is twofold. On the one hand we briefly revisit the problem of dynamic fracture to set the stage for the subsequent implementation of the cohesive model. On the other side, we introduce a basic form of rate-dependency into the cohesive relationship initially developed in [8] in a way to make it adjustable in run-time to account for variations in the dissipation power with the speed of the running crack. This requires in turn a suitable modification of the traction-separation law in which one includes, though in a rather implicit form, those additional dissipation mechanisms that come into the picture with kinetic energy.

<footnote n="1">¹ Università di Napoli Parthenope, Dipartimento per le Tecnologie, nunziante.valoroso@uniparthenope.it</footnote>

<footnote n="2">² EDF R&D, LaMSID - UMR EDF/CNRS/CEA 8193, {gilles.debruyne, jerome.laverne}@edf.fr</footnote>
In particular, one can show that, even in a simplified phenomenological model, to capture the essen-
tial information the traction-separation relationship has not to be just uniformly expanded or contracted
with varying velocity. On the contrary, in order to obtain a sufficiently smooth response, the velocity
toughening effect should be uniform along the interface. Numerical results and comparisons with true
experimental data are presented that also provide useful hints in view of the design and application of
advanced rate-dependent cohesive models to large-scale industrial simulations.

## 2 Cohesive damage model

The formulation of a cohesive-like model basically requires the definition of a relationship between
interface tractions and displacement discontinuities and the introduction of a criterion for damage to grow
and the process zone to advance. In particular, in the following we briefly recall the model proposed in [8]
and consider the simplest two-parameter version for the mode I case, that is governed by the following
equations :

$$
\begin{align}
t &= (1-D)k \langle [[u]] \rangle_+ + k^- \langle [[u]] \rangle_- \nonumber \\
Y &= \frac{1}{2} k \langle [[u]] \rangle_+^2 \nonumber \\
\phi &= Y - F(D) \le 0 \nonumber \\
\phi &\le 0; \quad \dot{D} \ge 0; \quad \dot{D} \phi = 0 \tag{1}
\end{align}
$$

where $t$ and $Y$ respectively denote the interface traction and the damage-driving force, $D \in [0, 1]$ is the scalar damage variable, $\llbracket u \rrbracket$ is the displacement jump in the direction normal to the interface while $k$ and $k^-$ are the undamaged interface stiffnesses in tension and compression, respectively. The impenetrability constraint is introduced in penalty form via the stiffness coefficient $k^-$ and by distinguishing between the positive $\langle \cdot \rangle_+$ and negative part $\langle \cdot \rangle_-$ of the displacement jump.

The positive function $F(D)$ represents the critical damage-driving force, i.e. a monotone non-decrea-
sing energy threshold whose value depends upon the current damage state. Typical forms of $F$ are power
laws or exponential functions, and their explicit expressions are constructed in a way to ensure that the
energy dissipated in the formation of a new unit traction-free surface equals the fracture energy $G_{Ic}$ of
the material. In particular, the two-parameter version of the exponential traction-separation relationship
presented in [8] is obtained using the following expression for the damage function :

$$
F(D) = -G_{Ic} \log(1 - D) \qquad (2)
$$

We recall that in the present model the computation of the damage state, is completely explicit ; in
particular, for damage loading ( ˙D > 0) at each time τ the damage variable is computed as :

$$
D(\tau) = \min \left( 1, \max_{(\sigma \le \tau)} \{F^{-1}(Y(\sigma))\} \right) \quad (3)
$$

see also [9] for a more detailed account.

3 A look at dynamic fracture

Dynamic fracture mechanics has been initially developed based on two fundamental hypotheses. The first, more explicit one, is that the material is brittle, i.e. linear elastic up to fracture. This statement is more properly traduced by the small-scale yielding (SSI) assumption, whereby all the dissipative processes accompanying the formation of stress-free surfaces occur in a region around the crack tip (the so-called process zone) which is small in size compared to the crack extension and the specimen dimensions and is dominated by the square-root singularity. The energy cost to pay for a new unit surface to form is the fracture energy and it collects contributions from all nonlinear mechanisms taking place in the process zone. The second assumption, which is so obvious that is seldom explicitly stated, is that there is only one single active crack and that so it remains during the entire process. In other words, no secondary or subscale fracturing process is admitted.
The fundamental reference on the subject of dynamic fracture is the monograph by LB Freund [2], where a comprehensive theory is developed. In particular, under the above assumptions the equation of motion for a straight semi-infinite crack in an elastic infinite medium is derived as :

$$ \dot{a} = c_R \left( 1 - \frac{G_{Ic}}{G(a)} \right) \qquad (4) $$

where a and $\dot{a}$ respectively denote the crack length and its velocity, G is the energy release rate and $c_R$ is the Rayleigh wave speed of the material, which represents the maximum theoretical speed of a propagating crack due to obvious constraints on energy transport by elastic waves.

Limitations to the use of the above equation are well-known and have been put forward by experiments. In particular, two effects are usually recorded, i.e. crack speeds that usually do not exceed 60% of the theoretical limit even in most brittle materials, see e.g. [6] and references therein, and significant rise of material toughness with increasing velocity of the crack front.

These apparent inconsistencies are the outcome of the hypothesis of absence of subscale processes, i.e. secondary cracks, which are usually observed for crack speeds exceeding one third of $c_R$, a limit beyond which eq. (4) is no longer valid. From the phenomenological standpoint the apparent raise of energy cost for crack progression is well explained by the appearance of secondary cracks aside of the main pre-existing one when crack speed increases beyond a critical threshold. Actually, the appearance of micro-cracks corresponds to an increase of stress-free surfaces, whereby a fracture energy rise with velocity can be introduced to model in a simplified manner all those secondary processes that are neglected in the quasi-static regime.

In the following we shall outline a rate-dependent model based on the hypothesis that the fracture energy is still a material property but no longer independent from the velocity of the propagating crack.

<figure><img src="image_1.png" /><figcaption>Figure 1 – Left : amplification function for increasing rate-sensitivity, see eqn. (5). Right : two possible modifications of the initial traction-separation relationship (red). Uniform expansion (green) and $G_{Ic}$ amplification at fixed strength (brown).</figcaption></figure>

## 3.1 A rate-dependent model

Based on the discussion of the previous section, we introduce in the cohesive model a dependence of the fracture energy from the macroscopic crack velocity that is typically measured during mechanical tests in order to account for experimentally observed behaviour at high crack speeds. Moreover, with a view towards subsequent parameter identification we keep this functional dependence as simple as possible and introduce only one additional material constant. Accordingly, the rate effects is described in the form of an amplification function via the following expression, see also Figure 1 :

$$ \frac{G_{Ic}^{dyn}(\dot{a})}{G_{Ic}} = 1 + \eta H(m) \qquad (5) $$
where $G_{Ic}^{dyn}$ denotes the speed-dependent dynamic fracture energy, $\eta$ is a non-dimensional rate-sensitivity parameter, $H$ is a monotone-non decreasing positive function and $m$ is the crack Mach number:

$$m = \frac{\dot{a}}{c_H} \qquad (6)$$

$c_H$ being a critical velocity value which is understood as a threshold beyond which the role of the amplification function becomes significant.

Equation (5) has been formulated in a way to depend from the macroscopic instantaneous crack velocity $\dot{a}$, which is a global measurable datum. In order to ease the implementation it can be useful to recall the relationship:

$$\dot{a} = - \left[ \left[ \dot{u} \right] \left( \frac{\partial \left[ u \right]}{\partial x} \right) \right]^{-1} \qquad (7)$$

which strictly holds under steady-state conditions since in this case the material time derivative of any quantity associated with a material point whose placement is $x$ is nihil.

In the adopted implementation the amplification function (5) is used to modify in run-time the damage mode (2) that characterizes the cohesive relationship. However, depending on the cohesive model and/or on the parametrization being used, an increment in the fracture energy $G_{Ic}$ can also result in a rise of the cohesive strength $t_{max}$. As shown in [9], this is indeed the case for the model described in Section 2 since these two material parameters are not independent.

In the present context a uniform amplification of the cohesive relationship seems to be inadequate in consideration of the physics of the phenomenon to be modeled, which manifests itself through a spread of the process zone due to the appearance of micro-cracks. Contrariwise, this effect is well described by allowing an increase in size of the cohesive zone with a limited increase (or no change at all) of strength. A pictorial representation of the two situations is given in Figure 1.

The velocity-dependent changes of the cohesive relationship in the present implementation freeze the peak stress and lower the initial stiffness in order to allow an increase of the fracture energy. It bears emphasis that at the present stage this change in the cohesive relationship is chosen to be the same for all points falling into the cohesive zone and is also kept constant within each load step in order to avoid spurious oscillations of the computed response. Conversely, the crack speed that is used to change the fracture energy can be computed in various ways, i.e. either as an average velocity of the cohesive zone, or as the velocity of a particular point that is conventionally designed as the crack tip, without significantly affecting the results.

## 4 Numerical example

In order to show the capabilities of the rate-sensitive cohesive model we consider as an example a Compact Tension (CT) test on a 18MND5 steel alloy for which experimental data on dynamic crack propagation are available from [10].

<figure><img src="image_2.png" /><figcaption>Figure 2 – Compact tension specimen. Geometry and experimental test (after [5]).</figcaption></figure>
The samples were first pre-cracked by fatigue to create a sharp crack starting from the preexisting notch. The initial length of the crack and the load was then estimated from the measured specimen compliance. Afterwards the specimens were cooled by nitrogen injection in an insulating enclosure up to the target temperature (–150°C for the case of interest) [5]. Loading was then applied via displacement control. In the above conditions the bulk material exhibits a very brittle behaviour and its constitution before fracture occurs is well described by linear elasticity. Once the crack is initiated it propagates abruptly to reach its final length in less than 100 microseconds [10]. Crack length and speed were recorded using electric wires bonded on both sides of the specimen, see Figure 2.

<table><caption>Table 1 – Compact Tension specimen. Material data</caption><tr><td>Bulk material</td><td>ρ = 7.80 (g/cm³)</td><td>E = 2.11e+05 (MPa)</td><td>ν = 0.3</td><td>c<sub>R</sub> ≈ 2960 m/s</td></tr><tr><td>Interface</td><td>G<sub>Ic</sub> = 1.85 (N/mm)</td><td>t<sub>max</sub> = 2.10e+03 (MPa)</td><td>η ∈ [0, 1.2]</td><td></td></tr></table>

The rate-dependent cohesive law described in previous sections has been implemented in a customized version of the finite element code FEAP [7] as the material constitution for zero-thickness cohesive elements. The mechanical test on the CT samples has been indeed modeled using a simplified geometrical and mechanical scheme by defining a priori the crack trajectory which is meshed using 4-noded interface elements whereas EAS quadrilaterals are used for the bulk material. The data set is summarized in Table 1.

In numerical simulations loading is prescribed via a two-stage process. In a first phase a vertical displacement is imposed in quasi-static manner up to a pre-determined value by keeping interface elements fully constrained so to inhibit crack propagation. Afterwards, the interface elements are released and the transient solution is computed using the implicit constant average acceleration Newmark integrator.

<figure><img src="image_3.png" /><figcaption>Figure 3 – Compact tension specimen. Experimental VS computed crack position (left) and speed (right). Rate-independent model.</figcaption></figure>

Simulations have been first carried out using the rate-independent model to ascertain its accuracy compared to experimental data. The main results of the computation are summarized in Figure 3, where are depicted the crack progression and the crack velocities curves. It is clearly seen that the agreement between numerical and experimental results is loosely qualitative since the simulation overpredicts both the crack length at arrest and the propagation speed. In particular, this last one is initially estimated to approach 2000 m/s, which is about twice the experimental value. This can well be attributed to the fact that in the rate-independent model all subscale fracture mechanisms are neglected by taking as the fracture energy of the material the (constant) value obtained in a quasi-static test.

The global beneficial effects subsequent to the introduction of rate-sensitivity can be appreciated from Figure 4, which refers to the case $η = 1.2$. As expected, the effect of rate-sensitivity is that of slowing down crack propagation and of reducing the distance ran by the crack.
In this case agreement with experimental data is significantly improved though the computed curves are somehow stiffer than the experimental ones. Moreover, a careful inspection of the experimental curves, which are not exactly coincident on the two faces of the specimen, suggests that the crack front is neither planar nor normal to the specimen surface. Hence, apart from other approximations (such as the straight crack path), the assumption of plane strain conditions that is used throughout computations does probably affect the results as well. Use of a fully three-dimensional model is therefore expectable to further mitigate differences between experimental and numerics.

<figure><img src="image_4.png" /><figcaption>Figure 4 – Compact tension specimen. Experimental VS computed crack position (left) and speed (right). Rate-dependent model.</figcaption></figure>

References

[1] G. Debruyne, J. Laverne, P.E. Dumouchel. *Dynamic crack growth: Analytical and numerical cohesive zone models approaches from basic tests to industrial structures*, Engineering Fracture Mechanics, 90, 1–29, 2012.

[2] L.B. Freund. *Dynamic Fracture Mechanics*. Cambridge University Press, 1990.

[3] N. Murphy, A. Ivankovic. *The prediction of dynamic fracture evolution in PMMA using a cohesive zone model*, Engineering Fracture Mechanics, 72, 861–875, 2005.

[4] A. Pandolfi, P.R. Guduru, M. Ortiz, A.J. Rosakis. *Three dimensional cohesive-element analysis and experiments of dynamic fracture in C300 steel*, International Journal of Solids and Structures, 37(27), 3733–3760, 2000.

[5] B. Prabel. *Modélisation avec la méthode X-FEM de la propagation dynamique et de l'arrêt de fissure de clivage dans un acier de cuve REP*, PhD Thesis, INSA, Lyon, 2007.

[6] K. Ravi-Chandar. *Dynamic fracture of nominally brittle materials*, International Journal of Fracture, 90, 83–102, 1998.

[7] R.L. Taylor. *FEAP - Programmer Manual*. University of California at Berkeley, http://www.ce.berkeley.edu/~rlt, 2010.

[8] N. Valoroso, L. Champaney. *A damage-mechanics-based approach for modelling decohesion in adhesively bonded assemblies*, Engineering Fracture Mechanics, 73(18), 2774–2801, 2006.

[9] N. Valoroso, R. Fedele. *Characterization of a cohesive-zone model describing damage and de-cohesion at bonded interfaces. sensitivity analysis and mode-I parameter identification*, International Journal of Solids and Structures, 47(13), 1666–1677, 2010.

[10] T. Yuritzinn, T. Le Grasse. *Essais isothermes de propagation dynamique de fissure sur CT et anneaux en acier 18MND5*», Rapport interne CEA n° SEMT/LISN/RT/08-023/A, 2008.

[11] F. Zhou, J.F. Molinari, T. Shioya. *A rate-dependent cohesive model for simulating dynamic crack propagation in brittle materials*, Engineering Fracture Mechanics, 72, 1383–1410, 2005.