# Dynamic crack propagation under mixed-mode loading - comparison between experiments and X-FEM simulations

David Grégoire, Hubert Maigre, Julien Réthoré, Alain Combescure

► To cite this version:

David Grégoire, Hubert Maigre, Julien Réthoré, Alain Combescure. Dynamic crack propagation under mixed-mode loading - comparison between experiments and X-FEM simulations. International Journal of Solids and Structures, 2007, 44, pp.6517-6534. 10.1016/j.ijsolstr.2007.02.044. hal-00938608

HAL Id: hal-00938608

https://hal.science/hal-00938608

Submitted on 4 Jun 2021

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.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Dynamic crack propagation under mixed-mode loading – Comparison between experiments and X-FEM simulations

D. Grégoire, H. Maigre, J. Réthoré, A. Combescure *

Laboratory of Mechanics and Solid Structures, LaMCoS, INSA-LYON, CNRS UMR5259, F69621 Villeurbanne, France

The objective of this paper is to describe a simple dynamic crack propagation experiment which reproduces two phenomena: mixed-mode propagation and crack stop and restart. This experiment is explained and interpreted using X-FEM simulations. We show that a simple fracture theory which consists in using a dynamic crack initiation toughness, a crack orientation along the maximum principal stress and a simple equation for the calculation of the crack speed is sufficient to explain what is observed experimentally.

**Keywords:** Dynamic crack propagation; Elastic; X-FEM; Experiments; Split Hopkinson Pressure Bar; PMMA; Mixed-mode loading

## 1. Introduction

The objective of this paper is to propose a methodology for assessing dynamic crack propagation laws under mixed-mode loading. Unlike quasi-static cases where the loading and the crack position can be easily established, in dynamic impact cases the loading conditions, the variation of the propagation parameters and the exact position of the crack are difficult to control. Thus, the determination of relevant constitutive crack propagation laws from dynamic crack propagation experiments is a rather challenging operation. Therefore, if one wants to evaluate dynamic crack propagation laws under mixed loading, one must perform numerical simulations and assess the quality of the laws by comparing numerical simulations using these laws with experimental results. This process requires that the quality of both the numerical simulations and the experiments be perfectly controlled. The interest of the tests presented in this paper is that they provide a well-controlled experimental basis against which numerical simulations can be tested.

The eXtended Finite Element Method (X-FEM) (Moës et al., 1999; Dolbow et al., 2000; Belytschko et al., 2001; Réthoré et al., 2005) is chosen for the simulations because the cracks are not described explicitly by the

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: alain.combescure@insa-lyon.fr (A. Combescure).</footnote>
mesh. Furthermore, the implicit description of a crack geometry is compatible with any crack path, even if this path is a priori unknown. This is typically the case of mixed-mode loading in dynamics: the evolution of the crack orientation during the propagation is unpredictable. It was also proven and verified numerically by Réthoré et al. (2005) that contrary to any other existing numerical method the X-FEM does not introduce or dissipate numerical energy during crack propagation. The objectives of this paper are, first, to check the applicability of a well-known model and, second, to identify the dynamic fracture parameters of this model for polymethylmethacrylate (PMMA) through comparison of simulations with experimental test results.

PMMA is chosen for the dynamic crack propagation experiments because this classical material has brittle fracture behavior under dynamic loading. In these tests, various means if improving the reliability of the simulations are studied. The propagation laws depend upon three parameters: an initiation parameter which controls the propagation, a mixity parameter which enables the prediction of the direction of the propagation, and a velocity equation which gives the velocity of the crack tip. These laws are used for the interpretation of the experimental tests.

Besides the propagation history, the loading conditions are also measured. Direct measurements of stress intensity factors are not performed. Only purely mechanical quantities are measured and the identification of the fracture parameters of the material is achieved by fitting the numerical simulations to the experimental results.

## 2. Dynamic crack propagation laws

The crack propagation laws are chosen according to the global macroscopic concept of stress intensity factor (Irwin, 1957) extended to elastodynamics within the framework of (Bui, 1978; Freund, 1990).

Since PMMA is brittle, the fracture phenomenon is assumed to be governed by the intensity of the hoop stress $\sigma_{00}$ near the crack tip, which is evaluated using a hoop stress intensity factor $k_{00}$

$$k_{00} = \lim_{r \to 0} \sqrt{2\pi r \sigma_{00}} \quad (1)$$

where $(r, \theta)$ are the local polar coordinates of the crack tip.

The maximum hoop stress intensity factor and the corresponding local polar angle are denoted $K^*$ and $\theta^*$

$$K^* = \max_{\theta \in [-\pi, \pi]} k_{00} = k_{\theta^* \theta^*} \quad (2)$$

No propagation initiates as long as the maximum hoop stress intensity factor remains less than a critical value, the dynamic crack initiation toughness, $K_{1d}$. As soon as this threshold is reached, the direction of the maximum hoop stress defines the critical direction $\theta_c$ of the incipient propagation (Maigre and Rittel, 1993). The crack initiation criterion is

$$ \begin{aligned} K^* &< K_{1d} && (\text{no initiation}) \\ K^* &= K_{1d}, \quad \theta^* = \theta_c && (\text{initiation}) \end{aligned} \quad (3) $$

The dynamic crack initiation toughness is a material property and must be evaluated from experiments. It may depend on the temperature and the loading rate. Since experiments are made at room temperature and at constant loading rate (Table 2), their influences are considered the same for every experiment and the possible dependencies are not studied.

