# A variational phase-field model for ductile fracture depending on hydrostatic stresses

Anne-Sophie Sur, Laura De Lorenzis, Corrado Maurini, Odd Sture Hopperstad

► To cite this version:

Anne-Sophie Sur, Laura De Lorenzis, Corrado Maurini, Odd Sture Hopperstad. A variational phase-field model for ductile fracture depending on hydrostatic stresses. 2025. hal-05105037

HAL Id: hal-05105037

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

Preprint submitted on 10 Jun 2025

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
# A variational phase-field model for ductile fracture depending on hydrostatic stresses

Anne-Sophie Sur*,¹,⁴; Laura De Lorenzis², Corrado Maurini³, Odd Sture Hopperstad*,¹

¹ Department of Structural Engineering, Norwegian University of Science and Technology (NTNU), Richard Birkelands vei 1A, 7034 Trondheim, Norway, odd.hopperstad@ntnu.no, anne.s.sur@ntnu.no

² Department of Mechanical and Process Engineering, ETH Zürich, Tannenstr. 3, 8092 Zürich, Switzerland

³ Institut Jean Le Rond d'Alembert, Sorbonne University and CNRS, UMR 7190, 75005 Paris, France

⁴ Institute of Numerical Mathematics, Technische Universität Dresden (TUD), Zellescher Weg 21-25a, 01069 Dresden, Germany

## Abstract

We model ductile fracture for geometrically linear deformations by coupling plasticity and phase-field fracture models in a variationally consistent framework. The main aim of the proposed model is to account for the effect of stress triaxiality, in order to accurately reproduce ductile fracture, in particular, the instant and location of fracture initiation. For this purpose, we couple the modified Cam-Clay plasticity model with a phase-field fracture model. We study the behaviour of the model analytically in terms of homogeneous material responses, and numerically on plane-strain and axisymmetric specimens under tension with different notches.

*keywords:* phase-field, ductile fracture, stress triaxiality, variational approach

## 1. Introduction

Ductile fracture is the failure of materials associated with plastic deformation, and its correct understanding and modelling are important for the safety of engineering structures involving metals, polymers and composites. Complex, interacting, material-specific processes influencing the fracture behaviour, such as nucleation, growth and coalescence of voids, the stress state, and plastic strain, need to be considered to model ductile fracture. Various local damage models have been proposed and extensively used to model ductile fracture in different elasto-plastic materials, including some models which account for stress triaxiality effects [1–8]. However, the localisation of damage results in the loss of ellipticity of the incremental equilibrium equations, which leads to undesirable mesh dependency in finite element (FE) simulations [9–12].

For brittle fracture, phase-field modelling has emerged as an elegant and robust approach to regularise local damage models and retrieve an energetic equivalence with sharp-interface fracture models. The introduction of variational phase-field models to describe fracture goes back to Bourdin et al. [13]; it is based on the variational approach to fracture [14], which formulates brittle fracture mechanics as an energy minimisation problem. In phase-field models a regularisation is introduced by the penalisation of the gradient of the internal variable in the energy functional, involving an internal length. For vanishing internal length, the global minimisers of the regularised models converge to the sharp-interface Griffith models in the sense of Γ-convergence.

Associative plasticity can also be formulated as an energy minimisation problem [15–17]. Thus, a few models have attempted to couple plasticity and phase-field modelling of fracture within the variational framework. Alessi et al. [18] extend the phase-field approach to ductile fracture by introducing an
energy functional describing elasto-plastic material behaviour coupled with a phase-field damage model in a one-dimensional context. Further works by the same authors [19, 20] formulate the model in the multi-dimensional setting for von Mises perfect plasticity, but limit the study of its response to a uniaxial tension test. For the anti-plane case, Dal Maso et al. [21] prove the $\Gamma$-convergence of these models to cohesive fracture models. Kuhn et al. [22] propose another variational model using again the von Mises yield criterion with a linear hardening term.

Several other phase-field models for ductile fracture have been proposed, which depart in one way or another from the variational framework. Ambati et al. [23] propose a model where the elastic energy is split into a compressive and a tensile part, and the phase-field variable is coupled to plasticity through the elastic degradation function. In [24], the model in [23] is extended to finite strains. Ulloa et al. [25] study the influence of a hardening parameter and the option to consider gradient plasticity with a separate length scale. Miehe et al. [26] consider finite strains in the logarithmic strain space and adapt a history variable to the material properties, including a plasticity threshold on principal tensile stresses. The work by Borden et al. [27] for finite strains comprises an elastic energy splitting, considers von Mises plasticity with a threshold value and shows how stress triaxiality can be introduced through the history variable. In [28], the effect of stress triaxiality is accounted for by incorporating the Gurson-Tvergaard-Needleman (GTN) plasticity model [3]. The phase-field variable evolves based on a history variable that depends on the material porosity and on a given threshold value. In [29], Dittmann et al. propose a phase-field formulation where the effective stress under plastic deformation is determined by the GTN model. In [30], Choo and Sun present a pressure-sensitive fracture model for geomaterials such as porous rock. Drucker-Prager type plasticity is smoothly connected to an elliptical yield surface at high pressure and coupled to the phase-field model through a history variable. Many further phase-field models for ductile fracture have been proposed, including, for example, dynamic loading, the addition of gradient plasticity, various underlying plasticity models, and the degradation of the fracture toughness. Here, we do not aim at covering a full literature review on phase-field modelling for ductile fracture. The review paper [31] by Alessi et al. summarises the underlying theory of phase-field models coupled with plasticity and compares the constitutive equations of existing phase-field models and their numerical results for a uniaxial tension case. The more recent review paper by Marengo and Perego [32] includes theoretical aspects such as the derivation of the strong form equations from the energy functional and from a micromorphic approach.

Summarising, the vast majority of the available models are of non-variational nature. Since their equilibrium and damage evolution equations do not stem from an energy functional, these models enjoy a greater flexibility to reproduce more complex material behaviour, but they give up the theoretical and practical advantages of the variational setting [33]. On the other hand, the few variationally consistent models [18–20, 22], based on von Mises plasticity, are not able to reproduce important phenomenological features of ductile fracture in metals, e.g. triaxiality effects. Even in a specimen loaded in uniaxial tension, when necking starts to take place during plastic yielding a three-dimensional stress state occurs, in which hydrostatic stresses are superposed to axial stresses [6]. The more pronounced the neck, the larger the hydrostatic stresses, which increases the resulting stress level in the stress-strain response [6, 34]. This highlights the importance of considering triaxiality effects for the accurate prediction of ductile fracture when necking occurs. These effects are even more pronounced for complex geometries with notches, where additional hydrostatic stresses already occur in the elastic stage of the loading process.

In this paper, we propose a variational phase-field model coupled with plasticity that can capture with some flexibility stress triaxiality effects relevant to ductile fracture initiation. Our plasticity formulation is based on the modified Cam-Clay plasticity model, that dates back to Roscoe and Burland [35] and is widely used in geomechanics, where it is applied to describe the compressive behaviour of fine-grained soils such as sand, silt and clay [36–40], but may also be used to describe failure of other materials, e.g. snow [41]. As we show in this paper, this plasticity model may also be relevant for other materials such as metals, where the stress triaxiality plays a role in ductile fracture behaviour due to void evolution during plastic deformation. Unlike in standard von Mises plasticity, the yield function of the modified Cam-Clay plasticity model depends on the deviatoric and the hydrostatic stresses, whereby the influence of the hydrostatic stresses can be adjusted with a parameter. Hence, we describe stress
triaxiality dependent ductile fracture by coupling modified Cam-Clay plasticity with a phase-field model for fracture. The proposed variational model is restricted to geometrically linear deformations. Varying the intensity of the influence of stress triaxiality is limited to adjusting a single parameter mentioned above, which is included in the formulation through the yield surface.

The paper is structured as follows. In section 2 we first recall the modified Cam-Clay plasticity model, and then formulate the proposed phase-field model coupled with plasticity and its numerical implementation in the commercial FE program Abaqus/Standard (2022). Analytical homogeneous material responses are presented in section 3. Section 4 illustrates and discusses the numerical results for two-dimensional plane-strain geometries and axisymmetric tensile specimens. Conclusions are drawn in section 5.

In terms of notation, we denote by $\dot{\mathbf{\epsilon}}(\cdot) = d(\cdot)/dt$ the derivative of a quantity with respect to time $t$ and by $\partial_x(\cdot) = \partial(\cdot)/\partial x$ the partial derivative with respect to any other variable $x$. Also, we denote with $\cdot$ the inner product between same-order tensors, as $\|\mathbf{\epsilon}\|^2 = \mathbf{\epsilon} \cdot \mathbf{\epsilon}$ the (squared) Frobenius norm of the second-order tensor $\mathbf{\epsilon}$, and as $I_3 \in \mathbb{R}^{3,3}$ the identity matrix. For better clarity, we use round brackets for functional dependencies and square brackets for mathematical operations.

## 2. Formulation of the phase-field model coupled with plasticity

Let us consider a body occupying a domain $\Omega \subseteq \mathbb{R}^3$ and denote by $\mathbf{u} : [0, T] \times \Omega \to \mathbb{R}^3$ its time-dependent displacement field. We assume the body to be mechanically loaded in quasi-static conditions by a hard device imposing a displacement $\mathbf{u} = \bar{\mathbf{u}}(t)$ on a part of the boundary, $\partial\Omega_D$. We adopt a geometrically linear deformation theory, denoting by $\mathbf{\epsilon}(\mathbf{u}) := \nabla_{\text{sym}} \mathbf{u}$ the linearised strain tensor, which is additively decomposed into

$$ \mathbf{\epsilon} = \mathbf{\epsilon}^e + \mathbf{\epsilon}^p, \qquad (1) $$

where $\mathbf{\epsilon}^e$ denotes the elastic and $\mathbf{\epsilon}^p$ the plastic contribution. The elastic strain energy density and the strain-stress relationships are

$$ \psi_e(\mathbf{\epsilon} - \mathbf{\epsilon}^p) = \frac{1}{2} \mathbb{C}[\mathbf{\epsilon} - \mathbf{\epsilon}^p] \cdot [\mathbf{\epsilon} - \mathbf{\epsilon}^p], \quad \mathbf{\sigma} = \partial_{\mathbf{\epsilon}^e} \psi_e(\mathbf{\epsilon}^e) = \mathbb{C} \mathbf{\epsilon}^e. \qquad (2) $$

with $\mathbb{C}$ as the fourth-order elasticity tensor.

In the following, we use the decomposition of the stress tensor $\mathbf{\sigma}$ in the volumetric part, $\mathbf{\sigma}_{\text{vol}} = \sigma_H I_3$, with $\sigma_H := \text{tr} \mathbf{\sigma}/3$ as the hydrostatic stress, and the deviatoric part, $\mathbf{\sigma}_{\text{dev}} = \mathbf{\sigma} - \mathbf{\sigma}_{\text{vol}}$.

### 2.1. Modified Cam-Clay plasticity criterion with linear hardening

We consider a modified Cam-Clay plasticity criterion where the admissible stresses $\mathbf{\sigma}$ are contained in a convex domain $K(p) := \{\mathbf{\sigma} | f_y(\mathbf{\sigma}, p) \le 0\}$ defined by a yield function [40, 42–44]

$$ f_y(\mathbf{\sigma}, p) = \sigma_{eq,cc}(\mathbf{\sigma}) - \sigma_y(p), \qquad (3) $$

where $p$ is the cumulated plastic strain (to be defined later),

$$ \sigma_{eq,cc}(\mathbf{\sigma}) := \sqrt{\frac{3}{2} \|\mathbf{\sigma}_{dev}\|^2 + M^2 \sigma_H^2}, \quad \sigma_y(p) = \sigma_0 + H p \qquad (4) $$

are respectively the equivalent stress and the yield stress considering linear isotropic hardening; $M > 0$ is the hydrostatic parameter, $\sigma_0$ is the initial yield stress and $H \ge 0$ is the hardening modulus. Thus, the yield function for modified Cam-Clay plasticity depends on the von Mises equivalent stress $\sigma_{eq,VM} = \sqrt{3/2} \|\mathbf{\sigma}_{dev}\|$ and on the hydrostatic stress $\sigma_H$; hence, it takes into account the stress triaxiality characterised by $\eta^* = \sigma_H/\sigma_{eq,VM}$.
<figure><img src="image_2.png" /><figcaption>Figure 1: a) Modified Cam-Clay initial yield surface for different values of M, plotted in the plane of the equivalent stress $\sigma_{eq,VM}$ and the hydrostatic stress $\sigma_H$, normalised by the initial yield stress $\sigma_0$; b) Comparison of the modified Cam-Clay and GTN initial yield surfaces for different porosity values f with the GTN parameters $q_1 = 1.74$, $q_2 = 1.0$, $q_3 = q_1^2$. The colour code of the yield surfaces for the modified Cam-Clay model corresponds to the legend in a).</figcaption></figure>

From (3) and (4), it is evident that the function $f_y$ describes an elliptic yield surface in the $\sigma_H$-||$\sigma_{dev}$||-plane with centre in the origin, see Fig. 1a). Unlike the standard von Mises yield surface, which describes a cylinder along the hydrostatic axis in the three-dimensional principal stress space and where the yielding is independent of the hydrostatic stresses, in the modified Cam-Clay plasticity model yielding can take place in a purely hydrostatic stress state. Modified Cam-Clay plasticity thus shows some ability to account for different stress states that arise with differently shaped geometries and different loading scenarios. The yield surface shows some similarity with the one of the well-known GTN plasticity model [2, 3, 45-47], see Fig. 1b). This work considers the same behaviour in compression and tension, as is the case for the GTN model. The radii of the ellipse on the normalised $\sigma_{eq,VM}$ and $\sigma_H$ axes are equal to 1 and $1/M$, respectively; hence, the hydrostatic parameter $M > 0$, specifying the degree of influence of the hydrostatic stress, affects the shape of the yield surface as shown in Fig. 1a). In the limit $M \to 0$, one recovers the classical von Mises plasticity model. In Fig. 1b), the yield surface for different values of M is compared to the yield surface of the GTN model that shrinks for a growing porosity f.

