# An improvement of Gurson-type models of porous materials by using Eshelby-like trial velocity fields

Vincent Monchiet, Eric Charkaluk, Djimedo Kondo

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

## ► To cite this version:

Vincent Monchiet, Eric Charkaluk, Djimedo Kondo. An improvement of Gurson-type models of porous materials by using Eshelby-like trial velocity fields. Comptes Rendus. Mécanique, 2007, 335 (1), pp.32-41. 10.1016/j.crme.2006.12.002. hal-00136688

HAL Id: hal-00136688

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

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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# An improvement of Gurson-type models of porous materials by using Eshelby-like trial velocity fields

Vincent Monchiet *, Eric Charkaluk, Djimedo Kondo

Laboratoire de mécanique de Lille-UMR CNRS 8107, université de sciences et technologies de Lille, cité scientifique,
59655 Villeneuve d'Ascq cedex, France

## Abstract

New expressions of the macroscopic criteria of perfectly plastic rigid matrix containing prolate and oblate cavities are presented. The proposed approach, derived in the framework of limit analysis, consists in the consideration of Eshelby-like trial velocity fields for the determination of the macroscopic dissipation. It is shown that the obtained results significantly improve existing criteria for ductile porous media. Moreover, for low stress triaxialities, these new results also agree perfectly with the (nonlinear) Hashin–Shtrikhman bound established by Ponte-Castañeda and Suquet.

## Résumé

Une amélioration des modèles de type Gurson pour les mileux poreux par utilisation des champs tests d'Eshelby. On présente de nouvelles expressions du critère macroscopique de milieux poreux constitués d'une matrice rigide parfaitement plastique contenant des cavités allongées ou aplaties. L'approche proposée, formulée dans le cadre de l'analyse limite, repose sur la considération de champs test de vitesse de type Eshelby pour la détermination de la dissipation macroscopique. On démontre que les résultats obtenus améliorent de manière significative les critères de milieux poreux ductiles existants. De plus, pour les faibles triaxialités de contrainte, ces nouveaux résultats s'accordent aussi parfaitement avec les bornes (non linéaires) d'Hashin–Shtrikhman établies par Ponte-Castañeda et par Suquet.

**Keywords:** Computational solid mechanics; Ductile porous metals; Generalized Gurson criterion; Prolate and oblate voids; Eshelby velocity fields; Limit analysis

**Mots-clés :** Mécanique des solides numériques ; Métaux ductiles poreux ; Critère de Gurson généralisé ; Cavités allongées et aplaties ; Champ de vitesse d'Eshelby ; Analyse limite

<footnote>* Corresponding author.
E-mail addresses: vincent.monchiet@ed.univ-lille1.fr (V. Monchiet), eric.charkaluk@univ-lille1.fr (E. Charkaluk), kondo@univ-lille1.fr (D. Kondo).</footnote>
# 1. Introduction

Since three decades, the modelling of the behavior of ductile porous media has been the subject of important researches in non linear mechanics of materials. In his pioneering work, Gurson [1] developed a limit analysis approach of a hollow sphere. The plastic matrix is assumed to obey the von Mises criterion:

$$f(\sigma) = \sigma_{\text{eq}} - \sigma_0 \le 0; \quad \text{with: } \sigma_{\text{eq}} = \sqrt{\frac{3}{2}\bar{\sigma} : \bar{\sigma}} \qquad (1)$$

where $\sigma$ denotes the microscopic stress field, $\bar{\sigma}$ its deviatoric part and $\sigma_{\text{eq}}$ the microscopic von Mises equivalent stress. $\sigma_0$ represents the yield stress in tension. More specifically, Gurson obtained a macroscopic criterion which, in the case of spherical voids, reads:

$$\frac{\Sigma_{\text{eq}}^2}{\sigma_0^2} + 2f \cosh\left\{\frac{3\Sigma_h}{2\sigma_0}\right\} - 1 - f^2 = 0 \qquad (2)$$

$\Sigma_h$ denotes the hydrostatic stress and $f$ the porosity. $\Sigma_{\text{eq}}$ is the macroscopic von Mises equivalent stress. It has been demonstrated (see for instance [2]), that the yield surface defined by (2) constitutes an upper bound for Hashin's well-known composite spheres assemblage and gives the exact result for purely hydrostatic macroscopic loading. Later, using variational techniques, Ponte-Castañeda [3] and Suquet [4] obtained a rigorous nonlinear Hashin-Shtrikhman upper bound which, for spherical voids, takes the following form<fnref n="1" />:

$$(1 + \frac{2}{3}f)\frac{\Sigma_{\text{eq}}^2}{\sigma_0^2} + \frac{9}{4}f\Sigma_h^2 - (1 - f)^2 = 0 \qquad (3)$$

An important observation is that the Gurson model (Eq. (2)) violates this upper bound for low values of the stress triaxiality $T = \Sigma_h / \Sigma_{\text{eq}}$. However, its predominance over the Hashin-Shtrikhman bound is still observed for high stress triaxialities. A possible method to improve the predictions of the original limit analysis approach of Gurson and in fact the subsequent models consists in considering refined trial velocity fields (see for instance [6] and [7]). Still, due to the limitation of trial velocity fields which have been explored, few significant improvements have been obtained in the past studies.

Therefore, the main objective of the present Note is to develop a limit analysis approach based on Eshelby-like velocity fields and to derive new expression of the yield function. The calculations will be performed in the general case of prolate and oblate voids (Sections 2 and 3). In this way, it is expected that the new approach will also provide an improved version and a generalization of the results obtained by [8] and [9]. In order to provide simple illustrations of the obtained approximate criterion, some specific cases such as spherical or cylindrical voids will be examined (Section 4).

# 2. Basic concepts and methodology

## 2.1. The studied cell

