# Cohesive model for fracture simulation of uranium dioxide at high temperature and strain rate

Clément Tosi, Rodrigue Largenton, Yann Monerie, Pierre-Guy Vincent

► To cite this version:

Clément Tosi, Rodrigue Largenton, Yann Monerie, Pierre-Guy Vincent. Cohesive model for fracture simulation of uranium dioxide at high temperature and strain rate. 2024. hal-04965661v1

HAL Id: hal-04965661

https://hal.science/hal-04965661v1

Preprint submitted on 7 Apr 2025 (v1), last revised 14 Jan 2026 (v2)

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
Cohesive model for fracture simulation of uranium
dioxide at high temperature and strain rate

Clément Tosi¹, Rodrigue Largenton², Yann Monerie³, and Pierre-Guy Vincent¹

¹ Institut de Radioprotection et de Sûreté Nucléaire (IRSN), PSN-RES/SEMIA/LSMA, B.P. 3, 13115 Saint-Paul-lez-Durance Cedex, France

² Electricité de France, EDF R&D Dpt MMC, Avenue des Renardières, Moret-Loing et Orvanne, Ecuelles, 77250, France

³ Univ. Montpellier, Laboratoire de Mécanique et Génie Civil (LMGC), CNRS, CC 048 Place Eugène Bataillon, 34095 Montpellier Cedex, France

A model known as the cohesive zone model is employed to accurately simulate the initiation and propagation of cracks within uranium dioxide fuel under high-temperature and high-strain-rate conditions. This model incorporates two energy-based failure criteria to distinguish between normal and tangential modes of crack propagation. Furthermore, it explicitly takes into account the dependencies on temperature, strain rate, and porosity. The effectiveness of this cohesive zone model approach is validated through three-dimensional simulations of a uniaxial compression test conducted under high-temperature and high-strain-rate conditions. Notably, the critical cohesive energy value identified for the normal mode of crack propagation is found to be lower than that for the tangential mode of crack propagation.

Keywords uranium dioxide, viscoplasticity, fracture, cohesive zone model, crack growth simulation, compressive test simulation

# 1 Introduction

Uranium dioxide and mixed oxide (MOX) fuels are extensively used in French Pressurized Water Reactors. Continual research is dedicated to understanding the mechanical behavior of these nuclear fuels, crucial for assessing reactor safety, particularly under design-based accident scenarios such as Reactivity Initiated Accidents (RIA) (Koo et al. 1997; Schmitz et al. 1999; Fuketa 2012). Our focus lies on the fracture modeling of uranium dioxide fuel under conditions of high temperature and high strain rate, representative of those encountered during a postulated RIA event. Above a specific temperature, uranium dioxide exhibits an elasto-viscoplastic behavior which depends on both temperature and strain rate (Labat et al. 2023).

On one hand, Salvo et al. (2015) recently derived a constitutive model for uranium dioxide creep which is suitable for a broad range of strain rates (10⁻⁹ - 10⁻¹/s) and temperatures (1100 - 1700°C). A hyperbolic sine model which accounts for porosity was proposed. This poro-viscoplastic model constitutes a damage model since the porosity can increase during the loading, inducing a softening behavior of the material. This model results from a comprehensive experimental campaign involving uniaxial compression tests on uranium dioxide fuel pellets conducted at high temperatures and strain rates, complemented by Scanning Electron Microscopy analyses, Finite Element simulations of these compression tests, and comparison with a database comprising over two hundred creep tests.

On the other hand, within the domain of fracture modeling for various polycrystalline materials encountered in nuclear reactors, that exhibit a ductile behavior at high temperatures, Nkoumbou Kaptchouang et al. (2021) proposed a cohesive-volumetric approach. This approach was developed specifically to simulate crack initiation and propagation in ductile materials experiencing void nucleation, growth, and coalescence, particularly in regions of the structure where plastic strain accumulates. The overall elastoplastic damage material response results in the coupling of both a hardening bulk constitutive law and a softening surface traction-separation law, effectively constituting a cohesive zone model, embedded between each mesh element of the finite element discretization.

The phenomenological approach of fracture through cohesive zone models is due to the pioneering works of Dugdale (1960) and Barenblatt (1959). It is now widely used to simulate
crack initiation and propagation in homogeneous or heterogeneous materials (see Schwalbe et al. (2012) among others). In this framework, the modeling introduces a process zone near the crack tip, avoiding the unrealistic unbounded stresses predicted by Linear Elastic Fracture Mechanics. Consequently, the stress remains finite at the crack tip and tends to vanish while the crack tip is opening. As delineated in Nkoumbou Kaptchouang et al. (2021), crack growth, within a standard cohesive surface theory of fracture, is modeled as the decohesion of two surfaces within a continuum material and the cohesive stress corresponds to the amount of stress which is needed to initiate the separation of the two surfaces. In such a modeling, the overall material behavior is governed by two constitutive relations, the bulk stress-strain constitutive law and the traction-separation law across the cohesive surface embedded in a three-dimensional continuum. From the general idea introduced by Tvergaard et al. (1992) in the context of ductile fracture, and then extended by Siegmund et al. (1999); Siegmund et al. (2000); Scheider (2009); Anvari et al. (2006); Anvari et al. (2007), the traction-separation law should represent the failure behavior of a unit material cell. In this context, Nkoumbou Kaptchouang et al. (2021) proposed to derive a cohesive zone model for ductile materials from a micromechanical approach. The proposed cohesive zone model of Nkoumbou Kaptchouang et al. (2021) is based on the widespread micromechanical Gurson-Tvergaard-Needleman (GTN) model (Gurson 1977; Koplik et al. 1988; Tvergaard 1990) and adapted to the reduced kinematics of a surface through the application of a scaling rule Suquet (1988); Michel et al. (1994). Although initially formulated within the specific framework of ductile damage and the GTN model, the approach outlined by Nkoumbou Kaptchouang et al. (2021) is quite general and can be used for other damage processes and models.

In this study, we adopt the methodology proposed by Nkoumbou Kaptchouang et al. (2021), departing from the viscoplastic damage model of Salvo et al. (2015) dedicated to uranium dioxide fuel, instead of the GTN plastic damage model. The resulting cohesive zone model is thus fully dedicated to uranium dioxide fuel and explicitly depends on temperature, strain rate, and porosity.

The structure of this paper is organized as follows. The second section deals with the modeling and the derivation of the cohesive zone model. The third section is devoted to the simulation of a uniaxial compression test of a uranium dioxide pellet up to fracture and its comparison to available experimental results.

## 2 Modeling

The fracture of uranium dioxide under high temperature and strain rate conditions is described by a cohesive zone model derived from the volume-to-surface approach of Nkoumbou Kaptchouang et al. (2021) applied to the bulk viscoplastic damage model of Salvo et al. (2015).