## 2.2. Incremental variational formulation of the modified Cam-Clay plasticity model with linear hardening

We consider an associated plasticity model. Following the framework of generalised standard materials [48], we introduce the plastic dissipation potential $\pi_{K(p)}(\dot{\boldsymbol{\epsilon}}^p)$ and the convex set $K(p) := \{\boldsymbol{\sigma} | f_y(\boldsymbol{\sigma}, p) \le 0\}$ of admissible stresses [49]. Considering the yield function in (3), for the modified Cam-Clay model the rate of plastic energy density is given by

$$ \pi_{K(p)}(\dot{\boldsymbol{\epsilon}}^p) := \sup_{\boldsymbol{\sigma} \in K(p)} \boldsymbol{\sigma} \cdot \dot{\boldsymbol{\epsilon}}^p. \qquad (5) $$

Introducing the following Cam-Clay M-dependent norm for the stress and the corresponding dual norm for the plastic strain rate as

$$ \|\boldsymbol{\sigma}\|_{M} := \sigma_{eq,cc}(\boldsymbol{\sigma}) = \sqrt{\frac{3}{2} \|\boldsymbol{\sigma}_{dev}\|^2 + M^2 \sigma_{H}^2}, \quad \|\dot{\boldsymbol{\epsilon}}^p\|_{M}^{*} := \sup_{\|\boldsymbol{\sigma}\|_{M} \le \sigma_{y}(p)} \frac{\boldsymbol{\sigma}}{\sigma_{y}(p)} \cdot \dot{\boldsymbol{\epsilon}}^p = \sup_{\|\tilde{\boldsymbol{\sigma}}\|_{M} \le 1} \tilde{\boldsymbol{\sigma}} \cdot \dot{\boldsymbol{\epsilon}}^p \quad (6) $$

with $\tilde{\boldsymbol{\sigma}} := \boldsymbol{\sigma}/\sigma_y(p)$, the rate of plastic energy density can be rewritten as follows

$$ \pi_{K(p)}(\dot{\boldsymbol{\epsilon}}^p) := \sup_{\|\boldsymbol{\sigma}\|_{M} \le \sigma_y(p)} \boldsymbol{\sigma} \cdot \dot{\boldsymbol{\epsilon}}^p = \sigma_y(p) \sup_{\|\tilde{\boldsymbol{\sigma}}\|_{M} \le 1} \tilde{\boldsymbol{\sigma}} \cdot \dot{\boldsymbol{\epsilon}}^p = (\sigma_0 + Hp) \dot{p}, \quad (7) $$
where

$$ \dot{p} = \|\dot{\boldsymbol{\varepsilon}}^p\|_{\mathrm{M}}^* = \sqrt{\frac{2}{3} \|\dot{\boldsymbol{\varepsilon}}_{\mathrm{dev}}^p\|^2 + \frac{\mathrm{tr}^2 \dot{\boldsymbol{\varepsilon}}^p}{M^2}} \quad (8) $$

is the cumulated plastic strain rate. The explicit computation of the dual norm requires the solution of the minimisation problem in the definition of the dual norm, see Appendix A. The potential $\pi_{\mathcal{K}(p)}$ in (7) is composed of the dissipation potential and the hardening term, considered as the rate of plastic free energy density. Another approach would be to derive $\pi_{\mathcal{K}(p)}$ as dissipation potential only, subtracting the plastic part of the Helmholtz free energy from the rate of plastic work [50]. In this case, $\pi_{\mathcal{K}(p)}$ is additionally a function of $\dot{p}$, where $p$ is introduced as an additional state variable, see [49], and the hardening term is added to the free elastic energy density. Both derivation approaches yield the same plastic energy given below, which is the relevant term for the proposed model.

The plastic energy density is the integral in time of the rate of plastic energy density. It is a function of the cumulated plastic strain $p$ only

$$ \psi_p(p) := \int_0^t \pi_{\mathcal{K}(p)}(\dot{\boldsymbol{\varepsilon}}^p(\tau)) \mathrm{d}\tau = \sigma_0 p + \frac{1}{2} H p^2, \quad p = \int_0^t \|\dot{\boldsymbol{\varepsilon}}^p(\tau)\|_{\mathrm{M}}^* \mathrm{d}\tau. \quad (9) $$

Considering a body $\Omega$ under quasi-static loading in the form of imposed displacements on the Dirichlet part of the boundary $\partial\Omega_{\mathrm{D}}$, the total energy of the body obeying the modified Cam-Clay plasticity model is defined as the sum of the elastic energy and the plastic energy

$$ \mathcal{E}(\boldsymbol{u}, \boldsymbol{\varepsilon}^p, p) = \int_{\Omega} \psi_e(\boldsymbol{\varepsilon}(\boldsymbol{u}) - \boldsymbol{\varepsilon}^p) \mathrm{d}\boldsymbol{x} + \int_{\Omega} \psi_p(p) \mathrm{d}\boldsymbol{x}, \quad (10) $$

which is a function of the state variables $\boldsymbol{u}$, $\boldsymbol{\varepsilon}^p$, and $p$.

In the time-discrete setting, given the solution $(\boldsymbol{u}_{i-1}, \boldsymbol{\varepsilon}_{i-1}^p, p_{i-1})$ at the time step $t_{i-1}$, the cumulated plastic strain at the next time step $t_i = t_{i-1} + \Delta t$ can be expressed using a backward Euler scheme as

$$ p_i = p_{i-1} + \|\Delta \boldsymbol{\varepsilon}_i^p\|_{\mathrm{M}}^*, \quad (11) $$

where $\Delta \boldsymbol{\varepsilon}_i^p := \boldsymbol{\varepsilon}_i^p - \boldsymbol{\varepsilon}_{i-1}^p$ is the increment of the plastic strain between the two steps. The solution of the time-discrete elasto-plastic evolution problem can be obtained by solving the following incremental variational problem for the displacement $\boldsymbol{u}_i$ and the plastic strain $\boldsymbol{\varepsilon}_i^p$ at the $i$-th time step

$$ (\boldsymbol{u}_i, \boldsymbol{\varepsilon}_i^p) \in \underset{\boldsymbol{u} \in V_i, \boldsymbol{\varepsilon}^p \in \mathcal{P}}{\operatorname{argmin}} \mathcal{E}_i(\boldsymbol{u}, \boldsymbol{\varepsilon}^p), \quad (12) $$

where the total energy at time $t_i$ is given by

$$ \mathcal{E}_i(\boldsymbol{u}, \boldsymbol{\varepsilon}^p) := \int_{\Omega} \psi_e(\boldsymbol{\varepsilon}(\boldsymbol{u}) - \boldsymbol{\varepsilon}^p) \mathrm{d}\boldsymbol{x} + \int_{\Omega} \psi_p(p_{i-1} + \|\boldsymbol{\varepsilon}^p - \boldsymbol{\varepsilon}_{i-1}^p\|_{\mathrm{M}}^*) \mathrm{d}\boldsymbol{x}, \quad (13) $$

and, assuming work hardening, the spaces of the admissible displacements and plastic strain tensors are

$$ V_i := \{\tilde{\boldsymbol{u}} \in H^1(\Omega, \mathbb{R}^3) \mid \tilde{\boldsymbol{u}} = \tilde{u}_i \text{ on } \partial_{\mathrm{D}}\Omega\}, \quad \mathcal{P} := \{\tilde{\boldsymbol{\varepsilon}}^p \in L^2(\Omega, \mathbb{R}^{3,3})\}, \quad (14) $$

where $L^2$ and $H^1$ denote respectively the spaces of square integrable functions and square integrable functions with square integrable weak derivatives [15–17, 51]. From now on, for notational simplicity, we omit the subscript $i$ for all quantities related to the current time step $t = t_i$.

The time-discrete elasto-plastic evolution equations at time $t_i$ are the first-order optimality conditions of the minimisation problem (12), see Appendix B.

* The stationarity conditions with respect to the displacement field $\boldsymbol{u}$ yield the strong form equilibrium equations and boundary conditions for the Cauchy stress $\boldsymbol{\sigma}$

$$ \mathrm{div} \, \boldsymbol{\sigma} = \mathbf{0} \quad \text{in} \, \Omega, \quad \boldsymbol{\sigma} \cdot \boldsymbol{\nu} = \mathbf{0} \quad \text{on} \, \partial\Omega_{\mathrm{N}}, \quad (15) $$

where the boundary decomposes as $\partial\Omega = \partial\Omega_{\mathrm{D}} \cup \partial\Omega_{\mathrm{N}}$ and $\boldsymbol{\nu} \in \mathbb{R}^3$ denotes the unit outer normal vector.
• The stationarity conditions with respect to the plastic strain $\boldsymbol{\varepsilon}^p$ imply (see Appendix B)

$$ \Delta \boldsymbol{\varepsilon}^p = \Delta \lambda \frac{\partial f_y}{\partial \boldsymbol{\sigma}}, \quad \Delta \lambda \ge 0, \quad f_y(\boldsymbol{\sigma}, p) \le 0, \quad f_y(\boldsymbol{\sigma}, p) \Delta \lambda = 0, \qquad (16) $$

with $\Delta \lambda = \|\Delta \boldsymbol{\varepsilon}^p\|_M^*$, i.e. the classical normality rule and Karush-Kuhn-Tucker conditions of elasto-plasticity. Moreover, defining the elastic predictor stress as

$$ \boldsymbol{\sigma}^{tr} := \boldsymbol{\sigma}_{i-1} + \mathbb{C}[\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}_{i-1}] = K \operatorname{tr}(\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}_{i-1}^p) \boldsymbol{I}_3 + 2\mu[\boldsymbol{\varepsilon}_{dev} - \boldsymbol{\varepsilon}_{dev,i-1}^p], $$

where $\mu = E/[2[1+\nu]]$ denotes the shear modulus, $K = E\nu/[[(1+\nu)[1-2\nu]] + 2/3\mu$ the bulk modulus and $\nu$ the Poisson's ratio, yields

- if $\|\boldsymbol{\sigma}^{tr}\|_M \le \sigma_y(p_{i-1})$, $\Delta \boldsymbol{\varepsilon}^p = 0$;

- if $\|\boldsymbol{\sigma}^{tr}\|_M > \sigma_y(p_{i-1})$, the incremental relations (16) can be integrated by solving the following nonlinear scalar equations for $\operatorname{tr}(\Delta \boldsymbol{\varepsilon}^p)$ and $\|\Delta \boldsymbol{\varepsilon}_{dev}^p\|$

$$ \left[ \frac{1}{M^2} \left[ \frac{\sigma_y(p_{i-1})}{\|\Delta \boldsymbol{\varepsilon}^p\|_M^*} + H \right] + K \right] \operatorname{tr}(\Delta \boldsymbol{\varepsilon}^p) = \frac{1}{3} \operatorname{tr}(\boldsymbol{\sigma}^{tr}) \qquad (17) $$

$$ 2 \left[ \frac{1}{3} \left[ \frac{\sigma_y(p_{i-1})}{\|\Delta \boldsymbol{\varepsilon}^p\|_M^*} + H \right] + \mu \right] \|\Delta \boldsymbol{\varepsilon}_{dev}^p\| = \|\boldsymbol{\sigma}_{dev}^{tr}\|. \qquad (18) $$

The tensors $\boldsymbol{\sigma}$, $\Delta \boldsymbol{\varepsilon}^p$, and $\boldsymbol{\sigma}^{tr}$ are co-axial.

This return mapping algorithm was already proposed by Aravas [52].

## 2.3. Phase-field model for ductile fracture based on modified Cam-Clay plasticity

To account for damage in the material, the modified Cam-Clay plasticity model is coupled with a phase-field fracture model. The phase-field variable representing the damage in the material is denoted by $d: [0, T] \times \Omega \to [0, 1]$, whereby $d=0$ and $d=1$ respectively describe the intact and the completely damaged material states. The incremental variational functional of the coupled model extends the formulation of the plasticity model in the previous section by adding the damage contribution, following [18]. The total energy of the body at time $t_i$ reads as

$$ \mathcal{E}_i(\mathbf{u}, \boldsymbol{\varepsilon}^p, d) = \int_{\Omega} g(d) \psi_e(\boldsymbol{\varepsilon}(\mathbf{u}) - \boldsymbol{\varepsilon}^p) \, d\mathbf{x} \\ + \int_{\Omega} \rho(d) \psi_p(p_{i-1} + \|\boldsymbol{\varepsilon}^p - \boldsymbol{\varepsilon}_{i-1}^p\|_M) \, d\mathbf{x} + \int_{\Omega} \frac{G_c}{c_w} \left[ \frac{w(d)}{\ell} + \ell \|\nabla d\|^2 \right] \, d\mathbf{x}, \quad (19) $$

where the elastic and plastic energies are coupled to the damage via the degradation functions $g$ and $\rho$, respectively, and the last integral represents the damage energy. In this work, we consider standard quadratic degradation functions [13, 32], i.e. $g, \rho: [0, 1] \to \mathbb{R}$ given by

$$ g(d) = (1-d)^2, \quad \rho(d) = (1-d)^2. \qquad (20) $$

The parameter $\ell > 0$ denotes the internal length that controls the width of the damage zone, $G_c$ represents the fracture toughness and $c_w$ is a scaling parameter given by [53]

$$ c_w := 4 \int_0^1 w(s)^{\frac{1}{2}} ds. \qquad (21) $$

In the following, we only consider the AT1-model [31, 54] as it allows for an initial purely elastic and plastic stage without damage evolution. Hence, the dissipation function is defined by

$$ w: [0, 1] \to \mathbb{R}_+, \quad w(d) = d, \qquad (22) $$
resulting in $c_w = 8/3$.

To ensure the irreversibility of the damage process, the condition $d_i(x) \ge d_{i-1}(x)$ must apply for all $x \in \Omega$ at all times. We assume $d$ is sufficiently regular such that at fixed time $t_i$

$$d \in W_i := \{ \tilde{d} \in H^1(\Omega) \mid \tilde{d} \ge d_{i-1}, 0 \le \tilde{d} \le 1 \}.$$

The solution of the time-discrete elasto-plastic fracture evolution problem can again be obtained by incremental energy minimisation. Given $u_{i-1}$, $\varepsilon_{i-1}^p$, $p_{i-1}$ and the phase-field variable $d_{i-1}$ at the previous time step $t_{i-1}$, the solution of the minimisation problem at time $t_i$ is given by

$$(u, \boldsymbol{\varepsilon}^p, d) \in \underset{u \in V_i, \boldsymbol{\varepsilon}^p \in P, d \in W_i}{\operatorname{argmin}} \mathcal{E}_i(u, \boldsymbol{\varepsilon}^p, d). \quad (23)$$

First-order optimality for the problem above gives the following necessary conditions

* $u$-optimality gives the equilibrium equations (15), where the stress is affected by the phase-field variable

$$\boldsymbol{\sigma} = g(d) \partial_{\boldsymbol{\varepsilon}^e} \psi_e(\boldsymbol{\varepsilon}^e) = g(d) \mathbb{C}[\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}_p]. \quad (24)$$