Following crack initiation, a different criterion must be determined in order to take the crack growth process into account. During the dynamic growth of a crack, the instantaneous maximum hoop stress intensity factor stays equal to the dynamic crack growth toughness $K_{1D}$, which can now depend on the velocity $\dot{a}$ of the crack tip

$$K^*(t, \dot{a}) \ge K_{1d} \Rightarrow K^*(t, \dot{a}) = K_{1D}(\dot{a}) \quad (\text{propagation}) \quad (4)$$

Again, possible dependencies on the loading rate will not be considered whereas the dependencies on the temperature are implicitly studied with the dependencies on the crack tip velocity and the dynamic crack growth toughness must be characterized through experiments. It is difficult to obtain a relevant estimate of $K_{1D}$ experimentally, but this quantity is necessary in numerical simulations in order to calculate the velocity of the
crack tip. Therefore, one uses a modified version of the expression given in Kanninen and Popelar (1985), in which the quasi-static toughness is replaced by the dynamic crack initiation toughness and the dynamic crack growth toughness is assumed to be

$$K_{1D}(\dot{a}) = \frac{K_{1d}}{1 - \left(\frac{\dot{a}}{c_R}\right)} \qquad (5)$$

where $c_R$, the velocity of the Rayleigh waves, is the theoretical maximum velocity of a crack in a homogeneous medium.

Some authors (as in Ravi-Chandar, 2004) introduce an additional dynamic crack arrest toughness in order to account for the observation that the dynamic crack growth toughness does not tend toward the dynamic crack initiation toughness when the crack tip comes to a stop. This additional parameter is not introduced here. The same critical value is chosen for crack initiation and crack arrest.

## 3. Experiments

### 3.1. Descriptions of the test rig and the specimen

The test rig is a Split Hopkinson Pressure Bar (SHPB) test developed by Kolsky (1949) and primarily used for the measurement of a material dynamic behavior. It is schematized in Fig. 1.

The test specimen is made of PMMA. The mechanical properties of this transparent, brittle material are summarized in Table 1. The bars must be chosen in order to guarantee a good waves transmission on their interfaces with the specimen and a measurable and elastic strain gage response. In our case, nylon is used for the bars because its mechanical properties are similar to those of PMMA. The material properties of nylon are summarized in Table 1.

There is no standard for dynamic fracture tests; therefore, we chose a simple, but distinctive, geometry designed to obtain separate values of the fracture criteria from the experiment. The specimen is shown in Fig. 2 and its dimensions are given in Fig. 3.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Characteristics of the Split Hopkinson Pressure Bar.</figcaption></figure>

<table><caption>Table 1
Mechanical properties of specimens and bars</caption><thead><tr><th>Mechanical properties</th><th>Symbols</th><th>(Units)</th><th>PMMA</th><th>nylon</th></tr></thead><tbody><tr><td>Density</td><td>ρ</td><td>(kg m<sup>-3</sup>)</td><td>1180</td><td>1145</td></tr><tr><td>Static Young's modulus</td><td>E<sub>S</sub></td><td>(GPa)</td><td>3.3</td><td>3.5</td></tr><tr><td>Dynamic Young's modulus</td><td>E<sub>D</sub></td><td>(GPa)</td><td>4.25</td><td>3.6</td></tr><tr><td>Poisson's ratio</td><td>ν</td><td>(-)</td><td>0.42</td><td>0.41</td></tr><tr><td>Lamé's constants</td><td>λ</td><td>(GPa)</td><td>7.8</td><td>5.8</td></tr><tr><td></td><td>μ</td><td>(GPa)</td><td>1.5</td><td>1.3</td></tr><tr><td>Velocity of compressive waves</td><td>c<sub>P</sub></td><td>(m s<sup>-1</sup>)</td><td>3032</td><td>2704</td></tr><tr><td>Velocity of shear waves</td><td>c<sub>S</sub></td><td>(m s<sup>-1</sup>)</td><td>1126</td><td>1056</td></tr><tr><td>Velocity of one-dimensional elastic waves</td><td>c<sub>L</sub></td><td>(m s<sup>-1</sup>)</td><td>1898</td><td>1773</td></tr><tr><td>Velocity of Rayleigh waves</td><td>c<sub>R</sub></td><td>(m s<sup>-1</sup>)</td><td>1064</td><td>996</td></tr><tr><td>One-dimensional impedance</td><td>√ρE</td><td>(10<sup>6</sup> kg m<sup>-2</sup> s<sup>-1</sup>)</td><td>2.25</td><td>2.03</td></tr></tbody></table>

(Karimzada and Maigre, 2000).
<figure><img src="image_3.png" /><figcaption>Fig. 2. Specimen with a hole.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Specimen geometry.</figcaption></figure>

In order to obtain fracture opening modes at the crack tip, the compressive waves must be converted into tensile waves. Therefore, we used a rectangular specimen with a circular hole to provide direct wave conversion. This configuration produces no friction at the specimen-bar interfaces, which would have generated perturbations in the gage signals upon impact of the specimen by a wedge (Karimzada and Maigre, 2000). In addition, crack arrest is provided by a circular hole as suggested by Karimzada and Maigre (2000): this interesting feature enables one to validate an arrest criterion.

In order to initiate the crack, an initial notch is machined at the border of the hole. Mixed-mode loading and crack orientation effects during the propagation are included by moving the initial notch upwards by a constant length $\Delta h = 7.5$ mm from the specimen axis of symmetry (Fig. 3).