### 2.1 Reduced kinematics of a cohesive zone model

In the framework proposed by Nkoumbou Kaptchouang et al. (2021), it was assumed that the failure of a material domain results from the strain localization within a narrow band, where plastic strain accumulates (illustrated in Figure 1). The thickness of this narrow band is denoted by $h$ and is assumed to be significantly smaller than a characteristic length scale of the material domain. This narrow band is thus summarized as a zero thickness cohesive surface $S$: the material domain is composed of two regions, $\Omega^-$ and $\Omega^+$, separated by this cohesive surface $S$. The displacement jump across $S$ is defined as:

$$[u] = u^+ - u^- = [u_n]n + [u_t]t + [u_s]s, \qquad (1)$$

where $(n, t, s)$ is the local normal-tangent basis with $n$ the unit normal vector to the surface $S$, and $t, s$ two unit vectors defining the tangent plane. The respective values of the displacement fields at the boundaries of $\Omega^+$ and $\Omega^-$ shared with $S$ are denoted by $u^+$ and $u^-$, The normal displacement jump reads $[u_n] = [u] \cdot n$ while tangent displacement jumps reads $[u_t] = [u] \cdot t$ and $[u_s] = [u] \cdot s$. Following Nkoumbou Kaptchouang et al. (2021), the strain rate over the
cohesive surface is approximated through the following enriched scaling rule:

$$ \dot{\epsilon} \approx \frac{[u] \otimes n + n \otimes [u]}{h} + \dot{\gamma}(i - n \otimes n) = \begin{bmatrix} \frac{[\dot{u}_n]}{h} & \frac{[\dot{u}_t]}{2h} & \frac{[\dot{u}_s]}{2h} \\ \frac{[\dot{u}_t]}{2h} & \dot{\gamma} & 0 \\ \frac{[\dot{u}_s]}{2h} & 0 & \dot{\gamma} \end{bmatrix}, \quad (2) $$

where $\otimes$ denotes the dyadic product and $i$ is the second order identity tensor. The strain rate $\dot{\gamma}$ is an additional in-plane component which depends on the loading state provided by a bulk finite element adjacent to the cohesive zone as detailed in the sequel (see 2.3.2).

<figure><img src="image_2.png" /><figcaption>Figure 1: The damage process occuring within a narrow band of thickness $h$ is summarized along a zero thickness cohesive surface $S$ with a local normal-tangent basis $(n, t, s)$.</figcaption></figure>

## 2.2 Bulk viscoplastic damage model of Salvo et al. (2015)

The damage viscoplastic constitutive model of Salvo et al. (2015) for uranium dioxide is defined by the help of a stress potential $\Psi$:

$$ \dot{\epsilon} = \frac{\partial \Psi}{\partial \sigma} \quad \text{with} \quad \Psi(\sigma_m, \sigma_{eq}, f) = \sigma_0 K \exp\left(-\frac{Q}{rT}\right) \cosh\left(\frac{\Sigma}{\sigma_0}\right), \quad (3) $$

where $\dot{\epsilon}$ is the viscoplastic strain rate tensor, $\sigma_m = 1/3 \text{tr}(\sigma)$ is the hydrostatic part of the Cauchy stress tensor $\sigma$, $\sigma_{eq} = \sqrt{3/2} \sigma^d : \sigma^d$ is the Von Mises stress where symbol ':' denotes the double dot product, $\sigma^d = \sigma - \sigma_m i$ is the deviatoric part of the stress, $i$ is the second order identity tensor, $f$ is the porosity, $\sigma_0$, $K$ and $Q$ are constants, $r$ is the perfect gas constant, $T$ is the temperature and $\Sigma$ is an effective stress defined by

$$ \Sigma = \sqrt{\alpha(f)\sigma_m^2 + \beta(f)\sigma_{eq}^2}, \quad (4) $$

where $\alpha(f)$ and $\beta(f)$ are two scalar functions derived from the hollow sphere model of Michel et al. (1992):

$$ \alpha(f) = \frac{9}{4} \left( n \left( f^{-\frac{1}{n}} - 1 \right) \right)^{-\frac{2n}{n+1}} \quad \text{and} \quad \beta(f) = \left( 1 + \frac{2}{3} f \right) \left( 1 - f \right)^{-\frac{2n}{n+1}}, \quad (5) $$

where $n$ is an exponent related a viscous power law. The viscoplastic strain rate tensor $\dot{\epsilon}$ is thus obtained by the derivation of (3):

$$ \dot{\epsilon} = K \exp\left(-\frac{Q}{rT}\right) \sinh\left(\frac{\Sigma}{\sigma_0}\right) \left(\frac{1}{3}\alpha(f)\frac{\sigma_m}{\Sigma}i + \frac{3}{2}\beta(f)\frac{1}{\Sigma}\sigma^d\right), \quad (6) $$

where is readily to see that $\sigma_0$ has the dimension of a stress, $K$ of a strain rate and $Q$ of an activation energy. Two primary characteristics of this constitutive model are the adoption of a hyperbolic sine law instead of a combination of power laws, and the incorporation of pore compressibility. This enables the presence of a hydrostatic component of the strain rate, and
consequently, the evolution of porosity as a damage parameter, as illustrated in (18). A simple derivation of the dual of Ψ reported in A gives the strain rate potential Φ:

$$ \sigma = \frac{\partial \Phi}{\partial \dot{\epsilon}} \quad \text{with} \quad \Phi(\dot{\epsilon}_m, \dot{\epsilon}_{eq}, f) = \sigma_0 \dot{E} \left( \sinh^{-1} \left( \frac{\dot{E}}{K \exp\left(-\frac{Q}{rT}\right)} \right) - \sqrt{1 + \left( \frac{K \exp\left(-\frac{Q}{rT}\right)}{\dot{E}} \right)^2} \right), \quad (7) $$

where $\dot{\epsilon}_m = 1/3 \text{tr}(\dot{\epsilon})$ and $\dot{\epsilon}_{eq} = \sqrt{2/3\dot{\epsilon}^d : \dot{\epsilon}^d}$ are the hydrostatic and equivalent parts of $\dot{\epsilon}$ respectively, $\dot{\epsilon}^d$ is the deviatoric part of $\dot{\epsilon}$ and $\dot{E}$ is an effective strain rate defined by:

$$ \dot{E} = \sqrt{\frac{(3\dot{\epsilon}_m)^2}{\alpha(f)} + \frac{(\dot{\epsilon}_{eq})^2}{\beta(f)}}. \quad (8) $$

The stress tensor is thus derived from (7):

$$ \sigma = \frac{\partial \Phi}{\partial \dot{\epsilon}} = \sigma_0 \sinh^{-1} \left( \frac{\dot{E}}{K \exp\left(-\frac{Q}{rT}\right)} \right) \frac{1}{\dot{E}} \left( \frac{3}{\alpha(f)} \dot{\epsilon}_m \dot{i} + \frac{2}{3\beta(f)} \dot{\epsilon}^d \right). \quad (9) $$