* $\boldsymbol{\varepsilon}^p$-optimality gives the plastic flow rule as in section 2.2 with a damage-dependent yield stress (see [19] and Appendix B)

- if $\|\boldsymbol{\sigma}^{tr}\|_M \le \sigma_y(p_{i-1}, d)$, $\Delta\boldsymbol{\varepsilon}^p = 0$;

- if $\|\boldsymbol{\sigma}^{tr}\|_M > \sigma_y(p_{i-1}, d)$, $\Delta\boldsymbol{\varepsilon}^p$ must verify

$$\left[ \frac{\rho(d)}{M^2} \left[ \frac{\sigma_y(p_{i-1})}{\|\Delta\boldsymbol{\varepsilon}^p\|_M^*} + H \right] + g(d) K \right] \operatorname{tr}(\Delta\boldsymbol{\varepsilon}^p) = \frac{1}{3} \operatorname{tr}(\boldsymbol{\sigma}^{tr}) \quad (25)$$

$$2 \left[ \frac{\rho(d)}{3} \left[ \frac{\sigma_y(p_{i-1})}{\|\Delta\boldsymbol{\varepsilon}^p\|_M^*} + H \right] + g(d) \mu \right] \|\Delta\boldsymbol{\varepsilon}_{dev}^p\| = \|\boldsymbol{\sigma}_{dev}^{tr}\| \quad (26)$$

where

$$\sigma_y(p, d) := \rho(d) \psi'_p(p) = \rho(d) \sigma_y(p) = \rho(d) [\sigma_y(p_{i-1}) + H \|\Delta\boldsymbol{\varepsilon}^p\|_M^*]. \quad (27)$$

The optimality conditions with respect to $\boldsymbol{\varepsilon}^p$ can be summarised as Karush-Kuhn-Tucker conditions as in (16) with the damage-dependent yield function given by

$$f_y(\boldsymbol{\sigma}, p, d) = \sigma_{eq,cc}(\boldsymbol{\sigma}) - \sigma_y(p, d) = \|\boldsymbol{\sigma}\|_M - \sigma_y(p, d). \quad (28)$$

* $d$-optimality, under the irreversibility condition, gives the phase-field evolution equations (see [53]) in the form of Karush-Kuhn-Tucker complementarity conditions

$$\begin{array}{l} f_d(\boldsymbol{\sigma}, p, d) \le 0, \quad d - d_{i-1} \ge 0, \quad f_d(\boldsymbol{\sigma}, p, d) [d - d_{i-1}] = 0 \quad \text{in } \Omega, \\ \nabla d \cdot \boldsymbol{\nu} \ge 0, \quad \nabla d \cdot \boldsymbol{\nu} [d - d_{i-1}] = 0 \quad \text{on } \partial\Omega, \end{array} \quad (29)$$

where

$$f_d(\boldsymbol{\sigma}, p, d) := s'(d) \psi_e^*(\boldsymbol{\sigma}) - \rho'(d) \psi_p(p) - \frac{G_c}{c_w} \left[ \frac{1}{\ell} - 2\ell \Delta d \right] \quad (30)$$

with $s(d) := 1/g(d)$ and the complementary elastic energy density

$$\psi_e^*(\boldsymbol{\sigma}) := \frac{1}{2} \mathbb{C}^{-1} \boldsymbol{\sigma} \cdot \boldsymbol{\sigma}. \quad (31)$$
## 2.4. Numerical implementation of the phase-field model coupled with modified Cam-Clay plasticity

The numerical implementation is performed in Abaqus/Standard (2022 version) by a time-implicit coupled temperature-displacement analysis [55]. To implement the proposed material model, as suggested in [56, 57], we deploy the user subroutines UMAT and HETVAL, which contain the plasticity return map and the strong form of the evolution equation for the phase-field variable, respectively. The built-in mechanical solver of Abaqus calculates the displacement and strain increments, and the phase-field increments are determined by the thermal solver using the analogy to the heat equation. To keep the formulation variationally consistent, as proposed in [58], the irreversibility of the phase-field variable at time $t_i$ is enforced approximately by adding the penalty term

$$ \frac{\gamma}{2} \int_{\Omega} \langle d - d_{i-1} \rangle_{-}^{2} \, d\boldsymbol{x}, \qquad (32) $$

to the energy functional $\mathcal{E}_i$ in (19), where $\langle \cdot \rangle_{-} := \min(\cdot, 0)$. At time step $t_i$, the minimisers of the unconstrained optimisation problem (23), which includes the additional penalty term in (32), are found by applying a staggered solution scheme, see Box 1. To compute the displacement and the phase-field variable at $t_i$, several iterations $k$ are performed.

At time $t_i$ and iteration $k$, the UMAT subroutine computes the trial stress $\boldsymbol{\sigma}^{tr,k} = \boldsymbol{\sigma}_{i-1} + g(d^{k-1})C\Delta\boldsymbol{\varepsilon}^k$ and the yield criterion is checked by computing the yield function in (28), $f_y(\boldsymbol{\sigma}^{tr,k}, p_{i-1}, d^{k-1}) = \|\boldsymbol{\sigma}^{tr,k}\|_M - \sigma_y(p_{i-1}, d^{k-1})$. If the yield criterion is violated, the return map algorithm [52] numerically solves (25), (26) for $\mathrm{tr}(\Delta\boldsymbol{\varepsilon}^{p,k})$, $\|\Delta\boldsymbol{\varepsilon}_{dev}^{p,k}\|$. To this end, the Newton-Raphson method is applied to iteratively solve the equivalent system of equations

$$ \left( \frac{3}{2} \| \boldsymbol{\sigma}_{dev}^k \| \mathrm{tr}(\Delta\boldsymbol{\varepsilon}^{p,k}) - M^2 \frac{\mathrm{tr}(\boldsymbol{\sigma}^k)}{3} \| \Delta\boldsymbol{\varepsilon}_{dev}^{p,k} \| \right) = \mathbf{0} \qquad (33) $$

for the plastic strain increments, where the first equation results from combining (25) and (26). As function of $\mathrm{tr}(\Delta\boldsymbol{\varepsilon}^{p,k})$ and $\|\Delta\boldsymbol{\varepsilon}_{dev}^{p,k}\|$, the plastic strain tensor at iteration $k$ is given by

$$ \boldsymbol{\varepsilon}^{p,k} = \boldsymbol{\varepsilon}_{i-1}^{p} + \frac{1}{3} \mathrm{tr}(\Delta\boldsymbol{\varepsilon}^{p,k}) \mathbf{I}_3 + \|\Delta\boldsymbol{\varepsilon}_{dev}^{p,k}\| \mathbf{n}^k, \qquad (34) $$

where we define the director $\mathbf{n}^k := \Delta\boldsymbol{\varepsilon}_{dev}^{p,k}/\|\Delta\boldsymbol{\varepsilon}_{dev}^{p,k}\|$. According to (11), the cumulated plastic strain is given by

$$ p^k = p_{i-1} + \|\Delta\boldsymbol{\varepsilon}^{p,k}\|_M^* \qquad (35) $$

and the stress tensor by

$$ \boldsymbol{\sigma}^k = \boldsymbol{\sigma}^{tr,k} - g(d^{k-1}) K \mathrm{tr}(\Delta\boldsymbol{\varepsilon}^{p,k}) \mathbf{I}_3 - g(d^{k-1}) 2 \mu \|\Delta\boldsymbol{\varepsilon}_{dev}^{p,k}\| \mathbf{n}^k. \qquad (36) $$

After solving (33), the consistent tangent operator $C^{t,k}$ is updated and passed to the internal solver.

The increment of the phase-field variable $\Delta d^k$ at iteration $k$ is determined by exploiting the analogy of the strong form evolution equation for the phase-field $f_d(\boldsymbol{\sigma}, p, d) = 0$, see (30), to the heat equation [56], [57]. We first determine the complementary elastic energy density in (31) and the plastic energy density in (9) as follows

$$ \psi_e^{*,k} = \frac{\boldsymbol{\sigma}_{dev}^k \cdot \boldsymbol{\sigma}_{dev}^k}{4 \mu} + \frac{(\boldsymbol{\sigma}_h^k)^2}{2 K}, \quad \psi_p^k = \sigma_0 p^k + \frac{1}{2} H (p^k)^2 \qquad (37) $$

and then use them to compute the input required by the subroutine HETVAL, i.e. the flux

$$ q := -\ell^2 \Delta d^{k-1} = -\frac{\ell}{2} \frac{c_w}{G_c} [-s'(d^{k-1}) \psi_e^{*,k} + \rho'(d^{k-1}) \psi_p^k] - \frac{1}{2} - S_{\mathrm{opt}} \langle d^{k-1} - d_{i-1} \rangle_{-}, \qquad (38) $$
corresponding to the discrete form of $f_d(\sigma, p, d) = 0$, with the constant [58]

$$S_{opt} = \frac{\ell}{2} \frac{c_w}{G_c} \gamma = \frac{9}{16 \text{TOL}_{ir}}. \qquad (39)$$

The tolerance threshold for the irreversibility of the phase-field variable is set to TOLir = 0.001 in order to obtain a sufficiently large penalty parameter. The thermal solver in Abaqus uses the flux in (38) to determine Δd^k.

**Increment t_{i-1} → t_i:** Given Δε^k, σ_{i-1}, p_{i-1}, d_{i-1}, Δd^{k-1}:

**Start iteration k:**

For d_{i-1} ∈ [0, 1], solve 1. - 3. in UMAT subroutine:

1. Update phase-field and degradation functions g(d^{k-1}), ρ(d^{k-1}) defined by (20),

$$d^{k-1} = d_{i-1} + \Delta d^{k-1} \in [0, 1], \quad \text{IF } k = 1 \quad \Delta d^{k-1} = 0$$

2. Plasticity return map

a) Elastic predictor with trial stress $\sigma^{tr,k} = \sigma_{i-1} + g(d^{k-1}) \circ \Delta\epsilon^k$,

b) Check for plastic deformation: IF $f_y(\sigma^{tr,k}, p_{i-1}, d^{k-1}) \le 0$ given in (28),

THEN Update $\sigma^k = \sigma^{tr,k}$, $p^k = p_{i-1}$ and continue with 3.

ELSE Go to (c).

c) Apply plasticity return map: Solve (33) iteratively for tr(Δε^{p,k}), ||Δε^{p,k}_{dev}|| by the Newton-Raphson method.

Update $\sigma^k$ by (36), $p^k$ by (35) and $\sigma_y(p^k, d^{k-1})$ iteratively until the return map has converged.

d) Update the consistent tangent operator

$$C^{t,k} = \frac{\partial \Delta \sigma^k}{\partial \Delta \epsilon^k}$$

3. Compute the elastic and plastic energy densities $\psi_e^{*,k}$, $\psi_p^k$ by (37)

4. HETVAL subroutine: Given $\psi_e^{*,k}$, $\psi_p^k$, $d_{i-1}$ and $d^{k-1}$, calculate the flux in (38)

5. IF $||R^u(d^{k-1})||/\bar{a}^u \le \delta$ and $||R^d(u^k)||/\bar{a}^d \le \delta$ (see (40)) THEN

update all variables $(\cdot)_i = (\cdot)_i^k$

ELSE

Convergence not achieved. Start with iteration k + 1.

Box 1: Numerical implementation of modified Cam-Clay plasticity coupled with the phase-field model in Abaqus/Standard (2022) using a staggered approach. The tolerance for convergence of the Newton-Raphson method in the return map is set to 10⁻⁶. The tolerance of the normalised residuals is set to δ = 0.005.

At every iteration k, we solve the equilibrium equations and the phase-field evolution equation as follows

$$K^{uu}(d^{k-1}) \Delta \bar{u}^k = -R^u(d^{k-1}), \quad K^{dd}(u^k) \Delta \bar{d}^k = -R^d(u^k), \qquad (40)$$