As in [6], consider a spheroidal (axisymmetric) prolate or oblate cavity with semi-axes $a_1$ (along $\mathbf{e}_3$), and $b_1$ (along $\mathbf{e}_1$ and $\mathbf{e}_2$) embedded in a cell which has the shape of a confocal spheroid with the semi-axes $a_2$ (along $\mathbf{e}_3$), and $b_2$ (along $\mathbf{e}_1$ and $\mathbf{e}_2$). $a_1 > b_1$ corresponds to a prolate cavity while $b_1 > a_1$ is associated to an oblate one. Let us denote $c$ the focal distance and $e_1$ the eccentricity defined by:

$$c = \sqrt{a_1^2 - b_1^2}; \quad e_1 = \frac{c}{a_1} \quad (\text{prolate}); \quad c = \sqrt{b_1^2 - a_1^2}; \quad e_1 = \frac{c}{b_1} \quad (\text{oblate}) \qquad (4)$$

It is convenient to introduce the spheroidal coordinates characterized by $\lambda, \beta, \theta$, and defined in the cylindrical frame (coordinates $\rho, \theta, z$) by $\rho = b \sin \beta, z = a \cos \beta$. The iso-$\lambda$ surface defines confocal spheroids, with semi-axes $a =$

<footnote n="1">¹ Similar results for fluid saturated porous media with a Drucker-Prager matrix can be found in [5].</footnote>
$c \cosh(\lambda)$, $b = c \sinh(\lambda)$ and eccentricity $e = c/a$, for prolate. An oblate spheroid is associated to semi-axes $a = c \sinh(\lambda)$, $b = c \cosh(\lambda)$ and eccentricity $c/b$. The unit vectors of the new base are:

$$ \underline{e}_{\lambda} = \frac{1}{L_{\lambda}} \{a \sin(\beta) \underline{e}_{\rho} + b \cos(\beta) \underline{e}_{z}\}; \quad \underline{e}_{\beta} = \frac{1}{L_{\lambda}} \{b \cos(\beta) \underline{e}_{\rho} - a \sin(\beta) \underline{e}_{z}\}; \quad \underline{e}_{\theta} = \underline{e}_{\theta} \quad (5) $$

with $L_{\lambda} = \sqrt{a^2 \sin^2(\beta) + b^2 \cos^2(\beta)}$, $\theta \in [0, 2\pi]$, $\beta \in [0, \pi]$ and $\underline{e}_{\rho} = \cos(\theta)\underline{e}_1 + \sin(\theta)\underline{e}_2$.

## 2.2. The Eshelby-like velocity field

The velocity field, $\underline{v}$, in the matrix, is classically decomposed into a uniform velocity field, $A \cdot \underline{x}$, and an heterogeneous velocity field, $\underline{v}^E$, as follows: $\underline{v} = A \cdot \underline{x} + \underline{v}^E$. For $\underline{v}^E$, we consider the exterior point Eshelby solution (see [10]) adapted here to an incompressible viscous fluid containing a spheroidal inclusion. For convenience, this solution (see also [11]) can be put in the form:

$$ \underline{v}^E = \sum_{r=1}^{r=6} \underline{v}^r d_r^* \quad (6) $$

where $d_r^* = \underline{d}^* : \underline{Q}_r / (\underline{Q}_r : \underline{Q}_r)$. Tensor $\underline{d}^*$ is an eigenstrain in the inhomogeneity and $\underline{Q}_r$ are defined by:

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

with the property $\underline{Q}_r : \underline{Q}_s = 0$ if $r \neq s$. The quantities $d_r^*$ are then related to the components of $\underline{d}^*$ by:

$$ \begin{aligned} d_1^* &= \frac{1}{3}(d_{11}^* + d_{22}^* + d_{33}^*); & d_2^* &= \frac{1}{3}\left[\frac{d_{11}^* + d_{22}^*}{2} - d_{33}^*\right] \\ d_3^* &= \frac{d_{22}^* - d_{11}^*}{2}; & d_4^* &= d_{12}^*; & d_5^* &= d_{13}^*; & d_6^* &= d_{23}^* \end{aligned} \quad (8) $$

The velocity fields $\underline{v}^r$ for $r = 1, 6$ are given in Appendix A, in the spheroidal frame.

Note that $\underline{v}^1$ and $\underline{v}^2$ are two axisymmetric velocity fields; they are independent of the coordinate $\theta$ and the component $v_{\theta} = 0$. The first field $\underline{v}^1$ is the one used by [8] and by [9] for the determination of the macroscopic yield function of a plastic matrix containing prolate and oblate voids respectively. Due to some limitations of $\underline{v}^1$, it has been proposed in [6] to incorporate a supplementary velocity field which does not appear in the Eshelby velocity field.

Garajeu and Suquet [7] have also derived an expression for the macroscopic yield function of the porous material in the case of prolate cavities by using a truncated expression of the velocity fields $\underline{v}^1$ and $\underline{v}^2$ which corresponds to $f_1(\lambda) = g_1(\lambda) = 0$ (see Appendix A). It is also interesting to mention that $\underline{v}^1$ and $\underline{v}^2$ are contained in the general axisymmetric velocity fields proposed by Lee and Mear [12]. However, the other fields $\underline{v}_r$ for $r = 3, 6$, considered in this study, are nonaxisymmetric and are not contained in the Lee and Mear fields.

The microscopic plastic strain rate is defined by:

$$ \underline{d} = A + \sum_{r=1}^{r=6} \underline{d}^r d_r^* $$

in which $\underline{d}^r = \nabla_s \underline{v}^r$. The microscopic dissipation reads $\pi(\underline{d}) = \sigma_0 d_{eq}$ where $d_{eq} = \sqrt{\frac{2}{3} \bar{\underline{d}} : \bar{\underline{d}}}$ is the equivalent plastic strain rate defined by:

$$ d_{\text{eq}}^2 = A_{\text{eq}}^2 + 2 \sum_{r=1}^{r=6} d_r^* A : \underline{d}^r + \sum_{r=1}^{r=6} \sum_{s=1}^{s=6} d_r^* d_s^* \underline{d}^r : \underline{d}^s \quad (9) $$