## 2.3 Cohesive zone model

### 2.3.1 Traction-separation law

Combining the surface reduced kinematics (2) with the constitutive equation (9) yields to the following expression for the stress vector $R([u]) = \sigma \cdot n$ over the cohesive surface S:

$$ R = R_n n + R_t t + R_s s \quad \text{with} \quad \left\{ \begin{array}{l} R_n = C \left( \frac{1}{\alpha(f)} \left( \frac{[u_n]}{h} + 2\gamma \right) + \frac{4}{9\beta(f)} \left( \frac{[u_n]}{h} - \gamma \right) \right), \\ R_t = C \left( \frac{1}{3\beta(f)} \frac{[u_t]}{h} \right), \\ R_s = C \left( \frac{1}{3\beta(f)} \frac{[u_s]}{h} \right), \end{array} \right. \quad (10) $$

where C is a prefactor

$$ C = \sigma_0 \sinh^{-1} \left( \frac{E}{tK \exp\left(-\frac{Q}{rT}\right)} \right) \frac{1}{E} \quad (11) $$

and $R_n = R \cdot n$ is the normal part of the cohesive stress vector, $R_t = R \cdot t$ and $R_s = R \cdot s$ are the tangential parts, E is the effective strain as a counterpart of (8)

$$ E = \sqrt{\frac{(3\epsilon_m)^2}{\alpha(f)} + \frac{(\epsilon_{eq})^2}{\beta(f)}}. \quad (12) $$

The traction-separation law (10) was obtained by imposing a monotonous and radial strain, such that:

$$ \epsilon(t) = \dot{\epsilon} t, \quad (13) $$

where t denotes the time and $\dot{\epsilon}$ is a constant strain rate.

### 2.3.2 Estimate of the in-plane strain

The in-plane strain $\gamma = \dot{\gamma}t$ is obtained by importing into the cohesive zone informations from the bulk close to the cohesive zone. Among several possibilities, it is chosen to compute $\gamma$ by keeping the equality between the stress triaxiality $\eta_\sigma = \sigma_m/\sigma_{eq}$ in the cohesive zone and the
stress triaxiality $\eta_{\sigma^{vol}} = \sigma_m^{vol}/\sigma_m^{vol}$ computed from the stress $\sigma^{vol}$ in a volumetric element adjacent to the cohesive zone:

$$ \eta_{\sigma} = \eta_{\sigma^{vol}}. \qquad (14) $$

From (6), the strain rate triaxiality in the cohesive zone, $\eta_{\dot{\epsilon}} = \dot{\epsilon}_m/\dot{\epsilon}_{eq}$, is expressed in terms of the stress triaxiality through (see A for details):

$$ \eta_{\dot{\epsilon}} = \frac{\alpha(f)}{3\beta(f)}\eta_{\sigma}. \qquad (15) $$

Using (14), (15) and the specific loading (13), the in-plane strain $\gamma$ is thus obtained by solving the following equation:

$$ \frac{\alpha(f)}{3\beta(f)}\eta_{\sigma^{vol}} = \frac{\frac{[u_n]}{h} + 2\gamma}{2\sqrt{\left(\frac{[u_n]}{h} - \gamma\right)^2 + \frac{3}{4}\frac{[u_t]^2}{h^2} + \frac{3}{4}\frac{[u_s]^2}{h^2}}}. \qquad (16) $$

Denoting $\eta = \eta_{\dot{\epsilon}}$ for sake of simplicity, a straightforward calculation gives two solutions for $\gamma$ depending $\eta$:

$$ \gamma_{\pm}(\eta) = \frac{2\eta^2 + 1}{2(\eta^2 - 1)} \frac{[u_n]}{h} \pm \frac{\sqrt{\eta^2 \left(9\frac{[u_n]^2}{h^2} + 3(1-\eta^2)\left(\frac{[u_t]^2}{h^2} + \frac{[u_s]^2}{h^2}\right)\right)}}{2(\eta^2 - 1)}. \qquad (17) $$

For a simple traction along the normal direction ($[u_n] > 0$, $[u_t] = [u_s] = 0$), the solution $\gamma_+(\eta)$ exhibits two discontinuities (for $\eta = \pm 1$), while $\gamma_-(\eta)$ remains a continuous function. The negative branch is thus selected in the following: $\gamma(\eta) = \gamma_-(\eta)$.

### 2.3.3 Evolution of the porosity

In Salvo et al. (2015), the evolution of porosity as a damage variable is given by the mass balance equation:

$$ \dot{f} = (1 - f) \text{tr}(\dot{\epsilon}). \qquad (18) $$

The time integration of this equation together with (2) gives:

$$ f = 1 - (1 - f_0) \exp\left(-\frac{[u_n]}{h} - 2\gamma\right), \qquad (19) $$

where $f_0$ is the initial porosity. Following GTN-type ductile damage models, e.g. Tvergaard (1990), an effective porosity $f^*$ is introduced in order to accelerate the porosity growth beyond a critical porosity $f_c$:

$$ f^* = \begin{cases} f & \text{if } f < f_c, \\ f_c + k(f - f_c) & \text{if } f \ge f_c, \end{cases} \qquad (20) $$

where $k$ is a large constant parameter. This effective porosity $f^*$ replaces the true porosity $f$ when computing the stress vector $R$. The stress triaxiality $\eta_{\sigma^{vol}}$ is transmitted to the cohesive zone until the initiation of coalescence (i.e. since $f \le f_c$), and remains constant at its final value during the post-coalescence phase (i.e. when $f > f_c$).

### 2.3.4 Unloading case

As stipulated in Nkoumbou Kaptchouang et al. (2021), local unloading may occur during crack propagation within a fully general structural computation. It is essential to precisely define the cohesive zone model's response during unloading, considering the irreversible nature of the damage process. First, a criterion must be established to distinguish between loading and unloading states. Following Nkoumbou Kaptchouang et al. (2021), unloading is addressed
separately along each direction. Unloading occurs along the normal and tangential directions respectively, if:

$$|[u_i]^{(t)}| < |[u_i]^{(t-\delta t)}| \quad \text{for } i \in \{n, t, s\}, \qquad (21)$$

where $t$ is the time and $\delta t$ is the time increment. Second, following again Nkoumbou Kaptchouang et al. (2021), the behavior of the cohesive zone during unloading is assumed to be elastic. As a result, the stress vector during unloading $\bar{R}$ is written as follows:

$$\bar{R}_i = (R_i^{int})_{\max} \frac{[u_i]}{[u_i]_{\max}}, \quad \text{for } i \in \{n, t, s\} \text{ and } [u_i]_{\max} = \max(|[u_i]^{(\tau)}|, 0 \le \tau \le t), \quad (22)$$

where $(R_i^{int})_{\max}$ is the component $i$ of the stress vector $R^{int}$ evaluated at the maximal opening $([u_n]_{\max}, [u_t]_{\max}, [u_s]_{\max})$.

These instances of unloading situations do no longer meet the requirement of a local monotonous and radial strain, assumption (13). Nevertheless, we will subsequently verify that this deviation does not result in significant unphysical behaviors in the practical scenarios under consideration.

### 2.3.5 Intrinsic formulation

The cohesive zone model is transformed into an intrinsic model with a finite initial slope in order to fulfill requirements of our cohesive-volumetric numerical implementation, see e.g. Perales et al. (2010), and also to tackle elastic cohesive reloading after a partial or total unloading. The stress vector $R$ is replaced by:

$$R_i^{int} = \min(C_i[u_i], R_i), \quad i \in \{n, t, s\}, \qquad (23)$$

where $C_n, C_t$, and $C_s$ are initial cohesive stiffnesses. The additional overall compliance induced by this intrinsic formulation is controlled using theoretical criteria proposed by Blal et al. (2013).

### 2.3.6 Energetic failure criteria

A key parameter of the resulting cohesive zone model is the critical porosity $f_c$ which characterizes the onset of void coalescence and thus the failure of the cohesive zone. We suggest relating this parameter to a fracture criterion based on energetic considerations. More precisely, two energetic criteria are considered independently along the normal and tangential directions:

$$G^I = \int_0^{[u_n]_{\max}} \bar{R}_n([\mathbf{u}]) d[u_n] > G_c^I \qquad (24)$$

and

$$G^II = \int_0^{[u_t]_{\max}} \bar{R}_t([\mathbf{u}]) d[u_t] + \int_0^{[u_s]_{\max}} \bar{R}_s([\mathbf{u}]) d[u_s] > G_c^II. \qquad (25)$$

These criteria introduce two critical values for the cohesive energies, denoted by $G_c^I$ and $G_c^II$, respectively for the normal and tangential modes of crack propagation. If one of these criteria is satisfied, then the critical porosity $f_c$ is set to the value of the current porosity $f$. Consequently, the effective porosity $f^*$ is highly accelerated following (20), leading to a sudden material failure. These critical cohesive energies can be regarded as critical values for energy release rate. They can be estimated for uranium dioxide by analyzing experimental data. Several studies have been conducted to assess the fracture properties of uranium dioxide. For non-irradiated UO$_2$, numerous results from macroscopic tests are available, see for example Evans et al. (1969) or Gatt et al. (2015) among others. However, for irradiated UO$_2$, obtaining specimens for macroscopic testing is challenging due to partial pellet fracturing. An alternative approach is to scale down the tests. Many results from micro and nano indentation tests can be found, for example, in Kutty et al. (1987), Matzke et al. (1980), and Spino et al. (1996). Nonetheless, these measurements are influenced by numerous parameters that are difficult to control (Quinn et al. 2007), such as grain size, porosity, and crack geometry. More recently, other methods have been developed to complement the results obtained from indentation tests on irradiated fuel (Hosemann et al. 2015;
Frazer 2018; Henry 2019; Henry et al. 2020). The concept involves replicating conventional macroscopic tests by downsizing them, while mitigating the microstructural defects present in the fuels (such as irradiation bubbles and microcracks). For instance, Henry et al. (2020) conducted bending tests on microbeams of irradiated fuel, positioned within the grain or at grain boundaries. Quantitatively, critical values for the energy release rate can be estimated from tests performed on macroscopic specimens or from indentation tests, ranging from 1 J m⁻² to 24 J m⁻² (up to about 40 J m⁻² in the high burn-up structure zone, from the results of Spino et al. (1996)). Values estimated from microbeam tests fall within the range of 7 to 24 J m⁻². Ultimately, these findings provide a ballpark figure for the energy release rate of uranium dioxide, with potential values ranging from a few J m⁻² to a few tens of J m⁻².

The criterion (24)–(25) has to be considered as a first attempt to differentiate between the normal and tangential modes of crack propagation. More complex criteria can be derived, see e.g. Shen et al. (1994). Our main idea underlying the material under consideration in the present paper is that it is significantly easier to initiate and propagate a crack in mode I compared to mode II.

### 2.3.7 Numerical implementation into a cohesive-volumetric finite element formulation

The cohesive zone model was then implemented in the XPER code (eXtended cohesive zone models and PERiodic homogenisation) dedicated to the simulation of fracture dynamics of heterogeneous materials (Perales et al. 2010). XPER is based on the Non-Smooth Contact Dynamics (NSCD) method (Jean 1999; Jean et al. 2001), wherein: 1) cohesive models are integrated as mixed boundary conditions between each volumetric finite element, 2) instabilities linked to non-convex surface damage behavior are circumvented through a dynamic approach, and 3) unilateral contact and friction on the crack lips are neither regularized nor penalized. In this approach, each element (or each set of elements) of the mesh can be considered as an independent body. A slight modification of Signorini-Coulomb conditions at the interface between bodies is used to introduce the cohesive zone model (Jean et al. 2001; Acary et al. 2006). A pivotal aspect of this approach is that the bulk behavior of each element is governed by a hardening constitutive law devoid of any damage, while the damage is solely encapsulated within the cohesive zone model at the interface between the bodies.