where $K^{uu}$, $K^{dd}$ denote the submatrices on the diagonal of the Jacobian matrix of the coupled damage-displacement problem, $R^u$ and $R^d$ the corresponding right-hand side (residual) vectors, and $\bar{u}$, $\bar{d}$ the vectors of the nodal displacement and phase-field variable values, respectively. Before proceeding to the next time increment, iterations are performed until the normalised nodal residuals reach a specified tolerance, i.e. $||R^u||/\bar{a}^u < \delta_1$ and $||R^d||/\bar{a}^d < \delta_2$. The normalisation constants $\bar{a}^u$ and $\bar{a}^d$ respectively depend on the spatial and time-averaged nodal forces and phase-field fluxes ($-\ell^2 \Delta d_i^k$). The tolerance is set to $\delta_1 = \delta_2 = 0.005$.
Choosing sufficiently small (pseudo) time steps is crucial for the plasticity return map to converge during the softening stage. In the following simulations, the time step is set to $\Delta t = 10^{-4}$ s, corresponding to $10^4$ increments. In addition, a sufficient number of staggered iterations must be allowed for, especially in the plastic and softening stages and at the beginning of the analysis to find the initial phase-field increment. In the following simulations, we need up to about 230 iterations per time increment, depending on the specific input parameters, time step and geometry.

# 3. Homogeneous material responses

In this section, analytical material responses are derived for the proposed model under the assumption of homogeneous damage, starting in one dimension and continuing with the multi-dimensional case.

## 3.1. One-dimensional analytical solutions

In the one-dimensional case, the yield function reads

$$f_y(\sigma, \varepsilon^p, d) = \sigma - \rho(d) [\sigma_0 + H \varepsilon^p], \quad (41)$$

where $\sigma = g(d) E (\varepsilon - \varepsilon^p)$ and $E > 0$ denotes the Young's modulus. We assume monotonic loading and thus directly consider the plastic strain $\varepsilon^p \in \mathbb{R}$ instead of the cumulated plastic strain. The hydrostatic parameter $M$ does not play a role in the one-dimensional case, and the model in this particular setting reduces to von Mises plasticity. Due to the condition $f_y \Delta\lambda = 0$ in (16), for evolving plastic strain $\Delta\lambda > 0$ it is $f_y = 0$. The corresponding value of stress, which we refer to as yield stress $\sigma_y$, is then given by

$$\sigma_y(\varepsilon^p, d) = \rho(d) [\sigma_0 + H \varepsilon^p]. \quad (42)$$

The damage yield function is given by (30) and is expressed for a homogeneous damage state as

$$f_d(\sigma, \varepsilon^p, d) = s'(d) \frac{1}{2} E^{-1} \sigma^2 - \rho'(d) \psi_p(\varepsilon^p) - \frac{G_c}{c_w} \frac{1}{\ell}. \quad (43)$$

Due to the condition $f_d [d - d_{i-1}] = 0$ in (29), for evolving damage $d > d_{i-1}$ it is $f_d = 0$. The corresponding value of stress, which we refer to as damage yield stress $\sigma_d$, is then given by

$$\sigma_d(\varepsilon^p, d) = \sqrt{E [1-d]^3 \left[ \rho'(d) \psi_p(\varepsilon^p) + \frac{G_c}{c_w \ell} \right]}. \quad (44)$$

For the AT1-model, $\sigma_d$ reaches its maximum value in the elastic stage; this value is given by

$$\sigma_d(0, 0) = \sqrt{\frac{3 E G_c}{8 \ell}}, \quad (45)$$

see also the derivations in [59]. Considering again the damage yield function in (43) and the quadratic degradation functions in (20), by solving $f_d(\sigma, \varepsilon^p, d) = 0$ for the phase-field variable we obtain

$$d(\varepsilon, \varepsilon^p) = 1 - \frac{\frac{G_c}{c_w \ell}}{E [\varepsilon - \varepsilon^p]^2 + 2 \psi_p(\varepsilon^p)}. \quad (46)$$

The initiation of damage depends on the constant

$$\zeta = \frac{1}{\sigma_0} \sqrt{\frac{3 E G_c}{8 \ell}} > 0. \quad (47)$$

For $\zeta < 1$, it is $\sigma_d(0, 0) < \sigma_y(0, 0)$, i.e. the damage yield criterion is reached before the plastic yield criterion during the loading process [20]. Instead, for $\zeta \ge 1$, it is $\sigma_y(0, 0) \le \sigma_d(0, 0)$, thus the plastic
yield criterion is fulfilled first. The behaviour of $\sigma_y$ and $\sigma_d$ is displayed in Fig. 2 for different values of $\zeta$. In the following, we focus on examples where $\zeta \ge 1$, as this reflects the experimentally observed ductile failure behaviour. In this case, a plastic stage is reached, where the plastic strain evolves before the phase-field variable starts to grow. Increasing $\zeta$ leads to a higher material ductility [60]. If all other parameters are kept fixed, a higher value of $\zeta$ is associated with a higher fracture toughness, which increases the initial maximum value of the damage stress $\sigma_d(0, 0)$ in (45) and hence extends the plastic stage. As a related effect, the higher fracture toughness also leads to a slower phase-field evolution in the softening stage, which becomes visible by a less rapid decrease of the stress in Fig. 2. The slope of the hardening branch appears small because the hardening parameter $H$ is relatively small. However, we use this value as it is of the same order of magnitude used for the simulations presented in section 4.

<figure><img src="image_3.png" /><figcaption>Figure 2: Damage yield stress $\sigma_d$ and plastic yield stress $\sigma_y$, normalised by the initial yield stress $\sigma_0$, as functions of the strain for different values of $\zeta$ given in (47). The material parameters are $E = 210$ GPa, $\sigma_0 = 730$ MPa, $H = 500$ MPa, $\ell = 0.15$ mm.</figcaption></figure>

Assuming $\zeta \ge 1$, the homogeneous ductile material response of the proposed model in the one-dimensional case can be divided into the following three stages. These stages are illustrated in Fig. 3, displaying the corresponding evolution of the plastic strain $\varepsilon^p$, the phase-field variable $d$ and the stress-strain curve.

(i) Linear-elastic stage (E), which is guaranteed by the AT1-model, as the damage only starts to evolve after a certain threshold value of the strain/stress. The maximum stress in this stage is given by the initial yield stress $\sigma_0$. The plastic yield stress $\sigma_y$ and damage yield stress $\sigma_d$ remain constant during the elastic stage, see Fig. 3b). We have

$$ \varepsilon \in [0, \frac{\sigma_0}{E}], \quad \varepsilon^p \equiv 0, \quad d \equiv 0, \quad \sigma = E \varepsilon; $$

(ii) Plastic stage with isotropic linear hardening (P) in a strain interval of

$$ \varepsilon \in \left[ \frac{\sigma_0}{E}, \left( \frac{1}{E} + \frac{1}{H} \right) \sqrt{ \frac{E \sigma_0^2 + E H \frac{G_c}{c_w \ell}}{H + E} - \frac{\sigma_0}{H} } \right]. $$

The plastic yield criterion $f_y(\sigma, \varepsilon^p, d) = 0$ is fulfilled in this stage, and plastic strain starts to evolve. Due to the plastic strain evolution, $\sigma_d$ decreases during the plastic stage, while the damage yield criterion is still not met since we consider the case $\zeta \ge 1$. The plastic strain, the phase-field variable and the stress are given by

$$ \varepsilon^p = \frac{E \varepsilon - \sigma_0}{E + H}, \quad d \equiv 0, \quad \sigma = \sigma_y(\varepsilon^p, d); $$
(iii) When the plastic yield stress $\sigma_y$ and the damage yield stress $\sigma_d$ intersect at a plastic strain $\hat{\epsilon}^p > 0$, i.e.

$$ \sigma_y(\hat{\epsilon}^p, 0) = \sigma_d(\hat{\epsilon}^p, 0), \qquad (48) $$

the phase-field variable $d$ starts to evolve and leads to material softening [20]. This stage occurs in a strain interval of

$$ \varepsilon \in \left[ \left( \frac{1}{E} + \frac{1}{H} \right) \sqrt{ \frac{E \sigma_0^2 + E H \frac{G_c}{c_w \ell}}{H + E} - \frac{\sigma_0}{H} } , \infty \right]. $$

The minimum of the strain interval is determined by $\varepsilon = \varepsilon^e + \varepsilon^p$ with the elastic strain given by $\varepsilon^e = \sigma_y(\varepsilon^p, 0)/E$ and the plastic strain given by the condition in (48). A distinction is made between damage with constant plastic strain (D) and coupled plasticity and damage (PD). With regard to practical applications with pronounced ductile material behaviour, in the examples of this work we focus on the (E-P-PD) response, where damage and plasticity both evolve after an elastic and a purely plastic stage. The choice of identical degradation functions $g(d) = \rho(d)$ precludes the (E-P-D) response.

In the case of the plastic-damage stage (PD), plasticity continues to evolve when the phase-field $d$ starts to evolve [20], since the yield criteria $f_y(\sigma, \varepsilon^p, d) = 0$ and $f_d(\sigma, \varepsilon^p, d) = 0$ both apply. We have

$$ \varepsilon^p = \frac{E \varepsilon - \sigma_0}{E + H}, \quad d(\varepsilon, \varepsilon^p) = 1 - \frac{\frac{G_c}{c_w \ell}}{\frac{\sigma_y(\varepsilon^p, 0)^2}{E} + 2 \psi_p(\varepsilon^p)}, \quad \sigma = \sigma_y(\varepsilon^p, d). $$

As shown in Fig. 3b), for the given parameters, the condition in (48) is fulfilled with $\hat{\epsilon}^p = 0.015$ at a total strain of $\varepsilon = 0.018$. That is, the damage yield stress $\sigma_d$ meets the plastic yield stress $\sigma_y$ at $\sigma/\sigma_0 = 1.01$, and this marks the beginning of softening, where $\sigma_d(\varepsilon^p, d) = \sigma_y(\varepsilon^p, d)$.

<figure><img src="image_4.png" /><figcaption>Figure 3: Example of an E-P-PD response, generated with the standard degradation functions $g(d) = (1-d)^2$, $\rho(d) = (1-d)^2$, $E = 210$ GPa, $\sigma_0 = 730$ MPa, $H = 500$ MPa, $\ell = 0.15$ mm, $G_c = 10$ N/mm; a) plastic strain and phase-field evolution with respect to strain, b) stress-strain response considering the normalised stress, damage yield stress and plastic yield stress.</figcaption></figure>

## 3.2. Multi-dimensional material behaviour

In the following, we analyse the stress evolution for the three-dimensional homogeneous case under uniaxial tension. The stress triaxiality is accordingly $\eta^* = 1/3$. The material constants are given by $E = 210$ GPa, $\nu = 0.3$, $M = 0.5$ and $G_c = 150$ N/mm. For a given constant strain increment, the stresses are computed numerically as described in section 2.4. The resulting material response is shown in Fig. 4, where we illustrate the stress evolution in the hydrostatic-deviatoric plane for a monotonic loading and the corresponding plastic and damage yield surfaces $f_y(\sigma, p, d)$ and $f_d(\sigma, p, d)$.
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4: Stress evolution for a given strain path with the corresponding plastic yield surface $f_y(\sigma, p, d)$, where the initial plastic yield surface $f_y(\sigma, 0, 0)$ is marked in blue, and damage yield surface $f_d(\sigma, p, d)$ in the hydrostatic-deviatoric plane. The evolution of the stress and the yield surfaces is shown in a) for a hardening parameter of $H = 500$ MPa and in b) for a hardening parameter of $H = 5000$ MPa. The numbers represent the sequence of stress evolution under the loading process.</figcaption></figure>
Fig. 4a) shows the material response for a hardening parameter $H = 500$ MPa, which is in the order of magnitude of the hardening parameter used for the simulations in section 4. To better visualise the initial enlargement of the plastic yield surface during the plastic stage, the hardening parameter is increased to $H = 5000$ MPa in Fig. 4b). The computed stress state is marked with the numbers next to it indicating the order in which these stresses occur during loading. The plastic yield surface $f_y(\sigma, p, d)$ first grows and then shrinks (corresponding to a development from black to white in Fig. 4), and the damage yield surface $f_d(\sigma, p, d)$ shrinks (development from orange to red) during the loading process.

As shown in Fig. 4, the stresses are initially located in the elastic domain, the corresponding initial yield surface is $f_y(\sigma, 0, 0)$ (marked in blue). Plastic strain starts to evolve as soon as the plastic yield criterion is met, i.e. the stress is on the yield surface, and the yield surface then grows with increasing plastic strain. When both the plasticity and the damage criterion are met, i.e. $f_y(\sigma, p, 0) = 0$, $f_d(\sigma, p, 0) = 0$, the stress state evolves along the intersection of the two shrinking yield surfaces. More precisely, damage and plasticity continue to evolve together, and the stress decreases. We notice a slow plastic evolution before the onset of damage, but a rapid decrease in stresses as soon as damage and plasticity evolve together, which slows down with a higher value of the phase-field variable $d$.

# 4. Numerical simulations

Simulation results for different geometries are presented using the algorithm described in section 2.4. In the following simulations, the Young's modulus is set to $E = 210$ GPa and the Poisson's ratio to $\nu = 0.3$. The plastic parameters typical for steel and the fracture toughness $G_c$ are specified individually for the geometries under consideration.

## 4.1. Axisymmetric specimen under tension

We first consider an axisymmetric specimen under tension. The specimen is clamped at one end and subjected to a uniform displacement $\bar{u} = (U, 0, 0)$ at the other end. The overall length is set to $L_0 = 1$ mm, and the internal length to $\ell = 0.15$ mm. The specimen is discretised by four-node bilinear axisymmetric elements with a uniform mesh size of $h_e = 5 \cdot 10^{-3}$ mm, which provides a sufficiently high resolution to display a smooth distribution of the phase-field variable and of the cumulated plastic strain over the length of the bar. The specimen has a small imperfection, namely, the diameter of the cross-section is decreased by $2 \cdot 10^{-4}$ mm in the centre and gradually increases to the full diameter of $D_0 = 0.1$ mm over a total length of 0.2 mm. The element size is adjusted accordingly in the neighbourhood of the small defect. The plastic parameters are defined by an initial yield stress of $\sigma_0 = 730$ MPa, a hydrostatic parameter $M = 0.5$ and a hardening parameter $H = 500$ MPa, and the fracture toughness is given by $G_c = 10$ N/mm.