The geometry was chosen to ensure a useful propagation zone large enough to produce interesting effects. In the mean time, the initial notch has to be short enough so that the crack is not subjected to boundary effects during most of its propagation (Fig. 2).

## 3.2. Measurements

The use of SHPB is attractive in our case because it provides both an accurate measurement of the applied loading and the global response of the test specimen during the transient experiment, thus enabling good control of the quality of the experimental tests. Reliable experimental data is necessary to ensure that the simulations are physically meaningful. This is a key to success in comparing numerical simulations with experiments.

Fig. 4 shows a sketch of the test rig. One can identify:

* the striker bar, input bar and output bar,

* the test specimen between the input and output bars,

* 3 strain gages connected to their amplifiers,

* 4 cameras Proxitronic HF-1 (640 × 288 pixels, 256 grey levels),
<figure><img src="image_5.png" /><figcaption>Fig. 4. Experimental test rig.</figcaption></figure>

* 1 standard flash unit,

* 1 delay line,

* 1 data acquisition adapter (4 channels, 1 MHz),

* 1 optical sensor connected to an oscilloscope,

* 2 computers.

One of our objectives is to obtain the history of the position of the crack tip. Since PMMA is transparent, optical tools are used. Four cameras with a very short exposure time (1 µs), connected through a delay line, provide four pictures (one picture per camera). The delay line time step adjustment is one microsecond. The results are checked to be highly reproducible. Therefore, by carrying out the same tests, repetitively, we are able to obtain precise paths and detailed histories. Different means of stress intensity factors measurement, like optical tools, are difficult to implement in this test rig because of lack of space. Therefore, only purely mechanical quantities are collected and stress intensity factors will be estimated through numerical simulations.

The loading is adjusted via several experimental parameters: the striker bar velocity, length and shape enabled us to control the amplitude, the duration and the shape of the loading.

Five experimental tests are performed by shifting the shooting time with the same specimen geometry under the experimental conditions described in Table 2. Unprocessed gage signals are shown in Fig. 5.

To limit wave superposition, the gages are generally located near the middle of the bars (Fig. 1). Consequently, the waves have to be shifted to the specimen-bar interfaces to obtain forces and velocities at the specimen faces (Fig. 6). Wave dispersion and geometry effects are taken into account in the shifting as in (Zhao and Gary, 1995) in order to obtain more accurate measurements.

### 3.3. Experimental results

The graphs shown in Fig. 6 confirm the reproducibility of the experimental tests. Since the duration of the propagation in a typical experiment is roughly 600 µs, the input and output velocities are drawn from 0 to 600 µs, as used for the numerical simulations.

<table><caption>Table 2
Output velocities of the striker bar for an air gun pressure of 1.2 bar</caption><thead><tr><th>Test</th><th>a</th><th>b</th><th>c</th><th>d</th><th>e</th></tr></thead><tbody><tr><td>Velocity (m s<sup>-1</sup>)</td><td>12.6</td><td>12.5</td><td>12.4</td><td>12.4</td><td>12.4</td></tr></tbody></table>

Uncertainty: 0.08 m s<sup>-1</sup>.
<figure><img src="image_6.png" /><figcaption>Fig. 5. Unprocessed gage signals.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Mechanical quantities at the specimen-bar interfaces.</figcaption></figure>

The input and output forces at the interfaces are fully represented in Fig. 6 to show that the loading agrees with the theoretical results of wave propagation through a cylindrical bar. The whole loading duration is approximately equal to twice the duration of the wave propagation through the length of the striker bar.

The crack evolutions for the five specimens are identical. Fig. 7 shows a post-mortem picture of a crack path. A crack arrest can be identified by the change in propagation direction. It has been checked that every

<figure><img src="image_8.png" /><figcaption>Fig. 7. Post-mortem view of a crack path.</figcaption></figure>
crack arrest in every post-mortem specimen appears in the same location. During the propagation, the crack turns towards the axis of symmetry (maximum stress zone).

The time when compressive waves reach the left-hand face of the specimen is chosen as the reference time. Fig. 8 shows a sequence of pictures of the cracked specimen taken by the cameras at the beginning of the propagation and before the crack arrest ($t = 216$ µs, $t = 226$ µs and $t = 246$ µs). There are only three pictures because one of the four cameras does not work sometimes. On each specimen, two marks are machined defining a reference length on each picture. In this way, a simple commercial software is needed to process the pictures and crack lengths are obtained by a rule of three.

The experimental crack length histories are collected in Fig. 9.

Their consistency inspires confidence in the quality and repeatability of the experiments.

Following the initiation at $t \approx 200$ µs, three different propagation phases can be easily observed:

* first, a propagation phase at constant horizontal velocity ($V_x \approx 211$ m s⁻¹), which corresponds to a maximum of the curvilinear velocity of the crack tip $\dot{a} \approx 260$ m s⁻¹,

* then, the crack stops for 50 µs,

* a second propagation stage at constant horizontal velocity 157 m/s ($\dot{a} = 160$ m/s) occurs until the final stop ($t \approx 500$ µs).

Graphs 6 and 9 reveal that initiation ($t \approx 200$ µs) occurs during a highly transient phase because the delay between the incident and transmitted signals is of the order of 100 µs. Therefore, the whole propagation

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

Uncertainty: Time: 1 µs Abscissa: 0.15 mm