Since the Eshelby-like velocity field introduced by (6) does not comply with the uniform strain on the cell boundary, we propose to use the classical average rule which relates the macroscopic strain rate tensor $\underline{D}$ to $\underline{d}$. This allows us to relate $\underline{d}^*$ to $\underline{D}$ by:

$$ \underline{D} = \frac{1}{|\Omega|} \int_{\Omega} \underline{d} \, dV = A + f \mathbb{S}(e_2) : \underline{d}^* \quad (10) $$
where $|\Omega| = 4\pi a_2 b_2^2/3$ denotes the volume of the studied cell (matrix + void) and $\mathcal{S}(e)$ is the Eshelby tensor corresponding to a spheroidal cavity of eccentricity $e$ embedded in an incompressible medium; the components of $\mathcal{S}(e)$ are expressed in Appendix B. By combining then (B.1) and (10), one obtains:

$$ 
\begin{aligned}
d_1^* &= \frac{D_h}{f}; & A_{33} &= D_{33}(1-\alpha_2) - (D_{11} + D_{22})\alpha_2 + \varepsilon f \frac{2\alpha_2 a_2^2 + (\alpha_2 - 1)b_2^2}{2c^2} d_2^* \\
A_{22} - A_{11} &= D_{22} - D_{11} - \varepsilon f \frac{2a_2^2 + 3b_2^2(\alpha_2 - 1)}{4c^2} d_3^*; & A_{12} &= D_{12} - \varepsilon f \frac{2a_2^2 + 3b_2^2(\alpha_2 - 1)}{8c^2} d_4^* \\
A_{13} &= D_{13} - \varepsilon f \frac{a_2^2 + b_2^2}{4c^2}(1 - 3\alpha_2)d_5^*; & A_{23} &= D_{23} - \varepsilon f \frac{a_2^2 + b_2^2}{4c^2}(1 - 3\alpha_2)d_6^*
\end{aligned}
\quad (11) $$

in which $D_h \mathbf{1}$ represents the spherical part of $\mathbf{D}$ and the convention $\varepsilon = 1$ for prolate voids and $\varepsilon = -1$ for oblate voids is adopted. $\alpha_2 = \alpha(e_2)$ where $\alpha(e)$ is function of eccentricity $e$ and is given by (A.2).

It is readily seen that the complete velocity field is defined by 11 parameters, the components of $\mathbf{A}$ and the scalar $d_r^*$ for $r = 1, 6$. Condition (10), detailed in (11) together with (A.2), gives six relations between the different parameters. Still, their remains five unknown parameters which have to be determined.<fnref n="2" /> These are the $d_r^*$ for $r = 2, 6$ or equivalently the components of $\bar{\mathbf{d}}^*$ (the deviatoric part of $\bar{\mathbf{d}}^*$). It is interesting to note that, by putting $d_r^* = 0$ for $r = 2, 6$, and choosing $\mathbf{A}$ axisymmetric in (11), the homogeneous boundary strain rate conditions obtained by [8,9] are recovered.

### 2.3. The macroscopic dissipation $\Pi(D)$ and the minimization principle

Due to the dependence of the velocity field $\mathbf{d}$ with the following unknown kinematic parameters $\bar{\mathbf{d}}^*$, let us now introduce $\tilde{\Pi}(\mathbf{D}, \bar{\mathbf{d}}^*)$ defined by:

$$ \tilde{\Pi}(\mathbf{D}, \bar{\mathbf{d}}^*) = \frac{\sigma_0}{|\Omega|} \int_{\Omega-\omega} d_{\text{eq}} dV \quad (12) $$

where $\omega$ denotes the volume of the void. It is readily understood that the searched expression of the macroscopic dissipation, $\Pi(D)$, is derived from a minimization procedure of $\tilde{\Pi}(D, \bar{d}^*)$ with respect to $\bar{d}^*$:

$$ \Pi(D) = \min_{\bar{d}^*} [\tilde{\Pi}(D, \bar{d}^*)] \quad (13) $$

The yield surface, related to the macroscopic dissipation, is then assumed to be given by:

$$ \Sigma = \frac{\partial \Pi}{\partial D} \quad (14) $$

In fact such consideration can be interpreted as the condition of the necessary coherence of the energy definition at the two scales.

Due to the difficulty of integrating (12), some approximations are needed either for prolate and oblate voids in order to determine the macroscopic dissipation and then the macroscopic yield function.

### 3. Approximate expression of $\tilde{\Pi}(D, \bar{d}^*)$ and determination of the macroscopic yield locus

We now aim to derive the expression of $\tilde{\Pi}(D, d_r^*)$ by using the same type of approximations as in [13]. Let us first recall the expression of $\tilde{\Pi}(D, d_r^*)$ in the spheroidal frame.

$$ \tilde{\Pi}(D, \bar{d}^*) = \frac{\sigma_0}{|\Omega|} \int_{\lambda=\lambda_1}^{\lambda=\lambda_2} \int_{\beta=0}^{\beta=\pi} \int_{\theta=0}^{\theta=2\pi} d_{\text{eq}} b L_{\lambda}^2 \sin \beta \, d\lambda \, d\beta \, d\theta \quad (15) $$

The following approximations are made:

<footnote n="2">² Note that the components of $\mathbf{A}$ could also be chosen as unknowns but this does not change the final results.</footnote>
A1: $d_{eq}$ is replaced by a mean value $\{\langle d_{eq}^2 \rangle \epsilon\}^{1/2}$ on each confocal ellipsoid $\mathcal{E}$.

It follows that (15) takes the following form:

$$ \tilde{\Pi}(\mathbf{D}, \bar{\mathbf{d}}^*) = \frac{\sigma_0}{a_2 b_2^2} \int_{\lambda=\lambda_1}^{\lambda=\lambda_2} \{\langle d_{eq}^2 \rangle \epsilon\}^{1/2} b(2a^2 + b^2) d\lambda \quad (16) $$

where $\langle d_{eq}^2 \rangle \epsilon$ is such that:

$$ \langle d_{eq}^2 \rangle \epsilon = \frac{3}{4\pi(2a^2 + b^2)} \int_{\beta=0}^{\beta=\pi} \int_{\theta=0}^{\theta=2\pi} d_{eq}^2 L_{\lambda}^2 \sin(\beta) d\beta d\theta \quad (17) $$

Then, from the definition (9), one has:

$$ \langle d_{eq}^2 \rangle \epsilon = A_{eq}^2 + 2 \sum_{r=1}^{r=6} d_r^* \mathbf{A} : \mathbf{W}_r(e) + \sum_{r=1}^{r=6} \sum_{s=1}^{s=6} d_r^* d_s^* P_{rs}(e) \quad (18) $$

in which $\mathbf{W}_r(e) = \frac{2}{3} \langle \mathbf{d}^r \rangle \epsilon$ and $P_{rs}(e) = \frac{2}{3} \langle \mathbf{d}^r : \mathbf{d}^s \rangle \epsilon$. Note that $\mathbf{W}_r(e)$ and $P_{rs}(e)$ also depend on $e_1$. Let us now introduce the following variables $x$ and $y$ defined by:

$$ x = \frac{a_1 b_1^2}{ab^2}; \quad y = \frac{\chi a_1 b_1^2}{c^3 + \chi a_1 b^2} \quad (19) $$

with $\chi = \frac{3}{4}\sqrt{\pi^2 + \frac{32}{3}}$. As a unified notation, let us introduce $u$ such that $u = x$ for prolate cavities and $u = y$ for oblate cavities.

A2: Following [6], the expressions of $\mathbf{W}_r(e)$ and $P_{rs}(e)$ are replaced by $\mathbf{w}_r u^2$ and $p_{rs} u^2$ respectively where $\mathbf{w}_r$ and $p_{rs}$ are constant.

As in [13], parameters $p_{rs}$ and tensors $\mathbf{w}_r$ are respectively determined as the mean values of $P_{rs}(e)/u^2$ and $\mathbf{W}_r(e)/u^2$ along the interval $[u_1, u_2]$:

$$ p_{rs} = \frac{1}{u_2 - u_1} \int_{u_1}^{u_2} P_{rs}(e) \frac{du}{u^2}; \quad \mathbf{w}_r = \frac{1}{u_2 - u_1} \int_{u_1}^{u_2} \mathbf{W}_r(e) \frac{du}{u^2}; \quad \text{for } r = 1, 6 \quad (20) $$

$u_1 = 1$ for a prolate cavity and $u_1 = \chi a_1 b_1^2 / (c^3 + \chi a_1 b_1^2)$ for an oblate one. This approximation is consistent only if the obtained values for $p_{rs}$ and $\mathbf{w}_r$ show low variations according to $e_1$ and $e_2$. The explicit determination of $p_{rs}$ and $\mathbf{w}_r$ is performed by using Maple software; the resulting expressions are given in Appendix C. The quantity $\langle d_{eq}^2 \rangle \epsilon$ can be rewritten as follows:

$$ \langle d_{eq}^2 \rangle \epsilon = A_{eq}^2 + 2 \sum_{r=1}^{r=6} \mathbf{A} : \mathbf{w}_r d_r^* u^2 + \sum_{r=1}^{r=6} \sum_{s=1}^{s=6} p_{rs} d_r^* d_s^* u^2 = A_{eq}^2 + [2 \mathbf{A} : \mathbf{W} : \mathbf{d}^* + \mathbf{d}^* : \mathbf{P} : \mathbf{d}^*] u^2 \quad (21) $$

where $\mathbf{W}$ and $\mathbf{P}$ are two fourth order tensors defined in Appendix C. Eq. (21) can also be expressed as follows:

$$ \langle d_{eq}^2 \rangle \epsilon = A_{eq}^2 - \mathbf{A} : \mathbf{W} : \mathbf{P}^{-1} : \mathbf{W}^T : \mathbf{A} u^2 + [\mathbf{d}^* + \mathbf{A} : \mathbf{W} : \mathbf{P}^{-1}] : \mathbf{P} : [\mathbf{d}^* + \mathbf{P}^{-1} : \mathbf{W}^T : \mathbf{A}] u^2 \quad (22) $$

where $\mathbf{W}^T$ represent the transpose of $\mathbf{W}$ such that $W_{ijkl}^T = W_{klij}$.

A3: Term $\mathbf{A} : \mathbf{W} : \mathbf{P}^{-1} : \mathbf{W}^T : \mathbf{A} u^2$ is neglected.
Although this term is not neglected in [9], [6] and [13], its effect on the macroscopic criterion is generally weak<fnref n="3" /> and motivates the approximation A3. Let us then introduce the following change of variable:

$$B^2 = \mathbf{B} : \mathbf{P} : \mathbf{B}; \quad \text{with: } \mathbf{B} = \mathbf{d}^* + \mathbf{P}^{-1} : \mathbf{W}^T : \mathbf{A} \qquad (23)$$

$\tilde{\Pi}(\mathbf{D}, \tilde{\mathbf{d}}^*)$ (Eq. (16)) takes then the following form:

$$\tilde{\Pi}(\mathbf{D}, \tilde{\mathbf{d}}^*) = \sigma_0 f \int_{u_2}^{u_1} \left\{ A_{\text{eq}}^2 + B^2 u^2 \right\}^{1/2} \frac{du}{u^2} = \sigma_0 f \left[ B \operatorname{arcsinh} \left\{ \frac{uB}{A_{\text{eq}}} \right\} - \frac{\sqrt{A_{\text{eq}}^2 + u^2 B^2}}{u} \right]_{u_2}^{u_1} \qquad (24)$$

As previously indicated, the macroscopic yield locus is defined by (14) in which $\Pi(\mathbf{D})$ corresponds to the minimum of $\tilde{\Pi}(\mathbf{D}, \tilde{\mathbf{d}}^*)$ with respect to $\tilde{\mathbf{d}}^*, d_1^*$ being completely identified and related to $D_h$ (see (11)). The macroscopic criterion is then obtained by resolution of:

$$\boldsymbol{\Sigma} = \frac{\partial \tilde{\Pi}(\mathbf{D}, \tilde{\mathbf{d}}^*)}{\partial \mathbf{D}} \quad \text{with:} \quad \frac{\partial \tilde{\Pi}(\mathbf{D}, \tilde{\mathbf{d}}^*)}{\partial \tilde{\mathbf{d}}^*} = 0 \qquad (25)$$

An analytical solution to this system, given in Appendix D, leads to the following generalized approximate criterion corresponding to prolate and oblate voids<fnref n="4" />:

$$\frac{\tilde{\Sigma}_{\text{eq}}^2}{\sigma_0^2} + 2(1 + g)(f + g) \cosh\left\{\frac{\tilde{\Sigma}_P}{\sigma_0}\right\} - (1 + g)^2 - (f + g)^2 = 0 \qquad (26)$$

where $g$ is conventionally taken as 0 for a prolate void and is given by $g = c^2/(\chi a_2 b_2^2)$ for an oblate one. $\tilde{\Sigma}_P$ and $\tilde{\Sigma}_{\text{eq}}$ which enters in the criterion (26) are defined by:

$$\tilde{\Sigma}_P^2 = \boldsymbol{\Sigma} : \mathcal{S}(e_2) : \mathbf{P}^{-1} : \mathcal{S}^T(e_2) : \boldsymbol{\Sigma}; \quad \tilde{\Sigma}_{\text{eq}}^2 = \frac{3}{2} \tilde{\boldsymbol{\Sigma}} : \tilde{\boldsymbol{\Sigma}}; \quad \text{with: } \tilde{\boldsymbol{\Sigma}} = \boldsymbol{\Sigma} - \mathbf{W} : \mathbf{P}^{-1} : \mathcal{S}(e_2) : \boldsymbol{\Sigma} \qquad (27)$$

where the components of $\mathbf{P}^{-1}$ are given in Appendix C. Eqs. (26) and (27) constitute the most important results of the present Note. Their new features and advantage over the existing Gurson-type criteria of porous media (see [8] and [9]) come mainly from the terms introduced by the components of $\tilde{\mathbf{d}}^*$. In particular, it is interesting to point out that, due to $\tilde{\Sigma}_P$, the new criterion introduces a coupling between pure shear stresses and the porosity $f$ which does not exist in the existing criteria. Note also that the nonlinear Hashin–Shtrikman bound exhibits a coupling of this type (see [14]).

## 4. The particular cases of spherical and cylindrical voids

Let us consider, as a first particular illustration, the case of the spherical cavity; this is obtained by putting $e_1 \to 0$ and $e_2 \to 0$. It follows that $\mathbf{W} = 0$, and then $\tilde{\Sigma}_{\text{eq}} = \tilde{\Sigma}_{\text{eq}}$. Coefficients $p_{rs}$ are $p_{11} = 4$, $p_{12} = 0$, $p_{22} = 3p_{33} = 3p_{55} = 24/25$, and then $\tilde{\Sigma}_P^2 = 9\tilde{\Sigma}_h^2/4 + 2\tilde{\Sigma}_{\text{eq}}^2/3$. The resulting macroscopic yield locus<fnref n="5" /> is given by:

$$\frac{\tilde{\Sigma}_{\text{eq}}^2}{\sigma_0^2} + 2f \cosh\left\{\frac{1}{\sigma_0}\sqrt{\frac{9}{4}\tilde{\Sigma}_h^2 + \frac{2}{3}\tilde{\Sigma}_{\text{eq}}^2}\right\} - 1 - f^2 = 0 \qquad (28)$$

In the particular loading case $\tilde{\Sigma}_{\text{eq}} = 0$, the proposed criterion (28) reduces to the Gurson one which gives the exact solution under purely hydrostatic pressure. Eq. (28) provides also the strongly remarkable property that the deviatoric equivalent stress $\tilde{\Sigma}_{\text{eq}}$ enters with the mean stress $\tilde{\Sigma}_h$ in the cosh term which controls the cavity growth process.

Consider now low values of $\tilde{\Sigma}_h$; it can be verified that in this case $\cosh\{\frac{1}{\sigma_0}\sqrt{\frac{9}{4}\tilde{\Sigma}_h^2 + \frac{2}{3}\tilde{\Sigma}_{\text{eq}}^2}\} \simeq 1 + \frac{9}{8}\frac{\tilde{\Sigma}_h^2}{\sigma_0^2} + \frac{\tilde{\Sigma}_{\text{eq}}^2}{\sigma_0^2}$. Reporting this result into (28), one recovers the Hashin–Shtrikman bound (see Eq. (3)). The predictions of the yield surface by the different models, in the case of spherical voids and for two different porosities, are represented on Fig. 1 and confirm the good performance of the proposed model.

<footnote n="3">³ However, note that the consideration of these terms does not rise an important computational difficulty in the determination of the macroscopic yield locus but gives more complex expressions.</footnote>

<footnote n="4">⁴ Obviously, it will be verified that this result generalizes also the Gurson criterion.</footnote>

<footnote n="5">⁵ Note that this result extends to fluid saturated porous media by replacing $\tilde{\Sigma}_h$ by $\tilde{\Sigma}_h + P$ where $P$ is the fluid pressure.</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1. Yield surface of the porous material—case of spherical void: comparison of the proposed model (28), the Gurson yield locus and the Ponte-Castañeda model. (a) Porosity $f = 0.01$; (b) porosity $f = 0.3$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 1. Surface de charge des matériaux poreux—cas d'une cavité sphérique : comparaison du modèle proposé (28) avec la surface de Gurson et avec le modèle de Ponte-Castañeda. (a) Porosité $f = 0.01$; (b) porosité $f = 0.3$.</figcaption></figure>

It must be recalled that, unlike the studies [1] and [8,9], the results obtained in the present study are based on a nonkinematical velocity field defined by the condition (10); it is this choice which allows the connection with the nonlinear homogenization bounds [14].