The stress $\sigma = F/A(D_0)$ in the axisymmetric specimen depending on the evolution of the strain $\varepsilon = U/L_0$ is shown in Fig. 5, where $A$ denotes the surface of the cross-section with a diameter of $D_0 = 0.1$ mm. The material response describes an elastic stage, a purely plastic stage with linear isotropic hardening and a softening stage, where the stress in the material decreases. The onset of softening at a strain of $\varepsilon = 0.014$ and a stress of $\sigma = 728.2$ MPa coincides with the onset of damage, i.e. $d > 0$, in the centre of the specimen. Compared to the one-dimensional example in section 3.1, here localisation of the plastic strain occurs and a three-dimensional stress state arises where the effect of the hydrostatic stress leads to the onset of softening at a lower strain. The evolving phase-field profile over the length of the bar is illustrated in Fig. 6, showing the finite support of the phase-field variable. The damage is regularised over a damage zone proportional to $\ell$ [20, 58], while plasticity is strongly localised in the centre of the damage zone where the imperfection is located, as illustrated in Fig. 7.
<figure><img src="image_7.png" /><figcaption>Figure 5: Stress-strain response of the axisymmetric specimen for $G_c = 10$ N/mm.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: Evolution of the distribution of the phase-field variable $d$ over the length of the axisymmetric specimen along the symmetry axis.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7: Evolution of the distribution of the cumulated plastic strain $p$ over the length of the axisymmetric specimen along the symmetry axis.</figcaption></figure>

## 4.2. Two-dimensional plane strain specimens

Two differently shaped two-dimensional plane strain specimens are considered to investigate the damage initialisation and evolution as well as its connection to the deviatoric and hydrostatic stress. The specimens are discretised by four-node bilinear plane strain elements. The geometry of the first specimen, similar to a geometry used in [60], is shown in Fig. 8. The specimen with a length of $L_0 = 6$ mm is clamped on the left-hand side, and a uniform displacement $U$ is applied on the right-hand side. The mesh is refined in the centre, where plastic strain and damage are expected to arise, with a uniform mesh size of $h_e = 0.02$ mm.

The internal length of the phase-field model is set to $\ell = 0.15$ mm. The material parameters are $\sigma_0 = 730$ MPa, $M = 0.5$, $H = 200$ MPa, and the fracture toughness is $G_c = 30$ N/mm.

Fig. 9 illustrates the evolution of the phase-field variable $d$ and the cumulated plastic strain $p$ at different average strain levels $\varepsilon = U/L_0$. Damage occurs first in the centre of the specimen, where the hydrostatic stress $\sigma_H$ concentrates, see Fig. 10. As shown in Fig. 11, the deviatoric part of the equivalent stress in (4), represented by $\sigma_{eq,VM}$, is distributed over the whole gauge area. The cumulated plastic strain $p$ in Fig. 9 b) concentrates in shear bands crossing in the centre of the specimen. The

<figure><img src="image_10.png" /><figcaption>Figure 8: Geometry, boundary conditions and mesh of a plane strain specimen with an element size of $h_e = 0.02$ mm in the centre. For visualisation purposes, the mesh size is increased by a factor of three.</figcaption></figure>
phase-field variable continues to evolve in one of the shear bands, in which $p$ concentrates and which becomes predominant over the other.

<figure><img src="image_11.png" /><figcaption>Figure 9: a) Evolution of the phase-field variable $d$ with the average strain $\varepsilon = U/L_0$, compared to b) the evolution of the cumulated plastic strain $p$. The images in the same row are related to the same average strain of $\varepsilon = 0.006$, $\varepsilon = 0.007$, $\varepsilon = 0.008$ and $\varepsilon = 0.012$, respectively.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 10: Distribution of the hydrostatic stress $\sigma_{\text{h}}$ at an average strain of $\varepsilon = 0.004$. The hydrostatic stress concentrates in the centre, where damage formation is triggered.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 11: Distribution of the deviatoric part of the stress tensor represented by the von Mises equivalent stress $\sigma_{\text{eq,VM}}$ at an average strain of $\varepsilon = 0.004$.</figcaption></figure>

The global behaviour of the specimen is shown in Fig. 12, where the nominal stress is expressed as a function of the average strain. The nominal stress is given by $\sigma = F/A(D_0)$, where $A$ denotes the surface of the centre cross-section with a width of $D_0 = 0.9$ mm and a thickness of 1 mm. We compare two curves corresponding to two different values of the hydrostatic parameter $M$. For a smaller $M$, a higher deviatoric stress $\sigma_{\text{dev}}$ is required to fulfil the yield criterion, resulting in a higher overall stress level. For a critical element in the centre of the specimen, the local phase-field evolution and the local equivalent stress $\sigma_{\text{eq,cc}}$ are illustrated in Fig. 13, showing that the onset of stress softening coincides with the onset of phase-field evolution. The equivalent stress in the purely plastic stage is independent of the choice of the parameter $M$ and is equal to the yield stress $\sigma_y$. However, for a smaller $M$ the onset of damage begins at a larger global strain, which also affects the evolution of the equivalent stress $\sigma_{\text{eq,cc}}$.
<figure><img src="image_14.png" /><figcaption>Figure 12: Global material response of the specimen shown in Fig. 8 for two different values of M.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 13: Local response of one critical element in the specimen centre. The equivalent stress $\sigma_{eq,cc}$ and the phase-field variable $d$ are shown versus the average strain for two different values of M.</figcaption></figure>

Although the plastic strain is localised in shear bands, the damage-plasticity problem is well-posed, as the damage is regularised by the phase-field. Even in the plastic limit case, when the plastic strain becomes a bounded Radon measure, the discretised measure is mesh-independent [61],[20]. In Fig. 14, the mesh in-dependency of the proposed model is verified. For this purpose, different mesh refinements in the centre of the specimen are considered, more precisely $h_e = 0.005$ mm, $h_e = 0.02$ mm, $h_e = 0.04$ mm, $h_e = 0.06$ mm, where the internal length remains fixed at $\ell = 0.15$ mm. The material response, i.e. the solutions for the displacement and the phase-field variable, converges with mesh refinement towards the exact solution, which is best approximated by the finest mesh, here given by $h_e = 0.005$ mm, which indicates mesh independency.

Another plane strain specimen with pre-defined notches in the centre is shown in Fig. 15. The initial

length is set to $L_0 = 6$ mm, and the constitutive parameters, length scale parameter, element type, and mesh size are the same as in the previous example. The specimen shows a similar ductile fracture behaviour as the one discussed above. Fig. 16 illustrates the distribution of the phase-field variable $d$ and the cumulated plastic strain $p$ at an average strain of $\varepsilon = 0.012$. As illustrated in Fig. 16 a), the maximum value of the phase-field variable is reached in the centre, where a non-straight regularised crack path can be observed between the two notches. The cumulated plastic strain $p$ in Fig. 16 b) localises in the centre and splits into two branches near the notches, while the phase-field variable concentrates in one of the branches. As illustrated in Fig. 16, damage is regularised by the phase-field approach, however, the plastic strain can accumulate locally and reach large values in the shear bands.

<figure><img src="image_16.png" /><figcaption>Figure 14: Global stress-strain curves of the specimen shown in Fig. 8 for different mesh refinements with corresponding element size $h_e$.</figcaption></figure>
<figure><img src="image_17.png" /></figure>

Figure 15: Geometry, boundary conditions and mesh of a notched plane strain specimen. The specimen has two notches with a radius of 0.5 mm and a clamped boundary on the left. The mesh size is set to $h_e = 0.02$ mm in the centre. For visualisation purposes, the mesh size is increased by a factor of four.

Figure 16: Distribution of a) the phase-field variable $d$ and b) the cumulated plastic strain $p$ at an average strain of $\varepsilon = 0.012$.

### 4.3. Axisymmetric notched specimens under tension

In this section, we apply the proposed phase-field model to two axisymmetric specimens with a different notch radius at the centre to demonstrate the influence of the stress triaxiality, i.e. the hydrostatic stresses, on the ductile fracture behaviour. The specimen with a notch radius of 1 mm is shown in Fig. 17 a) and the specimen with a notch radius of 2 mm in Fig. 17 b). The boundary on the left-hand side is clamped, while the specimens are subjected to a uniform displacement $U$ on the right-hand side. The total length is given by $L_0 = 35$ mm with a minimum diameter in the centre of 3 mm. The specimens are discretised by four-node bilinear axisymmetric elements. The mesh size in the centre is set to $h_e = 0.05$ mm, and the internal length to $\ell = 0.15$ mm. The hydrostatic parameter is defined by $M = 0.8$ to ensure a significant influence of the hydrostatic stress. The initial yield stress is given by $\sigma_0 = 730$ MPa, the hardening parameter by $H = 200$ MPa, and the fracture toughness is set to $G_c = 300$ N/mm.

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

Figure 17: Geometry, boundary conditions and mesh, mirrored on the symmetry axis, of an axisymmetric specimen with a) notch radius of 1 mm, b) notch radius of 2 mm and a diameter of 5 mm, reduced to 3 mm in the centre. For visualisation purposes, the mesh size is increased by a factor of three.

Figure 18: Stress triaxiality $\eta^*$ for the specimen with a) notch radius of 1 mm, b) notch radius of 2 mm at an average strain of $\varepsilon = 0.0025$.

With this example, we aim at showing that the proposed variational model can account for varying stress triaxialities that occur in the differently notched specimens. The stress triaxiality $\eta^* = \sigma_H/\sigma_{eq,VM}$ for the two specimens is shown in Fig. 18 at an average strain $\varepsilon = U/L_0 = 0.0025$. This strain corresponds to the plastic stage before the phase-field variable starts to evolve. In the critical element in the centre, the specimen with a notch radius of 1 mm in Fig. 18 a) shows a stress triaxiality of $\eta^* = 1.62$, whereas the specimen with the larger notch radius of 2 mm shows a lower stress triaxiality of $\eta^* = 1.52$. The inhomogeneous distribution of the hydrostatic stress has an influence on the location of the crack initiation [62], resulting in fracture initiation in the centre of the specimens, and the amount of stress triaxiality is expected to have an influence on the global stress-strain curve, in particular on
the onset of material softening.

Fig. 19 shows the stress-strain curves. The specimen with a notch radius of 1 mm shows an earlier material softening than the one with a notch radius of 2 mm. This is due to the incorporation of the hydrostatic stress in the proposed model coupling plasticity and the phase-field variable, which thus takes into account the different stress triaxiality of the two specimens for modelling fracture initiation. Fig. 20 illustrates the local hydrostatic stress $\sigma_H$ and the von Mises equivalent stress $\sigma_{eq,VM}$ in the critical element in the centre of both specimens during the purely elastic and plastic stage before damage initiation. Compared to the specimen with notch radius of 2 mm, the specimen with notch radius of 1 mm shows a higher hydrostatic stress $\sigma_H$ and a lower von Mises equivalent stress $\sigma_{eq,VM}$.

<figure><img src="image_19.png" /><figcaption>Figure 19: Global stress-strain curve for the specimen with notch radius of 1 mm (solid black) and with notch radius of 2 mm (dashed gray).</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 20: Local evolution of the hydrostatic stress $\sigma_H$ and the von Mises equivalent stress $\sigma_{eq,VM}$ in the critical element versus the average strain for the specimen with notch radius of 1 mm (solid line) and with notch radius of 2 mm (dashed line).</figcaption></figure>

The reason for the different fracture initiation and thus material softening of the two specimens is analysed in Fig. 21 and Fig. 22. The cumulated plastic strain $p$ and the phase-field $d$ in the critical element in the centre of both specimens are shown in Fig. 21. The kink in the evolution of the cumulated plastic strain occurs because, as the phase-field develops, plasticity becomes localised in two rows of elements where a significant increase in plastic strain is observed. The phase-field starts to evolve in the two specimens at a different strain; in the critical element of the specimen with a notch radius of 1 mm at a lower average strain of $U/L_0 = 0.0026$, in the critical element of the other specimen at an average strain of $U/L_0 = 0.0031$. Fig. 22 shows a detail of Fig. 21 considering the strain interval $U/L_0 = [0, 0.0035]$, as this shows the different evolution of $p$ in the specimens that triggers the damage initiation. The phase-field $d$ starts to develop in both specimens when the cumulated plastic strain reaches the value $p = 0.14$. Thus, the different growth of the cumulated plastic strain in the purely plastic stage leads to damage initiation at a different average strain.

The specific choice of $M$ plays an important role in weighting the hydrostatic stress or strain against the deviatoric components. It can be observed that a smaller $0 \le M < 1$ has a greater influence on the hydrostatic part of the plastic strain tensor in (9), while reducing the influence on the equivalent stress $\sigma_{eq,cc}(\sigma)$ defined in (4) at the same time. Since $M$ determines the influence of the hydrostatic stress on the equivalent stress during plastic deformation, it is the parameter that, in addition to the hydrostatic stress arising in the elastic stage, causes a different stress-strain behaviour with different triaxiality $\eta^*$. For $M$ close to zero, the hydrostatic stress hardly influences the equivalent stress $\sigma_{eq,cc}$ during plastic yielding. In this case, the different stress levels visible in the stress-strain curves of the two specimens are only due to triaxiality effects in the elastic stage caused by the predefined notches, resulting in different hydrostatic stresses and volumetric elastic strains. However, ductile fracture initiation cannot be influenced by triaxiality effects for $M$ close to zero.
<figure><img src="image_21.png" /><figcaption>Figure 21: Local evolution of the cumulated plastic strain *p* and the phase-field *d* versus the average strain in the critical element for the specimen with notch radius of 1 mm (solid line) and with notch radius of 2 mm (dashed line).</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 22: Detail of Fig. 21 illustrating the purely elastic and plastic stage until the onset of damage.</figcaption></figure>