<caption>Fig. 8. Pictures of Test *d* (time: 216, 226 and 246 µs).</caption>
<figure><img src="image_10.png" /><figcaption>Fig. 9. History of the X-coordinate of the crack tip.</figcaption></figure>

(initiation, arrest and restart until the final stop) is subject to transient stress waves evolutions in the specimen. These conditions provide a challenging test of the reliability of the numerical simulations.

Since there is a good impedance match between the output bar and the specimen, the crack arrest is certainly not due to wave reflection on the output interface. The study of different hole geometries (Karimzada and Maigre, 2000) indicates the fundamental influence of this parameter during the wave propagation and crack arrests.

# 4. Numerical simulations

## 4.1. The X-FEM method

In order to capture the discontinuity and singularity of the strain field at the crack tip, we use an enrichment to the classical finite element approximation defined through the partition of unity method developed in (Babuska and Melenk, 1997). This method, called X-FEM, was first developed for quasi-static analysis in (Moës et al., 1999). Here, the method is used for dynamic crack propagation as in (Réthoré et al., 2005).

Let us consider the two-dimensional problem of a homogeneous, elastic cracked body $\Omega$ in plane strain and assume that the only change in the material is due to crack propagation. Thus, the problem can be expressed in the framework of elastodynamics with small strains by adding an unknown $a(t)$, the position of the crack tip.

A classical finite element approximation of the displacement fields $\tilde{U}$ is

$$ \tilde{U}(t) = \sum_{i \in \mathcal{N}} N_i(x)u_i(t) \qquad (6) $$

where $\mathcal{N}$ is the set of the nodes which support the shape functions $\{N_i\}_{i \in \mathcal{N}}$ and $u_i(t)$ is the vector of the nodal degrees of freedom at time $t$.

The crack is represented by a set of straight segments. For the description of the crack to be independent of the mesh, discontinuous enrichments are added to the set of the nodes which have their support entirely cut by the crack, while the set of the nodes which contain the crack tip in their support is enriched by a singular set of functions, as shown in Fig. 10. When the crack propagates this set is modified as follows: all ancient enrichments are kept. New enrichments corresponding to the new position of the crack are simply added to the old set.
<figure><img src="image_11.png" /><figcaption>Fig. 10. Enrichment to the classical finite element approximation.</figcaption></figure>

The approximation of the displacement fields becomes

$$ \tilde{U}(t) = \sum_{i \in \mathcal{N}} N_i(x) u_i(t) + \sum_{i \in \mathcal{N}_{\text{cut}}(t)} N_i(x) \mathcal{H}(x) a_i^e(t) + \sum_{i \in \mathcal{N}_{\text{tip}}(t)} \sum_{j \in [1,4]} N_i(x) B_j^e b_{ij}^e(t) \quad (7) $$

where, at time $t$, $N_{\text{cut}}^e(t)$ is the set of the nodes which have had their support entirely cut by the crack, $a_i^e(t)$ are the nodal degrees of freedom corresponding to the step function $\mathcal{H}$, $N_{\text{tip}}^e(t)$ is the set of the nodes which have contained the crack tip in their support, $b_{ij}^e(t)$ are the nodal degrees of freedom corresponding to the functions $B_j^e$ which span the near-tip asymptotic fields.

According to Moës et al. (1999):

$$ \mathcal{H}(x) = \begin{cases} +1 & \text{if } x \text{ is above the crack} \\ -1 & \text{if } x \text{ is below the crack} \end{cases} \quad (8) $$

$$ \{B_j^e(r, \theta)\}_{j \in [1,4]} = \left\{ \sqrt{r} \cos \frac{\theta}{2}, \sqrt{r} \sin \frac{\theta}{2} \sin \theta, \sqrt{r} \sin \frac{\theta}{2}, \sqrt{r} \sin \frac{\theta}{2} \sin \theta \right\} \quad (9) $$

where $(r, \theta)$ are the local polar coordinates of the crack tip.

The displacements, velocities and accelerations are discretized through Eq. (7) and the mass and stiffness matrices are calculated using ad hoc integration techniques for X-FEM elements (Moës et al., 1999).

The discrete equation at time $t_n$ is

$$ M_n \ddot{U}_n + K_n U_n = F_n \quad (10) $$

where $M_n$ and $K_n$ are the mass and stiffness matrices at time $t_n$ and $F_n$ is the vector of external forces at time $t_n$.

## 4.2. Discretized energy conservation in the dynamic crack propagation case

During propagation, the new equilibrium at time $t_{n+1} = t_n + dt$ must be calculated by projecting the fields calculated at time $t_n$ with the geometry at time $t_n$ onto the new geometry at time $t_{n+1}$.

It was proved by Réthoré et al. (2005) that if all singular enrichments (Fig. 11) are kept and the new enrichments are initialized to zero (Eq. (11)), using a Newmark scheme, stability and exact conservation of the discretized energy are guaranteed. Hence, if $X_j^i$ is a field (displacement, velocity or acceleration) defined at time $t_j$ with the configuration at time $t_i$, there are two basic steps to be performed:
<figure><img src="image_12.png" /><figcaption>Fig. 11. Updating of the enrichment during propagation.</figcaption></figure>

$$ X_n^n \xrightarrow{\text{projection}} X_n^{n+1} \xrightarrow{\text{time stepping}} X_{n+1}^{n+1} \\ \text{With : } [X_n^{n+1}] = \begin{bmatrix} \frac{X_n^n}{0} \\ \vdots \\ 0 \end{bmatrix} \begin{array}{l} \text{old degrees of freedom} \\ \text{new degrees of freedom} \\ \text{with the crack growth} \end{array} \qquad (11) $$

