# Macroscopic yield criteria for ductile materials containing spheroidal voids: An Eshelby-like velocity fields approach

Vincent Monchiet, Eric Charkaluk, Djimedo Kondo

► To cite this version:

Vincent Monchiet, Eric Charkaluk, Djimedo Kondo. Macroscopic yield criteria for ductile materials containing spheroidal voids: An Eshelby-like velocity fields approach. Mechanics of Materials, 2014, 72, pp.1-18. 10.1016/j.mechmat.2013.05.006. hal-00955545

HAL Id: hal-00955545

https://hal.science/hal-00955545

Submitted on 26 Jul 2021

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Macroscopic yield criteria for ductile materials containing spheroidal voids: An Eshelby-like velocity fields approach

V. Monchiet<sup>a,*</sup>, E. Charkaluk<sup>b</sup>, D. Kondo<sup>c</sup>

Making use of limit analysis theory, we derive a new expression of the macroscopic yield function for a rigid ideal-plastic von Mises matrix containing spheroidal cavities (oblate or prolate). Key in the development of the new criterion is the consideration of Eshelby-like velocity fields which are built by taking advantage of the solution of the equivalent inclusion problem in which the eigenstrains rate are unknown for the plasticity problem. These heterogeneous trial velocity fields contain non-axisymmetric components which prove to be original in the context of limit analysis of hollow spheroid. After carefully computing the macroscopic plastic dissipation and implementing a minimization procedure required by the use of the Eshelby-like velocity fields, we derive, for the porous medium, a two-field estimate of the anisotropic yield criterion whose closed-form expression is provided. This estimate is compared to existing criteria based on limit analysis theory. Interestingly, in contrast to these criteria, the new results predict a significant effect of shear loadings in the particular case of ductile materials weakened by penny-shaped cracks.

# 1. Introduction

Since the pioneering works of McClintock (1968), Rice and Tracey (1969) and Gurson (1977), several studies have been performed in order to incorporate void shape in ductile damage analysis. To this end, Lee and Mear (1992) carried out numerical investigations of cavities shape effects on damage growth. In their analysis, a spheroidal cavity is embedded in an infinite plastic or viscoplastic material (obeying to the classical Norton law) subjected to an axisymmetric remote stress field. This study, in fact, extends the works of Rice and Tracey (1969) and Budiansky et al. (1982) to the case of spheroidal voids for which a class of axisymmetric trial velocity fields has been considered. This

family of velocity fields contains the exact solution for a spheroidal void embedded in an infinite linear viscous material subjected to an axisymmetric remote stress field. Moreover, it also includes the trial fields already used by Gurson (1977) in order to obtain coupled model in the case of spherical or cylindrical cavities. Later, Tvergaard (1981) observed that Gurson model is too stiff when compared with finite element unit-cell computations. This observation motivated the heuristic extension of the Gurson model which has been proposed by Tvergaard and Needleman (1984) (see also Tvergaard, 1981, 1990). Such extension, known as the GTN model, introduces three parameters, $q_1$, $q_2$ and $q_3$, which have to be determined. It is widely used in structural computations.

Later, Gologanu et al. (1993, 1994, 1997) and Garajeu et al. (2000) (see also Gârâjeu, 1995) have extended the Gurson approach by accounting for the cavities shape effects. They performed the limit analysis approach on a spheroidal unit cell containing a confocal spheroidal cavity.

<footnote><sup>a</sup> Université Paris-Est, Modélisation et Simulation Multi Echelle, MSME UMR 8208 CNRS, 5 bd Descartes, 77454 Marne-la-Vallée, France</footnote>

<footnote><sup>b</sup> Université de Lille I, Laboratoire de Mécanique de Lille, LML UMR 8107 CNRS, Cité Scientifique, Bd. Paul Langevin, 59655 Villeneuve d'Ascq, France <sup>c</sup></footnote>

<footnote>Université Pierre et Marie Curie (UPMC), Institut D'Alembert, UMR 7190 CNRS, 4, Place Jussieu, 75252 Paris Cedex 05, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: vincent.monchiet@univ-mlv.fr (V. Monchiet), eric.charkaluk@univ-lille1.fr (E. Charkaluk), djimedo.kondo@upmc.fr (D. Kondo).</footnote>
A two-field trial velocity, derived from Lee and Mear (1992) and complying with homogeneous strain rate boundary conditions, is used to obtain an anisotropic macroscopic yield function of the porous medium. It is worth noting that these criteria were derived assuming macroscopic axisymmetric loadings, this has motivated an heuristical extension in Gologanu et al. (1997) to arbitrary loadings.

Also in the framework of limit analysis, the effects of plastic anisotropy of the matrix have been also addressed in Benzerga and Besson (2001), Monchiet et al. (2008) and Keralavarma and Benzerga (2010), etc. In addition to these works, mention can be made of very recent studies devoted to porous materials with incompressible matrix exhibiting asymmetry between tension and compression (Cazacu and Stewart, 2009), or with matrix displaying pressure sensitive behavior (see for instance Guo et al., 2008; Shen et al., 2012, etc.) For a recent review on the topic of ductile fracture, the reader is referred to Benzerga and Leblond (2010).

Important contributions to the study of porous non-linear materials were also obtained in the framework of the non linear homogenization methods introduced by Ponte-Castañeda (1991) and Ponte-Castañeda and Zaidman (1994) (see also Ponte-Castañeda and Suquet, 1998). By assuming that the spatial distribution function of the voids is identical to that of their shape, these authors derived, in particular, Hashin-Shtrikman type upper bounds for ductile media containing ellipsoidal cavities. These approach, applied to a representative elementary volume for the porous medium, used the concept of a linear comparison material and provided an elliptic criterion, the so-called Hashin-Shtrikman non linear bound, which takes properly into account the void induced orthotropy.<fnref n="1" /> Later, Kailasam et al. (1997) extended the analysis of Ponte-Castañeda and Zaidman (1994) by considering a spatial distribution function of cavities which differs from the one characterizing the cavities shape. For completeness, it is convenient to note that Danas et al. (2008) recently derived accurate yield surfaces for the porous materials by using a second order homogenization method. The reader interested by the basic aspects of micromechanics of heterogeneous materials may refer for instance to the textbook of Nemat Nasser and Hori (1999) (see also Dormieux et al., 2006 for micromechanics of porous media).

Still in the context of the limit analysis approach, we propose to investigate the use of more refined velocity fields in order to improve the predictions of the Gurson-type models accounting for voids shape effects. Note that this type of study has been performed in the case of spherical voids by Gārājeu (1995) and Garajeu et al. (2000) who considered the exact solution of the elastic hollow sphere subjected to an arbitrary loading. More recently, still in the case of spherical cavities, Monchiet et al. (2011) considered Eshelby-based trial velocity field and successfully derived a new macroscopic criterion.

The main objective of the present paper is to present the limit analysis approach based on the Eshelby-like velocity

to the case of non spherical (oblate or prolate) voids. The paper is organized as follows. In Section 2, we first present the basic principle of the limit analysis of a spheroidal unit cell containing a confocal spheroidal void. Section 3.1 is devoted to a description of the new trial velocity fields inspired of the “exterior point” solution of the Eshelby equivalent inclusion problem (Eshelby, 1959) which, in the case of a spheroidal cavity, contains non-axisymmetric components. The conditions related to the chosen trial velocity fields are analyzed in Section 3.2. Section 4 is devoted to the derivation of an approximate expression of the macroscopic yield function of the porous medium. To this end, the macroscopic plastic dissipation is first computed and then used with an appropriate minimization procedure to derive the yield function. In order to assess the accuracy of the established results, we compare, in Section 5, the yield surface corresponding to the closed-form expression of the new macroscopic criterion with that derived from the “exact” two-field criterion which is obtained numerically without any approximation. Moreover, the predictions of the new criterion are compared with existing results. Finally, the particular case of ductile media containing penny-shaped cracks is analyzed in order to emphasize some key features of the new criterion. Note that the very technical details of the study are reported in appropriate appendices.

## 2. Basic principles and methodology

### 2.1. The studied cell

Following (Gologanu et al., 1997), let us first consider an axisymmetric ellipsoidal cavity of semi-axes $a_1$ and $b_1$ embedded in a confocal ellipsoid of semi-axes $a_2$ and $b_2$. Fig. 1 depicts the cell relatively to the $(x_1, x_2, x_3)$ cartesian coordinates system of orthonormal basis ($e_1, e_2, e_3$), $x_3$ being aligned with the axis of symmetry of the void. The shape of the cavity is defined through the aspect ratio $a_1/b_1$, with $a_1 > b_1$ corresponding to a prolate cavity while $b_1 > a_1$ is associated to an oblate void. The focal distance is denoted $c$; $e_1$ and $e_2$ represent the void eccentricity and the exterior boundary eccentricity respectively:

$$ 
\begin{aligned}
c &= \sqrt{a_1^2 - b_1^2} = \sqrt{a_2^2 - b_2^2} \quad e_1 = \frac{c}{a_1} \quad e_2 = \frac{c}{a_2} \quad \text{(prolate)} \\
c &= \sqrt{b_1^2 - a_1^2} = \sqrt{b_2^2 - a_2^2} \quad e_1 = \frac{c}{b_1} \quad e_2 = \frac{c}{b_2} \quad \text{(oblate)}
\end{aligned} 
\quad (1) $$

Both cylindrical coordinates $\rho$, $\theta$, $z$ (with ($e_\rho, e_\theta, e_z$) the associated orthonormal basis) and the spheroidal coordinates $\lambda$, $\varphi$, $\theta$ (with associated orthogonal basis ($e_\lambda, e_\varphi, e_\theta$) shown in Fig. 2) are considered:

$$ 
\left\{ 
\begin{array}{l}
x_1 = b \sin(\varphi) \cos(\theta) \\
x_2 = b \sin(\varphi) \sin(\theta) \\
x_3 = a \cos(\varphi)
\end{array} 
\right. 
\quad 
\left\{ 
\begin{array}{l}
\rho = b \sin(\varphi) \\
\theta = \varphi \\
z = x_3 = a \cos(\varphi)
\end{array} 
\right. 
\quad (2) $$

and

$$ 
\left\{ 
\begin{array}{l}
\underline{e}_{\lambda} = \frac{1}{L_{\lambda}} \{a \sin(\varphi) \underline{e}_{\rho} + b \cos(\varphi) \underline{e}_{3}\} \\
\underline{e}_{\varphi} = \frac{1}{L_{\varphi}} \{b \cos(\varphi) \underline{e}_{\rho} - a \sin(\varphi) \underline{e}_{3}\} \\
\underline{e}_{\theta} = \underline{e}_{\theta}
\end{array} 
\right. 
\quad (3) $$

<footnote n="1">¹ Note that, Qiu and Weng (1993) derived similar criterion of porous media in the particular case of spheroidal aligned voids.</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** The considered cell: (left) prolate and (right) oblate spheroidal void embedded in a confocal spheroid relative to a cartesian coordinates system ($x_1, x_2, x_3$).</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Illustration of the spheroidal coordinates in the case of a prolate void.</figcaption></figure>

in which $L_\lambda = \sqrt{a^2 \sin^2(\varphi) + b^2 \cos^2(\varphi)}$, $\varphi \in [0, \pi]$ and $e_\rho = \cos(\theta)e_1 + \sin(\theta)e_2$, $\theta \in [0, 2\pi]$.

In the above equations: $a = c \cosh(\lambda)$ and $b = c \sinh(\lambda)$ for a prolate void, while for the case of an oblate void $a = c \sinh(\lambda)$ and $b = c \cosh(\lambda)$ with $\lambda \in [0, +\infty]$. The iso-$\lambda$ surfaces define confocal spheroids with foci $c = \sqrt{a^2 - b^2}$ and eccentricity $e = c/a$, for a prolate void, $e = c/b$ for an oblate one.

The porosity $f$ and the void shape parameter $S$ are given by:

$$f = \frac{a_1 b_2^2}{a_2 b_1^2}; \quad S = \ln\left(\frac{a_1}{b_1}\right) \qquad (4)$$

Due to the confocality of the exterior and interior spheroids, the porosity $f$ can be equally expressed in terms of the eccentricities $e_1$ and $e_2$ as:

$$f = \begin{cases} \frac{e_2^2}{e_1^2} \frac{1 - e_1^2}{1 - e_2^2} & (\text{prolate void}) \\ \frac{e_2^2}{e_1^2} \frac{\sqrt{1 - e_1^2}}{\sqrt{1 - e_2^2}} & (\text{oblate void}) \end{cases} \qquad (5)$$

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

It must then be noticed that the porosity is fixed when the two eccentricities $e_1$ and $e_2$ are given. Conversely, when the porosity $f$ and the void shape $S$ are fixed, the shape of the unit cell and that of the void are determined and cannot been chosen arbitrarily.

## 2.2. Principle of determination of the macroscopic yield surface

The matrix of the spheroidal unit cell is made up of a rigid-ideal plastic material obeying to the von Mises yield criterion with an associated flow rule as required in limit analysis theory. The local plastic dissipation, denoted $\pi(\mathbf{d})$, is defined, for every traceless strain rate field $\mathbf{d}$, as:

$$\pi(\mathbf{d}) = \begin{cases} \sigma_0 d_{eq} & (\text{in the matrix}) \\ 0 & (\text{in the void}) \end{cases} \qquad (6)$$

where $\sigma_0$ represents the yield stress in uniaxial tension and $d_{eq}$ the von Mises equivalent local strain rate, $d_{eq} = \sqrt{\frac{2}{3} \mathbf{d} : \mathbf{d}}$. The following inequality then holds for all macroscopic stresses $\Sigma$ and macroscopic strain rates $\mathbf{D}$ (see for instance Suquet, 1982):

$$\Sigma : \mathbf{D} \le \Pi(\mathbf{D}) = \inf_{\nu, K, A} \left[ \frac{1}{|\Omega|} \int_{\Omega-\omega} \sigma_0 d_{eq} dV \right] \qquad (7)$$

In the above expression, $\Pi(\mathbf{D})$ represents the macroscopic dissipation, $\Omega$ denotes the domain of the unit cell, $|\Omega| = 4\pi a_2 b_2^2/3$, whereas $\omega$ corresponds to the domain occupied by the void, $|\omega| = 4\pi a_1 b_1^2/3$. The infimum is taken over all kinematically admissible (K.A) velocity fields, $\nu$.

The volume integral over the matrix (volume $\Omega - \omega$) can be expressed in the spheroidal coordinates system as:

$$\frac{1}{|\Omega|} \int_{\Omega-\omega} \sigma_0 d_{eq} dV = \frac{3\sigma_0}{4\pi a_2 b_2^2} \int_{\lambda=\lambda_1}^{\lambda=\lambda_2} \int_{\varphi=0}^{\varphi=\pi} \int_{\theta=0}^{\theta=2\pi} d_{eq} b L_{\lambda}^2 \\ \times \sin \varphi d\lambda d\varphi d\theta \qquad (8)$$

The macroscopic yield locus is classically shown to be deduced from $\Pi(\mathbf{D})$ in the form:

$$\Sigma = \frac{\partial \Pi}{\partial \mathbf{D}} \qquad (9)$$
As an example, in the Gurson approach Gurson (1977), instead of seeking the infimum in Eq. (7), $\Pi(D)$, is estimated by choosing a particular microscopic velocity field $v$ complying with uniform strain rate boundary conditions $v = D \cdot \mathbf{x}$. In the solid, the latter is defined as the sum of a linear part, involving a homogeneous traceless second order tensor $A$, and $v^E$ which is the solution to the problem of a hollow sphere loaded hydrostatically, i.e.:

$$v = A \cdot \mathbf{x} + v^E \qquad (10)$$

hence the local strain rate $d$ is obtained as:

$$d = A + d^E \qquad (11)$$

$d^E$ being the strain rate tensor associated to $v^E$.

Extensions of the Gurson limit analysis approach to materials containing spheroidal cavities have been proposed by several authors (see for instance Gologanu et al., 1993, 1994, 1997; Gologanu, 1997; Gârâjeu, 1995).<fnref n="2" /> These studies consider the same type of decomposition of the velocity field $v$ still complying with uniform strain rate boundary conditions. The choice for $v^E$ is inspired from the work of Lee and Mear (1992). These authors proposed a family of axisymmetric velocity fields to investigate the growth of a spheroidal cavity in an infinite incompressible viscoplastic matrix and subjected to a remote axisymmetric stress loading. As mentioned before, it is desirable to generalize the proposed extensions by introducing non-axisymmetric fields still in the context of limit analysis theory. This is the purpose of the next section.

### 3. The trial velocity field

Let us first indicate that, as previous authors, the trial velocity field is taken in the form given by (10). Our main proposal consists in considering the solution to the Eshelby equivalent inclusion problem (Eshelby, 1959) (see also Mura, 1987) for the heterogeneous part $v^E$ in (10).

#### 3.1. Expression of the Eshelby-like velocity fields

To this end, consider an infinite, incompressible, linear viscous medium in which a spheroidal subdomain (defined by its semi-axes $a_1$ and $b_1$) is subjected to a uniform arbitrary eigenstrain rate $d'$. The velocity field $v^E$ for an interior point (within the spheroidal subdomain) and an exterior point (outside the spheroidal subdomain) reads in the form given by Mura (1987) and reported in Appendix A (see Eq. (A.2)). In the case of spheroidal inclusions, $v^E$ reduces to:

$$v^E = \begin{cases} [T(\lambda) : d^*] \cdot \mathbf{x} + \frac{a_1 b_1^2}{a^2 b^3} (a^2 - a_1^2) L_\lambda d_{\lambda\lambda}^* e_\lambda & \text{for } \mathbf{x} \in \Omega - \omega \\ [T(\lambda_1) : d^*] \cdot \mathbf{x} & \text{for } \mathbf{x} \in \omega \end{cases} \qquad (12)$$

The corresponding computation is detailed in Appendix A.

In (12), the quantity $d_{\lambda\lambda}^*$ is defined by $d_{\lambda\lambda}^* = d^* : (e_\lambda \otimes e_\lambda)$. The fourth order tensor $T(\lambda)$ has not the major and first

minor symmetry ($T_{ijkl} \neq T_{klij}$ and $T_{ijkl} \neq T_{jikl}$ but $T_{ijkl} = T_{ijkl}$). Its components, which are only functions of the coordinate $\lambda$, are given by:

$$T_{1111} = T_{2222} = 3T_{1122} = 3T_{2211} = \frac{3a_1 b_1^4}{8ab^4} (3 - 3\alpha - \beta)$$

$$T_{1212} = T_{2112} = T_{1122}$$

$$T_{1133} = T_{2233} = \frac{a_1^3 b_1^2}{2a^3 b^2} \beta$$

$$T_{3311} = T_{3322} = \frac{a_1 b_1^4}{2a^3 b^2} \beta$$

$$T_{3333} = \frac{a_1^3 b_1^2}{a^3 b^2} (1 - \beta)$$

$$T_{1313} = T_{2323} = \frac{a_1 b_1^4}{2a^3 b^2} \beta$$

$$T_{3113} = T_{3223} = \frac{a_1^3 b_1^2}{2a^3 b^2} \beta$$

(13)

where $\alpha$ and $\beta$, function of $\lambda$, can be also be expressed in term of eccentricity $e$:

$$\alpha(e) = \begin{cases} \frac{1-e^2}{e^2} \operatorname{arctanh} \{e\} - \frac{1-e^2}{e^2} & \text{(prolate void)} \\ -\frac{\sqrt{1-e^2}}{e^3} \operatorname{arctan} \left\{\frac{e}{\sqrt{1-e^2}}\right\} + \frac{1}{e^2} & \text{(oblate void)} \end{cases} \qquad (14)$$

and

$$\beta(e) = \begin{cases} (1 - 3\alpha(e)) \frac{1}{e^2} & \text{(prolate void)} \\ -(1 - 3\alpha(e)) \frac{1 - e^2}{e^2} & \text{(oblate void)} \end{cases} \qquad (15)$$

Introduction of $\alpha(e)$ and $\beta(e)$ presents the advantage to derive various results using unified expressions for both oblate and prolate cavities. The variation of $\alpha(e)$ and $\beta(e)$ is shown on Figs. 3 and 4, respectively. We also adopt the following notations: $\alpha_1 = \alpha(e_1)$ and $\beta_1 = \beta(e_1)$ and similarly $\alpha_2 = \alpha(e_2)$ and $\beta_2 = \beta(e_2)$. Values of these quantities in the particular cases of special interest are the following. For cylindrical voids, one has $e_1 \to 1$, $e_2 \to 1$ and then $\alpha_1 = \alpha_2 = 0$, $\beta_1 = \beta_2 = 1$. In the case of a spherical cavity, $e_1 \to 0$, $e_2 \to 0$ we have $\alpha_1 = \alpha_2 = 1/3$, $\beta_1 = \beta_2 = 2/5$. Finally, for penny shaped cracks $e_1 \to 0$ and $e_2$ arbitrary, $\alpha_1 = 1$, $\beta_1 = 0$.

Note that the second term in the expression (12) for $\mathbf{x} \in \Omega - \omega$, is zero at $\lambda = \lambda_1$ which ensures the continuity of $v^E$ through the boundary of the spheroidal subdomain.

We now propose to compare the “exterior point” velocity field in (12) with the one that has been used by Gologanu et al. (1993, 1994) in the context of ductile porous media. For convenience, we decompose the eigenstrain rate $d^*$ appearing in (12) as follows:

$$d^* = \sum_{n=1}^{n=6} d_n^* Q_n \qquad (16)$$

where the $Q_n$ for $n = 1, \dots, 6$ are given by:

<footnote n="2">² In Gologanu (1997) and Gologanu et al. (1997), supplementary fields have been introduced in order to take into account the finiteness of the studied unit cell and to obtain a more refined numerical estimation of the yield surface.</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 3. Representation $\alpha$ as function of $e$ for a prolate and an oblate void.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. Representation $\beta$ as function of $e$ for a prolate and an oblate void.</figcaption></figure>

$$ 
\begin{aligned}
\mathbf{Q}_1 &= \mathbf{1}; \quad \mathbf{Q}_2 = \mathbf{1} - 3\mathbf{e}_3 \otimes \mathbf{e}_3; \\
\mathbf{Q}_3 &= \mathbf{e}_2 \otimes \mathbf{e}_2 - \mathbf{e}_1 \otimes \mathbf{e}_1 \quad \mathbf{Q}_4 = \mathbf{e}_1 \otimes \mathbf{e}_2 + \mathbf{e}_2 \otimes \mathbf{e}_1; \\
\mathbf{Q}_5 &= \mathbf{e}_1 \otimes \mathbf{e}_3 + \mathbf{e}_3 \otimes \mathbf{e}_1 \quad \mathbf{Q}_6 = \mathbf{e}_2 \otimes \mathbf{e}_3 + \mathbf{e}_3 \otimes \mathbf{e}_2
\end{aligned} 
\quad (17) $$

where $d_1^* = d_m^*$ represents the mean part of $\mathbf{d}^*$.

The components of $\nu^E$ according to the above decomposition are given in Appendix A (see Eqs. (A.10)–(A.12)). Note that the fields associated to $d_1^*$ and $d_2^*$ are axisymmetric while the ones corresponding to $d_3^*, d_4^*, d_5^*$ and $d_6^*$ are non axisymmetric. Indeed, the axisymmetric velocity field in (12) can be expressed as a combination of the three fields corresponding to coefficients $B_{00}$, $B_{20}$ and $B_{22}$ (using the Lee and Mear, 1992 notations). The expansions fields considered by Gologanu et al. (1993, 1994) and Gărațăjeu (1995), for their two-field estimate of the macroscopic yield function in the case of spheroidal cavities, correspond to the one proportional to $B_{00}$ and $B_{22}$ in Lee and Mear (1992). It must be emphasized that those fields are also contained in the exterior point solution in (12) and are associated to $\mathbf{d}^* = d_m^*\mathbf{1}$.

In order to reproduce numerical results in the case of oblate cavities and more specially for penny-shaped cracks, Gologanu (1997) and Gologanu et al. (1997) proposed to include two additional fields corresponding to coefficients $B_{20}$ and $B_{21}$. Note that the axisymmetric field associated to $B_{21}$ is not contained in (12).

The strain rate tensor associated to $\nu^E$ reads:

$$ \mathbf{d}^E = \mathbb{D} : \mathbf{d}^* \quad \text{with} : \quad \mathbb{D} = \begin{cases} \mathbb{T}^s(\lambda) + \mathbb{H}(\lambda, \varphi, \theta) & \text{for } \mathbf{x} \in \Omega - \omega \\ \mathfrak{S}(e_1) & \text{for } \mathbf{x} \in \omega \end{cases} \quad (18) $$

in which $\mathbb{T}^s(\lambda)$ is given by:

$$ T_{ijkl}^{s}(\lambda) = \frac{1}{2}(T_{ijkl}(\lambda) + T_{jikl}(\lambda)) \quad (19) $$

In Eq. (18), $\mathfrak{S}(e)$ denotes the classical fourth order Eshelby tensor for a linear viscous incompressible medium which is a function of the eccentricity $e$. Its components, in the cartesian coordinates system, are given by:

$$ 
\begin{aligned}
S_{1111} &= S_{2222} = 3S_{1122} = 3S_{2211} = 3S_{1212} = \frac{3}{8}(3 - 3\alpha - \beta) \\
S_{1133} &= S_{2233} = \frac{1}{2}\beta \\
S_{3311} &= S_{3322} = \frac{1}{2}(\beta + 3\alpha - 1) \\
S_{3333} &= 1 - \beta \\
S_{1313} &= S_{2323} = \frac{1}{4}(2\beta + 3\alpha - 1)
\end{aligned} 
\quad (20) $$

In (18), the components of the fourth order tensor $\mathbb{H}$ are function of $\lambda$, $\varphi$ and $\theta$ and are given by:
$$
\begin{align}
H_{\lambda\lambda\lambda\lambda} &= -\frac{a_1 b_1^2}{ab^2} \left[ \frac{2(a^2 - a_1^2)}{L_\lambda^2} + \frac{a_1^2}{a^2} \right]; & H_{\lambda\lambda\phi\phi} &= \frac{a_1 b_1^2}{ab^2} \frac{1}{L_\lambda^2} \left[ a^2 - a_1^2 - L_\lambda^2 \right] \nonumber \\
H_{\phi\phi\lambda\lambda} &= \frac{a_1 b_1^2}{ab^2} \frac{a^2 - a_1^2}{L_\lambda^2}; & H_{\lambda\lambda 00} &= -\frac{a_1 b_1^4}{ab^4} \nonumber \\
H_{00\lambda\lambda} &= \frac{a_1 b_1^2}{ab^4} (a^2 - a_1^2); & H_{\lambda\phi\lambda\phi} &= \frac{a_1 b_1^2}{ab^2} \frac{1}{2L_\lambda^2} \left[ 2(a^2 - a_1^2) - L_\lambda^2 \right] \nonumber \\
H_{\lambda 0\lambda 0} &= \frac{a_1 b_1^2}{2ab^4} \left[ b^2 - 2b_1^2 \right] \nonumber \\
H_{\lambda\lambda\lambda\phi} &= H_{\lambda\phi\lambda\lambda} = \frac{a_1 b_1^2}{a^2 b^3} \frac{\sin(2\phi)}{2L_\lambda^2} (a^2 - a_1^2)(a^2 - b^2) \tag{21}
\end{align}
$$

Note that the components of $\mathbf{D}$ in the cartesian coordinates system can be found in Mura (1987). $\mathbf{D}$ has minor symmetry but not the major one; moreover, due to the matrix incompressibility, one has $\mathbf{1} : \mathbf{D} = 0$.

## 3.2. Conditions associated to the Eshelby-like velocity field

Let us come back now to the problem of the ductile porous material, for which the inclusions correspond to spherical cavities. Substituting the Eshelby-like velocity field in (10), it appears that the trial velocity field $\underline{v}$ depends on 11 undefined kinematical parameters, the components of the uniform strain rate $\mathbf{A}$ (5 independent kinematical parameters due to $\text{tr}(\mathbf{A}) = 0$ which corresponds to the matrix incompressibility) and the components of the eigenstrain rate tensor $\mathbf{d}^*$ (6 independent components). The verification of uniform strain rate boundary conditions leads to:

$$ \mathbf{D} \cdot \underline{\mathbf{x}} = \mathbf{A} \cdot \underline{\mathbf{x}} + [\mathbb{T}(\lambda_2) : \mathbf{d}^*] \underline{\mathbf{x}} + f \frac{a_2^2 - a_1^2}{a_2 b_2} L_\lambda d_{\lambda\lambda}^* e_\lambda \quad (22) $$

Important observations are the followings:

* the first two terms on the right hand of the equality (22) being proportional to $\underline{x}$, they are then compatible with uniform strain rate boundary conditions,

* it can be shown that the last term on the right-hand side of (22) is proportional to $\underline{x}$ if and only if $\mathbf{d}^* = \mathbf{d}_m^* \mathbf{1}$.

Consequently, a necessary condition for the considered velocity field to comply with uniform strain rate boundary conditions is $\mathbf{d}^* = 0$, which leads to:

$$
\begin{align}
D_{11} &= A_{11} + f \frac{3}{2} (1 - \alpha_2) d_m^* \nonumber \\
D_{22} &= A_{22} + f \frac{3}{2} (1 - \alpha_2) d_m^* \tag{23} \\
D_{33} &= A_{33} + 3f \alpha_2 d_m^* \nonumber \\
D_{12} &= A_{12}; \quad D_{13} = A_{13}; \quad D_{23} = A_{23} \nonumber
\end{align}
$$

Under these conditions, the velocity field (10) reduces to the axisymmetric one used by Gologanu et al. (1993, 1994).

It follows that in the general case ($\mathbf{d}^* \neq 0$) the complete Eshelby-based velocity field does not comply with uniform strain rate boundary conditions for the finite hollow spheroid. However, the use of the complete Eshelby-based

velocity field in the kinematical limit analysis approach is rigorously justified for uniform stress boundary conditions, $\boldsymbol{\sigma} \cdot \mathbf{n} = \boldsymbol{\Sigma} \cdot \mathbf{n}$, ($\mathbf{n}$ being the outward normal unit vector taken on the external spheroid $\lambda = \lambda_2$). It follows from the Hill-Mandel lemma that the macroscopic strain rate $\mathbf{D}$ is defined as the volume average of the local strain rate $\mathbf{d}$:

$$ \mathbf{D} = \langle \mathbf{d} \rangle_{\Omega} = \frac{1}{\Omega} \int_{\Omega - \omega} \mathbf{d} \mathbf{d} V \quad (24) $$

Note that for uniform stress boundary conditions, the kinematical approach, may still provide an upper bound for the macroscopic yield locus.

Substituting (11) (with the definition (18) for $\mathbf{d}^E$) in Eq. (24) leads to:

$$ \mathbf{D} = \mathbf{A} + f \mathfrak{S}(e_2) : \mathbf{d}^* \quad (25) $$

where $\mathfrak{S}(e_2)$ is the fourth order Eshelby tensor given by (20) for an incompressible medium with a spheroidal subdomain of eccentricity $e_2$. Further, the components of $\mathbf{D}$ in (25) reads:

$$
\begin{align}
D_{11} &= A_{11} + f \left[ \frac{1}{8} (3 - 3\alpha_2 - \beta_2) (3d_{11}^* + d_{22}^*) + \frac{\beta_2}{2} d_{33}^* \right] \nonumber \\
D_{22} &= A_{22} + f \left[ \frac{1}{8} (3 - 3\alpha_2 - \beta_2) (d_{11}^* + 3d_{22}^*) + \frac{\beta_2}{2} d_{33}^* \right] \nonumber \\
D_{33} &= A_{33} + f \left[ \frac{1}{2} (3\alpha_2 + \beta_2 - 1) (d_{11}^* + d_{22}^*) + (1 - \beta_2) d_{33}^* \right] \nonumber \\
D_{12} &= A_{12} + f \frac{1}{4} (3 - 3\alpha_2 - \beta_2) d_{12}^* \nonumber \\
D_{13} &= A_{13} + f \frac{1}{2} (3\alpha_2 + 2\beta_2 - 1) d_{13}^* \nonumber \\
D_{23} &= A_{23} + f \frac{1}{2} (3\alpha_2 + 2\beta_2 - 1) d_{23}^* \tag{26}
\end{align}
$$

Relation (25), or equivalently (26), provides a link between the macroscopic strain rate tensor $\mathbf{D}$ and all the unknown kinematical parameters (components of $\mathbf{A}$ and $\mathbf{d}^*$) which enter into the definition of the velocity field defined by (10) together with (12). The mean part of $\mathbf{D}$ computed from (25) reads:

$$ D_m = \frac{d_m^*}{f} \quad (27) $$

It follows that, $d_m^*$ is known and explicitly given in terms of the macroscopic mean strain rate. Moreover, it is interesting to note that, by considering the particular case of an isotropic eigenstrain ($\mathbf{d}^* = 0$), (26) reduces to (23), that is the mapping between the average rule and the uniform strain boundary conditions.

It is readily seen that the complete velocity field is defined by 11 parameters. Conditions (25), detailed in (26) provide six relations between these parameters. Therefore, there remain five unknown parameters which have to be determined. These are the components of $\mathbf{d}^*$.<fnref n="3" />

<footnote n="3">³ It must be emphasized that the components of $\mathbf{A}$ could be also chosen as unknowns but this does not change the final results.</footnote>
## 4. Determination of a closed-form expression of the macroscopic yield locus

For the derivation of a closed-form expression of the macroscopic dissipation, we follow here a methodology already used in Monchiet et al. (2006, 2007, 2008). The latter is inspired of the one introduced by Gologanu et al. (1993, 1994), but it differs in several points due to the remaining unknown $\bar{\mathbf{d}}^*$ which characterize the non axisymmetry of the considered velocity fields. The macroscopic dissipation reads then:

$$ \Pi(\mathbf{D}) = \inf_{\bar{\mathbf{d}}^*} \Gamma(\mathbf{D}, \bar{\mathbf{d}}^*) \quad (28) $$

where $\Gamma(\mathbf{D}, \bar{\mathbf{d}}^*)$ is defined by:

$$ \begin{aligned} \Gamma(\mathbf{D}, \bar{\mathbf{d}}^*) &= \frac{1}{|\Omega|} \int_{\Omega \to \infty} \sigma_0 d_{eq} dV \\ &= \frac{3\sigma_0}{4\pi a_2 b_2^2} \int_{\lambda_1 \to \lambda_2} \int_{\varphi=0}^{\varphi=\pi} \int_{\theta=0}^{\theta=2\pi} d_{eq} b L_{\lambda}^2 \\ & \quad \times \sin \varphi d\lambda d\varphi d\theta \end{aligned} \quad (29) $$

in which $d_{eq}$ takes the form:

$$ d_{eq}^2 = A_{eq}^2 + \frac{4}{3} \mathbf{A} : \mathbf{D} : \mathbf{d}^* + \frac{2}{3} \mathbf{d}^* : \mathbf{D}^T : \mathbf{D} : \mathbf{d}^* \quad (30) $$

with $A_{eq}^2 = \frac{2}{3} \mathbf{A} : \mathbf{A}$. The fourth order tensor $\mathbf{D}$ has been introduced in the previous section (see Eqs. (18), (19) and (21)).

Finally, the macroscopic yield locus is obtained from (9) with the definition (28).

The approximate expression of $\Gamma(\mathbf{D}, \bar{\mathbf{d}}^*)$ is established in Appendix B, while the determination of the macroscopic yield function, after a minimization procedure, is presented in Section 4.1. One has (see Appendix B):

$$ \begin{aligned} \Gamma(\mathbf{D}, \bar{\mathbf{d}}^*) &= -\sigma_0 f \int_{u_1}^{u_2} \left\{ Y^2 + X^2 u^2 \right\}^{1/2} \frac{du}{u^2} \\ &= \sigma_0 f \left[ X \operatorname{arcsinh} \left\{ \frac{uX}{Y} \right\} - \frac{\sqrt{Y^2 + u^2 X^2}}{u} \right]_{u_1}^{u_2} \end{aligned} \quad (31) $$

with

$$ X^2 = \frac{2}{3} [\mathbf{d}^* + \mathbf{A} : \mathbf{Q} : \mathbf{P}^{-1}] : \mathbf{P} : [\mathbf{d}^* + \mathbf{A} : \mathbf{Q} : \mathbf{P}^{-1}] \quad (32) $$

$$ Y^2 = A_{eq}^2 - \frac{2f^2}{3(1+g)(g+f)} \mathbf{A} : \mathbf{Q} : \mathbf{P}^{-1} : \mathbf{Q}^T : \mathbf{A} $$

Note that the derivation of $\Gamma(\mathbf{D}, \bar{\mathbf{d}}^*)$ has required three approximations among which $\mathcal{A}2$ and $\mathcal{A}3$ that do not preserve the upper bound character of the approach. This will be assessed by means of numerical solutions.

### 4.1. The minimization procedure and derivation of the macroscopic yield function

We now derive a closed-form expression of the macroscopic criterion. To this aim, one has to compute the dissipation (see Eq. (7)) for which the minimum of $\Gamma(\mathbf{D}, \bar{\mathbf{d}}^*)$ with respect to $\bar{\mathbf{d}}^*$ is taken. Once this minimization is performed, the macroscopic criterion can be obtained from (9). These two steps are performed simultaneously as:

$$ \Sigma = \frac{\partial \Gamma(\mathbf{D}, \bar{\mathbf{d}}^*)}{\partial \mathbf{D}} \quad \text{with:} \quad \frac{\partial \Gamma(\mathbf{D}, \bar{\mathbf{d}}^*)}{\partial \bar{\mathbf{d}}^*} = 0 \quad (33) $$

where the approximate expression of $\Gamma(\mathbf{D}, \bar{\mathbf{d}}^*)$, given by (31), is used. To solve (33), it is convenient to introduce the following change of variable: $\Gamma(\mathbf{D}, \bar{\mathbf{d}}^*) \equiv \Gamma(\mathbf{X}, \mathbf{Y})$. Consequently, the macroscopic stress tensor, $\Sigma$, reads:

$$ \Sigma = f \frac{\partial X}{\partial \mathbf{D}} \Sigma_X + \frac{\partial Y}{\partial \mathbf{D}} \Sigma_Y \quad (34) $$

with

$$ \Sigma_X = \frac{1}{f} \frac{\partial \Gamma}{\partial \mathbf{X}}; \quad \Sigma_Y = \frac{\partial \Gamma}{\partial \mathbf{Y}} \quad (35) $$

On the other hand, the minimum of $\Gamma(\mathbf{D}, \bar{\mathbf{d}}^*)$ with respect to $\bar{\mathbf{d}}^*$ is given by:

$$ f \frac{\partial X}{\partial \bar{\mathbf{d}}^*} \Sigma_X + \frac{\partial Y}{\partial \bar{\mathbf{d}}^*} \Sigma_Y = 0 \quad (36) $$

By using now the approximate expression of $\Gamma(\mathbf{D}, \bar{\mathbf{d}}^*)$, given by (31), in relations (35), one gets the following expressions for $\Sigma_X$ and $\Sigma_Y$:

$$ \begin{aligned} \Sigma_X &= \sigma_0 [\operatorname{arcsinh}(u_2 \xi) - \operatorname{arcsinh}(u_1 \xi)] \\ \Sigma_Y &= -\sigma_0 f \left[ \frac{\sqrt{1 + u_2^2 \xi^2}}{u_2} - \frac{\sqrt{1 + u_1^2 \xi^2}}{u_1} \right] \end{aligned} \quad (37) $$

with $\xi = X/Y$. By eliminating $\xi$ between the last two relations, one obtains:

$$ \begin{aligned} \Phi(\Sigma) &= \frac{\Sigma_Y^2}{\sigma_0^2} + 2(1+g)(f+g) \cosh \left( \frac{\Sigma_X}{\sigma_0} \right) - (1+g)^2 \\ & \quad - (f+g)^2 = 0 \end{aligned} \quad (38) $$

In the expression (38), expressions of $\Sigma_X$ and $\Sigma_Y$ as function of $\Sigma$ are still unknown. Their expressions are determined from (34) together with (36), following a procedure which is detailed in Appendix D. This leads to:

$$ \begin{aligned} \Sigma_X^2 &= \frac{3}{2} \Sigma : \mathcal{S}(e_2) : \mathcal{P}^{-1} : \mathcal{S}^T(e_2) : \Sigma \\ \Sigma_Y^2 &= \Sigma_{eq}^2 - (1+g)(f+g)\Sigma_X^2 \end{aligned} \quad (39) $$

with:

$$ \tilde{\Sigma}_{eq}^2 = \Sigma_{eq}^2 + \frac{3f}{2} \Sigma : [\mathcal{S}(e_1) - f\mathcal{S}(e_2)] : [\mathcal{L}(e_1) - f\mathcal{L}(e_2)]^{-1} : \Sigma \quad (40) $$

Eq. (38), together with (39) constitute the most important results of the present paper. Tensors $\mathcal{S}$ which appears in (39) and (40) is given by (20), while $\mathcal{P}$ and $\mathcal{L}$ are provided in (C.5) and (C.6) depicted in Appendix C. Let us recall that the parameter $g$ has been taken as zero in the case of a prolate void and given by Eq. (B.10) for an oblate void. The coefficient $\chi$ which enters in the definition of (B.10) is chosen in order to obtain a better agreement with “exact” numerical solutions of the two-field criteria that are presented in Section 5. The expression $\chi = (3\alpha_1 - 1)/4$ will be used.

Let us introduce the following macroscopic stress components related to the transverse isotropy induced by the geometry of the spheroidal cavity:
$$
\begin{align}
\Sigma_m &= \frac{1}{3}(\Sigma_{11} + \Sigma_{22} + \Sigma_{33}); \quad \Sigma_q = \frac{1}{2}(\Sigma_{11} + \Sigma_{22}) - \Sigma_{33} \nonumber \\
\Sigma_s &= \sqrt{\frac{1}{4}(\Sigma_{22} - \Sigma_{11})^2 + \Sigma_{12}^2}; \quad \Sigma_t = \sqrt{\Sigma_{13}^2 + \Sigma_{23}^2} \tag{41}
\end{align}
$$

With these notations, $\tilde{\Sigma}_{eq}$ can be expressed as follows:

$$
\tilde{\Sigma}_{eq}^2 = (1 + \eta_1)\Sigma_q^2 + \eta_2\Sigma_m^2 + \eta_3\Sigma_m\Sigma_q + 3(1 + \eta_4)\Sigma_s^2 \\
\qquad + 3(1 + \eta_5)\Sigma_t^2 \tag{42}
$$

where coefficients $\eta_i$, are given by:

$$
\begin{align}
\eta_1 &= \frac{4}{3\Upsilon} - 1; \quad \eta_2 = \frac{3}{2\Upsilon}(3\tilde{\alpha} + 3\tilde{\beta} - 1); \quad \eta_3 = \frac{2}{\Upsilon}(1 - 3\tilde{\alpha}) \nonumber \\
\eta_4 &= \frac{3 - 3\tilde{\alpha} - \tilde{\beta}}{1 + 3\tilde{\alpha} + \tilde{\beta}}; \quad \eta_5 = \frac{3\tilde{\alpha} + 2\tilde{\beta} - 1}{3 - 3\tilde{\alpha} - 2\tilde{\beta}} \tag{43}
\end{align}
$$

while $\Upsilon$ is defined by $\Upsilon = 2\tilde{\beta} + 4\tilde{\alpha} - 3\tilde{\alpha}^2 - 1$. $\tilde{\alpha}$ and $\tilde{\beta}$ are:

$$
\tilde{\alpha} = \frac{\alpha_1 - f\alpha_2}{1 - f}; \quad \tilde{\beta} = \frac{\beta_1 - f\beta_2}{1 - f} \tag{44}
$$

The decomposition of $\Sigma_X$ leads to:

$$
\Sigma_X^2 = \kappa_1 \Sigma_p^2 + \kappa_2 \Sigma_q^2 + \kappa_3 \Sigma_p \Sigma_q + \kappa_4 \Sigma_s^2 + \kappa_5 \Sigma_t^2 \quad (45)
$$

with:

$$
\Sigma_p = \Sigma_m + \frac{1}{3}(1 - 3\alpha_2)\Sigma_q \quad (46)
$$

and:

$$
\begin{align}
\kappa_1 &= \frac{9p_{22}}{p_{11}p_{22} - p_{12}^2}; \quad \kappa_2 = \frac{9p_{11}(1 - \alpha_2 - \beta_2)^2}{p_{11}p_{22} - p_{12}^2} \nonumber \\
\kappa_3 &= \frac{18p_{12}(1 - \alpha_2 - \beta_2)}{p_{11}p_{22} - p_{12}^2} \tag{47} \\
\kappa_4 &= \frac{(3 - 3\alpha_2 - \beta_2)^2}{4p_{33}}; \quad \kappa_5 = \frac{(2\beta_2 + 3\alpha_2 - 1)^2}{p_{55}} \nonumber
\end{align}
$$

where the quantities $p_{rs}$ are given by (C.4), in which the last terms proportional to $(1 - a_1^2/a_2^2)(1 - b_1^2/b_2^2)$ are neglected.

**Remark.** Due to the presence of non axisymmetric veloc-
ity field considered in the present study, the macroscopic
criterion shows, through (42) and (45), some new cou-
plings between pure shear components (Σs, Σt) and the
porosity. This point will be particularly emphasized and
discussed later in the case of penny-shaped cracks.

**5. Validation and comparisons with existing criteria**

First, it is convenient to notice that the results described
in the above section extend to spheroidal voids the one
established and illustrated in Monchiet et al. (2011) for
spherical voids. The later can be obtained here as a partic-
ular case.

5.1. The case of a cylindrical cavity

The case of a hollow cylinder corresponds to the limits
$e_1 \to 1$ and $e_2 \to 1$, from which one obtains for (43), (44)
and (47):

$$
\frac{\Sigma_{eq}^2}{\sigma_0^2} + 2f \cosh \left\{ \frac{1}{\sigma_0} \sqrt{\frac{3}{2} \left( \Sigma_{11}^2 + \Sigma_{22}^2 \right) + 3\Sigma_{12}^2 + 3\Sigma_{13}^2 + 3\Sigma_{23}^2} \right\} - 1 - f^2 \le 0 \quad (48)
$$

For comparison purpose, let us recall the expression of the
macroscopic yield function obtained by Gurson (1977) in
the case of a cylindrical cavity:

$$
\frac{\Sigma_{eq}^2}{\sigma_0^2} + 2f \cosh \left\{ \frac{\sqrt{3}}{2\sigma_0} (\Sigma_{11} + \Sigma_{22}) \right\} - 1 - f^2 \le 0 \quad (49)
$$

Note that (48) coincides with (49) for macroscopic axisym-
metric loading (Σ₁₁ = Σ₂₂ and Σ₁₂ = Σ₁₃ = Σ₂₃ = 0) for
which the results are exact. The reason is that the "axisym-
metric" part of the trial velocity field used in the present
study coincides with the one considered by Gurson and
correspond to parameters d₁ and d₂ in (12) (see Eq. (16)
for the definition of the dₙ). Moreover, note that the
approximation A1 is also used by Gurson who replaced
the equivalent strain rate, dₐₗ, by the square root of the
mean value of d²ₐₗ over the unit circle. Approximations
A2 and A3 are exact when the hollow cylinder is subjected
to an axisymmetric loading. Nevertheless, for general load-
ings, comparison between Gurson (49) and the criterion
(48) shows that the latter contains additional terms which
appear in the hyperbolic cosine. Those terms introduce
couplings between the porosity and the macroscopic trans-
verse and longitudinal shear stresses that does not exist in
(49). It is worth noticing that the presence of those compo-
nents of the macroscopic stress within the hyperbolic
cosine in (48) are related to the consideration of the
non-axisymmetric components of the trial velocity fields
associated to d₃, d₄, d₅, d₆ in (12).

5.2. Prolate cavities

Let us first recall that the determination of a closed-
form expression of the macroscopic criterion for spheroidal
cavities has required a number of “uncontrolled” approxi-
mations. In order to check the validity of those approxima-
tions, we shall first evaluate the accuracy of (38) by
comparing with numerical exact two-field solution. Figs. 5
and 6 show various yield loci for prolate cavities, for two
aspect ratios, a₁/b₁ = 2 and a₁/b₁ = 5 respectively. An axi-
symmetric macroscopic loading is considered,
Σ₁₁ = Σ₂₂ ≠ 0, Σ₃₃ ≠ 0 and all other components of Σ are
taken as zero. The criterion is then plotted in the plane
Σₐₗ = Σ₃₃ - Σ₁₁ versus Σₘ = (2Σ₁₁ + Σ₂₂)/3. On each figure
three values of the porosity are considered, namely
f = 0.01, f = 0.05 and f = 0.2. The full line corresponds
to the approximate criterion (38) whereas the discrete
points (the circles) refer to the exact two-field numerical
solution. It is observed that, irrespective of the value of
the aspect ratio and the porosity, there is a good agreement
between the “closed form” and “exact” solutions; this
$$ (\Sigma_{33} - \Sigma_{11})/\sigma_0 $$

<figure><img src="image_7.png" /><figcaption>Fig. 5. Yield loci for a prolate cavity having an aspect ratio $a_1/b_1 = 2$. Comparison between the approximate two-field criterion (38) and the numerical two-field criterion.</figcaption></figure>

$$ (\Sigma_{33} - \Sigma_{11})/\sigma_0 $$

<figure><img src="image_8.png" /><figcaption>Fig. 6. Yield loci for a prolate cavity having an aspect ratio $a_1/b_1 = 5$. Comparison between the approximate two-field criterion (38) and the numerical two-field criterion.</figcaption></figure>

validates for the prolate case all the approximations used in Section 4.

We now present comparisons of the new criterion (38) with existing criteria. As already mentioned, extensions of the Gurson model to the case of spheroidal cavities have

$$ \Sigma_s/\sigma_0 $$

<figure><img src="image_9.png" /><figcaption>Fig. 7. Yield stress for $\Sigma_s$ as function of the porosity $f$ for a prolate cavity having an aspect ratio $a_1/b_1 = 5$. Comparison between the approximate two-field criterion (38) and the exact (numerical) two-field criterion and the GLPD model.</figcaption></figure>

been proposed by Gologanu et al. (1993, 1994) (GLD model), Gologanu et al. (1997) (GLPD model) and Gârâjeu (1995) and Garajeu et al. (2000).<fnref n="4" /> Those models were also obtained in the context of Limit Analysis of a hollow conical spheroidal unit cell subjected to an axisymmetric loading.<fnref n="5" /> The GLPD model is expressed as:

$$ \frac{C}{\sigma_0^2}(\Sigma_{33} - \Sigma_{11} + \eta\Sigma_h)^2 + 2(1+g)(f+g)\cosh\left(\frac{\kappa\Sigma_h}{\sigma_0}\right) - (1+g)^2 - (f+g)^2 = 0 \quad (50) $$

In this criterion, the definition of the parameter $g$ differs from our definition (B.10) which has a multiplicative factor $\chi$. The coefficients, $C$, $\eta$ and $\kappa$ are function of the two eccentricities $e_1$, $e_2$, the parameter $g$ and of the porosity $f$ (see Gologanu, 1997 for the corresponding expressions). Note that the term $\Sigma_h$ linearly depends of the macroscopic stress components $\Sigma_{11}$ and $\Sigma_{33}$; this constitutes a difference with (38) together with (39) and (40).

In order to highlight the improvement obtained with the consideration of the Eshelby-like velocity fields, we now provide results for non axisymmetric loadings. In Fig. 7 we represent the yield stress for the component $\Sigma_s$, defined in Eq. (41), as function of the porosity for an aspect ratio $a_1/b_1 = 5$. The predictions of the new criterion (38) are compared to the exact two-field criterion (obtained numerically) and the prediction of the GLPD model. Note that the later coincides with their earlier criterion

<footnote n="4">⁴ In Gârâjeu (1995) and Garajeu et al. (2000) only the case of a prolate cavity has been studied.</footnote>

<footnote n="5">⁵ Note also that in Gologanu (1997) an heuristical extension has been proposed for arbitrary loading.</footnote>
(Gologanu et al., 1993) for non axisymmetric loadings and reduce to $\Sigma_s = (1-f)\sigma_0$. Fig. 8 provide similar results for the yield stress $\Sigma_t$. It can be observed that the new criterion improves the prediction of the GLPD model for non axisymmetric loading case. Again these improvements come from the consideration of the non-axisymmetric velocity fields which are contained in the Eshelby-like trial velocity fields.

### 5.3. The case of oblate cavities

We now consider the case of an oblate spheroidal cavity and represent the new yield criterion for axisymmetric loading. Fig. 9 shows the macroscopic yield locus for an aspect ratio $b_1/a_1 = 2$ whereas in Fig. 10, the aspect ratio is $b_1/a_1 = 5$. For low values of the porosity ($f = 0.01$ and $f = 0.05$), a good agreement of the new criterion with the two-field numerical data is observed. However for larger values of the porosity $f$, it can be observed that the exact two-field solution gives an estimate of the macroscopic yield locus which is more interior than the one given by the approximate criterion. These differences result, obviously, from the approximations used in Section 4.

We now provide comparisons of the new criterion with numerical upper and lower bounds of the macroscopic criterion, recently obtained by Pastor et al. (2011) by means of finite element-based limit analysis combined with convex optimization. In this numerical limit analysis, the macroscopic criterion has been determined by using a kinematical and a static approach to the problem of a confocal spheroidal cell subjected to an axisymmetric loading. This approach leads then to an upper bound and a lower bound of the macroscopic yield surface. On Figs. 11 and 12, we

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Yield stress for $\Sigma_t$ as function of the porosity $f$ for a prolate cavity having an aspect ratio $a_1/b_1 = 5$. Comparison between the approximate two-field criterion (38) and the exact (numerical) two-field criterion and the GLPD model.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Yield loci for an oblate cavity having an aspect ratio $b_1/a_1 = 2$. Comparison between the approximate two-field criterion (38) and the numerical two-field criterion.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Yield loci for an oblate cavity having an aspect ratio $b_1/a_1 = 5$. Comparison between the approximate two-field criterion (38) and the numerical two-field criterion.</figcaption></figure>

represent the yield locus for an aspect ratio $a_1/b_1 = 1/2$ and for the porosities $f = 0.01$ and $f = 0.1$ respectively. On these figures, we compare the numerical lower bound (static) and upper bound (kinematic) of Pastor et al. (2011) with the new criterion (38). For completeness, we
$$ (\Sigma_{33} - \Sigma_{11})/\sigma_0 $$

<figure><img src="image_13.png" /><figcaption>Fig. 11. Yield loci for oblate cavities with an aspect ratio $a_1/b_1 = 1/2$, for a porosity $f = 0.01$. Comparison of the new criterion (38), the GLPD model (Gologanu et al., 1997) and the numerical bounds of Pastor et al. (2011).</figcaption></figure>

$$ (\Sigma_{33} - \Sigma_{11})/\sigma_0 $$

<figure><img src="image_14.png" /><figcaption>Fig. 13. Yield loci for oblate cavities with an aspect ratio $a_1/b_1 = 1/5$ and for the porosity $f = 0.01$. Comparison of the new criterion (38), the GLPD model (Gologanu et al., 1997) and the numerical bounds of Pastor et al. (2011).</figcaption></figure>

$$ (\Sigma_{33} - \Sigma_{11})/\sigma_0 $$

<figure><img src="image_15.png" /><figcaption>Fig. 12. Yield loci for oblate cavities with an aspect ratio $a_1/b_1 = 1/2$, for a porosity $f = 0.1$. Comparison of the new criterion (38), the GLPD model (Gologanu et al., 1997) and the numerical bounds of Pastor et al. (2011).</figcaption></figure>

$$ (\Sigma_{33} - \Sigma_{11})/\sigma_0 $$

<figure><img src="image_16.png" /><figcaption>Fig. 14. Yield loci for oblate cavities with an aspect ratio $a_1/b_1 = 1/5$ and for the porosity $f = 0.1$. Comparison of the approximate criterion (38), the GLPD model (Gologanu et al., 1997) and the numerical bounds of Pastor et al. (2011).</figcaption></figure>

still report on these figures the prediction of Gologanu et al. (1997). It can be noted that the numerical bounds are very closed to each other and thus gives an accurate estimate of the macroscopic yield surface. Figs. 13 and 14 display

similar results for the spheroidal cavity with an aspect ratio $a_1/b_1 = 1/5$. On these figures, the numerical bounds are less accurate for high values of the stress triaxiality. It is generally observed that our results as well as that of
the GLPD model are in a good agreement with the numerical solutions. Note however that the results provided on Figs. 12 and 14, show that the GLPD model slightly overestimates the yield stress for low stress triaxialities.

Consider now the case of non axisymmetric loadings. In Fig. 15, we represent the yield stress for $\Sigma_s$ as function of the porosity while, in Fig. 16, similar results are shown for the yield stress $\Sigma_t$. The aspect ratio of the oblate cavity is $a_1/b_1 = 1/5$. Again, the predictions of the GLPD model coincide with the criterion of Gologanu et al. (1994) for non axisymmetric loadings, which reduce to $\Sigma_s = (1-f)\sigma_0$ and $\Sigma_t = (1-f)\sigma_0$. Our results clearly show the modifications obtained for $\Sigma_t$. It is observed that, for the component $\Sigma_s$, the prediction are very closed to the ones obtained with the GLPD model.

### 5.4. Penny-shaped cracks

This subsection is devoted to the important case of a circular crack, which is obtained by considering an oblate spheroidal cavity with an aspect ratio $a_1/b_1 \to 0$ or equivalently with an eccentricity $e_1 \to 1$. The porosity $f \to 0$ and $g \to d/2$ where $d$ is the crack density parameter given by $d = b_1^3/(a_2b_2^2)$. This latter parameter was introduced first by Bristow (1960) and later considered by Budiansky and O'Connell (1976). The new macroscopic criterion (38) reduces then to:

$$ \frac{\Sigma_Y^2}{\sigma_0^2} + d \left(1 + \frac{d}{2}\right) \cosh \left\{ \frac{\Sigma_X}{\sigma_0} \right\} - \left(1 + \frac{d}{2}\right)^2 - \frac{d^2}{4} = 0 \quad (51) $$

where $\Sigma_Y$, defined by (39), reads for the case of the penny-shaped cracks:

$$ \Sigma_s / \sigma_0 $$

<figure><img src="image_17.png" /><figcaption>Fig. 15. Yield stress for $\Sigma_s$ as function of the porosity $f$ for an oblate cavity having an aspect ratio $a_1/b_1 = 1/5$. Comparison between the approximate two-field criterion (38) and the exact (numerical) two-field criterion and the GLPD model.</figcaption></figure>

$$ \Sigma_t / \sigma_0 $$

<figure><img src="image_18.png" /><figcaption>Fig. 16. Yield stress for $\Sigma_t$ as function of the porosity $f$ for an oblate cavity having an aspect ratio $a_1/b_1 = 1/5$. Comparison between the approximate two-field criterion (38) and the exact (numerical) two-field criterion and the GLPD model.</figcaption></figure>

$$ \Sigma_Y^2 = \bar{\Sigma}_{eq}^2 - \frac{1}{4}(2+d)d\Sigma_X^2 \quad (52) $$

For the definition of the quantity $\bar{\Sigma}_{eq}$ we refer to (40) together with (43). When the limit $e_1 \to 1$ is taken, $\bar{\Sigma}_{eq}$ reduces to:

$$ \bar{\Sigma}_{eq}^2 = \Sigma_{eq}^2 + \frac{3d\Sigma_{33}^2}{\pi + 2(\alpha_2 - \beta_2 - 1)d} + \frac{12d(\Sigma_{13}^2 + \Sigma_{23}^2)}{3\pi + 2(3\alpha_2 + 2\beta_2 - 3)d} \quad (53) $$

For the computation of the limit $e_1 \to 1$ of the coefficients $\eta_i$ defined in (43), it is useful to put the porosity in the form $f = dz_1$ with $z_1 = \sqrt{1 - e_1^2}$ (in the case of a penny-shaped cracks, $z_1 \to 0$) and to take into account the following series:

$$ \alpha_1 = 1 - \frac{\pi}{2}z_1 + o(z_1^2) \quad \beta_1 = o(z_1^2) \quad (54) $$

Quantity $\Sigma_X$ is given by (45) together with (47) and (C.4). When the limit $e_1 \to 1$, coefficients $p_{rs}$ which appear in the definition of $\kappa_i$ reads:

$$ \begin{aligned} p_{11} &= \frac{3(2+d)}{4} [2\pi + d(1+3\alpha_2)(1-\alpha_2)] \\ p_{22} &= \frac{3(2+d)}{4} [\pi - d(3\alpha_2 + 3\beta_2 - 1)(1-\alpha_2 - \beta_2)] \\ p_{12} &= -\frac{3(2+d)}{4} [\pi + d(1-3\alpha_2)(1-\alpha_2 - \beta_2)] \\ p_{33} &= p_{44} = \frac{2+d}{48} [6\pi - d(1+3\alpha_2 + \beta_2)(3-3\alpha_2 - \beta_2)] \\ p_{55} &= p_{66} = \frac{2+d}{12} [3\pi - d(1-3\alpha_2 - 2\beta_2)(3\alpha_2 + 2\beta_2 - 3)] \end{aligned} \quad (55) $$
As an illustration, on Fig. 17 is plotted the macroscopic criterion for the following values of the crack density parameter: $d = 0.05$, $d = 0.5$ and $d = 0.5$. Again, a good agreement between the exact two-field criterion (circles) and the approximate one (full line) is observed. On Figs. 18 and 19, we represent the variations of the yield stress for $\Sigma_s$ and $\Sigma_t$ as function of the crack density parameter, $d$. On these figures, we compare the predictions of the approximate criterion (51) with the exact two-field

<figure><img src="image_19.png" /><figcaption>Fig. 17. Yield loci for a penny-shaped crack. Comparison between the approximate two-field criterion (51) and the numerical two-field criterion.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 18. Variations of the yield stress $\Sigma_s$ as function of the crack density parameter for the penny-shaped crack.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 19. Variations of the yield stress $\Sigma_t$ as function of the crack density parameter for the penny-shaped crack.</figcaption></figure>

criterion and the GLPD model. For the latter, the damage parameter $d$ does not affect the yield stress while the new criterion shows an important effect.

## 6. Conclusion

A new analytic yield criterion for a rigid perfectly plastic porous medium containing spheroidal cavities has been derived using a limit analysis based micro-macro approach. In the analysis, trial velocity fields inspired from the exterior point solution to the Eshelby inclusion problem were considered. An important feature of these trial velocity fields is that they contain non-axisymmetric components. The methodology followed in the study has required a rigorous computation of the macroscopic dissipation combined with an appropriate minimization procedure from which we derived a closed-form expression of the yield function at macroscale. Although the criterion has a form similar to the one established by Gologanu et al. (1993, 1994) (see also Gologanu et al., 1997, GLPD model) or by Garajeu et al. (2000) (or Gārājeu, 1995), it differs significantly. In particular, owing to the presence of shear stress components in the argument of the hyperbolic cosine, the new criterion captures the yielding response for purely deviatoric loadings. Indeed, in the case of a cylindrical cavity and in the more general case of spheroidal voids, the new criterion exhibits couplings between the transverse and longitudinal shear stresses which were not properly accounted in previous limit analysis based models. The effect of these couplings on the yield surface has been clearly shown in the present study. Note also that in the particular case of a penny-shaped crack and purely shear loadings, the predictions of the GLPD model reduce to the von Mises criterion while the established criterion shows an influence of the shear stress
components. Assessment of the accuracy of the new criterion has been carefully done by comparing its predictions with the “exact” two-field numerical solutions and the numerical lower and upper bounds of Pastor et al. (2011). An overall good agreement has been observed.

The new criterion being analytical, it can be easily implemented in FE codes and can be readily used for efficient and fast computation of damage in a variety of engineering structural applications.

Finally, it is interesting to note that, based on Monchiet (2006) and Monchiet et al. (2007), Pastor and Kondo (2013) recently provided numerical results for the limit analysis problem of the voided spheroidal cell; this brings complementary validation of our theoretical results. The ultimate goal is to formulate, implement and validate the constitutive law resulting from the present study. Such validation will be carried out by means of comparisons to experimental data concerning the response of mechanical structures. A first step in this direction has been already provided in Lin et al. (2010) on the isotropic version of the model (case of spherical voids).

## Appendix A. Detailed expression of the Eshelby-based components of the velocity field

An ellipsoidal inclusion, of an infinite viscous elastic matrix of rigidity $C^0$ is submitted to the constant eigenstrain rate $d^t$. We introduce by $(x_1, x_2, x_3)$ the cartesian coordinates and by $c_1, c_2, c_3$ the radii of the ellipsoid along the three axis of cartesian frame. The inclusion is located at the origin and its volume is defined by:

$$ \frac{x_1^2}{c_1^2} + \frac{x_2^2}{c_2^2} + \frac{x_3^2}{c_3^2} - 1 \le 0 \qquad (A.1) $$

The velocity has the form (see Kachanov et al., 2003, p. 234, and Mura, 1987 chapter 2):

$$ \begin{aligned} \forall \underline{x} \in \omega : 8\pi(1 - v^0) \underline{v}_I^E(\underline{x}) &= 3d_{kk}^* x_I [c_K^2 I_{KL}(0) - (1 - 2v^0)I_L(0)] \\ &\quad + d_{ij}^* [(c_I^2 + c_J^2)I_{Ly}(0) + (1 - 2v^0)(I_L(0) \\ &\quad + I_J(0)) + \frac{1}{2}(I_J(0) - I_L(0))] x_J \\ \forall \underline{x} \in \Omega - \omega : 8\pi(1 - v^0) \underline{v}_I^F(\underline{x}) &= 3d_{kk}^* x_I [c_K^2 I_{KL}(\underline{\xi}) - (1 - 2v^0)I_L(\underline{\xi})] \\ &\quad + d_{ij}^* [(c_I^2 + c_J^2)I_{Ly}(\underline{\xi}) + (1 - 2v^0)(I_L(\underline{\xi}) \\ &\quad + I_J(\underline{\xi})) + \frac{1}{2}(I_J(\underline{\xi}) - I_L(\underline{\xi}))] x_J \\ &\quad + (c_K^2 I_{KL}(\underline{\xi}) - I_L(\underline{\xi}))_i d_{kl}^* x_k x_l \end{aligned} \qquad (A.2) $$

The following summation convention has been used: repeated lower case indices are summed from 1 to 3; upper case indices take on the same values as the corresponding lower case ones but are not summed. This summation convention has been used in Mura (1987). For example, in the monomial $a_1 a_i b_i$, the repeated index is $i$ and the upper case index, $l$, takes the same value as $i$; it gives: $a_1 a_i b_l = a_1^2 b_1 + a_2^2 b_2 + a_3^2 b_3$. In (A.2), variable $\xi$ is the higher-valued solution of the equation:

$$ \frac{x_1^2}{c_1^2 + \xi} + \frac{x_2^2}{c_2^2 + \xi} + \frac{x_3^2}{c_3^2 + \xi} = 1 \qquad (A.3) $$

Functions $I_i(\xi)$ and $I_{ij}(\xi)$ are defined by the following elliptic integrals:

$$ \begin{aligned} I_i(\xi) &= 2\pi c_1 c_2 c_3 \int_{\xi}^{+\infty} \frac{ds}{(c_1^2 + s)\Delta(s)} \\ I_{ij}(\xi) &= 2\pi c_1 c_2 c_3 \int_{\xi}^{+\infty} \frac{ds}{(c_1^2 + s)(c_2^2 + s)\Delta(s)} \\ \Delta(s) &= [(c_1^2 + s)(c_2^2 + s)(c_3^2 + s)]^{1/2} \end{aligned} \qquad (A.4) $$

The derivative of $c_K^2 I_{KL}(\underline{\xi}) - I_L(\underline{\xi})$ with respect to $x_i$ and which appears in (A.2), reads:

$$ (c_K^2 I_{KL}(\underline{\xi}) - I_L(\underline{\xi}))_i = 4\pi c_1 c_2 c_3 \frac{\underline{\xi}}{\Delta(\underline{\xi})} \frac{\phi_K \phi_L \phi_I x_I}{x_P x_P \phi_P^2} \qquad (A.5) $$

where $\phi_K$ is given by:

$$ \phi_K = \frac{1}{c_K^2 + \xi} \qquad (A.6) $$

Consider now the case of a spheroidal inclusion, for which $c_1 = c_2 = b_1$, $c_3 = a_1$ embedded in an incompressible infinite medium ($v_0 = 1/2$). The higher-valued solution of (A.3) is $\xi = a^2 - a_1^2 = b^2 - b_1^2$. Then $\Delta(\xi) = ab^2$ and functions $I_i(\xi)$ and $I_{ij}(\xi)$ read:

$$ \begin{aligned} I_1(\xi) &= I_2(\xi) = 4\pi \frac{a_1 b_1^2}{ab^2} \frac{1 - \alpha}{2}; & I_3(\xi) &= 4\pi \frac{a_1 b_1^2}{ab^2} \alpha \\ I_{11}(\xi) &= I_{12}(\xi) = I_{22}(\xi) = 4\pi \frac{a_1 b_1^2}{8ab^4} (3 - 3\alpha - \beta) \\ I_{13}(\xi) &= I_{23}(\xi) = 4\pi \frac{a_1 b_1^2}{2a^3 b^2} \beta; & I_{33}(\xi) &= 4\pi \frac{a_1 b_1^2}{3a^3 b^2} (1 - \beta) \end{aligned} \qquad (A.7) $$

Taking into account that $\phi_1 = \phi_2 = 1/b^2$ and $\phi_3 = 1/a^2$, one has:

$$ x_i \phi_I = \left[ \frac{L_\lambda}{ab} e_\lambda \right]_i \qquad (A.8) $$

Using relations (A.5)–(A.8), and after some algebraic manipulations, expression (A.2) can be put into the form (12). The following decomposition of the Eshelby velocity field is now used:

$$ \underline{v}^F = \sum_{r=1}^{r=6} d_r^t \underline{v}_r^F \qquad (A.9) $$

where the $d_r^t$ for $r = 1, \dots, 6$ has been introduced in (16) and the $\underline{v}_r^F$ are defined by:

$$ \begin{cases} \underline{v}_{1\lambda}^F = \frac{a_1 b_1^2}{bL_\lambda} \left[ 1 + \frac{1}{2}(1 - 3\alpha)(1 - 3\cos^2(\varphi)) \right] \\ \underline{v}_{1\varphi}^F = -\frac{3a_1 b_1^2 (a_1^2 - b_1^2)}{4ab^2 L_\lambda} (1 - \alpha - \beta) \sin(2\varphi) \\ \underline{v}_{1\theta}^F = 0 \\ \\ \underline{v}_{2\lambda}^F = \frac{a_1 b_1^2}{bL_\lambda} \left[ 1 - \frac{2a_1^2 + b_1^2}{2a^2} \beta \right] (1 - 3\cos^2(\varphi)) \\ \underline{v}_{2\varphi}^F = \frac{3a_1 b_1^2 (2a_1^2 + b_1^2)}{4ab^2 L_\lambda} (1 - \alpha - \beta) \sin(2\varphi) \\ \underline{v}_{2\theta}^F = 0 \end{cases} \qquad (A.10) $$
$$
\left\{
\begin{aligned}
v_{3\lambda}^E &= -\frac{a_1 b_1^2}{b L_\lambda} \left[ 1 - \frac{b_1^2}{4b^2} (1 + 3\alpha + \beta) \right] \sin^2(\varphi) \cos(2\theta) \\
v_{3\phi}^E &= -\frac{a_1 b_1^4}{8ab^2 L_\lambda} (3 - 3\alpha - \beta) \sin(2\varphi) \cos(2\theta) \\
v_{3\theta}^E &= \frac{a_1 b_1^4}{4ab^3} (3 - 3\alpha - \beta) \sin(\varphi) \sin(2\theta) \\
v_{4\lambda}^E &= \frac{a_1 b_1^2}{b L_\lambda} \left[ 1 - \frac{b_1^2}{4b^2} (1 + 3\alpha + \beta) \right] \sin^2(\varphi) \sin(2\theta) \\
v_{4\phi}^E &= \frac{a_1 b_1^4}{8ab^2 L_\lambda} (3 - 3\alpha - \beta) \sin(2\varphi) \sin(2\theta) \\
v_{4\theta}^E &= \frac{a_1 b_1^4}{4ab^3} (3 - 3\alpha - \beta) \sin(\varphi) \cos(2\theta)
\end{aligned}
\right.
\qquad (A.11)
$$
$$
\left\{
\begin{aligned}
v_{5\lambda}^E &= \frac{a_1 b_1^2}{4ab^2 L_\lambda} \left[ (3\alpha + 1)(1 - k)a^2 + 3(1 - \alpha)(1 + k)b^2 \right] \sin(2\varphi) \cos(\theta) \\
v_{5\phi}^E &= -\frac{a_1 b_1^2}{b L_\lambda} \frac{1}{2}(1 - 3\alpha)[1 - k\cos(2\varphi)]\cos(\theta) \\
v_{5\theta}^E &= \frac{a_1 b_1^2}{2b^2}(1 - 3\alpha)(1 - k)\cos(\varphi)\sin(\theta) \\
v_{6\lambda}^E &= \frac{a_1 b_1^2}{4ab^2 L_\lambda} \left[ (3\alpha + 1)(1 - k)a^2 + 3(1 - \alpha)(1 + k)b^2 \right] \sin(2\varphi)\sin(\theta) \\
v_{6\phi}^E &= -\frac{a_1 b_1^2}{b L_\lambda} V_5(\lambda)[1 - k\cos(2\varphi)]\sin(\theta) \\
v_{6\theta}^E &= -\frac{a_1 b_1^2}{2b^2}(1 - 3\alpha)(1 - k)\cos(\varphi)\cos(\theta)
\end{aligned}
\right.
\qquad (A.12)
$$

with

$$k = \frac{a_1^2 + b_1^2}{a_1^2 - b_1^2} \qquad (A.13)$$

## Appendix B. An approximate expression of $(D, \bar{d}')$

For any function $F(\varphi, \theta)$, its mean value over an $iso - \lambda$ spheroid is defined by:

$$\langle F(\varphi, \theta) \rangle_{\mathcal{E}(\lambda)} = \frac{3}{4\pi(2a^2 + b^2)} \int_{\varphi=0}^{\varphi=\pi} \int_{\theta=0}^{\theta=2\pi} F(\varphi, \theta) L_{\lambda}^2 \times \sin(\varphi) d\varphi d\theta \qquad (B.1)$$

$F(\varphi, \theta)$ being assumed positive, the following inequality holds:

$$\langle F(\varphi, \theta) \rangle_{\mathcal{E}(\lambda)} \le \sqrt{\langle F^2(\varphi, \theta) \rangle_{\mathcal{E}(\lambda)}} \qquad (B.2)$$

To begin with, the following first approximation is then considered:

A1: $d_{eq}$ is replaced by its mean value along each spheroid confocal to the cavity, $[\langle d_{eq}^2 \rangle_{\mathcal{E}(\lambda)}]^{1/2}$.

This approximation has the advantage to preserve the upper bound character of the approach. It follows that the expression (29) of $\Gamma(D, \bar{d}')$ reads:

$$\Gamma(D, \bar{d}') = \frac{\sigma_0}{a_2 b_2^2} \int_{\lambda_1}^{\lambda_2} \left[ \langle d_{eq}^2 \rangle_{\mathcal{E}(\lambda)} \right]^{1/2} b(2a^2 + b^2) d\lambda \qquad (B.3)$$

A being uniform on $\Omega - \omega$, one has:

$$\langle d_{eq}^2 \rangle_{\mathcal{E}(\lambda)} = A_{eq}^2 + \frac{4}{3} \mathbf{A} : \langle \mathcal{D} \rangle_{\mathcal{E}(\lambda)} : \mathbf{d}^* + \frac{2}{3} \mathbf{d}^* : \langle \mathcal{D}^T : \mathcal{D} \rangle_{\mathcal{E}(\lambda)} : \mathbf{d}^* \qquad (B.4)$$

The next step of the calculation of $\Gamma(D, \bar{d}')$, is to perform the integral with respect to $\lambda$. In order to express the new criterion in a Gurson-like form, we aim to replace (B.3) by an integral of the form:

$$\Gamma(D, \bar{d}') = C \int_{u_1}^{u_2} \sqrt{Y^2 + X^2 u^2} \frac{du}{u^2} \qquad (B.5)$$

in which C, X and Y are three constants (independent of u) which will be determined later.

Due to the complex dependence of $\langle \mathcal{D} \rangle_{\mathcal{E}(\lambda)}$ and $\langle \mathcal{D}^T : \mathcal{D} \rangle_{\mathcal{E}(\lambda)}$ on the coordinate $\lambda$, a second approximation is needed. To this end, let us consider the term $\langle \mathcal{D}^T : \mathcal{D} \rangle_{\mathcal{E}(\lambda)}$ which appears in (B.4) in the case of a prolate void. In the limiting case of a hollow sphere ($b = a$), $\langle \mathcal{D}^T : \mathcal{D} \rangle_{\mathcal{E}(\lambda)}$ is slightly proportional to $(a_1/a)^6$, whereas for a hollow cylinder ($a = a_1 \to +\infty$), it is approximatively “proportional” to $(b_1/b)^4$. Let us introduce a variable x defined by:

$$x = \frac{a_1 b_1^2}{ab^2} \qquad (B.6)$$

The study of the two limiting cases of a prolate cavity shows that $\mathcal{D}$ is slightly proportional to $x^2$. Consider now the crossed term $\langle \mathcal{D} \rangle_{\mathcal{E}(\lambda)}$ in expression (B.4), still in the case of a prolate cavity. Tensor $\langle \mathcal{D} \rangle_{\mathcal{E}(\lambda)}$ is null for both the case of spherical and cylindrical cavities and, as in Gologanu et al. (1993), it can be neglected for simplicity. However, we propose here a better approximation by replacing $\langle \mathcal{D} \rangle_{\mathcal{E}(\lambda)}$ by a term proportional to $x^2$.

The case of oblate cavities is somewhat different. Although for low values of the eccentricity, the fourth order tensor $\langle \mathcal{D}^T : \mathcal{D} \rangle_{\mathcal{E}(\lambda)}$ is still proportional to $x^2$. When $e \to 1$ (namely in the “sandwich” case, consisting of two plane and parallel layers separated by an empty space), the components of tensor $\langle \mathcal{D}^T : \mathcal{D} \rangle_{\mathcal{E}(\lambda)}$ tend to a finite limit. The approximation considered in the case of a prolate cavity appears then inappropriate for oblate cavity since $x \to \infty$ when $e \to 1$. The idea proposed by Gologanu et al. (1994) and also used in Monchiet et al. (2006, 2007, 2008) is to introduce a new variable y defined as:

$$y = \frac{a_1 b_1^2}{\chi c^3 + ab^2} \qquad (B.7)$$

where $\chi$ is a constant whose value will be specified later. Note that y is proportional to x for low values of e, but has a finite value when $e \to 1$. When $x \to +\infty$ one has $y \to 1/\chi$. The second advantage is that the change of variable $\lambda \to y$ preserves the form of the element of integration, that is $dx/x^2 = dy/y^2$.

Let us now introduce the variable u defined by $u = x$ for prolate cavities and by $u = y$ for oblate cavities. We consider the following approximation:

A2 Terms $\langle \mathcal{D} \rangle_{\mathcal{E}(\lambda)}$ and $\langle \mathcal{D}^T : \mathcal{D} \rangle_{\mathcal{E}(\lambda)}$ which appear in the expression of $d_{eq}$ (see Eq. (B.4)) are replaced by terms proportional to $u^2 : \langle \mathcal{D} \rangle_{\mathcal{E}(\lambda)} = \mathcal{Q}u^2$ and $\langle \mathcal{D}^T : \mathcal{D} \rangle_{\mathcal{E}(\lambda)} = \mathcal{P}u^2$ (B.4) reads then:

$$d_{eq}^2 = A_{eq}^2 + \frac{2}{3} [2\mathbf{A} : \mathcal{Q} : \mathbf{d}^* + \mathbf{d}^* : \mathcal{P} : \mathbf{d}^*] u^2 \qquad (B.8)$$

It remains now to establish the required formula for $\mathcal{Q}$ and $\mathcal{P}$. Since the variations of the components of tensors
$\langle \mathbb{D} \rangle_{\mathcal{E}(\lambda)}/u^2$ and $\langle \mathbb{D}^T : \mathbb{D} \rangle_{\mathcal{E}(\lambda)}/u^2$ with $u$ are low, they are replaced by their means values over $[u_1, u_2]$:

$$ 
\begin{aligned}
\mathbb{Q} &= \frac{1}{u_2 - u_1} \int_{u_1}^{u_2} \langle \mathbb{D} \rangle_{\mathcal{E}(\lambda)} \frac{du}{u^2} = \frac{(1+g)(f+g)}{f^2(1-f)} \langle \mathbb{D} \rangle_{\Omega-\omega} \\
\mathbb{P} &= \frac{1}{u_2 - u_1} \int_{u_1}^{u_2} \langle \mathbb{D}^T : \mathbb{D} \rangle_{\mathcal{E}(\lambda)} \frac{du}{u^2} = \frac{(1+g)(f+g)}{f^2(1-f)} \langle \mathbb{D}^T : \mathbb{D} \rangle_{\Omega-\omega}
\end{aligned}
\quad (B.9) $$

where $g = 0$ for a prolate cavity and is given by:

$$ g = \chi \frac{c^3}{a_2 b_2^2} \quad (B.10) $$

for an oblate one. The expressions of $\mathbb{P}$ and $\mathbb{Q}$ are given in appendix C.

While the integration $\Gamma(\mathbf{D}, \bar{\mathbf{d}}')$ with (B.8) can a priori be done, expression (B.8) is not in a symmetric form according to $\mathbf{A}$ and $\mathbf{d}'$, due to the presence of the crossed term. A rewriting of (B.8) as a "symmetric" expression is needed not only for the estimate of the macroscopic plastic dissipation but also for calculation of the stresses at yielding (see subsection (4.1)). To this end, we express (B.8) into the following equivalent form:

$$ 
\begin{aligned}
d_{eq}^2 = & \frac{2}{3} \mathbf{A} : [\mathbb{I} - u^2 \mathbb{Q} : \mathbb{P}^{-1} : \mathbb{Q}^T] \\
& + \frac{1}{3} [\mathbf{d}^* + \mathbf{A} : \mathbb{Q} : \mathbb{P}^{-1}] : \mathbb{P} \\
& + [\mathbf{d}^* + \mathbb{P}^{-1} : \mathbb{Q}^T : \mathbf{A}] u^2
\end{aligned}
\quad (B.11) $$

Finally, in order to get the form (B.5) for $\Gamma(\mathbf{D}, \bar{\mathbf{d}}')$, the following last approximation is used:

A3 The term $u^2$ which multiplies $\mathbb{Q} : \mathbb{P}^{-1} : \mathbb{Q}^T$ in (B.11) is replaced by the constant determined as the mean value of $u^2$ over $u \in [u_1, u_2]$, i.e. $u_2 u_1$.

It is convenient to introduce the following definitions (using also the fact that $(1+f)(f+g)u_1 u_2 = f^2$):

$$ X^2 = \frac{2}{3} [\mathbf{d}^* + \mathbf{A} : \mathbb{Q} : \mathbb{P}^{-1}] : \mathbb{P} : [\mathbf{d}^* + \mathbf{A} : \mathbb{Q} : \mathbb{P}^{-1}] \quad (B.12) $$

$$ Y^2 = A_{eq}^2 - \frac{2f^2}{3(1+g)(g+f)} \mathbf{A} : \mathbb{Q} : \mathbb{P}^{-1} : \mathbb{Q}^T : \mathbf{A} $$

$\Gamma(\mathbf{D}, \bar{\mathbf{d}}')$ is then given by:

$$ 
\begin{aligned}
\Gamma(\mathbf{D}, \bar{\mathbf{d}}') &= -\sigma_0 f \int_{u_1}^{u_2} \{Y^2 + X^2 u^2\}^{1/2} \frac{du}{u^2} \\
&= \sigma_0 f \left[ X \operatorname{arcsinh} \left( \frac{uX}{Y} \right) - \frac{\sqrt{Y^2 + u^2 X^2}}{u} \right]_{u_1}^{u_2}
\end{aligned}
\quad (B.13) $$

## Appendix C. Expressions of $\mathbb{Q}$, $\mathbb{P}$ and $\mathbb{L}$

We now provide the expressions of tensors $\mathbb{Q}$ and $\mathbb{P}$ introduced in relations (B.9). We first note that:

$$ 
\begin{aligned}
\langle \mathbf{d}^E \rangle_{\Omega-\omega} &= \frac{1}{|\Omega-\omega|} \int_{\partial\Omega} \frac{v^F \otimes e_j}{|\Omega-\omega|} dS - \frac{1}{|\Omega-\omega|} \int_{\partial\omega} \frac{v^F \otimes e_\lambda}{|\Omega-\omega|} dS \\
&= f(\mathfrak{S}(e_2) - \mathfrak{S}(e_1)) : \mathbf{d}^*
\end{aligned}
\quad (C.1) $$

Consequently, tensor $\mathbb{Q}$, defined in (B.9), takes the form:

$$ \mathbb{Q} = \frac{(1+g)(f+g)}{(1-f)f} (\mathfrak{S}(e_2) - \mathfrak{S}(e_1)) \quad (C.2) $$

We now come to the determination of $\mathbb{P}$. We propose to compute the quantity $\langle \mathbf{d}^E : \mathbf{d}^E \rangle_{\Omega-\omega}$ for which it is recalled that $\mathbf{d}^E = \mathbb{D} : \mathbf{d}'$. The decomposition (16) is used for the eigenstrain rate tensor and leads to:

$$ \langle \mathbf{d}^E : \mathbf{d}^E \rangle_{\Omega-\omega} = \sum_{rs} p_{rs} d_r^* d_s^* \quad (C.3) $$

Coefficients $p_{rs}$ for $r, s = 1, \dots, 6$ has been computed with MAPLE software; their exact expressions obtained from the above relation read:

$$ p_{11} = 3 \frac{(1+g)(f+g)}{f(1-f)} \{(1+3\alpha_1)(1-\alpha_1) - f(1+3\alpha_2)(1-\alpha_2)\} $$

$$ p_{12} = 3 \frac{(1+g)(f+g)}{f(1-f)} \{(1-3\alpha_1)(1-\alpha_1 - \beta_1) - f(1-3\alpha_2)(1-\alpha_2 - \beta_2)\} $$

$$ p_{22} = 3 \frac{(1+g)(f+g)}{f(1-f)} \left\{ \begin{aligned} & (3\alpha_1 + 3\beta_1 - 1)(1 - \alpha_1 - \beta_1) \\ & - f(3\alpha_2 + 3\beta_2 - 1)(1 - \alpha_2 - \beta_2) \\ & + 3f \left(1 - \frac{a_1^2}{a_2^2}\right) \left(1 - \frac{b_1^2}{b_2^2}\right) \left[\alpha_2(1 - 3\beta_2) - (1 - \beta_2)^2\right] \end{aligned} \right\} $$

$$ p_{33} = p_{44} = \frac{1}{12} \frac{(1+g)(f+g)}{f(1-f)} \left\{ \begin{aligned} & (1+3\alpha_1+\beta_1)(3-3\alpha_1-\beta_1) \\ & -f(1+3\alpha_2+\beta_2)(3-3\alpha_2-\beta_2) \\ & +f\left(1-\frac{a_1^2}{a_2^2}\right)\left(1-\frac{b_1^2}{b_2^2}\right)[3(\alpha_2-1)(2-\beta_2)-\beta_2^2-12] \end{aligned} \right\} $$

$$ 
\begin{aligned}
p_{55} = p_{66} = & \frac{1}{3} \frac{(1+g)(f+g)}{f(1-f)} \{(1-3\alpha_1-2\beta_1)(3\alpha_1+2\beta_1-3) \\ & -f(1-3\alpha_2-2\beta_2)(3\alpha_2+2\beta_2-3) \\ & +2f\left(1-\frac{a_1^2}{a_2^2}\right)\left(1-\frac{b_1^2}{b_2^2}\right)[-6\alpha_2\beta_2+3\beta_2-2\beta_2^2+3\alpha_2-3] \end{aligned}
\quad (C.4) $$

In the expressions of coefficients $p_{22}$, $p_{33}$, $p_{44}$, $p_{55}$ and $p_{66}$ the last terms proportional to $(1 - a_1^2/a_2^2)(1 - b_1^2/b_2^2)$ vanishes in the case of a hollow cylinder ($a_1 = a_2$) and in the "sandwich" case ($b_1 = b_2$) and are neglected in all cases.

With this last approximation, $\mathbb{P}$ can be put in the tensorial form:

$$ \mathbb{P} \approx \frac{(1+g)(f+g)}{f(1-f)} [\mathfrak{S}^T(e_1) : \mathbb{L}(e_1) - f\mathfrak{S}^T(e_2) : \mathbb{L}(e_2)] \quad (C.5) $$

in which the components of tensor $\mathbb{L}(e)$ are given by:

$$ L_{1111} = L_{2222} = \frac{1}{8}(3\beta + 13\alpha + 3); \quad L_{1122} = L_{2211} = \frac{1}{8}(\beta + 7\alpha + 1) $$

$$ L_{1133} = L_{3311} = \frac{1}{2}(2 - \beta - 2\alpha); \quad L_{3333} = 1 + \beta - \alpha $$

$$ L_{1313} = L_{2323} = \frac{1}{4}(3 - 3\alpha - 2\beta); \quad L_{1212} = \frac{1}{8}(1 + \beta + 3\alpha) \quad (C.6) $$

## Appendix D. The minimization procedure

The macroscopic criterion is defined by Eqs. (34) with (36):
$$ \left\{ \begin{array}{l} \Sigma = f \frac{\partial X}{\partial \mathbf{A}} \Sigma_X + \frac{\partial Y}{\partial \mathbf{D}} \Sigma_Y \\ f \frac{\partial X}{\partial \mathbf{d}^*} \Sigma_X + \frac{\partial Y}{\partial \mathbf{d}^*} \Sigma_Y = 0 \end{array} \right. \qquad (D.1) $$

The first step of the demonstration consists to apply, in the above relations, the following change of variable: $(\mathbf{D}, \bar{\mathbf{d}}^*) \equiv (\mathbf{A}, \bar{\mathbf{d}}^*)$. Using the following identities:

$$ \begin{gathered} \frac{\partial \mathbf{A}}{\partial \mathbf{D}} = \mathbb{I} - \mathfrak{S}(e_2) : \mathbb{J}, \qquad \frac{\partial \mathbf{A}}{\partial \bar{\mathbf{d}}^*} = -f \mathfrak{S}(e_2) : \mathbb{K} \\ \frac{\partial \mathbf{d}^*}{\partial \mathbf{D}} = \frac{1}{f} \mathbb{J}, \qquad \frac{\partial \mathbf{d}^*}{\partial \bar{\mathbf{d}}^*} = \mathbb{K} \end{gathered} \qquad (D.2) $$

It follows that the two relations (D.1) reads:

$$ \left\{ \begin{array}{l} \Sigma = \left[ f \frac{\partial X}{\partial \mathbf{A}} \Sigma_X + \frac{\partial Y}{\partial \mathbf{A}} \Sigma_Y \right] : \mathbb{I} - \mathfrak{S}(e_2) : \mathbb{J} + \Sigma_X \frac{\partial X}{\partial \mathbf{d}^*} : \mathbb{J} \\ - \left[ f \frac{\partial X}{\partial \mathbf{A}} \Sigma_X + \frac{\partial Y}{\partial \mathbf{A}} \Sigma_Y \right] : \mathfrak{S}(e_2) : \mathbb{K} + \Sigma_X \frac{\partial X}{\partial \mathbf{d}^*} : \mathbb{K} = 0 \end{array} \right. \qquad (D.3) $$

The second step of the demonstration consists to compute $\Sigma : \mathfrak{S}(e_2)$ with the first relation in (D.3), it leads to (using the fact that $\mathbb{J} : \mathfrak{S}(e_2) = \mathbb{J}$):

$$ \Sigma : \mathfrak{S}(e_2) = \left[ f \frac{\partial X}{\partial \mathbf{A}} \Sigma_X + \frac{\partial Y}{\partial \mathbf{A}} \Sigma_Y \right] : \mathfrak{S}(e_2) : \mathbb{K} + \Sigma_X \frac{\partial X}{\partial \mathbf{d}^*} : \mathbb{J} \qquad (D.4) $$

Combining this last relation with the second in (D.3), gives:

$$ \Sigma : \mathfrak{S}(e_2) = \Sigma_X \frac{\partial X}{\partial \mathbf{d}^*} \qquad (D.5) $$

Using the definition of $X$, see (B.12), it appears that:

$$ \Sigma_X^2 = \frac{3}{2} \Sigma : \mathfrak{S}(e_2) : \mathbb{P}^{-1} : \mathfrak{S}^T(e_2) : \Sigma \qquad (D.6) $$

We now look for the expression of $\Sigma_Y$. The sum of the two equations in (D.3) gives:

$$ \Sigma = \left[ f \frac{\partial X}{\partial \mathbf{A}} \Sigma_X + \frac{\partial Y}{\partial \mathbf{A}} \Sigma_Y \right] : \left[ \mathbb{I} - \mathfrak{S}(e_2) \right] + \Sigma_X \frac{\partial X}{\partial \mathbf{d}^*} \qquad (D.7) $$

In the above expression, the last term is given by (D.5). It follows that:

$$ \Sigma : [\mathbb{I} - \mathfrak{S}(e_2)] = \left[ f \frac{\partial X}{\partial \mathbf{A}} \Sigma_X + \frac{\partial Y}{\partial \mathbf{A}} \Sigma_Y \right] : [\mathbb{I} - \mathfrak{S}(e_2)] \qquad (D.8) $$

Which imply that:

$$ \bar{\Sigma} = f \frac{\partial X}{\partial \mathbf{A}} \Sigma_X + \frac{\partial Y}{\partial \mathbf{A}} \Sigma_Y \qquad (D.9) $$

Using now the property:

$$ \frac{\partial X}{\partial \mathbf{A}} = \frac{\partial X}{\partial \mathbf{d}^*} : \mathbb{P}^{-1} : \mathbb{Q}^T \qquad (D.10) $$

It follows, from (D.5), that:

$$ \Sigma_X \frac{\partial X}{\partial \mathbf{A}} = \Sigma : \mathfrak{S}(e_2) : \mathbb{P}^{-1} : \mathbb{Q}^T \qquad (D.11) $$

A combination of Eqs. (D.9) and (D.11) gives:

$$ \Sigma : [\mathbb{K} - f \mathfrak{S}(e_2) : \mathbb{P}^{-1} : \mathbb{Q}^T] = \frac{\partial Y}{\partial \mathbf{A}} \Sigma_Y \qquad (D.12) $$

Using the definition of $Y$, given by (B.12), one finally obtains:

$$ \begin{aligned} \Sigma_Y^2 = \frac{3}{2} \Sigma : [\mathbb{K} - f \mathfrak{S}(e_2) : \mathbb{P}^{-1} : \mathbb{Q}^T] \\ : \left[ \mathbb{I} - \frac{f}{\mu} \mathbb{Q} : \mathbb{P}^{-1} : \mathbb{Q}^T \right]^{-1} \\ : [\mathbb{K} - f \mathbb{Q} : \mathbb{P}^{-1} : \mathfrak{S}^T(e_2)] : \Sigma \end{aligned} \qquad (D.13) $$

where $\mu = (1+g)(f+g)/f$. A simplification of the above expression is possible. To do that, let us first introduce the fourth order tensor $\tilde{\mathfrak{S}}$ defined by:

$$ \tilde{\mathfrak{S}} = \frac{\mathfrak{S}(e_1) - f \mathfrak{S}(e_2)}{1 - f} \qquad (D.14) $$

It can be shown that $\mu \mathfrak{S}(e_2) = \mathbb{Q} + \mu \tilde{\mathfrak{S}}$ where it is re-called that: $\mu = (1+g)(f+g)/f$. Consequently, one has:

$$ \begin{aligned} [\mathbb{K} - f \mathbb{Q} : \mathbb{P}^{-1} : \mathfrak{S}^T(e_2)] : \Sigma &= \left[ \mathbb{I} - \frac{f}{\mu} \mathbb{Q} : \mathbb{P}^{-1} : \mathbb{Q}^T \right] \\ &: \bar{\Sigma} - f \mathbb{Q} : \mathbb{P}^{-1} : \tilde{\mathfrak{S}}^T : \Sigma \end{aligned} \qquad (D.15) $$

where the equality $\mathbb{Q}^T : \Sigma = \mathbb{Q}^T : \bar{\Sigma}$ has been used (since $\mathbb{Q}^T : \mathbb{J} = 0$). Using relation (D.15) in (D.13), one obtains for $\Sigma_Y$:

$$ \begin{aligned} \Sigma_Y^2 = \frac{3}{2} \Sigma : [\mathbb{K} - f \mathfrak{S}(e_2) : \mathbb{P}^{-1} : \mathbb{Q}^T] : \bar{\Sigma} &= \frac{3f}{2} \Sigma \\ : [\mathbb{Q} - f \mathfrak{S}(e_2) : \mathbb{P}^{-1} : \mathbb{Q}^T : \mathbb{Q}] : [\mathbb{P} - \frac{f}{\mu} \mathbb{Q}^T : \mathbb{Q}]^{-1} \\ &: \tilde{\mathfrak{S}}^T : \Sigma \end{aligned} \qquad (D.16) $$

Because $\mathbb{Q}$ can be decomposed into $\mathbb{Q} = \mu \mathfrak{S}(e_2) - \mu \tilde{\mathfrak{S}}$ we have:

$$ \mathbb{Q} - f \mathfrak{S}(e_2) : \mathbb{P}^{-1} : \mathbb{Q}^T : \mathbb{Q} = \mu \mathfrak{S}(e_2) : \mathbb{P}^{-1} : \left[ \mathbb{P} - \frac{f}{\mu} \mathbb{Q}^T : \mathbb{Q} \right] - \mu \tilde{\mathfrak{S}} \qquad (D.17) $$

Using the above relation in expression (D.16), we find:

$$ \begin{aligned} \Sigma_Y^2 = \Sigma_{eq} - \frac{3\mu f}{2} \Sigma : \mathfrak{S}(e_2) : \mathbb{P}^{-1} : \mathfrak{S}^T(e_2) : \Sigma + \frac{3\mu f}{2} \Sigma \\ : \tilde{\mathfrak{S}} : \left[ \mathbb{P} - \frac{f}{\mu} \mathbb{Q}^T : \mathbb{Q} \right]^{-1} : \tilde{\mathfrak{S}}^T : \Sigma \end{aligned} \qquad (D.18) $$

Finally, we use the property $\mathbb{P} - \frac{f}{\mu} \mathbb{Q}^T : \mathbb{Q} = \mu \tilde{\mathfrak{S}}^T : \tilde{\mathbb{I}}$ to obtain expression (39).

## References

* Benzerga, A., Besson, J., 2001. Plastic potentials for anisotropic porous solids. Eur. J. Mech. A. Solids 20, 397-434.

* Keralavarma, S.M., Benzerga, A.A., 2010. A constitutive model for plastically anisotropic solids with non-spherical voids. J. Mech. Phys. Solids 58, 874-901.

* Benzerga, A., Leblond, J.B., 2010. Ductile fracture by void growth to coalescence. Adv. Appl. Mech. 44, 169-305.

* Bristow, J.R., 1960. Microcracks and the static and dynamic elastic constants of annealed heavily cold-worked metals. Br. J. Appl. Phys. 11, 81-85.

* Budiansky, B., O'Connell, R.J., 1976. Elastic moduli of a cracked solid. Int. J. Solids Struct. 12 (2), 81-97.

* Budiansky, B., Hutchinson, J.W., Slutsky, S., 1982. Void growth and collapse in viscous solids. In: Hopkins, H.G., Sewell, M.J. (Eds.), Mechanics of Solids. Pergamon Press, Oxford, pp. 13-45.
Cazacu, O., Stewart, J.B., 2009. Plastic potentials for porous aggregates with the matrix exhibiting tension-compression asymmetry. J. Mech. Phys. Solids 57, 325-341.

Danas, K., Idiart, M.I., Ponte-Castañeda, P., 2008. A homogenization-based constitutive model for isotropic viscoplastic porous media. Int. J. Solids Struct. 45, 3392-3409.

Dormieux, L., Kondo, D., Ulm, F.-J., 2006. Microporomechanics. Wiley.

Eshelby, J.D., 1959. The elastic field outside an ellipsoidal inclusion. Proc. R. Soc. Lond. A 252, 561-569.

**Gârâjeu, M., 1995. Contribution à l'étude du comportement non linéaire de milieux poreux avec ou sans renfort. Thèse de doctorat de l'Université de Marseille.**

Garajeu, M., Michel, J.C., Suquet, P., 2000. A micromechanical approach of damage in viscoplastic materials by evolution in size, shape and distribution of voids. Comput. Methods Appl. Mech. Eng. 183, 223-246.

Gologanu, M., Leblond, J.-B., Devaux, J., 1993. Approximate models for ductile metals containing non-spherical voids - case of axisymmetric prolate ellipsoidal cavities. J. Mech. Phys. Solids 41 (11), 1723-1754.

Gologanu, M., Leblond, J.B., Perrin, G., Devaux, J., 1994. Approximate models for ductile metals containing non-spherical voids - case of axisymmetric oblate ellipsoidal cavities. J. Eng. Mater. Technol. 116, 290-297.

**Gologanu, M., 1997. Etude de quelques problèmes de rupture ductile des métaux. Thèse de doctorat de l'Université de Paris VI.**

Gologanu, M., Leblond, J.B., Perrin, G., Devaux, J., 1997. Recent extensions of Gurson's model for porous ductile metals. In: Suquet, P. (Ed.), Continuum Micromechanics. Springer Verlag.

Guo, T.F., Faleskog, J., Shih, C.F., 2008. Continuum modeling of a porous solid with pressure sensitive dilatant matrix. J. Mech. Phys. Solids 56, 2188-2212.

Gurson, A.L., 1977. Continuum theory of ductile rupture by void nucleation and growth. Part I: Yield criterion and flow rules for porous ductile media. J. Eng. Mater. Technol. 99, 2-15.

Kachanov, M., Shafiro, B., Tsukrov, I., 2003. Handbook of Elasticity Solutions. Kluwer Academic Publisher.

Kailasam, M., Ponte Casta neda, P., Willis, J., 1997. The effect of particle size, shape, distribution and their evolution on the constitutive response of nonlinear composites - I. Theory. Proc. R. Soc. 355, 1835-1852.

Lee, B.J., Mear, M.E., 1992. Axisymmetric deformation of power-law solids containing a dilute concentration of aligned spheroidal voids. J. Mech. Phys. Solids 40, 1805-1836.

Lin, J., Kanit, T., Monchiet, V., Shao, J.F., Kondo, D., 2010. Numerical implementation of a recent improved Gurson-type model and application to ductile fracture. Comput. Mater. Sci. 47, 901-906.

McClintock, F.A., 1968. A criterion for ductile fracture by the growth of holes. ASME J. Appl. Mech. 35, 363-371.

**Monchiet, V., 2006. Contributions à la modélisation micromécanique de l'endommagement et de la fatigue des métaux ductiles. Ph.D. Thesis, University of Sciences and Technologies of Lille.**

Monchiet, V., Gruescu, C., Charkaluk, E., Kondo, D., 2006. Approximate yield criteria for anisotropic metals with prolate or oblate voids. C.R. Mécanique 334 (7), 431-439.

Monchiet, V., Charkaluk, E., Kondo, D., 2007. An improvement of Gurson type models of porous materials by using Eshelby-like trial velocity fields. C.R. Mécanique 335 (1), 32-41.

Monchiet, V., Cazacu, O., Charkaluk, E., Kondo, D., 2008. Approximate criteria for anisotropic metals containing non spherical voids. Int. J. Plast. 24, 1158-1189.

Monchiet, V., Charkaluk, E., Kondo, D., 2011. A micromechanics-based modification of the Gurson criterion by using Eshelby-like velocity fields. Eur. J. Mech. A. Solids 30, 940-949.

Mura, T., 1987. Micromechanics of Defects in Solids. Martinus Nijhoff, Dordrecht.

**Nemat Nasser, S., Hori, M., 1999. Micromechanics: Overall Properties of Heterogeneous Materials. Elsevier, The Netherlands.**

Pastor, F., Kondo, D., Pastor, J., 2011. Numerical limit analysis bounds for ductile porous media with oblate voids. Mech. Res. Commun. 38, 350-354.

Pastor, F., Kondo, D., 2013. Assessment of hollow spheroid models for ductile failure prediction by limit analysis and conic programming. Eur. J. Mech. A. Solids 38, 100-114.

Ponte-Castañeda, P., 1991. The effective mechanical properties of nonlinear isotropic composites. J. Mech. Phys. Solids 39, 45-71.

Ponte-Castañeda, P., Zaidman, M., 1994. Constitutive models for porous materials with evolving microstructure. J. Mech. Phys. Solids 42, 1459-1492.

Ponte-Castañeda, P., Suquet, P., 1998. Nonlinear composites. Adv. Appl. Mech 34, 171-302.

Qiu, Y.P., Weng, G.J., 1993. Plastic potential and yield function of porous materials with aligned and randomly oriented spheroidal voids. Int. J. Plast. 9, 271-290.

Rice, J.R., Tracey, D.M., 1969. On a ductile enlargement of voids in triaxial stress fields. J. Mech. Phys. Solids 17, 201-217.

Shen, W.Q., Shao, J.F., Dormieux, L., Kondo, D., 2012. Approximate criteria for ductile porous materials having a Green type matrix: application to double porous media. Comput. Mater. Sci. 62, 189-194.

**Suquet, P., 1982. Plasticité et homogénéisation. Thèse de doctorat d'état, Université Pierre et Marie Curie, Paris VI.**

Tvergaard, V., 1981. Influence of voids on shear band instabilities under plane strain conditions. Int. J. Fract. 17, 389-407.

Tvergaard, V., Needleman, A., 1984. Analysis of cup-cone fracture in round tensile bar. Acta. Metall. 32, 157-169.

Tvergaard, V., 1990. Material failure by void growth and coalescence. Adv. Appl. Mech. 27, 83-151.