Fig. 23 shows the stress-strain curves of the geometry with a notch radius of 1 mm for different mesh sizes. The mesh size in the centre is set to $h_e = 0.05$ mm, $h_e = 0.1$ mm or $h_e = 0.2$ mm with an internal length of $\ell = 0.3$ mm, so that $h_e < \ell$ applies to all three simulations. The stress-strain curves underline again the mesh independency of the proposed model, see also section 4.2.

<figure><img src="image_23.png" /><figcaption>Figure 23: Stress-strain curves of the specimen with notch radius of 1 mm for different mesh refinements with corresponding element size $h_e$.</figcaption></figure>

# 5. Conclusion and perspectives

This paper demonstrates the potential of using the modified Cam-Clay plasticity model originating from soil mechanics to describe ductile material behaviour, for example of metals. The underlying physical processes, such as the formation of voids in metals, lead to ductile fracture which depends on triaxiality effects, in turn due to a three-dimensional stress state resulting from necking or notches. In this context, the magnitude and distribution of the hydrostatic stress is crucial for the fracture initiation. The yield surface of the modified Cam-Clay model is an ellipse along the hydrostatic axis in the principal stress space instead of a cylinder as for von Mises plasticity, so that the volumetric part of the stress and strain tensor is considered in addition to the deviatoric part during plastic deformation. This work represents a step towards a variationally consistent formulation that accounts for the triaxiality effects on ductile fracture initiation by coupling modified Cam-Clay plasticity with a regularising phase-field
model and ensures mesh independence.

We have shown that the proposed coupled model can reproduce the ductile fracture behaviour with its typical successive occurrence of elastic, plastic and softening phases through analytical calculations and various numerical simulations. The triaxiality effects are captured qualitatively so that the fracture initiation at different global strains can be modelled for geometries with different stress triaxiality. However, the model shows limited flexibility to adapt the magnitude of the triaxiality effects to the range of possible material responses. Therefore, it is necessary to further increase the influence of the stress triaxiality beyond the adjustment of the parameter *M*.

From the numerical perspective, an improvement of the convergence behaviour of the staggered algorithm using the implicit solver is desirable. This behaviour is currently sensitive to many interacting factors, such as the Newton-Raphson iterations in the plasticity return map, the penalty term, the loading process, the fracture toughness, the geometry and the chosen (pseudo) time step. Furthermore, we suggest to add geometric nonlinearity to the proposed model in order to predict the ductile fracture behaviour in real experiments.

To conclude, the limitations of the model are emphasised. The limited adaptability of the model to different material responses due to only one parameter, which enters the formulation through the yield surface and adjusts the effect of hydrostatic stresses and thus the onset of material softening, demands further work in this direction. In addition, a comparison of the simulation results with experimental tests for ductile fracture requires the extension of the variational model to large strains.

## Acknowledgments

Anne-Sophie Sur acknowledges the funding by the 'Norwegian CCS Research Centre' (NCCS), Centre for Environment-friendly Energy Research (Project No. 257579), the 'Centre for Advanced Structural Analysis' (CASA), Centre for Research-based Innovation (Project No. 237885), the financial support of 'Norges tekniske høgskoles fond' for a research stay, and wishes to thank David Morin (SIMLab, NTNU) for his support with the Abaqus implementation. Laura De Lorenzis acknowledges funding from the Swiss National Science Foundation through Grant No. 200021-219407 'Phase-field modeling of fracture and fatigue: from rigorous theory to fast predictive simulations'.

## References

[1] J. R. Rice, D. M. Tracey. 'On the ductile enlargement of voids in triaxial stress fields'. In: *Journal of the Mechanics and Physics of Solids* 17.3 (1969), pp. 201–217. DOI: 10.1016/0022-5096(69) 90033-7.

[2] A. L. Gurson. 'Continuum Theory of Ductile Rupture by Void Nucleation and Growth: Part I-Yield Criteria and Flow Rules for Porous Ductile Media'. In: *Journal of engineering materials and technology* 99.1 (1977), pp. 2–15. DOI: 10.1115/1.3443401.

[3] V. Tvergaard, A. Needleman. 'Analysis of the cup-cone fracture in a round tensile bar'. In: *Acta Metallurgica* 32.1 (1984), pp. 157–169. DOI: 10.1016/0001-6160(84)90213-X.

[4] G. R. Johnson, W. H. Cook. 'Fracture characteristics of three metals subjected to various strains, strain rates, temperatures and pressures'. In: *Engineering Fracture Mechanics* 21.1 (1985), pp. 31–48. DOI: 10.1016/0013-7944(85)90052-9.

[5] Y. Bai, T. Wierzbicki. 'Application of extended Mohr-Coulomb criterion to ductile fracture'. In: *International Journal of Fracture* 161.1 (2010), pp. 1–20. DOI: 10.1007/s10704-009-9422-8.

[6] M.G. Cockcroft, D.J. Latham. 'Ductility and the workability of metals'. In: *Journal of the Institute of Metals* 96 (1968), pp. 33–39.
[7] G. Gruben, O. S. Hopperstad, T. Børvik. 'Evaluation of uncoupled ductile fracture criteria for the dual-phase steel Docol 600DL'. In: *International Journal of Mechanical Sciences* 62.1 (2012), pp. 133–146. DOI: 10.1016/j.ijmecsci.2012.06.009.

[8] D. Mohr, S. J. Marcadet. 'Micromechanically-motivated phenomenological Hosford-Coulomb model for predicting ductile fracture initiation at low stress triaxialities'. In: *International Journal of Solids and Structures* 67-68 (2015), pp. 40–55. DOI: 10.1016/j.ijsolstr.2015.02.024.

[9] Davide Bigoni. *Nonlinear Solid Mechanics: Bifurcation Theory and Material Instability*. Cambridge: Cambridge University Press, 2012.

[10] A. Benallal, R. Billardon, G. Geymonat. 'Bifurcation and Localization in Rate-Independent Materials. Some General Considerations'. In: *Bifurcation and Stability of Dissipative Systems*. Ed. by Q. S. Nguyen. Vol. 327. International Centre for Mechanical Sciences. Vienna: Springer, 1993, pp. 1–44.

[11] J. Besson. 'Continuum Models of Ductile Fracture: A Review'. In: *International Journal of Damage Mechanics* 19.1 (2010), pp. 3–52. DOI: 10.1177/1056789509103482.

[12] R. de Borst, C.V. Verhoosel. 'Damage, Material Instabilities, and Failure'. In: *Encyclopedia of Computational Mechanics Second Edition: Solids and Structures, Part 2*. Ed. by R. Borst E. Stein, T.J.R. Hughes. Vol. 4. Chichester: Wiley, 2017, pp. 1–50. DOI: 10.1002/9781119176817.ecm2035.

[13] B. Bourdin, G.A. Francfort, J.-J. Marigo. 'Numerical experiments in revisited brittle fracture'. In: *Journal of the Mechanics and Physics of Solids* 48.4 (2000), pp. 797–826. DOI: 10.1016/S0022-5096(99)00028-9.

[14] G.A. Francfort, J.-J. Marigo. 'Revisiting brittle fracture as an energy minimization problem'. In: *Journal of the Mechanics and Physics of Solids* 46.8 (1998), pp. 1319–1342. DOI: 10.1016/S0022-5096(98)00034-9.

[15] G. Dal Maso, A. DeSimone, M. G. Mora. 'Quasistatic evolution problems for linearly elastic-perfectly plastic materials'. In: *Archive for rational mechanics and analysis* 180.2 (2006), pp. 237–291. DOI: 10.1007/s00205-005-0407-0.

[16] P.-M. Suquet. 'Sur un nouveau cadre fonctionnel pour les équations de la plasticité'. In: *CR Acad. Sci. Paris Sér. AB* 286.23 (1978).

[17] G. A. Francfort, A. Giacomini, J.-J. Marigo. 'The elasto-plastic exquisite corpse: A Suquet legacy'. In: *Journal of the Mechanics and Physics of Solids* 97. Supplement C (2016), pp. 125–139. DOI: 10.1016/j.jmps.2016.02.002.

[18] R. Alessi, J.-J. Marigo, S. Vidoli. 'Gradient damage models coupled with plasticity and nucleation of cohesive cracks'. In: *Ach. Rational Mech. Anal.* 214 (2014), pp. 575–615. DOI: 10.1007/s00205-014-0763-8.

[19] R. Alessi, J.-J. Marigo, S. Vidoli. 'Gradient damage models coupled with plasticity: Variational formulation and main properties'. In: *Mech. Mater.* 80, Part B (2015), pp. 351–367. DOI: 10.1016/j.mechmat.2013.12.005.

[20] R. Alessi, J.-J. Marigo, C. Maurini, S. Vidoli. 'Coupling damage and plasticity for a phase-field regularisation of brittle, cohesive and ductile fracture: One-dimensional examples'. In: *International Journal of Mechanical Sciences* 149 (2018), pp. 559–576. DOI: 10.1016/j.ijmecsci.2017.05.047.

[21] G. Dal Maso, G. Orlando, R. Toader. 'Fracture models for elasto-plastic materials as limits of gradient damage models coupled with plasticity: the antiplane case'. In: *Calculus of Variations and Partial Differential Equations* 55(3) (2016), pp. 1–39. DOI: 10.1007/s00526-016-0981-z.

[22] C. Kuhn, T. Noll, R. Müller. 'On phase field modeling of ductile fracture'. In: *GAMM-Mitteilungen* 39.1 (2016), pp. 35–54. DOI: 10.1007/s00466-015-1151-4.

[23] M. Ambati, T. Gerasimov, L. De Lorenzis. 'Phase-field modeling of ductile fracture'. In: *Computational Mechanics* 55(5) (2015), pp. 1017–1040. DOI: 10.1007/s00466-015-1151-4.

[24] M. Ambati, R. Kruse, L. De Lorenzis. 'A phase-field model for ductile fracture at finite strains and its experimental verification'. In: *Computational Mechanics* 57(1) (2016), pp. 149–167. DOI: 10.1007/s00466-015-1225-3.
[25] J. Ulloa, P. Rodríguez, E. Samaniego. 'On the modeling of dissipative mechanisms in a ductile softening bar'. In: *Journal of Mechanics of Materials and Structures* 11.4 (2016), pp. 463-490.
DOI: 10.2140/jomms.2016.11.463.

[26] C. Miehe, M. Hofacker, L.-M. Schänzel, F. Aldakheel. 'Phase field modeling of fracture in multi-physics problems. Part II. Coupled brittle-to-ductile failure criteria and crack propagation in thermo-elastic-plastic solids'. In: *Computer Methods in Applied Mechanics and Engineering* 294 (2015), pp. 486-522.
DOI: 10.1016/j.cma.2014.11.017.

[27] M. J. Borden, T. J. R. Hughes, C. M. Landis, A. Anvari, I. J. Lee. 'A phase-field formulation for fracture in ductile materials: Finite deformation balance law derivation, plastic degradation, and stress triaxiality effects'. In: *Computer Methods in Applied Mechanics and Engineering*. Phase Field Approaches to Fracture 312 (2016), pp. 130-166.
DOI: 10.1016/j.cma.2016.09.005.

[28] F. Aldakheel, P. Wriggers, C. Miehe. 'A modified Gurson-type plasticity model at finite strains: formulation, numerical analysis and phase-field coupling'. In: *Computational Mechanics* 62 (2018), pp. 815-833.
DOI: 10.1007/s00466-017-1530-0.

[29] M. Dittmann, F. Aldakheel, J. Schulte, F. Schmidt, M. Krüger, P. Wriggers, C. Hesch. 'Phase-field modeling of porous-ductile fracture in non-linear thermo-elasto-plastic solids'. In: *Computer Methods in Applied Mechanics and Engineering* 361 (2020), p. 112730.
DOI: 10.1016/j.cma.2019.112730.

[30] J. Choo, W. Sun. 'Coupled phase-field and plasticity modeling of geological materials: From brittle fracture to ductile flow'. In: *Computer Methods in Applied Mechanics and Engineering* 330 (2018), pp. 1-32.
DOI: 10.1016/j.cma.2017.10.009.

[31] R. Alessi, M. Ambati, T. Gerasimov, S. Vidoli, L. De Lorenzis. 'Comparison of Phase-Field Models of Fracture Coupled with Plasticity'. In: *Advances in Computational Plasticity: A Book in Honour of D. Roger, J. Owen*. Ed. by E. Oñate, D. Peric, E.A. De Souza Neto, M. Chiumenti. Vol. 46. Computational Methods in Applied Sciences. Cham: Springer, 2018, pp. 1-21.

[32] Al. Marengo, U. Perego. 'A concise review of small-strain phase-field modeling of ductile fracture'. In: *European Journal of Mechanics - A/Solids* 101 (2023), p. 105058.
DOI: 10.1016/j.euromechsol.2023.105058.

[33] L. De Lorenzis, C. Maurini. 'Nucleation under multi-axial loading in variational phase-field models of brittle fracture'. In: *International Journal of Fracture* 237 (2022), pp. 61-81.
DOI: 10.1007/s10704-021-00555-6.

[34] George E Dieter. *Mechanical metallurgy*. SI metric. McGraw-Hill series in materials science and engineering. London: McGraw-Hill, 1988.

[35] K.H. Roscoe, J.B. Burland. 'On the Generalised Stress-strain Behaviour of 'Wet' Clay.' In: *Engineering Plasticity: Papers for a Conference held in Cambridge March 1968*. Ed. by J. Heyman, F.A. Leckie. Cambridge: Cambridge University Press, 1968, pp. 535-609.

[36] G. Dal Maso, A. Desimone. 'Quasistatic ecolution for Cam-Clay plasticity: Examples of spatially homogeneous solutions'. In: *Mathematical Models and Methods in Applied Sciences* 19.09 (2009), pp. 1643-1711.
DOI: 10.1142/S0218202509003942.