Let us consider now, as a second illustration, the case of a cylindrical cavity; then $e_1 \to 1$, $e_2 \to 1$, $W = 0$, then $\tilde{\Sigma}_{eq} = \Sigma_{eq}$. Coefficients $p_{rs}$ are $p_{11} = 3$, $p_{12} = p_{22} = 0$ and $p_{rr} = 1/3$ for $r = 3, 6$, then $\tilde{\Sigma}_p^2 = (2\Sigma_{11}^2 + 2\Sigma_{22}^2 + 4\Sigma_{12}^2 + 4\Sigma_{13}^2 + 4\Sigma_{23}^2)/3$. The macroscopic yield surface is described by:

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

As for spherical voids, similar comments can be done for cylindrical cavities: the exact solution under purely hydrostatic pressure is recovered; for low values of $\Sigma_h/\sigma_0$ the criterion leads to the approximation:

$$ \frac{\Sigma_{eq}^2}{\sigma_0^2} + \frac{3f}{2} \frac{1}{\sigma_0^2} \{ \Sigma_{11}^2 + \Sigma_{22}^2 + 2\Sigma_{12}^2 + 2\Sigma_{13}^2 + 2\Sigma_{23}^2 \} - (1-f)^2 = 0 \quad (30) $$

which coincides with the result established by Suquet [4].

## Acknowledgements

The authors are grateful to Professor J.-B. Leblond for fruitful suggestions on Note [13] which are partly at the origin of the present study.

## Appendix A. Eshelby velocity field

The exterior point Eshelby solution associated to a spheroidal inclusion in an isotropic matrix, can be found in [10] or in the textbook [11]; the corresponding velocity field is given as a combination of the following velocities:

$$ \underline{v}^1 = \frac{a_1 b_1^2}{b L_\lambda} \left[ 1 + f_1(\lambda) \frac{1 - 3 \cos^2(\beta)}{2} \right] \underline{e}_\lambda + \frac{a_1 b_1^2}{b L_\lambda} g_1(\lambda) \sin(2\beta) \underline{e}_\beta \quad (A.1) $$

$$ \underline{v}^2 = \frac{a_1 b_1^2}{b L_\lambda} [1 - k_2 f_1(\lambda)] (1 - 3 \cos^2(\beta)) \underline{e}_\lambda - 2 \frac{a_1 b_1^2}{b L_\lambda} g_1(\lambda) k_2 \sin(2\beta) \underline{e}_\beta $$
with:

$$
\begin{align}
f_1(\lambda) &= 1 - 3\alpha; & g_1(\lambda) &= -\frac{3}{4ab}[2a^2\alpha + b^2(\alpha - 1)]; & k_2 &= \varepsilon \frac{2a_1^2 + b_1^2}{2c^2} \tag{A.2} \\
\alpha(e) &= \frac{ab^2}{c^3} \operatorname{arctanh}\left\{\frac{c}{a}\right\} - \frac{b^2}{c^2} & (\text{prolate}), & \alpha(e) &= -\frac{ab^2}{c^3} \operatorname{arctan}\left\{\frac{c}{a}\right\} + \frac{b^2}{c^2} & (\text{oblate})
\end{align}
$$

The two next velocity fields are given by:

$$
\begin{align}
\underline{v}^3 &= \frac{a_1 b_1^2}{b L_\lambda} \left[ f_3(\lambda) \sin^2(\beta) \underline{e}_\lambda + \frac{1}{2} g_3(\lambda) \sin(2\beta) \underline{e}_\beta \right] \cos(2\theta) - \frac{a_1 b_1^2}{b^2} g_3(\lambda) \sin(\beta) \sin(2\theta) \underline{e}_\theta \tag{A.3} \\
\underline{v}^4 &= \frac{a_1 b_1^2}{b L_\lambda} \left[ f_3(\lambda) \sin^2(\beta) \underline{e}_\lambda + \frac{1}{2} g_3(\lambda) \sin(2\beta) \underline{e}_\beta \right] \sin(2\theta) + \frac{a_1 b_1^2}{b^2} g_3(\lambda) \sin(\beta) \cos(2\theta) \underline{e}_\theta
\end{align}
$$

with:

$$
f_3(\lambda) = 1 + \frac{k_3}{b^2} [(3\alpha + 1)b^2 - 2a^2]; \quad g_3(\lambda) = \frac{k_3}{ab} [2a^2 + 3b^2(\alpha - 1)]; \quad k_3 = \varepsilon \frac{b_1^2}{4c^2}
$$

The last velocity fields are:

$$
\begin{equation}
\begin{aligned}
\underline{v}^5 = {}& \frac{a_1 b_1^2}{b L_\lambda} f_5(\lambda) \sin(2\beta) \cos(\theta) \underline{e}_\lambda + \frac{a_1 b_1^2}{b L_\lambda} g_5(\lambda) [1 + k_5 \cos(2\beta)] \cos(\theta) \underline{e}_\beta \\
& - \frac{a_1 b_1^2}{b^2} G_5(\lambda) (1 + k_5) \cos(\beta) \sin(\theta) \underline{e}_\theta
\end{aligned}
\tag{A.4}
\end{equation}
\begin{equation}
\begin{aligned}
\underline{v}^6 = {}& \frac{a_1 b_1^2}{b L_\lambda} f_5(\lambda) \sin(2\beta) \sin(\theta) \underline{e}_\lambda + \frac{a_1 b_1^2}{b L_\lambda} g_5(\lambda) [1 + k_5 \cos(2\beta)] \sin(\theta) \underline{e}_\beta \\
& + \frac{a_1 b_1^2}{b^2} g_5(\lambda) (1 + k_5) \cos(\beta) \cos(\theta) \underline{e}_\theta
\end{aligned}
\tag{A.5}
\end{equation}
$$

with:

$$
f_5(\lambda) = \frac{1}{4ab} [(3\alpha + 1)(1 + k_5)a^2 + 3(1 - \alpha)(1 - k_5)b^2]; \quad g_5(\lambda) = \frac{1}{2}(3\alpha - 1); \quad k_5 = -\varepsilon \frac{a_1^2 + b_1^2}{c^2}
$$