## 3 Application: simulation of a uniaxial compression test

In this section, the effectiveness of the proposed cohesive zone model approach in predicting crack initiation and propagation is demonstrated through an illustrative example that is of typical interest in modeling the mechanical behavior of uranium dioxide at high temperature and strain rate. Predicted results are compared with experimental results of Salvo et al. (2015). The compression of a uranium dioxide fuel pellet at 1100 °C and 0.1 s⁻¹ is considered in the present paper.

### 3.1 Experiment of Salvo et al. (2015)

The samples used in Salvo et al. (2015) are commercial uranium dioxide sintered cylindrical pellets of approximately $d_0 = 8.2$ mm in diameter and $H_0 = 12$ mm in height. The mean porosity of the samples is about 3.6% and the average grain size is close to 10 µm. Strain-rate driven tests are performed on the samples. Each sample is placed between two tungsten end-plates. The displacement rate of the crosshead is maintained constant until reaching final discharge, just prior to complete pellet failure. The crosshead speed is adjusted while monitoring the pellet’s continuous contraction. The true macroscopic strain $E_{zz}$ of the pellet is defined by:

$$E_{zz} = -\ln\left(\frac{H}{H_0}\right), \qquad (26)$$

where the deformed height of the pellet is denoted by $H$. The macroscopic axial stress $\Sigma_{zz}$ is defined from the measured force $F_z$:

$$\Sigma_{zz} = -\frac{4F_z}{\pi d_0^2}. \qquad (27)$$
The evolution of $\Sigma_{zz}$ with respect to $E_{zz}$ obtained at 1100 °C and 0.1 s⁻¹ is set as a reference curve in the sequel. In this test, the macroscopic axial strain reaches about 6%. Note that, due to friction at pellet ends, the sample exhibits a barrel shape with a noticeable increase in diameter at mid-height. The test was interrupted because of the onset of axial cracks (1 mm to 3 mm long). An example of a typical shape of such a cracked sample obtained by Salvo et al. (2015) is shown in Figure 2 (left). Note also that a similar compression test performed by Salvo et al. (2015) on a high density fuel pellet (8.2 mm in diameter, 15 mm in height, porosity close to 1.9%) at 1100 °C and 0.1 s⁻¹ leads to a complete failure of the sample. The resulting shape of the sample is shown in Figure 2 (right). As observed by Salvo et al. (2015), the cracks are mostly located at the center of the pellet and oriented along the loading direction (axial cracks). In this test, the remains of the pellet exhibit several axially oriented micro-columns, and the failure is attributed to shear cracking, characterized by oblique cracks with respect to the loading direction.

<figure><img src="image_3.png" /><figcaption>**Figure 2:** Typical shapes of cracked samples: commercial uranium dioxide pellet (left), high density uranium dioxide pellet (right), from Salvo (2014).</figcaption></figure>

## 3.2 Set up of the numerical simulation