[37] G. Bacquaert, S. Raude, V. Alves-Fernandes, F. Voldoire, D. Kondo. 'A standard thermodynamic-based extension of the Modified Cam-Clay soil model and its applications'. In: *European Journal of Mechanics - A/Solids* 103 (2024), p. 105122.
DOI: 10.1016/j.euromechsol.2023.105122.

[38] R. I. Borja, S. R. Lee. 'Cam-Clay plasticity, Part 1: Implicit integration of elasto-plastic constitutive relations'. In: *Computer Methods in Applied Mechanics and Engineering* 78.1 (1990), pp. 49-72.
DOI: 10.1016/0045-7825(90)90152-C.

[39] R.A. Saada. 'Effect of Some Parameters on the Elastoplastic Behavior of Green Sand'. In: *IUTAM Symposium on Multiscale Modeling and Characterization of Elastic-Inelastic Behavior of Engineering Materials*. Ed. by S. Ahzi, M. Cherkaoui, M. A. Khaleel, H. M. Zbib, M. A. Zikry, B. Lamatina. Vol. 114. Dordrecht: Springer Netherlands, 2004, pp. 83-88.
DOI: 10.1007/978-94-017-0483-0-11.
[40] N. Zouain, I. Pontes, J. Vaunat. 'Potentials for the modified Cam-Clay model'. In: *European Journal of Mechanics - A/Solids* 29.3 (2010), pp. 327-336. DOI: 10.1016/j.euromechsol.2009.11.008.

[41] Y. Surkutwar, C. Sandu, C. Untaroiu. 'Review of modeling methods of compressed snow-tire interaction'. In: *Journal of Terramechanics* 105 (2023), pp. 27-40. DOI: 10.1016/j.jterra.2022.10.004.

[42] E.A. De Souza Neto, D. Perić, D.R.J. Owen. *Computational methods for plasticity : theory and applications*. Chichester: Wiley, 2008.

[43] M Ortiz, A Pandolfi. 'A variational Cam-clay theory of plasticity'. In: *Computer Methods in Applied Mechanics and Engineering* 193.27 (2004), pp. 2645-2666. DOI: 10.1016/j.cma.2003.08.008.

[44] K. C. Bennett. 'An energy approach to Modified Cam-Clay plasticity and damage modeling of cohesive soils'. In: *Acta Geotechnica* 15.1 (2020), pp. 165-177. DOI: 10.1007/s11440-019-00880-0.

[45] V. Tvergaard. 'Influence of voids on shear band instabilities'. In: *International Journal of Fracture* 17.4 (1981), pp. 389-407. DOI: 10.1007/BF00036191.

[46] V. Tvergaard. 'On localization in ductile materials containing spherical voids'. In: *International Journal of Fracture* 18.4 (1982), pp. 237-252. DOI: 10.1007/BF00015686.

[47] C. Chu, A. Needleman. 'Void Nucleation Effects in Biaxially Stretched Sheets'. In: *Journal of Engineering Materials and Technology* 102.3 (1980), pp. 249-256. DOI: 10.1115/1.3224807.

[48] B. Halphen, Q.S. Nguyen. 'Sur Les Matériaux Standards Généralisés'. In: *Journal de Mécanique* 14.1 (1975), pp. 39-63.

[49] P. Steinmann, K. Runesson. *The Catalogue of Computational Material Models: Basic Geometrically Linear Models in 1D*. Cham: Springer, 2021.

[50] Han Yang, Hexiang Wang, Yuan Feng, Boris Jeremić. 'Plastic-Energy Dissipation in Pressure-Dependent Materials'. In: *Journal of Engineering Mechanics* 146.3 (2020), p. 04019139. DOI: 10.1061/(ASCE)EM.1943-7889.0001728.

[51] Martin Kružík, Tomáš Roubíček. *Mathematical Methods in Continuum Mechanics of Solids*. Interaction of Mechanics and Mathematics. Cham: Springer International Publishing, 2019.

[52] N. Aravas. 'On the numerical integration of a class of pressure-dependent plasticity'. In: *International Journal for Numerical Methods in Engineering* 24 (1987), pp. 1395-1416. DOI: 10.1002/nme.1620240713.

[53] K. Pham, H. Amor, J.-J. Marigo, C. Maurini. 'Gradient Damage Models and Their Use to Approximate Brittle Fracture'. In: *International Journal of Damage Mechanics* 20.4 (2011), pp. 618-652. DOI: 10.1177/1056789510386852.

[54] L. Ambrosio, V. M. Tortorelli. 'On the approximation of free discontinuity problems'. In: *Boll. Unione Matematica Italiana* 6-B(7) (1992), pp. 105-123.

[55] Dassault Systèmes Simulia Corp. *ABAQUS/Standard Documentation, Version 2022*. Dassault Systèmes, 2022.

[56] E. Azinpour, J. P. S. Ferreira, M. P. L. Parente, J. Cesar de Sa. 'A simple and unified implementation of phase field and gradient damage models'. In: *Advanced Modeling and Simulation in Engineering Sciences* 5.1 (2018). DOI: 10.1186/s40323-018-0106-7.

[57] Y. Navidtehrani, C. Betegón, E. Martínez-Pañeda. 'A simple and robust Abaqus implementation of the phase field fracture method'. In: *Applications in Engineering Science* 6 (2021), p. 100050. DOI: 10.1016/j.apples.2021.100050.

[58] T. Gerasimov, L. De Lorenzis. 'On penalization in variational phase-field models of brittle fracture'. In: *Computer Methods in Applied Mechanics and Engineering* 354 (2019), pp. 990-1026. DOI: 10.1016/j.cma.2019.05.038.
[59] E. Tanné, T. Li, B. Bourdin, J.-J. Marigo, C. Maurini. ‘Crack nucleation in variational phase-field models of brittle fracture’. In: *Journal of the Mechanics and Physics of Solids* 110 (2018), pp. 80–99. doi: 10.1016/j.jmps.2017.09.006.

[60] E. Tanné. ‘Variational phase-field models from brittle to ductile fracture: nucleation and propagation’. Materials and structures in mechanics. PhD thesis. Université Paris Saclay, 2017.

[61] Gianni Dal Maso, Antonio DeSimone, Maria Giovanna Mora. ‘Quasistatic Evolution Problems for Linearly Elastic-Perfectly Plastic Materials’. In: *Archive for Rational Mechanics and Analysis* 180.2 (2006), pp. 237–291. doi: 10.1007/s00205-005-0407-0.

[62] P.W. Bridgman. *Studies in large plastic flow and fracture: with special emphasis on the effects of hydrostatic pressure*. Metallurgy and metallurgical engineering series. New York: McGraw-Hill, 1952.

[63] S. Boyd, L. Vandenberghe. *Convex optimization*. Cambridge: Cambridge University Press, 2004.

[64] J.-J. Marigo, C. Maurini, K. Pham. ‘An overview of the modelling of fracture by gradient damage models’. In: *Meccanica* 51.12 (2016), pp. 3107–3128. doi: 10.1007/s11012-016-0538-4.
# Appendix

## A. Dual norm for the Cam-Clay criterion

We show as follows how to obtain the expression (6) for the dual strain norm $||\dot{\boldsymbol{\epsilon}}^p||_M^*$ of the Cam-Clay stress norm $||\boldsymbol{\sigma}||_M$. The following minimisation problem defines the dual norm

$$ ||\dot{\boldsymbol{\epsilon}}^p||_M^* := \sup_{||\boldsymbol{\sigma}||_M \le \sigma_y} \frac{\boldsymbol{\sigma}}{\sigma_y(p)} \cdot \dot{\boldsymbol{\epsilon}}^p = \sup_{||\tilde{\boldsymbol{\sigma}}||_M \le 1} \tilde{\boldsymbol{\sigma}} \cdot \dot{\boldsymbol{\epsilon}}^p \quad \text{with} \quad ||\tilde{\boldsymbol{\sigma}}||_M := \sqrt{\frac{3}{2} ||\tilde{\boldsymbol{\sigma}}_{dev}||^2 + M^2 \tilde{\sigma}_H^2} \quad (49) $$

and $\tilde{\boldsymbol{\sigma}} := \boldsymbol{\sigma}/\sigma_y$. Using the Lagrange multiplier method, the solution of the optimisation problem (7) including the convex constraint $||\tilde{\boldsymbol{\sigma}}||_M \le 1$ must be a stationary point of the Lagrangian

$$ \mathcal{L}(\tilde{\boldsymbol{\sigma}}_{dev}, \sigma_H, \tilde{\lambda}) = \tilde{\boldsymbol{\sigma}}_{dev} \cdot \dot{\boldsymbol{\epsilon}}_{dev}^p + \tilde{\sigma}_H \text{tr}(\dot{\boldsymbol{\epsilon}}^p) + \tilde{\lambda} \left[ 1 - \left[ \frac{3}{2} ||\tilde{\boldsymbol{\sigma}}_{dev}||^2 + M^2 \tilde{\sigma}_H^2 \right] \right] \quad (50) $$

for any deviatoric tensor $\tilde{\boldsymbol{\sigma}}_{dev}$, scalar $\tilde{\sigma}_H$, and Lagrange multiplier $\tilde{\lambda} \ge 0$ [63]. The first-order optimality conditions with respect to $\tilde{\boldsymbol{\sigma}}_{dev}$ and $\tilde{\sigma}_H$ give

$$ \dot{\boldsymbol{\epsilon}}_{dev}^p - 3\tilde{\lambda}\tilde{\boldsymbol{\sigma}}_{dev} = \mathbf{0}, \quad \text{tr}(\dot{\boldsymbol{\epsilon}}^p) - 2M^2\tilde{\lambda}\tilde{\sigma}_H = 0, \quad (51) $$

which imply in particular that $\tilde{\boldsymbol{\sigma}}_{dev}$ is co-axial to $\dot{\boldsymbol{\epsilon}}_{dev}^p$.

The unilateral optimality condition with respect to the non-negative Lagrange multiplier is in the form of a set of Karush-Kuhn-Tucker complementarity conditions [63]

$$ (i) \quad \tilde{\lambda} \ge 0, \quad (ii) \quad 1 - \left[ \frac{3}{2} ||\tilde{\boldsymbol{\sigma}}_{dev}||^2 + M^2 \tilde{\sigma}_H^2 \right] \ge 0, \quad (iii) \quad \left[ 1 - \left[ \frac{3}{2} ||\tilde{\boldsymbol{\sigma}}_{dev}||^2 + M^2 \tilde{\sigma}_H^2 \right] \right] \tilde{\lambda} = 0. \quad (52) $$

The case $\tilde{\lambda} = 0$ corresponds to the elastic regime in which $\dot{\boldsymbol{\epsilon}}_{dev}^p = \mathbf{0}$ and $\text{tr}(\dot{\boldsymbol{\epsilon}}^p) = 0$. In the plastic regime $\tilde{\lambda} > 0$ and the unilateral constraint (52ii) must be satisfied as an equality. Solving the optimality conditions (51) for the stress and replacing $\tilde{\boldsymbol{\sigma}}_{dev}$, $\tilde{\sigma}_H$ in the equality constraint gives the solution for the Lagrange multiplier

$$ \tilde{\lambda} = \frac{1}{2} \sqrt{\frac{2}{3} ||\dot{\boldsymbol{\epsilon}}_{dev}^p||^2 + \frac{\text{tr}^2(\dot{\boldsymbol{\epsilon}}^p)}{M^2}}. $$

Finally, replacing the optimal stresses from (51) in the expression for the dual norm in (49) gives the result

$$ ||\dot{\boldsymbol{\epsilon}}^p||_M^* = \sqrt{\frac{2}{3} ||\dot{\boldsymbol{\epsilon}}_{dev}^p||^2 + \frac{\text{tr}^2(\dot{\boldsymbol{\epsilon}}^p)}{M^2}}. $$

One can similarly prove that

$$ ||\tilde{\boldsymbol{\sigma}}||_M^{**} := \sup_{||\dot{\boldsymbol{\epsilon}}^p||_M^* \le 1} \tilde{\boldsymbol{\sigma}} \cdot \dot{\boldsymbol{\epsilon}}^p = ||\tilde{\boldsymbol{\sigma}}||_M. \quad (53) $$
## B. Variational derivation of the strong form plasticity equations and return mapping algorithm

We derive here the strong optimality conditions for the minimisation problem (12) of the elasto-plastic energy functional (13). Subsequently, we formulate a return mapping algorithm for time integration. The extension to the case including the phase-field damage variable is classical and it is not reproduced here. It can be found for example in [64] for the case without plasticity and in [18, 19] for the case with coupled damage and plasticity.

### A. Strong form equations

The solution $(u, \boldsymbol{\epsilon}^p) \in \mathcal{V}_i \times \mathcal{P}$ of the problem (12) must verify the condition of unilateral local directional minimality, i.e. for all test functions $(v, \eta) \in H_0^1(\Omega, \mathbb{R}^3) \times L^2(\Omega, \mathbb{R}^{3,3})$, there exists $\bar{h} > 0$ such that for all $h \in (0, \bar{h})$

$$ \mathcal{E}_i(u + h\mathbf{v}, \boldsymbol{\epsilon}^p + h\eta) - \mathcal{E}_i(u, \boldsymbol{\epsilon}^p) \geq 0, \quad (54) $$

where $\mathcal{V}_i$ and $\mathcal{P}$ are the space of admissible displacements at the $i$-th time step and the space of admissible plastic strain increments, respectively.

Taking $\eta = 0$ in (54) gives the standard variational formulation of the elastic sub-problem, which delivers the equilibrium equations for the stress tensor $\sigma$ along with the associated Neumann boundary conditions, see (15). The formal proof is classical and it is not reproduced here.

The case with $v = 0$ and $\eta \neq 0$ is less standard, because the presence of the norm function renders the energy non-smooth with respect to $\boldsymbol{\epsilon}^p$. The energy increment in this case is given by