The discretized dynamic Eqs. (10) and (14) are chosen such that they can be integrated numerically using Newmark's implicit mean acceleration scheme.

### 4.3. Boundary conditions

The numerical simulations are carried out using the input velocity collected from the experiments (Fig. 6) as a boundary condition at the input bar interface (Fig. 13).

The use of the experimental output velocity is difficult because the measurements are intrinsically imprecise. To achieve the same conditions during simulations and experiments would require perfect synchronization between input and output loading. Thus, one option would be to mesh the output bar entirely (Fig. 12): this would prevent the simulation of the specimen response from being affected by reflected waves, but at the expense of long calculation times.

Instead, we chose to model the output bar as in one-dimensional propagation. It leads the use of a simple impedance condition (Fig. 13).

At each interface node of the output bar interface, the stress and velocity are assumed to be related by

$$ \left\{ \begin{array}{l} \sigma.n = -z(v.n)n \quad \text{on } \partial\Omega_{\text{int}} \\ z = \rho^{\text{bar}}c_L^{\text{bar}} \end{array} \right. \qquad (12) $$

where $\sigma$ and $v$ are the stress and velocity on the specimen face, $z$ is the impedance, $\partial\Omega_{\text{int}}$ the interface and $n$ a normal vector, $\rho^{\text{bar}}$ and $c_L^{\text{bar}}$ are the density and the velocity of 1-D elastic waves in the bar.

<figure><img src="image_13.png" /><figcaption>Fig. 12. Case of a meshed output bar.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Fig. 13. Boundary conditions.</figcaption></figure>

Then, the impedance matrix is calculated by integrating Condition (12) at the interface

$$ (Z)_{ij} = \int_{\partial\Omega_{\text{int}}} z((N_i \cdot n)n) \cdot N_j dl \quad (13) $$

where Z is the impedance matrix, $N_i$ and $N_j$ are the shape functions of the displacement field (Eq. (7)).

Hence, the discrete form of Eq. (10) becomes

$$ M_n \ddot{U}_n + K_n U_n = F_n - Z_n \dot{U}_n \quad (14) $$

where $M_n$, $K_n$ and $Z_n$ are the mass, stiffness and impedance matrices at time $t_n$ and $F_n$ is the vector of external forces at time $t_n$.

The validity of this impedance condition was checked by numerical comparison with the calculation of a meshed output bar. The output velocity will be used further on to check the consistency between the experimental and numerical results (Fig. 21).

## 4.4. Discussion of the material properties of PMMA

PMMA is a viscoelastic material. Therefore, the loading rate affects its behavior and certainly has an influence on crack propagation. In compression, viscoelasticity tends to reduce the value of the dynamic Young's modulus as shown in Saad-Gouider et al. (2006). The same phenomenon appears in traction but in different proportion.

Our specimen is globally in compression, but the circular hole creates traction zones, which makes it difficult to characterize the mean effects of viscoelasticity and incorporate them into the model.

Finally, we retain elastic behavior, but choose a uniform dynamic Young's modulus taking into account the viscoelasticity of the material globally. Since the experimental tests provided input and output data, the global dynamic Young's modulus is chosen in order to have the correct delay between incident and transmitted signals in the simulations. A secant rather than dynamic Young's modulus is adopted.

The optimum value $E_D^{\text{secant}} = 2.4$ GPa is obtained and the resulting global mechanical properties of PMMA are shown in Table 3. These figures agree with the overall parameters of PMMA identified by Saad-Gouider et al. (2006).

## 4.5. Calculation of the fracture parameters

### 4.5.1. Dynamic stress intensity factors

In the vicinity of the crack tip, the dynamic stress intensity factors in Modes I and II are given (see Freund, 1990 for details) by

<table><caption>Table 3
The mechanical properties of PMMA chosen for the analysis</caption><thead><tr><th>ρ</th><th>1180</th><th>kg m<sup>-3</sup></th><th>λ</th><th>4.4</th><th>GPa</th><th>c<sub>S</sub></th><th>846</th><th>m s<sup>-1</sup></th></tr></thead><tbody><tr><th>E<sup>secant</sup><sub>D</sub></th><td>2.4</td><td>GPa</td><td>μ</td><td>0.8</td><td>GPa</td><td>c<sub>L</sub></td><td>1426</td><td>m s<sup>-1</sup></td></tr><tr><th>ν</th><td>0.42</td><td></td><td>c<sub>P</sub></td><td>2279</td><td>m s<sup>-1</sup></td><td>c<sub>R</sub></td><td>800</td><td>m s<sup>-1</sup></td></tr></tbody></table>
$$K_I^{\text{dyn}} = \lim_{r \to 0} \sqrt{2\pi r} \sigma_{22}(\theta=0) \quad \text{and} \quad K_{II}^{\text{dyn}} = \lim_{r \to 0} \sqrt{2\pi r} \sigma_{12}(\theta=0) \qquad (15)$$

where $(r, \theta)$ are the local coordinates of the crack tip, $(K_I^{\text{dyn}}, K_{II}^{\text{dyn}})$ are the dynamic stress intensity factors in Modes I and II.

A domain-independent integral $I^{\text{int}}$ (Réthoré et al., 2005; Réthoré, 2005) is used to calculate dynamic stress intensity factors. The mixed-mode separation is obtained via a two-field problem consisting of the actual fields $(u, \sigma)$ and the auxiliary fields $(u^{\text{aux}}, \sigma^{\text{aux}})$. Then, the Lagrangian conservation law for a virtual crack extension field $q$ leads to the following expression of $I^{\text{int}}$ (see Réthoré, 2005):

$$I^{\text{int}} = - \int_{A} q_{k,j} [(\sigma_{pq}^{\text{aux}} u_{p,q} - \rho \dot{u}_{p} \dot{u}_{p}^{\text{aux}}) \delta_{kj} - (\sigma_{ij}^{\text{aux}} u_{i,k} + \sigma_{ij} u_{i,k}^{\text{aux}})] dS \\ + 2 \int_{A} q_{k} [(\sigma_{ij,j}^{\text{aux}} u_{i,k} + \rho \ddot{u}_{i} u_{i,k}^{\text{aux}}) + (\rho \dot{u}_{i}^{\text{aux}} \dot{u}_{i,k} + \rho \dot{u}_{i} \dot{u}_{i,k}^{\text{aux}})] dS \qquad (16)$$

where $A$ is the area delimited by any contour $\partial A$ enclosing the crack tip, $q$ is pyramidal and compatible with the crack (tangent to the crack faces, with $\|q\| = 1$ at the tip and $\|q\| = 0$ on $\partial A$).

Then, a dynamic energy analysis provides an equivalent to Irwin's relation in plane strain

$$I^{\text{int}} = \frac{2(1 - v^2)}{E} (f_1(\dot{a}) K_I^{\text{dyn}} K_I^{\text{aux}} + f_2(\dot{a}) K_{II}^{\text{dyn}} K_{II}^{\text{aux}}) \qquad (17)$$

where $(K_I^{\text{aux}}, K_{II}^{\text{aux}})$ are the stress intensity factors of the auxiliary fields, $f_i$ are universal functions of the velocity of the crack tip

$$f_1(\dot{a}) = \frac{4\beta_1(1 - \beta_2^2)}{(\kappa + 1)D(\dot{a})} \quad \text{and} \quad f_2(\dot{a}) = \frac{4\beta_2(1 - \beta_2^2)}{(\kappa + 1)D(\dot{a})},$$

$$\beta_i = \sqrt{1 - \left(\frac{\dot{a}}{c_i}\right)^2} \quad \text{and} \quad D(\dot{a}) = 4\beta_1\beta_2 - (1 + \beta_2^2)^2,$$

$(c_1, c_2)$ are the velocities of the compressive and shear waves.

Eqs. (16) and (17) provide a value of the stress intensity factors $K_I^{\text{dyn}}$ and $K_{II}^{\text{dyn}}$ by choosing the Westergaard's exact asymptotic fields for auxiliary fields $u^{\text{aux}}$. $(K_I^{\text{dyn}}$ is evaluated through $(K_I^{\text{aux}}, K_{II}^{\text{aux}}) = (1, 0)$ and $K_{II}^{\text{dyn}}$ through $(K_I^{\text{aux}}, K_{II}^{\text{aux}}) = (0, 1))$.

In order to evaluate the stress intensity factors numerically, this interaction integral is calculated using a set of additional integration cells rich in Gauss points, called the $J$-domain (Gosz et al., 1998). This $J$-domain is independent of the mesh of the body and follows the crack tip and its orientation during the crack growth.

### 4.5.2. The crack growth criteria

A preferential direction $\theta^*$, defined as the direction of the maximum hoop stress as developed in Part 1 from Maigre and Rittel (1993), is assumed for the propagation of the crack. Thus, the preferential direction expression is obtained by solving

$$\frac{\partial \sigma_{00}}{\partial \theta} (\theta^*) = 0 \qquad (18)$$

where $(r, \theta)$ are the local polar coordinates and $\sigma_{00}$ and $\sigma_{r0}$ are the local components of the Westergaard asymptotic stress tensor.

In the initiation case ($\dot{a} = 0$), it leads to the equation

$$K_I^{\text{dyn}} \tan \left( \frac{\theta^*}{2} \right) + K_{II}^{\text{dyn}} \left( 1 - 2 \tan^2 \left( \frac{\theta^*}{2} \right) \right) = 0 \qquad (19)$$

and the analytical expression of this preferential direction is
$$ \theta^* = 2 \arctan \left( \frac{1}{4} \left[ \frac{K_I^{\text{dyn}}}{K_{II}^{\text{dyn}}} - \text{sign}(K_{II}^{\text{dyn}}) \sqrt{8 + \left( \frac{K_I^{\text{dyn}}}{K_{II}^{\text{dyn}}} \right)^2} \right] \right) \quad (20) $$

Then, the intensity of the loading near the crack tip is calculated using the maximum hoop stress intensity factor $K^*$ corresponding to $k_{\theta^*\theta^*}$, the hoop stress intensity factor in the preferential direction, Eq. (2)

$$ K^* = \cos^3 \frac{\theta^*}{2} <K_I^{\text{dyn}} > - \frac{3}{2} \cos \frac{\theta^*}{2} \sin \theta^* K_{II}^{\text{dyn}} \quad (21) $$

where $<K_I^{\text{dyn}}>$, the positive part of $K_I^{\text{dyn}}$, avoids any closure effect. Crack initiation occurs when this equivalent stress intensity factor $K^*$ reaches the dynamic crack initiation toughness $K_{1d}$.

Once the crack has been initiated, the propagating crack direction depends on the velocity of the crack tip (Eq. (4)) and there is no analytical solution of the nonlinear problem. Since the velocity of the crack tip is always small compared to that of Rayleigh waves and since substantial mixity takes place only at the beginning of the propagation, $\theta^*$ and $K^*$ are assumed to be given by Eqs. (20) and (21) throughout the propagation.

During the propagation, the velocity of the crack tip adapts itself in such a way that the equivalent stress intensity factor remains equal to the dynamic crack growth toughness (Eq. (4)). Hence, the velocity of the crack tip is obtained at each time step by solving Eq. (5)

$$ \dot{a} = \left( 1 - \frac{K_{1d}}{K^*} \right) * c_R \quad (22) $$

### 4.5.3. Dynamic crack initiation toughness: simulation with a fixed notch

The previous criteria provided the description of the dynamic propagation of the crack. All that remained to be found is an evaluation of the dynamic crack initiation toughness $K_{1d}$. The mixity at initiation is used to determine this material property through experiments. An initiation angle $\theta_c = -37, 5^\circ$ is measured on a post-mortem crack path (Fig. 14).

Then, the X-FEM numerical simulation of the dynamic response of the specimen with an initial fixed notch is performed. The mesh, shown in Fig. 13, consists of 1377 four-node elements with 4 integration points and the calculation requires 200 time steps with a step size chosen as $\Delta t = 1 \mu\text{s}$. Then, the evolutions of the maximum hoop stress direction $\theta^*$ and the maximum hoop stress intensity factor $K^*$ with time are calculated. The results are shown in Fig. 15.

The initiation angle ($37, 5^\circ$) determines the initiation time and the corresponding value of $K^*$, which is precisely the dynamic crack initiation toughness $K_{1d}$. For test $d$, the value of $K_{1d} = 1.47 \text{ MPa}\sqrt{\text{m}}$ is obtained with an uncertainty of $0.02 \text{ MPa}\sqrt{\text{m}}$ linked to the initiation angle uncertainty and the small gradient of the critical angle history (Fig. 15).

If the same method is applied to the five specimens, because of an initiation angle dispersion of $5^\circ$, the dynamic crack initiation toughness dispersion is $0.4 \text{ MPa}\sqrt{\text{m}}$.

Since simulations are performed with the loading corresponding to the test $d$, the dynamic crack initiation toughness value chosen for the subsequent simulations with propagation is $1.47 \pm 0.01 \text{ MPa}\sqrt{\text{m}}$.

<figure><img src="image_15.png" /><figcaption>Uncertainty: $0.25^\circ$

Fig. 14. Initiation angle.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. 15. Evaluation of $K_{1d}$.</figcaption></figure>

### 4.6. Simulations with propagation

The value of the dynamic crack initiation toughness $K_{1d}$ obtained above is used for the propagation criterion given by Eq. (22).

Simulations with propagation are performed using the same mesh (Fig. 13) with a total duration of 500 µs (100 steps with $\Delta t = 5$ µs).

Figs. 16 and 17 show, respectively, the comparison of the experimental and numerical results for the crack path and the history of the X-coordinate of the crack tip, taking into account the propagation of the crack ($K_{1d} = 1.47 \text{ MPa}\sqrt{\text{m}}$).

As shown in Fig. 16, while the crack path is correct in the initiation phase, the results are less satisfactory during the propagation phase. Moreover, the crack arrest obtained numerically occurs earlier and lasts longer than in the experiment.

<figure><img src="image_17.png" /><figcaption>Fig. 16. Comparison of crack paths, $K_{1d} = 1.47 \text{ MPa}\sqrt{\text{m}}$.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. 17. Comparison of crack tip X-coordinate histories, $K_{1d} = 1.47 \text{ MPa}\sqrt{\text{m}}$.</figcaption></figure>

Thus, the value $K_{1d} = 1.47 \text{ MPa}\sqrt{\text{m}}$ seems too high for the proper simulation of the crack propagation. A simple explanation is that the radius of the initial notch tip is larger than that of the propagating crack tip, as can be observed on a photograph of the propagating crack (Fig. 18). It has been shown by other means in Saad-Gouider et al. (2006) for PMMA (and also in (Akourri et al., 2000) for metals) that fracture toughness increases with an increase in notch radius.

Another explanation is that the crack initiation point is intrinsically not on the curve characterizing the dynamic crack growth criterion: $\lim_{\dot{a} \to 0} K_{1D}(\dot{a}) \neq K_{1d}$ in Eq. (5) as it has been noticed by Ravi-Chandar (2004).

Figs. 19–21 show, respectively, the comparison of the experimental and numerical results for the crack path, the history of the X-coordinate of the crack tip and the input and output velocities histories obtained with a lower value of the dynamic crack initiation toughness ($K_{1d} = 1.33 \text{ MPa}\sqrt{\text{m}}$).

On Fig. 19, there is a good matching of the crack paths and the crack arrests occurred at the same time and in the same location on Fig. 20. Thus, this lower value of the dynamic crack initiation toughness allows a

<figure><img src="image_19.png" /><figcaption>Fig. 18. Zoom on the initial notch and crack at initiation.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Fig. 19. Comparison of crack paths, $K_{1d} = 1.33 \text{ MPa}\sqrt{\text{m}}$.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 20. Comparison of crack tip X-coordinate histories, $K_{1d} = 1.33 \text{ MPa}\sqrt{\text{m}}$.</figcaption></figure>

better representation of the dynamic crack propagation. The propagation is slightly shorter in the experimental test, but this might be due to boundary effects which are not taken into account in the simulations.

On Fig. 21, oscillations of the numerical output velocity are noticeable. These are due to numerical noise in the calculations of the velocity and acceleration after propagation. Before the initiation, the numerical output velocity match the experimental output velocity. It is a check of the validity of the numerical impedance model of the output bar.

## 5. Conclusion

This work shows that a good combination of well-controlled experiments and refined X-FEM simulations enables one to explain the history of brittle dynamic crack propagation and arrest using simple ingredients.

The principal tensile stress model seems sufficient to explain the direction of the crack propagation, and the crack velocity seems to verify the set of equations proposed by Kanninen and Popelar (1985).
<figure><img src="image_22.png" /><figcaption>Fig. 21. Comparison of velocities histories, $K_{1d} = 1.33 \text{ MPa}\sqrt{\text{m}}$.</figcaption></figure>

Figs. 17 and 20 seem to indicate the need to distinguish between two different crack growth criterion to get an excellent prediction: a dynamic crack initiation toughness and a dynamic crack arrest toughness. The higher value ($K_{1d} = 1.47 \text{ MPa}\sqrt{\text{m}}$) seems a valid estimate to represent the initiation phase, whereas the lower value ($K_{1d} = 1.33 \text{ MPa}\sqrt{\text{m}}$) is better to represent the crack propagation and arrest.

Since the numerical model has only one parameter, the choice of the lower value is the best compromise to compute both the initiation and the propagation phases. If the numerical tool had two parameters (a dynamic crack initiation toughness and a dynamic crack arrest toughness), the value of the dynamic crack propagation could be perfectly simulated. Finally, the dynamic crack initiation toughness should be evaluated using the mixity at initiation as developed in the paper but it is necessary to carry out extra experimental tests in order to further evaluate the dependences of dynamic crack arrest toughness on temperature and loading rate.

## Acknowledgments

We are grateful to the *Laboratoire de Mécanique des Solides, École Polytechnique, France* where the experiments were performed.

This work was sponsored by ANR/CNRS Grant DYNRUPT(JC05_45254).

## References

* Akourri, O., Louah, M., Kifani, A., Gilgert, G., Pluvinage, G., 2000. The effect of notch radius on fracture toughness $j_{tc}$. *Engineering Fracture Mechanics* 65 (4), 491–505.

* Babuska, I., Melenk, J., 1997. The partition of unity method. *International Journal for Numerical Methods in Engineering* 40, 727–758.

* Belytschko, T., Moës, N., Usui, S., Parimi, C., 2001. Arbitrary discontinuities in finite elements. *International Journal for Numerical Methods in Engineering* 50, 993–1013.

* Bui, H.D., 1978. Mécanique de la rupture fragile. *Masson, Paris*.

* Dolbow, J., Moës, N., Belytschko, T., 2000. Discontinuous enrichment in finite elements with a partition of unity method. *International Journal for Numerical Methods in Engineering* 36, 235–260.

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

* Gosz, M., Dolbow, J., Moran, B., 1998. Domain integral formulation for stress intensity factor computation along curved three-dimensional interface cracks. *International Journal of Solids and Structures* 35 (15), 1763–1783.

* Irwin, G., 1957. Analysis of stresses and strains near the end of a crack traversing a plate. *Journal of Applied Mechanics* 24 (3), 361–364.

* Kanninen, M., Popelar, C.H., 1985. *Advanced Fracture Mechanics*. Oxford University Press, Oxford.

* Karimzada, T., Maigre, H., 2000. Modelisation of dynamic crack propagation criteria. In: *Proceedings of the 6th International Conference on Mechanical and Physical Behaviour of Materials under Dynamic Loading*, vol. 10. Krakow, Poland.

* Kolsky, H., 1949. An investigation of the mechanical properties of material at a very high rate of loading. In: *Proceedings of the Physical Society, vol. B* 62. pp. 676–700.
Maigre, H., Rittel, D., 1993. Mixed-mode quantification for dynamic fracture initiation: application to the compact compression specimen. *International Journal of Solids and Structures* 30 (23), 3233–3244.

Moës, N., Dolbow, J., Belytschko, T., 1999. A finite element method for crack growth without remeshing. *International Journal for Numerical Methods in Engineering* 46, 131–150.

Ravi-Chandar, K., 2004. *Dynamic Fracture*. Elsevier, Amsterdam.

Réthoré, J., 2005. *Méthodes éléments finis étendus en espace et en temps: application à la propagation dynamique des fissures*. Ph.D. thesis, Insa de Lyon.

Réthoré, J., Gravouil, A., Combescure, A., 2005. An energy conserving scheme for dynamic crack growth with the extended finite element method. *International Journal for Numerical Methods in Engineering* 63, 631–659.

Saad-Gouider, N., Estevez, R., Olagnon, C., Séguéla, R., 2006. Calibration of a viscoplastic cohesive zone for crazing in PMMA. *Engineering Fracture Mechanics* 73, 2503–2522.

Zhao, H., Gary, G., 1995. A three dimensional analytical solution of the longitudinal wave propagation in an infinite linear viscoelastic cylindrical bar. *Application to experimental techniques*. *Journal of Mechanics and Physics of Solids* 43 (8), 1335–1348.