The geometry, mesh, and boundary conditions are summarized in Figure 3. Linear tetrahedrons finite elements with size $l_{mesh} = 600$ µm are used. The pellet and the two end-plates are modeled as plain cylinders in contact. The displacement of the nodes on the bottom circular surface of the structure is set to zero along the z axis (i.e. the compression direction). A velocity along the z axis is imposed on the top circular surface of the structure (−1.6 × 10⁻³ m s⁻¹) in order to get a macroscopic uniaxial strain rate close to 0.1 s⁻¹, following the experimental procedure of Salvo et al. (2015). The displacements along the x and y directions are set to zero for the nodes located on the top and bottom external circles between the pellet and the end-plates. This specific boundary condition is introduced to approach the test conditions, in which the pellet is supported and centered by end-plates presenting a central gap with a diameter close to the diameter of the pellet. The end-plates behavior is assumed to be stiffer than the sample and governed by an isotropic elastic constitutive law (large Young modulus set to 2000 GPa, Poisson ratio set to 0.3 and density set to 10 500 kg m⁻³. The pellet behavior is governed by the isotropic elasto-viscoplastic constitutive law of Salvo et al. (2015) with no damage (no porosity evolution) keeping in mind that damage is entirely taken into account through the cohesive zones:

$$ \sigma = \mathbb{C} : (\varepsilon^{tot} - \varepsilon), \qquad (28) $$

where $\varepsilon^{tot}$ is the total strain, $\varepsilon$ is the viscoplastic strain whose evolution is given by (6), and $\mathbb{C}$ is the fourth order elasticity tensor. The Young modulus of the sample is set to 170 GPa according to Martin (1989) (for $f = 3.6\%$ and $T = 1100$ °C), the Poisson ratio is set to 0.3 and density to 10 500 kg m⁻³.

<table><caption>**Table 1:** Parameters for the bulk viscoplastic behavior and the cohesive zone model.</caption><thead><tr><th>K</th><th>σ<sub>0</sub></th><th>n</th><th>f<sub>0</sub></th><th>Q</th><th>r</th></tr></thead><tbody><tr><td>11 652 × 10<sup>4</sup> s<sup>-1</sup></td><td>5 MPa</td><td>6</td><td>3.6%</td><td>532.7 kJ mol<sup>-1</sup></td><td>8.314 J mol<sup>-1</sup> K<sup>-1</sup></td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>Figure 3: Geometry, mesh, and boundary conditions used for the simulation of the uniaxial compression test.</figcaption></figure>

The cohesive zone model of section 2 is adopted between each solid element of the sample and the thickness of the cohesive zone h is chosen to be equal to the mesh size

$$h = l_{\text{mesh}}.$$

All the parameters are reported in Table 1. The initial cohesive stiffness parameters $C_n$, $C_t$ and $C_s$ are set to $10^{17}$ Pa m$^{-1}$, ensuring the criterion of Blal et al. (2013). Frictionless conditions are introduced at the interfaces between the pellet and the plates. The resolution is performed in dynamics and the loading velocity is multiplied by 4000 in order to speed up the results. By doing so, the velocity along the z axis which is imposed on the top circular surface of the structure is set to $\dot{u}_z = 4000 \times (-1.6 \times 10^{-3} \text{ m s}^{-1})$. This increase in the loading velocity is balanced by an increase in the K parameter (i.e. the original K value specified in Salvo et al. (2015) is multiplied by 4000). It justifies the high value of K in the Table 1 with respect to the original K value specified in Salvo et al. (2015). Note that all the other parameters reported in Table 1 are coming from Salvo et al. (2015). The macroscopic axial stress and strain are computed from the simulation using (27) and:

$$E_{zz} \approx -\ln \left( 1 + \frac{u_z}{H_0} \right). \qquad (29)$$

A high value for $k = 6000$ is used in (20) in order to induce a strong increase in the effective porosity as soon as one of the fracture criteria is reached, see (24) and (25). By doing so, the cohesive stress is quickly set to zero, corresponding to a sudden local failure. To the end, all the simulation parameters are set except for the two critical values for the cohesive energies ($G_c^I$ and $G_c^II$).

### 3.3 Results and discussion

The axial cracks experimentally observed prior to the total pellet rupture suggest, as expected, that the critical cohesive energy $G_c^I$ for the normal mode of crack propagation has to be lower than the critical cohesive energy $G_c^II$ for the tangential mode. A large value for $G_c^II = 400 \text{ J m}^{-2}$ is thus imposed. To fix ideas, it is about 2 orders of magnitude higher than the smallest value measured by Evans et al. (1969) from 1000 °C to 1200 °C from three-point bending tests conducted on uranium dioxide specimens. The $G_c^I$ value ranges from $2 \text{ J m}^{-2}$ to $400 \text{ J m}^{-2}$. The simulated macroscopic stress-strain curves are shown in Figure 4. In this figure, the experimental curve of Salvo et al. (2015) has been shifted by -0.7% along the strain axis, as the initial slope reflects the loading of the compression test machine (gap closure between the plates and the specimen). The following observations can be made.
* The present numerical results have to be compared to the numerical simulation of Salvo et al. (2015) whose the present cohesive zone model is derived, they do not have to be compared to the experimental results of Salvo et al. (2015)

* The high oscillations observed at the beginning of the simulated stress-strain curves are due to the artificial increase of the loading velocity together with the dynamic resolution. Once these oscillations have been smoothed, the values of the stress remain close to that simulated by Salvo et al. (2015) in 2D using standard finite elements in quasi-static.

* All the simulated compression curves superimposed up to 5% of strain.

* As expected, the decrease of the stress due to damage and crack propagation is delayed when $G_c^I$ increases from $2 \text{ J m}^{-2}$ to $6 \text{ J m}^{-2}$, see Figure 5. The rupture is delayed with increased critical energy.

* More surprisingly, the rupture occurs slightly earlier when $G_c^I$ goes from $6 \text{ J m}^{-2}$ to $8 \text{ J m}^{-2}$, and then does not vary anymore between $8 \text{ J m}^{-2}$ to $400 \text{ J m}^{-2}$, see Figure 5.

<figure><img src="image_5.png" /><figcaption>**Figure 4:** Compression of a uranium dioxide fuel pellet at 1100 °C and 0.1 s⁻¹. Macroscopic stress-strain curves: simulations for different values of the critical cohesive energy $G_c^I$ (gray), experiment (dashed black) and simulation (black) of Salvo et al. (2015).</figcaption></figure>

From this numerical study, the value of $G_c^I = 6 \text{ J m}^{-2}$ is selected. Firstly, for this value, sample failure occurs at a strain exceeding 6%, consistent with the experimental findings of Salvo et al. (2015) (see Figure 4). Secondly, this value closely aligns with the critical energy release rate at 1100 °C estimated by Evans et al. (1969) from three-point bend tests conducted on notched bars of uranium dioxide (specimens with a higher grain size of approximately 25 µm compared to Salvo et al. (2015), but with a smaller porosity of approximately 3%). Thirdly, Figure 6 and Figure 7 show that, with this value, the outer boundaries of the pellet initially exhibit axial cracks, aligned with the loading direction, as observed by Salvo et al. (2015). Then, a total failure of the pellet occurs with multiple cracks, including oblique cracks with respect to the loading direction, consistent with the experimental observations of Salvo et al. (2015). Note that, in these figures, the damage variable $\beta = 1 - f^*$ is introduced to exhibit the cracks, $\beta = 0$ corresponds to a fully damaged cohesive zone.

Figure 6 shows that there is a change in the cracks observed on the pellet curved surface at the beginning of its failure, while varying the $G_c^I$ value. For $G_c^I = 400 \text{ J m}^{-2}$, one can mainly observe an oblique crack with respect to the direction of the loading. For $G_c^I = 20 \text{ J m}^{-2}$, one can mainly observe an oblique crack together with some axial cracks. For $G_c^I = 6 \text{ J m}^{-2}$ or $G_c^I = 4 \text{ J m}^{-2}$, one can mainly observe axial cracks (shorter cracks for $G_c^I = 6 \text{ J m}^{-2}$ than for $G_c^I = 4 \text{ J m}^{-2}$). It suggests that, during this test, there is a competition between normal and tangential crack propagation modes.

In Figure 6 and Figure 7, one can observe that the areas between the pellet and the two end plates appear in blue, i.e. $\beta = 0$. It does not represent any damage of the cohesive zone model. This is simply due to the sliding contact conditions imposed at the interfaces between the pellet and the plates. Nevertheless, one can also observe in Figure 6 that, just below or above this interface
<figure><img src="image_6.png" /><figcaption>**Figure 5:** Compression of a uranium dioxide fuel pellet at 1100 °C and 0.1 s⁻¹. Macroscopic stress-strain curves: simulations for different values of the critical cohesive energy $G_c^I$, zoom of Figure 4. Simulation Salvo et al. (2015) (black), present results (gray), highlighted values (blue).</figcaption></figure>

between the pellet and the end plates, a zone of the pellet curved surface exhibits a complete failure. It was observed that this local failure on the pellet curved surface appears just before the development of the axial or oblique cracks on the the pellet curved surface.

From these simulations, the selected values of the critical cohesive energies are $G_c^I = 6 \text{ J m}^{-2}$ and $G_c^{II} = 400 \text{ J m}^{-2}$. Note that such a large difference between $G_c^I$ and $G_c^{II}$ was also obtained by Shen et al. (1994) for gypsum ($G_c^I = 30 \text{ J m}^{-2}$ and $G_c^{II} = 550 \text{ J m}^{-2}$ for gypsum).

Using this set of parameters, the cohesive zone model is seen to be suitable for uranium dioxide fuel at high temperature (1100 °C) and strain rate (0.1 s⁻¹).

<figure><img src="image_7.png" /><figcaption>**Figure 6:** Compression of a uranium dioxide fuel pellet at 1100°C and 0.1/s. Beginning of the failure of the pellet curved surface ($E_{zz} \approx 5.79\%$) for different $G_c^I$ values: $\beta$ field (initial configuration).</figcaption></figure>

# 4 Conclusion

This paper introduces a novel cohesive zone model, derived from the volume-to-surface methodology outlined in Nkoumbou Kaptchouang et al. (2021) and applied to the viscoplastic damage model proposed by Salvo et al. (2015). Notably, energetic failure criteria are incorporated to
<figure><img src="image_8.png" /><figcaption>**Figure 7:** Compression of a uranium dioxide fuel pellet at 1100°C and 0.1/s. Total failure of the pellet at the end of the simulation ($E_{zz} \approx 6.69\%$) with $G_c^I = 6\text{J/m}^2$. Left: $\beta$ field (initial configuration, same legend as in Figure 6). Right: the center of each blue line indicates the failure of a contact point between two elements in the pellet (current configuration).</figcaption></figure>

differentiate between the normal and tangential modes of crack propagation. The resulting model is dedicated specifically for uranium dioxide fuel and explicitly accounts for temperature, strain rate, and porosity. The effectiveness of this cohesive zone model approach is demonstrated through 3-dimensional simulations of a uniaxial compression test conducted at high temperature and strain rate. Comparison with experimental results from Salvo et al. (2015) validates the model's predictive capabilities for crack initiation and propagation. Finally, the identified critical values for the cohesive energies for the normal and tangential modes of crack propagation are set to $6 \text{ J m}^{-2}$ and $400 \text{ J m}^{-2}$ respectively.

The general modeling proposed here is expected to be useful for a wide range of temperatures and strain rates. To show this, other compression tests in terms of temperatures and strain rates presented in Salvo et al. (2015) could also be simulated. This is left for future work.

## References

Acary, V. and Y. Monerie (2006). *Nonsmooth fracture dynamics using a cohesive zone approach*. Tech. rep. INRIA

Anvari, M., J. Liu, and C. Thaulow (2007). “Dynamic ductile fracture in aluminum round bars: experiments and simulations”. *International Journal of Fracture* 143.4, pp. 317–332

Anvari, M., I. Scheider, and C. Thaulow (2006). “Simulation of dynamic ductile crack growth using strain-rate and triaxiality-dependent cohesive elements”. *Engineering Fracture Mechanics* 73.15, pp. 2210–2228

Barenblatt, G. I. (1959). “The formation of equilibrium cracks during brittle fracture. General ideas and hypotheses. Axially-symmetric cracks”. *Journal of Applied Mathematics and Mechanics* 23.3, pp. 622–636

Blal, N., L. Daridon, Y. Monerie, and S. Pagano (2013). “Micromechanical-based criteria for the calibration of cohesive zone parameters”. *Journal of Computational and Applied Mathematics* 246, pp. 206–214

Dugdale, D. (1960). “Yielding of steel sheets containing slits”. *Journal of the Mechanics and Physics of Solids* 8.2, pp. 100–104

Evans, A. and R. Davidge (1969). “The strength and fracture of stoichiometric polycrystalline $UO_2$”. *Journal of Nuclear Materials* 33.3, pp. 249–260

Frazer, D. (2018). “Elevated Temperature Small Scale Mechanical Testing of Uranium Dioxide”. PhD thesis. Berkeley University

Fuketa, T. (2012). “Transient Response of LWR Fuels (RIA)”. *Comprehensive Nuclear Materials*. Ed. by R. J. Konings. Oxford: Elsevier, pp. 579–593

Gatt, J.-M., J. Sercombe, I. Aubrun, and J.-C. Ménard (2015). “Experimental and numerical study of fracture mechanisms in $UO_2$ nuclear fuel”. *Engineering Failure Analysis* 47, pp. 299–311
Gurson, A. (1977). "Continuum Theory of Ductile Rupture by Void Nucleation and Growth: Part I - Yield Criteria and Flow Rules for Porous Ductile Media". *Journal of Engineering Materials and Technology* 99, pp. 1-15

Henry, R. (2019). "Caractérisation locale des propriétés à la rupture du combustible nucléaire irradié". PhD thesis. INSA Lyon

Henry, R., I. Zacharie-Aubrun, T. Blay, N. Tarisien, S. Chalal, X. Iltis, J.-M. Gatt, C. Langlois, and S. Meille (2020). "Irradiation effects on the fracture properties of UO₂ fuels studied by micro-mechanical testing". *Journal of Nuclear Materials* 536, p. 152179

Hosemann, P., C. Shin, and D. Kiener (2015). "Small scale mechanical testing of irradiated materials". *Journal of Materials Research* 30, pp. 1231-1245

Jean, M. (1999). "The non-smooth contact dynamics method". *Computer Methods in Applied Mechanics and Engineering* 177 (3), pp. 235-257

Jean, M., V. Acary, and Y. Monerie (2001). "Non-smooth contact dynamics approach of cohesive materials". *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences* 359 (1789), pp. 2497-2518

Koo, Y.-H., D.-S. Sohn, and B. Volkov (1997). "A comparative analysis of UO₂ and MOX fuel behavior under reactivity initiated accident". *Annals of Nuclear Energy* 24.11, pp. 859-870

Koplik, J. and A. Needleman (1988). "Void growth and coalescence in porous plastic solids". *International Journal of Solids and Structures* 24, pp. 835-853

Kutty, T. R. G., K. N. Chandrasekharan, J. P. Panakkal, and J. K. Ghosh (1987). "Fracture toughness and fracture surface energy of sintered uranium dioxide fuel pellets". *Journal of Materials Science Letters* 6, pp. 260-262

Labat, J., R. Largenton, and J.-C. Michel (2023). "Multiscale identification and NTFA reduction of the elasto-viscoplastic polycrystalline behaviour of uranium dioxide (UO₂) for wide range of loading conditions". *Journal of Nuclear Materials* 582, p. 154471

Martin, D. (1989). "The elastic constants of polycrystalline UO₂ and (U, Pu) mixed oxides: a review and recommandations". *High Temperatures - High Pressures* 21, pp. 13-24

Matzke, H., T. Inoue, and R. Warren (1980). "The surface energy of UO₂ as determined by Hertzian indentation". *Journal of Nuclear Materials* 91, pp. 205-220

Michel, J.-C., P. Suquet, and F. Thébaud (1994). "Une modélisation du rôle des interfaces dans le comportement des composites à matrice métallique". *Revue Européenne des Éléments Finis* 3.4, pp. 573-595

Michel, J. and P. Suquet (1992). "The constitutive law of nonlinear viscous and porous materials". *Journal of the Mechanics and Physics of Solids* 40.4, pp. 783-812. DOI: https://doi.org/10.1016/0022-5096(92)90004-L

Nkoumbou Kaptchouang, N. B., Y. Monerie, F. Perales, and P.-G. Vincent (2021). "Cohesive GTN model for ductile fracture simulation". *Engineering Fracture Mechanics* 242, p. 107437

Perales, F., F. Dubois, Y. Monerie, B. Piar, and L. Stainier (2010). "A NonSmooth Contact Dynamics-based multi-domain solver". *European Journal of Computational Mechanics* 19 (4), pp. 389-417

Quinn, G. D. and R. C. Bradt (2007). "On the Vickers Indentation Fracture Toughness Test". *Journal of the American Ceramic Society* 90, pp. 673-680

Salvo, M. (2014). "Etude expérimentale et modélisation du comportement mécanique du combustible UO₂ en compression à haute température et forte vitesse de sollicitation". PhD thesis. Université d'Aix Marseille

Salvo, M., J. Sercombe, T. Helfer, P. Sornay, and T. Désoyer (2015a). "Experimental characterization and modeling of UO₂ grain boundary cracking at high temperatures and high strain rates". *Journal of Nuclear Materials* 460, pp. 184-199

Salvo, M., J. Sercombe, J.-C. Ménard, J. Julien, T. Helfer, and T. Désoyer (2015b). "Experimental characterization and modelling of UO₂ behavior at high temperatures and high strain rates". *Journal of Nuclear Materials* 456, pp. 54-67

Scheider, I. (2009). "Derivation of separation laws for cohesive models in the course of ductile fracture". *Engineering Fracture Mechanics*. MatModels 2007 76.10, pp. 1450-1459

Schmitz, F. and J. Papin (1999). "High burnup effects on fuel behaviour under accident conditions: the tests CABRI REP-Na". *Journal of Nuclear Materials* 270.1, pp. 55-64
Schwalbe, K.-H., I. Scheider, and A. Cornec (2012). Guidelines for Applying Cohesive Models to the Damage Behaviour of Engineering Materials and Structures. Springer Science & Business Media. Springer

Shen, B. and O. Stephansson (1994). "Modification of the G-criterion for crack propagation subjected to compression". Engineering Fracture Mechanics 47.2, pp. 177-189

Siegmund, T. and W. Brocks (1999). "Prediction of the Work of Separation and Implications to Modeling". International Journal of Fracture 99.1-2, pp. 97-116

Siegmund, T. and W. Brocks (2000). "A numerical study on the correlation between the work of separation and the dissipation rate in ductile fracture". Engineering Fracture Mechanics 67.2, pp. 139-154

Spino, J., K. Vennix, and M. Coquerelle (1996). "Detailed characterisation of the RIM microstructure in PWR fuels in the burn-up range 40 - 67 GWd/tM". Journal of Nuclear Materials 231.3, pp. 179-190

Suquet, P. M. (1988). "Discontinuities and Plasticity". Nonsmooth Mechanics and Applications. Ed. by J. J. Moreau and P. D. Panagiotopoulos. Vol. 27. Vienna: Springer Vienna, pp. 279-340

Tvergaard, V. (1990). "Material Failure by Void Growth to Coalescence". Advances in Applied Mechanics. Ed. by J. Hutchinson and T. Wu. Vol. 27. New-York: Academic Press, pp. 83-151

Tvergaard, V. and J. Hutchinson (1992). "The relation between crack growth resistance and fracture process parameters in elastic-plastic solids". Journal of the Mechanics and Physics of Solids 40.6, pp. 1377-1397

**Acknowledgements** The authors would like to thank Maxime Salvo from the IRSN for fruitful discussions on the experimental characterization of uranium dioxide in compression. The authors would like to thank Frederic Perales from the IRSN, and Frederic Dubois, Remy Mozul from the CNRS, for their help in numerical computations. The authors would like to thank Maxence Wangermez from the CEA for his assistance in the implementation of the constitutive model. The simulations presented here were made possible thanks to the resources of the CCRT (Centre de Calcul Recherche et Technologie) of the CEA/DIF.

**Data availability** The raw and processed data required to reproduce these findings are available from the corresponding author upon reasonable request.

## A Derivation of the strain rate potential Φ

Φ and Ψ are dual potentials such that:

$$ \Phi(\dot{\varepsilon}) = \sigma : \dot{\varepsilon} - \Psi(\sigma). \qquad (30) $$

To simplify the expressions, a new variable, denoted by $\tilde{K}$, is introduced:

$$ \tilde{K} = K \exp \left( -\frac{Q}{rT} \right). \qquad (31) $$

From (3) and (6), it leads to:

$$ \Phi = \tilde{K} \sinh \left( \frac{\Sigma}{\sigma_0} \right) \left( \alpha \frac{\sigma_m^2}{\Sigma} + \beta \frac{\sigma_{eq}^2}{\Sigma} \right) - \sigma_0 \tilde{K} \cosh \left( \frac{\Sigma}{\sigma_0} \right) \qquad (32) $$

and

$$ \dot{\varepsilon}_m = \tilde{K} \sinh \left( \frac{\Sigma}{\sigma_0} \right) \frac{\alpha}{3} \frac{\sigma_m}{\Sigma}, \quad \dot{\varepsilon}^d = \tilde{K} \sinh \left( \frac{\Sigma}{\sigma_0} \right) \frac{3\beta}{2} \frac{1}{\Sigma} \sigma^d, \quad \dot{\varepsilon}_{eq} = \tilde{K} \sinh \left( \frac{\Sigma}{\sigma_0} \right) \beta \frac{\sigma_{eq}}{\Sigma}. \qquad (33) $$

It appears that the strain rate triaxiality can be expressed in terms of the stress triaxiality through:

$$ \frac{\dot{\varepsilon}_m}{\dot{\varepsilon}_{eq}} = \frac{1}{3} \frac{\alpha}{\beta} \frac{\sigma_m}{\sigma_{eq}}. \qquad (34) $$
To simplify the expressions, a new variable, denoted by $\tilde{S}$, is introduced:

$$ \tilde{S} = \sinh\left(\frac{\Sigma}{\sigma_0}\right). \qquad (35) $$

Equations (33) write:

$$ \sigma_m = \frac{3\dot{\epsilon}_m\Sigma}{\tilde{K}\tilde{S}\tilde{\alpha}}, \quad \sigma^d = \frac{2\Sigma}{3\tilde{K}\tilde{S}\tilde{\beta}}\dot{\epsilon}^d, \quad \sigma_{eq} = \frac{\dot{\epsilon}_{eq}\Sigma}{\tilde{K}\tilde{S}\tilde{\beta}}. \qquad (36) $$

Using the definition of $\Sigma$ (4), it leads to:

$$ \Sigma^2 = \alpha\sigma_m^2 + \beta\sigma_{eq}^2 = \frac{(3\dot{\epsilon}_m)^2\Sigma^2}{\tilde{K}^2\tilde{S}^2\alpha} + \frac{\dot{\epsilon}_{eq}^2\Sigma^2}{\tilde{K}^2\tilde{S}^2\beta} \quad \text{and} \quad \tilde{S} = \sqrt{\frac{1}{\alpha}\left(\frac{3\dot{\epsilon}_m}{\tilde{K}}\right)^2 + \frac{1}{\beta}\left(\frac{\dot{\epsilon}_{eq}}{\tilde{K}}\right)^2}. \qquad (37) $$

Using the definition of $\dot{E}$ (8), it writes:

$$ \tilde{S} = \frac{1}{\tilde{K}}\dot{E}. \qquad (38) $$

Finally, the strain rate potential writes:

$$ \Phi = \frac{9\dot{\epsilon}_m\Sigma}{\tilde{K}\tilde{S}\tilde{\alpha}}\dot{\epsilon}_m + \frac{2\Sigma}{3\tilde{K}\tilde{S}\tilde{\beta}}\dot{\epsilon}^d : \dot{\epsilon}^d - \sigma_0\tilde{K}\cosh\left(\sinh^{-1}\left(\frac{1}{\tilde{K}}\dot{E}\right)\right) \qquad (39) $$

or in simpler form:

$$ \Phi = \frac{\Sigma}{\tilde{K}\tilde{S}} \left( \frac{9\dot{\epsilon}_m^2}{\alpha} + \frac{\dot{\epsilon}_{eq}^2}{\beta} \right) - \sigma_0 \sqrt{\tilde{K}^2 + \dot{E}^2} = \frac{\Sigma}{\tilde{K}\tilde{S}} \dot{E}^2 - \sigma_0 \sqrt{\tilde{K}^2 + \dot{E}^2}. \qquad (40) $$

From the definition (35) and the expression (38), one gets:

$$ \Sigma = \sigma_0 \sinh^{-1}\left(\frac{\dot{E}}{\tilde{K}}\right). \qquad (41) $$

From (38), (40), and (41), the strain rate potential writes:

$$ \Phi = \sigma_0 \left( \sinh^{-1}\left(\frac{\dot{E}}{\tilde{K}}\right) \dot{E} - \sqrt{\tilde{K}^2 + \dot{E}^2} \right). \qquad (42) $$

Using the definition of $\tilde{K}$ (31), it leads to the expression (7).