$$ \begin{aligned} \mathcal{E}_i(\mathbf{u}, \boldsymbol{\epsilon}^p + h\eta) - \mathcal{E}_i(\mathbf{u}, \boldsymbol{\epsilon}^p) &= \int_{\Omega} h \frac{\partial \psi_e}{\partial \boldsymbol{\epsilon}^p} (\boldsymbol{\epsilon}(\mathbf{u}) - \boldsymbol{\epsilon}^p) \cdot \eta \, d\mathbf{x} \\ &\quad + \int_{\Omega} [\psi_p (p_{i-1} + \|\Delta \boldsymbol{\epsilon}^p + h \eta\|_{\mathrm{M}}^*) - \psi_p (p_{i-1} + \|\Delta \boldsymbol{\epsilon}^p\|_{\mathrm{M}}^*)] \, d\mathbf{x} + o(h) \\ &= \int_{\Omega} -h \sigma(\boldsymbol{\epsilon}(\mathbf{u}) - \boldsymbol{\epsilon}^p) \cdot \eta \, d\mathbf{x} \\ &\quad + \int_{\Omega} [\psi_p (p_{i-1} + \|\Delta \boldsymbol{\epsilon}^p + h \eta\|_{\mathrm{M}}^*) - \psi_p (p_{i-1} + \|\Delta \boldsymbol{\epsilon}^p\|_{\mathrm{M}}^*)] \, d\mathbf{x} + o(h) \end{aligned} $$

where $\Delta \boldsymbol{\epsilon}^p = \boldsymbol{\epsilon}^p - \boldsymbol{\epsilon}_{i-1}^p$ is the plastic strain increment. In the expression above, we expand the derivative of the first term exploiting the smoothness of the elastic energy density with respect to the elastic strain and that

$$ \sigma = - \frac{\partial \psi_e}{\partial \boldsymbol{\epsilon}_p} (\boldsymbol{\epsilon} - \boldsymbol{\epsilon}^p) = K \operatorname{tr}(\boldsymbol{\epsilon} - \boldsymbol{\epsilon}^p) \mathbf{I}_3 + 2\mu [\boldsymbol{\epsilon}_{dev} - \boldsymbol{\epsilon}_{dev}^p]. $$

For the second term, we should separately consider the case of solution with $\|\Delta \boldsymbol{\epsilon}^p\|_{\mathrm{M}}^*$ equal to zero or not, as follows.

• If $\|\Delta \boldsymbol{\epsilon}^p\|_{\mathrm{M}}^* = 0$,

$$ \begin{aligned} \mathcal{E}_i(\mathbf{u}, \boldsymbol{\epsilon}^p + h\eta) - \mathcal{E}_i(\mathbf{u}, \boldsymbol{\epsilon}^p) &= h \int_{\Omega} [-\boldsymbol{\sigma} \cdot \eta + \psi'_p(p_{i-1}) \|\eta\|_{\mathrm{M}}^*] \, d\mathbf{x} + o(h) \\ &= h \int_{\Omega} [-\boldsymbol{\sigma} \cdot \eta + \sigma_y(p_{i-1}) \|\eta\|_{\mathrm{M}}^*] \, d\mathbf{x} + o(h) \end{aligned} $$

where we used $\psi'_p(p_{i-1}) = \sigma_y(p_{i-1}) = \sigma_0 + H p_{i-1}$ (note also that in this case $p = p_{i-1}$). Taking $\eta(\mathbf{x}) = \tilde{\eta}\eta(\mathbf{x})$, where $\eta \in L^2(\Omega)$ denotes an arbitrary scalar function and $\tilde{\eta} \in \mathbb{R}^{3,3}$ an arbitrary constant, yields the pointwise necessary optimality condition

$$ \sigma \cdot \tilde{\eta} \leq \sigma_y(p) \|\tilde{\eta}\|_{\mathrm{M}}^* \quad \text{in } \Omega $$
for all $\tilde{\eta} \in \mathbb{R}^{3,3}$. In particular

$$ \|\boldsymbol{\sigma}\|_{\mathrm{M}}^{**} := \sup_{\|\tilde{\eta}\|_{\mathrm{M}}^{*} \le 1} \boldsymbol{\sigma} \cdot \tilde{\eta} \le \sigma_{y}(p), $$

where $\|\boldsymbol{\sigma}\|_{\mathrm{M}}^{**} = \|\boldsymbol{\sigma}\|_{\mathrm{M}}$ is the dual of the dual Cam-Clay norm, coinciding with the Cam-Clay norm itself, see (53). Hence the optimality condition implies that

$$ f_y(\boldsymbol{\sigma}, p) = \|\boldsymbol{\sigma}\|_{\mathrm{M}} - \sigma_y(p) \le 0, $$

where $\|\boldsymbol{\sigma}\|_{\mathrm{M}} = \sigma_{eq,cc}(\boldsymbol{\sigma})$ according to (6).

* If $\|\Delta\boldsymbol{\varepsilon}^p\|_{\mathrm{M}}^* > 0$, the energy increment is given by

$$ \begin{aligned} \mathcal{E}_i(\mathbf{u}, \boldsymbol{\varepsilon}^p + h\boldsymbol{\eta}) - \mathcal{E}_i(\mathbf{u}, \boldsymbol{\varepsilon}^p) &= \int_{\Omega} [-h \boldsymbol{\sigma} \cdot \boldsymbol{\eta} + \psi'_p(p) (\|\Delta\boldsymbol{\varepsilon}^p + h\boldsymbol{\eta}\|_{\mathrm{M}}^* - \|\Delta\boldsymbol{\varepsilon}^p\|_{\mathrm{M}}^*)] d\mathbf{x} + o(h) \\ &= h \int_{\Omega} [-\boldsymbol{\sigma} \cdot \boldsymbol{\eta} + \sigma_y(p) \frac{\partial \|\Delta\boldsymbol{\varepsilon}^p\|_{\mathrm{M}}^*}{\partial \Delta\boldsymbol{\varepsilon}^p} \cdot \boldsymbol{\eta}] d\mathbf{x} + o(h) \end{aligned} $$

and the arbitrariness of $\boldsymbol{\eta}$ implies the necessary optimality condition

$$ \boldsymbol{\sigma} = \sigma_y(p) \frac{\partial \|\Delta\boldsymbol{\varepsilon}^p\|_{\mathrm{M}}^*}{\partial \Delta\boldsymbol{\varepsilon}^p}. \qquad (55) $$

Since the derivative of the norm function reads

$$ \frac{\partial \|\Delta \boldsymbol{\varepsilon}^p\|_{\mathrm{M}}^*}{\partial \Delta \boldsymbol{\varepsilon}^p} = \frac{\partial}{\partial \Delta \boldsymbol{\varepsilon}^p} \sqrt{\frac{2}{3} \|\Delta \boldsymbol{\varepsilon}_{dev}^p\|^2 + \frac{\mathrm{tr}(\Delta \boldsymbol{\varepsilon}^p)^2}{M^2}} = \frac{1}{\|\Delta \boldsymbol{\varepsilon}^p\|_{\mathrm{M}}^*} \left[ \frac{2}{3} \Delta \boldsymbol{\varepsilon}_{dev}^p + \frac{\mathrm{tr}(\Delta \boldsymbol{\varepsilon}^p)}{M^2} \mathbf{I}_3 \right], $$

from (55) we infer that $\boldsymbol{\sigma}$ and $\Delta\boldsymbol{\varepsilon}^p$ are co-axial. Also, by computing the Cam-Clay norm of $\boldsymbol{\sigma}$ in (55) and using the previous expression for the derivative of the norm function, it is immediate to show that the Cam-Clay criterion is verified as an equality, $\|\boldsymbol{\sigma}\|_{\mathrm{M}} = \sigma_y(p)$. Moreover, the derivative

$$ \frac{\partial f_y(\boldsymbol{\sigma}, p)}{\partial \boldsymbol{\sigma}} = \frac{\partial \|\boldsymbol{\sigma}\|_{\mathrm{M}}}{\partial \boldsymbol{\sigma}} = \frac{\partial}{\partial \boldsymbol{\sigma}} \sqrt{\frac{3}{2} \|\boldsymbol{\sigma}_{dev}\|^2 + \frac{M^2 \mathrm{tr}^2(\boldsymbol{\sigma})}{9}} = \frac{1}{\|\boldsymbol{\sigma}\|_{\mathrm{M}}} \left[ \frac{3}{2} \boldsymbol{\sigma}_{dev} + \frac{M^2 \mathrm{tr}(\boldsymbol{\sigma})}{9} \mathbf{I}_3 \right] $$

is also co-axial to $\boldsymbol{\sigma}$, hence to $\Delta\boldsymbol{\varepsilon}^p$, and it is straightforward to show that

$$ \Delta\boldsymbol{\varepsilon}^p = \Delta\lambda \frac{\partial f_y}{\partial \boldsymbol{\sigma}} \qquad (56) $$

with $\Delta\lambda := \|\Delta\boldsymbol{\varepsilon}^p\|_{\mathrm{M}}^*$, which is the incremental form of the classical normality rule of elasto-plasticity. For later convenience, let us define the director of the co-axial tensors as $\mathbf{n} := \boldsymbol{\sigma}_{dev}/\|\boldsymbol{\sigma}_{dev}\|$, so that it is also

$$ \Delta\boldsymbol{\varepsilon}_{dev}^p = \|\Delta\boldsymbol{\varepsilon}_{dev}^p\| \mathbf{n} \qquad (57) $$

In conclusion, the necessary condition for minimality of the total energy with respect to the plastic strain implies that

$$ f_y(\boldsymbol{\sigma}, p) = \|\boldsymbol{\sigma}\|_{\mathrm{M}} - \sigma_y(p) \le 0 $$

and

$$ \left\{ \begin{array}{ll} \Delta\boldsymbol{\varepsilon}^p = \mathbf{0}, & \text{if } \|\boldsymbol{\sigma}\|_{\mathrm{M}} < \sigma_y(p); \\ \Delta\boldsymbol{\varepsilon}^p \text{ is given by (56),} & \text{if } \|\boldsymbol{\sigma}\|_{\mathrm{M}} = \sigma_y(p). \end{array} \right. \qquad (58) $$

The conditions above can be also rewritten as follows

$$ \Delta\boldsymbol{\varepsilon}^p = \Delta\lambda \frac{\partial f_y}{\partial \boldsymbol{\sigma}}, \quad \Delta\lambda \ge 0, \quad f_y(\boldsymbol{\sigma}, p) \le 0, \quad f_y(\boldsymbol{\sigma}, p) \Delta\lambda = 0, \qquad (59) $$

which is the classical form of the Kuhn-Tucker conditions of elasto-plasticity.
## B. Return mapping algorithm

For integration of the incremental constitutive equations (59), it is useful to introduce the elastic predictor stress

$$ \sigma^{tr} := K \operatorname{tr}(\boldsymbol{\epsilon} - \boldsymbol{\epsilon}_{i-1}^p) \boldsymbol{I}_3 + 2 \mu [\boldsymbol{\epsilon}_{dev} - \boldsymbol{\epsilon}_{dev,i-1}^p] = \boldsymbol{\sigma} + K \operatorname{tr}(\Delta \boldsymbol{\epsilon}^p) \boldsymbol{I}_3 + 2 \mu \Delta \boldsymbol{\epsilon}_{dev}^p \quad (60) $$

which can be computed based on the current strain and on the quantities at the previous time step. The deviatoric component of the above equation yields

$$ \sigma_{dev}^{tr} = \sigma_{dev} + 2\mu\Delta\boldsymbol{\epsilon}_{dev}^{p}, $$

which shows that, since $\sigma_{dev}$ and $\Delta\boldsymbol{\epsilon}_{dev}^p$ are co-axial, also $\sigma_{dev}^{tr}$ must be, i.e.

$$ \mathbf{n}^{tr} := \sigma_{dev}^{tr} / \|\sigma_{dev}^{tr}\| = \mathbf{n}. \quad (61) $$

Using (60), the optimality condition (55) can be written as

$$ \boldsymbol{\sigma} = \boldsymbol{\sigma}^{tr} - K \operatorname{tr}(\Delta \boldsymbol{\epsilon}^p) \boldsymbol{I}_3 - 2 \mu \Delta \boldsymbol{\epsilon}_{dev}^p = \left[ \frac{\sigma_y(p_{i-1})}{\|\Delta \boldsymbol{\epsilon}^p\|_{\mathrm{M}}^*} + H \right] \left[ \frac{2}{3} \Delta \boldsymbol{\epsilon}_{dev}^p + \frac{\operatorname{tr}(\Delta \boldsymbol{\epsilon}^p)}{M^2} \boldsymbol{I}_3 \right], \quad (62) $$

where we also used

$$ \sigma_y(p) = \sigma_y(p_{i-1}) + H \Delta p = \sigma_y(p_{i-1}) + H \|\Delta \boldsymbol{\epsilon}^p\|_{\mathrm{M}}^* $$

stemming from (4). Taking the trace and the deviatoric part and recalling (57) and (61) results in

$$ \left[ \frac{1}{M^2} \left[ \frac{\sigma_y(p_{i-1})}{\|\Delta \boldsymbol{\epsilon}^p\|_{\mathrm{M}}^*} + H \right] + K \right] \operatorname{tr}(\Delta \boldsymbol{\epsilon}^p) = \frac{1}{3} \operatorname{tr}(\boldsymbol{\sigma}^{tr}), \quad (63) $$

$$ 2 \left[ \frac{1}{3} \left[ \frac{\sigma_y(p_{i-1})}{\|\Delta \boldsymbol{\epsilon}^p\|_{\mathrm{M}}^*} + H \right] + \mu \right] \|\Delta \boldsymbol{\epsilon}_{dev}^p\| = \|\boldsymbol{\sigma}_{dev}^{tr}\|, \quad (64) $$

which is a system of two scalar nonlinear equations for $\operatorname{tr}(\Delta \boldsymbol{\epsilon}^p)$ and $\|\Delta \boldsymbol{\epsilon}_{dev}^p\|$ given the elastic predictor.