**Appendix B. Components of S(e) for a spheroidal inclusion with axis along e<sub>3</sub>**

$$
\begin{align}
S_{1111}(e) &= S_{2222}(e) = 3S_{1122}(e) = 3S_{2211}(e) = 3S_{1212}(e) = 3\epsilon \frac{2a^2 + 3(\alpha - 1)b^2}{8c^2} \nonumber \\
S_{3333} &= \frac{\epsilon}{c^2}(3\alpha a^2 - b^2); \quad S_{1133}(e) = S_{2233}(e) = (1 - 3\alpha) \frac{\epsilon a_1^2}{2c^2} \tag{B.1} \\
S_{3311}(e) &= S_{3322}(e) = (1 - 3\alpha) \frac{\epsilon b_1^2}{2c^2}; \quad S_{2323}(e) = S_{1313}(e) = \epsilon(1 - 3\alpha) \frac{a_1^2 + b_1^2}{4c^2} \nonumber
\end{align}
$$

Let us recall that $\alpha$ is function of eccentricity $e$ and is given by (A.2).

**Appendix C. Expression of tensor W and P**

$$
Q_r : \mathbb{P} : Q_s = p_{rs} \quad \text{for } r=1,6; \text{ with: } p_{rs} = \frac{(1+g)(f+g)}{f(1-f)} [H_{rs}(e_1) - f H_{rs}(e_2)] \quad (C.1)
$$

where the nonnull $H_{rs}(e)$ are:

$$
\begin{align}
H_{11} &= 3(1 + 3\alpha)(1 - \alpha); & H_{12} &= 3(1 - 3\alpha)(1 - \alpha - \beta); & H_{22} &= 3(3\alpha + 3\beta - 1)(1 - \alpha - \beta) \nonumber \\
H_{33} &= \frac{1}{12}(1 + 3\alpha - \beta)(3 - 3\alpha + \beta); & H_{55} &= \frac{1}{3}(1 - 3\alpha - 2\beta)(3\alpha - 3 + 2\beta) \tag{C.2}
\end{align}
$$
Tensor $\mathbb{W}$ is given by:

$$ \mathbb{W} = \frac{2(1+g)(f+g)}{3(f(1-f))} [\mathbb{S}(e_2) - \mathbb{S}(e_1)] \qquad (C.3) $$

The nonzero components of tensor $\mathbb{P}$ are:

$$ \begin{aligned} P_{1111} &= P_{2222} = \frac{1}{36}(4p_{11} + p_{22} + 9p_{33} + 4p_{12}); & P_{3333} &= \frac{1}{9}(p_{11} + p_{22} - 2p_{12}) \\ P_{1122} &= \frac{1}{36}(4p_{11} + p_{22} - 9p_{33} + 4p_{12}); & P_{1133} &= P_{2233} = \frac{1}{18}(2p_{11} - p_{22} - p_{12}) \end{aligned} \qquad (C.4) $$

$$ \begin{aligned} P_{1212} &= \frac{1}{4}p_{33}; & P_{1313} &= P_{2323} = \frac{1}{4}p_{55} \end{aligned} $$

from which are easily obtained the components of $\mathbb{P}^{-1}$

## Appendix D. Determination of the macroscopic yield locus

Let us introduce the following change of variables $(D, \bar{d}^*) \Rightarrow (A, B)$ in (25):

$$ \Sigma = \frac{\partial \tilde{\Pi}}{\partial A} : \frac{\partial A}{\partial D} + \frac{\partial \tilde{\Pi}}{\partial B} : \frac{\partial B}{\partial D}; \qquad \frac{\partial \tilde{\Pi}}{\partial A} : \frac{\partial A}{\partial \bar{d}^*} + \frac{\partial \tilde{\Pi}}{\partial B} : \frac{\partial B}{\partial \bar{d}^*} = 0 \qquad (D.1) $$

The aim of the following demonstration is to put (D.1) in the form:

$$ \tilde{\Sigma}_{\text{eq}} = \frac{\partial \tilde{\Pi}}{\partial A_{\text{eq}}}; \qquad \tilde{\Sigma}_P = \frac{1}{f} \frac{\partial \tilde{\Pi}}{\partial B} \qquad (D.2) $$

where $\tilde{\Sigma}_{\text{eq}}$ and $\tilde{\Sigma}_P$ have to be determined.

Considering (24), $\tilde{\Sigma}_{\text{eq}}$ and $\tilde{\Sigma}_P$, defined by (D.2) read:

$$ \begin{aligned} \tilde{\Sigma}_{\text{eq}} &= \frac{\partial \tilde{\Pi}}{\partial A_{\text{eq}}} = \sigma_0 f \left[ \frac{\sqrt{1 + f^2 \xi^2}}{f} - \sqrt{1 + \xi^2} \right] \\ \tilde{\Sigma}_P &= \frac{1}{f} \frac{\partial \tilde{\Pi}}{\partial B} = -\sigma_0 \left[ \arcsinh(f\xi) - \arcsinh(\xi) \right] \end{aligned} \qquad (D.3) $$

with $\xi = B/A_{\text{eq}}$. Eliminating $\xi$ in relations (D.3), we find:

$$ \left( \frac{\tilde{\Sigma}_{\text{eq}}}{\sigma_0} \right)^2 + 2(1+g)(f+g) \cosh \left( \frac{\tilde{\Sigma}_P}{\sigma_0} \right) - (1+g)^2 - (g+f)^2 = 0 \qquad (D.4) $$

Recalling that the expression of tensor $\mathbf{A}$ is given by (10) and the one of tensor $\mathbf{B}$ by (23), one has:

$$ \mathbf{A} = \mathbf{D} - f\mathbb{S}(e_2) : \mathbf{d}^*; \quad \mathbf{B} = \mathbf{d}^* + \mathbb{P}^{-1} : \mathbb{W}^T : [\mathbf{D} - f\mathbb{S}(e_2) : \mathbf{d}^*]; \quad \text{with: } \mathbf{d}^* = \frac{D_h}{f} + \bar{\mathbf{d}}^* \qquad (D.5) $$

It follows that Eqs. (D.1) take the form:

$$ \left\{ \begin{aligned} \Sigma &= \frac{\partial \tilde{\Pi}}{\partial A} : [\mathbb{I} - \mathbb{S}(e_2) : \mathbb{J}] + \frac{1}{fB} \frac{\partial \tilde{\Pi}}{\partial B} \mathbf{B} : \mathbb{P} : \{\mathbb{J} + f\mathbb{P}^{-1} : \mathbb{W}^T : [\mathbb{I} - \mathbb{S}(e_2) : \mathbb{J}]\} \\ \frac{\partial \tilde{\Pi}}{\partial A} : \mathbb{S}(e_2) : \mathbb{K} &= \frac{1}{fB} \frac{\partial \tilde{\Pi}}{\partial B} \mathbf{B} : \mathbb{P} : [\mathbb{K} - f\mathbb{P}^{-1} : \mathbb{W}^T : \mathbb{S}(e_2) : \mathbb{K}] \end{aligned} \right. \qquad (D.6) $$

with $\mathbb{J} = \frac{1}{3} \mathbb{1} \otimes \mathbb{1}$ and $\mathbb{K} = \mathbb{I} - \mathbb{J}$, $\mathbb{I}$ being the symmetric fourth-order identity tensor and $\mathbb{1}$ the second-order identity tensor.

For the determination of $\tilde{\Sigma}_{\text{eq}}$ and $\tilde{\Sigma}_P$, it is convenient to expressed $\Sigma : \mathbb{S}(e_2)$ from the first relation into (D.6), or more precisely its deviatoric part and its hydrostatic part $\Sigma : \mathbb{S}(e_2)$. Using the following properties $\mathbb{1} : \mathbb{S}(e_2) : \mathbb{1} = 3$ and $\mathbb{J} : \mathbb{S}(e_2) : \mathbb{K} = 0$, one obtains:

$$ \left\{ \begin{aligned} \Sigma : \mathbb{S}(e_2) : \mathbb{1} &= \frac{1}{fB} \frac{\partial \tilde{\Pi}}{\partial B} \mathbf{B} : \mathbb{P} : \mathbb{1} \\ \Sigma : \mathbb{S}(e_2) : \mathbb{K} &= \mathbb{K} + \frac{1}{B} \frac{\partial \tilde{\Pi}}{\partial B} \mathbf{B} : \mathbb{W}^T : \mathbb{S}(e_2) : \mathbb{K} \end{aligned} \right. \qquad (D.7) $$
Combining the second relation into (D.6) and the second relation in (D.7), one obtains:

$$ \Sigma : \mathbb{S}(e_2) : \mathbb{K} = \frac{1}{fB} \frac{\partial \tilde{\Pi}}{\partial B} \mathbf{B} : \mathbb{P} : \mathbb{K} \qquad (D.8) $$

Finally from (D.8) and the first relation in (D.7), it is readily seen that:

$$ \Sigma : \mathbb{S}(e_2) = \frac{1}{fB} \frac{\partial \tilde{\Pi}}{\partial B} \mathbf{B} : \mathbb{P} \qquad (D.9) $$

which leads to the expression of $\tilde{\Sigma}_P$ (27). We now come to the determination of $\tilde{\Sigma}_{eq}$. Considering the deviatoric part of $\Sigma$ from (D.6), and noticing that $W^T : \mathbb{K} = W^T$, we have:

$$ \bar{\Sigma} = \frac{\partial \tilde{\Pi}}{\partial A} + \frac{1}{B} \frac{\partial \tilde{\Pi}}{\partial B} \mathbf{B} : W^T \qquad (D.10) $$

Replacing $\mathbf{B}$ by its expression deduced from (D.9), one obtains:

$$ \Sigma : [\mathbb{K} - f\mathbb{S}(e_2) : \mathbb{P}^{-1} : W^T] = \frac{\partial \tilde{\Pi}}{\partial A} \qquad (D.11) $$

and then the searched results (27).

## References

[1] A.L. Gurson, Continuum theory of ductile rupture by void nucleation and growth: Part I—Yield criterion and flow rules for porous ductile media, J. Engrg. Mat. Technol. 99 (1977) 2–15.

[2] G. Perrin, Contribution à l'Étude Théorique et Numérique de la Rupture Ductile des Métaux, PhD thesis, Ecole Polytechnique Palaiseau, France, 1992.

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

[4] P. Suquet, On bounds for the overall potential of power law materials containing voids with an arbitrary shape, Mech. Res. Commun. 19 (1992) 51–58.

[5] L. Dormieux, D. Kondo, F.-J. Ulm, Microporomechanics, Wiley, New York, 2006.

[6] M. Gologanu, J.-B. Leblond, G. Perrin, J. Devaux, Recent extensions of Gurson's model for porous ductile metals, in: P. Suquet (Ed.), Continuum Micromechanics, Springer-Verlag, Berlin/New York, 1997, pp. 61–130.

[7] M. Garajeu, J.C. Michel, P. Suquet, A micromechanical approach of damage in viscoplastic materials by evolution in size, shape and distribution of voids, Comput. Methods Appl. Mech. Engrg. 183 (2000) 223–246.

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

[9] M. Gologanu, J.-B. Leblond, J. Devaux, Approximate models for ductile metals containing non-spherical voids—case of axisymmetric oblate ellipsoidal cavities, J. Engrg. Mat. Technol. 116 (1994) 290–297.

[10] J.D. Eshelby, The determination of the elastic field of an ellipsoidal inclusion, and related problem, Proc. R. Soc. Lond. A 241 (1957) 376–396.

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

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

[13] V. Monchiet, C. Gruescu, E. Charkaluk, D. Kondo, Approximate yield criteria for anisotropic metals with non spherical voids, C. R. Mecanique 334 (2006) 431–439.

[14] P. Ponte-Castañeda, P. Suquet, Nonlinear composites, Adv. Appl. Mech. 34 (1998) 171–302.